随着近年来我国经济的快速发展，能源消耗日益增加，煤炭作为我国的重要能源之一,2005—2017年产量由21.5亿t剧增到35.2亿t[1]。煤炭开采产生了大量的采空区，引起坍塌、涌水等地质灾害，严重威胁煤炭的安全生产，如何快速、高效、准确地探测煤炭采空区成为亟待解决的问题。为此，2016年山西省煤炭地质物探测绘院联合吉林大学在沁水煤田晋城区块开展了山西省煤炭采空区地空电磁探测研究。

时域地空电磁系统，融合了时域地面电磁系统和时域航空电磁系统的优势，该系统采用接地长导线为发射天线，飞行器搭载接收装置的工作方式，具有横向分辨率高、探测效率高等优点，适用于沼泽、山地、海陆交互带以及森林覆盖区等恶劣地形条件下的资源快速勘察。MOGI等(2009)将时域地空电磁法应用于日本东北部磐梯山的火山地质结构调查中[2];ITO等(2011)在日本千叶县九十九里滨海岸开展了时域地空电磁探测，进一步了解了这一地区地下水文地质结构[3];ALLAH等(2013)采用时域地空电磁法在日本东南沿海的海陆交互地区开展了探测实验，实验结果反映出了该地区地下500 m的电性结构，并查明了该地区的海水入侵情况[4];嵇艳鞠等(2014)分别在江苏省南通市滨海地区和内蒙古巴彦宝丽格盆地开展了时域地空电磁探测实验，实验结果揭示了这两个地区的地下水资源分布情况[5];ALLAH等(2014)将时域地空电磁法应用于日本千叶县的地质结构调查中，探明了该地区的地下电阻率分布[6]。黎东升等(2017)采用时域地空系统在江苏省燕尾港地区进行了海水入侵的调查[7]。在开展时域地空电磁探测的同时，电磁数据解释方法的相关研究也随之展开。嵇艳鞠、阳贵红等(2012)在地面LOTEM的基础上对时域地空电磁数据进行了解释[8];嵇艳鞠、徐江等(2013)开展了基于神经网络电性源地空视电阻率解释方法的研究[9];李貅等(2015)开展了电性源瞬变电磁逆合成孔径成像方法研究[10]。

笔者采用时域电性源地空电磁系统在山西省沁水煤田开展煤炭采空研究。首先介绍了时域电性源地空电磁系统的组成和工作原理，然后针对时域地空电磁数据的处理和成像方法开展了细致的研究，最终探测结果准确的反映出地下采空区的位置及深度信息。本研究创新性地将时域电性源地空电磁系统应用于煤炭采空电磁探测。

1测区地质描述

晋城区块所处大地构造单元为沁水块坳南部，大型沁水含煤盆地沁水复式向斜南段东翼，晋(城)—获(鹿)褶断带西侧，区域上大致由南向北依次出露地层有奥陶系、石炭系、二叠系，第四系。

石炭系上统太原组和二叠系下统山西组为本区块的主要含煤地层，赋存主要可采的3，9，15号煤层。其中石炭系上统太原组由多层石灰岩、不同粒级的砂岩、灰黑色泥岩及煤层组成，为一套海陆交互相含煤沉积;二叠系下统山西组为一套灰白、灰黑色砂岩、泥岩和煤层组成的陆相含煤沉积[11-13]。本次实验主要探测测区内主采的3号煤层的采空区及采空积水区，埋深约为250～350 m，测区涵盖5个生产矿井及部分关闭小煤窑，3号煤层的采区占据了测区的绝大多数面积，除部分生产矿井已知采空区及采空积水区分布外，小窑采空区的位置几乎没有资料。

本区降水量充足，地表第四系覆盖层主要为黄褐色砂土，山西组中的砂岩往往为含水砂岩，因此煤层经采掘后发生塌陷，破坏了原有地层的水力关系，导致砂岩水通过裂隙进入采空区，甚至形成老空积水区，根据以往在该区开展的电磁勘探工作，本区采空区往往是表现为低阻特征的含水采空区。

