基于层源位场的重力及其梯度数据联合相关成像

林涛，曾昭发，于平，周帅，焦健

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

0 引言

重力勘探方法即根据地质体的密度差异测量其在地表引起的重力异常，解释地质体场源几何参数和物性参数的方法。重力勘探被广泛应用于区域地质调查、矿区普查、工程地质及水文勘察等领域[1-2]。随着重力梯度仪的不断发展，航空重力梯度测量技术由于其测量分辨率高、测量效率高等优势，在油气资源勘查、矿产资源勘探等方面逐渐发挥更大的作用。综合利用重力及其梯度数据可对地质体的场源参数反演提供更多的信息，提高现有解释方法的精度[3-4]。常规重力及其梯度异常反演是基于反演理论，设置目标函数，使用非线性或线性优化算法，使目标函数达到最低的方法[5-10]。但重力及其梯度反演数学上呈病态，多解性强，需要较多的约束条件，虽然能够对地下结构进行定量刻画，却很难在实际应用中取得良好效果[11]。

Patella[12]基于归一化互相关系数对自然电位测量数据进行地质目标体空间分布的快速成像。Mauriello和Patella[13-14]分别将这种方法推广应用于大地电磁和重力测量领域。此后，相关成像方法被广泛应用于地球物理数据处理中[15-18]。郭良辉等[19]基于重力梯度数据建立了相关成像公式，并采用异常分离方法提高了成像精度。陈召曦等[20]、侯振隆[21]实现了基于GPU并行的重力及其梯度相关成像，使得算法的运算效率大为提高。近年来，相关成像与其他地球物理数据和处理方法的结合成为趋势[22-26]。重力及其梯度相关成像技术主要用于定位目标体的中心埋深，由于完全忽略先验信息，相关成像方法的纵向分辨率较差，一般用于为精细反演提供初始模型[27]。

层源位场是通过异常分离方法得到的地下某一深度层产生的位场。程方道等[28]首先提出了可以分离出重力层源位场的插值-切割分离法。汪炳柱等[29]使用这种方法，将胜利油田某地区磁异常划分为浅、中、深3层的层源位场，并经过化极垂向一次导数和垂向二次导数的验证，证明该方法得到的浅部异常的可信度高。余海龙[30]首次将分离得到的层源位场使用积分-迭代法向下延拓至地层顶层，并将每一层分别进行反演，提高了反演的速度和精度。

笔者在前人工作的基础上，提出基于层源位场的重力及其梯度数据联合相关成像的方法，将重力及其梯度数据进行异常分离，并采用向下延拓技术得到某一深度地层引起的重力及其梯度异常，再对每一层的数据分别进行联合相关成像，通过合成数据和实测Humble盐丘重力数据的试验，验证了该方法能有效提升纵向分辨率。

1 方法原理

1.1 重力及其梯度联合三维相关成像原理

传统重力数据相关成像方法的基本原理是：对地下空间采用均匀网格进行剖分，对每个剖分的网格节点采用正演计算公式计算其引起的重力异常，将计算的重力异常与观测异常进行归一化互相关系数的计算，基于相关系数三维成像结果刻画地下场源的空间分布特征。

假设剖分网格的中心坐标为(xq,yq,zq)，体积为vq，剩余密度为σq的某点q在一个测点(x,y,z)处的重力异常Δgq(x,y,z)可表示为：

(1)

式中：G为万有引力常数，G=6.67×10-11m3/(kg·s2)。

计算第q网格引起的重力异常与实测重力异常的归一化互相关系数[19]：

(2)

式中：(xi,yi,zi)为第i个观测点的坐标，Δg(xi,yi,zi)为第i个观测点的重力异常，Δgq(xi,yi,zi)为q网格在第i个观测点位置引起的重力异常。

将式(1)代入式(2)中，有：

(3)

其中，

(4)

式中：Bq为地下q点的质量对于观测点(xi,yi,zi)的基函数，Cq取值范围在-1和1之间。相关系数Cq代表了观测的重力异常与某一地下网格节点的质量基函数的相关程度。

郭良辉等[19]将互相关系数推广到重力梯度，得到重力梯度的相关系数计算方法：

(5)

其中，

(6)

(7)

根据该理论，可以将重力梯度不同分量相关系数进行联合成像。

1.2 重力异常分离原理

重力及重力梯度测量在水平方向的分辨率较高，为提高相关成像的纵向分辨率，笔者引入徐世浙等[32]使用的插值-切割分离方法进行不同深度重力场的分离，该方法计算速度快，而且能将异常近似分离为不同地层深度引起的异常，其基本思路如下。

设位场异常Z(x,y)由区域场R(x,y)和局部场L(x,y)组成，即：

Z(x,y)=R(x,y)+L(x,y)

(8)

令A(x,y)代表与点(x,y)相距r的4个测点测得异常的平均值：

Z(x,y+r)+Z(x,y-r)]

(9)

式中：r为切割半径。

