郭建磊，高小伟，侯彦威

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

煤矿开采条件复杂，安全生产问题突出，其中，在掘进工作面发生重特大水灾事故占比达51.16%[1]，因此，提高超前预报的准确性对保障煤矿安全生产至关重要。煤矿超前探测中一般采用钻探和物探两种方法，其中，钻探方法结果可靠，但施工周期长、费用高，对巷道掘进生产影响较大。物探方法中，矿井瞬变电磁法[2]因无损、快速、信息丰富且对含水体反应灵敏等独特优点而被广泛应用，但受巷道复杂环境影响较大。巷-孔瞬变电磁超前探测方法[3]有效结合了矿井瞬变电磁法和超前钻孔的优势，该方法将发射线框布置在掘进工作面处，接收探头放入钻孔中，按间隔逐点采集三分量响应。孙怀凤等[4]通过物理模拟试验证明可以通过巷-孔瞬变电磁垂直分量信号在幅值和方向反号时间上存在的差异判断掘进工作面前方是否存在异常构造。陈丁等[5]采用积分方程法以煤层底板下方存在低阻体为模型研究了全空间条件下巷-孔瞬变电磁场的响应特征，结果表明，水平分量响应的延续时间优于垂向分量，垂向分量响应的幅值强于水平分量，在实际观测中不仅要观测垂向分量，也要观测水平分量。范涛[6-7]利用钻孔瞬变电磁垂直分量实现二维拟地震反演，并进行了三维数值模拟、水槽物理模拟和井下现场模拟试验；后续利用趋势面提取技术提取了水平分量纯异常场，通过异常场形态组合和幅值关系实现异常体中心方位角定位，并在此基础上通过聚类算法对水平分量异常形态进行自动分类，完成异常体空间角度定位[8]。该方法达到“一孔多用”的目的，有效克服了电磁和人文干扰，具有观测点距离异常体近、采集数据信噪比高等优点，在全国多个煤矿获得了广泛应用[9-10]。

巷-孔瞬变电磁法在资料处理解释过程中将地层假设为各向同性介质，但随着电磁勘探方法的发展，越来越多的学者认识到电导率各向异性在地下构造中普遍分布[11-13]，并对介质电导率各向异性特性与矢量电磁场关系做了大量研究。当前，关于频率域电磁法，如大地电磁法[14-18]、海洋可控源电磁法[19-21]、直流电法[22-24]、井间电磁法[25-26]和航空电磁法[27-28]等方法的各向异性正反演研究均取得了较为丰富的成果。关于时间域电磁法的各向异性研究主要针对半空间介质和电导率各向异性对垂直分量的影响。周建美等[29]基于有限体积法实现双轴各向异性瞬变电磁三维正演，刘亚军等[30]基于有限体积法实现任意各向异性瞬变电磁三维正演。在全空间介质方面，程久龙等[31]基于时域有限差分法研究了不同主轴各向异性介质对矿井瞬变电磁场垂直分量的影响特征，解释了探测结果存在偏差的原因，并结合实例进行说明。上述研究并未涉及电导率各向异性对三分量响应的影响。

目前，巷-孔瞬变电磁法从硬件到软件均实现了三分量信号的采集与解释，通过将地层假设为各向同性介质，依据水平分量形态组合及幅值关系判定孔旁异常体的分布象限，依据垂直分量反演计算异常体相关参数，最终实现巷-孔瞬变电磁立体成像解释。由于资料解释受制于对介质电导率各向同性的假设，导致在存在明显电导率各向异性的勘探区会产生较大的解释误差，因此，非常有必要研究全空间条件下电导率各向异性对巷-孔瞬变电磁响应的影响方式及程度，为进一步提高超前精度提供指导依据。本文研究采用时域有限差分算法(Finite Difference Time Domain，FDTD)[32]实现各向异性瞬变电磁三维正演，通过引入电导率各向异性张量构建控制方程，以差分代替微分对控制方程进行离散，空间离散采用后向差分，时间离散采用中心差分，以电流密度的形式加入各向异性Maxwell 方程组的安培环路定理实现任意电流源的加载，进而推导出有源和无源区域电磁场迭代表达式。最后，构建全空间模型、层状模型和三维块状模型进行正演并分析电导率轴向各向异性对巷-孔瞬变电磁三分量响应的影响程度与方式。

