我国地域广袤，地质条件复杂，矿产资源与水资源形势严峻，加大勘查力度已是迫在眉睫.航空瞬变电磁法(Airborne Transient Electromagnetics，简称ATEM)能够进入山区、沙漠、沼泽、森林等地形条件复杂地区进行矿产资源普查、地下水探查与构造地质填图等工作，以其速度快、成本低、可大面积覆盖等优势能够很好地适应我国山区大面积找矿、干旱地区大面积找水、区域地质填图等各方面的需求，在我国资源探测领域有着广泛的应用前景(雷栋等，2006；张昌达，2006).目前对于ATEM的数据的计算与解释大多是基于平坦地形的理想情况，然而在数据解释时若不考虑地形效应则容易将地形影响错误的解释为地下异常体并将掩盖地下真实目标体的响应，从而造成数据解释的较大误差，进一步导致勘探成本的增加和系统数据的浪费.因此，研究起伏地形条件下航空瞬变电磁系统的正演模拟对实际应用与数据解释都有着重要的指导意义.

由于航空电磁方法(特别是直升机航空电磁系统)大部分都是在地形复杂地区进行作业，所以对于地形效应的研究也越来越受到国内外学者的重视.Liu和Becker(1992) 利用边界元方法对于二维条件下的斜坡、山脊、山谷几种典型地形的频率域航空电磁(Airborne Electromagnetics，简称AEM)响应进行了计算；Newman和Alumbaugh(1995，1998)利用基于交错网格有限元方法实现了频率域AEM的三维正演，并模拟了简单二维地形条件下频率域AEM响应；Sasaki和Nakazato(2003) 采用有限差分方法对于二维地形条件下的频率域AEM数据进行了模拟，并指出山脊地形其在山顶处产生高阻异常而在山脚转折处产生低阻异常响应.王卫平等(2015) 利用频率域有限差分方法模拟了二维、三维起伏地形条件下的频率域AEM响应，指出地形对于频率域航空电磁响应具有较大影响，特别是在地形局部变化较大的地方. 张博，殷长春等(2016) 采用非结构有限元方法对山谷、山脊简单地形条件下三维频域与时域航空电磁模型进行了数值模拟，计算结果表明地形对于航空电磁数据的影响主要集中在高频段与时间道早期，且地形突变区域电磁响应也出现快速变化.从以上研究中不难发现，目前对于地形效应的模拟大多针对频率域航空电磁数据，而对于航空瞬变电磁数据的研究，特别是三维地形条件下航空瞬变电磁数据的研究相对较少；相比频率域航空电磁法，时间域航空电磁法发射功率大，勘探深度深，具有其独特的优势.因此，研究三维航空瞬变电磁地形效应的影响是十分必要的.

瞬变电磁的正演是研究瞬变电磁响应规律的有效途径，也是反演方法的前提与基础.随着社会经济的迅速发展和科学技术的全面进步，航空电磁法也逐步由一维模拟向三维正演发展.J.O.Barongo(1998) 、罗廷钟等(2004) 、李永兴等(2010) 对层状大地条件下的时间域航空电磁一维正演做了研究，并细致分析了航空瞬变电磁的响应规律，为时间域航空电磁系统的设计方案提供了理论依据；殷长春等(2013) 对ATEM系统全时瞬变响应多分量正演模拟做了相应研究，并将其推广到三维模型.目前关于瞬变电磁三维正演常用的方法包括有限单元法(Finite-element techniques)(Sugeng，1998; Borner，2008; Li，2011)，有限体积法(Finite volume techniques)(Haber et al.，2007; Oldenburg et al.，2008，2012; Yang D et al.，2012)，有限差分方法(Finite-difference techniques)(Wang 和 Hohmann，1993; Commer 和 Newman，2004; Newman和Commer，2005; 许洋铖等，2012; 孙怀凤等，2013)，积分方程法(Integral equation techniques)(Xiong，1992; Zhdanov et al.，2006; Cox et al.，2010).李建慧等(2011，2013)利用矢量有限单元法对大回线源激发的瞬变电磁场进行了三维数值模拟，这对三维航空电磁数据有限元正演模拟有着重要的指导意义.Oldenburg等(2012) 采用有限体积法对瞬变电磁三维正演计算进行了研究，实现了多激励源条件下的三维正演并适用于时间域航空电磁系统.Cox等(2010) 利用积分方程法实现了对航空电磁数据的三维正演模拟.Wang和Hohmann(1993) 提出了一种基于Du Fort-Frankel有限差分格式的三维瞬变电磁场数值模拟算法，利用该方法可以显式的得到任意时间点的电磁场分量；许洋铖等(2012) 将有源有限差分方程转换为无源有限差分方程，实现了发射电流为任意波形时，全波形时间域航空电磁响应的三维数值计算；孙怀凤等(2013) 等对Wang和Hohmann提出的方法进行改进，实现了矩形回线源瞬变电磁三维时域有限差分正演研究，并在激发电流源的计算中考虑关断时间.由于时间域电磁场的复杂性，以上研究工作多以简单模型为主，且基于平坦地形理想的情况，对于地形效应研究相对较少.现有的少数对于地形效应的研究，大多采用有限元法与边界元法(Liu和Becker，1992; Chen et al.，1998; Mitsuhata，2000; Franke et al.，2007; Li和Pek，2008; Nam et al.，2007).有限单元法网格剖分自由，对地形的模拟程度较高，然而对于三维正演特别是拥有海量数据的航空电磁法三维正演模拟，有限单元法因其对单元网格模拟需要占用大量的内存，计算效率大大受限.时域有限差分方法的计算效率与计算精度较高，并且具有天然的并行性，是航空瞬变电磁三维正演模拟较为理想的方法.

