2019, Vol. 62
随着陆地油气资源大规模勘探开采，资源量呈现逐渐衰减趋势.为了满足国家经济发展的迫切需求，人们将探索目标转向了海洋.海洋油气资源储量丰富，具有很好的开发潜力，但总体资源储量，尤其是深水区的油气资源分布尚未探明.因此，研究高精度和高效的深海油气资源勘探方法具有重要意义.海洋可控源电磁方法(Marine Controlled-Source Electromagnetic Method, MCSEM)(Constable, 2010; 殷长春等，2014)是一种钻前储存预测的有效方法.该方法通过船载拖曳发射(多为水平偶极子源)电磁信号，利用海底阵列接收机采集电磁信号，进而分析海底电性结构，预测油气靶区并分析其空间赋存特征.MCSEM通常采用的工作频率为0.01~10 Hz.统计表明，海洋可控源电磁技术可有效地提高海洋油气勘探的成功率，已经成为继海洋地震之后有效的海洋油气勘探技术(Hesthammer et al., 2010).
迄今为止，海洋可控源一、二维正反演理论体系已经很成熟，并已在实际生产工作中得到应用.目前，国际上可提供MCSEM数据采集服务的仅有挪威的EMGS公司.该公司自2002年成立至今已完成从北极到澳大利亚30~3500 m海水深度范围内700多次海底勘查任务，并且可提供MCSEM三维数据反演解释.
随着在实测资料数据处理过程中对精细程度的要求不断提高，一、二维MCSEM正反演理论体系已经不能满足逐步提高的实际工作要求，而电磁法三维正演理论体系逐步趋向于成熟，因此开展复杂地形条件下MCSEM三维反演方法研究工作已成为目前研究的总体趋势.Zhdanov等(2007)基于严格的积分方程形式，对正则化梯度类型反演算法进行了改进，随后实现了利用积分方程算子的重加权正则化共轭梯度算法并将其成功的应用于海洋电磁数据三维反演当中(Zhdanov et al., 2014).Jing等(2008)测试了各向异性成像对于MCSEM数据解释的重要性.Zach等(2008)实现了快速有限差分时域模拟(FDTD)，能够采用相对较少的迭代得到较好的反演结果.Newman等(2010)基于非线性共轭梯度和全波场数值模拟，提出了可处理横向异性问题的三维MCSEM成像算法并验证了其有效性.Schwarzbach等(2011)提出了几种并行计算方案，能够提高计算效率加速模拟过程.Sasaki(2013)基于高斯牛顿算法，采用交错网格有限差分技术进行正演计算，完成了对CSEM和MT数据的分别反演和联合反演算法.Amaya等(2016)提出了用于反演三维MCESM数据的Hessian近似方法，降低了数值模拟的复杂程度.贲放等(2016)提出了基于L-BFGS最优化算法的MCSEM三维反演，降低了传统方法对于计算机内存的要求.彭荣华(2016)实现了针对海洋和陆地不同CSEM勘探环境的三维并行反演程序并验证了其有效性.
由于非结构网格在建模时具有较强的灵活性，能够精确地拟合复杂海底地形和地下异常体，因此基于非结构网格剖分的海洋电磁正反演已经日益成为近几年研究的热点.在正演计算方面，杨军等(2015)实现了MCSEM三维矢量有限元正演；叶益信等(2016)采用体积约束对异常区域进行了局部加密，实现了MCSEM三维有限元正演；陈汉波等(2017)基于非结构网格实现了电导率呈各向异性条件下MCSEM三维有限元正演.在反演方面，Schwarzbach和Haber(2013)采用基于非结构四面体网格有限元技术进行正演计算，实现了CSEM动态正则化反演方案并与传统标准正则化反演方案进行了比较；Wang等(2018)基于非结构四面体网格有限单元法和高斯-牛顿反演技术实现了MCSEM数据三维各向异性反演.
相较于国外研究，国内MCSEM技术研究起步较晚，目前尚未实现效果较好、适用于复杂起伏海底地形的三维电磁数据反演方法.为此，本文开展基于非结构网格剖分的频率域海洋可控源电磁三维反演研究.我们首先介绍基于非结构网格剖分MCSEM三维矢量有限元正演理论并进行精度验证，然后提出本文反演所适用的目标函数、模型粗糙度计算方法和反演参数上下限约束技术，最后利用理论算例验证反演程序的正确性与实用性.
1 正演理论 1.1 正演方程设时谐因子为eiωt，由麦克斯韦方程微分形式可推导出三维电磁正演双旋度电场方程：
|
(1) |
式中，E为电场强度，ω=2πf是角频率，σ为电导率，Js表示发射源项，i是虚单位，ε为介电常数，μ为磁导率.通常假设真空中介电常数ε0=8.854×10-12F·m-1和磁导率μ0=4π×10-7H·m-1.
为了保证电磁场存在唯一解，假设模拟区域的外边界Γ距离发射源足够远，电场强度满足第一类狄里克雷边界条件：
|
(2) |
式中，n表示边界上的单位法向量.根据加权余量法，使权函数W和残差R的内积在计算区域Ω上的积分最小，并利用分部积分展开(1)式的第一项最终得到
|
(3) |
式中，
我们假设伽辽金法中权函数W选用矢量基函数N，并将单元内任意点电场表示为矢量基函数与棱边电场的内积，则依据广义变分原理对每个单元离散可得
|
(4) |
将(4)式表示成矩阵的形式，即
|
(5) |
其中，各单元矩阵的表达式为
|
(6) |
|
(7) |
|
(8) |
|
(9) |
对于发射源与网格单元相交的情况，我们将发射源位于四面体单元内的部分设置为一个电偶极子，并假设其位于该单元两个交点的中心.由此，利用脉冲函数的特性，上式中左侧源项可以化简为Se(i)=Nie(x0, y0, z0)·Jse，其中(x0, y0, z0)为上述两个交点中心的坐标.对每个单元的单元矩阵进行组装(Jahandari et al., 2014; 齐彦福等, 2017)，并施加边界条件可得到大型稀疏对称方程组：
|
(10) |
使用直接求解器MUMPS (Amestoy et al., 2006)求解式(10)，可得四面体单元棱边上的电场值，从而根据矢量插值函数插值得到电场，并依据法拉第定律计算磁场响应.
1.3 正演精度验证为了验证正演模拟精度，我们针对一维海底含高阻层的海洋电磁模型，分别利用一维半解析和本文三维正演方法计算了电场分量Ex.从图 1给出的对比结果可以看出，两种方法计算结果吻合很好，振幅最大误差小于4%，相位最大误差小于0.7°，说明本文算法能够满足实际计算的精度要求，可为后续反演工作提供可靠的数据.
|
图 1 MCSEM非结构有限元正演模拟精度验证 (a) Ex分量振幅结果对比; (b) Ex分量振幅相对误差; (c) Ex分量相位结果对比; (b) Ex分量相位相对误差. Fig. 1 Accuracy check for MCSEM using FEM method (a) Comparison of amplitude of Ex component; (b) Relative error of amplitude; (c) Comparison of phase of Ex component; (d) Relative error of phase. |
本文建立如下频率域海洋可控源电磁三维反演目标函数：
|
(11) |
式中，m为M维电导率模型参数矢量，d为Nd维观测数据矢量.Cd为对角矩阵，其元素由数据误差的协方差构成，f(m)为通过矢量有限元方法计算的正演响应.m0是先验或初始模型参数，λ为正则化因子，Cm为模型粗糙度算子，表示了梯度或更高阶差分算子的离散形式.通过将公式(11)两侧对模型参数求导可得目标函数的梯度为
|
(12) |
式中，r=d-f(mn).利用泰勒展开并忽略高阶项可得高斯-牛顿反演方程：
|
(13) |
利用共轭梯度法求解方程(13)获得模型更新量Δm，然后可由mn+1=mn+sΔm获得新的模型, 其中s由Haber(2014)提出的线性搜索方法进行选择.
2.2 灵敏度计算前文正演理论中推导的海洋电磁三维正演方程组为
|
(14) |
式中，K为大型稀疏矩阵，E为待求解电场分量，S为源项.正演数据d可由d=LE获得，其中L为空间插值算子.数据关于第k个四面体网格电导率的灵敏度可写为
|
(15) |
定义如下N×M，N=Nedge×nf维矩阵
|
(16) |
由于本文正演计算中采用总场方法，所以
|
(17) |
由公式(13)可知，求解反演方程时需要计算灵敏度矩阵与向量的乘积和灵敏度矩阵的转置与向量的乘积.为此，假设与灵敏度矩阵和其转置相乘的向量分别为v和u.将两者分别与(17)式中两个方程的左右端相乘，并将右端部分乘积结果记为b和l，则有
|
(18) |
进而，将矩阵K移到方程左边并利用其对称性质可得到如下伴随正演方程：
|
(19) |
对(19)式中的伴随正演进行求解方程后，可得灵敏度矩阵与向量的乘积及灵敏度矩阵的转置与向量的乘积：
|
(20) |
实际利用直接法求解正演问题时，我们可以将矩阵K的分解结果保存，在伴随正演过程中直接调用可大幅度减少计算量.
2.3 非结构网格的模型粗糙度计算常用的模型粗糙度算子定义为模型梯度的L2范数(Parker, 1994; Zhdanov, 2002)
|
(21) |
上式对于非结构网格往往难以直接进行积分离散.为此，本文提出采用加权平方和方法(weighted sum of squares)，其表达式为
|
(22) |
式中各参数表示如下：
|
(23) |
|
(24) |
|
(25) |
其中，Vi是第i个四面体中的体积，N(i)是与四面体i共顶点的所有四面体网格个数.公式(22)中右侧方括号表示了梯度的L2范数.单元之间的距离Δrij是由各单元的质心进行计算，如图 2所示.
|
图 2 非结构网格模型粗糙度计算方法(参考Key, 2016) 蓝色圆点为四面体质心，白色虚线为质心间距离，黄色四面体为第i个四面体，绿色四面体为与第i个参数四面体共顶点的其他四面体. Fig. 2 Roughness calculation for unstructured grids (refer to Key, 2016). The blue dots are the tetrahedral centroid; the white dashed lines denote the distance between centroids; the yellow tetrahedron is the ith tetrahedron; the green tetrahedron is the rest tetrahedrons that share the vertices with the ith tetrahedron. |
反演过程中施加电导率上下限约束能够减少虚假异常，使反演结果更加接近真实情况.本文基于对数参数的电导率转换函数将模型参数变化到转换域内(称为对数变换)，对反演结果进行电导率上下限约束.首先，我们将模型变换到转换域内，即
|
(26) |
则关于xk的灵敏度矩阵为
|
(27) |
其中，
|
(28) |
公式(26)—(28)中，mk为第k个四面体网格的模型参数，ak和bk分别为模型参数的下界和上界，xk为转换域内对应的第k个四面体网格的模型参数.通过求解高斯-牛顿反演方程获得xk的改变量后，真实模型的改变量为
|
(29) |
为验证本文提出的海洋可控源电磁三维反演方法的有效性，我们分别设计了单测线水平海底模型和面积性起伏海底模型.我们首先讨论单测线水平海底模型.
3.1 单测线水平海底面模型如图 3所示，计算中心精细剖分区域为8000 m×4000 m×2000 m，模型总体大小为40000 m×40000 m×40000 m.为了对比和分析原始模型与反演结果，图中仅显示部分海水层与地下含异常体区域，显示区域大小为8000 m×4000 m×3300 m.空气层厚度为4000 m，电导率为1.0×10-8S·m-1.海水深度为1000 m，电导率为3.33 S·m-1.异常体中心距海面2000 m，大小为3000 m×2000 m×200 m，电导率为0.02 S·m-1.海底电导率为1.43 S·m-1.另外，我们设计发射源距海底50 m，共7个，相距1000 m，沿单条测线分布.发射频率为0.25 Hz.测点位于海底，与发射源同线分布，相距100 m，共81个.正演模型共包含452145个四面体单元，反演模型包含424987个四面体单元.我们使用电场Ex分量的实部和虚部进行三维反演，并在反演数据中加入了5%的高斯噪声.为了使反演更加稳定，根据对数变换将反演电导率的最大值和最小值限制在0.005 S·m-1到2 S·m-1范围之内.在每次反演迭代中，正则化参数λ减少到先前值的0.8倍.所有的计算都在配有Inter(R) Xeon(R) 2.8 GHz CPU的工作站上进行，其中CPU共有8个内核，内存为32 GB，总运行时间为22.6 h.
|
图 3 水平海底海洋可控源电磁正反演模型 (a)正演模型剖分; (b)反演模型剖分. Fig. 3 Forward and inverse model of a flat seafloor (a) Forward model; (b) Inverse model. |
图 4给出真实模型与添加5%高斯噪声后数据反演结果.从图可以看出，本文的反演方法较好地恢复了异常体的位置和电阻率.在x、y和z三个方向上，反演模型与真实模型均比较接近.反演过程中各参数的变化特征如图 5所示.从图可以看出，RMS和目标函数值快速下降，模型改变量最终也趋于稳定，证明本文方法收敛性很好.
|
图 4 水平海底真实模型与三维反演结果 (a)—(c)分别为真实模型沿y方向、x方向和z方向切片; (d)—(f)分别为相同方向反演模型切片. Fig. 4 Theoretical and recovered model of a flat seafloor (a)—(c) are slices of the true model respectively in y-, x-, z-directions; (d)—(f) are slices of the recovered model in the same direction. |
|
图 5 水平海底模型反演参数曲线 (a)数据拟合差(RMS); (b)目标函数值; (c)模型改变量的模; (d)正则化参数. Fig. 5 Inversion parameters vs. iterations (a) RMS; (b) Object functional; (c) Model variation norm; (d) Trade-off parameter. |
图 6给出水平海底模型反演结果与实际数据曲线对比.由于受随机噪声的影响，数据中量级较小的虚部产生部分细微波动，然而总体曲线较为光滑.通过曲线拟合结果可知反演结果与实际数据具有较好的一致性，验证了本文方法的有效性.
|
图 6 水平海底模型反演结果与实际数据对比 (a)实部曲线对比; (b)虚部曲线对比. Fig. 6 Comparison of predicted and observed data for the flat seafloor model. (a) Real part of Ex; (b) Imaginary part of Ex. |
考虑到实际海底往往起伏不平，本节设计如图 7所示的包含起伏海底地形的海洋电磁模型.地形起伏不对称，发射源距海底面50 m，分别沿3条相距1000 m的测线布设.发射源间距为1000 m，共21个发射源.测点位于海底，共243个，分别沿与源相对应的3条测线，间距为100 m，其余参数与图 3中模型一相同.
|
图 7 起伏海底海洋电磁模型 (a)海底地形(红色和蓝色点分别表示发射源和接收点位置); (b)正演模型剖分; (c)反演模型剖分. Fig. 7 Forward and inverse model for a topographic seafloor (a) Seafloor topography (Red and blue dots are respetively for transmitters and receivers.); (b) Forward model; (c) Inverse model. |
本算例中正演模型离散为1219644个单元(图 7b)，反演模型离散为1081795个单元(图 7c).同样采用与上述相同的计算条件，本模型的反演总耗时为73.93 h.图 8给出模型数据添加5%高斯噪声后的反演结果.比较图中的真实模型和反演结果可以发现，尽管存在复杂的海底地形，高阻储层的位置和电阻率得到了很好的恢复.储层的顶部和底部边界清晰，厚度正确，表明反演算法具有一定的实用性，可用于复杂海底地形条件下MCSEM数据反演.图 9给出反演参数随迭代次数的变化特征.从图可以看出，迭代过程中数据拟合差RMS和目标函数缓慢地下降，RMS接近目标值1，说明反演算法稳定性.
|
图 8 起伏海底真实模型与三维反演结果 (a)—(c)分别为真实模型沿y方向、x方向和z方向切片; (d)—(f)分别为相同方向反演模型切片. Fig. 8 Theoretical and recovered model of a topographic seafloor (a)—(c) are slices of the true model respectively in y-, x-, z-directions; (d)—(f) are slices of the recovered model in the same direction. |
|
图 9 起伏海底模型反演参数曲线 (a)数据拟合差RMS; (b)目标函数值; (c)模型改变量的模; (d)正则化参数. Fig. 9 Inversion parameter vs. iterations (a) RMS; (b) Object functional; (c) Model variation norm; (d) Trade-off parameter. |
图 10给出起伏海底模型反演结果与实际数据的对比.从图可以看出，较之于拟合程度较好的实分量，虚分量曲线受到地形起伏和噪声污染的影响稍显严重.这是由于本文反演中数据噪声是基于Ex分量的幅值添加的，因此量级较小的虚部影响较大.综合图 9和10可以看出，虽然数据中存在的噪声在一定程度上影响数据拟合，但反演结果很好地恢复了起伏海底下的高阻异常体特征.
|
图 10 起伏海底模型反演结果与实际数据对比 (a)实部曲线对比; (b)虚部曲线对比. Fig. 10 Inversion results and observed data (a) Real part of Ex; (b) Imaginary part of Ex. |
本文基于非结构网格成功实现频率域海洋可控源电磁(MCSEM)数据三维正则化反演，并通过理论数据测试了反演效果.数据试验表明，无论是对水平还是起伏海底地形，本文提出的三维反演方法均能够有效地恢复海底高阻储油层的位置和电阻率信息，本文的反演算法稳定可靠.由于本文使用直接求解技术实现正演问题求解，反演模拟中我们可重复使用正演矩阵分解结果，大大提高了反演计算效率.另外，从数据的拟合情况可以看出，虽然虚部数据受噪声影响稍大，然而由于其量级比实部小，反演模型仍能很好地恢复海底异常体的分布特征.本文中我们仅展示初步研究成果，未来将在提高反演效率的同时，开展利用其他地球物理资料进行层位约束反演及考虑各向异性模型的反演研究.
Amaya M, Morten J P, Boman L. 2016. A low-rank approximation for large-scale 3D controlled-source electromagnetic Gauss-Newton inversion. Geophysics, 81(3): E211-E225. DOI:10.1190/geo2015-0079.1 |
Amestoy P R, Guermouche A, L'Excellent J Y, et al. 2006. Hybrid scheduling for the parallel solution of linear systems. Parallel Computing, 2006, 32(2): 136-156. http://www.sciencedirect.com/science/article/pii/S0167819105001328
|
Ben F. 2016. 3D Marine Controlled-source Electromagnetic Modeling and Inversion Theory and Research on the Mechanism of Anisotropic Effect.[Ph. D. thesis] (in Chinese). Changchun: Jilin University. http://cdmd.cnki.com.cn/Article/CDMD-10183-1016090176.htm
|
Chen H B, Li T L, Xiong B, et al. 2017. Finite-element modeling of 3D MCSEM in arbitrarily anisotropic medium using potentials on unstructured grids. Chinese Journal of Geophysics (in Chinese), 60(8): 3254-3263. DOI:10.6038/cjg20170830 |
Constable S C. 2010. Ten years of marine CSEM for hydrocarbon exploration. Geophysics, 75(5): 75A67-75A81. DOI:10.1190/1.3483451 |
Haber E. 2014. Computational Methods in Geophysical Electromagnetics. SIAM-Society for Industrial and Applied Mathematics.
|
Hesthammer J, Stefatos A, Boulaenko M, et al. 2010. CSEM Technology as a value driver for hydrocarbon exploration. Marine and Petroleum Geology, 27(9): 1872-1884. DOI:10.1016/j.marpetgeo.2010.08.001 |
Jahandari H, Farquharson C G. 2014. A finite-volume solution to the geophysical electromagnetic forward problem using unstructured grids. Geophysics, 79(6): E287-E302. DOI:10.1190/geo2013-0312.1 |
Jing C, Green K, Willen D. 2008. CSEM inversion: impact of anisotropy, data coverage, and initial models.//78th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 604-608.
|
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 |
Newman G A, Commer M, Carazzone J J. 2010. Imaging CSEM data in the presence of electrical anisotropy. Geophysics, 75(2): F51-F61. DOI:10.1190/1.3295883 |
Parker R L. 1994. Geophysical Inverse Theory. Princeton: Princeton University Press.
|
Peng R H. 2016. Three-dimensional forward modeling and inversion of frequency-domain CSEM data[Ph. D. thesis] (in Chinese). Wuhan: China University of Geosciences.
|
Qi Y F, Yin C C, Liu Y H, et al. 2017. 3D time-domain airborne EM full-wave forward modeling based on instantaneous current pulse. Chinese Journal of Geophysics (in Chinese), 60(1): 369-382. DOI:10.6038/cjg20170131 |
Sasaki Y. 2013. 3D inversion of marine CSEM and MT data:An approach to shallow-water problem. Geophysics, 78(1): E59-E65. DOI:10.1190/geo2012-0094.1 |
Schwarzbach C, Börner R U, Spitzer K. 2011. Three-dimensional adaptive higher order finite element simulation for geo-electromagnetics-a marine CSEM example. Geophysical Journal International, 187(1): 63-74. DOI:10.1111/j.1365-246X.2011.05127.x |
Schwarzbach C, Haber E. 2013. Finite element based inversion for time-harmonic electromagnetic problems. Geophysical Journal International, 193(2): 615-634. DOI:10.1093/gji/ggt006 |
Wang F Y, Morten J P, Spitzer K. 2018. Anisotropic three-dimensional inversion of CSEM data using finite-element techniques on unstructured grids. Geophysical Journal International, 213(2): 1056-1072. DOI:10.1093/gji/ggy029 |
Yang J, Liu Y, Wu X P. 2015. 3D simulation of marine CSEM using vector finite element method on unstructured grids. Chinese Journal of Geophysics (in Chinese), 58(8): 2827-2838. DOI:10.6038/cjg20150817 |
Ye Y X, Li Y G, Liu Y, et al. 2016. 3D finite element modeling of marine controlled-source electromagnetic fields using locally refined unstructured meshes. Chinese Journal of Geophysics (in Chinese), 59(12): 4747-4758. DOI:10.6038/cjg20161233 |
Yin C C, Ben F, Liu Y H, et al. 2014. MCSEM 3D modeling for arbitrarily anisotropic media. Chinese Journal of Geophysics (in Chinese), 57(12): 4110-4122. DOI:10.6038/cjg20141222 |
Zach J J, Bjørke A K, Støren T, et al. 2008. 3D inversion of marine CSEM data using a fast finite-difference time-domain forward code and approximate Hessian-based optimization.//78th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 614-618.
|
Zhdanov M S. 2002. Geophysical Inverse Theory and Regularization Problems. New York, NY, USA: Elsevier.
|
Zhdanov M S, Gribenko A, Čuma M. 2007. Regularized focusing inversion of marine CSEM data using minimum vertical support stabilizer.//77th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 579-583.
|
Zhdanov M S, Endo M, Cox L H, et al. 2014. Three-dimensional inversion of towed streamer electromagnetic data. Geophysical Prospecting, 62(3): 552-572. DOI:10.1111/1365-2478.12097 |
贲放. 2016.海洋可控源电磁法三维正反演理论与各向异性影响机制研究[博士论文].长春: 吉林大学. http://cdmd.cnki.com.cn/Article/CDMD-10183-1016090176.htm
|
陈汉波, 李桐林, 熊斌, 等. 2017. 基于Coulomb规范势的电导率呈任意各向异性海洋可控源电磁三维非结构化有限元数值模拟. 地球物理学报, 60(8): 3254-3263. DOI:10.6038/cjg20170830 |
彭荣华. 2016.频率域可控源电磁法三维正反演研究[博士论文].武汉: 中国地质大学. http://cdmd.cnki.com.cn/Article/CDMD-10183-1016083550.htm
|
齐彦福, 殷长春, 刘云鹤, 等. 2017. 基于瞬时电流脉冲的三维时间域航空电磁全波形正演模拟. 地球物理学报, 60(1): 369-382. DOI:10.6038/cjg20170131 |
杨军, 刘颖, 吴小平. 2015. 海洋可控源电磁三维非结构矢量有限元数值模拟. 地球物理学报, 58(8): 2827-2838. DOI:10.6038/cjg20150817 |
叶益信, 李予国, 刘颖, 等. 2016. 基于局部加密非结构网格的海洋可控源电磁法三维有限元正演. 地球物理学报, 59(12): 4747-4758. DOI:10.6038/cjg20161233 |
殷长春, 贲放, 刘云鹤, 等. 2014. 三维任意各向异性介质中海洋可控源电磁法正演研究. 地球物理学报, 57(12): 4110-4122. DOI:10.6038/cjg20141222 |