2016, Vol. 59
2. 金策工业综合大学资源探测工学系, 平壤 999093, 朝鲜
2. Resource Exploration and Engineering Department, Kimchaek University of Technology, Pyongyang 999093, Democratic People's Republic of Korea
在地球物理反演中，由于在空间和时间上的观测局限性等原因，始终存在多解性问题(Tarantola，2004)．减少多解性的通常方法是应用模型光滑约束(Constable et al.，1987)，其能够有效地消除不光滑构造信息，并得到可靠的反演结果(吴小平和徐果明，2000)．为了进一步减少反演结果不确定性和提高解释精度，需要结合多种地球物理方法．不同地球物理数据的结合是基于相应的物性之间的某些耦合关系．最简单的形式是对同一种物性的不同探测方法数据的联合解释(Vozoff and Jupp，1975；Sasaki，1989)．不同物性的耦合关系可以分成岩石物性关系和结构相似性关系．岩石物性关系因为建立了不同物性之间的函数关系(例如：在多孔介质中电阻率与地震速度间岩石物性特征(Tillmann and Stcker，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 交叉梯度联合反演方法假设参与联合反演的模型类别个数和观测方法类别个数都是任意的，而且每个观测数据类别对所有模型参数类别都敏感．令 m j为第j 种模型参数向量，引进扩展模型向量 m = m 1T，m 2T，…]T，第i种方法正演理论响应的表达形式为 f i(m)． 多种地球物理方法交叉梯度联合反演的目标函数是由模型粗糙度和数据拟合差构成，并在不同物性参数梯度的叉乘为零的条件下进行约束(Gallardo and Meju，2004；Gallardo et al.，2012)：
|
(1) |
约束 τ jl=0:j≠l，其中，m 0j为 m j的先验模型参数，C di和C mj为观测数据和模型向量的协方差，d i为第i种探测方法的观测数据，αj 为阻尼系数．我们忽略了不同物性参数之间的岩石物性关系，只考虑结构上相似程度．在二维的情况下交叉梯度函数
|
(2) |
只在走向方向含有非零成分．联合反演是对不同单位和不同量级的物性参数之间的耦合，因此，为了防止数值计算的不稳定性，在式(2)中加入系数κ来标准化模型向量．一般情况下设定κj为向量 m j 的预期变化幅度．联合反演之前，分别进行单独反演，然后单独反演结果的每一种模型参数的最大值和最小值的差即为对应的κ值．为了简单的描述公式推导，将式(1)以矩阵形式表示为
|
(3) |
约束 τ(m)=0，
其中
|
将式(3)的目标函数及约束项在当前模型附近进行泰勒级数展开，联合反演可简化为线性约束条件下的最小二乘问题．再用拉格朗日乘子方法(Tarantola，2004)，第二次模型更新式为
|
(4) |
其中
|
Δ m LS为无交叉梯度约束的正则化最小二乘模型更新； B 为交叉梯度函数的雅克比矩阵；拉格朗日乘子 Λ 为
|
(5) |
Δ m LS和其协方差矩阵 C LS为
|
(6) |
|
(7) |
其中
|
部分区域联合反演是指在多个模型的重叠区域里交叉梯度为零的约束条件下进行的联合反演．因此目标函数跟上面的方法一致，只是式(1)的约束条件跟上面方法有区别．部分区域联合反演的交叉梯度约束条件为
|
(8) |
其中
|
图 1 重叠模型区域内的交叉梯度约束 Fig. 1 Crossgradient constraint in the overlapped region of models |
|
(9) |
将式(9)在当前模型附近进行泰勒级数展开得到
|
其中
|
第j种和第l种模型间的交叉梯度对第j种模型参数的雅克比矩阵为
|
(10) |
将式(10)跟原方法的交叉梯度雅克比矩阵 B jlj比较，两个雅克比矩阵之间的关系为
|
(11) |
整个交叉梯度雅克比矩阵 B 为
|
(12) |
我们可以用式(12)及式(4)—(7)计算出第二次迭代的模型近似 m ′．式(11)中的 B和式(4)中的 B T Λ 只在对应重叠区域的模型单元内存在非零元素项，但在式(4)中，由于受到平滑模型约束的 C LS影响，C LS B T Λ 对不在重叠区域的模型单元内也存在 非零数值，所以部分区域交叉梯度约束不仅对重叠区域元素而且对其他区域的元素都产生影响．换句话说联合反演在其他区域能得到观测数据拟合差和平滑模型约束的模型，而在重叠区域不但能得到以上两种约束条件下的模型还能得出结构相似的交叉梯度约束的模型．因此部分区域交叉梯度联合反演方法不但利用所有观测数据而且使模型在重叠区域满足结构相似性约束条件，还能在电磁法和重磁法及初至波走时的联合反演时，不需要密度、磁化率及速度模型跟着电阻率模型扩展(图 2)．
|
图 2 在包含电磁法数据的联合反演中电阻率模型与其他物性模型的关系 Fig. 2 Relations between resistivity model and other geophysical models for joint inversions including electromagnetic data |
部分区域交叉梯度联合反演的算法如图 3所示由二重迭代组成．外部迭代以一样的比率逐渐减少 阻尼系数，内部迭代对固定的阻尼系数寻求满足交 叉梯度约束条件的解．在内部迭代为了提高计算效率不再计算模型响应 f i(m j)和雅克比矩阵 A ij(Gallardo et al.，2012)．
|
图 3 部分区域交叉梯度联合反演算法 Fig. 3 Joint inversion algorithm with sub-region crossgradient constraints |
如果观测数据类别和模型参数类别个数相同并且每一种观测数据只跟一种物性有关，如图 2所示，电阻率和大地电磁数据、密度和布格异常、磁化率和磁异常及纵波速度和初至波走时的联合反演时，因为 A ij=0(i≠j)所以矩阵 A 变成分块对角矩阵． 通过式(7)可以得到矩阵 C LS也变成分块对角矩阵
|
(14) |
其中Δ m LSi为被第i种观测数据和正则化约束决定的最小二乘模型．用式(5)计算拉格朗日乘子时，系数矩阵 BC LS B T维数较大，然而矩阵 C LS为分块 对角矩阵时可以用迭代解法(比如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 model The 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 |
|
图 6 理论模型响应 (a)TE模式视电阻率;(b)初至波走时. Fig. 6 Synthetic model responses (a)TE mode apparent resistivity;(b)First-arrival traveltime. |
反演时，为了比较本文方法和常规交叉梯度联合反演方法，将速度模型的范围扩充至电阻率模型的范围，构建了范围一致的电阻率模型和速度模型．再考虑FAST方法对网格剖分的均匀大小要求和电磁场边界条件的影响，需要将电阻率模型和速度模型同时沿着水平方向和深度方向进行适当的模型扩展(分别10个模型单元)．对同一区域模型及不同区域模型分别进行了单独反演和联合反演．
无论是同一区域模型还是不同区域模型，电阻率单独反演都很明显的恢复了两个异常体(图 7a和图 8a)，但由于大地电磁测点之间的距离较大，第二个高阻异常体的宽度小于实际模型宽度．速度模型的射线分布主要集中在水平距离5~8 km和深度300 m以上(图 5)，其他范围内没有射线分布，导致反演会在射线分布以外的区域敏感度降低．在速度单独反演模型中，由于速度异常体位于射线分布区域内，所以在地震单独反演剖面上可以很明显地看到高速异常体(图 7b和图 8b)．
|
图 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 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 |
在联合反演的电阻率剖面上，无论同一区域模型(图 7c)还是不同区域模型(图 8c)高阻异常体的 宽度都比单独反演更接近模型真实值，这说明电阻率模型受到了速度模型结构的约束．对低阻异常体而言，因为速度模型上没有异常体从而未受到结构约束，所以联合反演模型比单独反演的没有改善．联合反演的速度模型(图 7d和图 8d)跟单独反演进行对比，速度异常体与背景介质的分离更明显，这说明 了速度模型也受到了电阻率模型的高阻异常体结构的影响．在同一区域模型和不同区域模型的联合反演中，速度模型在深度1 km以下部分的速度分布差异(图 7d和图 8d)可以推测为不同模型范围的电阻率和速度灵敏度和模型平滑约束共同引起的．
总体上，常规交叉梯度联合反演方法和本文方法都比单独反演更好地恢复异常体．但两个方法中存在差别．常规方法中，电阻率模型为了符合速度模型网格剖分规律用均匀大小的单元扩展了模型，而本文方法中，电阻率模型随着水平距离和深度的增加用逐渐变大的单元扩展模型，这是更符合大地电 磁正反演精度的要求．在第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中的红色圆形)．
|
图 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． |
密度与磁化率模型的深度设定为20 km．对于电阻率模型，考虑电磁场的衰减，用对数间隔单元沿着水平方向和深度方向扩展了模型(图 10)．主要解释区域(图 10中虚线表示的区域)网格单元大小设定为宽1100 m，高200~4000 m．反演初始模型为300 Ωm、0 g·cm-3及0 SI的均匀介质．重磁正演及 雅克比矩阵计算根据多边形异常体公式(Singh，2002)．
|
图 10 联合反演模型网格剖分虚线为密度及磁化率模型的边缘，灰色网格部分为电阻率模型剖分，交叉梯度约束只适用于虚线区域内． Fig. 10 Mesh generation of joint inversion model Dashedline 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. |
由于此次研究主要关心深度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所示．
|
图 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 |
|
表 1 地质体物性范围及属性 Table 1 Property ranges of geological bodies and their attributes |
针对联合反演中不同地球物理模型范围不一致的问题，本文提出并实现了部分区域约束下的交叉梯度多重联合反演，并进行了模型试算与应用试验，得到了如下结论：
(1) 本文提出的算法仅在重叠区域内应用交叉梯度约束，克服了所有模型区域需要一致的限制，能对不相同模型范围进行结构耦合联合反演．
(2) 理论模型结果表明：当观测数据覆盖的区域不完全一致时，本文算法不但可以用到所有观测数据，而且能在重叠区域内很好地恢复结构耦合的结果．除此之外，由于平滑模型约束的影响，在非重叠区域与重叠区域的边界处仍然可以得到平滑变化的模型．
(3) 应用试验得到的电阻率、密度及磁化率模型结构上很相似，较好地反映了研究地区的深部构造，为确定深部地质体性质提供了有力证据．
| Abub akar 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–L19305. 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. , 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 , 63. | |
| 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. | |
| 王俊, 孟小红, 陈召曦, 等. 2013. 交叉梯度理论及其在地球物理联合反演中的应用. 地球物理学进展 , 28(4): 2094–2103. | |
| 吴小平, 徐果明. 2000. 利用共轭梯度法的电阻率三维反演研究. 地球物理学报 , 43(3): 420–427. | |
| 张宏嘉, 曹亮, 李福文, 等. 2013. 本溪-集安地区三维地质结构重磁电综合解释. 吉林地质 , 32(1): 88–94. | |