三维重力、重力梯度和大地电磁数据的平滑聚焦结构约束联合反演

何浩源， 李桐林， 张镕哲*， 朱威

1 吉林大学地球探测科学与技术学院， 长春 130012 2 中国地质科学院地球物理地球化学勘查研究所自然资源部地球物理电磁法探测技术重点实验室， 廊坊 065000

0 引言

重力、重力梯度和大地电磁方法都是地球物理勘探领域中的重要技术手段，无论是在理论还是实际应用方面均已得到国内外学者的广泛关注与研究.这些方法出于原理或实践上的差异，在分辨能力上表现出不同特点.例如，重力和重力梯度数据虽具有较高的横向分辨率，却因信号随深度增加的迅速衰减，使其纵向分辨率较低.重力梯度数据作为重力数据的一阶导数，含有更多高频信息，但低频信息量少于重力数据.因此，重力梯度数据相比重力数据对浅层区域的分辨能力更强，而对深层区域的分辨能力较弱.大地电磁数据则因包含大量频率信息而更具测深能力，但通常大地电磁法的野外实测间距较大，测点较稀疏，使该方法的横向分辨能力相对较弱.单一物探数据只能从单一物性角度对地下介质作出评价，考虑到重力、重力梯度和大地电磁法各自的优缺点，单一物探方法的解释结果经常有失准确性.但如果采用多种数据从多个角度对同一区域的地下介质进行研究，则可以更加全面地预测地质信息.因此，多种地球物理数据综合解释已成为地球物理勘探发展的重要趋势.而联合反演作为多种地球物理数据综合解释的有效方法之一，是降低反演多解性和改善单一方法局限性问题的重要手段.

20世纪70年代，Vozoff和Jupp(1975)首次提出地球物理联合反演的思想，并将其应用于直流电阻率法和大地电磁法的联合解释.随后的45年中，联合反演技术发展迅速，从一维联合反演发展到三维联合，从同一物性联合反演发展到不同物性联合，从物性关系联合反演发展到结构约束联合.根据物性参数间构建联合关系方式的不同，可以将联合反演分为两大类.第一类为岩石物性关系联合反演，该类方法需要提前构建不同物性参数间的经验关系式，通过经验关系式进行物性参数间的转换(Nielsen and Jacobsen，2000；Afnimar et al.，2002；Colombo and De Stefano，2007).然而不同地区的物性关系存在差异，且复杂地区的经验关系式难以建立，因此该方法具有一定局限性.第二类为结构相似性联合反演，该类方法假设同一区域内，不同物性参数具有相同或部分相同的地下结构.目前，由Gallardo和Meju(2003)提出的基于交叉梯度约束的联合反演已成为最受欢迎的结构相似性联合反演方法之一.该方法通过约束不同物性参数的梯度方向一致，来寻找结构相似的模型，既不依赖经验关系式，也不需要对初始模型做任何假设，是一种有效的地球物理综合解释手段(Gallardo and Meju，2004；彭淼等，2013；Zhang et al.，2019).

重力勘探和大地电磁勘探因具有覆盖面广、效率高和费用低等特点而获得广泛应用，两种方法间的交叉梯度联合反演也受到了地球物理学者们的关注和研究.其中，李桐林等(2016)实现了部分区域约束的重力、磁法和大地电磁的二维交叉梯度联合反演；Pak等(2017)开展了二维数据空间下的重力、磁法和大地电磁的交叉梯度联合反演；张镕哲等(2019)将交叉梯度约束应用于重力、磁法、大地电磁和地震初至波走时四种方法的二维联合反演中；闫政文等(2020)实现了三维重力、磁法和大地电磁法的交叉梯度联合反演.上述研究成果均表明，基于交叉梯度约束的联合反演对异常的分辨能力相比单独反演更强，但也存在一定的局限性:1)现有的重力和大地电磁联合反演更多集中在二维研究，三维研究较少.2)由于重力数据纵向分辨能力差，通过交叉梯度结构约束传递的不准确纵向密度信息通常会误导电阻率模型，导致大地电磁反演结果偏离真实模型.3)平滑正则化方法在联合反演中的应用会使结果表现出异常发散的现象，影响联合反演的分辨率.4)传统的连续式联合反演策略，常因交叉梯度的结构相似约束使反演模型向真实模型迭代的趋势受到限制，导致反演模型难以收敛.