本文基于交错网格的时域有限差分算法，实现了对三维地形起伏条件下航空瞬变电磁数据的正演模拟.通过计算典型地形模型与含异常体的复杂模型的ATEM响应，给出多测道曲线及视电阻率断面曲线图，通过分析地形效应的影响特征与变化规律，对地形的影响有更加直观的认识.由于野外实际作业时地形效应还受到多种因素的影响，本文以山脊、山谷模型为例，讨论了地形的尺寸参数、电性参数、飞行轨迹与飞行高度等参数变化对于ATEM数据的影响情况，并得出有价值的结论.以上研究成果，为有效实现航空电磁数据的地形校正提供有价值的理论依据，对提高航空瞬变电磁法的解释精度具有重要意义.

本文时域有限差分算法采用Yee氏晶胞格式对电磁场各分量进行空间离散，利用差分代替微分的方式实现空间域和时间域偏微分方程的数值求解.为了进一步减少网格数量提高计算效率，采用交错网格进行模型剖分，在地形起伏范围内加密网格，其余区域适当增大网格，从而能够更加准确的模拟地形起伏进一步提高计算精度与计算效率.对于瞬变电磁激励源的加入，采用孙怀凤等(2013) 的处理方式，将电流密度加入麦克斯韦方程组的安培环路定理方程.对于航空系统进行模拟时，模型中大部分范围内为空气介质，发射源放置于空气层介质中，计算时将空气层作为模型的一部分，修正发射源与接收点的高度，对模型的六个面均采用Dirichlet边界条件(第一类边界条件)进行边界处理；并采用基于OpenACC应用程序编写接口，利用GPU设备对计算程序进行了并行化设计，有效提高计算效率.

在瞬变电磁勘探中，一般忽略位移电流(朴化荣，1990)，因而在各向同性的均匀导电介质中，忽略位移电流后的无源媒质中的麦克斯韦方程组：

(1a)

(1b)

(1c)

(1d)

其中E为电场强度，B为磁感应强度，H为磁场强度，σ为电导率，t为时间.

为了构成可以适用于FDTD的显式的时间迭代格式，根据阻尼波动方程与扩散方程的相互转化关系，人为引入虚拟位移电流项来构建显示的时域有限差分格式(Oristaglio和Hohmann，1984)，由此方程(1b)变成

(2)

其中，γ称为虚拟介电常数，γ的取值需要满足一定的稳定性条件.

在有源媒质中，Maxwell方程组中的方程(1b)必须包含源的电流项，为了实现回线源的加入，将矩形回线源电流密度加入麦克斯韦方程组的安培环路定理方程，则将以上方程(2) 修改为(孙怀凤等，2013)

(3)

其中J_s表示源的电流密度.

本文采用FDTD计算中应用最为广泛的Yee晶胞格式进行网格离散，这样使得每一个磁场(电场)分量周围都有四个电场(磁场)分量所包围，形成电磁场的交替迭代模式，满足麦克斯韦方程组的差分计算要求.将电流密度施加到与电场水平分量重合的Yee元胞的棱边上，采用梯形波发射，进一步获得回线源情况下的电磁场值.

航空瞬变电磁系统为大面积扫面测量，采集范围广、模型尺寸较大.为了利用有限的计算资源对于较大尺寸的模型进行模拟，时域有限差分一般采用非均匀网格剖分，在有源区域及异常体分布区域加密网格，其余区域则适当增加网格尺寸，从而提高计算效率.非均匀剖分时电场计算公式中的电导率与均匀剖分时的计算不同，需要按照每个网格的体积来计算电场采样点处相邻四个网格的贡献值，由此得到电场迭代的FDTD形式如下(孙怀凤，2013)：

(4)

(5)

(6)

与地面瞬变电磁不同，航空瞬变电磁在建模时需要考虑包含空气层的全空间模型，本文将整个模型的外边界共6个面全部施加Dirichlet边界条件，每次迭代过程中将电场的切向分量和磁场的垂向分量强制赋零.在实际计算过程中为了避免由于电磁场分量强制赋零所产生的电磁波减弱的改变对于中心计算区域的影响，计算时必须保证模型网格剖分的规模足够大，可以忽略这种影响.

对于航空系统进行模拟时，模型中大部分范围内为空气介质，发射源放置于空气层介质中，计算时将空气层作为模型的一部分，修正发射源与接收点的高度.空气是一种近似绝缘的介质，电磁场在空气中转变为波动场.FDTD在波动场中的计算对时间离散步长要求非常短，如果以扩散场中对时间离散步长的要求来统一波场中的时间离散步长，将会导致整个计算过程非常耗时.根据已有的近似情况与计算经验(Um et al.，2012)，给定空气电导率为一个足够小的值而不赋零，这样既能够使计算过程中地质体与空气介质采用相同的迭代方程又能够保证计算精度.本文给定空气电导率为10^-6S·m^-1.

针对三维FDTD问题，其时间步长需要受到Courant-Friedrich-Levy(CFL)(Courant R. et al.，1967) 稳定性条件的限制，均匀剖分时，时间步长限制为如下形式：

(7)

