任朗宁，张凤旭，郝梦成，李银飞

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

0 引言

垂向导数作为常规的位场数据处理方法，不仅在分离叠加异常、突出局部异常和识别地质体边界等方面具有重要作用[1--3]，还是其他一些位场数据处理方法的重要环节。例如欧拉反褶积、向下延拓[4]和重力归一化总梯度等方法都需要先进行垂向导数的换算，并且导数换算的结果精度将直接影响到上述位场处理方法的最终计算精度。由于导数换算对噪声的敏感性，常常得不到较好的求导结果，影响后续异常的处理与解释。为此，许多改进位场垂向导数换算的方法被提出，它们大致可分为两类：一是空间域法；二是波数域法。空间域法包括量板法、边界单元法和样条函数法等[5--7]。波数域法包括常规FFT法、维纳滤波法[8]、补偿圆滑滤波法[9]和正则化方法[10]等。

为改善导数换算的不稳定，常规做法是通过引入低通滤波器，将其与垂向导数算子相结合，构建新的垂向导数算子，来压制垂向导数计算过程中噪声的放大作用。Butterworth低通滤波、Hanning窗滤波以及向上延拓滤波是3种常见的低通滤波，在导数计算中取得了较好的效果，但各方法的滤波参数常常需要反复尝试，降低了计算效率，笔者利用相关系数法来选择合适的滤波参数，实现位场数据3种方法快速垂向导数换算，降低了人为因素对导数结果的影响。

1 方法原理

1.1 垂向导数中常用的低通滤波因子

位场G(x,y)与其各阶垂向导数VDm(x,y)(m代表阶数)在波数域的关系式为[11]：

(1)

VDm(u,v)=φm·G(u,v)

(2)

由前人研究可知[12]，φm是一个放大型因子，求导过程会对高频噪声放大，导致垂向导数结果精度变差。通常的做法是在常规导数因子中附加一个低通滤波器来减小φm在高频部分的放大作用，保留低频有效信息。下面给出三种常规的低通滤波方法：

①Butterworth低通滤波器是一个在通带内拥有最大限度平坦度的滤波器，其频率响应为：

(3)

式中：ω为频率；ωb为Butterworth滤波器的截止频率;N为滤波器阶次。

②Hanning窗低通滤波是位场信号处理中常用的一类窗函数滤波因子,其表达式[13]为：

(4)

式中：ωh为Hanning窗滤波器的截止频率。

③向上延拓是将观测平面上的异常换算到观测平面之上，相当于一个低通滤波[14]，在波数域中，向上延拓因子是:

(5)

式中：h是上延高度。

由上述的式(2～5)可各自构造出新的求导公式：

VDB(u,v)=B(ω)·VDm(u,v)

(6)

VDH(u,v)=H(ω)·VDm(u,v)

(7)

VDY(u,v)=Y(h)·VDm(u,v)

(8)

将结果进行Fourier反变换，即得到垂向导数结果。

从上式及前人研究结果可知，这些方法的滤波效果主要与低通滤波器滤波参数(截止频率或上延高度)的选择有关。以截止频率的选择为例，当截止频率过小时，求导导致低频部分的有效异常损失；当截止频率过大时，求导过程中对噪声的压制作用较弱。这两种情况都会导致导数换算结果的精度不高，因此需要选择一个合适的截止频率来平衡这两种状况，即在保留有用信号的同时，最大限度的减少噪声干扰。

1.2 基于相关系数法的滤波参数选取方法

为了快速确定滤波参数，这里借鉴曾华霖确定最佳上延高度时采用的相关系数法[15]。即利用取不同滤波参数时，求导结果之间的相关性来确定最佳滤波参数。具体做法是：

①求原始异常的Fourier变换谱；

②在一定范围内，设置n个离散的滤波参数，并计算相应的垂向导数；

③利用下列公式计算相邻两个滤波参数求导结果的相关系数；

(9)

④当相关系数曲线开始趋于稳定的时候所对应的参数确定为最佳滤波参数。

2 模型试验

为检验滤波参数选取的相关系数法在垂向导数计算中的效果， 本文设计包含4个大小埋深各不相同的长方体的组合模型进行试验。 模型参数见表1。

表1 模型参数

在实际中，由于高频噪声干扰的始终存在，故对该组合模型产生的重力添加均值为0，标准差为0.1 g.u.的高斯白噪声，其含噪声重力异常Δg见图1a。从图1a中无法准确识别地质体的真实边界位置，需要对模型重力异常进行处理。图1b是该模型的理论垂向二阶导数。对比图1a和图1b，可以看出异常的垂向二阶导数换算，对地质体的异常信息进行突出。

对含噪声重力异常(图1a)采用式(6～8)进行垂向二阶导数计算。其结果Δgzz见图1d、图1f和图1h。图1c、图1e和图1g分别是Butterworth低通滤波、Hanning窗滤波以及向上延拓滤波进行低垂向二阶导数计算的滤波参数选取图，为验证相关系数法选取的滤波参数是最佳的，本文还求取了各方法导数计算结果与理论垂向二阶导数之间的均方根误差。从图1c、图1e和图1g中可以看出，Butterworth低通滤波、Hanning窗滤波以及向上延拓滤波的滤波参数分别是ωb=3.413 8，ωh=5.862 0和h=0.206 9，其所在位置与均方根误差最低的位置相对应，说明利用相关系数法确定滤波参数是合理的。对比图1d、1f和1h可以看出，经过本文选取合适滤波参数，各方法垂向导数计算中的高频噪声均得到了有效压制，图1d和1f中的结果相较于图1h中的处理结果更好，能够准确地看出异常体的位置和形态，与图1b中的结果相似。图1d和1f中的导数结果虽然都能够较好显出异常体的位置和形态，但是图1d异常的幅值与理论值(图1b)的差距较小，结果精度较高。