为解决上述问题，本文提出了一种基于交叉梯度约束的三维重力、重力梯度和大地电磁数据联合反演方法.首先，引入重力梯度数据能够为交叉梯度联合提供更为可靠的密度模型，有利于提高联合反演的准确性.其次，为了提高电阻率模型的自身分辨率，本文提出了一种平滑聚焦的正则化方法.该方法既能解决平滑反演结果异常体发散的问题，又能避免聚焦反演伴生的近地表假异常问题.同时，为聚焦和平滑聚焦正则化加入了聚焦强度阈值方法，以限制高物性异常区域的聚焦优势，避免过度聚焦问题.此外，借鉴Um等(2014)交替进行结构约束和数据约束的思想，本文提出了一套间断式联合的交叉梯度联合反演流程.整个反演流程仅经历一次收敛过程，且避免了某一方法对联合反演结果绝对主导的问题.在保证约束不同物性模型间结构相似的同时，又能一定程度上减小结构相似性约束对反演模型向真实模型迭代带来的限制.最后，通过理论模型测试，验证了本文所提出方法的有效性和准确性.

1 三维重力、重力梯度和大地电磁正反演理论

本文中，重力、重力梯度及大地电磁的正反演均在笛卡尔坐标系下进行计算，地下模型被划分为若干块物性均匀的规则长方体.

1.1 三维重力、重力梯度和大地电磁正演

根据Nagy(1966)，地下长方体在测点处产生的重力Vz及重力梯度Vzz的正演计算公式如下：

(1)

(2)

图1 长方体与测点坐标关系示意图Fig.1 Coordinate relationship between cuboid and measuring point

其中，(ξ1，ξ2)、(η1，η2)和(ζ1，ζ2)依次为长方体沿x、y、z轴的坐标范围，(x,y,z)为观测点P的坐标，空间关系如图1所示.G为万有引力常数，ρ为长方体密度，r为长方体中心与观测点P的距离.

假设地下模型被划分成M=Nx×Ny×Nz块小长方体，则有N=Nx×Ny个测点位于各长方体中心沿Z方向在水平地表的投影处.重力或重力梯度的观测数据d与模型物性m的正演关系可写成如下矩阵形式：

d=Jm,

(3)

式中，J是N×M的灵敏度矩阵.

大地电磁正演及计算灵敏度矩阵J的部分主要参考三维大地电磁反演程序WSINV3D(Siripunvaraporn et al.，2005).该程序采用FDE(Finite Difference Equations)方案(Yee，1966)进行有限差分交错网格划分，即电场分量定于网格棱边的中点，方向沿着棱边，磁场分量定于网格表面的中心，方向沿着面法向分量.该方案由麦克斯韦方程组出发，消去磁场，得到电场的波动方程.对波动方程离散化后使用QMR算法求解电场，再由电场计算磁场，最终计算出观测点处的四个阻抗分量Zxx、Zxy、Zyx和Zyy.

1.2 三维重力、重力梯度和大地电磁单独反演

本文中，重力、重力梯度及大地电磁均采用了模型空间下的正则化反演，其目标函数如下：

(4)

对式(4)求导，令求导结果等于0，可以得到模型空间下的反演模型迭代表达式：

(5)

在没有明确参考模型的情况下，每一次迭代时将当前模型mk设置为参考模型mref.使用共轭梯度法求解线性方程组即可得到每次迭代的模型更新量.

2 三维重力、重力梯度和大地电磁数据交叉梯度约束联合反演理论

