重力异常及其梯度张量DEXP定量解释方法的影响因素分析

邱峰，杜劲松，陈超

(1.中国地质大学(武汉) 地球物理与空间信息学院，湖北 武汉 30074; 2.中国地质大学(武汉) 地球内部多尺度成像湖北省重点实验室，湖北 武汉 30074; 3.中国地质大学(武汉) 地质过程与矿产资源国家重点实验室，湖北 武汉 430074)

0 引言

重力勘探在油气与矿产资源勘查、地质填图、考古探测、工程与环境等地质学科中发挥着重要作用。随着测量手段的多样化和测量数据精度的提高，欲在一定精度的数据基础上，达到理想的处理和解释效果，则需要与数据相对应的新方法与新技术。因此，新的处理和解释方法是地球物理学者长期以来一直致力于探索和研究的重要问题。伴随着全(或部分)重力梯度张量数据越来越多地应用于各个领域，其处理和解释的方法和技术研究也己经成为新的重要课题。

快速成像法是相对于传统常规反演方法而提出的[1-3]，由于不需要任何的迭代过程，因此能够快速地得到地下异常体的分布状况或相关参量(如深度、倾角及构造指数等)，避免了传统方法耗时长、计算量大等缺陷。Fedi和Pilkington于2012年总结并比较了位场成像方法[4]，包括Cribb[5]于1976年提出的Cribb成像法、Moreau等[6]于1999年阐述的连续小波变换法、Mauriello等[7]于2001年提出的概率成像法、Zhdanov等[8]在2011年提出的偏移成像法以及Fedi[9]于2007年提出的极值点深度估计法(DEXP)等。为了将这些方法应用于重力梯度张量数据解释，国内外学者们将许多方法均进行了一定程度的改进(如Guo等[10]; Zhdanov等[11]; Zhou等[12])。

DEXP方法是一种重要的成像方法，该方法由Fedi在2007年提出并对其进行了详细的阐述[9]，由于在计算过程中深度加权函数引入了构造指数，使得与其他方法相比具有许多优点。该方法的适用性比较广，不仅可以用于重力异常不同阶次导数以估计异常体的物性参数[9](如场源深度、剩余密度及构造指数等)，同时也是一种解释自然电位数据的可靠方法[13]。Abbas等[14]又在2014年提出了一种利用两个不同阶导数场之间的比值解释位场数据的自动DEXP方法，该方法在成像过程中不依赖于异常体的构造指数。此外，DEXP方法还可以与其他方法结合使用(Fedi等[15])。之后，Zhou等[12]又进一步将DEXP成像法用于重力梯度张量数据处理中，取得了较好的成像结果。笔者根据DEXP变换理论，通过设置多组理论模型及其模拟数据进行试验研究，结合实际应用，着重讨论了DEXP成像方法的抗噪能力以及数据的点距、计算范围和背景场对DEXP变换结果的影响。

1 方法原理

对于可视作质点或均匀球体的重力异常为：

(1)

其中：G为万有引力常数(G=6.67×10-11m3·kg-1·s-2)，r为坐标原点指向测量点P(x,y,z)的矢量，r0为坐标原点指向质点或球心Q(x0,y0,z0)的矢量。在笛卡尔直角坐标系下，假设质心或球心的位置为(0,0,z0)，测量点的水平位置为x=x0=0、y=y0=0以及z>0，令质量M=1，并且将万有引力常数进行归一化处理，则式(1)可简化为：

(2)

进一步地，式(2)的任意p阶导数fp可以写为：

(3)

事实上，Fedi[9]以及Abbas等[14]均已说明了DEXP方法的关键原理是尺度函数τp，即为位场函数的对数对深度z的对数的导数，故而根据式(3)，τp可表示为：

(4)

其中：z0为场源的中心埋深，Np为重力异常p阶垂向导数对应的构造指数，定义为Np=N+p。其中N为构造指数，是与场源几何形态相关的参数[16]。

根据式(4)，Fedi[9]将重力异常p阶垂向导数的DEXP变换函数Wp定义为：

Wp=zNp/2fp，

(5)

并且证明了DEXP变换函数的极值点可用于估算场源深度。

重力梯度张量的各个元素为重力位U在笛卡尔直角坐标系下三个坐标轴方向的二阶导数，即

(6)

将式(6)写成矩阵形式，即：

(7)

由于重力异常位在场源外部空间满足拉普拉斯方程，即Δ2U(r)=Txx+Tyy+Tzz=0，并且矩阵T为对称矩阵，故重力梯度张量仅有5个独立元素。

