部分区域约束下的交叉梯度多重地球物理数据联合反演

李桐林，张镕哲，朴英哲,*，金钢燮，金长鲜，催邴敏

1 吉林大学地球探测科学与技术学院，长春　130026，中国 2 金策工业综合大学资源探测工学系，平壤　999093，朝鲜



部分区域约束下的交叉梯度多重地球物理数据联合反演

李桐林1，张镕哲1，朴英哲1,2*，金钢燮2，金长鲜2，催邴敏2

1 吉林大学地球探测科学与技术学院，长春130026，中国 2 金策工业综合大学资源探测工学系，平壤999093，朝鲜

交叉梯度联合反演方法通过对多种地球物理模型实现结构耦合，在岩石物性关系不确定的情况下，既能提高反演结果的可靠性，又能减少反演的多解性，还能减少不同方法解释结果之间的矛盾．当不同的模型观测数据覆盖范围不一致时，交叉梯度联合反演通常需要取出重叠区域数据进行联合反演，并且建模时还要扩展一些模型范围．本文首先提出并实现了部分区域约束下的交叉梯度多重地球物理数据联合反演算法；接着进行了算法的模型试算；最后，我们将该反演算法用于本溪—集安深部地质调查重磁电综合地质地球物理解释中．结果表明：该算法不但能在重叠区域内很好地恢复结构相似的模型，而且在非重叠区域与重叠区域的边界处仍然可以得到平滑变化的模型；在本溪—集安10号剖面所获得的结构上相似的电阻率、密度及磁化率模型较好地反映了该区的深部地质结构，对于确定深部地质体的性质提供了有力的证据．

交叉梯度；联合反演；大地电磁测深；初至波走时；重磁法

1　引言

在地球物理反演中，由于在空间和时间上的观测局限性等原因，始终存在多解性问题(Tarantola，2004)．减少多解性的通常方法是应用模型光滑约束(Constable et al．，1987)，其能够有效地消除不光滑构造信息，并得到可靠的反演结果(吴小平和徐果明，2000)．为了进一步减少反演结果不确定性和提高解释精度，需要结合多种地球物理方法．不同地球物理数据的结合是基于相应的物性之间的某些耦合关系．最简单的形式是对同一种物性的不同探测方法数据的联合解释(Vozoff and Jupp，1975；Sasaki，1989)．不同物性的耦合关系可以分成岩石物性关系和结构相似性关系．岩石物性关系因为建立了不同物性之间的函数关系(例如：在多孔介质中电阻率与地震速度间岩石物性特征(Tillmann and Stöcker，2000))所以在联合反演解释中起到明显的作用(Moorkamp et al．，2011)，但是在很多情况下，这种参数直接关系是不确定的，而且其稳定性也是未知的．结构相似性关系反映的是不同物性的空间变化关系，其基础是多种地球物理方法都能反映相同的地下构造及地质体(Gallardo and Meju，2011)．结构耦合联合反演最初是通过模型曲率来约束结构相似性(Haber and Oldenburg，1997；Zhang and Morgan，1997)．基于模型曲率的结构约束适合于边缘检测问题，但这种方法并没有考虑模型梯度的方向．Gallardo和Meju(2003)提出的另一种方法是用不同物性梯度的叉乘来识别结构边界的相似性．这种交叉梯度联合反演方法已经被用到地球物理综合解释中(Gallardo and Meju，2004；Gallardo，2007；Gallardo et al．，2012；Abubakar et al．，2012；Fregoso and Gallardo，2009；Moorkamp et al．，2013；Linde et al．，2008；Um et al．，2014；彭淼等，2013；Hamdan et al．，2012)及水文问题(Lochbühler et al．，2013)解释中．

交叉梯度联合反演在结构相似约束方式上有两种思路，一种思路为寻找不同模型之间严密的结构相似性，另一种思路为通过最小化结构相似项的模来实现结构约束(Gallardo and Meju，2011)．目前提出的严密结构相似方法都只在相同模型范围内进行联合反演．这意味着有以下两个问题：

