易洪春， 袁永榜， 牟义， 常江浩

(1.瓦斯灾害监控与应急技术国家重点实验室， 重庆 400039； 2.中煤科工集团重庆研究院有限公司， 重庆 400037； 3.煤炭科学技术研究院有限公司安全分院， 北京 100013； 4.河北地质大学地球科学学院， 石家庄 050031)

陷落柱是一种埋藏型岩溶在动力作用下发生坍塌和剥落后形成的锥状塌陷体。陷落柱除自身常储聚大量地下水外，还常成为沟通其他水源的导水通道，故陷落柱的存在容易造成矿井突水事故。突水性强、规模大、破坏严重一直是陷落柱突水的突出特点。

瞬变电磁法(transient electromagnetic method，TEM)是利用不接地回线或接地线源向地下发射一次脉冲磁场，在一次脉冲磁场间歇期间利用线圈或接地电极观测地下介质中引起的二次感应涡流场，从而探测介质电阻率的一种方法，近年来被广泛应用到煤矿井下水害探测中，取得了大量精准有效的探测成果。而瞬变电磁场在含水陷落柱中的传播机理及响应特征一直是利用瞬变电磁法探测含水陷落柱需要弄清的问题。时域有限差分(finite-difference time-domain，FDTD)算法已被许多学者借鉴到电磁场的数值计算问题上来[1-3]，它不仅可以研究瞬变电磁场的产生机理，扩散过程及响应特征，同时也为瞬变电磁法探测含水低阻体提供了大量理论支撑。岳建华等[4]研究了层状介质中地下瞬变电磁场的全空间效应；杨海燕等[5]对矿用多匝小回线源瞬变电磁场的分布规律进行了详细分析；李建慧等[6]对基于电场方程的回线源产生的瞬变电磁场进行了三维正演研究；孟庆鑫等[7]对井中磁源瞬变电磁场进行了三维时域有限差分数值模拟；席振铢等[8]首次对瞬变电磁法三分量测量方法进行了研究；常江浩等[9]对矿井全空间磁源瞬变电磁场的扩散规律进行了详细研究；程久龙等[10]对矿井全空间三维主轴各向异性介质瞬变电磁场响应特征进行了深入研究，解释了探测结果存在偏差的原因。

目前，关于瞬变电磁场的研究主要集中于三分量中的某一分量，而对某两个或三个分量的合矢量的研究相对较少。鉴于此，在前人研究的基础之上，现采用基于非均匀网格的时域有限差分算法，研究全空间瞬变电磁场的扩散，并首次基于电场合矢量，详细分析电场的传播规律及响应特征，为瞬变电磁法超前探测掘进头前方的含水区提供理论支撑，并为用合适量来进行瞬变电磁数据处理解释奠定基础。最后通过数值模型算例以及实际应用验证瞬变电磁法探测含水陷落柱的有效性和适用性。

1 瞬变电磁场的FDTD算法

1.1 FDTD方程

时域有限差分法(FDTD)是求解三维时域电磁场问题的主要方法之一，其不需要求解大型方程，因此求解速度快且能够较好地实行并行运算。基于此，现采用时域有限差分法来进行瞬变电磁场三维模拟计算。FDTD算法是一种直接基于时域电磁场方程的数值算法，算法的核心思想是将电磁场的时域微分方程的求解转化为差分方程的迭代求解，近年来被广泛用于瞬变电磁场的数值计算。经常用到的Yee网格其磁场分量位于面的中心，而电场分量位于棱的中间[11-12]。 在计算瞬变电磁场时，使用准静态条件下的麦克斯韦方程[13]：

(1)

(2)

(3)

(4)

式中：H(r，t)为磁场强度，A/m；E(r，t)为电场强度，V/m；σ(r)为介质电导率，S/m；μ为介质磁导率，H/m；γ为位移介电常数；r为电磁波源点到场点的位置矢量；t为电磁波传播时间。电磁场的差分方程可以通过对微分方程[式(1)～式(4)]的离散来获得[14-15]。

1.2 源的加载

全空间条件下磁偶源中电流被突然切断所产生的瞬变磁场[16-17]为

(5)

1.3 稳定条件

电磁场时域有限差分法采用显式差分格式，其空间网格大小和时间步长需要满足一定的条件才能保证迭代的稳定。电磁场时域有限差分迭代的数值稳定性由Courant条件确定[18-19]。

(6)

式(6)中：ν为电磁波传播速度；Δx、Δy、Δz分别为x、y、z方向的空间最小网格间距。

1.4 边界条件

由于计算机模拟有限空间，截断边界会对计算结果产生影响，需采用合适的边界条件进行消除，本文研究中采用适合瞬变电磁场计算的卷积完全匹配层(CPML)边界进行求解。

2 低阻陷落柱影响下矿井全空间瞬变电磁场的扩散及响应特征

2.1 低阻陷落柱影响下的全空间瞬变电磁场的扩散

陷落柱是煤田广泛存在的一种地质构造，其导水性强，易于充水，已成为矿井水害的元凶之一。建立煤系地层中含水陷落柱的井下全空间地球物理模型，模拟陷落柱对瞬变电磁的响应特征及规律，把控瞬变电磁探测陷落柱的机理，可有效预防陷落柱突水事故的发生[20-21]。

