积分方程法模拟层状介质中薄板的瞬态响应

刘愿愿，席振铢，蒋 欢，王 鹤

(1.中南大学 地球科学与信息物理学院，湖南 长沙 410083；2.中南大学 海洋矿产探测技术与装备研究所，湖南 长沙 410083



积分方程法模拟层状介质中薄板的瞬态响应

刘愿愿1,2，席振铢1,2，蒋 欢1,2，王 鹤1,2

(1.中南大学 地球科学与信息物理学院，湖南 长沙 410083；2.中南大学 海洋矿产探测技术与装备研究所，湖南 长沙 410083

瞬变电磁法高维正反演技术发展缓慢，主要是由于正演方法不够成熟，计算速度慢，精度低，达不到反演的要求。文阐明了薄板的定义，采用积分方程法模拟层状介质中薄板的瞬态响应，简化了网格剖分，控制矩阵大小，采取直接解法求取系数矩阵，与均匀半空间的解析解对比验证了该算法和程序的可靠性，通过误差分析体现了该算法的高精度。同时，研究了覆盖层与基岩对薄板响应结果的影响，发现高阻覆盖层和低阻覆盖层对薄板的响应影响都很小，高阻基岩会扩大薄板的响应，低阻基岩会屏蔽薄板的响应。该研究为瞬变电磁二维反演奠定了基础，也为以后生产中薄板解释提供了理论借鉴。

瞬变电磁法；积分方程法；层状介质；薄板

1 引 言

瞬变电磁法在一次脉冲磁场的间歇期测量二次场，简化了对异常响应的研究，因此具有更高的分辨率和灵敏度，近年来已成为物探热门方法之一[1,2]。但是和MT，CSAMT等相比，瞬变电磁法正反演大多都还停留在一维的水平，高维正反演尚不成熟[3-5]。

正演速度是制约反演速度的关键因素[6]。有限单元法需要进行离散化，网格剖分数目较多，计算速度慢，对内存要求也很高[7,8]。电磁波速度变化范围比较大，因此时间步长一般比较小，有限差分法需要迭代很多次才能完成正演过程，速度同样较慢[9-11]。

积分方程法是最早实现瞬变电磁高维模拟的数值算法，因为它只计算小体积异常区的场，计算速度快。1985年，SanFilipo和Hohmann首先利用时间域积分方程法求出均匀半空间中的对称异常体的三维响应，源是矩形的发射回线[12]。1986年至1988年，Newman和Hohmann等学者利用积分方程求得频率域时间响应，然后通过傅式变换转到时间域的响应[13]。国内的学者主要致力于研究电磁张量格林函数的数值计算。2007年，魏宝君等人提出了并矢格林函数的递推算法，引发了人们对早先使用的计算公式的修正与校验[14]。2009年，陈桂波系统地推导了积分方程的计算公式，并应用于各向异性地层中进行电磁场三维模拟[15]。2014年，胡俊华利用积分方程法进行三维正演[16]，赖刘保等人利用积分方程法模拟接地线源近区的三维薄板的响应[17]，说明积分方程法在瞬变电磁响应计算中日益成熟。

为了提高正演速度，本文在阐明薄板概念的基础上，采用积分方程法模拟层状介质下任意位置任意形态薄板的瞬态响应，并采用LU分解法求解矩阵。通过与均匀半空间的对比验证了该程序的正确性和精确度，研究了层状介质中覆盖层和基岩对于薄板响应值的影响。

2 理论基础

薄板，从理论上来说是厚度为0的导电板，这意味着感应电流将会集中在薄板所在的那个面上。实际上，薄板的定义必须根据实际情况来定。5 m厚的薄板在早期或者高频阶段不可以被当做薄板计算，但是在晚期或者低频阶段，40 m厚的薄板也可以当做薄板来模拟，因为导电覆盖层会滤掉高频，从而扩大薄板的定义范围。

由电磁理论可知，任意点r的电磁场E或者H可表示为源在围岩中的背景场Eb(Hb)和薄板引起的散射场Es(Hs)的叠加：

E(r)=Eb(r)+Es(r)

(1-1)

H(r)=Hb(r)+Hs(r)

(1-2)

根据并矢格林函数理论，散射磁场可表示为：

(2)

(3)

Js(r)=Δσ(r) ·E(r)

(4)

其中，Δσ是围岩与薄板电导率之差。把式(2)和式(4)代入式(1)，可得：

(5)

将异常体所在区域进行网格剖分，并假定每个单元内电场与电导率等于单元中心的值，则式(5)可写为：

(6)

薄板是二维体，因此系数矩阵一般不会很大，采取LU分解法求解式(6)，可求出空间中任意位置的磁场。

3 正确性验证

电阻率为100欧姆的均匀半空间中磁通密度B随时间t变化的解析表达式为：

(7)

同样，在电阻率为100Ω的均匀半空间指定一块100m×100m的区域看作薄板，分别是X方向(-50m，50m)，Z方向是(-50m，-150m)，在地面布置测线一条，(-200m,0m)到(200m,0m)，一共11个点，发射线圈边长100m×100m，电流为5A，使用该算法进行正演模拟，模拟结果同解析解进行对比，对比结果见图1，具体误差见表1。

由图1可见，积分方程法计算的响应值和解析解的响应值得出的曲线形态完全一致，验证了该算法的正确性。由表1可见，除开早期的几个节点，该算法的误差基本控制在3%以内，确保了该算法的高精度，可以计算更为复杂的模型。

图1 均匀半空间解析解与积分方程法对比图Fig.1 Comparison between the analytical solution and the integral equation method for the homogeneous half space

时间/ms误差百分比/%0.30631.110.40422.920.50217.110.60012.810.6989.700.8455.191.0412.641.2371.091.4330.121.6290.612.3152.133.0991.773.4911.764.8631.785.6471.376.4311.147.2150.738.3911.539.9590.8911.5270.6014.6630.2917.0150.2520.1510.68

4 实例模拟

4.1 两层介质下直立薄板模拟

为了清晰直观地反应层状介质对薄板的响应，首先定义直立薄板的电阻率为1 Ω，沿X方向厚度为0，位于零点处，Y方向是(-200 m，200 m)，Z方向是(-25 m，-75 m)，在地表Z=0的平面布置测线一条，沿(-400 m，0 m)至(400 m，0 m)，一共41个测点，发射线圈边长为100 m×100 m，电流为5 A，然后设定三种基岩模型：①是电阻率为100 Ω的均匀半空间；②是电阻率为100 Ω的基岩上覆20 m厚电阻率为1 000 Ω的覆盖层；③是电阻率为1 000 Ω的基岩上覆80 m厚的电阻率为100 Ω的覆盖层。模型示意图见图2，Z方向和X方向的响应分别见图3～图6。

由图3至图6可见，模型1和模型2的Z方向和X方向的响应在0.5 ms和2 ms曲线形态完全一致，数值相差也很小，可得高阻覆盖层对于薄板响应影响很小，几乎可以忽略。但是模型3靠近薄板时响应值变大，曲线光滑性明显变差，极值点处的响应更加尖锐，且越靠近晚期，区别更加明显，可见低阻基岩会增大上方介质中的薄板响应值，并且增幅会随着关断时间的增加而增加。

图2 直立薄板模型示意图Fig.2 Sketch map of vertical thin plate

4.2 两层状介质下水平薄板模拟

定义水平薄板的电阻率为1 Ω，沿Z方向厚度为0，位于Z=-50 m处，Y方向是(-200 m，200 m)，X方向是(-25 m，25 m)，在地表Z=0的平面布置测线一条，沿(-400 m，0 m)至(400 m，0 m)，一共41个测点，发射线圈边长为100 m×100 m，电流为5 A，然后设定三种基岩模型：①是电阻率为1 000 Ω的均匀半空间；②是电阻率为1 000 Ω的基岩上覆20 m厚电阻率为100 Ω的覆盖层；③是电阻率为100 Ω的基岩上覆80 m厚的电阻率为1 000 Ω的覆盖层。模型示意图见图7，Z方向和X方向的响应分别见图8～图11。

图3 t=0.5 ms直立薄板Z分量响应图Fig.3 Z component response diagram of vertical thin plate at t=0.5 ms

图4 t=2 ms直立薄板Z分量响应图Fig.4 Z component response diagram of vertical thin plate at t=2 ms

图5 t=0.5 ms直立薄板X分量响应图Fig.5 X component response diagram of vertical thin plate at t=0.5 ms

图6 t=2 ms直立薄板X分量响应图Fig.6 X component response diagram of vertical thin plate at t=2 ms

图7 水平薄板模型示意图Fig.7 Sketch map of horizontal thin plate model

图8 t=0.5 ms水平薄板Z分量响应图Fig.8 Z component response diagram of horizontal thin plate at t=0.5 ms

图9 t=2 ms水平薄板Z分量响应图 Fig.9 Z component response diagram of horizontal thin plate at t=2 ms

图10 t=0.5 ms水平薄板X分量响应图Fig.10 X component response diagram of horizontal thin plate at t=0.5 ms

图11 t=2 ms水平薄板X分量响应图Fig.11 X component response diagram of horizontal thin plate at t=2 ms

由图8至图11可见，模型4和模型5的Z方向和X方向在0.5 ms和2 ms时曲线基本重合，响应值相差很小，说明低阻覆盖层对于薄板的影响可以忽略。但是模型6基岩是低阻，所以Z方向背景场较大，但是Z方向和X方向薄板的响应值很小，和模型4和模型5的响应值相差约1个数据级，并且随着关断时间的增加，背景场逐步趋于重合，但是薄板响应值的差距不会减少，依然十分明显。可见低阻基岩会部分屏蔽掉上方介质中薄板的响应，并且这种屏蔽效果不会随着关断时间的增加而减少。

5 结 论

1)积分方程法适合模拟层状介质中薄板的响应，计算速度快，计算结果准确，为反演奠定了扎实的基础。

2)无论是高阻覆盖层或者低阻覆盖层对于下方基岩中薄板的响应值影响很小，几乎可以忽略。

