耦合势有限体积法模拟海洋移动导线源三维频率域电磁响应

周建美 李 貅 戚志鹏(长安大学地质工程与测绘学院,陕西西安 710054)

1 引言

海洋可控源电磁法(MCSEM)是一种重要的海上油气勘探方法[1,2]。典型的MCSEM勘探采用水平导线源(长度一般为几百米)随船拖拽移动发射低频信号，在海底布设接收机阵列测量水平电场和三分量磁场响应[3]。在海洋电磁勘探设计及资料处理和解释过程中,均需要进行大量的数值模拟。三维正演算法研究一直都是电磁法的研究热点。MCSEM的三维正演相比其他电磁方法，有其特殊性：由于海水的滤波作用，MCSEM的工作频率非常低(0.1～10Hz)，导致控制方程呈病态，出现所谓的低感应数问题[4]； MCSEM的数值模拟需要同时计算近区和远区的电磁场响应(如发射源频率为0.25Hz时一般需要求解发射源附近±15km范围内的电磁场分布)，由于在大区域范围内电磁场的近场和远场的振幅相差很大(10个数量级以上),同时由于空气层的存在使得模型电导率差很大,导致离散化方程的条件数非常大[5]； MCSEM是一种典型的移动源电磁方法，需要对多次移动源进行三维模拟[6]。

目前主要采用有限差分法(FDM)[7-10]、有限体积法(FVM)[11]和有限元法(FEM)[12]模拟复杂三维地形的MCSEM响应。这三种方法的共同点是直接对Maxwell方程进行离散处理，得到一个大型的稀疏代数方程，借助于迭代法[13]或直接法[14]得到方程的数值解。直接法由于其对于矩阵的条件数不敏感,矩阵求解稳定,能有效处理移动源的三维正演问题[15]，广泛应用于MCSEM电磁响应的数值模拟[12,14]，但直接法的内存需求和计算时间随着未知数个数呈指数增长[16]，极大地制约了直接法用于复杂模型三维正演计算的效率。低感应数问题和大条件数的系数矩阵使常规的迭代算法很难有效求解MCSEM的离散线性方程[5]。利用电场的矢势和标势分解可有效克服低感应数问题[17]，结合有效的预处理能够实现MCSEM离散的大条件数线性方程的求解[5,13,18]。

在MCSEM三维正演计算中，一般将导线源近似为电偶极源。对于电偶极源采用直接离散源项的方法，但电偶极源的奇异性将导致源附近区域存在较大的离散误差。目前主要采用一次场/二次场方法利用等效源[13,14]或者差异场技术[15]克服源的奇异性问题。一次场/二次场方法能够有效模拟复杂场源形态的MCSEM三维响应[14]，但需要额外计算所有网格节点上的一次场分布；差异场技术虽然能够消除源附近的离散误差，但需要两次三维正演计算。对于实际MCSEM勘探中采用导线源的源项直接离散技术，这不同于电偶极源，在保证离散网格较导线源长度足够小时可得到源附近较为精确的电磁场分布。

MCSEM是一种典型的移动源电磁方法。对于移动源三维正演计算，常规迭代算法设置迭代初始解为零，总的计算时间随着移动源个数呈线性增长。设定移动发射源位置为网格中心位置，发射源移动，离散网格大小不变，但整体随着发射源移动。由于发射源移动前后地下电性分布变化不大，因此发射源移动前后的电磁场分布也非常接近。选择发射源移动前的电磁场分布作为发射源移动后电磁场的初始解，将有效减少迭代次数从而减少总的计算时间。

本文从Maxwell方程出发，首先通过引入电场的矢势标势分解，将Maxwell方程转化为一个关于矢势和标势的混合亥姆霍兹方程定解问题[17]；采用洛伦兹规范[16]，借助Yee氏交错网格和FVM对旋度和散度算子进行离散，得到稀疏对称系数矩阵的离散方程；采用直接加源技术实现导线源的离散；采用基于ILUT预处理的BICGSTAB迭代算法[17]实现离散方程的迭代求解。采用改进的迭代算法实现移动源的快速正演，即采用发射源移动前的电磁场响应作为发射源移动后的电磁场分布的初始解，从而提高迭代收敛速度，减少计算时间。

2 耦合势有限体积算法

2.1 控制方程

设电磁场随时间变化关系为eiω t，MCSEM三维正演对应的控制方程为

(1)

(2)

(3)

(4)

为了保证控制方程的完整性，对式(1)两边求散度，得到

(5)

将式(5)展开并利用洛伦兹规范进行简单的整理，得到

(6)

将式(4)与式(6)联立，表示为矩阵形式

(7)

(8)

由式(8)可知，该控制方程的系数矩阵是稀疏对称阵。采用x方向的线电流源，式(8)进一步简化为

D1Ax+Cxφ=-iωμ0Jx

