地球物理学进展  2016, Vol. 31 Issue (4): 1796-1802   PDF    
2.5维海洋可控源电磁反演算法及影响参数研究
陈光源1,2, 杜立彬1,2, 景建恩3, 吴凯4, 邓明3, 陈凯3     
1. 山东省科学院海洋仪器仪表研究所, 青岛 266001
2. 山东省海洋环境监测技术重点实验室, 青岛 266001
3. 中国地质大学(北京), 北京 100083
4. 中国科学院电子学研究所, 北京 100190
摘要: 海洋可控源电磁法作为一种较新的地球物理勘探方法,在国外已被成功应用于海洋油气资源与天然气水合物的探测中,而我国在这方面的研究开始较晚.本文编写了基于非线性共轭梯度算法的海洋可控源电磁2.5维反演程序,结合理论数据对反演代码的正确性进行验证并讨论了反演参数对反演运算速度及效果的影响.研究结果表明,1)本文编写的2.5维反演程序正确可靠,计算结果与理论模型一致,2)在非线性共轭梯度反演计算中,正则化因子、预条件矩阵、线性搜索及共轭梯度更新因子等参数对实际的反演速度及精度都存在一定影响,可根据不同的勘探需求,调整反演参数以达到较好的反演结果.
关键词海洋可控源电磁法     非线性共轭梯度     正则化因子     预条件矩阵     线性搜索     梯度更新方向    
2.5D inversion algorithm and parameters of marine controlled-source electromagnetic method
CHEN Guang-yuan1,2 , DU Li-bin1,2 , JING jian-en3 , WU Kai4 , DENG Ming3 , CHEN Kai3     
1. Institute of Oceanographic Instrumentation, Shandong Academy of Sciences, Qingdao 266001, China
2. Shandong Provincial Key Laboratory of Ocean Environmental Monitoring Technology, Qingdao 266001, China
3. China University of Geosciences, Beijing 100083, China
4. Institute of Electronics, Chinese Academy of Sciences, Beijing 100190, China
Abstract: As a new geophysical exploration method, Marine controlled-source electromagnetic method has been intensively studied and successfully applied into marine hydration and gas hydrate resource detection abroad. This study in China began later than the foreign countries and short for effective inversion method for 2.5-dimensional problem. This paper written a 2.5-dimensional code cored by the nonlinear conjugate gradients algorithm and the forward method adopted the adaptive finite element method. The synthetic and measured data were applied to prove the code's validity and affectivity. The conclusion was drawn that the 2.5-dimensional inversion code was correct and effective and the inversion result is consist with other geophysical exploration method. This code has merits as high calculation rate and inversion resolution that could provide with a good way to study the 2-dimensional conductivity distribution under the seafloor.
Key words: marine CESM     nonlinear conjugate gradient     regularization parameters     preconditioner     line searches     search directions    
0 引 言

海洋可控源电磁法(MCSEM)的研究经历了几十年的发展,已经成为人们普遍采用的一种地球物理勘探方法,应用范围从大洋洋壳的电导率分布研究扩展到火山、热液系统、近海油气探测以及天然气水合物的研究中.随着海洋可控源电磁探测活动的增加,需要解决的实际地球物理问题也越来越复杂(Wannamaker et al.,198419861997).一维模拟研究虽然已经十分成熟,但是对实际地球物理模型的描述水平有限.而三维数值模拟虽然对实际的地球结构模拟水平较高(殷长春等,2014; 韩波等,2015; 杨军等,2015),但目前在算法及模拟方面的研究还不够充分,很多问题还处于广泛的研究和讨论之中.因此,二维模拟就成为了最接近实际的地球物理结构且最具有可实现性的选择.

