地球物理学报  2021, Vol. 64 Issue (1): 327-342   PDF    
引用本文
张铭, Farquharson Colin G., 刘长胜. 2021. 2.5维起伏地表条件下频率域地空电磁正演模拟. 地球物理学报, 64(1): 327-342, doi: 10.6038/cjg2020N0182.  
Zhang M, Colin G. Farquharson, Liu C S. 2021. 2.5-D forward modeling of the frequency-domain ground-airborne electromagnetic response in areas with topographic relief. Chinese J. Geophys. (in Chinese), 64(1): 327-342, doi: 10.6038/cjg2020N0182.  
2.5维起伏地表条件下频率域地空电磁正演模拟
张铭1,2,3, Farquharson Colin G.3, 刘长胜1,2     
1. 吉林大学仪器科学与电气工程学院, 长春 130061;
2. 吉林大学地球信息探测仪器教育部重点实验室, 长春 130061;
3. Department of Earth Sciences, Memorial University of Newfoundland, St. John's, Canada A1B3S7
收稿日期 2020-02-14, 2020-07-28 收修定稿
基金项目: 国家自然科学基金面上项目（42074152）以及国家自然科学基金国家重大科研仪器研制项目（41827803），国家留学基金委（201706170173）联合资助
第一作者简介: 张铭, 女, 1989年生, 博士研究生, 主要从事电磁法正反演方法研究.E-mail:zhangm_sunflower@163.com
通讯作者: 刘长胜, 男, 1979年生, 教授, 主要从事电磁探测技术及仪器的研究.E-mail:liuchangsheng@jlu.edu.cn
摘要：频率域地空电磁探测方法是指在地面布设人工场源，在空中测量电磁场的一种高效的地球物理勘探技术.该方法具有大范围、高通过性、快速测量的优势，尤其适合崎岖山地、沙漠、沼泽、海陆交互带等复杂地貌区域的资源勘查.但是这些地区的地形起伏通常较大，因此分析地形对地空电磁响应的影响具有重要意义.本文利用有限元法对频率域地空电磁响应进行了正演计算，分析了起伏地表条件下的频率域地空电磁响应特征.首先利用傅里叶变换将2.5维问题转化成二维问题，利用伽辽金加权余量法推导了相应的离散有限元方程组.采用任意四边形单元对区域进行不均匀网格剖分，源和异常体附近网格加密处理，保证计算精度，远离目标区域网格逐渐稀疏，模拟无穷远边界，降低对计算资源的要求.在单元内进行插值，将有限元方程组变换为线性方程组，采用总场算法，利用具有一定面积的伪δ函数表达源电流分布，源项近似为分布在以电偶极源为中心的25个节点上.通过求解线性方程组得到波数域电磁响应，再对波数域电磁场响应进行反傅里叶变换从而获得空间域2.5维频率域电磁场值.通过对比2.5维正演结果与均匀半空间解析解，验证了本文算法的精度，同时本文还对地空电磁场与地面电磁场的响应特性进行了对比.
关键词: 地空电磁法      2.5维      正演      起伏地表      总场算法     
2.5-D forward modeling of the frequency-domain ground-airborne electromagnetic response in areas with topographic relief
ZHANG Ming1,2,3, Farquharson Colin G.3, LIU ChangSheng1,2     
1. College of Instrumentation and Electrical Engineering, Jilin University, Changchun 130061, China;
2. Key Lab of Geo-Exploration Instrumentation, Ministry of Education, Changchun 130061, China;
3. Department of Earth Sciences, Memorial University of Newfoundland, St. John's, A1 B3 S7, Canada
Abstract: Ground-airborne frequency-domain electromagnetic method is an effective and efficient geophysical prospecting method. It usually places the controlled source on the ground and measure the electromagnetic field in the air. This method can achieve fast measurement in large scale, is especially suitable for resource exploration in rugged mountains, deserts, swamps, sea-land interaction zones and other areas with complex topography relief. However, in these areas, the topography relief is serious, so the study of the topographic effect on ground-airborne electromagnetic response is meaningful. In this paper, the 2.5D frequency domain ground-airborne electromagnetic response of topography earth model was calculated based on the finite element method and the response characteristic is analyzed. Fourier transform was applied to the Maxwell's equations to transform a 2.5D problem into a 2D problem, and then the discrete finite element equation was deduced using the Galerkin weighted residual method. The surveyed region is subdivided into many arbitrary quadrilateral elements with different sizes, the ones region near the source and the target have smaller size to ensure the accuracy of numerical calculation, the ones away from the target region have bigger size to simulate infinite boundary and reduce the requirements for computing resource. Interpolation is conducted in the element to transform the finite element equations to the linear system of equations. Using the total field algorithm, the pseudo-delta function with a certain area is employed to approximate the distribution of the source current, the source terms are approximately distributed at 25 nodes centered on the electric dipole source. By solving the equation, the solution of electromagnetic field in wavenumber domain was obtained, then the solution of electromagnetic field in spatial domain can be calculated by the Invers Fourier Transformation. The accuracy of our algorithm was verified by the comparison of our 2.5D forward modeling results and the analytical results. We also compared the electromagnetic field response characteristics measured in the air and on the ground.
Keywords: Ground-airborne electromagnetic method    2.5D    Forward modeling    Topographic surface    Total field algorithm    
0 引言

