基于约束相干增幅扩散滤波的航空磁测
导读:
地球物理数据必须网格化才能进一步处理。然而，航测数据通常是在沿测线方向高频采集，测线间距较大，只有当网格尺寸接近样本间距时，网格化才能充分利用测线数据，否则，高频信息总是会丢失，从而产生混叠效应。例如，与测线呈锐角的线性趋势在图像中会呈现“牛眼”，类似于相交线处的串珠状异常。这些伪影的存在可能会掩盖目标异常，导致对目标体形态的错误解释。本文评估了一种名为“约束相干增幅扩散滤波”的方法，该方法仅在强各向异性的特定区域内平滑图像，本文用组合模型和实测数据进行了方法的测试。结果表明，该方法可以有效地以各向异性为基础沿着构造方向进行平滑。处理后的图像能够得到显著改善，例如垂直梯度图。施加的约束保证了原始测线有用信息得以保留。将此方法应用于实测数据，与测区传统的仅执行同等扩散平滑，而不考虑各向异性的处理结果进行了比较，证明本文提出的方法结果更优的同时压制了混叠效应。
0 引言
地球物理数据一般情况下是沿测线方向进行高频采集，测线间距相对点距较大，这种采样间距不一致的问题在航空勘查中尤为明显，导致测线数据进行网格化时在测线间进行插值，从而抑制小幅度异常信号，造成细微趋势信息的丢失。由于高频信息总是丢失，与测线呈锐角的线性趋势在测线交叉处往往会产生“牛眼”等伪影。根据涉及的波长，这些牛眼在图像中呈现为“串珠”或“阶梯”效应，对于这种狭长的磁异常，当测线与地质走向垂直时，几乎不产生影响，但随着测线方向与线性趋势之间角度的变化，这一问题会加重，在导出垂直梯度图或其他趋势增强图时，混叠效应会加剧，从而使地质结构的解释变得困难。
一般而言，混叠程度与测线间距、飞行高度或场源深度成正比，难以在保留原始数据中不同波数信息的前提下进行较好的图形处理。测量数据中横断面的混叠程度可以通过棱柱模型进行量化，三维测量对飞行高度与测线间距的比率做了严格限制，将混叠效应影响最小化。Reid得出结论，在区域总磁场测量中，测线间距不应超过飞行高度的两倍，而在梯度测量中不应超过一倍飞行高度，这些限制在二维尺度中可以酌情减小，但与测线方向成锐角的线性趋势仍然会呈现“牛眼”。
航磁数据处理中一般要求数据网格以测线数据为基准，这是因为经过适当处理后，能使测量精度高于±0.1nT，但是观察值相对精度可能大三到四个数量级。观察值相对精度是指航磁数据的网格算法在周边的±0.1nT内生成控制网格。对于大量航磁数据，有三种主流的算法满足维持原始数据有效性的前提下在测线之间产生平滑插值和高效的标准，包括最小曲率网格化、双向网格化和局部分段多项式网格化。最小曲率技术通过解决双调和差分方程生成具有最小总曲率的平滑表面，具有许多优良特性，因此广泛用于地球物理数据处理中。双向网格化方法沿线插值，然后横向重插值创建正方形网格，通过坐标系旋转增强趋势，然而，该技术不能充分处理多重或局部趋势。局部分段多项式网格化使用邻域点来拟合趋势面，并根据其与邻域中心的距离加权，这种插值会产生几个理想特性：在平滑区域中，网格具有较高精度；数据始终回归到中心，具有最佳收敛性；若存在不连续数据，不会出现线性插值的吉布斯现象。此外，等效层技术也可以用于数据网格化，但在网格化超大数据时效率太低。这些网格化方法都有一个共同的问题：与测线成锐角的线性趋势往往会在交叉处产生“牛眼”，它源于测线间数据的欠采样和网格化方法无法预测细微趋势特征。
增强趋势的方法尤为重要，前人已提出了几种用于增强航磁图中趋势的技术。使用梯度测量可以显著改善与测线呈锐角的尖锐特征和狭长磁异常的趋势。Yunxuan采用有限Radon变换（RT）应用于模拟重力数据去除带状伪影，但仅限于单一趋势方向。Sykes和Das进一步发展了这个想法，使用二维RT来增强地球物理图中的线性特征，但是RT在数据重构中维持数据精度的能力有限，限制了其在航磁数据中的应用。GoldenSoftwareSURFER8提供了SURFER各向异性的选项，提供灵活的搜索椭圆，可以进行权重和角度设置，并且可以集成于最小曲率网格化等网格化方法，能够显著改善网格化效果。最近，Guo等提出使用逆插值方法来网格化航磁数据，该方法在产生平滑图像的同时，一定程度上抑制混叠效应。
在网格化时引入了额外数据来增强趋势也是可行的方法。Fitzgerald等描述了一种根据最小方差原则插入额外数据点的方法，但该算法高度依赖所选参数，且在与测线夹角大于45°的特征上应用效果不佳。Keating提出通过使用邻近分析插入额外数据点的一种方法，能够有效去除了带状伪影，在某些情况下，当几个极大值无法连接或极大值和极小值的趋势线相交时，仍需进一步处理。Zhang等提出了一种网格化参数提取方法设置最佳插值节点，提高网格化效率，但是该方法面对大量航磁数据时运行较慢且不适用。
本文聚焦Smith和O’Connell的研究，通过相干增强扩散（CED）滤波改善成图效果。Smith和O’Connell使用约束各向异性扩散滤波（“TRUST”网格化）技术来增强航磁数据中的趋势。这种方法以测线数据为基准，能有效增强多个方向的趋势。然而，将网格转换为增强图（例如分析信号幅度）时，由于各向异性扩散滤波不仅作用于相干趋势，还作用于非相干结构，致使沿测线出现小的高频伪影。Smith和O’Connell改进了他们的方法，抑制沿测线生成的高频伪影，遗憾的是Smith和O’Connell所称的“改进”并未有详细记录。本文通过应用约束CED滤波弥补了Smith和O’Connell方法的不足，约束扩散仅在相干结构上作用，而不影响其他区域。
相干增强各向异性扩散是由Weickert首次提出，旨在增强流动状纹理的相干性的技术。这种技术已应用于多个领域，如指纹和医学图像的增强。本文将相干增强各向异性扩散应用于航磁数据的网格化中。
各向异性扩散滤波根据以下方程对初始图像（网格化的航磁数据）进行处理：
∂uσ/∂x和∂uσ/∂y分别是u(x,y)在x和y方向上的导数；而∇uσ表示高斯平滑后的u(x,y)的梯度。
Gρ由方程4定义，其中σ被ρ取代。积分尺度ρ决定网格的平均程度，局部尺度σ表示边缘检测忽略小于σ的细节噪声的尺度。两者是不同的平均距离，后续在参数选择部分详细讨论。
对称矩阵Jp具有正半定性，其特征值是实数，且不小于0。Jp的正交特征向量集{v1,v2}由以下方式给出：
Jp的特征值为结构的相干性提供了有效信息。最大的特征值μ1与特征向量v1相关联，后者对应与地质场源走向垂直的方向。相反，最小的特征值μ2与特征向量v2相关联，后者对应于地质场源走向的方向。一般而言，没有相干结构平滑区域的特征是μ1≈μ2≈0；结构边缘处μ1≫μ2≈0，而中间部分的特征是μ1≥μ2≫0。因此，两个特征值的差异可以作为相干性的度量，即实际的各向异性度量。各向异性系数k定义为：
随着特征值的显著变化，相关系数增加，在各向同性结构中趋于零。Weickert创新地提出在结构张量Jp上进行修改，使图像的最大滤波沿着趋势方向进行，最小滤波则垂直于趋势方向，通过重新分配两个特征值实现这一改进。重建的矩阵是Smith和O'Connell在TRUST方法中定义的扩散张量D，默认μ1=0和μ2=1。
Weickert建议，相干增强各向异性扩散的扩散张量D参考与结构张量Jp(∇uσ)相同的特征向量，且其特征值为：
引进公式中的指数函数是为了确保扩散张量随k的变化时保持平滑性，λ2不会超过1，那么扩散张量D是各向异性的。正参数α(α∈[0,1])是可交互数值，在扩散的计算中用以保证连贯性：即使结构变得各向同性(k→0)，具有轻微线性扩散的情况仍然存在。参数c也是可交互的数值。
为了在数值上逼近CED，导数使用有限差分替代，则方程1的最简单离散化表示为：
其中Δt是时间步长。Δt小于0.25时，时间步长算法有效。n=0,1,2,...表示时间步长。初始位场定义为u0。各向异性扩散量un+1通过上述方程获得。Weickert和Weickert和Scharr提供了更多实现相干增强各向异性扩散的详细资料。由于方程14中的扩散张量Dn的计算耗时较长，且Dn随时间步长的变化不大，因此在若干时间步后重新计算Dn，以提高计算效率，采用Smith和O'Connell的程序。
上述模型中的多个参数在实际应用中需要进行调试。为了满足增强航磁数据趋势的要求，在本节中提供了经验参数的选取指南。
方程5表明卷积的积分尺度ρ决定了梯度项的平均值。为了消除来自Jp的空间混叠效应，尺度ρ宜取测线间距的1.5倍，确保只有两条或更多条线上的特征得到增强。在方程3中，局部尺度σ反映噪声尺度，作用是使边缘检测器忽略小于σ的细节。σ的值通常比ρ小。实际上，取测线间距0.5到1.0倍之间的值会有更好效果。Smith和O'Connell使用Gρ=1和滤波后的数据网格u(x,y)进行二维卷积滤波，其宽度约为测线间距的1.5倍。鉴于各向异性的方向至关重要，局部尺度过大（例如，1.5倍的测线间距）可能导致计算错误，针对因地制宜地进行不同平均过程比Smith和O'Connell的做法更加灵活。
为了保留原始测线的有用信息，穿过网格节点的测线的值不能通过扩散滤波修改。
锚定或约束这些固定值是解决此问题的一种方法。在实际工作中，通过固定测线半个网格单元内的网格节点实现。一系列实验表明（未展示数据），即使α非常小（例如0.001），在平滑区域仍然会出现一些小的伪影，为了减少其数量，对结构张量D的特征值做修改：
这些特征值不包含相干性的任何信息。由于λ2始终等于1，即使各向异性不强的情况下，扩散也在全域得以应用。各向异性的方向至关重要，应该在高各向异性方向执行扩散步骤。所以，计算过程是基于5×5模板的不同有限算子（∇和∇·）执行。
图2 (a)对图1a中的图像应用约束CED滤波器（300个时间步）。
所有角度上的线性结构被平滑成像。
(b)从趋势增强网格(a)派生的垂直梯度图。
此图表明其线性特征的连续性得到了显著改善。
图3 (a)图1a和图2a数据之间的差异
扩散主要发生在线性特征和块体边缘处。
(b)脊线剖面上的原始数据（图1a中的实线）以及经过300个时间步长的约束CED处理后的数据（虚线）。星标处的数据保持不变。高频起伏已被消除。
其中m表示每条测线的测量点数，l表示测线编号。由于相干特征的增强导致磁场在趋势方向上的梯度趋于零，则δn的值在大约300个时间步长后趋于零，相应的扩散系数λ2也趋于零。针对这种情况，可以选择两个不同的标准来终止算法的迭代，一是最大迭代次数，二是当磁场的平均差异δn达到预定值时终止迭代。从图中可以发现，当扩散张量元素重新计算时，变化值会突然增加（图4）。
图4 约束CED经过0、30、...450次迭代后磁场平均差异δn
平均差异在大约300个时间步长后趋于零。
图5 由图2a中的趋势增强磁场图生成的特征值差异图
此算法比最初的计算复杂度要低：对于300次迭代并在每30个时间步长后重新计算扩散张量元素，101x81网格数据的执行时间为2秒，512x512网格数据为12秒，1024x1024网格数据为43秒，2048x2048网格数据为253秒，使用的是2.27GHzCPU的PC。
5 实测案例
图6展示了一个7000x7000m2的小型航磁测量项目，包含多个狭窄场源。测线方向为东西方向，测线间距为200米，每17米采样一次，距地面高度为80米。图6a的网格使用50米间隔的双向样条插值生成。区域的地质体趋势多为北西向，多个狭长的表现为 “串珠”状的北西向线性异常。将约束CED算法应用于图6a，C的值设定为2000；ρ和σ分别等于测线间距的1.5倍和0.5倍，并在300个时间步长后可以得到图6b，显然，图6b的趋势增强磁场图中的“串珠”状异常得到了改善，中心东-西区域的趋势也更加清晰。
通过减小各向异性阈值参数C可以凸显低振幅异常，但在本例中，采用较大的C值以减少各向异性扩散滤波，进而保留可能是排水或古排水造成的低振幅信息。此外，由于网格数据受到测线数据的约束，图6b中的趋势增强网格保持了与原始数据异常特征的一致性。
从垂直梯度图（图7b）可以看出，由趋势增强网格计算获得的线性异常，比从原始网格计算的（图7a）连续性更好，垂直梯度图也更易于解释，并且回避了断层或间断。增强磁场图（图6b）中的特征值差异图（图8）也凸显了主要异常边缘位置的重要信息，有助于地质解释。
本文将C=2000的算法与C=0的算法进行了比较，C=0的算法最接近Smith和O’Connell的方法。结果如图9a所示，相应的垂直梯度图如图9b所示。图9a中，使用TRUST网格化的线性特征得到了充分增强，但垂直梯度图测线交叉处仍然存在一些低振幅伪影。Smith和O’Connell认为这些低振幅伪影与测线附近地表变化的限制有关。除与测线相交的低振幅伪影外，图中的次生区域也存在低振幅伪影，可能是排水或古排水系统造成的。本文提出的约束CED滤波方法可以自动检测这些区域的各向异性并选择性地进行扩散滤波，从而保留低振幅信息，所有趋势增强方法都应满足的这一前提条件，确保获得高质量的航磁测量图。
实测数据测试成果斐然
相较TRUST方法，特征参数的不同选项作用于各向异性扩散滤波，使其仅应用于具有较强各向异性的区域。将约束相干增强各向异性扩散方法应用于组合模型和实测数据都得到良好的平滑效果。
相比将各向异性滤波应用所有区域的情况，约束相干增强各向异性扩散网格在垂直梯度图上有效压制伪影的产生，其次，新方法生成的特征值差异图还突出了地质体边缘的有效信息。因此，约束相干增强各向异性扩散能够更准确地从原始数据中提取信息，有利于后期地质解释。