在海洋可控源电磁法研究中,2.5维模拟具体表现为地球物理模型的二维性和场源发射的三维性,与二维场源发射模拟相比,能够进一步提高实际的探测效果.1971年,Coggon首先发表了利用有限元方法解决2.5维电磁问题的文章(Coggon,1971).1995年,Unsworth在有限元正演的基础上提出了一种子空间海洋可控源电磁的二维反演算法(Unsworth and Oldenburg,1995),并于2001年引入二维OCCAM算法,在汤加lau盆地得以成功应用.2008年Abubakar实现了一种更为快速可靠的二维反演算法(Abubakar et al.,2008).2011年之后,Kerry Key在自适应有限元(Key and Ovall,2011)的基础上发展了基于OCCAM的反演算法.而在我国,2.5维数值模拟起步于直流正演计算,周熙襄、钟本善(周熙襄等,1983)等分别利用有限差分、有限单元法实现了点源二维电阻率正演.随后许多学者也做了这方面的研究,探讨了如2.5维问题中傅氏变换的波数的选取、边值条件的选定以及反傅氏变换算法等内容(罗延钟和孟永良,1986;徐世浙,1988).

1 非线性共轭梯度反演方法 1.1 非线性共轭梯度算法

Fletcher和Reeves于1964年提出非线性共轭梯度(NLCG,Nonlinear Conjugate Gradient)算法,1987年由Tarantla将其引入地球物理非线性问题的求解中.根据Tikhonov和Arsenin提出的理论(Tikhonov and Arsenin,1977),本文采用正则化方法(Newman and Alumbaugh,2000Rodi and Mackie,2001)定义非线性共轭梯度的目标函数,公式为

(1)

其中正则化参数λ为正数,正定矩阵V-1表示误差向量e的方差,L定义为与三角网格剖分有关的二阶差分算子.利用公式(2)~(3),对关于m的目标函数最小化,得到最佳模型为(Fletcher and Reeves,1964).

(2)
(3)

公式(2)中的pl为模型空间的搜索方向,αl为该搜索方向上的最优步长.

搜索方向可通过公式(4)及(5)迭代求解.公式(5)中的左端第一项为最速下降方向,第二项用于修正最速下降方向.在线性共轭梯度求解中用于保证当前迭代方向与上一个迭代方向共轭,而在非线性共轭梯度求解中则不作此类要求,只需要共轭搜索方向与梯度满足公式(6)中给出的弱条件.公式为

(4)
(5)
(6)
1.2 反演程序正确性验证

为了验证本文编写的2.5维反演代码的正确性,利用Weitemeyer等(2010)的海洋二维高阻异常地电模型进行正演计算,模型及参数见图 1,设异常体沿Y方向无限延伸,接收机布置在x方向-8~8 km的海底面上,间距1 km.发射机在距离海底面50 m高度的海水中,从-8~8 km拖曳行进,点距500 m,发射频率0.25 Hz.在正演响应结果中添加5%的高斯噪声,采用本文中编写的2.5维反演程序及传统的OCCAM算法分别对上述模型进行2.5维反演计算,结果见图 2,图中彩色背景显示为海底介质的电阻率分布情况.将NLCG和OCCAM算法的反演结果进行对比可知,二者均能较好的分辨埋深1 km的高阻异常,NLCG的反演高阻异常体电阻率最大值为8.5 Ω·m,电阻率接近真实值,横向分辨率较好;而OCCAM反演得到的高阻异常的电阻率为 6.6 Ω·m,所得到的模型更平滑.

图 1 二维高阻地电模型 Figure 1 Twodimensional model with resistive anomaly
图 2 不同算法对二维高阻地电模型反演结果对比 Figure 2 Comparison of different algorithms on two dimensional high resistivity model in inversion results

因此,本文编写的海洋可控源2.5维NLCG反演程序能够正确反映海底地层中的高阻异常体;对比OCCAM算法,能够得到更接近真实异常值的反演电阻率分布.

2 反演程序影响因子讨论 2.1 共轭梯度方向更新因子

对于NLCG方向的更新因子β,前人各有不同的算法,比较常见的求解方式如公式(7)~(11)所示,分别由Hestenes和Stiefel于1952年,Fletcher和Reeves于1964,Daniel于1967年,Polak和Ribiere于1969年以及Dai和Yuan于1999年提出(Dai and Yuan,1996).不同的算子对反演计算结果的影响不同.公式为

(7)
(8)
(9)
(10)
(11)

为选取对海洋可控源电磁法非线性共轭梯度求解最有益的更新因子,建立如图 3所示二维浅层小范围高阻异常地电模型,利用公式(7)~(11)给出的五种非线性共轭梯度方向更新因子进行2.5维反演计算,对比分析不同更新因子的反演速度和效果.