地空电磁法最早起源于俄罗斯和西欧地区，随后在加拿大得到发展和应用.早在20世纪70年代，就诞生了Turair地空电磁探测系统并在安大略湖北部的金属硫化物矿体探测方面取得了成功应用(Pemberton et al., 1970).随后，在20世纪90年代，包括澳大利亚的FLAIRTEM系统、加拿大的TerraAir系统和日本的GREATEM系统在内的几种地空电磁探测系统也相继问世(Elliott, 1998; Mogi et al., 1998; Smith et al., 2001).近些年，国内地空电磁法也得到了较好的发展，林君、李貅、王绪本等教授带领的团队分别在地空电磁法的方法理论、仪器系统研制以及野外应用等方面进行了深入研究(张莹莹，2013；李肃义等，2013；许洋，2014；康利利，2019).地空电磁法，主要包括时间域地空电磁法和频率域地空电磁法(Thomson et al., 2007; Abdallah et al., 2017).时间域地空电磁法主要关注发射电流关断后的二次感应信号，信号较弱，为了保证测量信号的强度，通常在距离发射源较近的区域进行测量，在浅地表探测中可以发挥良好作用(Mogi et al., 2009; 嵇艳鞠等, 2013; 李貅等，2015; Wu et al., 2019).频率域地空电磁法通过地面布设发射源激励空间电磁场，空中接收不同发射频率下的电磁场响应获取不同深度下的大地电阻率信息，在整个采集过程中，电流是连续发射的，信号较强，可测量范围较广，有利于实现大范围、大深度的探测，从而达到地质填图、油气勘探、矿产普查及环境和工程调查等主要应用目的(Zhou et al., 2016；Lin et al., 2019).

地空电磁法中的发射源可以是电性源(如接地长导线源)或者是磁性源(如大回线源).长导线源由于其布设简单，实施方便，较为常用.在远区测量时，收发距远大于发射天线长度，长导线源可看作电偶极子源.由于地质结构是二维的，测线通常沿着地质构造倾向方向布设.地空电磁法需要测量Hz分量，当电偶极源沿着地质倾向方向布设，沿地质构造倾向方向的测线位于电偶极源的延长线上时，Hz分量幅度微弱.当电偶极源沿着地质走向方向布设，沿地质构造倾向方向的测线通过电偶极源的中点时，Hz分量的幅度较强，有利于实际测量(底青云和王若，2008).本文后续分析中，装置形式采用后者.

地空电磁法结合了地面大功率发射和空中快速非接触式接收的优势，为地表条件比较复杂、地面难以进入区域的地下电性信息快速获取提供了良好的解决方案，但是这些区域地表结构通常较为复杂，地形对电磁响应的影响不可忽略(殷长春等，2015).近年来，很多学者关注到了地形对电磁响应的影响，在地面电磁法和航空电磁法中，地形影响的分析已经比较深入，对地空电磁地形影响的分析还有待更深入的研究(Sasaki and Nakazato, 2003；底青云等，2004；Nam et al., 2007；张继峰等，2013).

多种数值模拟方法可以应用于地空电磁地形响应特征的研究，如边界元法、积分方程法、有限差分法、有限元法等.有限元方法比较适用于复杂电性结构的模拟，而且具有较高的精度，本文采用该方法(徐世浙，1994；Newman and Alumbaugh, 1995；金建铭，1998；Sasaki and Nakazato, 2003).利用有限元方法进行电磁场模拟有两种常用算法，一种是二次场算法，一种是总场算法，二次场算法无需对源附近网格进行加密剖分，可减少网格数量，节约计算时间，同时可以消除场源奇异性，但是由于其需要额外计算一次场，而复杂模型的一次场又较难模拟，因此在复杂电性结构模型的处理上有一定局限(Lu et al., 1999; Li and Key, 2007; Mitsuhata, 2000).本文采用总场算法对带地形的地空电磁响应进行了数值计算和响应特性分析.通过与一维解析解的对比验证了该方法的正确性.同时本文还将带地形的地空电磁响应与带地形的地面电磁响应进行了对比，从而评价地空电磁法的异常分辨能力和受地形影响程度.

1 2.5维正演算法及结果验证 1.1 正演算法

二维地电模型示意图如图 1所示，模型由空气、大地和异常体组成.y方向与地质构造走向方向平行，电导率、介电常数、磁导率等电性参数在y方向上不变，仅在x、z方向变化.电偶极源沿y方向，位于坐标原点.测线沿x方向布放.

图 1 地电模型与坐标示意图 Fig. 1 Schematic diagram of geo-electric model and coordinate

利用总场算法对电磁场进行求解，电磁场总场满足如下形式的麦克斯韦方程：

(1)

(2)

式中，w表示角频率，μ0表示真空磁导率，σ表示电导率，E表示电场，H表示磁场，Jc表示源电流密度.

对式(1)、(2)中E、H的每个分量及源项，沿构造走向y方向做傅里叶变换，从空间域转换到波数域，即可将三维的微分方程转换为二维的情况，从而得到波数域中, 两个相互关联的电磁场分量的耦合方程，式(3)和式(4).

(3)

(4)

式中, 和表示波数域中的电磁场分量，ky表示波数值.

求解上述方程组即可得到波数域中的电磁场分量.为了求解方程组，首先需要得到离散化的有限元方程.将整个计算区域剖分为多个小四边形单元，利用伽辽金加权余量法即可得到电偶极源沿地质构造走向方向布放时的离散化有限元方程：

(5)

(6)

求解有限元方程首先需要对研究区域进行网格剖分，本文采用任意四边形剖分，对任意四边形进行等参变换，将图 2a所示的x-z坐标系下的不规则四边形子单元转化成图 2a所示的ξ-η坐标系下的正方形母单元.

图 2 单元编号和坐标变换关系示意图 (a)子单元; (b)母单元. Fig. 2 Schematic diagram of cell numbering and coordinate conversion (a) Sub-element; (b) Parent element.

母单元的形函数由(7)式表示，子单元和母单元之间的坐标变换关系由(8)式表示.网格剖分时，采用不均匀网格剖分，对源和异常体附近网格进行加密，保证计算精度；远离目标区域网格逐渐稀疏，模拟无穷远边界，并降低对计算资源的要求.

(7)

(8)

(9)

对和进行插值，在每个单元内, 和满足(9)式.所以和可由(10)式表示.将单元离散变量(10)式代人(5)式和(6)式，可得到每个剖分小单元的矩阵形式的离散有限元方程，由(11)式表示.

(10)

(11)