第一，对于单独反演，地下模型的网格剖分规律取决于解释什么样的物理场．例如对直流电阻率法和大地电磁测深法反演，为了考虑边界条件的影响，二维问题时需要沿着4个模型边缘进行模型扩展来减少边界对主要解释区域的影响，然而初至波成像或重磁法反演不需要过于扩展模型．但在进行联合反演时，密度、磁化率或者纵波速度模型必须跟着电阻率模型进行扩展．

第二，当一种探测方法的测线可能是其他探测方法测线的一部分时，此时只能用重叠区域的观测数据进行联合反演．例如在Hamdan等(2012)的野外实测中，初至波测线是ERT测线的一部分，所以在进行联合反演时只用到ERT观测数据的一部分和全部的初至波走时数据．

本文首先介绍了基于严密结构相似的交叉梯度联合反演方法，然后提出了部分区域约束下的交叉梯度多重地球物理数据联合反演方法．最后通过模型试验和野外数据应用证明了该方法的合理性和实用性．

2　交叉梯度联合反演方法

(1)

约束τjl=0:∀j≠l，其中，m0j为mj的先验模型参数，Cdi和Cmj为观测数据和模型向量的协方差，di为第i种探测方法的观测数据，αj为阻尼系数．我们忽略了不同物性参数之间的岩石物性关系，只考虑结构上相似程度．在二维的情况下交叉梯度函数

(2)

只在走向方向含有非零成分．联合反演是对不同单位和不同量级的物性参数之间的耦合，因此，为了防止数值计算的不稳定性，在式(2)中加入系数κ来标准化模型向量．一般情况下设定κj为向量mj的预期变化幅度．联合反演之前，分别进行单独反演，然后单独反演结果的每一种模型参数的最大值和最小值的差即为对应的κ值．为了简单的描述公式推导，将式(1)以矩阵形式表示为

(3)

约束τ(m)=0，

其中

Cd=diag[Cd1,Cd2,…],

将式(3)的目标函数及约束项在当前模型附近进行泰勒级数展开，联合反演可简化为线性约束条件下的最小二乘问题．再用拉格朗日乘子方法(Tarantola，2004)，第二次模型更新式为

(4)

其中

ΔmLS为无交叉梯度约束的正则化最小二乘模型更新；B为交叉梯度函数的雅克比矩阵；拉格朗日乘子Λ为

(5)

ΔmLS和其协方差矩阵CLS为

(6)

(7)

其中

3　部分区域交叉梯度联合反演

3.1部分区域交叉梯度联合反演原理

部分区域联合反演是指在多个模型的重叠区域里交叉梯度为零的约束条件下进行的联合反演．因此目标函数跟上面的方法一致，只是式(1)的约束条件跟上面方法有区别．部分区域联合反演的交叉梯度约束条件为

(8)

图1　重叠模型区域内的交叉梯度约束Fig.1　Crossgradient constraint in the overlapped region of models

(9)

将式(9)在当前模型附近进行泰勒级数展开得到

其中

第j种和第l种模型间的交叉梯度对第j种模型参数的雅克比矩阵为

(10)

(11)

整个交叉梯度雅克比矩阵B为

(12)

图2　在包含电磁法数据的联合反演中电阻率模型与其他物性模型的关系Fig.2　Relations between resistivity model and other geophysical models for joint inversions including electromagnetic data

图3　部分区域交叉梯度联合反演算法Fig.3　Joint inversion algorithm with sub-region crossgradient constraints

