航空电磁法作为一种经济、高效的勘探手段在地质填图、油气勘探、矿产普查及环境和工程调查等方面得到了广泛应用.航空电磁解释通常使用水平层状模型(Yin and Hodges, 2007; Siemon et al., 2009; Vallee and Smith, 2009;蔡晶等,2014).近年来,地形对航空电磁响应的影响吸引了世界范围内众多学者的关注. Liu和Becker(1992)使用边界元方法计算了带地形频率域航空电磁响应;Newman和Alumbaugh(1995)使用有限差分法模拟了二维和三维简单地形的频率域航空电磁响应;Avdeev等(1998)使用积分方程法模拟了三维航空电磁响应;Sasaki和Nakazato(2003)计算了二维和三维带地形频率域航空电磁响应.然而,所有这些研究都是针对频率域问题,对于带地形时间域航空电磁正演模拟研究较少.随着时间域航空电磁法越来越广泛地应用于地球物理勘查,各种解释手段也得到迅速发展,地形对电磁响应的影响不容忽视(Annetts et al., 1998;Liu and Becker, 1992; Mitsuhata,2000; Baba and Seama, 2002; Sasaki and Nakazato, 2003; Nam et al., 2007; 刘云鹤和殷长春,2013;蔡晶等,2014).
带地形航空电磁响应可以利用各种数学方法进行模拟.边界元法虽然可以模拟起伏地形,但不适用于复杂模型;有限差分法模拟带地形航空电磁响应时必须使用规则的楼梯来模拟起伏地形.为了更好地拟合地形,必须进行精细网格划分,从而导致了计算资源的较大浪费.为解决以上问题,本文提出了一种基于非结构化网格的有限元方法模拟带地形时间域航空电磁响应.为了把2.5维问题转化为二维问题,对麦克斯韦方程进行傅里叶变换并在波数域中对其进行求解.我们使用伽辽金方法(金建铭,1998;张继锋等,2013)将波数域中的物理方程离散化,并利用MUMPS求解器求解离散后的有限元方程.通过与半空间模型的解析结果以及和已发表文献模型结果进行对比验证本文的算法精度.最后,我们以水平共面(HCP)和直立共轴(VCA)装置系统为例,分析了山峰和山谷地形对航空电磁响应的影响特征.
2 正演算法为了避免奇异性,我们将电磁场分解成背景场和二次场进行求解.二次场满足如下形式的麦克斯韦方程:
式中σ表示电导率,σs=σ-σp(σp表示背景电导率),E p表示背景电场(本文使用源在全空间的电磁响应作为背景场).将式(1)和(2)按分量展开并假设x轴沿走向方向,对 Es,H s各分量以及(1)(2)式中的右端源项在x方向进行傅里叶变换,可得到波数域中麦克斯韦方程组:
对(3)—(8)式进行化简,消去y和z方向的电、磁场分量,得到波数域中电磁场x分量满足的耦合方程 式中,^表示波数域场,kx表示沿x方向的波数,=iωμ0,ke2=kx2-σ.为了得到离散化有限元方程,我们将整个计算区间划分成小单元,并使用伽辽金方法得到每个小三角形单元对应(9)和(10)式的弱形式
式中Ne为小单元对应的基函数(徐世浙,1994).(11)和(12)式的矩阵形式如下: 其中 K ie是3×3的小矩阵,B ie表示源向量.将所有小单元对应的系数矩阵 K ie和右端的源向量组合起来,并假设当外边界足够远时,场值为零(=0,=0),得到
式中 K 表示整体刚度矩阵,e表示待求解的波数域中的电磁场,b 表示源项.利用MUMPS求解器在求解方程时,将系数矩阵的分解与方程求解分开处理.当方程系数矩阵不变而仅源项变化时可以快速重复求解方程,因而特别适合求解航空电磁正反演问题.使用MUMPS求解器求解方程(15),可以得到各个节点上的和,其他场分量可以由这两个分量的导数求得(Mitsuhata,2000;范翠松等,2014).利用反傅里叶变换可以得到空间域的电磁场值(徐世浙,1988;Yin et al., 2008; 范翠松,2013; 底青云等,2004),进而利用正弦变换(等价于0.5阶汉克尔变换),可得到时间域航空电磁系统响应 其中,J表示贝塞尔函数. 3 结果验证
下面我们将本文模拟结果与Sasaki和Nakazato(2003)使用有限差分计算的二维带地形频率域航空电磁响应进行对比,以验证本文算法精度;对于时间域航空电磁系统,我们将本文的结果与半空间模型 的半解析解进行对比.图 1a给出Sasaki和Nakazato(2003)计算的梯形山峰模型(简称模型一).山峰高 50 m,顶部和底部分别宽20 m和220 m;模型的电阻率为100 Ωm;收发装置采用水平共面装置,收发 距10 m,发射频率为16 kHz,飞行高度固定为30 m.图 1b给出了本文模拟结果与Sasaki和Nakazato(2003)给出结果之间的比较,图 1c给出两者之间的相对误差.从图中可以看出,两者吻合较好,最大相对误差不超过3%.
由于带地形模型时间域航空电磁响应正演模拟结果发表的较少,本文将计算结果与半空间模型的半解析解(殷长春等,2013)进行比较.收发装置采用水平共面装置,收发距为10 m,地下电阻率为10 Ωm,飞行高度为30 m.发射波形为阶跃波.正演模拟结果对比如图 2a和2c所示,相对误差如图 2b和2d所示.从图中可以看出,两组模拟结果的最大相对误差不超过1.5%.由此得出结论:本文的正演模拟方法具有较高的计算精度.
下面我们通过模拟山峰和山谷地形频率域和时间域航空电磁响应,分析地形对航空电磁响应的影响特征.山峰和山谷地形模型及其网格划分如图 3a和3b所示(简称模型二和模型三).山峰高度和山谷深度均为20 m,在顶部和底部的宽度分别为60 m和220 m.模型的电阻率为100 Ωm.对于频率域航空电磁系统,本文对水平共面和直立共轴两种装置 计算了380 Hz,1600 Hz,6300 Hz,25 kHz,120 kHz 五个频率的航空电磁响应.对于时间域航空电磁系统,本文计算了垂直磁偶极源发射条件下磁场Bz和磁感应dBz/dt.发射和接收机高度均为30 m,收发距为10 m.发射波形为阶跃波,我们在10-5~10-2s之间选取了12个时间道.图 3给出频率域正演模拟结果,图 4给出时间域正演模拟结果.
从上面的模拟结果和响应特征分析可以看出:
(1)地形对频率域和时间域航空电磁响应均有 较大影响,航空电磁数据处理解释不能忽略地形效应;
(2)频率越高或时间道越早,电磁响应受地形的影响越大.这是因为高频段电磁信号穿透能力较弱,反映的信息主要来自距离航空电磁收发装置较近的地表介质;而低频段电磁信号穿透能力较强,反映的信息主要来自距离航空电磁收发装置较远的深部地下介质.由此可以理解地形主要对高频段电磁信号影响大.与之相对应,早期道电磁信号反映的信息来自距离航空电磁收发装置较近的地表介质,而晚期道信号反映的信息主要来自离航空电磁收发装置较远的深部地下介质.因此地形对早期时间道电磁信号影响更大;
(3)带地形航空电磁响应的曲线形态与地形之间存在镜像关系.在地形的拐点位置,电磁响应显示出剧烈的变化.
5 地形对异常体响应的影响为分析地形对异常体响应的影响,我们设计了山峰和山谷地形下埋藏低阻异常体的模型(简称模型四和模型五).如图 5a和5b,山峰和山谷模型参数与图 3中的模型参数相同.异常体埋在距地表 30 m的位置,大小40 m×30 m,电阻率为10 Ωm.围岩电阻率为100 Ωm.图 5和6分别给出频率域和时间域正演模拟结果.为便于分析,我们同时给出了水平地表半空间中埋有相同异常体的模型响应.
图 7a和7b分别给出了山峰和山谷地形下埋有两个异常体的模型(简称模型六和模型七).山峰和山谷模型参数与图 3中的模型参数相同.两个异常体的大小均为40 m×30 m.左侧的异常体埋深较浅,位于斜坡中点下方20 m处;右侧的异常体埋深较深,位于斜坡中点下方30 m处.两个异常体的电阻率均为10 Ωm.围岩电阻率为100 Ωm.图 7和8分别给出这两个模型的频率域和时间域正演模拟结果.为便于分析,我们同时给出了水平地表半空间中埋有相同异常体的模型响应.
对比频率域/时间域响应可以发现,纯地形与带地形异常体的响应曲线在高频段和早期时间道差异较小;然而,随着频率的降低或时间的延长,两者差异逐渐变大,地下电性不均匀体的响应突出表现出来.这从物理上很容易解释.事实上,在高频段和早期时间道,电磁信号穿透浅,主要反映地表信息,此时地形响应占主导地位,而异常体的响应微弱、难以识别;随着频率降低,或时间向晚期道推移,电磁信号逐步穿透地下介质,地表地形响应减弱,而异常体的响应增强.由此,可以得出结论,地形效应主要体现在频域电磁信号的高频段或时域电磁信号的早期道,而有一定埋深的目标体的电磁响应主要体现在中低频和中晚期时间道(取决于良导体导电性和埋深).由此,通过分析不同频段或不同时间道电磁信号,我们可有效识别地形效应和异常体响应.
为了便于分析异常体响应相对整个模型响应的大小,了解地形对航空电磁响应的影响程度,进而确定地形响应和异常体响应之间的量化关系,我们对异常体响应和纯地形影响进行计算,并和整个模型响应进行对比.针对模型四至七,图 9给出了频率域异常体响应相对于不含异常体模型响应的比值曲线;图 10给出了频率域纯地形响应相对于不含地形响应的比值曲线;图 11给出了时间域异常体响应相对于不含异常体模型响应的比值曲线;图 12给出了时间域纯地形响应相对于不含地形响应的比值曲线.
从图 5,7,9和10给出的频率域正演模拟结果、地形相对影响曲线和异常体相对响应曲线,可以看出:在380 Hz、1600 Hz等低频信号中,异常体相对响应较大;而在25 kHz、120 kHz等高频信号中异 常体相对响应较小;地形的相对影响与异常体响应相反,在25 kHz、120 kHz等高频信号,地形相对影响较大;而在380 Hz、1600 Hz等低频信号地形相对影响较小.这是因为地表地形距离航空电磁收发装置更近,高频段电磁信号穿透能力较弱,主要反映地形影响;随着频率降低,低频段电磁场穿透深度增大,有一定埋深的异常体响应得到明显反映.
从图 6,8,11和12给出的时间域正演模拟结果、地形相对影响曲线和异常体相对响应曲线,可以看出:(1)异常体在第三和第四时间道得到较好反应.这是由于更早时间道电磁信号尚未传播到异常体,而更晚期时间道则已穿透异常体.(2)地形在早期时间道得到良好的反映,随时间延迟的增加地形影响减小.这是由于早期信号(对应于高频段)反映距离航空电磁收发装置较近的地表介质-地形信息;随着时间的延迟,电磁信号反映离航空电磁收发装置较远的地下介质信息,地形影响减弱.
比较异常体的相对响应和地形的相对影响可以进一步发现:地形对电磁信号的影响与异常体的响应在相同数量级上,因此,地形对航空电磁信号的影响是不可忽视.
6 结论本文利用非结构有限元法很好地模拟了带地形模型的航空电磁系统响应.通过与均匀半空间的半解析解及已发表模型结果的对比,验证了本文算法具有较高的计算精度.
本文计算结果表明,地形对航空电磁响应影响很大,在地形复杂地区从事航空电磁观测,地形效应不容忽视,尤其是在频率域系统的高频段或时间域系统的早期时间道.对山脊和山谷地形正演模拟结果表明,在地形变化的突变点位置,电磁响应出现快速的变化.在这些地形突变区进行航空电磁数据解释时,应对地形影响给予高度重视.比较频率域和时间域系统的正演结果,可以进一步发现地形对频率域航空电磁响应主要发生在高频段,而对时间域航空电磁系统的影响主要发生在早期时间道.相比之下,具有一定埋深的地下良导体在高频段响应反应较弱,主要出现在中低频和中晚期时间道(取决于良导体导电性和埋深).这在一定程度上为我们有效识别地形效应和异常体响应提供可能.然而,由于地形响应和异常体响应处于相同数量级,为有效实现航空电磁数据反演解释,对地形影响进行模拟和校正十分必要.
致谢 我们向贲放、黄威、齐彦福、任秀艳、黄鑫和邱长凯在模型计算和论文准备过程中给予的帮助表示感谢.特别对审稿人提出宝贵的修改意见表示感谢,最后感谢编辑们的辛苦工作.[1] | Annetts D, Sugeng F, Raiche A. 1998. The effect of topography on airborne electromagnetic response. SEG Expanded Abstracts. |
[2] | Avdeev D B, Kuvshinov A V, Pankratov O V, et al. 1998. Three-dimensional frequency-domain modeling of airborne electromagnetic responses. Exploration Geophysics, 29(2):111-119. |
[3] | Baba K, Seama N. 2002. A new technique for the incorporation of seafloor topography in electromagnetic modelling. Geophysical Journal International, 150(2):392-402. |
[4] | Cai J, Qi Y F, Yin C C. 2014. Weighted Laterally-constrained inversion of frequency-domain airborne EM data. Chinese Journal of Geophysics (in Chinese), 57(3):953-960, doi:10.6038/cjg20140324. |
[5] | 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. |
[6] | Fan C. 2013. Research on complex resistivity forward and inversion with finite element method and its application[Ph. D. thsis](in Chinese). Jilin:College of Geo-exploration Sciences and Technology of Jilin University. |
[7] | Fan C, Li T L, Wang D Y, et al. 2014. Research on 2.5D SIP inversion with topography. Acta Geologica Sinica (in Chinese), 88(4):755-762. |
[8] | Jin J M. 1998. Electromagnetic Finite-Element Method (in Chinese). Xi'an:Xidian University Press. |
[9] | Liu G M, Becker A. 1992. Evaluation of terrain effects in AEM surveys using the boundary element method. Geophysics, 57(2):272-278. |
[10] | Liu Y H, Yin C C. 2013. 3D inversion for frequency-domain HEM data. Chinese Journal of Geophysics (in Chinese), 56(12):4278-4287, doi:10.6038/cjg20131230. |
[11] | Mitsuhata Y. 2000. 2-D electromagnetic modeling by finite-element method with a dipole source and topography. Geophysics, 65(2):465-475. |
[12] | Nam M J, Kim H J, Song Y, et al. 2007. 3D magnetotelluric modelling including surface topography. Geophysical Prospecting, 55(2):277-287. |
[13] | Newman G A, Alumbaugh D L. 1995. Frequency-domain modelling of airborne electromagnetic responses using staggered finite differences. Geophysical Prospecting, 43(8):1021-1042. |
[14] | Sasaki Y, Nakazato H. 2003. Topographic effects in frequency-domain helicopter-borne electromagnetics. Exploration Geophysics, 34(2):24-28. |
[15] | Siemon B, Auken E, Christiansen A V. 2009. Laterally constrained inversion of helicopter-borne frequency-domain electromagnetic data. Journal of Applied Geophysics, 67(3):259-268. |
[16] | Vallee M A, Smith R S. 2009. Inversion of airborne time-domain electromagnetic data to a 1D structure using lateral constraints. Near Surface Geophysics, 7:63-71. |
[17] | Xu S Z. 1988. Selection of wave number K in Fourier inverse transformation for point source and 2-D electric field problem. Computing Techniques for Geophysical and Geochemical Exploration (in Chinese), 10(3):235-239. |
[18] | Xu S Z. 1994. The Finite-element Method in Geophysics (in Chinese). Beijing:Science Press. |
[19] | Yin C C, Hodges G. 2007. Simulated annealing for airborne EM inversion. Geophysics, 72(4):F189-F195. |
[20] | Yin C, Smith R S, Hodges G, et al. 2008. Modeling results of on-and off-time B and dB/dt for time-domain airborne EM systems. Extended Abstract, 70th Annual EAGE Conference and Exhibition, Rome, 1-4. |
[21] | Yin C C, Huang W, Ben F. 2013. The full-time electromagnetic modeling for time-domain electromagnetic systems. Chinese Journal of Geophysics (in Chinese), 56(9):3153-3162, doi:10.6038/cjg20130928. |
[22] | 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. |
[23] | 蔡晶, 齐彦福, 殷长春. 2014. 频率域航空电磁数据的加权横向约束反演. 地球物理学报, 57(3): 953-960, doi: 10.6038/cjg20140324. |
[24] | 底青云, Unsworth M, 王妙月. 2004. 复杂介质有限元法2.5维可控源音频大地电磁法数值模拟. 地球物理学报, 47(4): 723-730. |
[25] | 范翠松. 2013. 基于有限元法的复电阻率正反演研究及应用[博士论文]. 长春: 吉林大学地球探测科学与技术学院, 69-99. |
[26] | 范翠松, 李桐林, 王大勇等. 2014. 起伏地形下复电阻率法2.5维反演研究. 地质学报, 88(4): 755-762. |
[27] | 金建铭. 1998. 电磁场有限元方法. 西安: 西安电子科技大学出版社. |
[28] | 刘云鹤, 殷长春. 2013. 三维频率域航空电磁反演研究. 地球物理学报, 56(12): 4278-4287, doi: 10.6038/cjg20131230. |
[29] | 徐世浙. 1988. 点电源二维电场问题中付氏反变换的波数k的选择. 物探化探计算技术, 10(3): 235-239. |
[30] | 徐世浙. 1994. 地球物理中的有限单元法. 北京: 科学出版社. |
[31] | 殷长春, 黄威, 贲放. 2013. 时间域航空电磁系统瞬变全时响应正演模拟. 地球物理学报, 56(9): 3153-3162, doi: 10.6038/cjg20130928. |
[32] | 张继锋, 智庆全, 李貅等. 2013. 起伏地形电偶源2.5维有限元数值模拟. 中国有色金属学报, 23(9): 2540-2550. |