其中δ表示网格尺寸，由于虚拟介电常数γ是人为引入的参量，其取值范围也有一定的限制，将上式(7) 进行改写后即可得到虚拟介电常数的限制

(8)

对于一次场的模拟，Wang和Hohmann(1993) 将γ的值给定为真实的介电常数，然后采用等间隔剖分依据公式(7) 给出时间迭代步长.然而，一次场在关断之前的大部分时间内都处于稳定状态，电磁场值变化较小，因此对于有限差分模拟来说并不需要如此小的时间迭代步长，时间迭代步长过小将导致一次场计算时间过长，严重降低计算效率.为了进一步提高计算效率，减短计算时间，本文中对于一次场的计算，采用变化时间迭代步长方法进行计算.以梯形波为例，对于梯形波的上升沿与下降沿阶段，电磁场值变化剧烈，给定其时间迭代步长为1 ns；对于梯形波的持续时间阶段，通过给定递增因子与递减因子的方式，先增大再缩小的给定时间迭代步长.为了进一步避免时间迭代步长无限制的增大，需要给定一个一次场最大时间迭代步长，通过大量的计算发现，最大的时间迭代步长取值可以在公式(7) 的基础上适当增大量级，也就是

(9)

此时γ取值为真空介电常数，k的取值范围为50~100，可根据计算需要精度范围进行选取.

现有的航空瞬变电磁观测系统，进行的是大面积扫面测量，这就要求三维航空电磁模拟需要随着发射源的移动，不断改变源的位置单独做一次正演计算，这样的计算工作量是十分巨大的，需要占用大量计算机内存与计算时间.然而时域有限差分方法具有天然的并行性，采用并行算法能够极大地提高计算效率.本文借助OpenACC应用程序编写接口，采用CPU+GPU异构并行策略对程序进行了并行化处理，有效提高计算效率，大大缩短计算时间.瞬变电磁三维FDTD程序的时间主要花费在六个电磁场分量的更新迭代计算上，利用CPU+GPU异构并行模式可以使程序中计算密集的迭代部分在GPU中进行，而逻辑控制部分和数据的读写部分则由CPU处理，每次由CPU向GPU提交一组循环计算，GPU计算完成后将结果返回给CPU，然后再次接受CPU分配的下一组循环计算.(文中GPU采用 NVIDIA Tesla K20工作站，卡内配置2688个CUDA计算核心).

由麦克斯韦方程组可以导出均匀半空间情况下回线源航空瞬变电磁的磁场垂直分量表达式如下(Nabighian，1988)

(10)

式中I₀为电流强度，a为发射线圈半径，μ₀为真空磁导率，J₀与J₁分别为零阶与一阶第一类贝塞尔函数.采用D. Guptasarma和 B. Sing(1997) 提出的改进的线性数字滤波方法求取以上的汉克尔型积分问题，并采用其给出的滤波系数计算(10) 式中的汉克尔积分，进一步提高计算精度.对于时间域的电磁场，可以由频率域电磁场通过正、余弦变换得到

(11)

为了验证本文所给出的航空瞬变电磁三维FDTD程序的正确性，以均匀半空间模型为例，与一维滤波解进行对比.采用中心回线方式，发射磁矩为9×10⁵Am²，发射高度与接收高度均为30 m，均匀大地电阻率为100 Ωm.计算结果如图 1所示，FDTD解与滤波解的磁感应强度曲线吻合较好，在晚期稍有差异，最大误差在3.5%左右，整体平均误差在0.6%以下.

图 1

Fig. 1

图 1 均匀半空间模型FDTD解与滤波解计算结果对比曲线 (a)磁感应强度衰减曲线;(b)磁感应强度相对误差. Fig. 1 Comparison between the FDTD solution and the analytical soltion of a homogeneous half-space model (a)Decay curves of BZ;(b)Relative error of BZ decay curves.

为了进一步验证本文算法对于起伏地形模型模拟的正确性，将本文计算结果与张博，殷长春等(2016) 使用有限元方法计算的三维山脊模型的时间域航空电磁响应进行对比.模型示意图如图 2a所示，山脊高20 m，顶部与底部宽度分别为60 m和220 m，模型电阻率为100 Ωm，收发装置采用水平共面装置，收发距为10 m，飞行高度为30 m.选取两个时间道(t=1.6×10^-5s、2.51×10^-3s)进行对比，图 2b为磁感应强度曲线对比图，图 2c为相对误差曲线.对比结果表明，本文FDTD算法与前人计算结果吻合较好，最大误差不超过2.1%，整体平均误差在0.8%以下，由此证明本文方法对于地形效应的模拟是正确可靠的.

图 2

Fig. 2

图 2 山脊模型FDTD解与殷长春等(2016) 模拟结果对比曲线 (a)模型示意图;(b)计算结果对比曲线;(c)两时间道相对误差曲线. Fig. 2 Comparison between the FDTD solution and the analytical soltion of a homogeneous half-space model (a)A trapezoid hill model;(b)Comparison of the FDTD solution with YIN et al.，(2015) ;(c)Relative errors.