图 3 二维浅层高阻地电模型 Figure 3 Twodimensional model with resistive shallow buried

采用上述五种迭代更新引子进行反演计算,均能快速收敛,但反演效果却存在一定差别.根据图 4中给出的浅层异常体的2.5维反演结果可知,利用βKDY因子计算得到的高阻最大电阻率值为5.65 Ω·m,最接近真实电阻率值,但在高阻异常右下方得出一个电阻率约为0.44 Ω·m的低阻假异常,影响对地下介质的真实反映;以βKFR作为更新因子的反演结果中,高阻异常体的最大电阻率较小,且存在较明显的低阻假异常;βKHSβKD对高阻异常的反演效果优于上述两种因子,但仍有假异常的存在,且βKD因子在求解过程中需要显式计算并存储海森矩阵,因而需要耗费大量的计算时间和存储空间;迭代更新因子为βKPRP时,高阻异常最大电阻率值为3.71 Ω·m,且在高阻异常下方没有低阻假异常显示,但对异常纵向边界分辨能力不强.

图 4 不同非线性共轭梯度更新因子对高阻异常反演效果影响 (a)更新因子为βDY K ;(b)更新因子为βFR K ;(c)更新因子为βHS K ;(d)更新因子为βDK;(e)更新因子为βPKRP . Figure 4 Inversion results of various update primer of the vector of conjugate gradient (a)Update primer is βDY K;(b)Update primer is βFKR; (c)Update primer is βHS K ;(d)Update primer is βDK ;(e)Update primer is βPKRP
2.2 线性搜索对反演结果的影响

对公式(2)的最小化需要通过迭代求解.由于该函数只包含一个未知量α,所以可以通过全局最小化方法来进行求解,即寻找关于α的ψ最小值(Moré and Thuente,1994).而高斯牛顿方法却无法区分函数的局部最小和全局最小.但是,全局最小化也会导致每次NLCG迭代都要进行很多次正演计算.为此,本文中不采用全局线性最小化的方法,而是对其进行一定的修正,在固定的共轭方向上进行线性搜索,是高斯牛顿方法的一元函数版本.定义关于α的一元函数Φl,其高斯牛顿近似值为

(12)

通过线性搜索,生成一组模型为

(13)

其中αl,0=0,

(14)

由于ψ是关于m的二次函数,Φl是关于α的二次函数,由上述公式的最小化过程,可以得到线性搜索的更新步长为

(15)

本文采取Rodi和Mackie(2001)中给出的对高斯牛顿公式的改正方式,为保证在线性搜索过程中得到最优模型,首先定义

(16)

在迭代过程中,如果Φl的值增大,即Φll,k)>Φll,k-1),则将线性搜索步长减半,步长更新公式变为

(17)

在经过线性搜索的第一次迭代后,如果当前模型与前一个模型之间存在使目标函数达到最小值的点,即满足Φl(αl,best)Φl(αl,k)<0,则通过三次样条插值得到更新的步长.当通过插值得到的目标函数的估计值在一定误差范围内与目标函数保持一致,则认为线性搜索收敛.在高斯牛顿迭代中,收敛条件一般如下公式(18)所示,其中τ<<1.如果在设定的最大搜索次数内,迭代无法收敛,或是Φl(αl,k+1)>1.5Φl(αl,best),则认为线性搜索失败.不论线性迭代成功与否,都以最终迭代的结果为最优模型,如公式(19)所示.如果线性搜索收敛,那么就通过公式(5)计算得到更新的共轭梯度方向pl+1,如果不收敛,就将公式中的第一项作为最速下降方向,此时新的搜索方向与上一次的搜索方向不满足共轭条件.公式为

(18)
(19)

为了说明线性搜索对NLCG反演的效果和作用,在图 5中给出在不进行线性搜索的情况下,NLCG方法的迭代收敛情况.从曲线图中可知,在不进行线性搜索的情况下,NLCG反演的RMS虽然仍然能够保证一直减小,即搜索方向保证向着目标函数最小的方向进行修正,但是收敛的速度极为缓慢,在图 5所示的迭代过程中,RMS从10.04降至9.1需要近40次的迭代反演,不但耗费机器内存,更加浪费时间.因此,在NLCG算法中加入线性搜索的部分是非常有必要的,可以加快迭代的收敛速度,节省内存空间及运算时间.