本文利用具有一定面积的伪δ函数表达(11)式右端的源电流分布，源项近似为分布在以电偶极源为中心的25个节点上.源沿着地质构造走向y方向布放时，波数域的电流密度可用(12)式表示，伪δ函数可用(13)、(14)式表示.将(12)式、(13)式、(14)式代入(5)式的右端，即可得到每个小单元的源项.

(12)

(13)

(14)

将整个剖分单元的系数矩阵及源项矩阵按节点的网格排列顺序进行扩展，得到大型线性方程组(15).求解该方程组即可得到各个节点上的波数域电磁场值.再对波数域电磁场值做傅里叶逆变换即可得到空间域电磁场值.

(15)

1.2 结果验证

为了验证上述算法及程序计算的正确性，将本文的2.5维正演计算结果与均匀半空间的解析解进行了对比.100 Ωm均匀大地，地面y=0 m和y=500 m测线上的电磁场分量及卡尼亚视电阻率在100 Hz和10 Hz两个频率下的结果如图 3所示.本文给出的结果均为单位电矩激励下的电磁场值.由于本文实际关注的区域位于距离发射源较远的地方，平面波理论在该区适用，因此可利用卡尼亚视电阻率(简称视电阻率)计算公式计算关注区域的视电阻率值(Cargniard, 1953).图中实线为均匀半空间解析解，圆点为本文计算结果.从图中可以看出，y=0 m测线和y=500 m测线上，较高频率(100 Hz)和较低频率(10 Hz)下，Ey, Hx, Hz三个分量的计算结果均与解析解对应良好，视电阻率也能够准确对应模型的实际电阻率，验证了本文计算的正确性.

图 3 2.5维有限元计算结果与解析解对比图 (a) Ey幅值(y=0 m, 10 Hz)；(b) Hx幅值(y=0 m, 10 Hz)；(c) Hz幅值(y=0 m, 10 Hz)；(d)视电阻率(y=0 m, 10 Hz); (e) Ey幅值(y=0 m, 100 Hz)；(f) Hx分量幅值(y=0 m, 100 Hz)；(g) Hz幅值(y=0 m, 100 Hz)；(h)视电阻率(y=0 m, 100 Hz); (i) Ey幅值(y=500 m, 10 Hz)；(j) Hx幅值(y=500 m, 10 Hz)；(k) Hz分量幅值(y=500 m, 10 Hz)；(l)视电阻率(y=500 m, 10 Hz); (m) Ey幅值(y=500 m, 100 Hz)；(n) Hx幅值(y=500 m, 100 Hz)；(o) Hz幅值(y=500 m, 100 Hz)；(p)视电阻率(y=500 m, 100 Hz). Fig. 3 Comparisons between analytical results and 2.5D FEM (a) Amplitude of Ey (y=0 m, 10 Hz); (b) Amplitude of Hx (y=0 m, 10 Hz); (c) Amplitude of Hz (y=0 m, 10 Hz); (d) Cagniard resistivity (y=0 m, 10 Hz); (e) Amplitude of Ey (y=0 m, 100 Hz); (f) Amplitude of Hx (y=0 m, 100 Hz); (g) Amplitude of Hz (y=0 m, 100 Hz); (h) Cagniard resistivity (y=0 m, 100 Hz); (i) Amplitude of Ey (y=500 m, 10 Hz); (j) Amplitude of Hx (y=500 m, 10 Hz); (k) Amplitude of Hz (y=500 m, 10 Hz); (l) Cagniard resistivity (y=500 m, 10 Hz); (m) Amplitude of Ey (y=500 m, 100 Hz); (n) Amplitude of Hx (y=500m, 100 Hz); (o) Amplitude of Hz (y=500 m, 100 Hz); (p) Cagniard resistivity (y=500 m, 100 Hz).
2 无地形异常体模型的电磁场响应特性

为了分析2.5维无地形异常体模型的响应特性，并对比地面电磁法和地空电磁法对异常体的响应能力，对图 4a所示的无地形低阻体模型(简称低阻体模型)和图 4b所示的无地形高阻体模型(简称高阻体模型)在地面和空中的电磁场响应进行了计算，并对比了地面测量和空中测量时电磁场幅度的相对异常，异常体模型参数如表 1所示.地面测量时，计算了测线y=0 m的电磁场及视电阻率响应，空中测量时，计算了两个不同高度(z=50 m和z=100 m)下y=0 m测线的磁场响应，测量频率为100 Hz.无地形异常体模型电磁场响应曲线如图 5所示，每个子图中的小尺寸细节图为异常体上方的局部异常特征.

图 4 无地形异常体模型及网格剖分正视图 (a)低阻体模型正视图; (b)高阻体模型正视图. Fig. 4 Front view of target earth model without topography (a) Conductive target earth model; (b) Resistive target earth model.
表 1 异常体模型参数 Table 1 Parameters of target earth model
图 5 无地形异常体模型电磁场响应曲线 (a)地面Ey幅值(低阻体)；(b)地面Hx幅值(低阻体)；(c)地面Hz幅值(低阻体)；(d)地面视电阻率(低阻体)；(e)地面Ey幅值(高阻体)；(f)地面Hx幅值(高阻体)；(g)地面Hz幅值(高阻体)；(h)地面视电阻率(高阻体)；(i)空中50 m Hx幅值(低阻体)；(j)空中50 m Hz幅值(低阻体)；(k)空中50 m Hx幅值(高阻体)；(l)空中50 m Hz幅值(高阻体)；(m)空中100 m Hx幅值(低阻体)；(n)空中100 m Hz幅值(低阻体)；(o)空中100 m Hx幅值(高阻体)；(p)空中100 m Hz幅值(高阻体)；(q) Hx相对异常(低阻体)；(r) Hz相对异常(低阻体)；(s) Hx相对异常(高阻体)；(t) Hz相对异常(高阻体). Fig. 5 Electromagnetic response curves of target earth model without topography (a) Amplitude of ground measured Ey (conductive target earth model); (b) Amplitude of ground measured Hx (conductive target earth model); (c) Amplitude of ground measured Hz (conductive target earth model); (d) Cagniard resistivity (conductive target earth model); (e) Amplitude of ground measured Ey (resistive target earth model); (f) Amplitude of ground measured Hx (resistive target earth model); (g) Amplitude of ground measured Hz (resistive target earth model); (h) Cagniard resistivity (resistive target earth model); (i) Amplitude of 50 m-air-measured Hx (conductive target earth model); (j) Amplitude of 50 m-air-measured Hz (conductive target earth model); (k) Amplitude of 50 m-air-measured Hx (resistive target earth model); (l) Amplitude of 50 m-air-measured Hz (resistive target earth model); (m) Amplitude of 100 m-air-measured Hx (conductive target earth model); (n) Amplitude of 100 m-air-measured Hz (conductive target earth model); (o) Amplitude of 100 m-air-measured Hx (resistive target earth model); (p) Amplitude of 100 m-air-measured Hz (resistive target earth model); (q) Relative anomaly of Hx (conductive target earth model); (r) Relative anomaly of Hz (conductive target earth model); (s) Relative anomaly of Hx (resistive target earth model); (t) Relative anomaly of Hz (resistive target earth model).