(9a)

D1Ay+Cyφ=0

(9b)

D1Az+Czφ=0

(9c)

(9d)

边界处采用简单的截断边界条件[9]

(A,φ)|∂Ω=0

(10)

式中Ω为求解区域。

2.2 有限体积离散

式(9a)对应于x方向的稀疏矩阵方程的有限体积离散，引入x方向对应的半整数控制体积单元

(11)

(12)

定义控制体积单元内的场值是均一的，则有

(13)

(14)

式(9b)～式(9d)左边项的有限体积离散形式与式(9a)类似。

2.3 源项直接离散

源项离散包括式(9a)和式(9d)右边项的有限体积离散。

对于式(12)右边项的离散，本文采用直接加源方式。为了简单起见，设定网格与发射源位置重合，发射源的两个端点分别为(xa,yc,zd)和(xb,yc,zd)，其中xa和xb分别对应x方向相邻的两个网格节点坐标，yc和zd分别为y方向和z方向某一网格节点坐标，则沿x方向长度为xb-xa的线电流源可以定义为

Jx=Iδ(y-yc)δ(z-zd)×

[H(x-xb)-H(x-xa)]

(15)

(16)

对于式(9d)右边源项的有限体积离散，对应的控制体积单元为

(17)

对于x方向的发射源

[δ(x-xb)-δ(x-xa)]

(18)

则对应的有限体积离散为

(19)

2.4 移动源控制方程的快速求解

对离散结果进行整理,得到关于矢势A和标势φ的方程组

(20)

3 算法验证

3.1 精度验证

通过与一维模型的半解析计算结果进行对比，验证FVM的计算精度； 通过对比偶极源与线源的计算精度，验证源项直接离散方法的有效性。一维模型如图1所示，空气层电阻率设为1.0×106Ω·m，海水层电阻率为0.25Ω·m，海水层厚度为1000m，海底地层电阻率为1.0Ω·m。分别计算发射源为x方向的水平电偶极源(长度为1m)和线源(长度为400m)两种情况，发射源极矩均为1A·m，发射源位于海底上方100m，发射频率为0.25Hz。接收机位于海底，接收沿测线方向的电场Ex。剖分网格单元数为121×121×101，最小网格长度为100m，即线源为4个网格长度，偶极源为1个网格长度，总的网格剖分范围如图1所示。分别采用基于ILUT预处理的BICGSTAB算法和基于Jacobi预处理的QMR算法[14]求解该一维模型。图2为采用电偶极源时两种算法的收敛曲线对比，得到小于1.0×10-8的归一化残差，基于Jacobi预处理的QMR算法需要3316次迭代，计算时间为918s；而基于ILUT预处理的BICGSTAB算法只需要336次迭代，计算时间为416s， 可见后者计算速度更快，本文在后面的数值计算中都采用该算法。本文的FVM计算结果与开源软件Dipole1D[19]的半解析计算结果对比如图3所示,图4为线源与偶极源响应结果对比。

图1 海洋一维模型

图2 两种迭代求解算法的收敛曲线对比

对比图4a与图4b可知，线源与偶极源的响应只是在发射源附近区域(偏移距小于2500m)存在差异(图3中的红色曲线清晰地显示了在发射源附近区域两种发射源响应的差别)，在偏移距大于2500m的区域，两种场源的响应基本相同，与文献[20]的结论一致。对于MCSEM，异常一般都出现在大偏移距区域，因此MCSEM采用偶极源近似线源是合理的。由图4c可知，在偏移距大于2500m的区域， FVM计算的线源与电偶极源的振幅的相对误差基本一致，均小于3%；在发射源附近，FVM计算的电偶极源振幅的相对误差较大，达到了23%，线源振幅的最大相对误差小于7%。由图4d可知，不论是电偶极源还是线源，在整个求解区域(包括源所在位置)的相位差均小于2°； 线源发射的情况下，源区域的相位差更小。综合图4c和图4d可知：①对电偶极源采用直接加源技术，会导致源附近区域的振幅存在较大误差；②对线源采用直接加源技术，并保证源附近区域网格长度相比发射源长度足够小(本文中发射源长度为400m，源附近网格长度为100m)，能够得到发射源附近区域较为精确的响应。

图3 海洋一维模型正演响应(a)电偶极源的Ex振幅； (b)电偶极源的Ex相位； (c)线源的Ex振幅; (d)线源的Ex相位

图4 FVM计算的电偶极源与线源计算结果对比(a)线源与偶极源Ex振幅比值； (b)线源与偶极源Ex相位比值; (c)电偶极源和线源相

3.2 单异常体模型