2时域地空电磁系统

日本北海道大学的茂木教授最早开展了时域电性源地空电磁系统的理论研究(1998)，并多次将地空电磁系统应用于火山地质结构的探测(2009至今)。从2009年开始，吉林大学、长安大学、成都理工大学和中科院等多家院校及科研单位也开展了地空电磁系统的研究，虽然国内相关研究开展的较晚，但近年来发展迅速。

本次探测中采用的是吉林大学开发的时间域地空电磁系统，其结构如图1所示。该系统采用接地长导线(1～5 km)作为发射天线，车载JDD-100(Jilin Double Domain 100 kW)型大功率发射机为天线提供双极性方波电流，关断时间小于200 μs，发射机最大功率100 kW，最大输出电压1 500 V，最大输出电流60 A，可提供多组发射基频(3.125，6.25，12.5，25 Hz)，系统中采用板状电极代替传统的接地电极，增大电极与地面接触面积。相比于磁性源(不接地回线)，电性源发射可以提供更大的磁矩和收发距，进而提升电磁探测深度。电磁接收机和接收线圈搭载在无人机上，接收机采用24位AD转换器，可以在30 kHz高速采样率下完成30通道电磁信号的模数转换，接收线圈有效面积2 160 m2，动态范围大于160 dB。探测时，通过GPS天线保证发射机和接收机之间的时钟同步。飞行搭载平台为DG654型六旋翼无人机，最大载荷10 kg，最大航速10 m/s。

图1 时域地空电磁系统结构示意
Fig.1 Structural sketch of time-domain geo-spatial electromagnetic system

一次场引起地下围岩和低阻目标体中产生涡流，当发射电流关断后，涡流形成的二次场向上传播被接收线圈采集并记录在接收机中。二次场中携带大量地下介质的电阻率信息，通过对二次场的分析计算可以得到地下的电阻率分布，实现对地下目标体视电阻率和视深度探测的目的。

3时域地空电磁噪声抑制

时域地空系统开展电磁探测时，信号容易受到运动噪声、背景噪声以及天电噪声等影响，降低电磁数据信噪比和成像质量。因此，在计算视电阻率和视深度前对噪声进行抑制是必要的。

3.1基线校正

相比于背景噪声和天电噪声，基线漂移的频率较低(<10 Hz)且没有奇异点。通过小波变换理论可知，低频信号的能量主要集中在小波变换后的高尺度近似分量中，细节分量中几乎不包含低频信号。因此，本研究采用数据小波分解后的高尺度近似分量对基线漂移进行估计，最终实现消除基线校正的目的。

本次探测中，电磁接收机的采样频率为30 kHz，根据奈奎斯特采样定律可知电磁数据中包含的频率信息为0～15 kHz。小波分解实际是对信号在中心频率进行截断，低频信息表征信号的近似部分(尺度系数)，而高频信息表征信号的细节部分(小波系数)。根据电磁数据的频率范围和基线漂移的频率确定小波分解的级数为11级，通过11级分解后的尺度系数对基线漂移进行估计，并从实测电磁数据中减去估计出的基线，得到校正后的数据，基线校正如图2所示。

图2 地空瞬变电磁数据基线校正效果
Fig.2 Base line correction of GREATEM data

3.2噪声抑制

在电磁勘探系统的应用中，电磁数据的信噪比较低会引起成像结果中的假异常，因此抑制噪声是数据处理的关键步骤。在实测GREATEM数据中，经常包含背景噪声、工频噪声及天电噪声等多种干扰，这些噪声可以总结为两类:近似白噪声和电磁尖峰。本文在前人研究的基础上[13-16]，针对GREATEM数据中的多种噪声设计了综合小波去噪方法。

为了表示实测GREATEM数据中的不同分量，笔者将原始数据描述为

G(n)=s(n)+w(n)+e(n)

(1)

式中，G(n)为实测数据;s(n)为纯的GREATEM响应;w(n)为近似白噪声;e(n)为电磁尖峰;n为离散采样点个数。

