轨道交通采空区探测中的航空电磁全域视电阻率成像算法研究

贾亚娟　郑建波　周红芳

摘 要：目前采空区探测的方法主要有地震勘探、地面瞬变电磁法，高密度电法等，然而在地形地貌复杂的地区开展传统地面勘探时，工作难度大、效率低、成本高，与上述方法相比，时间域航空电磁法可以快速地、大面积地进行采空区的勘查，工作效率大大提高，特别是在地形条件复杂的山区等，更能凸显其优势。

关键词：采空区 航空电磁 视电阻率 成像

1 引言

近年来，我国的轨道交通有了日新月异的发展，在几条贯穿全国的主干线基础上已经建成了四通八达的交通网。为了使人们的出行更加快速、方便，开展城际高速铁路和城际轨道交通建设成为当前的工作重点。然而部分城市在过去很长一段时间对地下煤炭和矿产资源的开采留下了大量采空区，这对城际轨道交通建设埋下了巨大安全隐患。对采空区进行精细探测成为保障轨道交通安全的关键。然而全域视电阻率成像方法并未应用于时间域航空电磁数据成像当中，与地面和井中电磁法相比，航空电磁法的发射波形十分复杂，不再是简单的阶跃电流波形，而当前的全域视电阻率成像方法并未考虑实际发射电流波形的影响，这将导致时间域航空电磁数据成像精度受到很大影响。为此，本文开展任意发射波形时间域航空电磁数据全域视电阻率成像算法研究，对提高时间域航空的电磁数据解释水平具有积极的促进作用。

2 时间域航空电磁响应模拟

针对当前轨道交通采空区勘探过程中航空电磁观测数据成像技术的不足，开展任意发射波形时间域航空电磁快速成像算法研究，利用时间域褶积算法模拟任意电流发射波形时间域航空电磁响应，基于反函数定理开发考虑发射波形的全域视电阻率成像方法和数据处理软件，最终实现时间域航空电磁数据快速、高精度成像，进而推动轨道交通采空区航空电磁勘探技术的发展。

2.1 时间域航空电磁阶跃响应快速数值模拟

从频率域麦克斯韦方程组出发推导一维介质模型频率域航空电磁响应表达形式，利用汉克尔变换算法获得频率航空电磁响应。采用基于正弦变换的频时转换技术将频率域航空电磁响应转换到时间域，完成阶跃发射波形时间域航空电磁响应快速数值模拟。在阶跃响应模拟过程中考虑航空电磁观测系统的装置形式和飞行高度的影响。

2.2 基于时间域褶积算法的任意电流发射波形时间域航空电磁响应模拟

时间域航空电磁系统实际发射波形十分复杂，而且发射波形对观测响应影响严重，要获得准确的航空电磁成像结果必须考虑发射波形的影响，为此本项目采用时间域阶跃响应与发射电流时间导数褶积的方法模拟任意发射波形时间域航空电磁响应。开发任意发射波形时间域航空电磁响应模拟算法是实现考虑发射波形的全域视电阻率成像算法的基础。

2.3 基于反函数定理的全域视电阻率成像方法

通过分析时间域航空电磁响应随地下介质电性的变化规律，可知地下介质电阻率与时间域航空电磁观测响应间存在一一对应关系。基于这一单调特性，利用反函数定理和泰勒展开公式将积分形式转化为级数形式，进而建立一种基于迭代格式的全域视电阻率成像方法。由于在航空电磁响应数值模拟过程中考虑了航空电磁系统的装置形式、飞行高度和发射电流波形，因此利用其定义的全域视电阻率不再受到装置形式、飞行高度和发射电流波形的影响，从而可以实现快速、高精度成像技术。开发考虑电流发射波形的全域视电阻率定义方法是对时间域航空电磁数据进行高精度成像的关键。

3 技術方案及算法研究

3.1 阶跃波形时间域航空电磁数值模拟方法

从频率域麦克斯韦方程组出发

其中E和H分别为电场和磁场，ω是角频率，σ、μ、ε分别是地下介质电导率、磁导率和介电常数。通过引入谢坤诺夫势和格林函数可以获得一维频率域航空电磁响应：

其中m为发射源磁矩，h为系统飞行高度，z是接收点坐标，rTE是地层对电磁波的反射系数，与地层厚度和电阻率有关，u0是关于频率的函数，J0是0阶贝塞尔函数。利用汉克尔变换技术可以将（3）式写成数值累加形式，通过合理设置汉克尔变换系数即可完成频率域航空电磁响应数值模拟。

