航空瞬变电磁法三维各向异性数值模拟和响应特征研究_管清辰.pdf

2023 年第 38 卷 第1期2023，38(1):0449-0465地球物理学进展Progress in Geophysicshttp:/wwwprogeophyscnISSN 1004-2903CN 11-2982/P管清辰，谭捍东，林昌洪2023 航空瞬变电磁法三维各向异性数值模拟和响应特征研究 地球物理学进展，38(1):0449-0465，doi:10 6038/pg2023GG0177GUAN QingChen，TAN HanDong，LIN ChangHong 2023 Three-dimensional modeling and response characteristics of airborne transient electromagneticdata for anisotropic medium Progress in Geophysics(in Chinese)，38(1):0449-0465，doi:106038/pg2023GG0177航空瞬变电磁法三维各向异性数值模拟和响应特征研究Three-dimensional modeling and response characteristics of airbornetransient electromagnetic data for anisotropic medium管清辰，谭捍东*，林昌洪GUAN QingChen，TAN HanDong*，LIN ChangHong收稿日期2022-06-16;修回日期2022-11-05投稿网址http:/www progeophys cn基金项目中国自然资源航空物探遥感中心青年创新基金课题(2020YFL05)、国家自然科学基金项目(41830429)和山西省重点研发计划项目(202102080301001)联合资助第一作者简介管清辰，男，博士研究生，主要从事地球物理电法的研究工作 E-mail:285230064 qq com*通讯作者谭捍东，男，教授，主要从事电法勘探理论和应用研究 E-mail:thd cugb edu cn中国地质大学(北京)地球物理与信息技术学院，北京100083School of Geophysics and Information Technology，China University of Geosciences，Beijing 100083，China摘要随着瞬变电磁法的快速发展，三维任意各向异性介质的数值模拟成为研究热点 本文从时间域的麦克斯韦方程组出发，采用时域交错采样有限差分法，推导了时域电磁场的时间分步迭代公式，实现了任意各向异性介质的航空瞬变电磁三维正演 设计地电模型，与已有的软件计算结果进行对比，检验了三维正演算法的计算精度 设计了典型的三维各向异性地电模型，改变各向异性参数计算了航空瞬变电磁响应，分析了各向异性参数对航空瞬变电磁响应的影响 开发的正演算法为研究瞬变电磁法各向异性响应特征和三维反演提供了重要技术支撑关键词航空瞬变电磁法;各向异性;三维正演;时域有限差分中图分类号P631文献标识码Adoi:10 6038/pg2023GG0177AbstractWith the rapid development of transientelectromagnetic method，numerical simulation of three-dimensional anisotropic media has become a hot researchtopicBased on Maxwells equation，we derived thediscrete iterative formula of electromagnetic field by usingtime domain finite difference method and realized thethree-dimensional forward modeling of airborne transientelectromagnetic data for anisotropic medium The accuracyof the forward algorithm is verified by comparing the resultswith those of other software We calculated the response ofa typical three-dimensional anisotropic geoelectric modeland analyzed the effect of anisotropic parameters on theairborne transient electromagnetic response The forwardalgorithm provides an important technical support for theresearch of anisotropic response characteristics and three-dimensional inversion of transient electromagnetic methodKeywordsAirborne transient electromagnetic sounding;Three-dimensional modeling;Anisotropy;Finite-differencetime method0引言由于环境恶劣、工作条件特殊，地面勘探方法在沙漠、森林覆盖区、湖泊、沼泽等地区的勘探效率不高 而航空电磁法采用飞行平台，无需地面人员布置仪器，可以实现快速勘探，在资源勘探领域拥有广阔的应用前景 航空瞬变电磁法是在地面瞬变电磁法的基础上发展起来的一种探勘方法，是航空电磁法地球物理学进展www progeophys cn2023，38(1)的一个重要分支 作为一种高效的时间域航空电磁法，该方法已经在矿产资源勘查，水资源勘探和环境检测等领域发挥重要作用，具有速度快、成本低、观测面积大等优点关于航空瞬变电磁正演方面的研究有很多，但大多是基于各向同性介质的成果(aiche，1998;罗延钟等，2003;殷长春等，2013)大量的基础理论研究和野外观测结果表明，电性各向异性现象在地球内部是普遍存在的 由于电各向异性对电磁数据的影响严重，关于各向异性的电磁法研究成为当前电磁领域的热点 已有的研究结果(Avdeev et al，1998;Yin and Fraser，2004;Liu and Yin，2014;殷长春等，2015)表明:航空电磁信号对地下介质各向异性相当敏感 利用传统的各向同性方法去解释本质上是各向异性的航空电磁资料，会导致不容忽视的误差，甚至错误的解释结果(Yin and Maurer，2001;Weiss and Newman，2002;Yin et al，2016)航空瞬变电磁资料的处理技术上，目前主要使用已经成熟的一维正反演技术，二维、三维的正反演仍在探索中，而考虑各向异性介质的航空瞬变电磁正反演研究很少，近几年才出现相关研究(Yin etal，2016;Huang et al，2017;齐彦福等，2020)电法勘探中考虑介质各向异性有助于建立更符合实际地质情况的模型，为电法数据的解释工作提供新的思路 正演是反演的基础和关键部分，本文在传统的时域有限差分算法的基础上，引入了各向异性电导率张量，推导了考虑电导率各向异性的时步迭代公式，并采用延拓的思想计算空中场值，最终实现了任意各向异性介质的航空瞬变电磁法三维正演算法，并设计模型进行了计算和分析1任意各向异性介质航空瞬变电磁法三维有限差分正演1 1介质的各向异性各向同性电导率为标量，而各向异性电导率为张量 电导率张量由三个主轴电导率 x，y，z以及走向角 S，倾向角 D，偏角 L三个欧拉旋转角共同决定 本算法采用的旋转变换如图 1 所示(Pek andSantos，2002)，先绕 z 轴旋转角度 S，然后绕 x 轴旋转角度 D，最后再绕 z 轴旋转角度 L 需要注意的是，当 D=0时，另外两个旋转角 S和 L的旋转效果是相同的，所以后边的讨论部分没有涉及旋转角 L设*为主轴电导率张量，其表达式如下:图 1坐标转换示意图(Pek and Santos，2002)Fig 1Transformation of the coordinate system(Pek and Santos，2002)*=x000y000z，(1)设 x和 z为关于 x 轴和 z 轴的旋转矩阵，其表达式如下:x()=1000cossin0 sincos ，z()=cossin0 sincos0 001，(2)考虑各向异性情况下的电导率张量具体可以通过以下计算求得，它包含九个元素:=xxxyxzyxyyyzzxzyzz=z(S)x(D)z(L)*z(L)x(D)z(S)(3)1 2任意各向异性介质航空瞬变电磁法电磁场满足的方程本文采用时域有限差分法(FDTD)直接模拟航空瞬变电磁响应 忽略位移电流，时间域麦克斯韦方程组如下(Wang and Hohmann，1993):b(t)t=e(t)，(4)js(t)+e(t)=h(t)，(5)b(t)=0，(6)j(t)=0，(7)并且:b(t)=h(t)，(8)其中 b(t)为磁感应强度，js(t)为源电流密度，h(t)为磁场强度，e(t)为电场强度，为磁导率，为电导率张量为了构成 FDTD 计算所需要的时间步进格式，0542023，38(1)管清辰，等:航空瞬变电磁法三维各向异性数值模拟和响应特征研究(www progeophys cn)图 2电磁场交错采样网格示意图Fig 2The staggered grid and distribution of electromagnetic field使用 Dufort-Frankel 方法对麦克斯韦方程组进行修正，在关断时间后源电流密度项 js(t)为零，然后式(5)左端加入一虚拟位移电流项后变成如下形式:e(t)t+e(t)=h(t)，(9)考虑各向异性情况下，电导率为张量，式(9)展开后形式如下:ext+xxex+xyey+xzez=hzyhyz，(10)eyt+yxex+yyey+yzez=hxzhzx，(11)ezt+zxex+zyey+zzez=hyxhxy(12)将式(4)、(6)展开后形式如下:bxt=ezyeyz，(13)byt=exzezx，(14)bzz=bxxbyy(15)1 3时间分步有限差分法采用时域交错采样有限差分法对求解区域进行网格剖分，将连续问题离散化，用有限节点代替连续空间，电场和磁场使用交错网格离散化 如图 2 所示，电场采样点都选取在棱柱边缘的中点，而磁场采样点均选取在棱柱面的中心 每一个电场分量的周围有四个磁场分量环绕;同样，每一个磁场分量的周围有四个电场分量环绕 如图 3 所示，电场和磁场在时间顺序上交替抽样，抽样时间间隔彼此相差半个时间步长，使麦克斯韦旋度方程离散以后构成显示差分方程，从而可以在时间上迭代求解 为了同时兼顾计算效率和计算精度，网格剖分采用的是非均匀的网格剖分方案 测区中异常体所在的核心区域采用较为密集的均匀网格剖分，随着距离的增加，网格边长按一定的步长扩大系数逐渐增大图 3电磁场时间采样示意图Fig 3Time sampling of electromagnetic field首先采用中心差分法将电场以及电场和磁场对时间的偏导数进行近似:en+1/2=en+1+en2，(16)e()tn+1/2=en+1 entn，(17)b()tn=bn+1/2 bn 1/2(tn 1+tn)/2，(18)然后将式(16)、(17)代入到式(10)(12)求解，即可得到如下的电场强度的迭代公式:en+1x=83enx+82(thn+1/2zyhn+1/2y)z 4t(2hn+1/2xzhn+1/2z)xxy(+hn+1/2yxhn+1/2x)yxz154地球物理学进展www progeophys cn2023，38(1)(hn+1/2zyhn+1/2y)zyy(hn+1/2zyhn+1/2y)zzz42tenx(xx yy zz)82t(xyeny+xzenz)+2t(3hn+1/2yxhn+1/2x)yxyyz(hn+1/2yxhn+1/2x)yxzyy(hn+1/2xzhn+1/2z)xxyzz(+hn+1/2xzhn+1/2z)xxzzy(+hn+1/2zyhn+1/2y)zyyzz(hn+1/2zyhn+1/2y)zyzzy2t2enx(xxyy xyyx+xxzz xzzx yyzz+yzzy)4t2eny(xyzz xzzy)4t2enz(xzyy xyyz)t3enx(xxyyzz+xxyxzy+xyyxzz+xyyzzx+xzyxzy+xzyyzx)/83+42t(xx+yy+zz)+2t2(xxyyxyyx+xxzz xzzx+yyzz yzzy)+t3(xxyyzz xxyxzy xyyxzz+xyyzzx+xzyxzy xzyyzx)，(19)en+1y=83eny+8(2hn+1/2xzhn+1/2z)xt+4t(2hn+1/2xzhn+1/2z)xxx(hn+1/2yxhn+1/2x)yyz(+hn+1/2xzhn+1/2z)xzz(hn+1/2zyhn+1/2y)zyz+42teny(xx yy+zz)82t(yxenx+yzenz)2t(3hn+1/2yxhn+1/2x)yxxyz(hn+1/2yxhn+1/2x)yxzyx(hn+1/2xzhn+1/2z)xxxzz(+hn+1/2xzhn+1/2z)xxzzx(+hn+1/2zyhn+1/2y)zyxzz(hn+1/2zyhn+1/2y)zyzzx2t2eny(xxyy xyyx xxzz+xzzx+yyzz yzzy)4t2enz(xxyz xzyx)4t2enx(yxzz yzzx)t3eny(xxyyzz xxyzzyxyyxzz+xyyzzx+xzyxzy xzyyzx)/83+42t(xx+yy+zz)+2t2(xxyyxyyx+xxzz xzzx+yyzz yzzy)+t3(xxyyzz xxyxzy xyyxzz+xyyzzx+xzyxzy xzyyzx)，(20)en+1z=83enz+8(2hn+1/2yxhn+1/2x)yt+4t(2hn+1/2yxhn+1/2x)yxx(+hn+1/2yxhn+1/2x)yyy(hn+1/2xzhn+1/2z)xzy(hn+1/2zyhn+1/2y)zzx+42tenz(xx+yy zz)82t(zxenx+zyeny)+2t(3hn+1/2yxhn+1/2x)yxxyy(hn+1/2yxhn+1/2x)yxyyx(hn+1/2xzhn+1/2z)xxxzy(+hn+1/2xzhn+1/2z)xxyzx(+hn+1/2zyhn+1/2y)zyxzy(hn+1/2zyhn+1/2y)zyyzx+2t2enz(xxyy xyyx xxzz+xzzx yyzz+yzzy)4t2eny(xxzy xyzx)+4t2enx(yxzy yyzx)t3enz(xxyyzz xxyzzyxyyxzz+xyyzzx+xzyxzy xzyyzx)/83+42t(xx+yy+zz)+2t2(xxyyxyyx+xxzz xzzx+yyzz yzzy)+t3(xxyyzz xxyxzy xyyxzz+xyyzzx+xzyxzy xzyyzx)(21)将式(18)代入方程对式(13)(15)求解，即可得到如下的磁感应强度迭代公式:bn+1/2x=bn 1/2x+tn 1+tn2enzyeny()z，(22)bn+1/2y=bn 1/2y+tn 1+tn2enxzenz()x，(23)bn+1/2z(zk)=bn+1/2z(zk+1)+zbn+1/2xx+bn+1/2y()y(24)对于场值对空间的偏导数项，我们同样采用中心差分法进行近似处理 之后，给定相应电磁问题的初始值，就可以逐步推进地求得以后各个时刻电磁场的空间分布1 4初始场和边界条件根据毕奥-萨伐尔定律，电流元在空间任意点产生的磁感应强度是可以计算的:db=04Idl rr3=04Idlsinr2，(25)对式(25)进行积分可以得到有限长导线产生的磁感应强度:2542023，38(1)管清辰，等:航空瞬变电磁法三维各向异性数值模拟和响应特征研究(www progeophys cn)b=04baIdlsinr2，(26)将多边形回线的边看作有限长导线，分别计算其产生的磁场，最后相加即可得到多边形回线产生的初始场由于空气的电导率太小，如果将空气层一起纳入进行剖分计算，将会导致迭代不稳定，所以本文只对地表以下部分进行剖分 采用将地表的场值向上延拓的方法计算空中场值，近似模拟地空边界条件地下其他边界，均采用第一类边界条件，将边界处的切向电场设置为零1 5空中场值的计算首先从自由空间的磁感应强度满足拉普拉斯方程出发:2b=0，(27)对其进行傅里叶变换将空间域的值转化到波数域，可得到如下公式(Wang and Hohmann，1993):Bx(u，v，z=0)=iuu2+v2Bz(u，v，z=0)，(28)By(u，v，z=0)=ivu2+v2Bz(u，v，z=0)(29)利用地表的 Bz通过计算可以得到地表的 Bx和By，然后空中的场值可以由地表的场值向上延拓得到:Bx(u，v，z=h)=exp(hu2+v2)Bx(u，v，z=0)，(30)By(u，v，z=h)=exp(hu2+v2)By(u，v，z=0)，(31)Bz(u，v，z=h)=exp(hu2+v2)Bz(u，v，z=0)(32)最后再通过反傅里叶变换，将波数域的值转换到空间域 采用公式(30)、(31)计算地空边界条件采用公式(32)计算空中的磁感应强度 bz，然后用相邻时间点的磁感应强度的差值除以时间间隔，即可得到磁感应强度对时间的导数 dbz/dt2三维正演算法计算精度检验为了验证航空瞬变电磁法任意各向异性介质三维正演程序的正确性，需要将本程序与其他软件的计算结果进行对比采用 80 80 55 的非均匀网格剖分，最小网格边长为 10 m，在 x 方向和 y 向上，以原点为中心，均匀网格的数量为 40 个，然后按 1 2 的扩大系数逐渐增大网格边长 在 z 方向上，均匀网格的数量为 30个，同样按 1 2 的扩大系数逐渐增大网格边长 剖分区域水平方向最远距离到 2500 m，垂直方向最深到地下 6000 m发射线圈的大小为 20 m 20 m 时间范围从从 105s 到 102s，供电电流为 1 A，发射波形为阶跃波2 1均匀半空间模型空中场值对比地下为电导率 0 01 S m1均匀半空间，发射线圈位于空中 30 m 处，接收点位于发射线圈中心 与丹麦 Aarhus 大学开发的 AarhusInv 软件(Auken etal，2015)计算结果进行对比，如图 4 所示，结果吻合的很好，误差都在 5%以内图 4与 AarhusInv 软件计算结果对比图Fig 4Comparison with results of AarhusInv software2 2任意各向异性半空间模型结果对比地下为任意各向异性半空间，三个主轴电导率分别为 x=0 03 S m1，y=0 006 S m1，z=001 S m1，三个旋转角 S=0，D=30，L=0，采用中心回线装置，发射和接收均位于地表 计算结果与殷长春(Yin and Weidelt，1999)的瞬变电磁各向异性一维正演软件计算结果进行对比，如图 5 所示，结果吻合的比较好354地球物理学进展www progeophys cn2023，38(1)图 5与一维各向异性软件计算结果对比图Fig 5Comparison with the results of one-dimensionalanisotropy software2 3二维主轴各向异性模型结果对比地下存在一沿 y 方向无限延伸的棱柱体，如图6 所示，围岩背景电阻率设为 0 01 S m1，异常体顶部埋深100 m，x 方向宽度为200 m，z 方向厚度100 m，y 方向无限延伸 棱柱体三个主轴电导率分别为 x=0 01 S m1，y=0 02 Sm1，z=0 02 Sm1，采用中心回线装置，发射和接收均位于地表 与本人之前用有限单元法开发的二维主轴各向异性瞬变电磁正演软件(管清辰，2019)计算结果进行对比 如图 7 所示，结果吻合的比较好图 6二维模型示意图Fig 6Sketch of two-dimensional model图 7与二维主轴各向异性软件计算结果对比图Fig 7Comparison with the results of two-dimensionalanisotropy software图 8三维模型侧视图Fig 8Side view of 3d model图 9三维模型俯视图Fig 9Top view of 3d model4542023，38(1)管清辰，等:航空瞬变电磁法三维各向异性数值模拟和响应特征研究(www progeophys cn)图 10各向异性介质模型示意图(x=zy)Fig 10Sketch of anisotropic medium while x=zy图 11y取不同值的航空瞬变电磁响应Fig 11Airborne transient electromagnetic response with different y通过与其他软件对各向异性半空间模型和二维主轴各向异性模型计算结果进行对比，检验了本文开发的三维任意各向异性正演算法的计算精度3三维各向异性地电模型响应计算和特征分析设计了典型的三维各向异性地电模型，用所实现的航空瞬变电磁三维各向异性正演程序计算响应，分析不同的各向异性参数对航空瞬变电磁响应的影响554地球物理学进展www progeophys cn2023，38(1)图 13z取不同值的航空瞬变电磁响应Fig 13Airborne transient electromagnetic response with different z图 12各向异性介质模型示意图(x=yz)Fig 12 Sketch of anisotropic medium while x=yz地下异常体为埋深 100 m，长、宽、高均为 200 m的立方体，测点点距和测线线距均为 20 m 发射线圈边长 20 m，接收点在线圈中心，高度距离地面30 m，测区范围 400 m 400 m，21 条测线，每条测线21 个测点，共 441 个测点 其他模型剖分参数和装置参数与前面正演验证部分相同 图 8 和图 9 为三维模型示意图，其中黑框为异常体，黑点为测点 异常体设置为规则立方体是为了消除因异常体三个方向构造尺度不同而产生的影响，使我们专注于研究各向异性参数对响应的影响6542023，38(1)管清辰，等:航空瞬变电磁法三维各向异性数值模拟和响应特征研究(www progeophys cn)图 14各向异性介质模型示意图(x=zy，S0)Fig 14Sketch of anisotropic medium while x=zyand S0754地球物理学进展www progeophys cn2023，38(1)图 15S取不同值的航空瞬变电磁响应Fig 15Airborne transient electromagnetic response with different S图 16各向异性介质模型示意图(x=zy，D0)Fig 16Sketch of anisotropic medium(x=zy，D0)3 1异常体为各向异性，背景围岩为各向同性3 1 1水平方向电导率的影响围岩电导率设置为 0 01 Sm1，异常体垂直方向电导率 z=0 02 Sm1，水平方向电导率 x=0 02 Sm1，y取 0 1 Sm1，0 02 Sm1，0 005 S m1 图 10 为该各向异性异常体的物理模型示意图，垂直方向电导率与其中一水平方向的电导率相同，而与另一水平方向的电导率不同 自然界中存在与该各向异性模型对应的地质情况，比如水平层状的沉积岩由于后期构造运动产生类似的直立岩层，或者各种岩体在应力作用下发生断裂，形成垂直方向的节理 由于异常体模型是规则立方体，异常体水平方向电导率 x和 y对响应的影响是对称的，我们只需展示 y变化的响应，x变化的响应只是方向不同响应如图 11 所示，我们选取了两个时间道 t=0 1 ms 和 t=0 4 ms，绘制了磁感应强度垂直分量对时间的导数 dbz/dt 在水平面上的扩散图，黑色虚线为异常体在地表投影的边界 当 y=0 1 S m1时，由于 x 方向的电导率比 y 方向的电导率小，响应产生沿 x 方向的拉伸，在早期 t=0 1 ms 时，响应极大值位 于 异 常 体 正 上 方，随 着 时 间 推 移 到 t=0 4 ms，沿响应拉伸的方向产生了两个对称的极大值 当 y=0 02 S m1，此时 x 方向和 y 方向的电导率相等，所以响应沿 x 方向和 y 方向的扩散速度相8542023，38(1)管清辰，等:航空瞬变电磁法三维各向异性数值模拟和响应特征研究(www progeophys cn)图 17D取不同值的航空瞬变电磁响应Fig 17Airborne transient electromagnetic response with different D等，呈现中心对称的形状 当 y=0 01 S m1时，由于 y 方向的电导率比 x 方向的电导率小，响应产生沿 y 方向的拉伸，由于电导率的减小，响应扩散速度增大，所以在 t=0 1 ms 时，响应已经延拉伸方向出现了两个极大值，随着时间推移到 t=0 4 ms，两个极大值的距离增加3 1 2垂直方向电导率的影响异常体水平方向电导率 x=y=0 02 S m1，异常体垂直方向电导率 z取01 S m1，002 S m1，0 005 Sm1，图 12 为该各向异性异常体的物理模型示意图 水平层状的沉积岩或者各种岩体在应力作用下发生断裂，形成水平方向的节理后，都会导致水平方向的电导率与垂直方向的电导率不同响应如图 13 所示，随着垂直方向的电导率 z的减小，响应衰减速度增大，同一时刻响应的幅值有轻微的减小 与之前水平方向电导率变化带来的影954地球物理学进展www progeophys cn2023，38(1)图 18y取不同值的航空瞬变电磁响应Fig 18Airborne transient electromagnetic response with different y响相比，垂直方向电导率变化对响应的影响较小3 1 3旋转角 S的影响围岩电导率 0 01 S m1，异常体三个主轴电导率分别为 x=0 02 Sm1，y=0 1 Sm1，z=0 02 S m1，旋转角 D=L=0，旋转角 S分别取0，30，60，90 图 14 为该各向异性异常体的物理模型示意图，相当于 3 1 1 节的各向异性模型沿 z轴发生了旋转，S不为零时，电导率主轴方向和坐标轴不再重合响应如图15 所示，在早期时刻 t=0 1 ms 时，响应极大值位于异常体正上方，由于 x 方向的电导率小于 y 方向的电导率，响应沿着电导率较小的 x 方向产生了拉伸 随着时间推移到 t=0 4 ms，响应逐渐扩散，沿早期拉伸的方向出现了两个对称的极大值 旋转角 S由0变化到30，60，90时，响应拉伸方向在水平面上进行了相应角度的旋转0642023，38(1)管清辰，等:航空瞬变电磁法三维各向异性数值模拟和响应特征研究(www progeophys cn)图 19z取不同值的航空瞬变电磁响应Fig 19Airborne transient electromagnetic response with different z3 1 4旋转角 D的影响围岩电导率 0 01 S m1，异常体三个主轴电导率分别为 x=0 02 Sm1，y=0 1 Sm1，z=0 02 S m1，旋转角 S=L=0，旋转角 D分别取0，30，60，90 图 16 为该各向异性异常体的物理模型示意图，相当于 3 1 1 节的各向异性模型沿 x轴发生了旋转，D不为零时，电导率主轴方向和坐标轴不再重合响应如图 17 所示，旋转角 D=0时与前面 S=0的情况是一样的 随着旋转角 D取到 30和60时，异常范围沿 y 方向增大 这是因为旋转角 D时是绕 x 轴旋转的，D=0时，异常体 y 方向的电导率等于 y D取 30和 60时，沿 y 方向的电导率将受到 y和 z共同作用，效果介于两者之间 D取30和 60时，响应极值沿 y 方向出现了轻微偏移，不再沿 x 轴完全对称 D=90时，此时 y 方向的电导164地球物理学进展www progeophys cn2023，38(1)图 20S取不同值的航空瞬变电磁响应Fig 20Airborne transient electromagnetic response with different S2642023，38(1)管清辰，等:航空瞬变电磁法三维各向异性数值模拟和响应特征研究(www progeophys cn)图 21D取不同值的航空瞬变电磁响应Fig 21Airborne transient electromagnetic response with different D364地球物理学进展www progeophys cn2023，38(1)率等于 z，跟 x 方向电导率 x相等，所以响应呈现出中心对称的形状3 2异常体为各向同性，背景围岩为各向异性3 2 1水平方向电导率的影响异常体电导率为01 S m1，围岩垂直方向电导率z=001 S m1，水平方向电导率 x=0 01 S m1，y取 0 03 Sm1，0 01 Sm1，0 003 Sm1，各向异性围岩的物理模型图如图 10 所示响应如图 18 所示，当 y=0 03 S m1时，由于y 方向的电导率比 x 方向的电导率大，响应产生沿 y方向的拉伸，在早期 t=0 1 ms 时，响应极大值位于异常体正上方，随着时间推移到 t=0 4 ms，沿响应拉伸的方向产生了两个对称的极大值 与前边异常体为各向异性时的情况不同，围岩为各向异性时，响应沿电导率更大的方向拉伸了，而不是电导率更小的方向 当 y=0 01 S m1，此时 x 方向和 y 方向的电导率相等，所以响应沿 x 方向和 y 方向的扩散速度相等，呈现中心对称的形状 当 y=0 003 S m1时，由于 x 方向的电导率比 y 方向的电导率大，响应产生沿 x 方向的拉伸，但在观测时间范围内，随着时间的推移，并没有出现两个极值点3 2 2垂直方向电导率的影响异常体电导率为 0 1 Sm1，围岩水平方向电导率 x=y=0 01 S m1，围岩垂直方向电导率 z取0 03 S m1，0 01 S m1，0 003 S m1，各向异性围岩的物理模型图如图 12 所示响应如图 19 所示，随着垂直方向的电导率 z的减小，响应衰减速度增大，同一时刻响应的幅值减小 由于围岩在构造尺度上远大于异常体，所以围岩垂直方向电导率比之前异常体垂直方向电导率对响应的影响更明显3 2 3旋转角 S的影响异常体电导率为 0 1 Sm1，围岩三个主轴电导率分别为 x=0 01 Sm1，y=0 03 Sm1，z=0 01 S m1 旋转角 D=L=0，旋转角 S分别取 0，30，60，90，各向异性围岩的物理模型图如图 14 所示响应如图20 所示，在早期时刻 t=0 1 ms 时，响应极大值位于异常体正上方，随着时间推移到 t=0 4 ms，响应逐渐扩散，沿早期拉伸的方向出现了两个对称的极大值 旋转角 S由 0变化到 30，60，90时，响应拉伸方向在水平面上进行了相应角度的旋转3 2 4旋转角 D的影响异常体电导率为 0 1 Sm1，围岩三个主轴电导率分别为 x=0 01 Sm1，y=0 03 Sm1，z=0 01 S m1 旋转角 S=L=0，旋转角 D分别取 0，30，60，90，各向异性围岩的物理模型图如图 16 所示响应如图 21 所示，旋转角 D=0时与前面 S=0的情况是一样的 旋转角 D取 30和 60时，可以看到响应在 y 方向出现了明显的偏移，不再沿 x轴对称 t=0 4 ms，D=30时，下边的极大值明显小于上边的极大值 D=60时，下边极值完全消失，只有上边存在一个极值 D=90时，此时 y 方向的电导率等于 z，跟 x 方向电导率 x相等，响应平面扩散图又呈现出中心对称的形状4结论本文从时间域的麦克斯韦方程出发，采用时域交错采样有限差分法，推导地下电磁场满足的离散的时间分步迭代方程 用 Fortran 语言编写程序实现了航空瞬变电磁法任意各向异性介质三维正演 与已有的一维和二维的软件进行对比，验证了程序的正确性 然后设计典型的各向异性三维理论地电模型，对响应的异常特征进行分析，探究了三维地电模型的各向异性参数对航空瞬变电磁响应的影响规律异常体为各向异性时，响应沿水平电导率更小的方向产生拉伸;围岩为各向异性时，响应沿水平电导率更大的方向拉伸 随着垂直方向的电导率的减小只会导致响应衰减速度增大，不会造成响应拉伸或者偏移 绕 z 轴的旋转角 S变化会导致响应拉伸方向在水平面上进行相应角度的旋转 绕 x 轴的旋转角 D变化会导致响应在 y 方向出现偏移，不再沿x 轴对称 可以看到，介质为各向异性时与各向同性介质的航空瞬变电磁响应有着明显的不同，不同的各向异性参数组合可能会导致航空瞬变电磁响应产生相应的拉伸，旋转以及偏移 这给航空瞬变电磁资料的解释提供了新的视角，也为航空瞬变电磁法任意各向异性介质三维反演工作打下了基础致谢感谢审稿专家提出的修改意见和编辑部的大力支持!eferencesAuken E，Christiansen A V，Kirkegaard C，et al 2015 An overview of4642023，38(1)管清辰，等:航空瞬变电磁法三维各向异性数值模拟和响应特征研究(www progeophys cn)a highly versatile forward and stable inverse algorithm for airborne，ground-basedandboreholeelectromagneticandelectricdataExploration Geophysics，46(3):223-235Avdeev D B，Kuvshinov A V，Pankratov O V，et al 1998 Three-dimensional frequency-domain modeling of airborne electromagneticresponses Exploration Geophysics，29(1-2):111-119Guan Q C 2019 25 dimensional modeling of transient electromagneticdata for anisotropic media Master s thesis(in Chinese)Beijing:China University of Geosciences(Beijing)Huang W，Ben F，Yin C C，et al 2017 Three-dimensional arbitrarilyanisotropicmodelingfortime-domainairborneelectromagneticsurveys Applied Geophysics，14(3):431-440Liu Y H，Yin C C 2014 3D a