我们可以用式(12)及式(4)—(7)计算出第二次迭代的模型近似m′．式(11)中的B和式(4)中的BTΛ只在对应重叠区域的模型单元内存在非零元素项，但在式(4)中，由于受到平滑模型约束的CLS影响，CLSBTΛ对不在重叠区域的模型单元内也存在非零数值，所以部分区域交叉梯度约束不仅对重叠区域元素而且对其他区域的元素都产生影响．换句话说联合反演在其他区域能得到观测数据拟合差和平滑模型约束的模型，而在重叠区域不但能得到以上两种约束条件下的模型还能得出结构相似的交叉梯度约束的模型．因此部分区域交叉梯度联合反演方法不但利用所有观测数据而且使模型在重叠区域满足结构相似性约束条件，还能在电磁法和重磁法及初至波走时的联合反演时，不需要密度、磁化率及速度模型跟着电阻率模型扩展(图2)．

3.2部分区域交叉梯度联合反演算法

部分区域交叉梯度联合反演的算法如图3所示由二重迭代组成．外部迭代以一样的比率逐渐减少阻尼系数，内部迭代对固定的阻尼系数寻求满足交叉梯度约束条件的解．在内部迭代为了提高计算效率不再计算模型响应fi(mj)和雅克比矩阵Aij(Gallardo et al．，2012)．

如果观测数据类别和模型参数类别个数相同并且每一种观测数据只跟一种物性有关，如图2所示，电阻率和大地电磁数据、密度和布格异常、磁化率和磁异常及纵波速度和初至波走时的联合反演时，因为Aij=0(∀i≠j)所以矩阵A变成分块对角矩阵．通过式(7)可以得到矩阵CLS也变成分块对角矩阵

(13)

为此式(6)可改写为

(14)

其中ΔmLSi为被第i种观测数据和正则化约束决定的最小二乘模型．用式(5)计算拉格朗日乘子时，系数矩阵BCLSBT维数较大，然而矩阵CLS为分块对角矩阵时可以用迭代解法(比如Minimal Residual Method(Barrett et al．，1994))避免构成大型逆矩阵．

做联合反演之前，对每一种探测方法做单独反演，然后从单独反演中得出阻尼系数的初值．阻尼系数的减少比率是通过阻尼系数在外部循环的最后迭代达到最小值得到的．根据我们的经验，αi最大值与最小值之间比率通常为100～1000，外部循环和内部循环迭代次数是5～20、3～5次合适．

4　模型试验与应用

4.1理论模型(音频大地电磁法和初至波走时数据的联合反演)

为了展示本文方法的优势，构建了与文献Hamdan等(2012)的野外数据例子相似的模型(图4)．电阻率模型为长度 8 km(水平距离1～9 km)、深度2 km，速度模型在位于水平距离约5～8 km的地方，深度为1.5 km．背景介质设定为电阻率300 Ωm，速度随着深度的变化线性增加的模型(图5)．两个矩形异常体的电阻率为30 Ωm及10000 Ωm，顶面埋深为100 m，底部埋深为300 m，其中第一个异常体不在速度模型区域里(图4)．在速度模型区域里的第二个异常体的速度设定为6000 m·s-1．大地电磁测量数据为TE、TM两种极化模式下的视电阻率和阻抗相位(9个测点(图4中黑色菱形)、15个频点(10～5000 Hz的对数间隔)，图6a所示的其中TE模式视电阻率的拟断面)，地震数据是由10个地震发射器(图5中倒三角)和30个接收器(图5中正三角)组成的观测系统记录的初至波走时(图6b)，并在所有数据中加入了2%的随机噪声．反演中初始电阻率模型为300 Ωm的均匀半空间，对于初始速度模型，为了保证折射波传播，设定为速度沿着深度方向递增的介质．用于大地电磁及初至波走时正演与雅克比矩阵计算的程序为Occam2DMT(deGroot-Hedlin and Constable，1990)和FAST(Zelt and Barton，1998)．

图4　理论模型电阻率网格剖分为灰色网格部分，虚线外部为扩展的电阻率单元，黑色实线为速度模型边缘，黑色为异常体的边缘．Fig.4　Synthetic modelThe mesh of resistivity model is plotted with gray lines.Outside of dashed lines are padded resistivity cells.Black solid lines are the edge of velocity model.The designed anomalies are outlined with black rectangles．