第1步:利用小波阈值法从实测数据中去除近似白噪声。考虑瞬变电磁法数据的衰减特性和奇异性，我们选择了启发式阈值和硬阈值收缩函数。Sym8小波基具有短支撑、正交和高阶消去矩等重要特性，因此被选用。通过试验确定了分解层数为5。

去除近似白噪声后，将GREATEM数据重写为

(2)

式中，为去除近似白噪声后的电磁数据。

第2步:采用平稳小波变换(SWT)及其逆小波变换(ISWT)进行电磁尖峰抑制。时域地空电磁信号与低通滤波器HP和高通滤波器GP进行卷积(上角标P为滤波器Pass的标识)，被分解为低频近似分量Aj和高频细节分量Dj。高频细节分量被保留，低频近似分量作为下一级频带分解的输入，经过变换后信号被划分为不同的频率成分。

Aj+1=HP*Aj

Dj+1=GP*Dj(j=0,1,…,J)

(3)

平稳小波变换的滤波器HP和GP可以通过修改正交离散小波变换算法获得。将正交离散小波变换的滤波系数表示为hj和gj，在每一层上对hj和gj插值补零得到平稳小波变换的滤波系数和从而获得平稳小波变换的滤波器HP和GP。经过平稳小波变换后，近似分量Aj和细节分量Dj与信号A0等长度，并且在时域上相对应。分解后的数据也可表示为近似系数和细节系数之和，J为最高分解层数,即

(4)

由前面的分析可知，近似分量表征信号的低频部分，而电磁尖峰具有迅速衰减、幅值变化较快的特点，其能量主要集中于高频部分，因此可以通过对细节分量的处理达到抑制电磁尖峰的目的。数据经过平稳小波变换后，关键问题在于识别各层细节分量上的噪声。本文采用基于能量的尖峰信号识别方法[17]:

(5)

其中，m为信号能量二阶矩，n为离散点，通过能量计算和设置合适的判别阈值实现GREATEM数据中的尖峰信号的识别。MALLAT曾指出尖峰信号的影响域在各层细节分量中呈锥形分布，首先计算GREATEM数据中的电磁尖峰长度，随着分解层数的增加，尖峰噪声的长度以1.5倍增长，各层细节分量中的电磁尖峰也随之被识别。将综合小波去噪法应用于GREATEM衰减曲线中，如图3所示，近似白噪声和电磁尖峰均得到有效的抑制。

图3 地空瞬变电磁数据噪声抑制效果
Fig.3 Noise suppression of GREATEM data

4视电阻率及视深度计算

为了将地下介质的电阻率分布表示成视电阻率-视深度图的形式，需要通过电磁数据计算视电阻率和视深度，KAUFAMN(1983)和 WEIR(1980)给出了电偶极源在均匀半空间下Z方向磁场响应[18-19]，即

(6)

式中，Hz为垂直磁场强度;y为垂直方向坐标;I为电流强度;ds为电偶极长度;t为采样时刻;r为接收机到电偶极的距离;θ为电导率σ的函数；erf为误差函数。

(7)

将接收线圈面积S和Bz=μHz(Bz为磁感应强度;μ为磁导率)代入到式(6)，并将多个电偶极的响应进行叠加，得到接地长导线Z方向电磁响应的表达式，如式(8)所示。通过对式(7)，(8)的求解，就可以得到电导率σ。

(8)

式中，ri为接收机到第i个电偶极的距离；xi为第i个电偶极与长导线源中点O之间的距离。

采用扩散深度公式(9)对视深度进行计算，其中t0为发射电流关断时间;δ(t)为视深度，本次试验使用地空系统的电流关断时间为200 μs。

(9)