由于进行有限差分正演的单元网格均为正交的矩形，当模拟地形起伏时，难以直接模拟弯曲或者倾斜的地形表面，所以一般采用“台阶状”网格逼近地形的起伏变化(如图 3)，在地形起伏区域加密网格，可以有效提高对于地形界面的刻画精度.由于直接观测磁感应强度B比测量∂B/∂t有诸多优点(Raiche，1983；Spies和Eggers，1986；白登海等，2003)，文中对磁感应强度垂直分量B_z进行计算，并利用基于B_z的全域视电阻率定义方法(赵越等，2015)对ATEM视电阻率进行计算，该方法有效实现了全时域、全空域的视电阻率定义、并且避免了利用∂B_z/∂t定义在早期出现undershoot/overshoot的假极值干扰(Raiche，1983)，从而能够更加准确反映起伏地形对于ATEM数据的影响.

图 3

Fig. 3

图 3 网格剖分示意图(X-Z平面) Fig. 3 Cross-section view in the X-Z plane at Y=0

为了考察地形起伏对于航空瞬变电磁数据的影响，以下选取几种典型地形模型进行分析处理.模型计算所需变量如表 1所示，计算采用中心回线装置(以下模型若无特殊说明均采用以下参数进行计算).

表 1(Table 1)

表 1 计算参数设置 Table 1 Parameter Settings 参数 值(单位) 参数 值(单位) 模型总网格数： 221×221×200 最小网格长度 10(m) 测线长度 800(m) 测点数 21 发射回线边长： 30(m) 发射电流： 1000(A) 采样时间 [10-5，10-2](s) 采样数目 40

表 1 计算参数设置 Table 1 Parameter Settings

为了考察山脊地形对于航空瞬变电磁场的影响情况，给出模型如图 4a所示，山脊的顶部面积为20 m×20 m，底部矩形面积为200 m×200 m，山脊高度为100 m，山脊与地下介质电阻率为100 Ωm，发射与接收高度均为20 m，发射磁矩为9×10⁵ Am²，采样时间为10^-5s~10^-2s，采用对数等间隔选取40个采样点，测线布置图如4b所示，主剖面为Y=0测线.

图 4

Fig. 4

图 4 山脊地形正演模型 (a)山脊模型示意图;(b)山脊模型XY平面视图与测线布置图;(c)山脊模型磁感应强度多测道曲线图(主剖面); (d)山脊模型视电阻率断面图(主剖面). Fig. 4 3-D trapezoidal hill model (a)Sketch map of the 3D trapezoidal hill model;(b)X-Y plane view and survey line arrangement;(c)Bz response for the 3-D trapezoidal hill model(Y=0) ;(d)Apparent resistivity contours of the 3-D trapezoidal hill(Y=0) .

图 4中结果表明，山脊地形的起伏对于ATEM探测存在显著的影响.由图 4c中多测道曲线图可以看出，山脊地形的影响主要集中在时间道早期，并且山脊地形顶部响应曲线下凹，但在山脊底部转折处曲线则成上凸趋势；随着时间的不断推移，曲线逐渐变为平缓，地形影响逐渐减弱.转换为视电阻率断面图可知，山脊地形在顶部产生高阻异常影响，但在底部转折处则产生低阻异常，这与Sasaki和Nakazato(2003) 计算所得结论一致.由此可知，地形对于航空瞬变电磁响应具有显著的影响，且影响主要集中于时间道早期.

为了进一步分析山脊地形对于地下异常体的影响情况，设计山脊地形条件下赋存一个低阻异常体的模型(如图 5a)，山脊地形参数与图 4a相同，低阻异常体的埋深80 m，异常体尺寸150 m×150 m×50 m，电阻率为1 Ωm.图 5b、c分别为带异常体情况下山脊地形的ATEM正演响应结果及视电阻率断面，图中结果显示，B_z响应曲线在时间道早期与纯山脊地形响应基本一致，主要反映地表起伏的影响，随着时间的进一步推移，地下异常体响应逐渐反映出来，地形影响进一步减弱.从视电阻率断面图中可以清晰的看到山脊地形顶部所产生的高阻异常与山脚处的低阻响应，地下的低阻异常体由于受到地形的影响，异常形态由原本的单个异常体变为两个独立异常体，偏离实际模型.由此说明，山脊地形对于地下异常体影响显著，若不考虑地形影响，直接将起伏地形按照平坦地面进行解释将造成较大误差.

图 5

Fig. 5

图 5 山脊地形带异常体正演模型 (a)山脊带异常体模型示意图;(b)磁感应强度多测道曲线图;(c)视电阻率断面图. Fig. 5 3-D trapezoidal hill model with an abnormal body (a)Sketch map of the 3D trapezoidal hill model with an abnormal body;(b)Bz response for the 3-D trapezoidal hill model with an abnormal body;(c)Apparent resistivity contours(Y=0) .

给出山谷模型如图 6a、b所示，山谷的顶部矩形面积为200 m×200 m，底部矩形面积为20 m×20 m，山谷深度为100 m，地下介质电阻率为100 Ωm，发射与接收高度均为50 m，发射磁矩为9×10⁵Am².

图 6

Fig. 6

图 6 山谷地形正演模型 (a)山谷模型示意图;(b)山谷模型XY平面视图与测线布置图;(c)山谷模型磁感应强度多测道曲线图(主剖面); (d)山谷模型视电阻率断面图(主剖面). Fig. 6 3-D trapezoidal valley model (a)Sketch map of the 3D trapezoidal valley model;(b)X-Y plane view and survey line arrangement;(c)Bz response for the 3-D trapezoidal valley model(Y=0) ;(d)Apparent resistivity contours of the 3-D trapezoidal valley(Y=0) .