图 5 线性搜索对非线性共轭梯度反演收敛速度影响 Figure 5 The influence of linear search on the convergence rate of the nonlinear conjugate gradient inversion
2.3 预条件矩阵

公式(5)中,-Clgl是预条件的最速下降方向,用于最小化目标函数ψ对ml的方向偏导数.由于预条件矩阵对共轭梯度的求解有很大的影响,因此在选择预条件矩阵时,要注意因增加预条件因子导致的计算量增加以及其将梯度向量转化为线性搜索方向的有效性等问题(林昌洪等,2008Lin et al.,2009).本文中采用的预条件矩阵为

(20)

其中γl为基于均一介质的雅可比矩阵的常量,与l无关,使预条件矩阵尽可能的接近海森矩阵的逆,从而避免显式求解并存储海森矩阵.

针对图 3中所示浅层小范围高阻异常模型,分别采取单位矩阵和-Cl为预条件矩阵进行NLCG反演,对比这两种预条件矩阵,对2.5维非线性反演问题求解速度和精度的影响.二者迭代速度曲线如图 6a所示,其中红色曲线代表-Cl为预条件矩阵时的迭代收敛情况,经过三次迭代就可达到收敛,RMS小于1;而蓝色曲线代表的利用单位阵进行预条件处理的求解方法,则需要8次迭代才能达到相同迭代水平.同时,对比二者反演结果,如图 6b所示,采用预条件处理方法的程序经迭代收敛后,能够正确得到高阻异常显示,见图 4;而没有进行预条件处理的计算方法,虽然经过较多的迭代仍能够达到目标的拟合差,但是却无法反映高阻异常的所在,见图 6b.因此,预条件矩阵的应用不但能够大大提高NLCG的收敛速度,而且有助于二维高阻异常体的反演成像.

图 6 预条件矩阵对2. 5 维共轭梯度反演的影响 (a)预条件矩阵对反演计算加速效果;(b)高阻异常在无预条件矩阵约束下的反演结果 Figure 6 The impaction of precondition matrix on the 2. 5D NLCG inversion (a)Comparison of inversion with and without precondition matrix; (b)The inversion result of 2. 5D NLCG without precondition matrix.
2.4 正则化因子

正则化方法由Tihonov提出,用于平衡对目标函数最小化求解过程中,数据偏差度与模型粗糙度之间的关系.根据公式(9)给出的NLCG的目标函数可知,正则化因子数值越大,模型越平滑,对数据的拟合度越低;正则化因子值越小,模型越粗糙,对数据的拟合度越高.本文中采用使用的固定因子方法,设定几个正则化因子值,通过对同一理论模型进行反演试算,选取反演效果最好的λ,作为固定的正则化因子使用.采用图 1所示二维地电模型,对正演理论计算结果添加5%的高斯噪声,分别选择正则化因子为0.3、 3和30,进行2.5维的非线性共轭梯度反演,分析不同正则化因子取值对反演效果的影响,并选取使反演结果最佳的正则化因子.

图 7中给出了正则化因子不同取值时,NLCG反演的迭代收敛速度及RMS值.结合图 8中a~c对应的三个正则化因子下的反演结果可知,对于同一地电模型,其他条件不变的情况下,上述三种正则化因子的取值均能保证反演迭代的快速收敛,且正则化因子取值越小,迭代收敛的越快;正则化因子越大,迭代收敛的越慢.但是当正则化因子取值较小时,该算法的反演分辨率教差.正则化因子取0.3时,最大高阻异常值达到31.1 Ω·m,已经无法正确反映海底地层下高阻异常的形态,且存在较多的假异常.正则化因子取3时,最大高阻异常值约为20.5 Ω·m,最接近正演模型电阻率值,但是对异常体形态的反映不够理想;正则化因子取30时,最大高阻异常值为8.5Ω·m,能够正确反映海底高阻异常体的形态及其在海底地层中所处位置.