为了验证视电阻率计算的准确性和飞行高度对于视电阻率计算的影响，文中模拟了层状模型下(两层模型，第1层100 Ω·m，第2层1 Ω·m)3个飞行高度的时域地空电磁响应半解析解(图4(a))，并利用模拟出的电磁响应进行视电阻率计算(图4(b))，结果表明模型的电阻率分布与计算出的视电阻率基本一致，而在地空飞行探测中(<50 m)飞行高度对于视电阻率计算结果影响较小，因此可以在视电阻率计算中忽略飞行高度的影响。

5结果分析

2016年，山西省煤炭地质勘查研究院委托山西省煤炭地质物探测绘院联合吉林大学在沁水煤田晋城区块开展23 km2的时域地空电磁探测，圈定了该区域采空积水区域。根据已知地质资料,结合勘查区内的钻孔资料,在本区300测线附近绘制了该区地层岩性分布示意图(图5)，虽然该测线与地空电磁飞行测线并不完全重合，有约15°的夹角，但是基本反映了该区煤系地层的岩性分布特征。

部分地空测线的视电阻率断面如图6所示，本次探测的主要目标为山西组3号煤层，由视电阻率断面图可以看出，沿着3号煤层走向，出现多处较明显的低阻区域，如图6中3 750～4 950 m处、图7中0～1 800 m处和2 400～5 500 m处,低阻异常区的标高深度为200～300 m。结合之前文中对测区的地质描述可知，这些低阻区域与煤层走向及标高吻合度较高，并且与已掌握的地质资料情况一致，可以认定为煤炭釆空后形成的积水区。

图4 两层模型不同飞行高度条件下地空瞬变电磁场与视电阻率曲线特征
Fig.4 Characteristic of incidental electromagnetic field and apparent resistivity curves of two-layer model at different flight heights

图5 300线附近已知剖面地层岩性分布示意
Fig.5 Schematic map of lithologic distribution of known sections near line 300

图6 0300,3300线时间域地空电磁视电阻率断面
Fig.6 Time domain GREATEM apparent resistivity profile of line 0300 and 3300

6结 论

(1)在地面利用接地长导线源进行大功率发射、在空中利用无人机搭载电磁接收系统的工作方式，获得了较高信噪比的时域地空电磁数据。

(2)在数据预处理中，先后开展了基线校正和电磁噪声抑制方法研究，对采集到的电磁数据进行了进一步的优化处理，提升了数据质量便于进行下一步计算处理。

(3)研究了视电阻率及视深度计算方法，通过计算时域地空电磁视电阻率和视深度，获取研究区域地下电阻率分布特性。

(4)利用时域地空电磁探测系统准确地揭示该区域地下采空积水的情况，表明GREATEM系统对釆空积水的探测具有效性和实用性。

致谢感谢山西省煤炭地质物探测绘院和吉林大学为本文撰写提供的技术支持。

参考文献 :

[1] 秦荣军.我国煤炭开采现状及政策研究[J].煤炭经济研究,2019,39(1):57-61.

QIN Rongjun.Research on current status and policy of coal mining in China[J].Coal Economic Research,2019,39(1):57-61.

[2] MOGI T,TANAKA Y,KUSUNOKI K,et al.Development of grounded electrical source airborne transient EM(GREATEM)[J].Exploration Geophysics,1998(29):61-64.

[3] ITO H,MOGI T,JOMORI A,et al.Further investigations of underground resistivity structures in coastal areas using grounded-source airborne electromagnetics[J].Earth,Planets and Space,2011(63):e9-e12.

[4] ALLAH S A,MOGI T,ITO H,et al.Three-dimensional resistivity characterization of a coastal area:Application of Grounded Electrical-Source Airborne Transient Electromagnetic (GREATEM) survey data from Kujukuri Beach[J].Journal of Applied Geophysics,2013(99):1-11.

[5] 嵇艳鞠,王远,徐江，等.无人飞艇长导线源时域地空电磁勘探系统及应用[J].地球物理学报,2013,56(11):3640-3650.

JI Yanju,WANG Yuan,XU Jiang,et al.Development and applicatition of the grounded long wire source airborne electromagnetic exploration system based on an unmanned airship[J].Chinese Journal Geophys,2013,56(11):3640-3650.