图5　速度模型射线分布Fig.5　Ray distribution in the velocity model

反演时，为了比较本文方法和常规交叉梯度联合反演方法，将速度模型的范围扩充至电阻率模型的范围，构建了范围一致的电阻率模型和速度模型．再考虑FAST方法对网格剖分的均匀大小要求和电磁场边界条件的影响，需要将电阻率模型和速度模型同时沿着水平方向和深度方向进行适当的模型扩展(分别10个模型单元)．对同一区域模型及不同区域模型分别进行了单独反演和联合反演．

无论是同一区域模型还是不同区域模型，电阻率单独反演都很明显的恢复了两个异常体(图7a和图8a)，但由于大地电磁测点之间的距离较大，第二个高阻异常体的宽度小于实际模型宽度．速度模型的射线分布主要集中在水平距离5～8 km和深度300 m以上(图5)，其他范围内没有射线分布，导致反演会在射线分布以外的区域敏感度降低．在速度单独反演模型中，由于速度异常体位于射线分布区域内，所以在地震单独反演剖面上可以很明显地看到高速异常体(图7b和图8b)．

在联合反演的电阻率剖面上，无论同一区域模型(图7c)还是不同区域模型(图8c)高阻异常体的宽度都比单独反演更接近模型真实值，这说明电阻率模型受到了速度模型结构的约束．对低阻异常体而言，因为速度模型上没有异常体从而未受到结构约束，所以联合反演模型比单独反演的没有改善．联合反演的速度模型(图7d和图8d)跟单独反演进行对比，速度异常体与背景介质的分离更明显，这说明了速度模型也受到了电阻率模型的高阻异常体结构的影响．在同一区域模型和不同区域模型的联合反演中，速度模型在深度1 km以下部分的速度分布差异(图7d和图8d)可以推测为不同模型范围的电阻率和速度灵敏度和模型平滑约束共同引起的．

图6　理论模型响应(a)TE模式视电阻率;(b)初至波走时.Fig.6　Synthetic model responses(a)TE mode apparent resistivity;(b)First-arrival traveltime.

图7　扩充模型的单独反演(a和b)与联合反演(c和d)结果．(a和c)电阻率反演模型，(b和d)速度反演模型Fig.7　Separate (a and b)and joint (c and d)inversion results of extended models.(a and c)resistivity models,(b and d)velocity models

图8　图4模型的单独反演(a和b)与联合反演(c和d)结果．(a和c)电阻率反演模型，(b和d)速度反演模型Fig.8　Separate (a and b)and joint (c and d)inversion results of Fig.4 models.(a and c)resistivity models,(b and d)velocity models

图9　本溪—吉安10号测线观测数据及反演模型响应(a)观测TM视电阻率；(b)图9a中电阻率模型的TM视电阻率响应；(c)和(d)分别为观测布格异常和磁异常及图9a中密度、磁化率模型的正演响应.Fig.9　Observations and model responses of the survey line No.10 in Benxi-Ji′an area(a)Observed TM mode apparent resistivity；(b)Response from the resistivity model of Fig.9a；(c)and(d)Observed Bouguer and TMI (Total Magnetic Intensity)anomaly,and responses from the density and susceptibility models of Fig.9a．

总体上，常规交叉梯度联合反演方法和本文方法都比单独反演更好地恢复异常体．但两个方法中存在差别．常规方法中，电阻率模型为了符合速度模型网格剖分规律用均匀大小的单元扩展了模型，而本文方法中，电阻率模型随着水平距离和深度的增加用逐渐变大的单元扩展模型，这是更符合大地电磁正反演精度的要求．在第11次迭代时，两种反演方法的拟合差都达到期望值1以下，反演时间分别为15分52秒，12分14秒(Intel Core2Duo E7200 2.53 GHz)，时间的增大是由速度模型扩展引起的.