为了进一步检验各方法的垂向二阶导数效果，对导数结果进行二次积分反算，并将积分结果与原异常求残差，其结果分别见图2a、2c、2e和图2b、2d、2f。对比图2a、2c、2e和图1a，可以看出图2a中Butterworth低通滤波的二次积分结果与原异常最相近。再对比图2b、2d、2f，可以看出图2b中的残差结果里包含的地质体异常信息最少，而Hanning窗滤波和向上延拓处理的结果中异常信息较多，说明这两种方法在垂向二阶导数计算中对异常的转换不够彻底。综上所述，使用相关系数法选取最佳的滤波参数的做法是合适的，利用确定的滤波参数获得垂向二阶导数结果中，Butterworth滤波的导数精度最高。

3 实际应用

为了检验上述方法在实际数据处理中的效果，对位于大兴安岭外围的松辽盆地镇赉地区的布格重力异常数据(图4a)进行了垂向二阶导数计算。数据网格点距0.5 km，网格点数124×231。首先，利用相关系数法选取本文各方法最佳的滤波参数分别为ωb=0.9，ωh=0.9，h=1.8(图3a、3b、3c)。再将选取的滤波参数代入到各个滤波器中进行垂向二阶导数计算，结果见图4b、4c、4d。对比上述结果可以看出，Butterworth低通滤波的垂向二阶导数(图4b)处理效果较好， 其异常曲线圆滑，对噪声的压制作用较好，同时对浅部异常的突出作用也较强；Hanning窗滤波的垂向二阶导数结果(图4c)精度较差，浅部异常较宽缓，说明相邻异常的分离不够彻底。而向上延拓滤波的垂向二阶导数(图4d)虽然突出了大部分的浅部异常特征，但是异常曲线明显存在畸变现象，说明噪声干扰被压制的不够彻底。为了进一步检验求导结果的准确性，对各方法垂向二阶导数结果求取二次积分并与原异常做差，结果分别见图5a、5c、5e和图5b、5d、5f。二次积分的结果与该区原异常进行对比，可以看出Butterworth低通滤波的结果(图5a)与原异常形态相似，幅值相近，而其他两种方法的异常形态与原异常存在较大差异，幅值也有不同。从图5b、5d、5f的残差结果中可看出，图5b中的异常特征最少，幅值最低，说明Butterworth低通滤波求取垂向二阶导数过程中异常转换地较为彻底，间接证明了其垂向二阶导数的计算结果比其他两种方法的结果精度高，这与模型试验的效果相似。

(a)含噪声重力异常；(b)理论垂向二阶导数；(c)Butterworth滤波参数选择；(d)垂向二阶导数(Butterworth滤波器)；(e)Hanning窗滤波参数选择；(f)垂向二阶导数(Hanning窗滤波)；(g)向上延拓滤波参数选择；(h)垂向二阶导数(向上延拓)。图1 模型重力异常及各方法截止频率处理效果对比分析图Fig.1 Comparative analysis of model gravity anomaly and processing effect of cutoff frequency of each method

(a)Butterworth二次积分结果；(b)二次积分与原异常残差(Butterworth滤波器)；(c)Hanning窗滤波二次积分结果；(d)二次积分与原异常残差(Hanning窗滤波)；(e)向上延拓二次积分结果；(f)二次积分与原异常残差(向上延拓)。图2 各方法求导结果二次积分及残差图Fig.2 Secondary integration and residual plot of each method

(a)Butterworth截止频率选取；(b)Hanning窗截止频率选取；(c)向上延拓上延高度选取。图3 松辽盆地镇赉地区各方法截止频率选取图Fig.3 Selection chart of cutoff frequency for each method of Zhenlai area, Songliao Basin

(a)西部斜坡布格重力异常；(b)Butterworth滤波器垂向二阶导数；(c)Hanning窗滤波垂向二阶导数;(d)向上延拓滤波器垂向二阶导数。图4 松辽盆地镇赉地区重力异常及垂向导数处理图Fig.4 Gravity anomaly and vertical derivatives processing map in Zhenlai area, Songliao Basin

(a)Butterworth二次积分结果；(b)二次积分与原异常残差(Butterworth滤波器)；(c)Hanning窗滤波二次积分结果；(d)二次积分与原异常残差(Hanning窗滤波)；(e)向上延拓二次积分结果；(f)二次积分与原异常残差(向上延拓)。图5 松辽盆地镇赉地区各方法求导结果二次积分及残差图Fig.5 Secondary integration and residual map of results of Zhenlai area, Songliao Basin

4 结论

(1)使用相关系数法确定了Butterworth低通滤波、Hanning窗滤波和向上延拓滤波的最佳滤波参数，实现了3种方法位场数据垂向导数的快速计算，降低了人为因素对计算结果的影响。

(2)模型试验结果表明，使用相关系数法选取最佳的滤波参数的做法是合理的，利用确定的滤波参数计算得到的垂向二阶导数结果中，Butterworth低通滤波法的导数结果与其他两种方法的结果相比，精度较高。

(3)松辽盆地镇赉地区的实测布格重力异常数据处理检验了3种垂向导数换算方法的实际效果，获得了与模型试验相似的效果，为进一步数据处理与解释提供了支持。