图 7 不同正则化因子对非线性共轭梯度反演迭代速度影响 Figure 7 The inversion speed of 2. 5D NLCG with different regularization parameters
图 8 不同正则化因子对非线性共轭梯度反演结果影响 (a)正则化因子为30;(b)正则化因子为3;(c)正则化因子为0. 3. Figure 8 The inversion results of 2. 5D NLCG with various regularization factors (a)Regularization parameter is 30;(b)Regularization parameter is 3;(c)Regularization parameter is 0. 3.
3 结 论

3.1 本文编写了基于非线性共轭梯度算法的海洋可控源电磁2.5维反演程序,实现了海洋可控源电磁2.5维反演模拟计算,经理论数据证明,该反演程序正确有效.同时,建立针对海底浅层高阻异常的地电模型,依次讨论梯度方向更新因子、线性搜索、预条件矩阵及正则化因子等反演参数对海洋可控源电磁2.5维非线性共轭梯度反演算法的计算精度和速度的影响,并得出如下结论:

1) 梯度方向更新因子:对于相同的地电模型,βKDYβKFR的反演结果中,假异常较明显;βKD因子在求解过程中,需要显示计算并存储海森矩阵,对计算时间和存储空间的需求较大;βKPRP因子对异常体的纵向边界分辨能力有限;βKHS对高阻的横向分辨能力较差,对水平方向异常区域的反映大于正演模型.

2) 线性搜索:在非线性共轭梯度算法中加入线性搜索的部分是非常有必要的,可以加快迭代的收敛速度,节省内存空间及运算时间.

3) 正则化因子:在NLCG反演计算中,正则化因子数值越大,模型越平滑,但模型对数据的拟合度越低.否则,反之.

4) 预条件矩阵:在NLCG的反演计算中增加预条件矩阵,能够有效提高反演的收敛速度及成像精度.

3.2 然而,在本文的研究中,所使用的理论模型相对简单,对复杂模型的讨论稍显不足.通过对理论模型响应数据的反演计算发现,非线性共轭梯度算法在处理复杂形态异常和带地形的反演问题时,迭代收敛速度会有所下降.在后续的研究中,建议增加此类模型的探讨研究并增加实测数据处理实例.另外,考虑二维模型对实际的三维结构存在近似误差,计划在以后的工作中,开展三维海洋可控源电磁法非线性共轭梯度反演研究.

