一种新型的航空重力梯度测量数据调平方法

孙 勇，周 帅

（吉林大学地球探测科学与技术学院,长春130026）

0 引言

航空地球物理测量的特点之一是测点间距远远小于测线间距，出于避免插值失真的考虑，网格化间距通常被定为测线间距的1/4左右。因此，测线与测线之间彼此的约束小，沿测线的低频噪声就会以条带状的形式显示在网格化后的图上。调平是用来去除数据中的低频干扰的常用方法，如移除线性或非线性的漂移、降低机械噪声的干扰等，最终使测量数据接近真实的场值。调平已经成为航空地球物理测量数据预处理工艺中的标准步骤之一，可以为后续的解释工作提供高质量的地球物理数据。

在具体的操作方法上，已有大量学者做出过研究。已有的调平方法可分为统计调平（切割线调平）与经验调平（微调平）。在统计调平这一范畴内，Foster等[1]用一组正交基来拟合测线与切割线上系统误差，假设系统误差缓慢变化，用最小二乘拟合原始数据与用正交基拟合的数据，以此来分离出这种意义下的 “真实”解。Harold等[2]在处理航空磁测数据时，通过最小化交叉点上的值来去除日变场的干扰，并且事先不需要地面基站的测量值。Bandy等[3]在Foster的研究基础之上推导了当测线与切割线在水平面上的投影为直线时的调平方法，简化了计算步骤，其研究结果可被视为Forster提出的最小二乘统计调平的一个特例。Huang等[4-5]在研究航空电磁场数据与航空磁测数据的调平时提出了一种基于每条测线逐步调平的方法。以误差干扰小的切割线为基准，将切割线与邻近的测线结合起来，用参数方程拟合了测线上的需要调平的误差，再利用最小二乘准则求取了切割线与测线在二范数意义上的最小，求得了方程中的参数。Gleave等[6]在研究航空重力数据的调平时发现，在用滤波切除噪声干扰时，容易将与噪声等波长的信号也切除掉，因此其建议将区域场先进行分离再做调平处理。统计调平的优点是可以对大范围的测线进行某种规则下的统一处理，缺点是在统计的意义下难以准确区分真实信号与噪声信号。Davydenko等[7]在处理航空电磁数据的调平时，引入了一种主成分分析法，来设计方向滤波器。这种方法假设切割线中的噪声差异大于测线，从而对方向滤波器的参数做出了设定。Zhang等[8]在Davydenko的主成分分析法上做出了改进，进一步降低了对切割线的依赖，使方向滤波器在没有切割线的地区同样可以应用。微调平专门用来处理统计调平所不能妥善处理的测线，微调平是一种经验调平，针对个别测线可用特定参数的滤波器进行处理，滤波器的参数则根据处理人员对该地区的了解而定。Nuady等[9]设计了一种非线性滤波器，该滤波器后被我国国土资源航遥中心进行了引进及改良。Minty[10]根据航磁数据调平误差与飞行方向的相关性设计出了一种方向滤波器。对统计调平与微调平的详细技术描述，也可以参见商用软件Oasis Montaji Geosoft中关于调平的帮助文档[11]。

本文在上述研究的基础上，针对航空重力梯度数据，提出了一种组合调平方法。重力梯度数据与航磁数据或航空电磁数据相比，受到的交叉测线位置不同、高度不同，飞行方向与载体质量干扰更大，相同之处在于同样包括仪器的漂移干扰。用商用软件Geosoft中的调平模块直接对数据进行处理，往往得不到良好的效果。因此，本文在统计调平的基础之上增加了3个处理步骤，分别是原始数据的去趋势、融合过密的测线与增加方向滤波器，3个处理步骤均在统计调平之前进行。经过实测数据验证，这种组合调平的方法相比单纯进行低通滤波或只进行统计调平可以取得更佳的效果。

1 组合调平的原理

传统的滤波手段会将噪声及与噪声波长相等的信号一起滤除，而组合调平利用沿测线的滤波器保留高频信息，再用网格化的数据进行低通滤波，以保留低频信息。最后，重新采样到原测线，与高频数据合并以得到最终结果。数据流程图如图1所示。

图1 组合调平的数据流程图Fig.1 Data flow diagram of combined leveling

图1左侧分支主要的步骤为去趋势，去趋势可以在空间域与频率域进行。在空间域进行计算时，除了早期的图像法与多项式拟合法之外，多数要进行褶积运算，计算复杂度高，时间成本高。因此，出于方便计算的考虑，大部分算法均在频率域中进行。在频率域中，常用的背景场去除方法包括了向上延拓法、维纳滤波、小波变换等。针对不同测区数据的特点，可以采用不同的方法。在本节的调平处理中，由于不需要准确地分离出区域场值，分离处理是为了避免在处理高频的噪声干扰时损失一部分信号。因此，可选用较为方便的空间域多项式拟合方法来分离趋势场。

以测线l为基本单位，用n阶多项式来拟合的区域场为Test。由于越高阶的项在趋势场中的占比越小，因此权因子b可以构造为[0，1]区间内的形式

式（1）中，x、y代表测点坐标，每一个测点都被单独赋予了一个权值b。

单条测线的重力梯度趋势场拟合值Test可表达为

式（2）中，B为用坐标信息求出的系数矩阵，l为未知数，n为多项式的阶数，m为单条测线上测点的个数。设趋势场的真值为Treal，在二范数最小的前提下求取最小二乘解，可以得到