由图 5可以看出，在地面测量时，无地形异常体模型的Hz分量在异常体附近产生了最明显的异常响应，Ey分量次之，Hx分量最不明显.在异常体边界附近Hz的异常达到峰值，在异常体中心附近Ey分量和Hx分量的异常达到峰值.在x=7000~9000 m范围内，低阻体模型对应的视电阻率整体呈现低阻响应，值约为80 Ωm，高阻体模型对应的视电阻率整体呈现高阻响应，值约为130 Ωm，其他位置上视电阻率值为100 Ωm，视电阻率变化趋势与图 4模型的电阻率的变化趋势一致.在空中测量时，无法获取电场分量，所以只需关注磁场分量.在空中不同高度测量时，其结果均与地面测量结果类似，Hx分量的异常并不明显，对异常体的响应不灵敏，Hz分量在低阻体和高阻体上方均出现了明显的异常，响应形态与地面的Hz分量相同，但响应幅度略有不同.

为了进一步定量分析空中测量和地面测量时，磁场分量对于异常体的响应灵敏程度差异，以评价空中测量的可行性，本文计算了三种测量方式下的磁场分量相对异常，结果如图 5q—t所示.从图中可以看出，对于Hx分量，在低阻体上方，地面和空中测量的Hx分量均呈现正异常，在高阻体上方，地面和空中测量的Hx分量均呈现负异常，异常峰值位置位于异常体中心位置.而对于Hz分量，在低阻体上方，地面和空中测量的Hz分量均先呈现正异常后呈现负异常，在高阻体上方，呈现相反的趋势，异常峰值位置约位于异常体边界位置.同时可以看出，虽然Hz分量的总场幅度小于Hx分量，但是其相对异常却远远大于Hx分量.而且地面测量的相对异常与空中50 m测得的相对异常几乎相同, 空中100 m测得的相对异常仅有小幅度减小，这说明空中测量仍然可以有效识别异常体，异常分辨率并没有显著降低.

3 带地形无异常体模型电磁场响应特性

为了分析带地形无异常体模型的响应特性，并对比地面电磁法和地空电磁法受地形影响的程度，对图 6a所示的凹地形无异常体模型简称凹地形模型和图 6b所示的凸地形无异常体模型简称凸地形模型在地面和空中的电磁场响应进行了计算.带地形无异常体模型参数如表 2所示.测量频率为100 Hz.地面测量时，测线随地形起伏.空中测量时，可以随地形测量也可以在同一高度进行测量.本文对以下三种测量方式下的电磁场响应进行了计算和分析：

图 6 带地形无异常体模型正视图 (a)凸地形模型正视图; (b)凹地形模型正视图. Fig. 6 Front view of topography earth model (a) Convex earth model; (b) Concave earth model.
表 2 带地形无异常体模型参数 Table 2 Parameters of topography earth model without target

第一种：地面测量时，测线沿地表随地形起伏，测量示意图如图 7a所示.

图 7 测线布置方式示意图 (a)地面测量; (b)空中随地形测量; (c)空中同一高度测量. Fig. 7 Schematic diagram of survey line configuration (a) Measured on ground surface along the topography; (b) Measured in the air along the topography; (c) Measured in the air at constant height.

第二种：空中随地形测量时，测线随地形起伏，位于地表上方50 m，测量示意图如图 7b所示.

第三种：空中同一高度测量时，测线不随地形起伏，固定位于同一海拔高度100 m，测量示意图如图 7c所示.

图 7中灰色虚线表示测线，黑色实线表示起伏地表.

图 7a、图 7b、图 7c所对应的三种测量方式下的电磁场响应曲线分别如图 8、图 9、图 10所示，每个子图中的小尺寸细节图为地形上方的局部异常特征.