由模型处理结果图 6c、d可以看出，山谷地形的起伏对于ATEM探测存在显著的影响，且该影响主要体现在早期.图 6c中的早期时刻，磁感应强度曲线形态为中心上凸而两侧地形转折处下凹的曲线，与山脊地形响应相反，随着时间的推移，多测道曲线趋于平缓，为近似的水平直线，表现为地下均匀大地响应与模型设计一致.视电阻率断面图在早期中心位置呈现出一个低阻异常，但在两侧地形转折处为高阻异常响应，与山脊视电阻率响应恰恰相反；随着二次场的向下传播，电阻率逐渐增大并趋于均匀大地电阻率.

同样设计山谷地形条件下含低阻异常体的模型(如图 7a)，山谷模型参数不变，低阻异常体的埋深50 m，异常体尺寸150 m×150 m×50 m，电阻率为1 Ωm.计算结果如图 7b、c所示，山谷地形影响主要集中在早期，随着时间的不断推移，山谷地形影响逐渐减弱，地下异常则逐渐显现，B_z响应曲线在晚期曲线明显上凸，对应地下低阻异常体的响应；视电阻率断面中也可以明显看出在时间道晚期存在一个低阻异常影响，这与模型设计一致.

图 7

Fig. 7

图 7 山谷地形带异常体正演模型 (a)山谷带异常体模型示意图;(b)磁感应强度多测道曲线图;(c)视电阻率断面图. Fig. 7 3-D trapezoidal valley model with an abnormal body (a)Sketch map of the 3D trapezoidal hill model with an abnormal body;(b)Bz response for the 3-D trapezoidal valley model with an abnormal body;(c)Apparent resistivity contours(Y=0) .

实际地形是复杂多变的，我们探测的目的是埋藏于起伏地形下的目标异常体，为了探讨地形对于地下含异常体的ATEM响应影响，设计了峰谷联合起伏地形下的复杂模型如图 8所示.模型参数见表 2.

图 8

Fig. 8

图 8 峰谷联合纯地形正演模型 (a)峰谷联合模型示意图;(b)XY平面视图与测线布置图;(c)磁感应强度多测道曲线图;(d)视电阻率断面图. Fig. 8 3-D trapezoidal complex terrain model (a)Sketch map of the 3D complex model;(b)X-Y plane view and survey line arrangement;(c)BZ response for the 3-D complex terrain model(Y=0) ;(d)Apparent resistivity contours of the 3-D compelx teerain mode(Y=0).

表 2(Table 2)

表 2 复杂模型参数 Table 2 Parameters of the complex model 模型 尺寸(长×宽×高)(m) 电性参数(Ωm) 山峰 150×150×70 100 山谷 150×150×70 106 山谷底部异常体 150×200×70(埋深170) 2 山脊底部异常体 100×100×60(埋深70) 0.3

表 2 复杂模型参数 Table 2 Parameters of the complex model

图 8显示的为“两山夹一谷”复杂地形模型正演结果，接收高度与发射高度均为30 m，得到磁场感应强度变化曲线及视电阻率断面如图 8c，d所示.由图 8c中的多测道曲线图可以清楚的看到ATEM电磁响应随着时间推移的变化过程，在时间早期，B_z曲线主要受到地形的影响，呈现出起伏状曲线形态，山峰所在位置曲线下凹、山谷对应位置曲线则上凸，曲线形态与地表起伏形状相反.随着时间的推移，地形影响逐渐减弱，多测道曲线逐渐趋于水平，指示地下均匀介质.图 8d视电阻率断面图中对于地形效应的影响则反映的更加直观，图中可得山谷处呈现低阻异常响应而山脊处则表现为高阻异常，并且在山峰的山脚转折处为低阻异常.由此可见峰谷联合模型地形对于ATEM响应的影响很大，且复杂地形模型的整体响应结果可以看做是各简单地形的组合.

为了进一步考察复杂地形对于地下异常体的影响情况，设计如图 9a模型，网格剖分示意图如图 9b所示，得到三维正演结果如图 10所示.

图 9

Fig. 9

图 9 三维航空正演复杂模型 (a)复杂模型工作示意图;(b)复杂模型网格剖分示意图. Fig. 9 3D terrain ground with three abnormal body for ATEM modeling (a)Sketch map of the complex model;(b)The grid subdivision schemes of the complex model.

图 10

Fig. 10

图 10 山谷地形正演模型 (a)复杂模型磁感应强度多测道曲线图;(b)复杂模型视电阻率断面图; (c)复杂模型三维电阻率切片图;(d)X-Y方向三维电阻率切片图. Fig. 10 Synthetic model of square frustum vally (a)Sketch map of ATEM system and the synthetic model;(b)Apparent resistivity contours of the complex model; (c)Slice map of the complex model;(d)Slice map of the 3-D complex model in X-Y dircetion.

从图 10a中的多测道曲线图可以清楚的看到ATEM电磁响应随着时间推移的变化过程，在早期，B_z曲线主要受到地形的影响，呈现出起伏状曲线形态，曲线起伏形态与纯地形形态一致，进一步说明地形影响体现在时间道早期，随着时间的推移，底部异常体逐渐显现，曲线呈现上凸趋势.图 10b的视电阻率断面图早期变化为地形起伏的影响，山峰主要显示为一高阻异常，山谷则主要为一低阻异常，然而由于受到山峰模型的压制，山峰底部低阻异常体在电阻率断面上一分为二，显示为两个尺寸更小的低阻异常；山谷底部异常受到地形与周边异常影响在电阻率图中的表现与实际大小与位置并不相符.由此可知地形对于地下介质影响明显，对于野外资料的处理解释地形效应不可忽略.图 10c、d为视电阻率三维立体切片图，图中反映出复杂模型的地下形态以及地形效应的影响过程，对地形效应的影响过程表现更加清晰与形象.