1 正演理论基础

1.1 控制方程

电导率各向异性、无源媒质中的Maxwell 方程组为：

式中：E为电场强度，V/m；H为磁场强度，A/m；为张量电导率，S/m；μ为磁导率，Wb/(A·m)；ε为介电常数，F/m；t为时间，s。

地球物理低频电磁勘探中一般忽略位移电流，为了满足三维时域有限差分计算要求，因此，加入虚拟介电常数构成显式的时间迭代格式，式(1b)变为：

式中：γ为虚拟介电常数。

采用经典Yee 氏网格对计算区域进行剖分，电场分量位于棱边中心点，磁场分量位于面中心点，每一个电场(磁场)分量均被相应的4 个磁场(电场)分量包围，电场分量和磁场分量在时间轴上交替采样。

采用非均匀网格形式[32]进行模型剖分，将式(2)在t(n+1/2)时刻进行空间和时间离散，以差分代替微分，其中，空间上采用后向差分、时间上采用中心差分，无源区域得到如下轴向各向异性的电场强度迭代公式：

在有源媒质区域，式(2)必须包含源电流密度项，将式(2)修改为：

式中：Js为源电流密度。

巷-孔瞬变电磁法采用回线源进行发射，坐标系满足“右手定则”，假定，巷道方向为z方向，上下为x方向，左右为y方向，发射源将被放置在xoy平面，因此，式(6)中仅存在x、y方向上的源电流密度项(Jsx、Jsy)。考虑低频近似后，在直角坐标系中将式(6)展开得到下式，按照差分格式正常离散，可以得到有源区域网格电场迭代公式。磁场迭代公式不包含电导率项，因此，各向异性条件下磁场迭代公式与各向同性条件下磁场迭代公式相同[32]。

1.2 算法验证

构建均匀全空间各向同性模型进行正演模拟，通过将数值解与解析解进行对比验证算法精度。全空间介质电导率为0.01 S/m，采用笛卡尔坐标系，回线源中心位于坐标系原点，边长为3 m×3 m，电流为1 A，发射源上升沿和下降沿均为1 μs，持续时间为10 ms，二次场采样时间为10 ms。最小网格尺寸为1 m，相邻网格放大系数为1.05，剖分网格数为221×221×200，接收点位于(0.5 m，0.5 m，0 m)，在接收点处接收∂Bx/∂t、∂By/∂t、∂Bz/∂t响应。

数值解与解析解对比结果如图1 所示。由图1 可知，数值解和解析解的三分量响应曲线基本重合，当时间达到0.003 ms 时误差控制在5%以内，其中，∂Bx/∂t和∂By/∂t的相对误差稍大于∂Bz/∂t的相对误差，由于接收点位于坐标系对角线上，故∂Bx/∂t与∂By/∂t响应曲线重合且相对误差一致。故本文算法满足计算精度要求。

图1 数值解与解析解对比结果Fig.1 Comparison of numerical solutions and analytical solutions

2 各向异性巷-孔瞬变电磁响应特征

2.1 各向异性全空间模型

巷-孔瞬变电磁法勘探环境为全空间介质，因此，建立如图2 所示模型研究全空间介质电导率各向异性对巷-孔瞬变电磁三分量响应的影响方式与程度。如图2 所示，采用上节的坐标系系统，巷道内空气电导率为0.000 1 S/m，空腔范围为(-3～3 m,-3～3 m,0～-∞ m)。发射源平贴放置于巷道工作面且中心位于坐标原点，测点位于超前钻孔中，钻孔孔口位于(0.5 m,0.5 m,0 m)，钻孔平行于z轴，长度范围为0～80 m，按照1 m 的点间距从孔底至孔口逐点采集三分量响应数据。模型a为各向同性，电导率为σx=σy=σz=0.01 S/m；模型b 为x轴各向异性，电导率为σx=0.1 S/m、σy=σz=0.01 S/m；模型c 为y轴各向异性，电导率为σy=0.1 S/m、σx=σz=0.01 S/m；模型d 为z轴各向异性，电导率为σz=0.1 S/m、σx=σy=0.01 S/m。模型剖分与发射源参数与1.2 节一致。