图 8 带地形无异常体模型，地面测量时的电磁场响应曲线 (a) Ey幅值(凸地形)；(b) Hx幅值(凸地形)；(c) Hz幅值(凸地形)；(d)视电阻率(凸地形)；(e) Ey相对异常(凸地形)；(f) Hx相对异常(凸地形)；(g) Hz相对异常(凸地形)；(h)视电阻率相对异常(凸地形)；(i) Ey幅值(凹地形)；(j) Hx幅值(凹地形)；(k) Hz幅值(凹地形)；(l)视电阻率(凹地形)；(m) Ey相对异常(凹地形)；(n) Hx相对异常(凹地形)；(o) Hz相对异常(凹地形)；(p)视电阻率相对异常(凹地形). Fig. 8 Electromagnetic response curves of topography earth model without target when measured on ground surface (a) Amplitude of Ey (convex earth model); (b) Amplitude of Hx (convex earth model); (c) Amplitude of Hz (convex earth model); (d) Cagniard resistivity (convex earth model); (e) Relative anomaly of Ey (convex earth model); (f) Relative anomaly of Hx (convex earth model); (g) Relative anomaly of Hz (convex earth model); (h) Relative anomaly of Cagniard resistivity (convex earth model); (i) Amplitude of Ey (concave earth model); (j) Amplitude of Hx (concave earth model); (k) Amplitude of Hz (concave earth model); (l) Cagniard resistivity (concave earth model); (m) Relative anomaly of Ey (concave earth model); (n) Relative anomaly of Hx (concave earth model); (o) Relative anomaly of Hz (concave earth model); (p) Relative anomaly of Cagniard resistivity (concave earth model).
图 9 带地形无异常体模型，空中随地形测量时的电磁场响应曲线 (a) Hx幅值(凸地形)；(b) Hz幅值(凸地形)；(c) Hx幅值(凹地形)；(d) Hz幅值(凹地形)；(e) Hx相对异常(凸地形)；(f) Hz相对异常(凸地形)；(g) Hx相对异常(凹地形)；(h) Hz相对异常(凹地形). Fig. 9 Electromagnetic response curves of topography earth model without target when measured in the air along topography (a) Amplitude of Hx (convex earth model); (b) Amplitude of Hz (convex earth model); (c) Amplitude of Hx (concave earth model); (d) Amplitude of Hz (concave earth model); (e) Relative anomaly of Hx (convex earth model); (f) Relative anomaly of Hz (convex earth model); (g) Relative anomaly of Hx (concave earth model); (h) Relative anomaly of Hz (concave earth model).
图 10 带地形无异常体模型，空中同一高度测量时的电磁场响应曲线 (a) Hx幅值(凸地形)；(b) Hz幅值(凸地形)；(c) Hx幅值(凹地形)；(d) Hz幅值(凹地形)；(e) Hx相对异常(凸地形)；(f) Hz相对异常(凸地形)；(g) Hx相对异常(凹地形)；(h) Hz相对异常(凹地形). Fig. 10 Electromagnetic response curves of topography earth model without target when measured in the air at constant height (a) Amplitude of Hx (convex earth model); (b) Amplitude of Hz (convex earth model); (c) Amplitude of Hx (concave earth model); (d) Amplitude of Hz (concave earth model); (e) Relative anomaly of Hx (convex earth model); (f) Relative anomaly of Hz (convex earth model); (g) Relative anomaly of Hx (concave earth model); (h) Relative anomaly of Hz (concave earth model).

通过图 8d和图 8l可以看出，地面测量时，带地形无异常体模型简称带地形模型和无地形无异常体模型简称无地形模型的视电阻率响应曲线在平坦区域几乎重合，但是在地形上方，即x=7600~8600 m范围内，这两个曲线出现了分离，凸地形上方，带地形模型的视电阻率呈现小幅度升高趋势，值约为103 Ωm，凹地形上方，带地形模型的视电阻率呈现小幅度下降趋势，值约为96 Ωm，其他位置上视电阻率值为100 Ωm.通过图 8a—c和图 8i—k中的小尺寸细节图也可以看出，带地形模型和无地形模型的电磁场响应曲线在地形上方也出现了分离.结合图 8e—g和图 8m—o所示的电磁场相对异常曲线可以看出，Ey和Hx分量在地形上方的异常幅度比较平稳，相对异常幅度较小，凸地形上方，Ey和Hx分量的相对异常呈现正异常，凹地形上方，其相对异常呈现负异常，地形以外区域，其相对异常趋近于0，Hz分量在地形上方存在异常且异常变化较快，其相对异常极值出现在地形边界处，幅度较大，地形以外区域，其相对异常趋近于0.

通过对比图 8b、9a、10a，对比图 8j、9c、10c可以看出，空中随地形测量及空中同一高度测量时，无论凸地形模型与凹地形模型，其Hx响应曲线在地形上方的形态与地面测量时相同，同时对比图 8f、9e、10e，对比图 8n、9g、10g可以看出，三种测量方式下，其相对异常曲线的幅度也比较接近.同理，通过对比图 8c、9b、10b，对比8k、9d、10d可以看出，空中随地形测量及空中同一高度测量时，无论凸地形模型与凹地形模型，其Hz响应曲线在地形上方的形态与地面测量时相同，同时对比图 8g、9f、10f，图 8o、9h、10h可以看出，三种测量方式下，其对异常响应曲线的幅度同样比较接近.可见，空中随地形测量及空中同一高度测量时，地形对磁场响应曲线所产生的影响与地面测量时相似，并没有本质上的明显差异.

4 带地形异常体模型电磁场响应特性

为了分析带地形异常体模型的地面和空中电磁场响应特性，对图 11所示的带地形异常体模型在地面和空中的电磁场响应进行了计算，模型中的异常体参数设置同表 1，地形参数设置同表 2，测线布置方式如图 7.测量频率为100 Hz.地面随地形测量时，凸地形异常体模型的电磁场响应曲线和凹地形异常体模型的电磁场响应曲线分别如图 12、图 13所示，空中随地形测量时和空中同一高度测量时，带地形异常体模型的磁场响应曲线分别如图 14、图 15所示，每个子图中的小尺寸细节图为异常体上方的局部异常特征.