4.2野外数据应用

本溪—集安地区横跨两个不同地质构造单元，一是太古代古陆核，由太古代变质岩组成，其控矿作用明显；二是元古代辽吉古裂谷带，其内矿产资源十分丰富．研究区太古代的含矿建造除了在研究区的北部出露地表外，大部都处于隐伏状态，元古代辽吉古裂谷带沉积巨厚并且被古生代、中新生代沉积盖层覆盖，且无论是太古代含矿建造还是辽吉古裂谷带都经历了多期构造与岩浆侵入活动．查明该区深部地质结构与构造，将会对研究区找矿具有非常重要的意义(张宏嘉等，2013)．为此在研究区共完成了十几条重、磁、大地电磁综合剖面．使用的仪器是加拿大凤凰公司V5-2000系列电磁仪、美国BURRIS型高精度金属弹簧重力仪和加拿大GEM公司GSM-19T微机质子磁力仪．此外对研究区涉及的主要地质体的物性也进行了测量．

用本文方法处理了其中10号测线数据(图9)．测线方向大致从北西到南东，测线长度为97 km．在测线上共有大地电磁法48个测点，测点之间平均距离为2 km，测量数据中采取了15个频点(360.0～0.02 Hz的对数间隔)的TM极化模式视电阻率(图9a)．从重磁观测数据得到布格异常和磁异常，其中采取189个测点数据(图9c和图9d中的红色圆形)．

密度与磁化率模型的深度设定为20 km．对于电阻率模型，考虑电磁场的衰减，用对数间隔单元沿着水平方向和深度方向扩展了模型(图10)．主要解释区域(图10中虚线表示的区域)网格单元大小设定为宽1100 m，高200～4000 m．反演初始模型为300 Ωm、0 g·cm-3及0 SI的均匀介质．重磁正演及雅克比矩阵计算根据多边形异常体公式(Singh，2002)．

由于此次研究主要关心深度10 km以内的地质结构，因此图11a中只画出了深度到10 km的剖面．联合反演过程在第6次迭代时以1.52(电)、0.81(重)、1.38(磁)的数据拟合差结束．最终模型响应分别表示为图9b的拟断面(TM视电阻率)和图9c(布格异常)、图9d(磁异常)的蓝色实线．

联合反演获得了结构上相似的电阻率、密度及磁化率模型(图11a)．根据物性范围，可清晰地划分出不同的地质体．结合地表地质情况(图11b)和物性测量数据作进一步解释，得到每个地质体的位置和物性范围，推断了各地质体的属性以及断层构造．综合解释剖面图(图11c)与10号测线附近地表地质(图11b)和联合反演剖面(图11a)基本吻合，其中水平位置72～76 km、深度1～6 km的地质体没有出现地表地质图上，但是根据其物性范围和在地质图外部浪子山组邻接的地层及侵入岩属性，推测该地质体为古元古代盖县组．10号测线剖面可以分出两个地质单元，第一段为0～50 km区间的桓仁地块，主要岩性由早白垩纪梨树沟组、小岭组火山沉积层、古元古代大石桥组组成；第二段为50～95 km区间的集安地块，主要岩性为早白垩纪小岭组火山沉积层、古元古代盖县组、浪子山组、早白垩纪侵入岩．测线下伏岩体为古元古代时期生成的辽吉花岗岩．根据综合解释剖面和联合反演结果推测的每个地质体的物性范围如表1所示．

图10　联合反演模型网格剖分虚线为密度及磁化率模型的边缘，灰色网格部分为电阻率模型剖分，交叉梯度约束只适用于虚线区域内．Fig.10　Mesh generation of joint inversion modelDashedline is the edge of density and susceptibility model.Gray line is the mesh of resistivity model.Crossgradient constraint is only applied to the dashed line rectangle.