由式(9)可以得出区域场是Z(x,y)和A(x,y)的加权平均：

R(x,y)=aA(x,y)+bZ(x,y)

(10)

式中：a和b是加权系数，且a+b=1。

令：

(11)

(12)

ΔBx=Z(x+r,y)-Z(x-r,y)

(13)

ΔBy=Z(x,y+r)-Z(x,y-r)

(14)

(15)

(16)

e=c+d(0≤e≤2)

(17)

取加权系数：

(18)

由此得到一次切割场R1(x,y)。对R1(x,y)重复以上步骤，不断迭代下去，有：

(19)

于是有：

(20)

式中：R(x,y)为切割半径为r的区域场。由式(8)可以得到最终局部场：

L(x,y)=Z(x,y)-R(x,y)

(21)

根据该方法在实际应用中的经验，可以近似地认为，切割半径为r时得到的局部场是深度r以上的地层产生的异常场。在此基础上，余海龙[30]提出了层切割技术，即用不同的两个切割半径对原始异常进行切割，用切割结果的局部场做减法，其结果便可近似地表征为地下某具有上下底面的层位在地表引起的层源异常。利用这个技术，将重力数据划分为不同深度的地层引起的多组重力异常，从而在数据处理过程中提高纵向分辨率。

1.3 重力异常向下延拓法原理

用切割法分离异常之后所得的异常场是地面测得的地下不同深度层的异常场。若直接采用分层相关成像的方法，深部底层距测点仍较远，会产生不必要的误差。所以在切割分离之后，采用向下延拓的方法，将测量面下延至切割层的顶层，此时再采用相关成像的方法，可以更加精确地分辨层内密度特征。

本文采用的向下延拓方法为徐世浙[33]提出的迭代法，其原理如下：

假设存在两个平面ΓA和ΓB，用u(x,y,0)代表平面的位，用U(kx,ky,0)=F[u(x,y,0)]表示u(x,y,0)的Fourier变换，则z=h(ΓB)平面的位可表示为：

(22)

式中：Γ-1是Fourier反变换。

(23)

式中：s是步长，一般为1。

重复进行以上过程，可以得到迭代公式：

(24)

一般迭代20～50次即可。这种方法将不稳定的向下延拓问题转换为稳定的向上延拓问题，使得延拓的结果更加稳定，并且能延拓更大的深度。经过异常分离和向下延拓后，可以将某一深度范围内产生的重力及重力梯度异常进行联合相关成像，得到更为精确和易于解释的结果。

2 模型验证

2.1 传统联合相关成像

设计由11个不同大小、不同密度的棱柱体组成的地质模型，其具体参数如表1所示。选择y=8 000 m测线为研究对象，地质体大致分为3层，A1、A2层埋深为1 500～3 000 m；B1、B2、B3层埋深为400～800 m；C1、C2层为100～200 m。其俯瞰图和y=8 000 m截面图见图1。

表1 复合模型的几何参数和剩余密度

对组合体的Δg、Gxz、Gyz、Gzz分别进行网格为20 000 m×20 000 m×3 000 m，步长为m的三维相关成像，选取y=8 000 m测线作剖面图。为使图像具有一致性，文中的相关系数结果均经过归一化处理。由图2可以看出，关于重力异常的相关成像可以突出底部较大地质体的异常A1、A2，而关于重力梯度的相关成像方法对中等深度的地质体B1、B2、B3有异常反应，但它们对浅层地质体C1、C2的信号都不明显。

对Gxz、Gyz、Gzz三个分量进行联合相关成像，得到图3。图像明显收敛，B3产生的异常虽然较小，但精确地反映了其位置，仍存在以下问题：

①A1的质心比A2要高，推测是测线在A1的中间部分而在A2的边缘部分，使得A2的异常值降低更大。

②受A1、A2影响比较大的B1、B2的异常消失了。

图2 y=8 000 m测线单分量相关成像Fig.2 Single-component probability tomography of y=8 000 m line

图3 y=8 000 m测线联合相关成像Fig.3 Joint probability tomography of y=8 000 m line

可见联合相关成像能够突出主要异常，但会忽略较小地质体的异常。为突出浅部较小地质体的异常，可以先进行异常分离，将分离后的数据再进行相关成像。图4是进行3×3均值滤波分离异常之后的多分量相关成像结果，中层B1、B2、B3和浅层C1、C2的地质体异常更加明显，虽然能根据最大值大致判断其质心的横纵向位置关系，但对中浅层地质体的深度分布却很模糊，纵向分辨率并不高。

此时再对Gxz、Gyz、Gzz三个分量进行联合成像，所得到的结果如图5所示，虽然浅层地质体C1、C2的异常值减弱了，但通过最大值位置的判断，可以清晰地分辨出B1、B2、B3以及C1、C2的质心深度的相对关系，并能推断其形态，但仍存在如下问题：

①对C1、C2的质心深度判断并不准确，推测是由于测点距为200 m，而C层位的分布深度在100～200 m导致精度不够。