图 11 带地形异常体模型正视图 (a)凹地形低阻体模型正视图; (b)凸地形低阻体模型正视图；(c)凹地形高阻体模型正视图; (d)凸地形高阻体模型正视图. Fig. 11 Front view of topography earth model with target (a) Front view of concave earth model with conductive target; (b) Front view of convex earth model with conductive target; (c) Front view of concave earth model with resistive target; (d) Front view of convex earth model with resistive target.
图 12 地面测量时，凸地形异常体模型的电磁场响应曲线 (a) Ey幅值(低阻体)；(b) Hx幅值(低阻体)；(c) Hz幅值(低阻体)；(d)视电阻率(低阻体)；(e) Ey相对异常(低阻体)；(f) Hx相对异常(低阻体)；(g)Hz相对异常(低阻体)；(h)视电阻率相对异常(低阻体)；(i) Ey幅值(高阻体)；(j) Hx幅值(高阻体)；(k) Hz幅值(高阻体)；(l)视电阻率(高阻体)；(m) Ey相对异常(高阻体)；(n) Hx相对异常(高阻体)；(o) Hz相对异常(高阻体)；(p)视电阻率相对异常(高阻体). Fig. 12 Electromagnetic response curves of convex topography target earth model when measured on ground surface (a) Amplitude of Ey (conductive target); (b) Amplitude of Hx (conductive target); (c) Amplitude of Hz (conductive target); (d) Cagniard resistivity (conductive target); (e) Relative anomaly of Ey (conductive target); (f) Relative anomaly of Hx (conductive target); (g) Relative anomaly of Hz (conductive target); (h) Relative anomaly of Cagniard resistivity (conductive target); (i) Amplitude of Ey (resistive target); (j) Amplitude of Hx (resistive target); (k) Amplitude of Hz (resistive target); (l) Cagniard resistivity (resistive target); (m) Relative anomaly of Ey (resistive target); (n) Relative anomaly of Hx (resistive target); (o) Relative anomaly of Hz(resistive target); (p) Relative anomaly of Cagniard resistivity (resistive target).
图 13 地面测量时，凹地形异常体模型的电磁场响应曲线 (a) Ey幅值(低阻体)；(b) Hx幅值(低阻体)；(c) Hz幅值(低阻体)；(d)视电阻率(低阻体)；(e) Ey相对异常(低阻体)；(f) Hx相对异常(低阻体)；(g)Hz相对异常(低阻体)；(h)视电阻率相对异常(低阻体)；(i) Ey幅值(高阻体)；(j) Hx幅值(高阻体)；(k) Hz幅值(高阻体)；(l)视电阻率(高阻体)；(m) Ey相对异常(高阻体)；(n) Hx相对异常(高阻体)；(o) Hz相对异常(高阻体)；(p)视电阻率相对异常(高阻体). Fig. 13 Electromagnetic response curves of concave topography target earth model when measured on ground surface (a) Amplitude of Ey (conductive target); (b) Amplitude of Hx (conductive target); (c) Amplitude of Hz (conductive target); (d) Cagniard resistivity (conductive target); (e) Relative anomaly of Ey (conductive target); (f) Relative anomaly of Hx (conductive target); (g) Relative anomaly of Hz (conductive target); (h) Relative anomaly of Cagniard resistivity (conductive target); (i) Amplitude of Ey (resistive target); (j) Amplitude of Hx (resistive target); (k) Amplitude of Hz (resistive target); (l) Cagniard resistivity (resistive target); (m) Relative anomaly of Ey (resistive target); (n) Relative anomaly of Hx (resistive target); (o) Relative anomaly of Hz(resistive target); (p) Relative anomaly of Cagniard resistivity (resistive target).
图 14 空中随地形测量时，带地形异常体模型的磁场响应曲线 (a) Hx幅值(凸地形低阻体); (b) Hz幅值(凸地形低阻体); (c) Hx幅值(凸地形高阻体); (d) Hz幅值(凸地形高阻体); (e) Hx相对异常(凸地形低阻体); (f) Hz相对异常(凸地形低阻体); (g) Hx相对异常(凸地形高阻体); (h) Hz相对异常(凸地形高阻体); (i) Hx幅值(凹地形低阻体); (j) Hz幅值(凹地形低阻体); (k) Hx幅值(凹地形高阻体); (l) Hz幅值(凹地形高阻体); (m) Hx相对异常(凹地形低阻体); (n) Hz相对异常(凹地形低阻体); (o) Hx相对异常(凹地形高阻体); (p) Hz相对异常(凹地形高阻体). Fig. 14 Magnetic field curves of topography target earth model when measured in the air along the topography (a) Amplitude of Hx (convex earth model with conductive target); (b) Amplitude of Hz (convex earth model with conductive target); (c) Amplitude of Hx (convex earth model with resistive target); (d) Amplitude of Hz (convex earth model with resistive target); (e) Relative anomaly of Hx (convex earth model with conductive target); (f) Relative anomaly of Hz (convex earth model with conductive target); (g) Relative anomaly of Hx (convex earth model with resistive target); (h) Relative anomaly of Hz (convex earth model with resistive target); (i) Amplitude of Hx (concave earth model with conductive target); (j) Amplitude of Hz (concave earth model with conductive target); (k) Amplitude of Hx (concave earth model with resistive target); (l) Amplitude of Hz (concave earth model with resistive target); (m) Relative anomaly of Hx (concave earth model with conductive target); (n) Relative anomaly of Hz (concave earth model with conductive target); (o) Relative anomaly of Hx (concave earth model with resistive target); (p) Relative anomaly of Hz (concave earth model with resistive target).
图 15 空中同一高度测量时，带地形异常体模型的磁场响应曲线 (a) Hx幅值(凸地形低阻体); (b) Hz幅值(凸地形低阻体); (c) Hx幅值(凸地形高阻体); (d) Hz幅值(凸地形高阻体); (e) Hx相对异常(凸地形低阻体); (f) Hz相对异常(凸地形低阻体); (g) Hx相对异常(凸地形高阻体); (h) Hz相对异常(凸地形高阻体); (i) Hx幅值(凹地形低阻体); (j) Hz幅值(凹地形低阻体); (k) Hx幅值(凹地形高阻体); (l) Hz幅值(凹地形高阻体); (m) Hx相对异常(凹地形低阻体); (n) Hz相对异常(凹地形低阻体); (o) Hx相对异常(凹地形高阻体); (p) Hz相对异常(凹地形高阻体). Fig. 15 Magnetic field curves of topography target earth model when measured in the air at constant height (a) Amplitude of Hx (convex earth model with conductive target); (b) Amplitude of Hz (convex earth model with conductive target); (c) Amplitude of Hx (convex earth model with resistive target); (d) Amplitude of Hz (convex earth model with resistive target); (e) Relative anomaly of Hx (convex earth model with conductive target); (f) Relative anomaly of Hz (convex earth model with conductive target); (g) Relative anomaly of Hx (convex earth model with resistive target); (h) Relative anomaly of Hz (convex earth model with resistive target); (i) Amplitude of Hx (concave earth model with conductive target); (j) Amplitude of Hz (concave earth model with conductive target); (k) Amplitude of Hx (concave earth model with resistive target); (l) Amplitude of Hz (concave earth model with resistive target); (m) Relative anomaly of Hx (concave earth model with conductive target); (n) Relative anomaly of Hz (concave earth model with conductive target); (o) Relative anomaly of Hx (concave earth model with resistive target); (p) Relative anomaly of Hz (concave earth model with resistive target).