重力、重力梯度和大地电磁数据联合反演的主要联合关系为：重力与重力梯度直接联立方程组，联合求解剩余密度模型(本文中为了方便理清关系，将此联合称为密度单独反演)；剩余密度模型与大地电磁对应的电阻率模型进行交叉梯度约束联合.联合关系框架如图2所示.

图2 重力、重力梯度和大地电磁联合关系示意图Fig.2 Joint relation diagram of gravity， gravity gradient and magnetotelluric

2.1 重力和重力梯度数据直接联合反演

重力和重力梯度的联合反演(即密度单独反演)，为密度模型不同观测数据之间的联合反演，直接联立它们的方程组求解即可.综合考虑反演的计算量和结果的准确性，本文只采用Vz和Vzz分量进行联合反演(秦朋波，2016).令：

(6)

(7)

(8)

将式(6)—(8)代入式(5)中，得到重力和重力梯度联合反演的模型迭代表达式：

×[dgra-Jgramk],

(9)

A·Δm=b，

(10)

(11)

求解出WmΔm，进而求得模型更新量.

2.2 交叉梯度约束与联合网格剖分

根据Gallardo和Meju(2003)定义的三维交叉梯度函数：

(12)

其中，mρ为剩余密度，mσ为对数电阻率.交叉梯度的三个分量分别为

(13)

(14)

(15)

采用向前差分的方式对式(13)—(15)进行离散(Fregoso and Gallardo，2009)：

-mρ d(mσ f-mσ c)],

(16)

-mρ r(mσ d-mσ c)],

(17)

-mρ f(mσ r-mσ c)],

(18)

其中，mc为当前块物性，mr为X方向下一块的物性，mf为Y方向下一块的物性，md为Z方向下一块的物性.

将交叉梯度项加入式(4)，得到交叉梯度联合反演的目标函数：

+βρtT(mρ,mσ)t(mρ,mσ),

(19)

+βσtT(mρ,mσ)t(mρ,mσ),

(20)

其中，β为交叉梯度约束项权重系数.令式(19)和式(20)的导数等于0，得到模型迭代表达式：

(21)

(22)

式(21)和式(22)的等号右端项对应第k次迭代的数据.B为交叉梯度的离散导数，其求取方式可参考Fregoso和Gallardo(2009)的文章，本文限于篇幅不作展开.

图3 联合网格划分示意图Fig.3 Diagram of joint grid generation

大地电磁反演通常要求纵向网格间距随着深度增加而逐渐增加，而重力和重力梯度对网格划分的要求则较为宽松.因此，如图3所示，联合网格的公共部分(中心区)采用X、Y方向均匀剖分，Z方向网格距随深度增加而增加的剖分方式.为满足大地电磁边界条件，还需要在中心区外设置扩边区.扩边区内，越远离中心区的网格，X、Y、Z方向的间距越大.

考虑到计算交叉梯度要求两种物性模型网格一致，在为大地电磁联合反演计算交叉梯度之前，需要先将剩余密度模型扩边，扩边区剩余密度按照背景值设定.而在为重力和重力梯度联合反演计算交叉梯度之前，则需要先将电阻率模型缩边，只保留中心区.

2.3 平滑聚焦的模型正则化方法

平滑约束和聚焦约束的有效性无论在理论还是应用上均已得到充分证明(Constable et al.，1987；Zhdanov，2009；刘小军，2007).本文结合两种传统约束方法，为大地电磁反演提出了一种平滑聚焦的模型正则化方法，并将其应用于联合反演中.该正则化方法既能避免平滑反演存在的异常体发散问题，又能解决聚焦反演伴生的近地表假异常问题，具体方法如下.

平滑聚焦正则化的平滑部分如下：

(23)

(24)

(25)

(26)

其中，Rx、Ry和Rz分别用于约束模型X、Y、Z方向的一阶导数.对于Rx，若当前块非X方向的最后一块，则用X方向的下一块减当前块作为一阶导数；若当前块是X方向的最后一块，则用当前块减X方向的上一块作为一阶导数.Ry和Rz同理.

