海洋电磁是海洋油气勘探的一种重要方法.该方法通常采用船载拖曳式水平电偶源发射电磁信号,海底阵列接收机或船载拖曳式接收机采集电磁信号,通过对采集的电磁信号分析处理可以获得地下的电性结构信息,从而识别出油气构造中含水还是含油,减少干井率.
近年来,随着人们对海底电性结构精细探测要求的不断提高,现有一、二维正反演方法已不能满足实际工作的要求.另一方面,计算设备的快速发展使得海洋电磁三维正反演方法变得可行,因此三维海洋电磁正反演技术在过去十几年里得到了快速发展.正演方面,Avdeeva等(2007)使用有限差分算法对频率域和时间域海洋可控源方法进行研究,发现时间域方法对于浅海油气储层具有更好的探测能力,Swidinsky和Edwards(2010)利用积分方程三维正演算法研究了海底浅层高阻体对背斜高阻油气藏的时间域海洋电磁响应的影响特征,而Um等(2010)利用非结构有限元对时间域海洋电磁三维问题进行正演模拟,实现了晚期道电磁响应的精确求解.Zhang等(2018)进行了双船拖曳方式带地形的频率域海洋电磁响应三维正演模拟研究,发现海洋电磁数据解释过程中海底地形的影响不能忽略.反演方面,Gribenko和Zhdanov(2007)进行了基于积分方程的三维海洋可控源数据反演,并采用正则化聚焦方法获得了尖锐的油层边界,而Commer和Newman(2007)采用有限差分法实现了三维频率域海洋电磁数据带地形反演.Zach等(2008)使用基于近似海森矩阵的最优化方法反演海洋可控源电磁数据,并通过较少的迭代次数得到好的反演结果,Sasaki(2013)使用高斯牛顿最优化方法进行了海洋可控源数据、海洋大地电磁数据的联合反演,分辨率得到显著提高,而Zhdanov等(2014)使用共轭梯度法对拖曳式海洋电磁数据三维反演进行了研究,验证了拖曳式海洋电磁观测方式的可行性.贲放(2016)基于交错网格有限差分和L-BFGS方法进行了海洋可控源数据的三维反演,验证了L-BFGS进行海洋电磁三维反演的有效性.然而,到目前为止,三维海洋电磁反演均在频率域进行.由于频率域海洋电磁受到空气波影响严重(殷长春等,2012),因此频率域海洋电磁在浅水区难以获得很好的应用效果.时间域海洋电磁方法可根据空气波和地层波到达接收机的时间不同将它们分离,因此,时间域方法相比于频率域方法在解决浅海问题时具有很大的优势(Weiss,2007).Avdeeva等(2007)、Connell和Key(2013)对浅海环境下频率域和时间域海洋电磁方法进行了对比分析得出一致结论:时间域方法在浅海环境下应用效果更好.为此,本文中我们开展时间域海洋电磁三维反演研究.考虑到海水的高导特性,下阶跃波的建场时间过长,当前海洋发射设备难以实现长时间供电,为此我们选用实际勘探中使用的方波进行模拟.首先介绍基于非结构网格时间域海洋电磁三维矢量有限元正演理论和算法,计算下阶跃波电场响应并进行精度验证;进而,计算针对典型高阻异常体模型的电磁响应,并分析验证该电磁响应对高阻异常体的识别能力.考虑到有限内存L-BFGS方法收敛速度快、内存需求较小,本文选择L-BFGS方法进行时间域海洋电磁三维反演.我们通过理论模型反演验证反演方法的有效性.
1 正演理论 1.1 时间域有限元方法时间域海洋电磁的电场扩散方程可写为
(1) |
式中,e(r, t)为t时刻位置r上的电场,μ为介质的磁导率,σ为电导率,js为源电流密度.
对公式(1)采用基于棱边的非结构矢量有限元法进行求解.为此,在每一个小单元中,我们将电场表示为
(2) |
式中,nie(r)和eie(t)分别为小单元e内的矢量基函数和第i个棱边上的电场(金建铭, 1998).
由公式(1)和(2)并应用伽辽金方法,可得到关于小单元e的有限元方程:
(3) |
其中,Ae、Be、Se分别为单元e的质量矩阵、刚度矩阵和源项.矩阵Ae和Be的第(i, j)个元素、矩阵Se的第i个元素分别表示为
(4) |
(5) |
(6) |
为了对(3)式中的时间导数进行离散,我们采用针对任意时间步长无条件稳定的后推欧拉方法(齐彦福等,2017)近似公式(3)中电场对时间的导数,其形式为
(7) |
其中,下角标i,i+1和i+2分别表示第i、i+1和i+2个时间道,Δt表示时间步长.将(7)式代入(3)式并组装所有单元,可得
(8) |
为了求解(8)式,我们施加第一类边界条件,即
(9) |
其中,n为计算域外边界Γ的法向量.
对所有时间道的公式(8)进行组装,可得到如下矩阵方程:
(10) |
其中,Cn=3A+2ΔtnB,sn=2ΔtnSn.公式(10)可简写为
(11) |
其中,K为全时正演系数矩阵.
1.2 初始电场从公式(10)中可以看出,时间域电磁模拟需要初始时刻各棱边上的电场值.当采用电性源发射上阶跃波或方波时,初始电场为0;然而,当采用电性源发射下阶跃波时,初始电场为
(12) |
其中,eDC为直流电场,LDC为插值函数,φ为节点上的电位.节点上的电位可通过节点有限元求解泊松方程获得,即
(13) |
而边界条件采用Dirichlet边界条件:φ|Γ=0.为节省篇幅,具体数值算法本文不作叙述,可参考任政勇等(2017)和Um(2011).
1.3 精度验证为了验证本文的正演模拟精度,我们设计了一个海底均匀半空间模型,图 1为该模型的示意图.海水深度为200 m,电阻率0.3 Ωm,海底沉积层电阻率0.7 Ωm.发射源为一个长100 m的水平电偶源,放置于距海底50 m处.发射波形为下阶跃波,发射电流为1 A.接收点位于距海底30 m处,观测方式为同线观测, 收发距(T-R)1000~7000 m.图 2给出了发射波形为下阶跃波时本文计算结果与一维半解析解的对比.从图可以看出本文的计算结果与一维半解析解拟合较好,最大相对百分比误差小于5%,验证了本文算法的精度.
本文选择方波作为时间域海洋电磁发射波形.为了得到方波电磁响应,我们可通过两个阶跃响应反向移位叠加得到(齐彦福等,2015).
为了分析浅海环境中方波电磁响应对高阻异常体的识别能力,我们设计如图 3所示的高阻模型:海水深度200 m,海水电阻率0.3 Ωm,沉积层电阻率0.7 Ωm,异常体电阻率100 Ωm,异常体大小为4000 m×1500 m×300 m,其顶部埋深为800 m.以异常体的中心在海平面的投影为坐标原点,异常体长轴沿x轴方向,短轴沿y轴方向,z轴垂直向下为正.发射源长100 m,沿x轴方向,距海底50 m,发射波形为宽度1 s的方波,发射源中心点坐标为(-3000 m,0 m,150 m),我们设计的7个测点均距海底30 m,其位置坐标为x=-2000~4000 m,y=0 m.图 4给出了存在高阻异常体的电场响应和没有高阻异常体时背景响应的对比.从图可以看出,对于所有收发距电场Ex均在1 s处存在峰值.对于较大收发距,电场曲线出现第二个峰值,1 s处的峰值是由于在0 s开始供电电场逐渐增大,到1 s开始断电电场开始衰减,对于小收发距空气波和地层波到达时间均较早,两者被掩盖在了供电过程中的电场中,对于大收发距空气波到达较早被掩盖在了供电过程中的电场中,但地层波到达较晚,与供电中的电场和空气波在时间上分开,形成第二个峰值.所以时间域方法在浅海条件下不受空气波的影响,可用于浅海环境.对于存在高阻异常体时,地层波在高阻中传播较快,到达时间早于背景模型,与背景电场相比存在较大的异常,所以可用于反演.
针对时间域海洋电磁三维反演,我们采用如下正则化目标函数:
(14) |
其中,dobs和d分别为Nd×Nc维观测数据和模拟数据,Nd为接收点的个数,Nc为观测的时间道个数.m和m0分别为Nm维电导率模型参数向量和参考模型参数向量.λ为正则化因子,Cd为数据加权算子,其元素为反演数据噪声标准差的倒数.Cm为模型粗糙度算子.将公式(14)两端分别对m求偏导数得到目标函数的梯度为
(15) |
其中, J为灵敏度矩阵,rd=CdTCd(dobs-d),rm=CmTCm(m-m0).
2.2 灵敏度矩阵的计算为求得目标函数的梯度,需要先计算灵敏度矩阵.本文通过伴随正演隐式求解灵敏度矩阵(即求解灵敏度矩阵与向量的乘积).该方法避免直接存储灵敏度矩阵,可大大减少计算过程中的内存消耗.灵敏度矩阵J可表示为
(16) |
式中,插值算子L可表示为L=LtLsLd,其中Ld将计算的棱边上的电场值插值到接收点处,Ls将接收点处的电场合成为方波电场响应,而Lt将计算时间道的电场插值到各观测时间道.
将式(11)左右两端对模型参数求导可得
(17) |
则JTrd可表示为
(18) |
假设
(19) |
则有
(20) |
由于公式(10)中的A和C均是对称矩阵,KT可写为
(21) |
从公式(21)可知,求解u的伴随正演中我们需要一个逆向的递推过程.在完成伴随正演后,将u代入(18)式中,可得
(22) |
为了得到灵敏度矩阵和向量的乘积,我们还需计算
(23) |
必须指出,当使用下阶跃波进行计算时,初始电场中包含了地下电导率信息,因此公式(10)(11)中的另一项
(24) |
其中,解析项
我们采用Key(2016)在二维海洋电磁反演中提出的方法,在三维时间域海洋电磁反演中的模型粗糙度(model roughness)可写为
(25) |
(26) |
(27) |
(28) |
其中,Vi为四面体i的体积,Ni为与四面体i共顶点的四面体个数,Δrij为四面体i的中心到与其共顶点的四面体j的中心的距离,ωj为四面体j的体积与其共顶点的四面体的体积和的比值.
3 理论数据反演为了验证本文时间域反演算法的有效性,我们分别设计了水平海底、起伏海底和一个复杂海底构造模型.对于这三个模型,均采用x方向的电场Ex进行反演,并在反演数据中加入3%的高斯噪声.考虑到时间域海洋电磁三维反演计算时间和内存消耗均较大,本文选择有限内存L-BFGS方法进行数据反演(Nocedal and Wright, 1999).所有计算均使用处理器3.20 GHz、运行内存128 GB的Dell工作站.
3.1 水平海底模型水平海底模型与图 3相同.本次反演采用3条测线,y方向测线范围为-1000~1000 m,线距为1000 m,每条测线采用7个发射源,发射源中点的位置坐标为x=-3000~3000 m.每条测线上81个测点,测点的位置坐标为x=-4000~4000 m,y=0 m,间距为100 m.每一次观测发射船固定在一个发射位置,接收船记录该测线所有接收点位置处的电场.对于一次发射观测时间范围为0.9~45 s,时间道数为43道.图 5a、5b分别为正演和反演使用的网格.正演四面体网格数为818097,反演网格数为784529.初始正则化参数为λ=0.1.如果数据拟合差RMS有明显下降,则正则化参数保持不变;否则,正则化参数降为原来的十分之一.反演迭代次数共计39次,用时190小时.反演过程中的各参数变化如图 6所示.从图 6a、6b和6c可以看出,数据拟合差(RMS)和目标函数下降较快,RMS最终趋近于1,而模型粗糙度趋于稳定,说明本文时间域反演方法对于浅海模型数据的收敛性很好.图 7给出了x、y、z三个方向的反演结果与真实模型的对比.从图可以看出,本文反演得到的模型与真实模型吻合较好,地下高阻异常体的位置和电阻率得到很好的恢复.图 8a给出了发射源位于x=-3000 m,y=0 m处,测线y=0 m上单一测点不同时间道的初始模型响应、预测的电场响应数据与反演数据的对比,而图 8b和8c分别给出发射源分别位于测线y=0上7个不同位置时,同测线上不同测点t=2 s和t=10 s时的初始模型响应、预测的电场响应数据与反演数据的对比.从图 8可以看出,在存在噪声影响的情况下,预测数据仍很好地拟合了反演数据,验证了本文反演方法的有效性.
实际海底地形往往是起伏的,而起伏地形的存在会使电磁响应曲线形态发生变化(殷长春等,2019),为此我们设置了一个起伏海底模型来验证本文算法的可行性.本文设计的起伏海底地形模型参见图 9.海底最大起伏高度为100 m,其他计算参数与水平海底模型相同.图 10a、10b分别给出正演和反演网格,正演四面体网格数为784296,反演四面体网格数为748353.初始正则化参数为λ=0.1,其改变方式与水平海底模型计算时相同.反演迭代次数为41次,总共用时193小时.反演过程中的各参数变化如图 11所示,从图 11a—11c可以看出,RMS和目标函数下降较快且RMS最终趋近于1,模型粗糙度趋于稳定,说明本文的时间域方法对于浅海起伏地形电磁数据反演的收敛性也很好.图 12给出了x、y和z三个方向的反演结果与真实模型的对比.从图可以看出,本文反演算法得到的模型与真实模型十分接近,说明本文反演算法对浅海起伏地形环境的适应性.图 13a给出了发射源位于x=-3000 m,y=0 m处,测线y=0 m上单一测点不同时间道的初始模型数据、预测数据与反演数据的对比,图 13b和13c分别给出发射源位于测线y=0 m上7个不同位置,观测时间为t=2 s和t=10 s时y=0测线上不同测点的初始模型数据、预测数据与反演数据的对比.从图 13可以看出,虽然反演数据受噪声影响变得不光滑,但预测数据仍然可以很好地对其进行拟合,验证了本文算法对起伏海底地形电磁数据反演的有效性.
为了更进一步说明本文时间域反演算法的有效性,我们设计了一个更为复杂的海洋地电模型.如图 14所示,该模型含有起伏海底地形(参见图 9)和多个构造体.其中,A为4000 m×1500 m×300 m、电阻率为100 Ωm的拱形储油层,B为电阻率1000 Ωm的高阻岩体,C为高阻基底,电阻率为1000 Ωm,D、E和F为电阻率10 Ωm的局部异常体.电磁系统参数及计算参数与水平海底模型中的参数一致.我们使用三条测线数据进行反演,图 15给出反演结果.从图可以看出,本文的时间域反演算法很好地恢复了高阻储油层A和高阻基底C的位置,然而高阻储油层的形状没能很好地恢复.对于高阻岩体B,虽然其位置可以得到恢复,但刻画效果较差.这主要是由于高阻岩体B形态较为直立,电磁波在其中传播的横向传播距离较短,电磁导波特征不明显,因此观测异常中无法得到清晰体现.
本文基于非结构有限元法和后推欧拉方法,进而利用伴随算法并结合有限内存L-BFGS方法,成功实现了起伏海底条件下时间域海洋电磁三维反演.针对水平海底、起伏海底及复杂海底构造模型的合成数据反演,表明本文算法具有很高的稳定性,反演结果较好地刻画了海底水平高阻异常体的分布特征.特别是在浅水区域,由于频率域海洋电磁法受空气波影响,难以对海底高阻储油构造进行有效刻画的情况下,本文提供的算例表明,基于非结构有限元的模拟算法,结合有限内存的BFGS反演算法可以对浅水区海底油气构造的分布特征进行较为精细的刻画.由此我们可以得出结论:借助于本文的三维反演算法,时间域海洋电磁勘查技术可以实现浅水区海底油气勘查.由于目前时间域海洋电磁系统仍在研发之中,尚缺乏实测数据.利用本文算法进行实测海洋数据的反演将是我们未来的研究目标.
致谢 作者向本文审稿人及编辑表示由衷感谢.另外,也向吉林大学殷长春电磁研究团队的其他成员表示感谢.
Avdeeva A, Commer M, Newman G A. 2007. Hydrocarbon reservoir detectability study for marine CSEM methods: time-domain versus frequency-domain.//77th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts.
|
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.
|
Cai H Z, Xiong B, Zhdanov M S. 2015. Three-dimensional marine controlled-source electromagnetic modelling in anisotropic medium using finite element method. Chinese Journal of Geophysics (in Chinese), 58(8): 2839-2850. DOI:10.6038/cjg20150818 |
Commer M, Newman G A. 2007. 3D CSEM modeling and inversion for hydrocarbon reservoir mapping: The bathymetry problem.//77th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 452-456.
|
Connell D, Key K. 2013. A numerical comparison of time and frequency-domain marine electromagnetic methods for hydrocarbon exploration in shallow water. Geophysical Prospecting, 61(1): 187-199. DOI:10.1111/j.1365-2478.2012.01037.x |
Gribenko A, Zhdanov M. 2007. Special Section-Marine Controlled-Source Electromagnetic Methods-Rigorous 3D inversion of marine CSEM data based on the integral equation method. Geophysics, 72(2). |
Jin J M. 1998. Electromagnetic Finite-Element Method (in Chinese). Wang J G trans. Xi'an: Xidian University Press.
|
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 |
Liu Y H, Yin C C, Qiu C K, et al. 2019. 3-D inversion of transient EM data with topography using unstructured tetrahedral grids. Geophysical Journal International, 217(1): 301-318. DOI:10.1093/gji/ggz014 |
Nocedal J, Wright S J. 1999. Numerical Optimization. New York: Springer-Verlag.
|
Qi Y F, Yin C C, Wang R, et al. 2015. Multi-transient EM full-time forward modeling and inversion of m-sequences. Chinese Journal of Geophysics (in Chinese), 58(7): 2566-2577. DOI:10.6038/cjg20150731 |
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 |
Ren Z Y, Tang J T, Pan K J, et al. 2017. The FEM Simulation of DC Resistivity (in Chinese). Beijing: Science Press.
|
Sasaki Y. 2013. 3D inversion of marine CSEM and MT data:An approach to shallow-water problem. Geophysics, 78(1): E59-E65. |
Swidinsky A, Edwards R N. 2010. The transient electromagnetic response of a resistive sheet:an extension to three dimensions. Geophysical Journal International, 182(2): 663-674. DOI:10.1111/j.1365-246X.2010.04649.x |
Um S E, Harris J M, Alumbaugh D L. 2010. 3D time-domain simulation of electromagnetic diffusion phenomena:A finite-element electric-field approach. Geophysics, 75(4): F115-F126. |
Um S E. 2011. Three-dimensional finite-element time-domain modeling of the marine controlled-source electromagnetic method[Ph. D. thesis]. Stanford: Stanford University.
|
Weiss C J. 2007. The fallacy of the "shallow-water problem" in marine CSEM exploration. Geophysics, 72(6): A93--A97. |
Yin C C, Hui Z J, Zhang B, et al. 2019. 3D adaptive forward modeling for time-domain marine CSEM over topographic seafloor. Chinese Journal of Geophysics (in Chinese), 62(5): 1942-1953. DOI:10.6038/cjg2019L0798 |
Yin C C, Liu Y H, Weng A H, et al. 2012. Research on marine control-source electromagnetic method airwave. Journal of Jilin University (Earth Science Edition) (in Chinese), 42(5): 1506-1620. DOI:10.13278/j.cnki.jjuese.2012.05.001 |
Zach J J, Bjørke A, 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.
|
Zhang B, Yin C C, Liu Y H, et al. 2018. 3D forward modeling and response analysis for marine CSEMs towed by two ships. Applied Geophysics, 15(1): 11-25. |
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.海洋可控源电磁法三维正反演理论与各向异性影响机制研究[博士论文].长春: 吉林大学.
|
蔡红柱, 熊斌, Zhdanov M S. 2015. 电导率各向异性的海洋电磁三维有限单元法正演. 地球物理学报, 58(8): 2839-2850. DOI:10.6038/cjg20150818 |
金建铭. 1998.电磁场有限元方法.王建国译.西安: 西安电子科技大学出版社.
|
齐彦福, 殷长春, 王若, 等. 2015. 多通道瞬变电磁m序列全时正演模拟与反演. 地球物理学报, 58(7): 2566-2577. DOI:10.6038/cjg20150731 |
齐彦福, 殷长春, 刘云鹤, 等. 2017. 基于瞬时电流脉冲的三维时间域航空电磁全波形正演模拟. 地球物理学报, 60(1): 369-382. DOI:10.6038/cjg20170131 |
任政勇, 汤井田, 潘克家, 等. 2017. 直流电阻率有限单元法及进展. 北京: 科学出版社.
|
殷长春, 惠哲剑, 张博, 等. 2019. 起伏海底地形时间域海洋电磁三维自适应正演模拟. 地球物理学报, 62(5): 1942-1953. DOI:10.6038/cjg2019L0798 |
殷长春, 刘云鹤, 翁爱华, 等. 2012. 海洋可控源电磁法空气波研究现状及展望. 吉林大学学报(地球科学版), 42(5): 1506-1620. DOI:10.13278/j.cnki.jjuese.2012.05.001 |