图11　10号测线联合反演结果(a)、地质图(b)及联合反演综合解释剖面图(c)Fig.11　(a)Joint inversion results of the survey line No.10，(b)geological map and (c)comprehensively interpreted section of the joint inversion result

地质体电阻率log10Ωm密度异常g·cm-3磁化率异常SI地质属性K1l2.0～2.3-0.03～0.01-0.017～0.000早白垩纪梨树沟组Jkx1.7～2.2-0.12～0.000.000～0.030晚侏罗纪—早白垩纪小岭组Pt1gx3.7～3.9-0.02～-0.01-0.001～0.000古元古代盖县组Pt1d1.8～2.50.02～0.14-0.029～-0.004古元古代大石桥组Pt1l3.6～3.90.00～0.050.000～0.007古元古代浪子山组K1Fh3.5～3.6-0.02～-0.010.003～0.025早白垩纪花岗岩K1Lj3.5～3.8-0.01～0.030.000～0.015早白垩纪连江辉石闪长岩Pt2κγ3.3～3.5-0.06～-0.01-0.021～0.001中元古代碱长花岗岩Pt1Hrγ3.3～3.8-0.03～-0.02-0.035～-0.013古元古代桓仁花岗杂岩

5　结论

针对联合反演中不同地球物理模型范围不一致的问题，本文提出并实现了部分区域约束下的交叉梯度多重联合反演，并进行了模型试算与应用试验，得到了如下结论：

(1)本文提出的算法仅在重叠区域内应用交叉梯度约束，克服了所有模型区域需要一致的限制，能对不相同模型范围进行结构耦合联合反演．

(2)理论模型结果表明：当观测数据覆盖的区域不完全一致时，本文算法不但可以用到所有观测数据，而且能在重叠区域内很好地恢复结构耦合的结果．除此之外，由于平滑模型约束的影响，在非重叠区域与重叠区域的边界处仍然可以得到平滑变化的模型．

(3)应用试验得到的电阻率、密度及磁化率模型结构上很相似，较好地反映了研究地区的深部构造，为确定深部地质体性质提供了有力证据．

致谢感谢匿名评审专家为本文修改提出的细致而宝贵的意见．

References

Abubakar A，Gao G，Habashy T M，et al．2012．Joint inversion approaches for geophysicalelectromagnetic and elastic full-waveform data．Inverse Problems，28(5)：055016，doi：10.1088/0266-5611/28/5/055016．

Barrett R，Berry M，Chan T F，et al．1994．Templates for the Solution of Linear Systems:Building Blocks for Iterative Methods．Philadelphia：SIAM．

Constable S C，Parker R L，Constable C G．1987．Occam′s Inversion：A practical algorithm for generating smooth models from electromagnetic sounding data．Geophysics，52(3)：289-300，doi:10.1190/1.1442303．

deGroot-Hedlin C D，Constable S C．1990．Occam′s Inversion to Generate Smooth Two Dimensional Models from Magnetotelluric Data．Geophysics，55(12)：1613-1624，doi：10.1190/1.1442813．

Fregoso E，Gallardo L A．2009．Cross-gradients joint 3D inversion with applicationsto gravity and magnetic data．Geophysics，74(4)：L31-L42，doi：10.1190/1.3119263．Gallardo L A，Meju M A．2003．Characterization of heterogeneousnear-surface materials by joint 2D inversion of DC resistivity and seismicdata．Geophys．Res．Lett．，30(13)，doi：10.1029/2003GL017370．Gallardo L A，Meju M A．2004．Joint two-dimensional DC resistivityand seismic travel time inversion with cross-gradients constraints．J．Geophys．Res．，109(B3)：B03311，doi：10.1029/2003JB002716．