图2 全空间各向异性巷-孔瞬变电磁超前探测模型Fig.2 Full-space anisotropic tunnel-hole transient electromagnetic advance detection model

图3 为4 个模型的三分量响应多测道图，时间道的时间为0.06～0.21 ms。由图3 可知，4 个模型的∂Bx/∂t、∂By/∂t多测道图形态呈反“C”型、∂Bz/∂t多测道图形态呈半“C”型；模型a 的∂Bx/∂t、∂By/∂t多测道图形态和幅值一致，最大幅值位于25 m；∂Bz/∂t多测道图最大幅值位于0 m，随深度增加逐渐衰减；与模型a 相比，模型b 的∂Bx/∂t多测道图形态和幅值在发射源处发生微弱变化，∂By/∂t多测道图幅值整体增大且最大幅值位于10 m、至40 m 趋于0，∂Bz/∂t多测道图幅值整体增大，至40 m 趋于0；与模型a 相比，模型c 的∂Bx/∂t多测道图幅值整体增大且最大幅值位于10 m、至40 m 趋于0，∂By/∂t多测道图形态和幅值在发射源处发生微弱变化；∂Bz/∂t多测道图幅值整体增大、至40 m处趋于0。模型d 的三分量响应多测道图特征与模型a 的基本一致。

图3 三分量响应多测道图Fig.3 Three-component responses multi-channel diagram

将4 个模型的0.06 ms 时刻的三分量响应曲线进行对比(图4)发现，x轴各向异性主要影响∂By/∂t和∂Bz/∂t响应，y轴各向异性主要影响∂Bx/∂t和∂Bz/∂t响应，z轴各向异性对三分量响应基本不产生影响。根据“烟圈效应”，回线源产生的感应电流主要在水平方向，三分量响应受到水平方向电导率的影响较大，而受到垂直方向的电导率影响很小；x方向感应电流垂直穿过y轴电导率，y方向感应电流垂直穿过x轴电导率，因此，x轴电导率主要影响∂Bx/∂t响应，y轴电导率主要影响∂By/∂t响应。

图4 三分量响应曲线对比(t=0.06 ms)Fig.4 Comparison of three-component responses (t=0.06 ms)

2.2 各向异性层状模型

对于水平钻孔而言，结合煤炭固有的地层特征，巷-孔瞬变电磁法勘探环境基本处于层状介质。当回线源位置固定后，地下不同深度地层受到一次场激发的程度不同，地层电导率各向异性对三分量响应的影响程度也不一样。

以各向同性全空间模型为基础，分别将距回线源不同距离的地层设置为各向异性(图5)。如图5 所示，模型a 为上部地层(40～∞ m)电导率分别为各向同性和各向异性、模型b 为中间地层(-40～40 m)电导率分别为各向同性和各向异性、模型c 为下部地层(-40～-∞ m)电导率分别为各向同性和各向异性。当地层为各向同性时，电导率σx=σy=σz=0.01 S/m；当x轴为各向异性时，电导率σx=0.1 S/m、σy=σz=0.01 S/m；当y轴为各向异性时，电导率σy=0.1 S/m、σx=σz=0.01 S/m；当z轴为各向异性时，电导率σz=0.1 S/m、σx=σy=0.01 S/m。模型剖分与发射源参数与1.2 节一致。

图5 层状模型Fig.5 Schematic diagram of the layered models