将L带回到式（2）中，可以得到趋势场的最小二乘计算结果。

图1所示的右侧分支包括网格化、方向滤波与统计调平，其算法流程图如图2所示。

图2的核心步骤是选择合适的滑动窗口对相邻的测线进行平均处理，窗口大小由窗口内测线的条数来定义。相邻测线包含的噪声相类似，因此，可以取经过平滑处理的相邻测线与原始测线差值的平均值作为噪声的近似值，再用原始数据减去这个近似值，以作为输出。在图2中，统计调平模块的编制参考了胥值礼等[12]对航磁数据调平的研究成果，其核心思路是使切割线与测线的值相接近。其步骤是首先调切割线，其次用调平后的切割线调平测线。设Te为修正前切割线的重力梯度值，Tl为修正前测线的重力梯度值，Te′为修正后切割线的重力梯度值，Tl′为修正后测线的重力梯度值。切割线的调整方法为

图2 图1右侧分支的算法流程图Fig.2 Algorithm flow chart of the right branch in Fig.1

在式（3）中，i为切割线的编号，j为测线的编号，Tcei，k表示第i条切割线上与测线相交的第k个点的重力梯度值，p1为第i条切割线上与测线的交点的个数，Tclj，k为k点处测线上的重力梯度值。因此，式（3）中ΔTei的物理意义为第i条切割线在与测线相交的各点处差值的算数平均值。测线的调整方法为

式（4）中，ΔTlj为第j条测线与切割线相交处差值的算数平均值，Tce′i，k为在第k个交点上调整后切割线上的重力梯度值，Tclj，k为第k个交点上测线的重力梯度原值。经过式（3）与式（4）的处理，可以拉近切割线与测线的差异，最终完成统计调平处理。

2 实测数据验证

以2003年在加拿大西北部的Diavik钻石矿区进行的航空全张量重力梯度测量为例，来验证调平方法的有效性。图3（a）展示了以金伯利岩为主的矿区，图3（b）展示了测线轨迹，切割线共有13条。

图3 矿区位置图Fig.3 Map of mine location

图 4（a）～图 4（f）展示了调平之前I1～I3、C1～C3数据的形态。可以看出，在未经调平处理时，数据包含了很大的与测线相关的方向噪声。为了方便说明调平的效果，在图4中，用红色线框圈出了相应的计算区域。

在图4中，I1～I3、C1～C3的表达式分别为

图4 各分量在调平之前的形态Fig.4 Shape of each component before leveling

在式（5）中，u、v、w为伞形坐标系下三个坐标轴的方向。为与本节提出的调平方法所得到的结果作出对比，首先需对数据进行常规的低通滤波处理。以图 5（a）所示的分量C1为例，经过图5（b）所示的截至波数为 1.84（单位 1/km）的空间域低通滤波，可得到如图5（c）所示的滤波结果，其与原数据的差如图5（d）所示。可以看到，大尺度的低通滤波可以去除相当多的干扰，滤波前后信号的最大值从165.8E降低到了146.9E。但对过滤之后的图像进行分析，可以看到滤波过程也滤掉了一部分有效的信号，即滤掉的部分不完全是条带状的干扰。

下面用本文提出的组合调平方法对C1进行处理，趋势场的分离结果如图6所示。

从图6可以看出，经过去趋势的处理，大部分沿测线方向的噪声都已经被去除，但在测区的北部及南北向的切割线上，仍能看到沿测线方向的强干扰的噪声。

图5 对C1进行低通滤波的结果与相应的滤波器参数设置Fig.5 Results of low-pass filtering for C1and the corresponding filter parameter setting

图6 C1去趋势前后的对比Fig.6 Comparison of C1before and after detrending

用此方法进行处理的关键之处在于相邻测线的框定范围，即选定合适的窗口大小。在垂直于测线方向的窗口分别为15条测线、10条测线、5条测线这3种不同的条件下，对C1的处理结果如图7所示。

图7 不同窗口大小对C1进行滑动的平均结果对比Fig.7 Comparison of sliding average results for C1with different window sizes

由图7（a）～图 7（c）可知，当窗口为 5 时，可以取得最佳的效果。此时，仍能观测到沿测线方向的噪声，因此需要做进一步的统计调平处理。对I1～I3、C1～C3做出相同大小的测线数据集与切割线数据集，之后统计交叉点上的信息，用式（3）与式（4）进行计算，得到的处理结果如图8所示。

图8 经过统计调平的各道数据集Fig.8 After statistical adjustment of each channel data set

图9 滤波与组合调平的结果对比Fig.9 Comparison of the results of filtering and combined leveling

滤波与组合调平的结果对比如图9所示，组合调平的结果上异常形态更加完整，沿测线方向的噪声也得到了更好的压制。组合调平的结果异常幅值最大值为63.7E，这说明简单进行低通滤波处理的结果不仅没有保留完整的异常形态，而且噪声也没有被完全去除。

3 结论

航空地球物理测量的高动态性使其无法避免被与航向有关的噪声干扰，如何在去除干扰的同时保留尽可能多的异常细节是处理过程中的难点，特别是对于高分辨率的航空重力梯度测量数据而言。与低通滤波处理或传统的切割线调平处理不同，本文提出了一种组合调平方法。经过实测数据验证，该方式可以明显提高原始数据的质量。从图9可以看出，在经过处理后，仍有少量非地质体信号的成分存在。因此，建议在组合调平后进行逐条测线的微调平处理，以进一步提高数据的质量。