以上通过几种典型地形模型的正演计算，对于地形起伏变化对于ATEM数据的影响有了进一步的认识，然而实际地形情况是复杂而多变的，更为复杂的地形起伏情况下的响应规律还有待做进一步深入研究.

上文针对几种典型地形的三维瞬变响应进行模拟，并给出其全域视电阻率曲线，从而对地形效应的影响有了更加直观的认识.然而在ATEM的实际测量中，地形效应还受到多种因素的影响，其中包括地形的尺寸参数变化的影响、地形的电性参数变化的影响(例如风化、冲蚀等作用使基岩出露造成地表的电性不均匀、雨后地表的电性变化等)、飞行的高度与飞行轨迹变化的影响等等.由此，下文通过设计大量模型进一步分析不同参数变化对于地形效应的影响情况.

由于单元网格为正交的矩形，为了阐述形象与计算方便，以三维直角坐标系作为基准，研究地形在X，Y，Z三个方向的地形尺寸变化对于ATEM响应的影响.以山谷地形为例(如图 11)，令坐标的原点位于模型中心，模型沿坐标轴X，Y，Z方向的尺寸参数分别为lx、ly、lz，测线沿X轴方向.以山谷模型为例研究X，Y，Z方向尺寸变化的影响(主剖面为Y=0测线)，山谷的顶部矩形面积为200 m×200 m，底部矩形面积为20 m×20 m，山谷深度为100 m，地下介质电阻率为100 Ωm，发射与接收高度均为50 m，发射磁矩为9×10⁵ Am².

图 11

Fig. 11

图 11 山谷模型三维示意图 Fig. 11 3D view of the trapezoidal valley model

保持Y、Z方向尺寸不变，将X方向尺寸lx的取值由200 m变化到400 m(步长为100 m)，绘制B_z曲线对比图及视电阻率断面对比图(图 12).计算结果表明，随着X方向尺寸的增大，地形对ATEM响应的影响范围也在逐渐的增大，但主要变化都集中在地形起伏区域.多测道曲线图中，不同宽度的山脊地形的磁感应强度曲线随时间的变化趋势是一致的，且地形影响主要集中在早期，晚期各曲线逐渐重合，表现为地下均匀大地的响应.多测道曲线图中B_z曲线形态与地形变化一致，随着X方向尺寸的不断增大，B_z曲线也逐渐增宽，但极值基本保持一致.视电阻率断面对比图中，地形尺寸越宽，低阻异常横向影响的范围也越大.由此说明X方向尺寸变化与ATEM电磁响应呈正相关(以下测量剖面均为Y=0主剖面).

图 12

Fig. 12

图 12 X方向尺寸变化对比图 (a)Bz响应对比图;(b)视电阻率断面对比图. Fig. 12 Comparison of X direction size change (a)Comparison of Bz-field;(b)Comparison of apparent resistivity contours.

保持X、Z方向尺寸不变，将Y方向尺寸ly的取值由200 m变化到400 m(步长为100 m)，绘制B_z曲线对比图及视电阻率断面对比图(图 13).图 13结果显示，Y方向尺寸的变化对于ATEM电磁响应影响较小，不同宽度B_z响应差别较小且几乎重合，仅在极值点处有微小差异；视电阻率断面对比图中不同地形尺寸的视电阻率断面图形态基本一致.由于测线沿X方向布置，故推断地形影响与测线布置位置有关，为了进一步比较X、Y方向尺寸变化的影响程度，将同样尺寸不同方向的响应曲线进行对比如图 14所示，由图中可以看出，两曲线存在明显差异，且主要集中在早期.图 14中沿测线X方向的曲线(红色)变化范围比Y方向(绿色)更宽，且尺寸越大两曲线B_z响应差异越大，由此说明地形宽度变化对于ATEM响应的影响主要受到测线方向的控制，沿测线方向的地形尺寸变化与ATEM响应成正相关.据此，实际工作中可根据实际地貌与目标选取合适的观测方向.

图 13

Fig. 13

图 13 Y方向尺寸变化对比图 (a)Bz响应对比图;(b)视电阻率断面图对比图. Fig. 13 Comparison of Y direction size change (a)Comparison of Bz-field;(b)Comparison of apparent resistivity contours.

图 14

Fig. 14

图 14 X、Y方向宽度变化对比图 (a)lx=ly=300;(b)lx=ly=400. Fig. 14 Comparison of size change in X and Y direction

保持X、Y方向尺寸不变，将Z方向尺寸lz的取值由100 m变化到300 m(步长为100 m)，绘制B_z曲线对比图及视电阻率断面对比图(图 15).图 15结果表明，随着Z方向尺寸的增大，地形对于ATEM数据的影响强度与影响范围也逐渐增大.在多测道曲线图中，随着地形尺寸的不断加深，地形对场的整体影响增强，异常幅度也相应增大.视电阻率断面对比图中，地形尺寸越深，中心低阻异常影响范围也越深，相应的电阻率值也越低，转折处电阻率则越高，这与地形整体的变化趋势是一致的.可见Z方向尺寸变化与ATEM电磁响应也呈正相关.