Fedi[9]指出DEXP变换方法适用于重力异常任意p阶导数，并且已给出了重力异常gz以及重力异常垂向导数Tzz的DEXP变换公式，故而DEXP变换同样可以应用于重力梯度张量[12]。由于重力梯度张量的各个元素为重力异常各分量在x、y以及z方向的一阶导数，所以Np=N+p可以写为N1=N+1，因此重力梯度张量各分量的DEXP变换函数为：

(8)

对于重力梯度数据，不同分量反映了地下异常体的不同特征信息，比如Txz、Tyz主要反映的是异常体的边界信息，而Tzz则可以描述场源中心位置信息。因此，重力梯度张量的DEXP变换也具有相似的特征。Zhou[12]指出Wxx、Wyy和Wzz可以用来描述地下异常体的中心埋深、Wxz和Wyz则可以根据不同极值点深度位置确定场源的边界信息，同时也指出DEXP变换的总水平导数Wthd能够反映异常体的总体边界信息，其计算公式为：

(9)

2 基于理论模型的影响因素分析

首先令均匀直立棱柱体模型的长(a)、宽(b)与高(c)分别为30、20与10 m，其中心埋深为15 m，剩余密度Δρ为0.5 g/cm3。首先， 通过正演计算可以得到该模型在地表即z=0处产生的重力异常以及重力梯度张量异常；其次，分别将模拟数据在波数域进行不同高度的向上延拓处理以得到3D数据体；然后，根据DEXP变换理论，分别可得模型重力异常以及重力梯度张量异常的DEXP变换结果，分别如图1与图2所示。其中，模拟数据的测点点距为2 m，向上延拓总高度为40 m，步长为1 m。图中黑色矩形描述了直立棱柱体模型的边界位置，白色实线为DEXP变换后所选取的剖面位置，白色圆点为DEXP变换对应的极值点。

图1 直立棱柱体模型的重力异常(a)及其DEXP变换结果(b)Fig.1 Gravity anomaly (a) of the rectangular prism model and its DEXP transform results (b)

图2 直立棱柱体模型的重力梯度张量异常(a)及其DEXP变换结果(b)Fig.2 Gravity gradient tensor components of the rectangular prism model (a) and its DEXP transform results (b)

由图1与图2可知，虽然重力异常以及重力异常梯度均能反映理论模型的中心深度，但是重力垂向梯度的DEXP变换结果更好。在重力梯度张量各分量的DEXP变换结果中，Wxx、Wyy和Wzz的极值点所在位置均可以描述模型的中心深度信息，而Wxz、Wyz、Wxy与Wthd的极值点主要反映模型的边界信息，Wthd极值点反映的场源边界信息更为准确。

在一定精度的数据基础及相同的数据处理方法上，如何减小其他影响因素对极值点深度成像结果的影响尤为重要。为了简化分析，笔者以重力垂向梯度的DEXP成像结果为例，分析数据的误差(噪声影响)、点距、计算范围和背景场对DEXP成像结果的影响。

2.1 DEXP变换的抗噪能力分析

令直立棱柱体模型的长a、宽b与高c分别为20、20与6 m，其中心埋深为15 m，剩余密度Δρ为1 g/cm3，首先计算得到模型的重力垂向梯度异常，然后分别向理论模型加入4组平均值为零但方差不同的高斯噪声进行对比研究，其方差分别为理论重力垂向梯度最大幅值的1%、2%、5%和10%，对应的重力垂向梯度异常及其极值点深度成像结果分别如图3与图4所示。其中，模拟计算测点的点距为2 m，向上延拓总高度为30 m，延拓步长为1 m。

由图4加入了不同噪声的Tzz分量DEXP变换成像结果可知，其深度估计结果均与理论Tzz异常DEXP变换结果一致。这说明DEXP变换对噪声具有一定的压制作用，对加入较高比例的高斯噪声，其DEXP变换结果仍较准确。这是由于该算法主要基于向上延拓，而向上延拓属于低通滤波，因此受数据噪声的干扰很小。

(a) 理论Tzz异常；(b)、(c)、(d)、(e)分别为加入1%、2%、5%、10%最大幅值的高斯噪声之后的Tzz异常 (a) Tzz component with 0% Gauss noise; (b) Tzz component with 1% Gauss noise; (c) Tzz component with 2% Gauss noise; (d) Tzz component with 5% Gauss noise; (e) Tzz component with 10% Gauss noise图3 加入不同高斯噪声的Tzz分量异常Fig.3 Tzz component anomaly of the rectangular prism model with different Gauss noise

图4 与图3对应的Tzz分量DEXP变换结果Fig.4 DEXP transform results of Tzz component corresponding to Fig. 3.

2.2 数据点距对DEXP变换成像的影响