图1 陷落柱模型Fig.1 Collapse column model

图1为煤系地层含水陷落柱的井下全空间地球物理模型。含水陷落柱发育于煤层底板以下，其顶部距煤层20 m，底部距离煤层70 m，长20 m、宽20 m，顶部中心坐标为(60，60，-20)，电阻率设为10 Ω·m，巷道位于煤层中，其截面大小为6 m×6 m。煤层厚10 m，其电阻率及围岩的电阻率为100 Ω·m。发射回线大小为2 m×2 m，位于煤层中，回线的中心为坐标系原点，其法向沿x轴方向，满足右手坐标系，发射电流为1 A。

图2为柱状低阻陷落柱影响下的不同时刻x方向感应磁场分量Hx在xy平面(z=-50 m)的等值线平面图，单位A/m。如图2所示，磁场的最大值始终保持在xy平面的中心位置(z轴上)。随着感应磁场从发射回线附近不断向周围空间扩散，其扩散范围到达陷落柱位置时，感应磁场的分布受到陷落柱的影响。

图3为柱状低阻陷落柱影响下的不同时刻y方向感应电场分量Ey在yz平面(x=60 m)的等值线剖面图，单位V/m。如图3所示，感应电场的正负最大值分别位于z=0平面的上下两侧，随着感应电场向上下空间不断地扩散，感应电场的正负最大值也向上下空间移动，当电场扩散到陷落柱体内，电场分布受到陷落柱的显著影响，而电场等值线畸变的范围恰好与陷落柱的位置和深度吻合(图3中虚线框所示)。当t=500 μs时，感应电场正负最大值已经扩散到周围较远空间，此时电场等值线近似水平分布。

2.2 基于合矢量的瞬变电磁场的响应特征

某一方向的场强(电场或磁场)分量的扩散过程，一定程度上反映了介质及异常体的特性，但很难体现出异常体的整体形态，而实际上，场强是x、y、z三个方向的合矢量，用合矢量来描述电磁场的扩散过程，更能体现低阻异常体的电性、形态、位置等响应特征。下面以电场为例，进行电磁场扩散特性的描述。

图4 感应电场分量Ey、Ez及其合矢量Ey+z的等值线分布图Fig.4 Contour distribution of induced electric field components Ey，Ez and their resultant vector Ey+z

图5 感应电场分量Ex、Ey及其合矢量Ex+y的等值线分布图Fig.5 Contour distribution of induced electric field components Ex，Ey and their resultant vector Ex+y

图6为不同深度(从上至下分别为z=0 m、z=-50 m和z=-100 m)的x和y方向的电场合矢量Ex+y的在xy平面的等值线平面图，图7为不同横向位置(从左到右分别为x=40 m、x=60 m和x=80 m)的y和z方向电场合矢量Ey+z在yz剖面的等值线剖面图。如图6所示，在初始阶段，感应电场主要分布在z=0平面附近，其分布只受发射回线周围介质电阻率的影响。随着时间的推移，感应电场的扩散范围慢慢扩大，其分布受到低阻陷落柱的影响。在t=500 μs时，感应电场已经进入陷落柱体内，它的等值线分布受到陷落柱的显著影响，等值线在陷落柱内发生畸变。在不同深度的畸变的范围与陷落柱模型的位置十分吻合，如图6(b)所示，特别是在z=-50 m的平面上，能清晰地看到等值线畸变的范围恰好是陷落柱在z=-50 m平面所在的位置。在不同横向位置的yz剖面的电场合矢量Ey+z的等值线分布的畸变区域的位置和高度也恰好与陷落柱模型的位置和深度十分吻合，如图7所示。当陷落柱完全包含于磁场畸变的中心时，此时的陷落柱已经成为一个“二次源”，向空间辐射电磁场。此时，观测到的信号主要反映了瞬变电磁场在低阻体影响下的衰减情况。

通过对电磁场的扩散过程进行数值模拟，表明瞬变电磁场的分布会受到低阻体的显著影响，并且可以通过电磁场的变化来确定低阻体的形态和位置。因此，可以通过观测瞬变电磁场的变化来探测低阻异常体。

图6 不同深度的x和y方向的电场合矢量Ex+y的等值线平面图Fig.6 Plan views of contours of the vector Ex+y in the electrical field in the x and y directions at different depths

3 应用实例

基于前文研究的理论基础，将瞬变电磁法应用到煤矿含水陷落柱的实际探测中，通过应用实例，进一步验证了瞬变电磁法探测低阻陷落柱的有效性。

3.1 实例1

山西官地矿23511工作面掘进巷道沿3#煤层布置，距2#煤底板平均5 m，平均煤厚3.0 m，煤层倾角平均4°。顶板为泥岩、细砂岩及砂质泥岩互层，底板为泥岩。工作面整体为褶皱构造，陷落柱较发育。据2#煤采掘资料，该巷道预计揭露正断层4条，其中断距大于1 m的2条，均不存在导水现象；工作面内存在1个陷落柱，积水情况不明。为保障巷道的安全掘进，需提前查明前方100 m范围内的陷落柱富含水情况。