本文采用迭代法计算移动源模型的三维正演，利用解的初值重置技术减少迭代次数从而减少计算时间。为了分析算法的有效性，设计如图5所示的三维模型。在海底下方1000m处存在一个直径4km、厚度100m、电阻率为100Ω·m的圆形高阻油气藏。发射源为沿x方向、长度为400m的线电流源，置于海底上方100m，发射频率为1.0Hz。发射源沿x方向移动，设定初始发射源与油气藏中心在x方向的距离为2500m，每隔500m移动一次，共移动11次。接收机位于海底，接收沿测线方向的电场Ex。采用FVM算法计算该移动源模型，剖分网格单元数为141×121×112，最小网格长度为100m(图6)。

图5 单异常体海洋三维模型

图6 移动源海洋3D模拟网格剖分示意图(a)发射源移动前； (b)发射源移动后

图7 FVM计算的海洋移动源三维模型正演响应(a)Ex振幅； (b)Ex相位

图8 移动源正演计算迭代收敛曲线

图9 移动源正演响应计算时间对比

3.3 双异常体模型

双异常体三维模型如图10所示。在海底下方1000m深度处存在两个直径为2km、厚度为100m、电阻率为100Ω·m的高阻油气藏。发射源为沿x方向、长度为400m的线电流源，位于海底上方100m，发射频率为1.0Hz, 发射源沿x方向移动。设初始发射源距离第一个目标体边界在x方向上的距离为500m，每隔1000m移动一次，共移动7次。接收机位于海底，接收沿测线方向的电场Ex。剖分网格单元数为201×121×112，最小网格长度为100m，网格以发射源为中心。

图10 双异常体模型

图11 双异常体模型正演响应(a)Ex振幅； (b)Ex相位

图12 双异常体模型移动源正演计算迭代收敛曲线

图13 双异常体模型移动源正演响应计算时间对比

3.4 非均匀覆盖层模型

为了进一步研究本文算法的适用性，设计如图14所示的非均匀覆盖层的三维模型。 在海底下方1000m深度处存在一个直径为4km、厚度为100m、电阻率为100Ω·m的圆形高阻油气藏。海水下方覆盖层中有两个电阻率分别为10Ω·m和5Ω·m的三维异常体。发射源为沿x方向、长度为400m的线电流源，发射源始终保持在海底上方100m处，发射频率为1.0Hz，发射源沿x方向移动。设定初始发射源距离高阻油气藏异常体边界在x方向上的距离为500m，每隔500m移动一次，共移动11次。接收机位于海底沿地形布设，接收沿测线方向的电场Ex。网格剖分单元为141×121×112，最小网格长度为100m，网格以发射源为中心。

图14 海洋非均匀覆盖层三维模型(a)xoz平面图； (b)xoy平面图

图15 非均匀覆盖层三维模型响应曲线(a)Ex振幅; (b)Ex相位

图16 非均匀层覆盖层模型的正演计算迭代收敛曲线对比

图中实线为移动源11次正演常规算法(即迭代初始解X0=0)的迭代收敛曲线，虚线为本文改进算法(即迭代初始解X0=X1)计算的后10次正演(不包括第1次正演)的迭代收敛曲线

图17 非均匀覆盖层三维模型正演响应计算时间对比

4 结束语

本文采用耦合势有限体积法建立了一套快速计算海洋移动导线源电磁响应的三维数值模拟技术：基于电场矢势和标势分解将Maxwell方程转换为关于矢势与标势的亥姆霍兹方程；采用洛伦兹规范整理得到对称形式的离散系数矩阵；通过基于ILUT预处理的BICGSTAB算法迭代求解离散线性方程；采用直接加源技术实现导线源的精确快速模拟；针对移动源正演响应计算，提出采用解的初值重置技术加快正演速度，即采用发射源移动前的电磁场响应作为发射源移动后电磁场分布的初始解，提高迭代收敛速度，减少迭代次数。

数值实例计算结果表明，采用直接加源技术，通过将导线源表示为多个电偶极源的叠加，能够显著减少源附近区域的计算误差。对于移动源的正演计算，采用解的初值重置技术，选择发射源移动前的电磁场分布作为移动后控制方程迭代计算的初始解，当发射源移动前后空间中电磁场分布变化较小时，能显著地提高正演速度；当发射源移动前后空间中电磁场分布变化较大时，正演速度提高有限。

参考文献

[1] Constable S.Ten years of marine CSEM for hydrocarbon exploration.Geophysics,2010，75(5):A67-A81.

[2] 何展翔,孙卫斌,孔繁恕等.海洋电磁法.石油地球物理勘探,2006,41(4):451-457.

He Zhanxiang,Sun Weibin,Kong Fanshu et al.Marine electromagnetic approach.OGP,2006,41(4):451-457.

[3] Constable S and Srnka L.An introduction to marine controlled-source electromagnetic methods for hydrocarbon exploration.Geophysics,2007,72(2):WA3-WA12.