致谢 感谢山东省科学院海洋仪器仪表研究所对本研究给予的支持,感谢中国地质大学(北京)地下信息探测技术与仪器教育部重点实验室提供的良好平台,感谢邓明教授及景建恩老师的帮助和学术指导.
参考文献
[1] Abubakar A, Habashy T M, Druskin V L, et al.2008. 2.5D forward and inverse modeling for interpreting low-frequency electromagnetic measurements[J]. Geophysics, 73 (4) : F165–F177. DOI:10.1190/1.2937466
[2] Coggon J H.1971. Electromagnetic and electrical modeling by the finite element method[J]. Geophysics, 36 (1) : 132–155. DOI:10.1190/1.1440151
[3] Dai Y H, Yuan Y.1996. Convergence properties of the Fletcher-Reeves method[J]. IMA Journal of Numerical Analysis, 16 (2) : 155–164. DOI:10.1093/imanum/16.2.155
[4] Fletcher R, Reeves C M.1964. Function minimization by conjugate gradients[J]. The Computer Journal, 7 (2) : 149–154. DOI:10.1093/comjnl/7.2.149
[5] Han B, Hu X Y, Schultz A, et al.2015. Three-dimensional forward modeling of the marine controlled-source electromagnetic field with complex source geometries[J]. Chinese Journal of Geophysics (in Chinese), 58 (3) : 1059–1071. DOI:10.6038/cjg20150330
[6] Key K, Ovall J.2011. A parallel goal-oriented adaptive finite element method for 2[J]. 5-D electromagnetic modelling[J]. Geophysical Journal International, 186 (1) : 137–154. DOI:10.1111/j.1365-246X.2011.05025.X
[7] Lin C H, Tan H D, Tong T.2009. Parallel rapid relaxation inversion of 3D magnetotelluric data[J]. Applied Geophysics, 6 (1) : 77–83. DOI:10.1007/s11770-009-0010-5
[8] Luo Y Z, Meng Y L.1986. Some problems on resistivity modeling for two-dimensional structures by the finite element method[J]. Acta Geophysica Sinica (in Chinese), 29 (6) : 613–621.
[9] Moré J J, Thuente D J.1994. On line search algorithms with guaranteed decrease[J]. ACM Transactions on Mathematical Software, 20 (3) : 286–307. DOI:10.1145/192115.192132
[10] Newman G A, Alumbaugh D L.2000. Three-dimensional magnetotelluric inversion using non-linear conjugate gradients[J]. Geophysical Journal International, 140 (2) : 410–424. DOI:10.1046/j.1365-246x.2000.00007.x
[11] Rodi W, Mackie R L.2001. Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion[J]. Geophysics, 66 (1) : 174–187. DOI:10.1190/1.1444893
[12] Tikhonov A N, Arsenin V Y.1977. Solutions of Ill-Posed Problems. Washington: V. H. Winston and Sons[M]. .
[13] Unsworth M, Oldenburg D.1995. Subspace inversion of electromagnetic data: Application to mid-ocean-ridge exploration[J]. Geophysical Journal International, 123 (1) : 161–168. DOI:10.1111/j.1365-246x.1995.tb06668.x
[14] Wannamaker P E, Hohmann G W, SanFilipo W A.1984. Electromagnetic modeling of three-dimensional bodies in layered earths using integral equations[J]. Geophysics, 49 (1) : 60–74. DOI:10.1190/1.1441562
[15] Wannamaker P E, Stodt J A, Rijo L.1986. Two-dimensional topographic responses in magnetotellurics modeled using finite elements[J]. Geophysics, 51 (11) : 2131–2144. DOI:10.1190/1.1442065
[16] Wannamaker P E, Johnston J J, Stodt J A, et al.1997. Anatomy of the southern Cordilleran hingeline, Utah and Nevada, from deep electrical resistivity profiling[J]. Geophysics, 62 (4) : 1069–1086. DOI:10.1190/1.1444208
[17] Weitemeyer K, Gao G Z, Constable S, et al.2010. The practical application of 2D inversion to marine controlled-source electromagnetic data[J]. Geophysics, 75 (6) . DOI:10.1190/1.3506004
[18] Yang J, Liu Y, Wu X P.2015. 3D simulation of marine CSEM using vector finite element method on unstructured grids[J]. Chinese Journal of Geophysics (in Chinese), 58 (8) : 2827–2838. DOI:10.6038/cjg20150817
[19] Yin C C, Ben F, Liu Y H, et al.2014. MCSEM 3D modeling for arbitrarily anisotropic media[J]. Chinese Journal of Geophysics (in Chinese), 57 (12) : 4110–4122. DOI:10.6038/cjg20141222
[20] 韩波, 胡祥云, SchultzA, 等.2015. 复杂场源形态的海洋可控源电磁三维正演[J]. 地球物理学报, 58 (3) : 1059–1071. DOI:10.6038/cjg20150330
[21] 林昌洪, 谭捍东, 佟拓. 2008. 大地电磁三维共轭梯度反演研究[C]. //中国地球物理学会第二十四届年会论文集. 北京: 中国大地出版社, 526.
[22] 罗延钟, 孟永良.1986. 关于用有限单元法对二维构造作电阻率法模拟的几个问题[J]. 地球物理学报, 29 (6) : 613–621.
[23] 徐世浙.1988. 点源二维各向异性地电断面的直流电场有限元解法[J]. 山东海洋学院学报, 18 (1) : 81–90.
[24] 杨军, 刘颖, 吴小平.2015. 海洋可控源电磁三维非结构矢量有限元数值模拟[J]. 地球物理学报, 58 (8) : 2827–2838. DOI:10.6038/cjg20150817
[25] 殷长春, 贲放, 刘云鹤, 等.2014. 三维任意各向异性介质中海洋可控源电磁法正演研究[J]. 地球物理学报, 57 (12) : 4110–4112. DOI:10.6038/cjg20141222
[26] 周熙襄, 钟本善, 江玉乐.1983. 点源二维电阻率法有限差分法正演计算[J]. 物探化探计算技术 (2) : 1–9.