3)基岩对于上方覆盖层中薄板的响应值有着显著的影响。高阻基岩会增大薄板的响应值，且增幅随着关断时间的增加而增加；低阻基岩会抑制薄板的响应值，这种影响不会随着关断时间的变化而变化。

致 谢 感谢湖南五维地质科技有限公司提供的帮助。

[1]牛之链.时间域电磁法原理[M].长沙:中南大学出版社,2007.

[2]李貅.瞬变电磁测深的理论与应用[M].西安:陕西科学技术出版社,2002.

[3]薛国强,李貅,底青云.瞬变电磁法正反演问题研究进展[J].地球物理学进展,2008,23(4):1 165-1 172.

[4]高阳,陈志军,熊华山.矩形框瞬变电磁法一维正演研究[J].工程地球物理学报,2015,12(6):726-731.

[5]李建慧,朱自强,曾思红,等.瞬变电磁法正演计算进展[J].地球物理学进展,2012,27(4):1 393-1 400.

[6]米萨克 N.纳比吉安.勘探地球物理电磁法[M].赵经祥,王艳君译.北京:地质出版社,1992.

[7]M.I.Epov, E.P.Shurina, O.V.Nechaev. 3D forward modeling of vector field for induction logging problems[J]. Russian Geology and Geophysics,2007,48(9):770-774.