平滑聚焦正则化的聚焦部分如下(Zhdanov，2015；孙思源，2019)：

(27)

Re=diag(1/me),

(28)

其中，常数e为聚焦因子，e越小，迭代模型越易聚焦.在没有明确参考模型的情况下，通常令此处的mref等于0，即均匀半空间.me是M×1的矩阵，Re为模型聚焦项，是以1/me为主对角线的M×M矩阵.

但该聚焦方式经常存在过度聚焦的问题.由式(27)和式(28)，当前模型物性参数与均匀半空间相差越大的块，越容易聚焦，越容易在下次迭代中产生更大的物性差异.这种聚焦优势经常会随着模型迭代逐渐增大，最终导致反演结果的物性异常聚集在比真实模型边界更小的范围内.仅根据实际资料在迭代过程中对物性参数值的上下限进行约束尚不足以充分解决此问题.因此，本文对模型聚焦项进行了改进，通过设置阈值et对聚焦强度进行限制.由式(27)得到me后，将me中所有高于et的值降至et，再将新得到的me代入后续计算.根据模型测试经验，et通常可取原me最大值的1/2至2/3.

模型正则化的本质是约束各模型块物性参数间的数学关系，平滑约束令相邻模型块间的物性参数变化极小，而聚焦约束对物性参数与均匀半空间相差大的块约束小，对相差小的块约束大.同时利用两种约束关系，丰富模型限制条件，有利于减少反演的多解性.于是为平滑聚焦正则化方法设计如下的模型协方差矩阵：

(29)

在迭代初期，由于物性起伏较小，Re的聚焦作用并不明显，平滑聚焦约束整体偏向平滑.随着迭代继续，物性起伏愈发明显，Re的聚焦作用也愈发突出，平滑聚焦约束逐渐偏向聚焦.考虑到如果任一物性在某处的梯度值为0，则该处两种物性的交叉梯度也为0.那么在交叉梯度联合反演的迭代初期，通过偏向平滑的正则化使物性异常发散，能够得到更大范围的有效交叉梯度约束，有利于提高两种物性模型结构变化的同步性.

(30)

(31)

其中，Rs为基于灵敏度矩阵的深度加权项(Mehanee et al.，1998)，是一个M×M的对角阵.

2.4 重力、重力梯度和大地电磁数据联合反演流程

图4 重力、重力梯度和大地电磁联合反演流程图Fig.4 Flow chart of gravity，gravity gradient and magnetotelluric joint inversion

如图4所示，在交叉梯度联合之前，先对两种物性的初始模型进行几次单独反演迭代，使物性模型初显结构.考虑到交叉梯度联合反演具有约束两种物性结构相似的特性，而迭代过程中常常存在两种物性结构已近似一致，但均与真实模型存在偏差的情况.在传统的连续式联合反演中，由于两种物性间的结构约束是分步、交替且连续进行的，反而使反演模型向真实模型迭代的趋势受到限制，导致联合反演收敛困难，反演结果发散.

Um等(2014)为了克服地震和电磁联合反演收敛困难的问题，提出了一种交替进行结构约束和数据约束的反演策略.其反演流程中，先分别令地震和电磁数据单独反演至收敛，然后固定所得速度模型对所得电阻率模型进行仅含模型约束项和交叉梯度约束项的纯结构约束反演，再将结构约束后的电阻率结果作为新的初始模型进行电磁单独反演直至收敛.一次纯结构约束反演与一次电磁单独反演被合称为一次电磁精化过程(EM refinement process).重复数次电磁精化过程直至连续两次精化的电阻率结果基本相近，再固定电阻率模型对地震施行类似的精化过程，如此往复，最终得到联合反演结果.该策略虽有效改善了联合反演的收敛性，但也明显存在结构约束过强的问题.尤其是在纯结构约束反演中缺少数据项约束，而多次电磁精化后才反轮到电阻率结构约束地震反演，导致联合反演的结果几乎完全取决于地震反演所得结构.此外，由于每次精化都意味着一个完整的反演收敛过程，反复地精化使联合反演流程过于冗长.

