近年来,由于时域航空电磁探测具有快速、便捷、经济高效等优点,在加拿大、澳大利亚、美国等国家开展了相关的电磁理论和数据处理研究,并进行了广泛应用[1-3],在实际探测方面上趋于成熟.在国内则起步较晚,我国于近年才开展基于直升机的吊舱式时间域航空电磁探测系统的研究.
时间域航空电磁系统,由于发射线圈及发射系统电感的影响,发射电流不能阶跃关断.因此,在数据处理时,主要对Off-time期间数据进行电导率-深度成像和反演[4-8].并且多数系统最早时间取样道从0.2 ms开始,如VTEM 系统,这导致浅部信息丢失.而在三维正演中,主要采用阶跃波替代发射波进行计算,这导致在早中期计算结果与实际产生较大误差.
全波形时间域航空电磁测量,不仅能实现浅层地质体的探测,而且具有解释精度高,对良导体有高分辨率等优点.2008 年,AeroTEM 系统利用On-time期间数据,基于高导电薄板模型计算视电导率[9],在Montcalm Deposit沉积层中成功发现了高导镍矿.Yin对半正弦和梯形波发射时,采用时域卷积方法,计算了均匀半空间模型的全波形航空电磁响应[10].白登海等对梯形波和指数波发射时,采用阶跃波离散方法,计算了均匀半空间模型的地面电磁响应[11].嵇艳鞠等通过电流环等效原理,研究了导电薄板和覆盖层存在时,全波形时间域航空电磁探测的分辨率[12].以上的全波形航电响应研究都是基于一维大地等简单模型下的计算,而对含有三维体的复杂大地模型下的全波形响应研究却相对较少.
对于时间域电磁响应的三维正演计算,主要有积分方程法[13]、有限元法[14-15]和有限差分法.积分方程法只能计算简单的三维体(如球体、板状体等规则模型),对于复杂的不规则大地模型(如垂直接触带模型等),则很难计算;有限元法则主要应用在频率域,在计算时间域电磁探测响应时需要将响应从频域变换到时域;而三维有限差分法,直接基于时间域和空间域离散,因此本文采用三维有限差分法进行时间域电磁计算.
综上所述,针对目前时间域航空电磁探测在正演计算和解释存在的不足,结合全波形探测的优势,在国内外阶跃波发射三维无源有限差分算法基础上,研究了任意波形发射全波形航空时间域响应三维有限差分(Finite-Differencesolutionin Time-Domain, FDTD)数值算法:结合矩形波无源三维有限差分算法[16-18],在Switch-off-time时段对发射电流进行场源离散,等效为若干小梯形脉冲元的和,并将小梯形脉冲元在时间上的电磁场转换成空间中初始电磁场,分时引入迭代方程中,逐渐将有源有限差分方程转化为无源有限差分方程.最后通过不同模型下的数值计算,与积分方程法(IntegralEquationSolution, IES)的对比,验证了全波形三维有限差分计算的正确性.
2 全波形三维有限差分计算航空时间域系统发射电流波形如图 1所示.若发射机导通时间足够长,则发射电流上升沿对下降沿产生感应电动势的影响可以忽略不计.因此,本文忽略发射电流上升沿的影响,研究发射电流关断后的电磁响应,将计算区间分为如图所示Switch-off-time时段和Off-time时段.
由于时间域航空电磁系统接收信号中只有二次场感应电动势V2(t)含有地下有用信息:
(1) |
dHr2(t)/dt为发射电流Ir(t)产生二次场的变化率,μ0 为自由空间磁导率.因此,只需研究针对二次场变化率的全波形三维有限差分数值算法,再由(1)式即可得出V2(t).
在Switch-off-time时段,Ir(t)≠0,计算空间为有源空间,迭代方程为有源有限差分方程.如果不对Ir(t)进行处理,直接计算有源方程,则会使迭代过程十分复杂.因此,为了简化差分方程,本文结合地面阶跃波三维无源有限差分算法,采用分段叠加法,将Ir (t)离散为若干小梯形脉冲元.将Ir (t)在Switch-off-time时段产生的二次场响应,等效为若干小梯形脉冲元完全关断后,在无源空间中二次场响应与若干梯形脉冲在有源空间中二次场响应的和.随着计算时间的增加,逐渐将有源有限差分方程转换为类似阶跃波关断后的无源有限差分方程.对于无源空间的二次场响应,将梯形脉冲元以初始场的形式分时代入无源有限差分方程中,迭代计算;对于梯形脉冲在有源空间中的二次场响应,则通过积分方程法计算.
在Off-time时段,Ir(t)=0,有限差分方程完全变为无源有限差分方程.只需将Switch-off-time时段所有小梯形脉冲元二次场响应叠加,作为Off-time时段初始场,代入无源有限差分方程,即可迭代计算出Off-time时段Ir(t)产生的二次场响应.
3 Switch-off-time二次场的计算 3.1 Switch-off-time场源离散设发射电流Switch-off-time 时段波形如图 2所示,关断时间为Toff.为了将有源有限差分方程转换为无源有限差分方程,需要将Ir(t)进行离散.
在Tn时刻,发射电流等效为
(2) |
(3) |
(4) |
(5) |
(6) |
其中:
Ir (t)在Tn时刻产生的二次场变化率dHr2(Tn)/dt等效为ΔI1(t)、…、ΔIn(t)、ΔIn_ini(t)、…、ΔIn_m(t)和ΔIn_off(t)在Tn时刻产生二次场变化率之和.ΔTn越小,近似程度越高.
(7) |
在Tn时刻,ΔIn_(n+2)(t)…ΔIn_m(t)和ΔIn_off(t)式中第一项为0,只剩下常数项,因此:
(8) |
由(8)式可得:
(9) |
在T1 时刻,ΔI1(t)完全关断.对于脉冲元ΔI1(t)产生的二次场响应dH12(T1)/dt,空间为无源空间,可以将dH12(T1)/dt以初始场的形式,在T1 时刻代入无源三维有限差分方程中迭代计算.
同阶跃波发射类似,dH12(T1)/dt需要用积分方程法来计算.因为在ΔI1 (t)关断以后Δtini内,dH12(t)/dt变化非常剧烈,如果在ΔI1(t)完全关断时刻t1,将dH12(t1)/dt代入有限差分方程进行迭代,为了使迭代稳定不发散,时间步长会非常小.这不但使迭代时间变长,而且对计算机硬件要求非常高.
因此,为了保证迭代稳定,在ΔI1(t)完全关断以后,延时Δtini, 再将dH12(T1)/dt作为初始场代入迭代方程中参与计算.如果ΔI1(t)关断时间Δt1 和Δtini足够小,则可以在均匀半空间假设下,用积分方程法计算dH12(T1)/dt.
为了满足均匀半空间的假设,且不产生数值色散.根据Hohmann的理论:
(10) |
(11) |
其中,μ1、σ1 和Δd1 分别为空间Yee氏网格中第一层大地网格的磁导率、电导率和最小网格边长.由(11)式可知,可以通过Δd1 来控制计算精度.
对于ΔI1_ini(T1)、ΔI1_2(T1),计算空间为有源空间.因为Δtini 与Δt2 足够小,dH1_ini2 (T1 )/dt和dH1_22(T1)/dt可以通过在均匀半空间条件下的积分方程法计算.
同理,在Tn时刻,对于ΔI1(t)、…、ΔIn-1(t)为无源空间,二次场变化率dH12(Tn)/dt、…、dH(n-1)2(Tn)/dt,可以通过无源三维有限差分方程迭代计算.dHn2(Tn)/dt、dHn_ini2(Tn)/dt和dHn_(n+1)2(Tn)/dt通过积分方程法计算.随着计算时间的增加,逐渐将有源有限差分方程转换为无源有限差分方程.
3.3 初始二次场的计算在Switch-off-time时段,在Tn时刻,ΔIn(t)产生的初始二次场dHn2(Tn)/dt为
(12) |
其中,Hi2(t)为系统单位冲击信号产生的二次磁场响应,Hs2(t)为系统单位阶跃信号产生的二次磁场响应.
同理,ΔIn_ini(t)和ΔIn_(n+1)(t)在Tn时刻产生的二次场响应为
(13) |
(14) |
根据频率域下,航空垂直回线源在均匀半空间中脉冲响应磁场表达式,利用汉克尔变换[19]、Guptasarma数字滤波方法[20]以及贝塞尔函数展开法与小宗量近似结合[21],可以用数值方法计算出Hs2(t),再代入(12)、(13)、(14)式则可计算出dHn2(Tn)/dt、dHn_ini2(Tn)/dt和dHn_(n+1)2(Tn)/dt.
3.4 Switch-off-time二次场的迭代计算在转换为无源有限差分方程以后,计算流程同阶跃波三维无源有限差分方程相同,只需再将由前述计算的初始场分时代入,即可迭代计算出Switch-off-time时段dHr2(t)/dt.无源有限差分方程及其计算流程见参考文献[16-18],这里就不多加叙述.
4 Off-time二次场的计算在Off-time时段,Tn时刻,当Toff ≤Tn≤Toff+Δtini 时,(9)式变为
(15) |
其中:
(16) |
(17) |
其中:
由(16)式得:
(18) |
dH12(Tn)/dt、…、dH(n-1)2(Tn)/dt,由迭代计算,dHn2(Tn)/dt、dHn_ini2(Tn)/dt由积分方程法计算.
当Toff+ΔTini ≤Tn时:
(19) |
至此,dHr2(t)/dt完全转换为若干小梯形脉冲在无源空间二次场响应的和,以初始场的形式叠加后,代入无源有限差分方程中迭代计算.
5 模型分析为了计算方便,下面通过梯形波发射,对比积分方程法,在均匀半空间模型下验证全波形航空时间域三维FDTD 的正确性.由于均匀半空间条件下,采用积分方程法计算时,航空时间域电磁响应理论解没有解析式,只能通过数值方法计算,只有地面时间域电磁响应理论解可以由解析式计算.因此,本文首先通过地面时间域响应解析解与数值解的对比,来验证数值方法的准确性,再通过全波形航空三维FDTD 与积分方程法数值解的对比,间接验证全波形航空三维FDTD 的正确性.然后由阶跃波发射与不同关断时间梯形波发射电磁响应的对比来说明全波形FDTD 的重要性.通过地下垂直接触带模型的航空电磁响应三维分布来验证算法对于复杂大地条件下任意源航空时间域二次场响应计算的有效性.最后根据高低阻三维体模型,研究全波形响应信号在Switch-off-time和Off-time时段对地下异常的分辨问题.表 1是计算参数表.
FDTD 迭代网格采用变间距步长,网格数为101×101×50=510050 个,计算空间为2740 m×2740m×2350m.
5.1 均匀半空间图 3、图 4和图 5分别是均匀半空间条件下,采用解析解和数值解计算,发射线圈中心点地面时间域二次场感应电动势衰减曲线和主剖面感应电动势分布.均匀半空间电阻率10Ωm, 发射电流关断时间300μs.
由图 3、图 4和图 5可以看出,采用积分方程法计算时,发射线圈中心点二次场感应电动势衰减曲线和主剖面感应电动势分布数值解同解析解的误差都小于1%,可以采用数值解代替解析解来验证航空全波形FDTD 的正确性.
图 6、图 7和图 8分别是在均匀半空间条件下,采用IES数值解和FDTD 计算,航空发射线圈中心点二次场感应电动势衰减曲线和主剖面感应电动势分布.
由图 6中可以看出,梯形波发射时,在Switch-off-time时段,发射线圈中心点二次场感应电动势为负;在Off-time时段,发射线圈中心点二次场感应电动势为正,与Switch-off-time时段相反.在计算时间内,FDTD计算结果与IES计算误差小于1%.
由图 7 和图 8 中可以看出,在Switch-off-time时段,“烟圈"内二次场感应电动势为负;在Off-time时段,“烟圈"内二次场感应电动势为正.随着时间的增加,“烟圈"逐渐向外扩散.在剖面中心部分,FDTD 与IES计算的二次感应电动势能很好重合.但在网格边界处,由于边界处的截断误差,随着计算时间的增加,响应存在一定误差.这可以通过加入完全匹配层等边界条件改善.
图 9为不同关断时间均匀半空间航空发射线圈中心点感应电动势衰减曲线.toff等于0时为阶跃波发射衰减曲线.由图中可以看出,在Switch-off-time期间,阶跃波发射与梯形波发射两者的二次场不仅符号相反,而且阶跃波发射的二次场比梯形波发射的二次场幅值大;在Off-time期间,梯形波发射二次场幅值比阶跃波发射二次场幅值大,在晚期时二者才相互重合.这也导致在数据解释时,如果采用阶跃波发射做全程视电阻率解释,无论是在Switch-off-time期间还是在Off-time早期,视电阻率都会产生误差;只有在晚期,视电阻率才和实际相同.
为了进一步说明FDTD 对任意复杂大地航空时间域响应计算的有效性,计算了采用IES 不容易计算的,在低阻覆盖层下,垂直接触带中三维低阻体模型下的二次场感应电动势.
图 10为发射线圈中心点感应二次电动势衰减曲线与均匀半空间条件下,发射线圈中心点二次感应电动势衰减曲线的对比.其中均匀半空间电阻率为10Ωm.由图中可以看出,无论是在Switch-off-time期间还是Off-time期间,早期,由于低阻覆盖层的影响,感应电动势衰减较慢;晚期,由于高阻背景的影响,发射线圈中心点感应电动势衰减较快.
图 11(a—c)分别是不同时刻发射线圈平面二次感应电动势分布图.由图中可以明显看出,无论是在Switch-off-time期间还是Off-time期间,由于低阻体的存在,“烟圈"在低阻体附近扩散较慢,导致低阻体附近感应电动势幅值低于其他地方.各个时段早期,由于低阻覆盖层的影响,“烟圈"扩散较慢;晚期由于高阻围岩的影响,“烟圈"扩散较快,感应电动势对低阻异常的反应也相对较慢,符合实际情况.
综上所述,无论是一维均匀大地还是三维复杂大地,无论是在Switch-off-time期间还是Off-time期间,全波形时间域航空电磁响应三维有限差分数值计算都能有效地模拟二次场信号在空间中的传播和随时间的变化.
5.3 均匀半空间中有低阻异常图 12是均匀半空间中含有低阻体时,发射线圈中心点感应电动势衰减曲线与均匀半空间感应电动势的对比.由图中可以看出,无论是在Switch-off-time时段还是Off-time时段,瞬变电磁响应对低阻异常的反应比较明显.在Off-time期间,随着时间的增大,衰减曲线对低阻异常的反映由小变大,然后变小,最后与均匀半空间响应相重合.因此,对于均匀半空间中低阻异常的检测,如地下水以及金属矿等低阻体的探测等,接收机只有在异常体附近,在Switch-off-time和Off-time期间的某一个时间段内可以分辨出低阻异常.
图 13a是均匀半空间中含有高阻体时,发射线圈中心点二次场感应电动势与均匀半空间感应电动势的对比.由图中可以看出,衰减曲线对高阻体的反应在Switch-off-time期间比较明显.
图 13b 为在Off-time 期间,衰减曲线的放大图.由图中可以看出,与均匀半空间衰减曲线相比只有较小偏差.所以,对于均匀半空间中高阻异常的检测,如果只处理Off-time以后的数据,对高阻异常的分辨不如Switch-off-time时段明显.因此,探测高阻空洞时,在Switch-off-time时段探测比在Off-time时段容易.
6 结论本文将发射源进行离散,通过初始场的形式将有源有限差分方程转化为无源有限差分方程,得出任意源全波形三维有限差分方程.通过不同模型的计算分析可知:
(1) 全波形三维有限差分算法,对于任意源发射,任意复杂大地,可以计算空间中任意点,在Switch-off-time和Off-time期间的感应电动势.在梯形波发射,均匀半空间模型下,计算稳定时间大于10ms, 发射线圈中心点二次感应电动势衰减曲线平均误差小于1%.
(2) 如果用阶跃波发射代替梯形波发射进行正演计算和数据解释,在Switch-off-time期间,二者符号相反,前者幅值比后者大;在Off-time期间,两者符号相同,前者幅值比后者幅值小.如果做全程视电阻率解释,采用阶跃波近似,则无论是在Switch-off-time期间还是在Off-time早期,视电阻率都会产生误差,只有到晚期才和实际重合.由此也证明了全波形FDTD 的重要性.
(3) 对于含有低阻体均匀半空间,二次场感应电动势在Switch-off-time和Off-time期间,在异常体附近,对异常的反应都较为明显.随着时间的增加,仪器对异常的分辨率由小变大,然后再变小.只有在某个时段内,才可以比较容易地分辨出低阻异常.
(4) 对于含有高阻体的均匀半空间,在Switch-off-time期间,感应电动势对异常的反应比在Off-time期间明显.因此,对于高阻体的检测,可以通过在Switch-off-time期间,剔除接收机接收信号中的一次场来分辨出高阻体[22].
(5) 全波形FDTD 在计算梯形波、三角波等关断沿线性变化的航空时间域响应时,计算时间较短,精度较高;在计算非线性关断源(正弦波等)航空电磁响应时,在Switch-off-time期间,非线性关断源二次场响应采用梯形脉冲二次场响应等效,产生的失真可以根据非线性源与线性源的误差来计算.ΔTn越小,误差越小,同时迭代计算时间也越长.ΔTn的选择可以根据关断源具体计算.
(6) 对于发射信号上升沿对下降沿二次场响应的影响不能忽略时,只需将场源离散应用于上升沿,即可计算出上升沿产生的二次场响应.
全波形时间域航空电磁响应三维有限差分数值计算的实现,可以更加实际地研究任意波形发射时,时间域航空电磁响应信号在大地中的传播特点,为仪器系统的设计提供理论支持,为将来研究任意波形三维全波形反演奠定基础.
致谢感谢吉林大学地球信息探测仪器教育部重点实验室为本文提供了研究平台,感谢吉林大学时间域航空电磁组全体成员对本研究的支持和帮助,感谢审稿者对本文的快速认真审阅.
[1] | Fountain D. Airborne electromagnetic systems-50years of development. Exploration Geophysics , 1998, 29(2): 1-11. DOI:10.1071/EG998001 |
[2] | 张昌达. 航空时间域电磁法测量系统: 回顾与前瞻. 工程地球物理学报 , 2006, 3(4): 265–273. Zhang C D. Airborne time domain electromagnetics system: look back and ahead. Chinese Journal of Engineering Geophysics (in Chinese) , 2006, 3(4): 265-273. |
[3] | Fountain D. 60 years of airborne EM-focus on the last decade.//AEM2008-5th International Conference on Airborne Electromagnetics. Haikko Manor. Finland:, 2008. |
[4] | Macnae J C, Xiong Z. Block modelling as a check on the interpretation of stitched CDI sections from AEM data. Exploration Geophysics , 1998, 29(2): 191-194. DOI:10.1071/EG998191 |
[5] | Eaton P A. Application of an improved technique for interpreting transient electromagnetic data. Exploration Geophysics , 1998, 29(2): 175-183. DOI:10.1071/EG998175 |
[6] | Huang H P, Rudd J. Conductivity-depth imaging of helicopter-borne TEM data based on a pseudolayer half-space model. Geophysics , 2008, 73(3): F115-F120. DOI:10.1190/1.2904984 |
[7] | Sattel D. Conductivity information in three dimensions. Exploration Geophysics , 1998, 29(2): 157-162. DOI:10.1071/EG998157 |
[8] | Sattel D. Inverting airborne electromagnetic (AEM) data with Zohdy's method. Geophysics , 2005, 70(4): G77-G85. DOI:10.1190/1.1990217 |
[9] | Walker S E, Rudd J. Advanced processing of helicopter TEM Data. 5th International Conference on Airborne Electromagnetics, Extended Abstract,AEM2008,2008. |
[10] | Yin C, Smith R S, Hodges G, et al. Modeling results of on-and off-time B and dB/dt for time-domain airborne EM system. 70th Annual Conference and Exhibition, Extended Abstract, EAGE, 2008. H048. |
[11] | 白登海, MejuM. 瞬变电磁法中两种关断电流对响应函数的影响及其应对策略. 地震地质 , 2001, 23(2): 245–251. Bai D H, Meju M. The effect of two types of turn-off current on TEM responses and the correction techniques. Seismology and Geology (in Chinese) , 2001, 23(2): 245-251. |
[12] | 嵇艳鞠, 栾卉, 李肃义, 等. 全波形时间域航空电磁探测分辨率. 吉林大学学报 (地球科学版) , 2011, 41(3): 885–891. Ji Y J, Luan H, Li S Y, et al. Resolution study of full-waveform detection in airborne TEM. Journal of Jilin University (Earth Science Edition) (in Chinese) , 2011, 41(3): 885-891. |
[13] | 殷长春, 刘斌. 瞬变电磁法三维问题正演及激电效应特征研究. 地球物理学报 , 1994, 37(Suppl. Ⅱ): 486–492. Yin C C, Liu B. The research on the 3D TDEM modeling and IP effect. Chinese J. Geophys. (in Chinese) , 1994, 37(Suppl. Ⅱ): 486-492. |
[14] | 王华军, 罗延钟. 中心回线瞬变电磁法2. 5维有限单元算法. 地球物理学报 , 2003, 46(6): 855–862. Wang H J, Luo Y Z. Algorithm of 2. 5 dimensional finite element method for transient electromagnetic with a central loop. Chinese J. Geophys. (in Chinese) , 2003, 46(6): 855-862. |
[15] | 严良俊, 徐世浙, 陈小斌, 等. 线源二维瞬变电磁场的正演计算新方法. 煤田地质与勘探 , 2004, 32(5): 58–61. Yan L J, Xu S Z, Chen X B, et al. Forward modeling method of 2D transient electromagnetic field induced by line source. Coal Geology & Exploration (in Chinese) , 2004, 32(5): 58-61. |
[16] | Wang T, Hohmann G W. A finite-difference, time-domain solution for three-dimensional electromagnetic modeling. Geophysics , 1993, 58(6): 797-809. DOI:10.1190/1.1443465 |
[17] | 闫述, 陈明生, 傅君眉. 瞬变电磁场的直接时域数值分析. 地球物理学报 , 2002, 45(2): 275–283. Yan S, Chen M S, Fu J M. Direct time-domain numerical analysis of transient electromagnetic fields. Chinese J. Geophys. (in Chinese) , 2002, 45(2): 275-283. |
[18] | 宋维琪, 仝兆歧. 3D瞬变电磁场的有限差分正演计算. 石油地球物理勘探 , 2000, 35(6): 751–757. Song W Q, Tong Z Q. Forward finite differential calculation for 3-D transient electromagnetic field. Oil Geophysical Prospecting (in Chinese) , 2000, 35(6): 751-757. |
[19] | Guptasarma D, Singh B. New digital linear filters for Hankel J0 and J1 transforms. Geophysical Prospecting , 1997, 45(5): 745-762. DOI:10.1046/j.1365-2478.1997.500292.x |
[20] | Guptasarma D. Computation of the time-domain response of a polarizable ground. Geophysics , 1982, 47(11): 953-963. |
[21] | 许洋铖, 林君, 嵇艳鞠, 等. 航空时间域电磁法回线源有限差分初始场计算. 电波科学学报 , 2010, 25(2): 259–264. Xu Y C, Lin J, Ji Y J, et al. Calculation of initial field for loop source in airborne time-domain electromagnetic by finite-difference approach. Chinese Journal of Radio Science (in Chinese) , 2010, 25(2): 259-264. |
[22] | 嵇艳鞠, 林君, 于生宝, 等. ATTEM系统中电流关断期间瞬变电磁场响应求解的研究. 地球物理学报 , 2006, 49(6): 1884–1890. Ji Y J, Lin J, Yu S B, et al. A study on solution of transient electromagnetic response during transmitting current turn-off in the ATTEM system. Chinese J. Geophys. (in Chinese) , 2006, 49(6): 1884-1890. |