地空瞬变电磁在积水采空区探测中的应用

侯彦威，高小伟，李雄伟，郭建磊，姜 涛

(中煤科工集团西安研究院有限公司，陕西 西安 710077)

地空瞬变电磁法[1]结合了地面瞬变电磁法和航空瞬变电磁法各自的特点，相比地面瞬变电磁法能有效提高野外工作效率；地空瞬变电磁法采用地面激发、无人机搭载接收探头空中接收的装置形式，有效解决了地表情况复杂地区的瞬变电磁勘探施工问题，目前已逐步被用于环境地质调查[2]、煤田采空区勘查[3，4]等方面。地空瞬变电磁法是Nabighian[5](1988年)基于水平电偶源的工作方式而提出；Smith等[6]早在2001年的时候在加拿大某矿区对地面、地空和航空瞬变电磁法的探测能力进行了对比分析，证明了同等条件下地空瞬变电磁法同样具有勘探深度大、数据信噪比高等优点；Mogi等[7]研制了电性源作为发射源、直升机搭载磁场接收探头的数据采集方式，即电性源地空瞬变电磁；吉林大学[8]研制的无人飞艇长导线源时域地-空电磁勘探系统采集的是磁场的时间导数，促进了地空电磁勘探的进步。在地空瞬变电磁正演及数据处理方面也有大量研究，如李贺和孙怀凤[9，10]分别实现了三维有限元和三维时域有限差分正演；李肃义等[11]研究了小波分析压制噪声；李貅等[12-14]从解释方面提出了全域视电阻率定义、逆合成孔径成像、拟地震成像及联合解释等技术；张澎等[15]将一维自适应正则化反演方法应用于半航空时间域电磁数据反演解释中，并实现了该反演算法的并行计算；王振荣等[4]利用地空时间域电磁系统在哈拉沟煤矿开展采空区勘查工作，采用视电阻率和视深度进行绘图。

综上所述，随着研究的深入，地空瞬变电磁逐步被推广应用，数据反演计算方法也逐步被探索应用。本文采用正反演相结合，通过模拟计算的数据，分析认识所建模型的地空瞬变电磁响应特征，尝试采用约束反演对数据进行反演计算，验证可行，并将其用于实际采空区探测数据处理中，解释成果符合地质电性规律，低阻异常区也获得了验证。

1 地空瞬变电磁工作原理

地空瞬变电磁工作原理与地面瞬变电磁近似，电性源地空瞬变电磁是在地面布置一定长度的导线源，采用大功率设备沿导线向地下发射阶跃变化的电磁波信号，此信号即为激发场源(称为“一次场”)，根据法拉第电磁感应定律，当发射电流被关断时，大地或探测目标体在一次场的作用下，其内部会感生具有空间和时间特性的涡流，其大小与诸多因素有关，如目标体的空间特征和电性特征、激发场的特征等，而且因为热损耗的缘故会逐渐减弱直至消失。通过搭载在无人机上的接收设备接收这种涡流产生的电磁场(称为“二次场”)，二次场具有强弱之分，也具有空间和时间特性，研究这些特性可分析地层或地下目标体的电性，并据此完成探测。

地空瞬变电磁探测装置的飞行高度可灵活调整，一次布设可在短时间内完成大面积勘探。

2 正演理论与模拟验证

在均匀、各向同性、有耗、非磁性、无源媒质中，Maxwell方程组为：

(1)

式中，E表示电场强度；B表示磁感应强度；H表示磁场强度；σ表示介质电导率；ε表示介电常数；t表示时间。

瞬变电磁勘探计算中一般忽略位移电流，而为了获得显式的时间迭代格式，需在上方程组中的式(1b)中加入虚拟介电常数项，在有源区域加入电流密度项，式(1b)即为：

(2)

式中，γ表示具有介电常数的量纲；Js表示源电流密度。

地空瞬变电磁法的激发源一般采用电性源[1]，为便于采用时间域有限差分法进行正演计算，常在位于直角坐标系中晶胞网格的Y轴方向加载电性源。

在地空瞬变电磁正演过程中仍需考虑其在低频情况下正演结果的正确性，因此需对磁场各分量采用低频近似处理，将方程组(1)中的式(a)、(d)和式(2)在直角坐标系下展开，得分项式为：

(3)

(4)