借鉴Um等交替进行结构约束和数据约束的思路，本文设计了一套间断式联合的反演流程，即根据当前迭代次数间隔式地令交叉梯度约束项权重系数β=0，相当于让联合反演和单独反演交替进行.根据参与交叉梯度联合的物性方法在可信度上的差异，又可为流程分出两种间断方案：若两种物性方法可信度相近，则采用如图5a的方案，在奇数次迭代中对方法一进行交叉梯度约束，令方法二单独反演，在偶数次迭代中对方法二进行交叉梯度约束，令方法一单独反演，整体上两种方法的物性结构互相影响；若方法二的可信度明显高于方法一，则采用如图5b的方案，在奇数次迭代中对两种方法都进行交叉梯度约束，在偶数次迭代中都不约束，整体上偏向于以方法二的物性结构为主导.对于第二种方案，虽然迭代过程中总是先由方法二的物性结构约束方法一，但由于约束策略的变动精细到每次迭代，被约束的物性结构不会立即与约束方高度一致，使得方法一的物性结构仍对方法二的物性结构有一定的约束作用，避免了某一方法因结构约束过强而对联合反演结果绝对主导的问题.具体到本文的重力、重力梯度和大地电磁联合反演，可通过单独反演测试来选择合适的间断方案.若大地电磁单独反演结果呈现出明显的纵向电阻率结构，则选择偏向以大地电磁为主导的方案.若物性无明显纵向结构，则可以选择互相约束的方案，抑或为突出横向物性变化而选择偏向以重力和重力梯度为主导的方案.这种间断式联合的交叉梯度联合反演策略在整个反演流程中仅经历一次收敛过程，相比于现有的连续式联合，既能保证不同物性参数间的结构相似性约束，又能让迭代模型的响应向观测数据充分拟合.

2.5 联合反演参数取值方案

本节对模型约束项权重系数λ、交叉梯度约束项权重系数β和交叉梯度联合前单独反演次数的取值方案进行说明.λ和β的计算方法如下：

图5 间断式联合方法示意图 (a) 方案一； (b) 方案二.Fig.5 Diagram of discontinuous joint method (a) Scheme 1; (b) Scheme 2.

(32)

(33)

通过比值方法为系数λ和β赋值，能够保证数据约束项、模型约束项和交叉梯度约束项在数量级上相近.对于重力和重力梯度反演，由式(32)为λ赋初值，并在迭代过程中采用与张镕哲(2020)相同的阻尼因子方法对λ进行修正.而大地电磁反演的λ及所有反演的β则均采用随迭代次数线性减小λ0和β0的方法进行取值：

(34)

(35)

其中，k为当前迭代次数，Nit为最大迭代次数，a和b为常数.

3 模型试验

采用双长方体理论模型对提出的算法进行模型试验.如图6a和图6b，左侧长方体大小为2000 m×1500 m×2000 m，顶面埋深850 m，剩余密度为1 g·cm-3，电阻率为1000 Ωm；右侧长方体大小为2000 m×1500 m×1900 m，顶面埋深1500 m，剩余密度为1 g·cm-3，电阻率为10 Ωm；背景剩余密度为0，电阻率为100 Ωm.中心区采用22×22×14的三维剖分，X、Y方向为均匀剖分，网格间距为500 m.Z方向网格间距随着深度增加逐渐从200 m增加到2000 m.中心区外，在X、Y的正负方向及Z的正方向各扩边4个网格，间距依次为4000、8000、1600和32000 m，地下网格共30×30×18块.重力和重力梯度的测点数均为22×22个，分布于中心区地表各网格的中心处.大地电磁共8×8个测点，如图6c所示布设于地表.大地电磁采用9个频率观测，依次为100、10、5、3 Hz、2、1.5、1、0.5 Hz和0.25 Hz.为重力、重力梯度和大地电磁的合成数据分别加入3%随机噪声.