通过图 12d、图 12l及13d、图 13l可以看出，地面测量时，带地形异常体模型和带地形无异常体模型的视电阻率曲线在平坦区域几乎重合，在地形边界位置附近形态幅度也较相近，但是在异常体上方，这两个曲线出现了明显的分离，低阻体上方，带地形异常体模型的视电阻率值明显下降；高阻体上方，带地形异常体模型的视电阻率值明显升高.通过图 12a—c、图 12i—k及13a—c、图 13i—k中的小尺寸细节图也可以看出，带地形异常体模型和带地形无异常体模型的电磁场响应曲线在异常体上方出现了明显分离，结合图 12e—g、图 12m—o及图 13e—g、图 13m—o所示的电磁场相对异常曲线可以看出，在低阻体上方，Ey呈现负异常，Hx呈现正异常，Hz先呈现正异常后呈现负异常，在高阻体上方，Ey呈现正异常，Hx呈现负异常，Hz先呈现负异常后呈现正异常.

同理，通过图 14和图 15可以看出，当空中随地形测量及空中同一高度测量时，带地形异常体模型和带地形无异常体模型的磁场响应曲线与地面随地形测量时特点相似，均为平坦区域几乎重合，地形边界位置附近形态幅度相近，异常体上方曲线明显分离.

综上所述，存在地形情况下，地面测量、空中随地形测量及空中同一高度测量所得的磁场响应曲线均能够有效反映异常体信息，地形并没有明显影响异常响应能力.

5 总结

本文利用有限元方法计算了带地形的地面和地空电磁响应，并对比了其响应特性.通过对比2.5维正演结果与均匀半空间的解析解，验证了本文算法的精度.通过对比无地形异常体模型的地面和地空电磁场响应发现，地面测量时，Hz的异常响应能力最强，Ey次之，Hx最弱；在低阻体上方，Ey呈现负异常，Hx呈现正异常，Hz分量先呈现正异常后呈现负异常；在高阻体上方，呈现相反的趋势.在空中不同高度测量时，其磁场响应特性与地面测量结果相似，磁场响应曲线与地面测量时的形态相同，异常幅度相近，并没有明显降低，这说明空中测量仍然可以有效识别异常体.本文还对比了带地形无异常体模型的地面和地空电磁场响应特性，分析了地形对地面和地空电磁场响应的影响，结果表明，地面测量时，地形对Hz的影响较大，对Ey和Hx的影响较小；Ey和Hx在凸地形上方呈现正异常，在凹地形上方呈现负异常，其异常变化比较平稳，相对异常幅度较小；Hz在地形边界位置附近异常变化较快，相对异常极值出现在地形边界处，幅度较大.空中随地形测量及空中同一高度测量时，其磁场响应特性与地面测量结果相似，变化趋势及幅度相近，这说明地形对地空磁场响应曲线所产生的影响与地面测量时相比，并没有本质上的明显差异.本文还分析了带地形异常体模型地面和地空电磁响应特性，地面测量时，带地形异常体模型和带地形无异常体模型的电磁场响应曲线在平坦区域几乎重合，在地形边界位置附近形态幅度相近，在异常体上方曲线明显分离；在低阻体上方，Ey呈现负异常，Hx呈现正异常，Hz先呈现正异常后呈现负异常；在高阻体上方，呈现相反的趋势.当空中随地形测量及空中同一高度测量时，带地形异常体模型和带地形无异常体模型的磁场响应曲线与地面随地形测量时特点相似，这说明，在存在地形的情况下，地面测量、空中随地形测量及空中同一高度测量所得的磁场响应曲线均能够有效反映异常体信息，地形并没有明显影响磁场的异常响应能力.