为了研究数据点距对异常体DEXP变换成像的影响，分别设置了3组不同点距(1、2以及5 m)进行试验研究。令直立棱柱体模型长a、宽b与高c分别为30、20与10 m，其中心埋深为30 m，剩余密度Δρ为1 g/cm3，计算范围为-70～70 m，向上延拓步长为1 m，总延拓高度为70 m。根据DEXP变换理论，即可得模型Tzz分量的深度成像结果，如图5所示。

由成像结果可知，当计算范围保持不变时，点距越小，所估计的模型埋深与理论中心埋深之间的差异越小。随着点距越大，所估计的深度结果更加偏离模型的中心深度位置，而往模型的上边界方向移动，并且其成像的整体形态越往异常体中心靠拢。

2.3 数据范围对DEXP变换成像的影响

采用与图5相同的模型，为了研究数据范围对模型DEXP变换成像的影响，分别设置了5组不同计算范围(-50～50 m、-60～60 m、-70～70 m、-80～80 m和-100～100 m)进行试验研究。其中，数据点距为2 m，向上延拓步长为1 m，总延拓高度为70 m。根据DEXP变换理论，即可得模型重力垂向梯度的DEXP成像结果，如图6所示。

由成像结果可知：当点距保持不变时，若计算范围小于一定范围，由于变换过程中没有足够的场源信息，其成像形态过于向两侧分散，导致模型的DEXP变换成像结果在计算范围内不存在极值点；随着计算范围的增大，模型的DEXP变换结果异常形态开始向中心收敛，其极值点开始向模型中心位置移动，当到某个计算范围时，其极值点与模型中心位置重合，即DEXP变化所估计的深度与理论埋深一致、误差为零；若继续增大计算范围， DEXP变换结果对应极值点则会向场源的上边界移动，当增大到某一个计算范围之后，极值点位置不再发生移动，其位置收敛在场源的上边界附近。

图5 数据点距为1 m(a)、2 m(b)和5 m(c)情况下的Tzz分量DEXP变换成像结果Fig.5 DEXP transform results of Tzz component with data spaces of 1 m (a), 2 m (b) and 5 m (c)

图6 计算范围为-50～50 m(a)、-60～60 m(b)、-70～70 m(c)、-80～80 m(d)、-100～100 m(e)时Tzz分量DEXP变换结果Fig.6 DEXP transform results of Tzz component with data spatial ranges of -50～50 m (a), -60～60 m (b), -70～70 m (c), -80～80 m (d) and -100～100 m (e)

2.4 背景场对DEXP变换结果的影响

采用与图5相同的模型进行试验，分别在模型正演中加入不同大小的背景场(分别为-1E、-0.5E、-0.2E、0E、0.2E、0.5E)，以研究背景场对DEXP变换成像结果的影响。模拟数据的测点点距为2 m，向上延拓步长为1 m，总延拓高度为70 m。图7为加入不同背景场后模型Tzz分量DEXP变换结果成像图。由图7可知，背景场对DEXP变换结果具有较大的影响，当背景场的绝对值过大时，DEXP变换成像形态会发生很大的改变，甚至可能导致变换结果没有极值点。因此，在实际数据资料处理时，应该首先采用适当的位场分离技术，再对局部异常进行DEXP变换。

3 基于实际数据的影响因素分析

研究区位于路易斯安那州西南部，靠近德克萨斯州边界的Vinton盐丘区。区内为沉积岩相，地下发育盐丘构造。Vinton盐丘区全张量重力梯度(FTG)数据由Bell Geospace公司测得，采用南北向飞行测线，间隔为150 m，平均海拔75 m。测区范围为北纬30.07°～30.23°、西经93.53°～93.66°，总测线长度为1 087.5 km，测区面积为192.6 km2。图8a即为Vinton盐丘地区的重力异常数据。在进行DEXP变换时，Wp中的深度加权函数zNp/2包含了地质体几何形状的构造指数信息，因此在DEXP成像过程中，需要预先指定异常体的构造指数。Zhou等[12]通过计算得出了盐丘地区梯度数据的构造指数(Np)大致为2.2，故而可得重力异常的构造指数为1.2，将其应用到式(8)中，即可得到Vinton盐丘地区重力异常DEXP变换结果，如图8b所示。由结果可得，重力异常数据的DEXP变换结果受到背景场的干扰很大，其成像结果成层状并且没有极值点。同时Fedi也指出重力梯度数据比重力异常的DEXP变换受附近异常体的干扰小[9]，故本文又将DEXP变换成像法应用于实测重力梯度张量数据处理中。图9为全张量重力梯度异常分布(ΔTxx、ΔTxy、ΔTxz、ΔTyy、ΔTyz和ΔTzz)。根据DEXP变换成像理论，即可以得到Vinton盐丘区全张量重力梯度数据的DEXP变换结果成像，如图10所示。图中DEXP变换的总水平导数W_thd对应的极值点位置可以大致描述盐丘的边界位置，Wxx、Wyy和Wzz变换结果的极值点位置可以反映盐丘的中心埋藏深度，极值点的平均位置约为335 m，故减去飞行高度计算出的平均深度约为260 m，其结果与Zhou等人[12]计算结果一致。