对比分析模型a-模型c 的不同电导率参数下0.06 ms 时刻的三分量响应曲线(图6)。由图6a 可知，对于回线源上部的地层，y轴各向异性对∂Bx/∂t响应影响明显，且使∂Bx/∂t幅值由负变正，对∂By/∂t、∂Bz/∂t响应几乎没有影响；同时，x轴、z轴各向异性与各向同性的三分量响应曲线几乎重合。由图6c 可知，y轴各向异性对∂Bx/∂t响应影响明显，与模型a 的∂Bx/∂t响应反映不同的是，模型c 的∂Bx/∂t响应幅值正负不变且幅值增大。因此，通过水平分量形态和幅值的相互关系可以判断各向异性介质所在方位及电导率主轴方向，这主要是因为模型a 的地层位于回线源上方，模型c的地层位于回线源下方，巷-孔瞬变电磁三分量响应具有方向性。由图6b 可知，x轴、z轴各向异性和各向同性的∂Bx/∂t响应曲线重合且偏离y轴各向异性的∂Bx/∂t响应曲线；y轴、z轴各向异性和各向同性的∂By/∂t响应曲线重合且偏离x轴各向异性的∂By/∂t响应曲线；x轴和y轴各向异性的∂Bz/∂t响应曲线重合且与z轴各向异性和各向同性的∂Bz/∂t响应曲线偏离较大，说明距离发射源较近地层的电导率各向异性对三分量响应的影响较大，也再次说明巷-孔瞬变电磁三分量响应主要受水平方向电导率的影响。

图6 三分量响应曲线对比(t=0.06 ms)Fig.6 Comparison of three-component responses (t=0.06 ms)

2.3 各向异性块状模型

煤炭资源广泛分布于沉积岩中，沉积岩地区由于层理发育导致地下介质电阻率随电流方向发生变化[33]；煤炭开采产生不同形态、高度和密度的采空裂隙和裂缝，裂隙和裂缝充水后形成裂隙含水体[34]，上述情况均会导致地层介质具有电导率各向异性特性。因此，本节主要研究各向同性地层中含有各向异性块状模型和各向异性地层中含有各向同性块体模型。

1)各向同性地层中含有各向异性块状模型

建立如图7 所示的模型，异常体1 的规模为20 m×20 m×20 m，异常体2 的规模为30 m×20 m×20 m，异常体3 的规模为20 m×30 m×20 m，3 个异常体的中心均位于(20 m，20 m，20 m)。将3 个异常体分别设置为各向同性和各向异性，为各向同性时电导率为σx=σy=σz=0.1 S/m；为x轴各向 异性时电导率 为σx=1 S/m、σy=σz=0.1 S/m；为y轴各向异性时电导率为σy=1 S/m、σx=σz=0.1 S/m；为z轴各向异性时电导率为σz=1 S/m、σx=σy=0.1 S/m。模型剖分与发射源参数与1.2 节一致。

图7 块状异常体模型俯视图Fig.7 Top view of block abnormal model

将3 个异常体的0.06 ms 时刻的三分量响应曲线进行对比(图8)。由图8 可知，当3 个异常体均为各向同性时，∂Bx/∂t响应大小为异常体3＞异常体2＞异常体1，说明y方向边长改变影响∂Bx/∂t响应；∂By/∂t响应的大小为异常体2＞异常体3＞异常体1，说明改变x方向边长改变影响∂By/∂t响应；∂Bz/∂t响应大小为异常体2=异常体3 ＞异常体1。当3 个异常体均为z轴各向异性时，其三分量响应形态、规律与各向同性时相同，但其水平分量响应幅度大于各向同性；当异常体均为x轴和y轴各向异性时，其∂Bx/∂t和∂By/∂t响应形态、规律与各向同性的相同，不同之处在于，当3 个异常体均为x轴各向异性时，异常体2 的∂Bz/∂t响应＞异常体3，当3 个异常体均为y轴各向异性时，异常体3 的∂Bz/∂t响应＞异常体2。

图8 三分量响应曲线对比(t=0.06 ms)Fig.8 Comparison of three-component responses (t=0.06 ms)