[4] Newman G A and Alumbaugh D L.Three-dimensio-nal induction logging problems,Part 2:A finite-diffe-rence solution.Geophysics,2002,67(2):484-491.

[5] Puzyrev V.A parallel finite-element method for three-dimensional controlled-source electromagnetic forward modelling.Geophysical Journal International,2003,193(2):678-693.

[6] Plessix R E,Darnet M and Mulder W A.An approach for 3D multisource,multi frequency CSEM modeling.Geophysics,2007,72(5):SM177-SM184.

[7] 殷长春,贲放,刘云鹤等.三维任意各向异性介质中海洋可控源电磁法正演研究.地球物理学报,2014,57(12):4110-4122.

Yin Changchun,Ben Fang,Liu Yunhe et al.MCSEM 3D modeling for arbitrarily anisotropic media.Chinese Journal of Geophysics,2014,57(12):4110-4122.

[8] 付长民,底青云,王妙月.海洋可控源电磁法三维数值模拟.石油地球物理勘探,2009,44(3):358-363.

Fu Changmin,Di Qingyun,Wang Miaoyue.3D nume-ric simulation of marine controlled source electromagnetics (MCSEM).OGP,2009,44(3):358-363.

[9] 严波,李予国,韩波等.任意方位电偶源的MCSEM电磁场三维正演.石油地球物理勘探,2017,52(4):859-868.

Yan Bo,Li Yuguo,Han Bo et al.3D marine controlled-source electromagnetic forward modeling with arbitrarily orientated dipole source.OGP,2017,52(4):859-868.

[10] 朱成,李桐林,杨海斌等.带地形频率域可控源电磁法三维反演研究.石油地球物理勘探,2016,51(5):1031-1039.

Zhu Cheng,Li Tonglin,Yang Haibin et al.3D controlled source electromagnetic inversion with topography in the frequency domain.OGP,2016,51(5):1031-1039.

[11] 杨波,徐义贤,何展翔等.考虑海底地形的三维频率域可控源电磁响应有限体积法模拟.地球物理学报,2012,55(4):1390-1399.

Yang Bo,Xu Yixian,He Zhanxiang et al.3D frequency-domain modeling of marine controlled source electromagnetic responses with topography using finite volume method.Chinese Journal of Geophysics,2012,55(4):1390-1399.

[12] 杨军,刘颖,吴小平.海洋可控源电磁三维非结构矢量有限元数值模拟.地球物理学报,2015,58(8):2827-2838.

Yang Jun,Liu Ying,Wu Xiaoping.3D simulation of marine CSEM using vector finite element method on unstructured grids.Chinese Journal of Geophysics,2015,58(8):2827-2838.

[13] 蔡红柱,熊彬,Michael Zhdanov.电导率各向异性的海洋电磁三维有限单元法正演.地球物理学报,2015,58(8):2839-2850.

Cai Hongzhu,Xiong Bin,Zhdanov M.Three-dimensional marine controlled-source electromagnetic mode-lling in anisotropic medium using finite element method.Chinese Journal of Geophysics,2015,58(8):2839-2850.

[14] 韩波,胡祥云,黄一凡等.基于并行化直接解法的频率域可控源电磁三维正演.地球物理学报,2015,58(8):2812-2826.

Han Bo,Hu Xiangyun,Huang Yifan et al.3-D frequency-domain CSEM modeling using a parallel direct solver.Chinese Journal of Geophysics,2015,58(8):2812-2826.

[15] 周建美,张烨,汪宏年等.耦合势有限体积法高效模拟各向异性地层中海洋可控源的三维电磁响应.物理学报,2014,63(15):159101-1.

Zhou Jianmei,Zhang Ye,Wang Hongnian et al.Efficient simulation of three-dimensional marine controlled-source electromagnetic response in anisotropic formation by means of coupled potential finite volume method.Acta Physica Sinica,2014,63(15):159101-1.

[16] Streich R.3D finite-difference frequency-domain mode-ling of controlled-source electromagnetic data:Direct solution and optimization for high accuracy.Geophy-sics,2009,74(5):F95.

[17] Haber E,Ascher U,Aruliah D et al.Fast simulation of 3D electromagnetic problems using potentials.Journal of Computational Physics,2000,163(1):150-171.

[18] Weiss C J. Project APhiD:a Lorenz-gauged A-decomposition for parallelized computation of ultra-broadband electromagnetic induction in a fully heterogeneous Earth.Computers & Geosciences,2013,58(1):40-52.

[19] Key K.1D inversion of multi component,multi frequency marine CSEM data:Methodology and synthetic studies for resolving thin resistive layers.Geophysics,2009，74(2):F9-F20.

[20] Streich R and Becken M.Electromagnetic fields ge-nerated by finite-length wire sources:comparison with point dipole solutions.Geophysical Prospecting,2011，59(2):361-374.