对式(3)和式(4)用差分代替微分，由于Euler前向差分对离散时间步的要求比较严格，故空间离散采用后向差分，均匀网格剖分中电/磁场的时间采样恰好在两相邻磁/电场采样时刻的中心，故时间离散采用中心差分，可以得到各个分量的迭代公式。其中，有源区域中的Ez分量迭代公式为：

(5)

采用第一类边界条件，为保证电磁场传播的因果关系，将需要计算的整个模型剖分得很大且在计算过程中时间域和空间域的稳定性条件(式6)满足时，在边界处可将电场与磁场进行强制性赋零。

(6)

在均匀半空间条件下构建电性发射源与接收试验模型，基于三维时域有限差分算法进行地空瞬变电磁正演计算，将计算结果与解析解进行对比以验证程序的可靠性。发射源与接收点相对位置如图1所示，模型参数为：均匀半空间地层电阻率100Ω·m，电性发射源长度500m，接收点坐标为(x=300m，y=-50m)，接收点高度50m，电流方向沿Y轴正向，电流1A，二次场采样时间10ms，最小网格尺寸为10m且总剖分网格数为221×221×200。将三维时域有限差分模拟计算结果和一维解析解通过直接与误差统计的方式进行对比，绘制对比图(图2)，其中图2(a)为直接对比图，响应曲线形态、特征基本一致；图2(b)为误差对比曲线图，经统计二者的最大误差为2.176%，误差小于1%的占比89.47%，可见三维时域有限差分模拟计算的结果满足精度要求；同时由图2(a)可见，地空瞬变电磁的感应电动势衰减曲线与地面瞬变电磁衰减曲线规律近似，因此可在地面瞬变电磁数据反演方法的基础上对地空瞬变电磁数据进行反演试算。

图1 发射源与接收点相对位置俯视图

图2 瞬变电磁正演模拟与验证对比

3 约束反演方法

地空瞬变电磁数据反演采用自适应的约束反演方法，此方法在反演的过程中采用自主调整正则化因子的方案，反演矩阵条件数比较大时增大正则化因子，当次反演拟合差大于上次反演的拟合差时减小正则化因子；为减小反演过程的病态性，反演的层厚按照等对数间隔进行离散；同时对每一层电阻率进行上下限约束，如果超过给定的上下限，那么就减小反演模型的步长，直到满足阈值[16]。反演目标函数及正则化因子为：

(7)

[(WJ1)TWJ1+μATA]Δmk=(WJ1)TWΔdk

(8)

(9)

4 地质模型地空瞬变电磁响应特征

为进一步对比与验证上述约束反演方法的适用性，特建立地质模型分别如下：

1)均匀半空间模型：电阻率为100Ω·m。

2)与陕北侏罗系煤田纵向整体电性特征相似的K型地质模型(图3a)：第一层模拟第四系地层，电阻率为100Ω·m，厚度50m；第二层假设为侏罗系含煤地层，电阻率为1000Ω·m，厚度为50m；第三层电阻率为100Ω·m。

3)含低阻异常体的K型模型：在侏罗系含煤地层中部加入电阻率为10Ω·m的低阻异常体(假设为积水采空区，见图3a)，低阻异常体体积为100m×100m×40m。

图3 地空瞬变电磁K型模型与计算结果对比

模拟的电性发射源长度500m，接收点高度50m，接收点位于低阻异常体中心上方，电流1A，二次场采样时间10ms，最小网格尺寸为10m且总剖分网格数为221×221×200。

对上述模型进行三维时域有限差分模拟计算，计算结果如图3(b)所示，图中三条衰减曲线分别为均匀半空间模型感应电动势衰减曲线、K型地质模型感应电动势衰减曲线和K型地质模型中间层加入低阻异常体后的感应电动势衰减曲线。可见，三条衰减曲线随时间均呈逐渐衰减趋势，在早期略有交叉；中期K型地质模型感应电动势值较低，加入低阻异常体后的感应电动势值有所增加但仍小于均匀半空间模型感应电动势值；晚期有重合趋势，但数值仍有差异，K型模型模拟计算值始终小于均匀半空间模型值。以上衰减曲线的不同特征反映了不同地质模型的电性分布规律，为反演取得好的结果奠定了基础。

将图3(b)中K型地质模型及其加入低阻异常体后模型的感应电动势数据分别进行约束反演计算，其结果如图3(c)所示。可见，从上而下两种模型的反演曲线在浅部均接近模型第一层电阻率值，随深度增加，电阻率值也在逐渐增大，且逐步跨进第二层高阻层，直至达到极值后又逐渐降低，约在-170m附近接近第三层电阻率值，整体由浅至深均呈“低阻-高阻-低阻”的变化趋势，与K型地质模型三层电性规律吻合较好。加入低阻异常体后的反演曲线所示电阻率值更低，可见该探测方法及约束反演对低阻异常体反应明显，以此通过横向对比可实现积水采空区的探测。