GallardoLA．2007．Multiple cross-gradient joint inversion for geospectralimaging．Geophys．Res．Lett．，34(19)：L19301-5，doi：10.1029/2007GL030409．Gallardo LA，MejuMA．2011．Structure-coupled multiphysics imagingin geophysical sciences．Rev．Geophys．，49(1)：RG1003，doi：10.1029/2010RG000330．Gallardo LA，FontesSL，MejuMA，et al．2012．Robust geophysical integration through structure-coupled joint inversionand multispectral fusion of seismic reflection，magnetotelluric，magnetic，and gravity images:Example from Santos Basin，offshore Brazil．Geophysics，77(5)：B237-B251，doi：10.1190/GEO2011-0394.1．Haber E，Oldenburg D．1997．Joint inversion:A structuralapproach．Inverse Problems，13(1)：63-77，doi：10.1088/0266-5611/13/1/006．

Hamdan H，Economou N，Kritikakis G，et al．2012．2D and 3D imaging of the metamorphic carbonates at Omalos plateau/polje，Crete，Greece by employing independent and joint inversion on resistivity and seismic data．International Journal of Speleology，41(2)：199-209，doi：10.5038/1827-806X.41.2.7．

Linde N，Tryggvason A，Peterson J E，et al．2008．Joint inversion of crosshole radar and seismic traveltimesacquired at the South Oyster Bacterial Transport Site．Geophysics，73(4)：G29-G37，doi：10.1190/1.2937467．

Lochbühler T，Doetsch J，Brauchler R，et al．2013．Structure-coupled joint inversion of geophysical and hydrological data．Geophysics，78(3)：ID1-ID14，doi：10.1190/GEO2012-0460.1．

Moorkamp M，Heincke B，Jegen M，et al．2011．A framework for 3D joint inversion of MT，gravity and seismicrefraction data．Geophys．J．Int．，184(1)：477-493，doi：10.1111/j.1365-246X.2010.04856.x．

Moorkamp M，Roberts A W，Jegen M，et al．2013．Verification of velocity-resistivity relationships derived fromstructural joint inversion with borehole data．Geophys．Res．Lett．，40(14)：3596-3601，doi：10.1002/grl.50696．

Peng M，Tan H D，Jiang M，et al．2013．Three-dimensional joint inversion of magnetotelluric and seismic traveltime data with cross-gradient constraints．Chinese J．Geophys．(in Chinese)，56(8)：2728-2738，doi：10.6038/cjg20130821．

Sasaki Y．1989．Two-dimensional joint inversion of magnetotelluric and dipole-dipole resistivity data．Geophysics，54(2)：254-262，doi：10.1190/1.1442649．

Singh B．2002．Simultaneous computation of gravity and magnetic anomalies resulting from a 2-D object．Geophysics，67(3)：801-806，doi：10.1190/1.1484524．

TarantolaA．2004．Inverse Problem Theory and Methods for Model Parameter Estimation．Philadelphia：SIAM．

Tillmann A，Stöcker T．2000．A new approach for the joint inversion of seismic and geoelectric data．//Extended Abstracts，63rdEAGE Conference and Technical Exhibition．Amsterdam：EAGE．

Um E S，Commer M，Newman G A．2014．A strategy for coupled 3D imaging of large-scale seismic andelectromagnetic data sets:Application to subsalt imaging．Geophysics，79(3)：ID1-ID13，doi：10.1190/GEO2013-0053.1．

Vozoff K，Jupp DLB．1975．Joint inversion of geophysical data．Geophys．J．Int.，42(3)：977-991，doi：10.1111/j.1365-246X.1975.tb06462.x．

Wang J,Meng X H,Chen Z X.2013.The theory of cross-gradient and its application in geophsical joint inversion.Chinese J.Geophys.(inChinese),28(4):2094-2103,doi:10.6038/pg20130454.Wu X P，Xu G M．2000．Study on 3-D resistivity inversion using conjugate gradient method．Chinese J．Geophys．(in Chinese)，43(3)：420-427．