图4 y=8 000 m测线基于均值分离的单分量相关成像Fig.4 Single-component probability tomography of y=8 000 m line based on mean separation

图5 y=8 000 m测线基于均值分离的联合相关成像Fig.5 Joint probability tomography of y=8 000 m line based on mean separation

②C2的位置比实际向下，推测是B2和A2的干扰所致。

尽管如此，联合相关成像的精度还是显著提高了。

根据联合相关成像在组合体模型中的实验，可以看出，相对于传统相关成像方法，联合相关成像的方法可以压制干扰，使得成像结果更加收敛，更加清楚地分辨地质体的质心深度。同时，该方法也会压制一些实际存在的较小异常。

2.2 基于层源位场的联合相关成像

对上述组合体模型采用切割分离的方法，将位场异常分离为不同深度的地质体引起的异常，再向下延拓至地层顶面后进行相关成像。研究切割厚度较大时的图像特征，将地层切割为两层，0～1 000 m为一层，1 000～2 000 m为一层，分别对其进行传统相关成像和联合相关成像，得到图6和图7。

可以看出，在切割范围较大的情况下，传统成像方法受干扰较大，虽然其异常范围可以与模型相对应，但假影较多，难以解释。而联合相关成像方法的结果清晰。图6中可以看到B层的埋深在600 m±，且其存在于0～1 000 m地层的下部，A1、A2的埋深在1 500 m和2 000 m，分布于1 000～2 000 m地层的下部，与模型较为符合。

图6 0～1 000 m地层y=8 000 m测线相关成像Fig.6 Probability tomography of y=8 000 m line within 0～1 000 m strata

图7 1 000～2 000 m地层y=8 000 m测线相关成像Fig.7 Probability tomography of y=8 000 m line within 1 000～2 000 m strata

图8 y=8 000 m测线多地层联合相关成像Fig.8 Multi-strata joint probability tomography of y=8 000 m line

在此基础上，为得到更精确的层源信息，将地层进一步精细切割，并对切割的结果进行联合相关成像(图8)。由此可推断出地质体的顶底面埋深：

①0～200 m地层中，只有C层位引起异常。由于顶界面埋深小于测点距，无法得出C层地质体在该地层的分布位置。记C层埋深为0～200 m。

②200～400 m地层中，异常主要由B层引起，可以明显看出该层异常存在于地层下部，可记其顶面埋深为300 m。

③400～800 m地层中，异常主要由B层地质体引起。

④800～1 200 m地层中，没有主要异常，故可将B层地质体的底面埋深记为800 m。

⑤1 200～1 600 m的地层中，A1的异常出现，记A1的顶面埋深为1 200 m。

⑥1 600～2 400 m的地层中，A2的异常出现，记其顶面埋深为1 600 m。

由于切割成像和向下延拓法均有深度限制，为保证切割数据的精确性，只将地层切割至2 400 m的深度。

由此可以推断出地质体顶底面埋深：

浅层：0～200 m(实际100～200 m)。

中层：300～800 m(实际400～800 m)。

深层(顶部埋深)：A1为1 200 m；A2为1 600 m(实际均为1 500 m)。

由于在数据处理过程中明确了层位信息，基于层源位场的重力及其梯度联合成像方法在模型实验中，精度高于传统相关成像方法。由于异常分离和延拓的过程产生噪音，单独数据的相关反演方法并不能在层源位场数据取得较好效果。

3 Humble盐丘实测重力数据应用

将本文方法应用到美国Houston地区Humble盐丘数据中。图9是盐丘的重力及重力梯度异常平面图，采样间隔为100 m。为研究盐丘的纵向分布特征，对该测区的重力梯度数据分别按1 000～3 000 m，3 000～5 000 m，5 000～7 000 m进行异常分离与延拓，将3个地层数据分别进行步长为100×100×100的多分量联合相关成像，并画出y=6 000 m截面的剖面图(图10)。可以看出，该盐丘的中心埋深为4 km ，顶面在3 km。表2为本方法得到的盐丘埋深与其他方法的对比。可以看出，该方法得到的埋深与其他方法相近，说明重力梯度多分量联合相关成像的反演结果是可靠的。

4 结论

(1)相较于单分量相关成像，联合相关成像可以有效压制假影，使成像结果更为聚焦，同时也会压制一些真实影像。

表2 重力及其梯度联合相关成像与其他方法的对比

图9 Humble盐丘重力及重力梯度异常平面图Fig.9 Gravity and gravity gradiometry anomalies in Humble salt dome

图10 Humble盐丘测区y=6 000 m测线多地层联合相关成像剖面图Fig.10 Multi-strata joint probability tomography of y=6 000 m line within Humble salt dome survey area

(2)单分量相关成像在处理层源位场数据时会产生大量假影，联合相关成像依然聚焦。

(3)对层源位场数据进行联合相关成像，可以在提升纵向分辨率的同时，抑制假影的产生，得到更准确的结果。