[8]殷长春,张博,刘云鹤,等.2.5维起伏地表条件下时间域航空电磁正演模拟[J].地球物理学报,2015,58(4):1 411-1 424.

[9]Michael Commer, Gregory Newman. A parallel finite-difference approach for 3D transient electromagnetic modeling with galvanicsources[J]. Geophysics,2004,69(5):1 192-1 202.

[10]孙怀凤,李貅,李术才,等.考虑关断时间的回线源激发TEM三维时域有限差分正演[J].地球物理学报,2013,56(3):1 049-1 064.

[11]徐正玉,杨海燕,邓居智,等.回线源三维地—井瞬变电磁法FDTD数值模拟[J].工程地球物理学报,2015,12(3):327-332.

[12]Wiliam A Sanfilipo, Gerald W Hohmann. Integral equation solution for the transient electromagnetic response of a three dimensional body in a conductive half-space[J]. Geophysics, 1985,50(5):798-809.

[13] Gregory A Newman, Gerald W Hohmann. Transient electromagnetic responses of high-contrast prisms in a layered earth[J]. Geophysics,1988,53(5):691-706.

[14]魏宝君，王甜甜，王颖.用磁流源并矢Green函数的递推矩阵方法计算层状各向异性地层中多分量感应测井响应.地球物理学报[J],2009,52(11):2 920-2 928.

[15]陈桂波.各向异性地层中电磁场三维数值模拟的积分方程算法及其应用[D].长春:吉林大学,2009.

[16]胡俊华.瞬变电磁积分方程法正演模拟研究[D].武汉:中国地质大学,2014.

[17]赖刘保.接地线源近区三维薄板瞬变电磁响应积分方程法数值模拟[D].长沙:中南大学,2014.

The Simulation of Transient Response of Thin Plate in Layered Media by Integral Equation Method Modeling

Liu Yuanyuan1,2,Xi Zhenzhu1,2,Jiang Huan1,2,Wang He1,2

(1.SchoolofGeosciencesandInfo-Physics,CentralSouthUniversity,ChangshaHunan410083,China;2.MarineMarineExplorationandEquationresearchInstitute,CentralSouthUniversity,ChangshaHunan410083,China)

Transient electromagnetic method multidimensional forward and inversion technology development is slow, because its forward modeling method is not mature, the calculated speed is slow and the accuracy is low, which can not be up to the requirements of inversion. This paper clarifies the definition of sheet, uses integral equation method to simulate the transient response of sheets in layered medium, simplifies meshing and controls the size of matrix to take direct method to solve coefficient matrix. The reliability of this algorithm and this program is verified by the comparison of analytical solution of the homogeneous half space, and the accuracy of the algorithm is demonstrated by the error analysis. At the same time, the effect of the overburden and bedrock on the response of the thin plate is studied. It is found that the response of the high resistance and low resistance to thin plate is very small, and that high resistance bedrock would expand the response of thin plate, and that low resistivity bedrock would shield the response of thin plate. This study lays the foundation for 2-dimensional inversion of transient electromagnetic, and also provides a theoretical reference for the production of thin plate in the future.

transient electromagnetic method; integral equation; layered media; thin plate

1672—7940(2016)02—0143—06

10.3969/j.issn.1672-7940.2016.02.001

国家自然科学基金项目(编号:41304090)；深圳市战略新兴产业发展专项资金项目(编号：CXZZ20120618165608947)

刘愿愿(1991-)，男，硕士研究生，主要从事于电磁法正反演研究。E-mail:csu_lyy@163.com

席振铢(1966-)，男，教授，博士生导师，主要从事瞬变电磁方法与技术研究。E-mail:xizhenzhu@163.com

P631.3

A

2015-10-09