References
Abdallah S, Mogi T, Kim H J. 2017. Three-dimensional inversion of GREATEM data:application to GREATEM survey data from Kujukuri Beach, Japan. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 10(10): 4321-4327.
Cagniard L. 1953. Basic theory of the magneto-telluric method of geophysical prospecting. Geophysics, 18(3): 605-635.
Di Q Y, Unsworth M, Wang M Y. 2004. 2.5-D CSAMT modeling with the finite element method over 2-D complex earth media. Chinese Journal of Geophysics (in Chinese), 47(4): 723-730.
Di Q Y, Wang R. 2008. Controlled Source Audio-Frequency Magnetotellurics (in Chinese). Beijing: Science Press.
Elliott P. 1998. The principles and practice of FLAIRTEM. Exploration Geophysics, 29(2): 58-60.
Ji Y J, Wang Y, Xu J, et al. 2013. Development and application of the grounded long wire source airborne electromagnetic exploration system based on an unmanned airship. Chinese Journal of Geophysics (in Chinese), 56(11): 3640-3650. DOI:10.6038/cjg20131105
Jin J M. 1998. Electromagnetic Finite-Element Method (in Chinese). Wang J G trans. Xi'an: Xidian University Press.
Kang L L. 2019. Research on suppression technology of motion-induce noise in ground-airborne frequency-domain electromagnetic system[Ph. D. thesis] (in Chinese). Changchun: Jilin University.
Li S Y, Lin J, Yang G H, et al. 2013. Ground-Airborne electromagnetic signals de-noising using a combined wavelet transform algorithm. Chinese Journal of Geophysics (in Chinese), 56(9): 3145-3152. DOI:10.6038/cjg20130927
Li X, Zhang Y Y, Lu X S. 2015. Inverse synthetic aperture imaging of ground-airborne transient electromagnetic method with a galvanic source. Chinese Journal of Geophysics (in Chinese), 58(1): 277-288. DOI:10.6038/cjg20150125
Li Y G, Key K. 2007. 2D marine controlled-source electromagnetic modeling:part 1-an adaptive finite-element algorithm. Geophysics, 72(2): WA51-WA62.
Lin J, Kang L L, Liu C S, et al. 2019. The frequency-domain airborne electromagnetic method with a grounded electrical source. Geophysics, 84(4): E269-E280.
Lu X Y, Unsworth M, Booker J. 1999. Rapid relaxation inversion of CSAMT data. Geophysical Journal International, 138(2): 381-392.
Mitsuhata Y. 2000. 2-D electromagnetic modeling by finite-element method with a dipole source and topography. Geophysics, 65(2): 465-475.
Mogi T, Tanaka Y, Kusunoki K I, et al. 1998. Development of grounded electrical source airborne transient EM (GREATEM). Exploration Geophysics, 29(1-2): 61-64.
Mogi T, Kusunoki K I, Kaieda H, et al. 2009. Grounded electrical-source airborne transient electromagnetic (GREATEM) survey of Mount Bandai, north-eastern Japan. Exploration Geophysics, 40(1): 1-7.
Nam M J, Kim H J, Song Y, et al. 2007. 3D magnetotelluric modelling including surface topography. Geophysical Prospecting, 55(2): 277-287.
Newman G A, Alumbaugh D L. 1995. Frequency-domain modelling of airborne electromagnetic responses using staggered finite differences. Geophysical Prospecting, 43(8): 1021-1042.
Pemberton R H, Seige H O, Bosschart R A. 1970. Deep penetrating airborne Turair EM system. Society of Exploration Geophysicists, 40th Annual International Meeting. New Orleans, LA, United States.
Sasaki Y, Nakazato H. 2003. Topographic effects in frequency-domain helicopter-borne electromagnetics. Exploration Geophysics, 34(1-2): 24-28.
Smith R S, Annan A P, McGowan P D. 2001. A comparison of data from airborne, semi-airborne, and ground electromagnetic systems. Geophysics, 66(5): 1379-1385.
Thomson S, Fountain D, Watts T. 2007. Airborne geophysics-evolution and revolution.//5th International Conference on Mineral Exploration. September, Toronto, Canada.
Wu X, Xue G Q, Fang G Y, et al. 2019. The development and applications of the semi-airborne electromagnetic system in China. IEEE Access, 7(99): 104956-104966.
Xu S Z. 1994. The Finite Element Method in Geophysics (in Chinese). Beijing: Science Press.
Xu Y. 2014. Study about 1D forward and inversion of SAEM[Master's thesis] (in Chinese). Chengdu: Chengdu University of Technology.
Yin C C, Zhang B, Liu Y H, et al. 2015. 2.5-D forward modeling of the time-domain airborne EM system in areas with topographic relief. Chinese Journal of Geophysics (in Chinese), 58(4): 1411-1424. DOI:10.6038/cjg20150427
Zhang J F, Zhi Q Q, Li X, et al. 2013. 2.5D finite element numerical simulation for electric dipole source on ridge terrain. The Chinese Journal of Nonferrous Metals (in Chinese), 23(9): 2540-2550.
Zhang Y Y. 2013. Study on multi-component interpretation and full field apparent resistivity of semi-airborne transient electromagnetic method with electric dipole on the surface[Ph. D. thesis] (in Chinese). Xi'an: Chang'an University.
Zhou H G, Lin J, Liu C S, et al. 2016. Interaction between two adjacent grounded sources in frequency domain semi-airborne electromagnetic survey. Review of Scientific Instruments, 87(3): 034503.
底青云, Unsworth M, 王妙月. 2004. 复杂介质有限元法2.5维可控源音频大地电磁法数值模拟. 地球物理学报, 47(4): 723-730.
底青云, 王若. 2008. 可控源音频大地电磁数据正反演及方法应用. 北京: 科学出版社.
嵇艳鞠, 王远, 徐江, 等. 2013. 无人飞艇长导线源时域地空电磁勘探系统及其应用. 地球物理学报, 56(11): 3640-3650. DOI:10.6038/cjg20131105
金建铭. 1998.电磁场有限元方法.王建国译.西安: 西安电子科技大学出版社.
康利利. 2019.地空频率域电磁探测系统运动噪声抑制技术研究[博士论文].长春: 吉林大学.
李肃义, 林君, 阳贵红, 等. 2013. 电性源时域地空电磁数据小波去噪方法研究. 地球物理学报, 56(9): 3145-3152. DOI:10.6038/cjg20130927
李貅, 张莹莹, 卢绪山, 等. 2015. 电性源瞬变电磁地空逆合成孔径成像. 地球物理学报, 58(1): 277-288. DOI:10.6038/cjg20150125
徐世浙. 1994. 地球物理中的有限单元法. 北京: 科学出版社.
许洋. 2014.半航空电磁一维正反演研究[硕士论文].成都: 成都理工大学.
殷长春, 张博, 刘云鹤, 等. 2015. 2.5维起伏地表条件下时间域航空电磁正演模拟. 地球物理学报, 58(4): 1411-1424. DOI:10.6038/cjg20150427
张继锋, 智庆全, 李貅, 等. 2013. 起伏地形电偶源2.5维有限元数值模拟. 中国有色金属学报, 23(9): 2540-2550.
张莹莹. 2013.水平电偶源地空系统瞬变电磁法多分量解释技术及全域视电阻率定义研究[硕士论文].西安: 长安大学.