图6 双长方体模型设计示意图 (a) 剩余密度模型三维透视图； (b) 对数电阻率模型三维透视图； (c) 大地电磁测点分布图.Fig.6 Design diagram of double cuboid model (a) 3D perspective of residual density model； (b) 3D perspective of log resistivity model； (c) Distribution of magnetotelluric measuring points.

各组模型试验均使用均匀半空间作为初始模型，并保证充分迭代.使用数据拟合差RMS和模型拟合差RMSE对反演结果进行评价：

(36)

(37)

其中，dobs为观测数据，dinv为反演结果的正演响应，mtrue为真实模型，minv为反演结果模型，err为观测误差.RMS越小表示数据拟合程度更高，RMSE越小则表示模型还原程度更高.

模型试验3.1—3.3节中的所有二维模型切片均在x=5500 m处，为方便对比，使用黑色粗线框在反演结果切片中标记了异常体的真实边界位置.试验中提及的所有“密度”均指剩余密度，“电阻率”均指对数电阻率.

3.1 正则化方法对比试算

3.2 不同数据的交叉梯度联合反演对比试算

图7 (a) 真实电阻率模型切片； (b) 大地电磁平滑反演结果切片； (c) 大地电磁聚焦反演结果切片； (d) 大地电磁平滑聚焦反演结果切片Fig.7 (a) Real resistivity slice； (b) Slice of magnetotelluric smooth inversion； (c) Slice of magnetotelluric focusing inversion； (d) Slice of magnetotelluric smooth-focusing inversion

图8 不同正则化方法大地电磁单独反演RMS曲线Fig.8 RMS curve of magnetotelluric inversion by different regularization methods

表1 不同正则化方法大地电磁单独反演结果的 RMS和RMSE值Table 1 RMS and RMSE of magnetotelluric single inversion results by different regularization methods

图9a和图9b分别为真实密度模型切片和真实电阻率模型切片；图9c和图9d分别为重力单独反演结果切片和大地电磁单独反演结果切片；图9e和图9f分别为重力和大地电磁交叉梯度联合反演结果的密度切片和电阻率切片；图9g和图9h分别为重力、重力梯度和大地电磁交叉梯度联合反演结果的密度切片和电阻率切片.图10为上述反演迭代的RMS曲线，各组反演结果的RMS和RMSE如表2所示.由图9(a—f)发现，重力纵向分辨能力弱，其单独反演结果的纵向密度结构与真实模型相差较大，无法为交叉梯度联合反演提供可靠的密度模型.由于交叉梯度联合方法的结构相似性约束，不可靠的纵向密度模型影响了本来更为可靠的纵向电阻率模型，导致联合反演效果差.而如图9g和图9h，异常体的位置、几何形态和物性参数幅值明显更还原真实模型，可见重力梯度信息的引入极大地提高了密度模型的可靠性，从而改善了联合反演效果.图11对比了重力、重力梯度和大地电磁联合反演中，剩余密度和对数电阻率间交叉梯度的模的变化.图11a为单独反演阶段结束时的交叉梯度结果切片，图11b则为联合反演最终结果的交叉梯度切片.交叉梯度模的减小说明了交叉梯度联合能够使两种物性反演结果具有更高的结构一致性.

表2 不同数据交叉梯度联合反演结果的RMS和RMSE值

3.3 聚焦强度阈值和间断式联合算法试算

沿用上一节的反演参数，基于重力、重力梯度和大地电磁交叉梯度联合反演对聚焦强度阈值和间断式联合算法进行模型试算分析.图12a和图12b为无聚焦强度阈值的交叉梯度联合反演结果的密度模型切片和电阻率模型切片；图12c和图12d为连续式联合的交叉梯度联合反演结果的密度模型切片和电阻率模型切片；图12e和图12f(同图9g和图9h)为加入聚焦强度阈值且间断式联合的交叉梯度联合反演结果的密度切片和电阻率切片，其中聚焦强度阈值设置为原最大值的0.6倍.各组反演迭代的RMS曲线如图13所示，反演结果的RMS和RMSE如表3所示.

