尽管单一物探方法的采集、处理和解释运用已经相对较为成熟,但是随着科学技术的提高,要求观测精度也随之提高,并且空间和时间观测存在的局限性,观测误差的不可避免性,导致反演过程中始终存在着多解性的问题.为了降低地球物理数据反演的多解性,可采用对反演的目标函数引入正则化约束的反演(Rodi and Mackie, 2001; Key, 2016),使得每次迭代过程中目标函数的极小值唯一,从而达到降低多解性的目的.同时,由于不同地球物理方法干扰来源不同、受干扰的影响程度也不同,可以采用多种地球物理数据进行联合反演,进而起到压制噪声和减少反演结果多解性的作用.
联合反演作为地球物理反演的新思路,其具体含义是利用不同类型的地球物理数据约束反演,提高反演结果的准确性的目的.不同地球物理数据的噪声来源以及对噪声的敏感程度不尽相同.虽然利用多种方法的数据对同一地质体的综合解释已经取得不错的效果(高德章等,2004;李德春等,2012;姜枚等,2012),但是采用更多的物性联合反演与增加同一物性的数据量进行反演相比,更有利于改善反演效果.随后诸多学者展开了广泛并且深入的研究,Sasaki(1989)基于相同的物性参数,进行了大地电磁法和偶极-偶极直流电法的联合反演.
在国际上比较流行的是采用模型结构耦合的联合反演,具有代表性的是Gallardo和Meju(2003)提出的基于交叉梯度约束的耦合方法, 利用交叉梯度函数约束直流电测深数据和地震折射数据之间的联合反演,在处理实测数据方面取得了较好的效果. Colombo和De Stefano(2007)基于不同的物性参数,完成了地震和重力、地震和大地电磁数据的同步联合反演,改善了地震反演的速度模型;Jegen等(2009)将联合反演应用于海洋勘探领域并取得了较好的效果.彭淼(2013)实现了基于交叉梯度约束的大地电磁和地震走时资料的三维联合反演,他的成果表明在联合反演对电阻率异常的边界刻画和异常值恢复都具有了较大的改善.另外,Wang等(2017)实现了基于交叉梯度约束的CSAMT和磁法二维联合反演,研究结果显示联合反演对于磁化率和电阻率的反演改善取得良好的效果.
本文在前人提出的交叉梯度的基础上,着重完善和论述了三维交叉梯度函数离散化方式,开发了基于交叉梯度约束的重力、磁法和大地电磁三维联合反演算法,并实现MPI并行加速算法,大大改善了三维大地电磁反演速度慢、时间长的问题.大地电磁的勘探深度较深能够弥补重力和磁法在深度上呈现指数衰减导致勘探深度浅的情况,而重、磁可以补充并丰富大地电磁在浅部勘探的地质构造信息,利用不同观测数据能够弥补彼此的敏感度不足的缺点,取长补短.通过理论模型的测试,探讨了联合反演算法的可行性和适用性.
1 重力、磁法和大地电磁三维正反演方法 1.1 重、磁三维正反演方法重力正演问题存在解析解公式(曾华霖, 2005);在给定地磁场强度、地磁场倾角、断面磁倾角、地磁场与断面夹角等参数后,磁法正演公式可简化为(管志宁, 2005)
(1) |
其中,d为观测数据向量(Δg重力异常,ΔT磁异常),G为核函数矩阵,m为模型向量(密度ρ,磁化率κ).重、磁的正演问题就是已知G和m求取d的过程.
由于密度、磁化率的变化范围较大并且反演过程中可能出现负值,本文使用Li和Oldenburg(1996)提出来的对数障碍法,选择以“10”为底模型向量的对数mk作为磁法反演的物性参数:
(2) |
由此,重、磁的正演公式变为
(3) |
本文采取的重、磁三维反演方法是Li和Oldenburg(1996, 2003)提出的基于对数障碍法的正则化反演算法,在考虑数据拟合的基础上加入了模型光滑约束、深度加权和对数障碍法.其中,模型光滑约束解决了反演结果的趋肤效应问题;深度权函数抵消了模型光滑矩阵在目标函数极小化过程中造成的物性结构上移;对数障碍约束解决了反演结果超出合理范围的问题.读者可参考Li和Oldenburg等相关参考文献,本文不再赘述.
1.2 大地电磁的三维正反演方法求解大地电磁正演问题所常使用的数值模拟方法分为有限差分法、有限单元法和积分方程法(柳建新等, 2012).本文联合反演中大地电磁的正演采用的是三维交错网格有限差分算法(谭捍东等,2003).
大地电磁反演方法采用了L-BFGS方法(Nocedal and Wright, 2006),即有限内存拟牛顿法.L-BFGS方法在BFGS利用迭代时目标函数的梯度巧妙地逼近当前海森矩阵逆的基础上进行了改善,采用了前几次(3~20)次迭代的目标函数梯度逼近,具有高效快速的优点,大大改善了三维大地电磁反演时间长的问题(Wang et al., 2017).
根据Egbert和Kelbert(2012)的模型协方差矩阵和修订的目标函数进行反演,修改后的大地电磁目标函数形式如下:
(4) |
Egbert提出了一种简单而且效率颇高的方法来光滑原始参数梯度的仿射线变换,将上述的目标函数进一步修正, 而且在王堃鹏的博士论文(2017)里也对这种变换进行了详细论述:
(5) |
修正后的目标函数:
(6) |
修正后的梯度表达式:
(7) |
对于真实的模型参数可以由下式获得:
(8) |
L-BFGS反演流程(图 1):
① 选择初始模型m0,设定最大的迭代次数IterMax,最小拟合差Rms;②根据LBFGS理论及公式计算拟海森矩阵的逆
③ 设定初始步长,并从初始步长开始进行一维线搜索,以寻找当前迭代相对最佳步长α;mk+1=mk+αpk.
④ 计算当前模型的拟合差是否小于设定的最小拟合差,若小于则停止当前迭代,若迭代次数大于设定的最大迭代次数,同样停止迭代.反之,则跳转到第二步进行反演的迭代.
2 基于交叉梯度耦合的三维联合反演算法 2.1 三维交叉梯度函数以两种物性ms,mr的交叉梯度函数为例,交叉梯度函数定义式如下:
其中tx, ty, tz具体形式如下:
(9) |
(10) |
(11) |
交叉梯度值Φcg表达式:
(12) |
彭淼(2012)推导了三维交叉梯度函数,将tx, ty, tz的表达式进行离散化,并未对不同网格离散化进行详细讨论,为了寻求更为合理并且精度较高的离散化方式,本文依据中心网格所在的位置我们分成以下四类网格进行离散化:非边界网格、边界面网格、边界线网格和边界点网格.
其中,m的下角标r,s代表两种物性,上角标C, R, L, N, P, B, U代表离散化过程中网格位置,具体含义如图 2所示.
① 非边界网格(如图 3所示)
(13) |
(14) |
(15) |
② 边界面网格(以边界面x=1为例,如图 4所示)
(16) |
(17) |
(18) |
③ 边界线网格(以边界线x=1,z=1为例,如图 5所示)
(19) |
(20) |
(21) |
④ 边界点网格(以边界点x=1,y=1,z=1为例,如图 6所示)
(22) |
(23) |
(24) |
为了获得稳定而且足够精确的正演解,不同方法对正演模拟的网格剖分具有不同的要求.考虑到大地电磁方法的分辨能力、场的衰减特点以及正演计算的效率,大地电磁的网格一般按照不等间隔进行剖分,在Z方向随着深度的增加网格间距逐渐变大,在X、Y方向上,远离研究区域的外围网格也逐渐变得稀疏.而对于重、磁正反演过程中网格没有特殊要求,一般是将网格均等剖分进行计算.但在联合反演过程中,为了对不同模型的物性参数进行耦合计算,这就要求不同方法的网格在我们关心的研究区域内需要一致,以此方便进行交叉梯度的计算(如图 7所示).
将加入交叉梯度后的重力、磁法和大地电磁的目标函数定义如下:
(25) |
(26) |
(27) |
式中,σ、κ、ρ分别为剩余密度、剩余磁化率和电阻率.G为灵敏度矩阵,Cm为模型协方差,Cd为数据协方差,λ为拉格朗日乘子,η1、η2、η3、η4、η5、η6为交叉梯度因子.
由此可见,在加入交叉梯度之前,各个方法之间没有相互联系.通过交叉梯度函数耦合三种物性,在单方法反演过程中受到其他两种方法反演结果的影响,增强了各方法之间反演结果在结构上的一致性.
2.4 联合反演目标函数的梯度表达式周丽芬(2012)给出了利用泰勒级数展开求取交叉梯度项关于模型向量偏导数的方法,本文依旧采用这种方案,当反演的初始模型为均匀半空间时,交叉梯度的偏导数矩阵可以与模型协方差矩阵进行合并,以联合反演重力目标函数的梯度表达式为例,具体的推导过程如下:
(28) |
将目标表达式对模型求取偏导得到梯度表达式:
(29) |
其中,
设Bκ为交叉梯度t(σ, κ)对密度模型σ的偏导数矩阵:
(30) |
当反演的初始模型为均匀半空间时,t0=0,
(31) |
同理分析:
(32) |
综上所述,联合反演重力目标函数的梯度表达式如下:
(33) |
同理所证,联合反演磁法目标函数的梯度表达式如下:
(34) |
联合反演大地电磁目标函数如下:
(35) |
修正后的梯度表达式:
(36) |
(37) |
同理所证:
(38) |
综上所述,联合反演大地电磁目标函数的梯度表达式如下:
(39) |
基于以上反演及交叉梯度的理论和算法,采用Fortran语言开发了重力、磁法和大地电磁三维联合反演算法,并实现MPI并行加速计算.为了检验联合反演算法的有效性和优越性,本文进行了理论模型合成数据的反演研究,分别设计了单棱柱体模型算例和组合棱柱体模型算例.
3.1 单一棱柱体模型算例异常体长宽高为4 km×4 km×4 km,顶面埋深2 km;异常体相对密度为2 g·cm-3,相对磁化率为0.25SI,电阻率为10 Ωm;异常体埋藏于相对密度为0 g·cm-3,相对磁化率为0 SI,电阻率为100 Ωm的均匀半空间,异常体的结构空间分布如图 8所示.
大地电磁采用全部波阻抗张量元素,反演的观测数据源自100、50、10、5、1、0.5、0.1、0.05和0.01 Hz9个频率正演计算加入3%随机高斯误差合成数据;重、磁反演的观测数据来自重磁三维理论正演计算加入3%随机高斯误差合成数据,对模型的合成数据分别进行了单独反演和联合反演,反演的初始模型采用均匀半空间模型,其物性参数为:0 g·cm-3、0 SI和100 Ωm.磁法正反演过程所假设的磁性参数为:地球磁场强度50000 nT与剖面夹角0°、剖面内磁倾角90°、地磁场倾角90°.
由以上三种方法的单一反演结果与联合反演结果的对比(如图 9所示),我们可以看出联合反演相比于单一反演,在数值的恢复上更能接近异常体的理论数值,而且在边界的分辨能力上也具有了较大的提升.
在第2.1节交叉梯度函数中详细地给出了非边界网格、边界面网格、边界线网格和边界点网格的tx、ty、tz离散计算方式,通过以上计算可以计算出每个物性单元交叉梯度值xyz三个方向的分量tx、ty、tz,我们这样定义交叉梯度值Φcg表达式:
(40) |
经过这样计算,我们可以得到交叉梯度值在空间的分布(如图 10所示).
为了进一步分析对比联合反演和单独反演的准确性,本文对单一反演和联合反演的最终反演结果进行了交叉梯度值的计算,计算结果如表 1所示,联合反演的交叉梯度值明显小于单独反演1~2个数量级,说明了联合反演结果的相似性程度高于单独反演的结果,从而可以得到更为一致性的模型进行更准确、更简单的地质解释.以下交叉梯度值在空间分布的切片深度与位置均与反演结果的切片深度和位置一致.
大地电磁停止迭代的条件是拟合差小于1或者大于反演最大迭代次数;重力和磁法的反演终止条件是反演的均方根误差达到收敛标准或者大于反演最大迭代次数.在本文的所有模型算例中,各方法的反演结果均在合理条件下获得的,反演拟合差迭代曲线如图 11所示.
根据联合反演和单独反演的模型结果,绘制了单独反演图和联合反演图中电阻率、磁化率和密度的物性交会图(如图 12所示).为了更加明显体现出联合反演的效果,我们选取单棱柱异常体,其中红色点代表异常体,蓝色点代表围岩;空间三维坐标轴:X轴-磁化率,Y轴-电阻率,Z轴-密度,各个点的三维坐标代表物性值.我们可以看到单独反演模型的物性参数交会图分散比较严重,而且物性值较真值具有一定差距;联合反演模型的物性参数交会图较为凝聚,数值的恢复程度较好.
为了测试复杂棱柱体模型是否会对联合反演结果产生虚假异常,我们设计了以下复杂棱柱体模型算例.三个棱柱体模型几何参数如下:棱柱体(左)长3 km,宽3 km,高3 km,顶面埋深2 km;棱柱体(中)长2.5 km,宽3 km,高3 km,顶面埋深1.5 km;棱柱体(右)长2 km,宽3 km,高2.5 km,顶面埋深1 km.异常体模型的物性参数如下:棱柱体(左)相对磁化率0.2 SI,相对密度2.5 g·cm-3,电阻率20 Ωm;棱柱体(中)相对磁化率0.2 SI,相对密度2.0 g·cm-3,电阻率10 Ωm;棱柱体(右)相对磁化率0.25 SI,相对密度2.5 g·cm-3,电阻率15 Ωm;围岩相对磁化率0.0 SI,相对密度0.0 g·cm-3,电阻率100 Ωm.异常体的空间结构分布如图 13所示.
大地电磁采用全部波阻抗张量元素,反演的观测数据源自100、50、10、5、1、0.5、0.1、0.05 Hz和0.01 Hz 9个频率正演计算加入3%随机高斯误差合成数据;重、磁反演的观测数据来自重磁三维理论正演计算加入3%随机高斯误差合成数据,对模型的合成数据分别进行了单独反演和联合反演,反演的初始模型采用均匀半空间模型,其物性参数为:0 g·cm-3、0 SI和100 Ωm.磁法正反演过程所假设的磁性参数为:地球磁场强度50000 nT与剖面夹角0°、剖面内磁倾角90°、地磁场倾角90°.组合模型的反演结果如图 14所示,交叉梯度值的空间分布如图 15所示,反演迭代曲线如图 16所示.
为了进一步分析对比联合反演和单独反演的准确性,本文对单一反演和联合反演的最终反演结果进行了交叉梯度值的计算,计算结果如表 2所示(括号内的数值表示为联合反演中每对物性参数之间的交叉梯度RMS值).
由于考虑到三维大地电磁反演速率较慢,花费时间较长,并且各个频率都要进行正演和“拟正演”求解过程互相独立,所以在正演和“拟正演”时利用MPI将其多个频率分配到不同的进程中进行计算,改善了三维大地电磁反演速度慢的问题.
对表 3进行分析,随着进程数的增加,并行加速比随之提高;当使用12个进程的时候,即每个进程计算1个频率,并行加速比有较大的提高;由于节点数增加时,进程之间的通信会占用一定的时间,导致并行效率有所下降.从表 3的统计结果来看,运用MPI并行加速相对于串行程序反演速率加速效果比较明显.
联合反演能够克服单一方法的局限性,减少地球物理反演的多解性,是今后地球物理发展的必然方向.本文实现了基于交叉梯度约束的重力、磁法和大地电磁三维联合反演算法,并设计了单棱柱体模型和组合棱柱体模型的合成数据对联合反演程序进行了测试,并将联合反演算法进行了MPI并行加速计算.
通过以上模型反演算例对比结果可以看出:
(1) 联合反演相对于单一反演不论在数值的恢复能力上还是在异常体边界的刻画能力上都具有较大程度上的提高.
(2) 在反演过程中,我们并不知道密度、磁化率和电阻率之间是否存在确定的关系式,但基于交叉梯度约束的联合反演依旧能够同时改善重力、磁法和大地电磁单独反演结果,这也进一步说明了交叉梯度联合反演的特点和优势:并不依赖于不同物性间的特定关系式,而是依靠地下结构的相似性进行约束反演,因而具有更加普遍的适用性.
(3) 联合反演程序中实现了MPI并行加速运算,很好地解决了反演计算速度慢的问题.
对实际重磁电资料应用中要注意的问题主要有:
(1) 网格的匹配问题.三种方法在实际工作采集数据时可能因为探测目标、精度、环境等因素的影响,使其在反演时候采用的反演网格有所差异,可以参考国外学者在处理网格兼容性问题的解决方案,利用多对一网格映射算法来解决.
(2) 交叉梯度权重因子的选取.需要借助模型试验、测区的地质资料和经验综合进行合理选取;如有该区钻探资料,则在权重因子的选取方面应考虑到钻探资料提供的先验信息.
(3) 计算效率.本文采用MPI并行计算对反演算法进行加速,对于实际重磁电资料的处理效率上往往还得取决于数据量、地下网格剖分等一系列因素.
致谢 感谢评审专家在百忙之中为本文的修改提出细致而宝贵的意见.
Colombo D, De Stefano M. 2007. Geophysical modeling via simultaneous joint inversion of seismic, gravity, and electromagnetic data:Application to prestack depth imaging. The Leading Edge, 26(3): 326-331. DOI:10.1190/1.2715057 |
Egbert G D, Kelbert A. 2012. Computational recipes for electromagnetic inverse problem. Geophysical Journal International, 189(1): 251-267. DOI:10.1111/j.1365-246X.2011.05347.x |
Gallardo L A, Meju M A. 2003. Characterization of heterogeneous near-surface materials by joint 2D inversion of DC resistivity and seismic data. Geophysical Research Letters, 30(13): 1658. DOI:10.1029/2003GL017370 |
Gao D Z, Zhao J H, Bo Y L, et al. 2004. A profile study of gravitative-magnetic and seismic comprehensive survey in the East China Sea. Chinese Journal of Geophysics (in Chinese), 47(5): 853-861. |
Guan Z N. 2005. Geomagnetic Field and Magnetic Exploration (in Chinese). Beijing: Geological Publishing House: 107-110.
|
Jegen M D, Hobbs R W, Tarits P, et al. 2009. Joint inversion of marine magnetotelluric and gravity data incorporating seismic constraints:Preliminary results of sub-basalt imaging off the Faroe Shelf. Earth and Planetary Science Letters, 282(1-4): 47-55. DOI:10.1016/j.epsl.2009.02.018 |
Jiang M, Peng M, Wang Y X, et al. 2012. Geophysical evidence for deep subduction of Indian Lithospheric plate beneath Eastern Himalayan Syntaxis. Acta Petrologica Sinica (in Chinese), 28(6): 1755-1764. |
Key K. 2016. MARE2DEM:a 2-D inversion code for controlled-source electromagnetic and magnetotelluric data. Geophysical Journal International, 207(1): 571-588. DOI:10.1093/gji/ggw290 |
Li D C, Yang S J, Hu Z Z, et al. 2012. Integrated interpretation of 3D gravity, magnetic, electromagnetic and seismic data; a case study of conglomerate mass investigation in piedmont area of Kuche Depression. Oil Geophysical Prospecting (in Chinese), 47(2): 353-359. |
Li Y G, Oldenburg D W. 1996. 3-D inversion of magnetic data. Geophysics, 61(2): 394-408. DOI:10.1190/1.1443968 |
Li Y G, Oldenburg D W. 2003. Fast inversion of large-scale magnetic data using wavelet transforms and a logarithmic barrier method. Geophysical Journal International, 152(2): 251-265. DOI:10.1046/j.1365-246X.2003.01766.x |
Liu J X, Tong X Z, Guo R W, et al. 2012. Magnetotelluric Sounding Method Exploration Data Processing, Inversion and Interpretation (in Chinese). Beijing: Science Press: 41-68.
|
Nocedal J, Wright S J. 2006. Numerical Optimization. 2nd ed. New York: Springer-Verlag: 135-160.
|
Peng M. 2012. Joint Inversion of Magnetotelluric and Teleseismic Data[Ph. D. Thesis] (in Chinese). Beijing: China University of Geosciences, Beijing.
|
Peng M, Tan H D, Jiang M, et al. 2013. Three-dimensional joint inversion of magnetotelluric and seismic travel time data with cross-gradient constraints. Chinese Journal of Geophysics (in Chinese), 56(8): 2728-2738. DOI:10.6038/cjg20130821 |
Rodi W, Mackie R L. 2001. Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion. Geophysics, 66(1): 174-187. DOI:10.1190/1.1444893 |
Sasaki Y. 1989. Two-dimensional joint inversion of magnetotelluric and dipole-dipole resistivity data. Geophysics, 54(2): 254-262. DOI:10.1190/1.1442649 |
Tan H D, Yu Q F, Booker J, et al. 2003. Magnetotelluric three-dimensional modeling using the staggered-grid finite difference method. Chinese Journal of Geophysics (in Chinese), 46(5): 705-711. |
Wang K P, Tan D H, Wang T. 2017. 2D joint inversion of CSAMT and magnetic data based on cross-gradient theory. Applied Geophysics, 14(2): 279-290. DOI:10.1007/s11770-017-0615-z |
Wang Y P. 2017. Research on Forward Modeling and Inversion of Tensor CSAMT in Three-Dimensional Axial Anisotropic Media[Ph. D. Thesis] (in Chinese). Beijing: China University of Geosciences, Beijing.
|
Zeng H L. 2005. Gravity Field and Gravity Exploration (in Chinese). Beijing: Geological Publishing House: 4-29.
|
Zhou L F. 2012. Two dimensional joint inversion of MT and seismic data[Master's thesis] (in Chinese). Beijing: China University of Geosciences, Beijing.
|
高德章, 赵金海, 薄玉玲, 等. 2004. 东海重磁地震综合探测剖面研究. 地球物理学报, 47(5): 853-861. DOI:10.3321/j.issn:0001-5733.2004.05.017 |
管志宁. 2005. 地磁场与磁力勘探. 北京: 地质出版社: 107-110.
|
姜枚, 彭淼, 王有学, 等. 2012. 喜马拉雅东构造结岩石圈板片深俯冲的地球物理证据. 岩石学报, 28(6): 1755-1764. |
李德春, 杨书江, 胡祖志, 等. 2012. 三维重磁电震资料的联合解释——以库车大北地区山前砾石层为例. 石油地球物理勘探, 47(2): 353-359. |
柳建新, 童孝忠, 郭荣文, 等. 2012. 大地电磁测深法勘探-资料处理、反演与解释. 北京: 科学出版社: 41-68.
|
彭淼. 2012.大地电磁与天然地震数据联合反演研究[博士论文].北京: 中国地质大学(北京).
|
彭淼, 谭捍东, 姜枚, 等. 2013. 基于交叉梯度耦合的大地电磁与地震走时资料三维联合反演. 地球物理学报, 56(8): 2728-2738. DOI:10.6038/cjg20130821 |
谭捍东, 余钦范, Booker J, 等. 2003. 大地电磁法三维交错采样有限差分数值模拟. 地球物理学报, 46(5): 705-711. DOI:10.3321/j.issn:0001-5733.2003.05.019 |
王堃鹏. 2017.张量CSAMT三维主轴各向异性正反演研究[博士论文].北京: 中国地质大学(北京).
|
曾华霖. 2005. 重力场与重力勘探. 北京: 地质出版社: 4-29.
|
周丽芬. 2012.大地电磁与地震数据二维联合反演研究[硕士论文].北京: 中国地质大学(北京).
|