[6] ALLAH S A,MOGI T,ITO H,et al.Three-dimensional resistivity modelling of grounded electrical-source airborne transient electromagnetic (GREATEM) survey data from the Nojima Fault,Awaji Island,south-east Japan[J].Exploration Geophysics,2014(45):49-61.

[7] LI Dongsheng,WANG Yuan,LIN Jun,et al.Electromagnetic noise reduction in grounded electrical source airborne transient electromagnetic signal using a stationary wavelet-based denoising algorithm[J].Near Surface Geophysics,2017,15,163-173.

[8] JI Yanju,YANG Guihong,GUAN Shanshan,et al.Interpretation research on electrical source of time domain ground-airborne electromagnetic data[J].IEEE Computer Society,2012,1-4.

[9] 嵇艳鞠,徐江,吴琼.基于神经网络电性源半航空视电阻率反演研究[J].电波科学学报,2014,29(5):937-980.

JI Yanju,XU Jiang,WU Qiong,et al.Apparent resistivity inversion of electrical source semi-airborne electromagnetic data based on neural network[J].Chinese Journal of Radio Science,2014,29(5):973-980.

[10] 李貅,张莹莹,卢绪山,等.电性源瞬变电磁地空逆合成孔径成像[J].地球物理学报,2015,58(1):277-288.

LI Xiu,ZHANG Yingying,LU Xushan,et al.Inverse Synthetic Aperture imaging of ground-airborne transient electromagnetic method with a galvanic source[J].Chinese Journal of Geophysics,2015,58(1):277-288.

[11] 王绪本,陈进超,郭全仕,等.沁水盆地北部煤层气富集区CSAMT勘探试验研究[J].地球物理学报,2013,56(12):4310-4323.

WANG Xuben,CHEN Jinchao,GUO Quanshi,el al.Research of the CSAMT exploration mode and experiment for the coalbed methane enrichment region in the north Qinshui Basin[J].Chinese Journal of Geophysics,2013,56(12):4310-4323.

[12] XUE Guoqiang,CHENG Jiulong,ZHOU Nannan,et al.Detection and monitoring of water-filled voids using transient electromagnetic method:A case study in Shanxi,China[J].Environmental Earth Sciences,2013,70,2263-2270.

[13] QIU Ning,CHANG Janjun,HE Zhanxiang.Comparison of several threshold selection rules of the wavelet denoising on the transient electromagnetic response[A].2ndWuhan International Conference on Environmental and Engineering Geophysics[C].Wuhan:2006.

[14] 李貅,宋建平,马宇,等.基于小波分析的TEM信号提取[J].煤田地质与勘探,2005,33(2):72-75.

LI Xiu,SONG Jianping,MA Yu,et al.The abstract of TEM signal based on the wavelet analysis[J].Coal Geology & Exploration,2005,33(2):72-75

[15] 李肃义,林君,阳贵红,等.电性源时域地空电磁数据小波去噪方法研究[J].地球物理学报,2013,56(9):3145-3152.

LI Suyi,LIN Jun,YANG Guihong,et al.Ground-airborne electromagnetic signals de-noising using a combined wavelet transform algorithm[J].Chinese Journal of Geophysics,2013,56(9):3145-3152.

[16] BOUCHEDDA A,CHOUTEAU M,KEATING P,et al.Sferics noise reduction in time-domain electromagnetic systems Application to MegaTEMⅡ signal enhancement[J].Exploration Geophysics,2010,41,225-239.

[17] JI Yanju,LI Dongsheng,YUAN Guiyang,et al.Noise reduction of time-domain electromagnetic data:Application of a combined wavelet de-noising method[J].Radio Science,2016,51(6):680-689.

[18] KAUFMAN A A,KELLER G V.Frequency and transient Soundings[M].London:Elsevier Applied Science Publishers,Ltd.,2014.

[19] WEIR G J.Transient electromagnetic fields about an infinitesimally long grounded horizontal electric dipole on the surface of a uniform half-space[J].Geophysics,1980,61,41-56.