对比图12a和图12e，发现图12a的反演结果存在过度聚焦问题，高密度块聚集在比真实模型边界更小的范围内，而加入了聚焦强度阈值的图12e物性参数聚焦地更为均匀，更吻合真实模型.对比图12b和图12f也能得到同样的结论，即聚焦强度阈值能有效改善聚焦反演和平滑聚焦反演中存在的过度聚焦问题.如图12c和图12d，连续式联合的每次迭代都进行交叉梯度约束，当密度和电阻率的结构已近似一致时，两种交替迭代的物性模型反而受结构相似性约束的限制，难以进一步向真实模型收敛.而间断式联合的方法，如图12e和图12f，明显改善了这一问题.

图9 (a) 真实密度模型切片； (b) 真实电阻率模型切片； (c) 重力单独反演结果切片； (d) 大地电磁单独反演结果切片； (e) 重力和大地电磁联合反演结果密度切片； (f) 重力和大地电磁联合反演结果电阻率切片； (g) 重力、重力梯度和 大地电磁联合反演结果密度切片； (h) 重力、重力梯度和大地电磁联合反演结果电阻率切片Fig.9 (a) Real density slice； (b) Real resistivity slice； (c) Slice of gravity single inversion； (d) Slice of magnetotelluric single inversion； (e) Density slice of gravity and magnetotelluric joint inversion； (f) Resistivity slice of gravity and magnetotelluric joint inversion； (g) Density slice of gravity， gravity gradient and magnetotelluric joint inversion； (h) Resistivity slice of gravity， gravity gradient and magnetotelluric joint inversion

3.4 复杂理论模型试算

为进一步验证本文联合反演算法的有效性，本节对复杂理论模型进行试算.沿用与前一模型相同的网格划分、测点布设及大地电磁观测频率组，同样为重力、重力梯度和大地电磁的合成数据分别加入3%随机噪声.如图14所示，模型由深浅不一的两块低阻高密度长方体和一块高阻高密度长方体组成，另有一条高阻高密度台阶斜搭在其中一块长方体上.

图10 不同数据交叉梯度联合反演RMS曲线 (a) 重力数据RMS； (b) 重力梯度数据RMS； (c) 大地电磁数据RMS.Fig.10 RMS curve of cross-gradient joint inversion with different data (a) RMS of gravity data; (b) RMS of gravity gradient data; (c) RMS of magnetotelluric data.

图11 重力、重力梯度和大地电磁联合反演交叉梯度变化对比图 (a) 单独反演阶段结果交叉梯度模切片； (b) 联合反演最终结果交叉梯度模切片.Fig.11 Contrast of cross-gradient for gravity，gravity gradient and magnetotelluric joint inversion (a) Modulus slice of cross-gradient for results in single inversion phase； (b) Modulus slice of cross-gradient for joint inversion final result.

图12 (a) 无聚焦强度阈值的联合反演结果密度切片； (b) 无聚焦强度阈值的联合反演结果电阻率切片； (c) 连续式联合反演结果密度切片； (d) 连续式联合反演结果电阻率切片； (e) 加入聚焦强度阈值且间断式联合的反演结果密度切片； (f) 加入聚焦强度阈值且间断式联合的反演结果电阻率切片Fig.12 (a) Density slice of joint inversion without focusing intensity threshold； (b) Resistivity slice of joint inversion without focusing intensity threshold； (c) Density slice of continuous joint inversion； (d) Resistivity slice of continuous joint inversion； (e) Density slice of discontinuous joint inversion with focusing intensity threshold； (f) Resistivity slice of discontinuous joint inversion with focusing intensity threshold

图13 不同方式交叉梯度联合反演RMS曲线 (a) 重力数据RMS； (b) 重力梯度数据RMS； (c) 大地电磁数据RMS.Fig.13 RMS curve of cross-gradient joint inversion by different methods (a) RMS of gravity data; (b) RMS of gravity gradient data; (c) RMS of magnetotelluric data.