利用频率域响应实部和虚部的奇偶性，由傅里叶逆变换可以将频率域航空电磁响应转换到间域：

其中Im[ ]表示取虚部。利用正弦函数和贝塞尔函数之间的关系可以得到时间域航空电磁阶跃响应：

J1/2为0.5阶贝塞尔函数。在次利用汉克尔变换即可完成时间域航空电磁阶跃响应计算。

3.2 任意电流发射波形时间域航空电磁响应模拟

对于任意波形发射电流的时间域电磁响应一般采用输入电流与系统响应进行褶积求得。然而Smith和Lee（2004）通过对航空电磁法的研究发现脉冲响应是不稳定的，当时间t趋近于0时，半空间的脉冲响应趋近于t-1/2，所以脉冲响应在t=0处存在一个奇点，这对于瞬变电磁响应B（t）及其对时间的导数dB/dt在断电瞬间的稳定性有很大影响。为此本项目采用阶跃响应与电流波形时间导数褶积的方法实现任意电流发射波形时间域航空电磁响应模拟。利用脉冲响应与阶跃响应间的微积分关系可得任意电流发射波形时间域航空电磁响应：

上式中是发射电流，是时间域阶跃响应。可以利用高斯积分方法完成褶积计算。

3.3 基于反函数定理的全域视电阻率成像方法

在实际工作过程中，探测结束后观测装置和飞行高度已经成为已知量，此时航空电磁响应可以简写为。通过分析时间域航空电磁响应随地下介质电性的变化关系可知，由于响应是关于电阻率ρ的单调函数，根据反函数定理知，必然存在一个电阻率值ρ唯一地对应着一个值，为此，采用泰勒展开的办法，将磁感应强度的积分表达式展成级数形式，取其线性主部，建立起迭代关系的视电阻率定义式。在初始电阻率的邻域内对进行泰勒展开：

取其线性主部建立迭代关系的视电阻率定义式完成视电阻率计算：

其中。通过设定迭代终止条件：

其中是实测航空电磁数据，是电阻率为的半空间模型航空电磁理论响应，ε是迭代终止误差。

基于反函数定理和泰勒展开公式建立一种基于迭代格式的全域视电阻率成像方法，在此过程中考虑实际航空电磁观测系统工作过程中的装置形式和飞行高度，可以实现时间域航空电磁数据快速成像。在时间域航空电磁全域视电阻率成像过程中，考虑实际电流发射波形的影响，通过时间域褶积算法模拟实际电流发射波形时间域航空电磁响应，实现任意发射波形航空电磁响应视电阻率成像，有效提高了航空电磁成像精度。

4 结论

本文研究的航空电磁全域视电阻率成像算法不受装置形式、飞行高度的影响，同时考虑了实际发射电流波形，可以实现对时间域航空电磁数据快速、高精度成像。该成果不但可以为城际铁路建设过程中的采空区探测提供更加准确的数据解释结果，同样可以应用于矿产勘查、环境地质调查等诸多领域，因此具有比较广泛的应用和推广前景。

陕西省教育厅科研计划项目资助（项目编号：21JK0746）。

参考文献：

[1]Huang H， Rudd J. 2008. Conductivity-depth imaging of helicopter-borne TEM data based on a pseudolayer half-space model[J].Geophysics，73（3）：F115-F120.

[2]曾友強，谢国胜，智庆全等. 2019. 小回线源TEM法全域视电阻率成像在某金矿水害中的应用[J]. 地球物理学进展，34（2）：0826-0831.

[3]雷栋，胡祥云，张素芳.2006.航空电磁法的发展现状[J].地质找矿论丛，21（1）： 40-44.

[4]武军杰，李貅，智庆全等.2017. 电性源地-井瞬变电磁全域视电阻率定义[J].. 地球物理学报， 60（4）：1595-1605.

[5]殷长春，张博，刘云鹤，等. 2015. 航空电磁勘查技术发展现状及展望[J]. 地球物理学报，58（8）：2637-2653.

[6]张莹莹，李貅，姚伟华，智庆全，李佳.多辐射场源地空瞬变电磁法多分量全域视电阻率定义[J].地球物理学报，2015，58（8）：2745-2758.

[7]赵越，李貅，王祎鹏. 2015. 大回线源瞬变电磁全域视电阻率定义[J].地球物理学进展，30（4）：1856-1863.