Zelt C A，Barton P J．1998．Three-dimensional seismic refraction tomography：A comparison of two methods applied to datafrom the Faeroe Basin．J．Geophys．Res．，103(B4)：7187-7210，doi：10.1029/97JB03536．

Zhang H J，Cao L，Li F W，et al．2013．GME comprehensive explanation of 3D geological structure of Benxi-Ji′an region．Jilin Geology (in Chinese)，32(1)：88-94．

Zhang J，Morgan FD．1997．Joint seismic and electrical tomography.//Symposium on the Applications of Geophysicsto Engineering and Environmental Problems．Keystone，Colo：EEGS，391-396．

附中文参考文献

彭淼，谭捍东，姜枚等．2013．基于交叉梯度耦合的大地电磁与地震走时资料三维联合反演．地球物理学报，56(8)：2728-2738，doi：10.6038/cjg20130821．

王俊,孟小红,陈召曦等.2013.交叉梯度理论及其在地球物理联合反演中的应用.地球物理学进展,28(4):2094-2103,doi:10.6038/pg20130454.

吴小平，徐果明．2000．利用共轭梯度法的电阻率三维反演研究．地球物理学报，43(3)：420-427．

张宏嘉，曹亮，李福文等．2013．本溪—集安地区三维地质结构重磁电综合解释．吉林地质，32(1)：88-94．

(本文编辑刘少华)

Multiple joint inversion of geophysical data with sub-region crossgradient constraints

LI Tong-Lin1，ZHANG Rong-Zhe1，PAK Yong-Chol1,2*，KIM Gang-Sop2，KIM Jang-Son2，CHOE Byong-Min2

1 College of Geo-Exploration Sciences and Technology,Jilin University,Changchun 130026,China 2 Resource Exploration and Engineering Department,Kimchaek University of Technology,Pyongyang 999093, Democratic People′s Republic of Korea

Crossgradient joint inversions,achieving structure coupling of multiple geophysical models,improve the reliability of inversion results,reduce the uncertainty of inversions,and reduce conflicts between results of different methods,in the case of unknown petrophysical relationships.When coverage of different observations is inconsistent,the crossgradient joint inversion usually need extract data corresponding to the overlapped region,and also expands the scope of some models.This paper proposed and implemented a multiple joint inversion of geophysical data with sub-region crossgradient constraints and followed by a synthetic model example.We have applied this method to comprehensive geological and geophysical interpretation of MT,gravity and magnetic field data sets for deep geological survey of Benxi-Ji′an.The results show that this algorithm not only can well restore similar structure model in the overlapped region,but the model parameters are varied smoothly on the boundary between overlapped and non-overlapped regions.The resistivity,density and susceptibility models,obtained from the joint inversion of Benxi-Ji′an survey line No.10,are very similar in their structures,well reflect the deep geological structure of this area,and provide strong evidences to determine the nature of deep geological bodies.

Crossgradient；Joint inversion；Magnetotellurics；First-arrival traveltime；Gravity and magnetic method

李桐林，张镕哲，朴英哲等.2016.部分区域约束下的交叉梯度多重地球物理数据联合反演.地球物理学报，59(8):2979-2988,

10.6038/cjg20160821.

Li T L,Zhang R Z,Pak Y C,et al.2016.Multiple joint inversion of geophysical data with sub-region crossgradient constraints.Chinese J.Geophys.(in Chinese),59(8):2979-2988,doi:10.6038/cjg20160821.

中国地质调查局项目(12120113098400)，国家重大科学仪器设备开发专项(2011YQ05006009).

李桐林，男，1962年生，教授，主要从事电磁法理论与应用教学和研究.E-mail:lilaoshizh@163.com

朴英哲，男，朝鲜人，1980年生，博士研究生，主要从事电磁正反演理论与应用研究，E-mail:piaoyz525@163.com

10.6038/cjg20160821

P318，P319

2015-12-12，2016-05-15收修定稿