5 积水采空区探测

5.1 测区地形与地空瞬变电磁工作布置

某煤矿地处陕北黄土高原，地形起伏较大、地表沟壑纵横。探测区域内海拔标高最高点达1320m，最低点约为1160m，高差达160m(如图4所示)，地面施工难度较大，所以采用地空瞬变电磁进行勘探。

图4 探测区地形

探测区域与工作布置如图5所示，煤矿生产巷道由西南向东北进入该探测区域(面积约0.7km2)，探测区域前方有小煤窑采空区分布，积水情况不明，采用地空瞬变电磁探测，设计线距40m，点距采用10m，沿设计线路迂回折返飞行，飞行高度为探测区域内最高点以上50m；电性发射源位于测区东南向3km处，长1.5km。

图5 探测工作布置

5.2 地层与地球物理特征

煤矿地层、相应岩性与地球物理特征见表1。

表1 地层岩性及其地球物理特征

由表1可见，探测区域由浅至深地层电性整体呈“低阻—高阻—低阻”的变化特征。当煤层被开采采空区未积水呈更高阻的特征；当采空区有积水时，水体的浸润使得周围介质电阻率明显降低，其电阻率值将明显降低，一般会小于30Ω·m，使得积水采空区与周围介质产生明显的电性差异，上述导电性差异为地空瞬变电磁勘探提供了条件。

5.3 探测效果分析

按照本文前述反演方法对实测数据进行计算，其结果可反映原生地层和被采动地层两种情况，故本文举例两条典型断面如图6所示，图6(a)为反映原生地层电性的3线反演电阻率断面，图6(b)为反映原生地层被采空区破坏且存在低阻异常区的7线反演电阻率断面。断面图横向为横向水平距离，纵向为高程，高程1190m以浅为低阻层，3-3煤层附近为高阻层，深部1110m以下为低阻层，详见断面图右侧柱状电阻率标尺。

图6 地空瞬变电磁探测反演电阻率断面图

由图6可见，反演电阻率横向上整体呈似层状分布，浅部为低阻层，对应新近系以浅地层；中部为高阻层，对应侏罗系含煤地层，探测目标层3-3煤层恰位于本层中部；深部为低阻层，对应三叠系地层。此种由浅至深的电性变化规律与前述地球物理特征一致，表明反演结果反映地层实际。由图6(a)可见含煤地层横向上连续性好，电性分布均匀，表明地层为未受采动破坏的原生地层；图6(b)中部高阻层横向上的连续性在420～800m处被打破，电阻率值突然降低，推断为低阻异常区，同时也表明此处地层应受到了采动破坏且局部积水。

依据3-3煤层底板高程和全区反演电阻率数据，抽取本层电阻率数据并横向取异常，如图7所示，图7中有颜色填充处为低阻区域。

图7 地空瞬变电磁探测3-3煤层顺层平面图

由图7可见，探测区内大部分区域反演电阻率等值线较为稀疏，表明该目标地层大部分区域电性值较为稳定，可以此间接推测等值线稀疏区域目标煤层未受采动破坏；同时也发现一大一小两处低阻异常区，经与后期提供的小煤窑老巷道对比，此两处低阻异常区均与老巷道有对应关系，推测其应为小煤窑沿巷道向周围开采煤层形成的采空区积水所致。从图7中煤层底板等高线分布趋势可见，大范围的低阻异常区恰位于局部低洼区域，在有采空空间存在的情况下，利于水的积聚，所以此大范围低阻异常区推测为采空积水区的依据充分，且经钻探(TC1号钻孔)验证，确为积水采空区。

6 结 论

1)尝试使用地面瞬变电磁数据反演方法对地空瞬变电磁数据进行反演计算，探测结果证明该约束反演算法可以有效的实现地空瞬变电磁数据的处理工作，为地空瞬变电磁数据处理提供一条新的思路。

2)利用地空瞬变电磁法实现煤矿积水采空区探测，将探测结果与采掘巷道资料及钻孔资料对比，验证了地空瞬变电磁法对于煤矿积水采空区探测的可行性与高效性，可为后期煤矿开采提供有效的参考资料。