图 15

Fig. 15

图 15 Z方向尺寸变化对比图 （ａ）Bz响应对比图；（ｂ）视电阻率断面图对比图. Fig. 15 Comparison of Z direction size change (a)Comparison of Bz-field;(b)Comparison of apparent resistivity contours.

为了进一步考察地形的电性参数变化对于ATEM响应的影响情况，以山脊模型为例(图 16a)，保持均匀大地电阻率100 Ωm不变，计算山脊电阻率ρ=10 Ωm、ρ=100 Ωm、ρ=1000 Ωm时的三维ATEM瞬变响应，计算结果如图 16(b、c)所示.

图 16

Fig. 16

图 16 地形电性变化对比图 (a)山脊模型示意图;(b)Bz响应对比图;(c)视电阻率断面图对比图. Fig. 16 Comparison of terrain apparent resistivity change (a)Sketch map of the 3D trapezoidal hill model with apparent resistivity change; (b) Comparison of Bz-field; (c) Comparison of apparent resistivity contours.

计算结果表明，地形电阻率变化对于ATEM响应有着显著的影响，当山脊电阻率为ρ=10 Ωm与ρ=100 Ωm时，ATEM响应截然不同.由上文中的计算结果可知，当地形电阻率与地下介质电阻率一致时，山脊地形在顶部产生一个高阻异常，但在其底部地形转折处产生的是低阻异常；然而当山脊地形电阻率发生改变时，这两种异常区域的范围与强度也发生改变.由图 16b中多测道对比图可以看出，当山脊电阻率减小为10 Ωm时，由于山脊电阻率与空气层及地表的电性差异增大，导致山脊底部转折处B_z响应进一步增强，而山顶处B_z响应则进一步减弱，由此转换到视电阻率断面(图 16c)中则表现为低阻异常的范围与强度明显增强，顶部高阻异常受到压制，整体地形模型呈现低阻异常影响.同理当山脊电阻率增大到1000 Ωm时，底部转折处的低阻异常范围与强度减小，顶部高阻异常进一步增强，整体模型主要体现为顶部的高阻异常.

由此可知，地形电阻率对于ATEM响应影响显著，同一地形模型不同电阻率情况下的ATEM响应可能完全不同，实际测量中应结合实际地质资料进行处理解释，不能仅凭电阻率表征或者地质地貌对地形效应进行判断.

在野外实际测量中，地形的电性还受到多种地质作用及自然因素的影响，例如风化、冲蚀等作用使基岩出露造成地表的电性不均匀、雨后地表的电性发生变化等等，这些因素都有可能造成地形中地表的电性发生改变，为了进一步了解地表的电性变化对于ATEM数据的影响，以山脊模型为例，通过改变山脊表层电阻率及表层电阻率的厚度，分析表层的电性变化对于ATEM数据的影响情况.

图 17a为山脊模型表层电阻率变化示意图，计算模型时按照台阶状起伏给出均匀表层厚度Z₁，表层介质电阻率为ρ，山脊模型内部电阻率与均匀大地电阻率一致，均为100 Ωm.以下分别给出表层电阻率为ρ=10 Ωm与ρ=1000 Ωm时，不同表层厚度情况下ATEM的地形效应.

图 17

Fig. 17

图 17 地形表层电阻率变化对比图 (a)山脊模型表层电阻率变化示意图;(b)不同表层厚度Bz响应对比图(ρ=10 Ωm);(c)视电阻率断面图对比图(ρ=10 Ωm); (d)不同表层厚度Bz响应对比图(ρ=1000Ωm);(e)视电阻率断面图对比图(ρ=1000 Ωm). Fig. 17 Comparison of terrain surface layer apparent resistivity change (a)Sketch map of the 3D trapezoidal hill model with surface layer apparent resistivity change;(b)Comparison of Bz-field with surface layer thicknes change(ρ=10 Ωm);(c)Comparison of apparent resistivity contours(ρ=10 Ωm);(d)Comparison of Bz-field with surface layer thicknes change(ρ=1000 Ωm);(e)Comparison of apparent resistivity contours(ρ=1000 Ωm).

图 17(b，c)为表层电阻率为ρ=10 Ωm时，表层厚度分别为Z₁=10 m、30 m、50 m与100 m时磁感应强度多测道及视电阻率断面对比图.图中结果显示，随着表层厚度的不断增加，ATEM响应也相应发生变化，但当表层厚度增加为30 m时，磁感应曲线值与厚度为50 m、100 m时基本一致，转换到视电阻率断面图中也一样，厚度30 m之后的视电阻率断面基本一致，由此可以说明，当表层为低阻时，地形效应主要与表层介质的电性有关，受表层介质影响大.

图 17(d，e)为表层电阻率为ρ=1000 Ωm时，层厚分别为Z₁=10 m、30 m、50 m与100 m时磁感应强度多测道及视电阻率断面对比图.与图(b，c)相比，由于表层是高阻介质，磁感应强度曲线变化更加和缓，不同厚度B_z曲线变化并不明显.当厚度为50 m时，磁感应强度曲线值与100时大致一致.由此表明，地形效应主要受到表层介质的电性影响，并且低阻表层影响比高阻表层影响更大.

在实际测量过程中，由于受到气流与天气等因素的影响，飞行高度并不是一个恒定值，由此将导致测量时的接收高度发生变化.因此，分析飞行高度对于地形效应的影响是很有必要的.以下同样以山脊模型为例，通过改变飞机的均匀离地高度，进一步讨论飞行高度对于ATEM地形效应的影响情况.