图14 复杂模型设计示意图 (a) 剩余密度模型三维透视图； (b) 对数电阻率模型三维透视图.Fig.14 Design diagram of complex model (a) 3D perspective of residual density model； (b) 3D perspective of log resistivity model.

图15 复杂模型的重力、重力梯度和大地电磁交叉梯度联合反演RMS曲线 (a) 重力数据RMS； (b) 重力梯度数据RMS； (c) 大地电磁数据RMS.Fig.15 RMS curve of gravity, gravity gradient and magnetotelluric cross-gradient joint inversion for complex model (a) RMS of gravity data; (b) RMS of gravity gradient data; (c) RMS of magnetotelluric data.

图16 (a) 真实密度模型切片1(x=3000 m)； (b) 真实密度模型切片2(y=5500 m)； (c) 真实电阻率模型切片1(x=3000 m)； (d) 真实电阻率模型切片2(y=5500 m)； (e) 联合反演结果密度切片1(x=3000 m)； (f) 联合反演结果密度 切片2(y=5500 m)； (g) 联合反演结果电阻率切片1(x=3000 m)； (h) 联合反演结果电阻率切片2(y=5500 m)Fig.16 (a) Real density slice 1 (x=3000 m)； (b) Real density slice 2 (y=5500 m)； (c) Real resistivity slice 1 (x=3000 m)； (d) Real resistivity slice 2 (y=5500 m)； (e) Density slice 1 of joint inversion (x=3000 m)； (f) Density slice 2 of joint inversion (y=5500 m)； (g) Resistivity slice 1 of joint inversion (x=3000 m)； (h) Resistivity slice 2 of joint inversion (y=5500 m)

表3 不同方式交叉梯度联合反演结果的RMS和RMSE值Table 3 RMS and RMSE of cross-gradient joint inversion results by different methods

4 结论

本文提出了一种基于平滑聚焦正则化法的三维重力、重力梯度和大地电磁数据交叉梯度联合反演算法，并通过理论模型试算对比分析得出以下结论：

(1)将交叉梯度结构相似约束应用于联合反演的重要前提是，参与联合的两种方法单独反演均能得到具有一定可靠性的物性模型.由于重力反演纵向分辨能力弱，其不可靠的纵向密度模型会因为结构相似性约束影响本来更可靠的纵向电阻率模型，导致联合反演结果准确性降低.本文通过引入重力梯度数据有效地提高了密度模型的可靠性，从而提高了联合反演的分辨率.

(2)本文提出的平滑聚焦正则化方法可以为大地电磁反演恢复出更真实的异常体几何形态和物性参数幅值，既能避免平滑反演异常体发散、边界模糊的问题，又能有效解决聚焦反演伴生的近地表假异常问题.同时，结合提出的聚焦强度阈值方法，可以显著地改善聚焦反演或平滑聚焦反演中易出现的过度聚焦问题，更好地恢复异常体特征.

(3)在交叉梯度联合反演中，若两种物性结构已近似一致，但均与真实模型有所偏差.此时，由于两种物性反演的分步性和交替性，连续的交叉梯度结构相似约束反而会限制反演模型向真实模型迭代的趋势.而本文提出的间断式联合方法，通过交替进行单独反演和交叉梯度联合反演，既能保证约束物性结构相似，又能允许迭代模型响应向观测数据充分拟合，可以获得更加收敛的联合反演结果.

(4)由于重力、重力梯度和大地电磁方法存在分辨率差异，交叉梯度联合之前，应使单独反演得到的密度异常和电阻率异常在形态上尽量相近.通过灵活调整反演参数来控制单独反演结果的异常发散程度，得到尽量大范围的有效交叉梯度，有利于改善最终联合反演结果.

本文限于条件仅采用了理论模型对算法进行研究，目前尚未将其应用于实测数据反演，接下来将进一步对此进行研究.

致谢感谢评审专家在百忙之中为本文修改提出的宝贵意见.