图7 不同横向位置的y和z方向的电场合矢量Ey+z的等值线剖面图Fig.7 Section views of contours of the vector Ey+z in the electric field in the y and z directions at different positions

为了能够准较直观地反映巷道掘进头前方低阻异常体的方位和距离，常使用扇形探测技术[18]来描述前方介质的地电特性，即在视电阻率等值线图上将每个探测方向上的视电阻率值按照时深转换结果绘制到该方向对应的位置上，形成二维的扇形图。其中由感应电动势求取视电阻率参考文献[19]，并通过时深转换来获得对应深度[20]。

如图8所示，根据探测任务，每隔15°进行一个方向探测，共13个方向，形成从左帮到右帮180°范围的扇形扫描。扫描结果经过视电阻率计算，时深转换和坐标换算，形成的二维扇形视电阻率等值线拟断面图，如图9所示。据图9可知，在探测位置前方30～90 m、左帮10～80 m范围视电阻率值相对正前方及右帮区域明显偏低，且为明显的视电阻率等值线低阻“圈闭”异常，结合水文地质资料解释此区域为陷落柱含水造成的低阻异常。经矿方钻探验证为陷落柱含水水区，出水点位于掘进前方33 m处顶板上方4.5 m，积水范围自巷道至左侧帮75 m，出水量约20 m3/h，采取防治水措施后安全掘过了此区域。

图8 瞬变电磁超前探测扫描方向示意图Fig.8 Schematic diagram of the scanning direction of transient electromagnetic advance detection

图9 瞬变电磁超前探测视电阻率等值线拟断面图Fig.9 Pseudo-section diagram of the apparent resistivity contour of the transient electromagnetic advanced detection

3.2 实例2

山西晋城矿区王台铺矿地质条件较复杂，为保障巷道的安全掘进，需要准确查明前方采空区及陷落柱含水情况。探测方式除同实例1中水平方向扫描巷道迎头左帮到右帮180°范围外(每15°一个方向，共13个方向)，如图10(a)所示，还在剖面方向探测顶板向上+90°、+60°、+30°、水平0°、底板向下-30°、-60°、-90°共7个方向，如图10(b)所示。

图10 探测方式布置示意图Fig.10 Schematic diagram of detection mode layout

探测结果经过视电阻率计算，时深转换和坐标转换，并利用Voxler成像软件形成三维超前探测成果图。图中坐标原点为探测位置，+x方向沿巷道右帮方向，+y方向沿巷道掘进方向，+z方向为巷顶底板方向，形成的三维数据体立体图和三维等值面图分别如图11和图12所示，图11和图12中从蓝色到红色(冷色调到暖色调)表示视电阻率不断升高。可以看出图中存在3处较大的异常范围，其中在迎头正前方存在1处高阻范围，推测为不含水区；在迎头两帮各存在1处低阻范围；推断为含水陷落柱，而在巷道两帮都存在低阻响应推测为瞬变电磁的全空间效应引起的。

后经过矿方打钻验证，正前方50 m之内的高阻异常区，为无水区；结合其他地质资料，向右侧帮的低阻异常区打钻，于+x方向35 m处出水，水量较大，揭露为含水陷落柱。探测结果说明，瞬变电磁法对探测含水陷落柱十分有效。

图11 三维数据体立体图图Fig.11 Three-dimensional view of a 3D data body

图12 三维数据体等值面图Fig.12 Isosurface map of 3D data body

4 结论

通过全空间三维时域有限差分数值算法，模拟了在低阻陷落柱影响下的全空间瞬变电磁场的扩散过程，并对电磁场某一方向上的分量及其合矢量进行了详细分析，最后将结果应用于实际探测。得出如下结论。

(1)FDTD法是模拟瞬变电磁场特征规律的有效方法，该方法在求解瞬变电磁场微分方程时将时空都离散化，把空间划分为网格，时间上按步进的方式迭代，在计算过程中的每一时间步都对全场所有的网格进行一次计算，因而能够展示出场随时间的波动在空间的传播过程，并给出清晰的物理图像。这种模拟过程可以为瞬变电磁法探测低阻陷落柱提供理论支撑。

(2)低阻体会严重影响瞬变电磁场的扩散过程。某一方向的场强分量的扩散过程，一定程度上反映了介质及异常体的特性，但很难体现出异常体的整体形态，而用合矢量来描述电磁场的扩散过程，更能体现异常体的电性、形态、位置等特征。

(3)多方向探测方法为掘进头前方的目标定位提供了一种解决方案，该方法不仅可以检测到低阻体，还可以确定其方位和距离。数值模拟和实际应用都证明了瞬变电磁法探测含水陷落柱的有效性。

(4)对于瞬变电磁的分支——钻孔瞬变电磁或者地井瞬变电磁，目前常规的分析方法是针对不同单一分量进行曲线分析，确定异常方位和形状，该方法分析起来十分繁琐，而用合矢量来分析异常响应特征可能会找到一种简便的数据处理解释方法，这需要进一步研究。