图 18(b，c)分别表示山脊地形情况下，均匀离地高度为h=0 m、20 m和50 m时，磁感应强度对比图及视电阻率断面对比图.由图 18b可知，不同飞行高度条件下的B_z响应曲线的形态是一致的，只是响应幅值存在差异；由于空气层的影响，飞行高度越高，B_z响应幅值越小.图 18c视电阻率断面图对比图中，不同的飞行高度产生的视电阻率断面结果也不尽相同，随着飞行高度的不断增大，山脊地形响应的影响强度也在不断增强，顶部的高阻异常区域与底部低阻异常区的范围与强度都在不断增大，地形效应随着飞行高度的增加而更加明显.

图 18

Fig. 18

图 18 地形飞行高度变化对比图 (a)山脊模型飞行高度变化示意图;(b)Bz响应对比图;(c)视电阻率断面图对比图. Fig. 18 Comparison of flight height change (a)Sketch map of the 3D trapezoidal hill model with flight height change;(b)Comparison of Bz-field with flight height change; (c)Comparison of apparent resistivity contours.

航空瞬变电磁的野外测量，一般采用沿地形飞行的飞行方式，但是这种飞行方式对于飞行器与飞行员自身的驾驶能力都具有一定的要求，当以上要求不满足时，则采用沿水平面飞行的飞行轨迹.因此，下文中将针对不同飞行轨迹测量条件下对于ATEM地形效应的影响情况进行分析.

图 19(b，c)分别表示飞行轨迹不同时，山脊地形所产生的视电阻率多测道曲线及视电阻率断面对比图.图中表明不同飞行轨迹条件下所产生的视电阻率曲线变化形态是基本一致的，但仍存在明显区别.当飞行轨迹为水平飞行时，在时间早期，由于地形的影响，地形周围的水平大地响应较真实电阻率更低，随着时间的推移，逐渐增大到真实电阻率值.由此说明，水平飞行情况下地形对于ATEM的影响范围更广，作用更强，沿地形起伏的飞行模式可以相对减小地形效应的影响程度.

图 19

Fig. 19

图 19 地形飞行轨迹变化对比图 (a)山脊飞行轨迹变化模型示意图;(b)视电阻率对比图;(c)视电阻率断面图对比图. Fig. 19 Comparison of flight path change (a)Sketch map of the 3D trapezoidal hill model with flight path change;(b)Comparison of apparent resistivity with different flight path; (c)Comparison of apparent resistivity contours.

本文基于时域有限差分算法实现了带地形条件下回线源航空瞬变电磁的三维正演，与前人计算结果进行对比，验证了算法的正确性.通过一系列三维地形的理论模型正演计算，对于带地形条件下ATEM的电磁响应有了初步认识，并得到了一些有价值的结论：

(1) 地形起伏对于ATEM数据有着显著的影响且主要集中在早期，随时间推移影响逐渐减弱；山脊地形在ATEM响应中为在山顶处产生高阻异常，而在山脚地形转折处为低阻异常，山谷则呈现完全相反的电阻率形态.

(2) ATEM的瞬变响应与地形尺寸参数的变化正相关，随着地形尺寸的增大，地形对ATEM数据的影响强度与影响范围也在逐渐的增大；ATEM的测量结果还主要受到测线方向的影响，沿测线方向的地形尺寸变化对于ATEM响应的影响更为显著.

(3) 地形的电性参数变化对于ATEM的瞬变响应具有显著影响.以山脊地形为例，当地形电阻率小于地下介质电阻率时，由于山脊电阻率与空气层及地表的电性差异增大，导致山脊底部转折处B_z响应进一步增强，从而导致视电阻率断面中低阻异常的范围与强度明显增强，顶部高阻异常受到压制，整体地形模型呈现出低阻异常影响.反之，当山脊电阻率大于地下介质电阻率时，底部转折处的低阻异常范围与强度减小，顶部高阻异常进一步增强，整体模型主要体现为顶部的高阻异常.并且地形效应主要受到表层介质的电性影响，低阻表层影响比高阻表层影响更大.

(4) ATEM的地形效应还受到飞行高度与飞行轨迹因素的影响.不同飞行高度所产生的ATEM响应结果也不尽相同，随着飞行高度的不断增大，地形响应的影响强度与影响范围都在不断增大.采用不同飞行轨迹进行测量所得到的ATEM响应的结果也存在区别，沿地形飞行情况更能够比水平飞行更好的减小地形效应的影响.

综上所述，地形的起伏变化对与ATEM数据影响显著，并且不同的地形类型对于ATEM数据的影响不同.通过文中多个模型算例结果可知，地形对于地下异常体的影响不容忽视，特别是将其按照平坦地表情况进行解释处理时，会导致计算结果偏离真实模型，造成解释的误差.本文通过大量模型进一步分析多种因素对于ATEM地形效应的影响情况，并总结规律，以上结果对于地形效应的识别与后期处理解释提供重要的理论依据，为下一步对地形效应进行校正、提高航空电磁的数据解释精度提供有效的理论基础.

作者向孙怀凤老师及卢绪山博士对于本文中的三维FDTD算法所提供的帮助表示感谢，特别对审稿专家对于本文所提出的宝贵意见表示衷心的感谢，最后感谢编辑老师的辛苦工作.