图7 加入-1E(a)、-0.5E(b)、-0.2E(c)、0E(d)、0.2E(e)和0.5E(f)的背景场之后的Tzz分量DEXP变换结果Fig.7 DEXP transform results of Tzz component with background fields of -1E(a)、-0.5E(b)、-0.2E(c)、0E(d)、0.2E(e)和0.5E(f)

图8 Vinton盐丘地区重力异常(a)及其DEXP变换结果(b)Fig.8 Gravity anomaly (a) in Vinton Salt Dome region and its DEXP transform results (b)

图9 Vinton盐丘地区的全张量重力梯度异常数据Fig.9 The full gravity gradient tensor fields in Vinton Salt Dome region

图10 Vinton盐丘地区全张量重力梯度数据的DEXP变换结果Fig.10 DEXP transform results of the full gravity gradient tensor field data in Vinton Salt Dome region

但是，根据前文的影响因素分析，结合图10中Wxx、Wyy和Wzz的成像形态，可以发现计算结果存在背景场的影响。因此，为了提高深度估计的精度，本文将ΔTxx、ΔTyy和ΔTzz分量进行匹配滤波分场处理，再将分场得到的局部场(如图11所示)进行DEXP变换成像，成像结果如图12所示。由变换结果计算得到的平均深度约为330 m，与前人计算的深度估算结果200～400 m一致[17-25]。

根据式(4)可以首先计算ΔTzz分量的尺度函数τp，再由Np=-2τp(z=-z0)即可估算得出Vinton盐丘区重力垂向梯度数据的构造指数Np。若不对ΔTzz分量进行分场处理，则根据DEXP变换极值点估计的构造指数为2.377 4，如图13a所示。图13b为利用ΔTzz分量局部场估计的构造指数结果，故可计算Vinton盐丘区重力垂向梯度数据的构造指数约为2.198，结果与Zhou等[12]估算结果一致，从而可得盐丘的构造指数为1.198，可知盐丘的几何形态比较接近于柱体状。

图11 Vinton盐丘地区ΔTxx、ΔTyy和ΔTzz分量分场之后的局部场Fig.11 The residual fields of ΔTxx, ΔTyy and ΔTzz components after anomaly separation in Vinton Salt Dome region

图12 Vinton盐丘地区ΔTxx、ΔTyy和ΔTzz分量的局部场DEXP变换成像结果Fig.12 DEXP transform results of the residual fields of ΔTxx, ΔTyy and ΔTzz components after anomaly separation in Vinton Salt Dome region

图13 构造指数Np估计结果Fig.13 The results of the corresponding structural index Np

4 结论

笔者基于DEXP变换理论，将其应用于重力异常以及重力梯度张量极值点深度成像中，结果表明重力梯度张量的成像结果较好，并且重力梯度张量极值点深度成像不仅能够反映地下异常体的中心埋深，也能大致描述其边界位置。

同时也设置了多组直立棱柱体模型进行试验研究，以讨论不同因素(数据误差、点距、计算范围和背景场)对极值点深度成像结果的影响。结果表明：极值点深度成像法具有较好的抗噪能力，具有计算稳定性和准确性特点；当计算范围保持不变时，点距越小，所估计的模型埋深与理论中心埋深之间的差异越小；当点距保持不变时，计算范围由小变大，从不存在极值点到极值点开始向模型上边界方向移动直至模型上边界附近的某一深度值后不再变化；背景场对DEXP变换结果具有较大的影响，当其绝对值过大时，DEXP变换成像形态会发生很大的改变，甚至可能导致变换结果没有极值点。因此，在实际数据资料处理时，应该综合考虑相关因素对DEXP变换结果的影响，并且进行相关预处理(尤其是合适的局部异常提取)以提高深度估计结果的精度。

将DEXP变换应用于Vinton盐丘地区的全张量重力梯度数据成像中，首先对ΔTxx、ΔTyy和ΔTzz分量进行了分场处理，再将分场得到的局部场进行DEXP变换成像，得到了盐丘较为精确的中心埋深与构造指数。鉴于DEXP成像方法具有抗噪性、准确性以及便捷性，而且可以同时得到场源埋深与构造指数，也可以将其成像结果作为自约束信息引入三维物性反演之中，因而具有较好的应用前景，值得推广。

致谢:感谢Bell Geospace公司提供的全张量重力梯度数据！