将水平分量响应曲线进行进一步对比(图9)。当3 个异常体均为x轴各向异性时，异常体1 的∂By/∂t响应＞∂Bx/∂t响应，异常体3 的∂Bx/∂t(∂By/∂t)响应＞异常体2 的∂By/∂t(∂Bx/∂t)响应。当3 个异常体均为y轴各向异性时，异常体1 的∂Bx/∂t响应＞∂By/∂t响应，异常体2 的∂By/∂t(∂Bx/∂t)响应＞异常体3 的∂Bx/∂t(∂By/∂t)响应。当3 个异常体均为z轴各向异性时，异常体1 的∂Bx/∂t响应和∂By/∂t响应曲线重合，异常体2 的∂By/∂t(∂Bx/∂t)响应与异常体3 的∂Bx/∂t(∂By/∂t)响应曲线重合，且异常体2 的∂By/∂t(异常体3 的∂Bx/∂t)响应＞异常体2 的∂Bx/∂t(异常体3 的∂By/∂t)响应。上述情况说明异常体边长不同对水平分量响应的影响强弱不同，与模型长轴平行方向的电导率改变对三分量响应的影响较大。

图9 水平分量响应对比(t=0.06 ms)Fig.9 Comparison diagram of horizontal component responses (t=0.06 ms)

2)各向异性地层中含有各向同性块体模型

基于图7 中的异常体1 模型，将异常体1 设置为各向同性，电导率为σx=σy=σz=0.1 S/m，将全空间介质分别设置为各向同性和各向异性，为各向同性时电导率为σx=σy=σz=0.01 S/m；为x轴各向异性时电导率为σx=0.1 S/m、σy=σz=0.01 S/m；为y轴各向异性时电导率为σy=0.1 S/m、σx=σz=0.01 S/m；为z轴各向异性时电导率为σz=0.1 S/m、σx=σy=0.01 S/m。

将全空间介质为各向同性和各向异性时0.06 ms时刻的三分量响应曲线进行对比(图10)。由图10 可知，全空间介质为各向同性与z轴各向异性的三分量响应几乎重合，∂Bx/∂t和∂By/∂t响应曲线在距离发射源40 m 后出现偏离；改变全空间介质x轴和y轴电导率对异常体的三分量响应均影响较大，使∂Bx/∂t和∂By/∂t异常幅值由正变负，且改变y轴电导率对∂Bx/∂t响应的影响＞改变x轴电导率、改变x轴电导率对∂By/∂t响应的影响＞改变y轴电导率，改变x轴和y轴电导率对∂Bz/∂t响应的影响一样，这与上述小节的规律一致。全空间介质为各向异性的三分量异常响应大于各向同性的三分量异常响应，说明异常体产生的异常响应被淹没在全空间介质电导率各向异性产生的异常响应中，如果此时不考虑介质的各向异性特性会对巷-孔瞬变电磁解释造成较大误差。

图10 三分量响应对比(t=0.06 ms)Fig.10 Comparison of three-component responses (t=0.06 ms)

3 结论

a.介质距回线源的距离不同，其各向异性特性对巷-孔瞬变电磁三分量响应的影响强弱和方式不同，可以通过其三分量响应的形态和幅值的相互关系辨别各向异性介质所在方位及主轴电导率方向，同时由于三维异常体的边长不同对其三分量响应的影响强弱不同，进一步增加了各向异性条件下三分量解释的复杂性。

b.回线源产生的感应电流主要是水平方向，并且x方向感应电流垂直穿过y轴电导率，y方向感应电流垂直穿过x轴电导率，因此，z轴各向异性对三分量响应基本没有影响，x轴和y轴各向异性对巷-孔瞬变电磁三分量响应有较大影响，同时∂Bx/∂t响应主要受y轴电导率影响，∂By/∂t响应主要受x轴电导率影响。

c.当异常体所处介质呈电导率各向异性时，异常体产生的异常响应被淹没在全空间介质电导率各向异性产生的异常响应中，介质的电导率各向异性特征在巷-孔瞬变电磁法解释过程中不可忽略，因此，本文研究内容为进一步研究各向异性提供一定的基础，也为后续进行各向异性反演提供参考。

d.未来的工作包括引入更复杂的电导率张量、进行亚网格剖分、添加吸收边界条件等，实现任意各向异性的不规则复杂模型三维快速正演，为实际生产等提供重要的参考意义。