随着中国煤矿开采深度及开采强度不断加大，煤层气勘探工作不断深入，采空区、采动裂隙、煤岩层中的煤层气资源，已成为煤层气富集矿井的重要资源^[1]。煤层气主要是以吸附状态赋存于煤层微孔隙中，由于不同的沉积环境、煤层埋藏和热演化史、构造作用、温压条件和水文条件，导致煤层气在煤层中的吸附解吸动态平衡发生变化，使煤层气分布不均匀，进而形成煤层气的相对富集区^[2-4]。煤矿采动过程会加速破坏解吸动态平衡，使煤层气加速向采空区和采动裂隙中扩展^[5-6]，因此，精准探测采空区及采动裂隙范围，是矿井采空区煤层气资源开发的前提。

在煤层气富集矿井中，采空区及采动裂隙中一般存在一定量的积水，形成煤层气与水共存的地质情况。采空区内的积水与水基压裂液作为低阻目标体，与煤层气位于同一采空范围之内，因此可通过对其进行探测确定采空区及采动裂隙发育范围，进而确定煤层气富集区。在众多积水采空区探测方法中，矿井瞬变电磁法对低阻体反应敏感、多匝小回线装置轻便高效，成为最常用的探测富水采空区及采动裂隙方法之一^[7]，也成为间接探测煤层气的有效方法。在煤层气富集矿井瞬变电磁法探测方面，葛燕燕等^[8]以潘庄区块为研究对象，利用瞬变电磁法对排采井附近的地下水进行探测，分析煤层气排采井地下水的运移规律，对煤层气排采期间水文地质动态变化规律进行监测。范涛等^[9]提取了瞬变电磁信号中的波场特性，基于速度分析反演成像方法，通过对比压裂前后的数据对煤层气水力压裂的效果进行了评价。段建华等^[10]采用瞬变电磁法和微震监测技术联合对井下煤层气水力压裂效果及影响范围进行监测。

综合以上研究成果可以看出，目前利用瞬变电磁法进行煤层气富集矿井低阻目标体的探测多为基于传统数据处理与解释方法，受体积效应影响，视电阻率断面探测成果可反映低阻区域，但对于最为重要的采空区及采动裂隙边界位置则难以区分，且基于井下全空间瞬变电磁法的含水采空区及采动裂隙电性界面探测技术研究较少，未能实现采空区及采动裂隙范围的精准探测。

笔者从波场变换理论积分形式入手，采用非线性阻尼最小二乘法求解波场变换积分方程，通过Fletcher策略进行迭代计算，针对井下探测进行全空间响应校正与波场数据处理，重点对采空区探测波场合成孔径成像方法进行了阐述，实现了矿井瞬变电磁扩散场到波动场的稳定变换计算和数据处理，开展了数值模拟及井下超前与工作面探测应用试验，对波场变换成像方法的有效性进行了综合分析及验证。

1 全空间波场变换 1.1 波场变换原理

在瞬变电磁波场变换研究方向，Weidelt^[11]、Kunetz^[12]、Levy等^[13]发现了地震波场和大地电磁场之间的内在关系，Lee等^[14-15]基于麦克斯韦方程组推导了瞬变电磁扩散场与对应波场之间的数学关系，并证明了这种变换关系适用于任意源和矢量电磁场，瞬变电磁场分量与其相对应的波场之间的数学积分形式为

$ \boldsymbol{P}(r, t)=\frac{1}{2 \sqrt{{\rm{ \mathsf{ π} }} t^{3}}} \int\limits_{0}^{\infty} q \;\; \mathrm{e}^{-\frac{q^{2}}{4 t}} \boldsymbol{U}(r, q) \mathrm{d} q, $

(1)

式中：P(r, t)为瞬变电磁场分量，U(r, q)为波场分量，t为瞬变电磁场时间，q为波场类时间变量，以s^1/2为单位，r为各场分量对应时间的场值，e为自然对数的底数。通过计算该积分方程，可以将瞬变电磁场随时间t的响应转换为波场随类时间q的响应，便于利用地震数据处理方法对转换波场的物性与几何参数进行求解，提取原瞬变电磁扩散场数据中包含的电性界面信息，有利于降低方法的体积效应，提高分辨率。

将式(1)写成矩阵形式为

$ \boldsymbol{A} \boldsymbol{U}=\boldsymbol{P}。$

(2)

利用梯形公式对式(1)的积分方程进行离散，得到式(3)：

$ \boldsymbol{U}=\left(u_{1}, \cdots, u_{j}, \cdots, u_{n}\right)^{\mathrm{T}}, \boldsymbol{P}=\left(p_{1}, \cdots, p_{i}, \cdots, p_{m}\right)^{\mathrm{T}} \quad(i=1,2, \cdots, m, j=1,2, \cdots, n)。$

(3)

式(2)中，A为m×n阶矩阵，U为转换波场数据，P为离散化的瞬变电磁场数据。求解式(2)，即可得到变换后的波场响应，矩阵方程式(2)的求解为不适定反问题，可采用正则化方法求解^[16]。波场变换计算的关键在于正确求解式(2)，本次采用阻尼最小二乘法求解该不适定反问题，并利用Fletcher策略进行迭代计算^[17]，根据目标函数的非线性程度对阻尼因子大小进行调整，保证稳定收敛，并加快收敛速度，主要计算步骤如下：

1) 建立目标函数ϕ=f^Tf，f=AU-P；

2) 给定初值U₀，允许误差ε，初始阻尼因子λ₀>0，令迭代次数k=0；

3) 计算方程(H_k+λ_kI)δ_k=g_k，g_k=-J^Tf_k，H_k为海森矩阵，J为雅可比矩阵，可求解δ_k，当δ_k < ε，停止，得最优解U^*=U_k-1+δ_k；

4) 令U_k+1=U_k+δ_k，根据Fletcher策略修改λ值，令k=k+1，返回步骤3，迭代计算直至计算精度满足要求。

1.2 全空间响应校正

在波场变换之前必须对全空间响应进行校正，由于煤矿井下全空间探测条件下，探测的位置位于井下巷道，与地面半空间探测不同，此时线圈的响应为线圈平面两侧的响应，直接利用井下采集的全空间信号进行波场变换数值计算，显然是错误的，Kaufman等^[18]经过研究得出在均匀介质全空间与半空间的瞬变响应存在倍数关系，即

$ B_{z}^{\mathrm{W}}=k B_{z}^{\mathrm{H}}, $

(4)

式中：B_z^W为全空间瞬变电磁响应，B_z^H为半空间瞬变电磁响应。如满足式(5)条件，式(4)中k取值为2.5^[19]：

$ \varphi=2 {\rm{ \mathsf{ π} }} \sqrt{\frac{2 \rho t}{\mu}}, \varphi / a>15, $

(5)

式中：φ为函数表达式，ρ为介质电阻率，t为采样延迟时间，μ为介质磁导率，a为等效发射线圈半径。在井下瞬变电磁探测情况下，小回线源发射线圈边长不大于2 m，煤岩地层电阻率一般大于20 Ω·m，可得出在矿井瞬变电磁观测时间段内，表达式φ值满足式(4)，因此，可得到全空间与半空间响应磁场的关系，便于进行全空间波场变换计算及后续处理。

1.3 波场数据处理

在矿井瞬变电磁探测方法中，式(2)中P通过仪器数据采集预处理测得，由于瞬变电磁场产生的误差会对变换后的波场造成较大的影响，特别是式(2)中系数矩阵条件数较大时，这种影响会更大，所以在波场变换之前，采用小波变换法进行滤波，利用式(4)进行响应校正，再通过式(2)求解U，在完成波场变换后，采用最小平方反褶积对转换波场数据进行处理^[20-21]，提高变换波场数据的分辨率，得到单位测点压缩后的虚拟波场值，再进行合成孔径相关分析处理，关于虚拟波场传播速度的计算研究，可利用基于“烟圈”理论的最简化反演，由瞬变电磁扩散场数据获取等效涡流环的传播速度^[22-23]，也可利用基于虚拟波场波形特征的速度分析反演成像方法进行速度分析^[24]，这里采用第一种方法计算虚拟波场波速，进而完成时深转换计算过程。

2 采空区探测波场成像方法 2.1 地面瞬变电磁合成孔径原理

地面瞬变电磁波场成像中的相关合成借鉴了雷达中的合成孔径技术，在波场变换后，各个测点的转换波场信号可等效于地震勘探中自激自收的波动场信号，且该信号具有相关叠加性^[25-26]，图 1为地表瞬变电磁合成孔径成像原理示意，线圈在测点所在测线做定向移动, 各个测点测得的信号值相当于小孔径信号值，可利用相关分析的方法将小孔径的信号值合成为等效大孔径的信号值，其公式为：

$ U^{\prime}\left(r_{i}, t_{j}\right)=\sum\limits_{i=-N}^{N} W^{m}\left(r_{i}, \tau_{i}^{m}\right) U\left(r_{i}, t_{j}-\tau_{i}^{m}\right) \quad(j=1,2, \cdots, m) 。$

(6)

式中：左端U′(r_i, t_j)为第i点t_j时刻合成孔径波场值，右端U(r_i, t_j)为t_j时刻第i点的波场值；合成孔径总长度为2N+1，-N~N为线圈定向移动路径；第i点到-N~N内某点的距离为r_i，τ为相对时移，τ_i^m为最佳时移，W^m(r_i, τ_i^m)为最佳时移对应的最大归一化互相关系数。合成孔径示意见图 1，通过相关叠加，异常体对应的地表测点接收到的异常响应可得到增强。

2.2 矿井瞬变电磁合成孔径方法

在井下工作面区域探测工作中，线圈在顺槽巷道内测线做定向移动, 将井下工作面测线每一个测点等效为小孔径，与地面瞬变电磁合成孔径成像公式原理相似，可按照式(6)进行相关合成。图 2为波场成像方法示意图，在井下超前探测工作中，将井下超前探测各个测点角度信号值等效为小孔径，对其各个角度的小孔径信号进行相关合成，其合成孔径成像公式为：

$ U^{\prime}\left(r_{M+N}, t_{j}\right)=\sum\limits_{i=m}^{M+2 N} \sum\limits_{k=1}^{l} W_{i k} U\left(r_{i}, t_{j k}-\tau_{i j k}\right), $

(7)

式中：左端U′(r_i, t_j)为X_i测点t_j时刻合成孔径后波场值，右端U(r_i, t_k)为X_i测点即r_i探测方向测点第k次测量时的波场值，M为起始探测角度编号，M+2N为终止探测角度编号，l指单一角度重复测量次数，τ_ijk为相关系数取极值时的时移量，相关系数为W_ik，按式(8)计算：

$ W_{i k}=\frac{\sum\nolimits_{j=1}^{n} U\left(r_{i}, t_{j k}\right) U\left(r_{M+N}, t_{j k}-\tau_{i j k}\right)}{\left\{\sum\nolimits_{j=1}^{n}\left[U\left(r_{i}, t_{j k}\right)\right]^{2} \cdot \sum\nolimits_{j=1}^{n}\left[U\left(r_{M+N}, t_{j k}\right)\right]^{2}\right\}^{\frac{1}{2}}}。$

(8)

由雷达及地面瞬变电磁合成孔径成像原理，合成孔径成像方法是进行不同测点数据间的相关叠加，井下常见的2种工作方式中，不同测点或不同方向X处的探测结果均可等效为小孔径信号值，当相邻测点的波场信号均对异常体存在反应，那么根据合成公式对各测点的小孔径波场信号值进行等效相关合成，由合成孔径原理可知，合成后异常体信号值可得到加强。

需要注意的是，岩石电阻率具有各向异性，煤岩层沿层理方向和垂直层理方向的导电性不同，地表瞬变电磁采空区探测为垂向探测，探测的是煤层层位采空区及采动裂隙横向上的电性差异，采用大回线装置，测点间距、探测距离范围一般较大，井下瞬变电磁探测探测为煤岩层顺层方向探测，采用多匝小回线装置，测点间距、探测距离范围一般小于地面，在异常体大小相近的情况下，可以适当增大合成孔径的长度，有利于提高合成孔径的效果。

3 理论模型计算

为了验证采空区波场成像探测方法的有效性，设计了全空间中包含低阻块状异常体的三维模型，采用中心回线激发瞬变电磁场，采用有限元法进行数值计算获取正演数据^[27]，并用本文的算法进行波场变换成像。三维模型示意图如图 3(a)所示，发射边长为2 m矩形重叠回线，收发线圈中心点在z=0平面沿x方向作测线参数化扫描，图中共11个测点，测点间距为50 m，围岩电阻率ρ_w为100 Ω·m，低阻块状异常体模拟含水采空低阻区域，位于测线正下方，埋深50 m，大小为200 m×200 m×50 m，电阻率ρ_s为1 Ω·m，图 3(b)为主测线剖面成像结果，红色矩形实线为低阻块状异常，在测线剖面图中，可以反映纵向两个电性界面的存在，深部界面振幅相对较小，相邻深部存在异常体的测点波场曲线相关性较强，震动较为明显，在相关性较弱的异常体的侧边缘，波场曲线震动减弱，证明了方法探测采空区边界的有效性。

4 工程应用 4.1 煤层气富集矿井采空区地球物理特征分析

煤层气富集矿井采空区及采动裂隙内的积水与煤层气共存，积水对采空区冒裂带、采动裂隙内的游离气的赋存、运移产生重要影响。采空区内的积水与水基压裂液作为低阻目标体，与煤层气在同一采空范围之内，使采空区域与周围岩体之间存在明显的电性差异，由前文所述方法原理与理论模型的计算结果可知，波场变换与成像可以反映含水采空区、裂隙带与围岩之间的电性差异，与地震勘探自激自收剖面类似，在地层中电性界面位置附近，波场曲线会出现较为明显的幅值跳动，电性差异越大，跳动越明显，可以综合视电阻率探测结果与异常电性界面波场成像结果对采空区及采动裂隙的电性边界进行判断，更为精确地确定采空区的范围，进而为采空区岩体孔隙、积水量、含水饱和度等模型的计算与煤层气富集矿井资源量的评价工作提供帮助。

晋城矿区煤系地层为石炭系、二叠系地层，即二叠系下统山西组(P₁s)、石炭-二叠系太原组(C₂-P₁)t，主要可采煤层为3号、9号及15号煤层，各煤层煤变质程度已达到无烟煤至中级阶段；3层煤层的底板均为泥岩，渗透率低，各岩层透气性较弱，围岩条件作为煤层的盖层，以及采空区的存在，均为煤层气的赋存提供了良好条件。将该技术在晋城矿区煤层气富集矿井进行了工程应用，并结合巷道超前探测与工作面探测的应用实例对方法的有效性进行了分析。

4.2 巷道掘进超前探测

在晋城矿区某矿井53093巷迎头1 638 m处采用矿井瞬变电磁法进行巷道掘进超前探测，探测采用重叠回线装置，收发线圈边长2 m，发射电流3.5 A，在巷道掘进方向顺煤层扇形断面布置测点，测点角度间隔15°，共13个测点，顺时针方向依次为1~13号测点，探测结果如图 4所示。从图 4可以看出，在巷道掘进正前方偏左发现1处低阻异常区域，位于1~4号测点探测方向。超前探测波场结果如图 5所示，从图 5可以看出，经过波场变换与合成孔径成像，在左侧1~4号测点波场曲线震动较为明显，在该区域存在比较明显的电性界面，深度位置为30~45 m，结合图 4低阻异常响应方向为1~4号测点探测方向，综合分析确定电性界面的位置如图中红色实线圈定位置所示，界面区域采动裂隙较为发育，存在裂隙富水，成像结果对巷道掘进前方的采动裂隙与电性界面反映较为明显，经过钻探验证，3号测点方向孔深33 m处发现采空区且钻孔出水，与综合分析结果吻合较好。

4.3 工作面区域探测

晋城矿区某矿3号煤层局部存在采空区及废弃巷道，为了查明研究矿井3308工作面可能存在的采空区及采动裂隙的影响范围，采用矿井瞬变电磁法进行工作面区域探测，并采用波场变换与合成孔径成像技术对工作面进行成像，探测采用重叠回线装置，收发线圈边长2 m，发射电流4 A，在工作面轨道巷布置测线，测点间距5 m，巷道测线长度500 m，共101个测点。图 6为该工作面探测视电阻率等值线断面图，由传统视电阻率等值线断面图可知，在距离巷道开口0~60，210~305，325~385，400~495 m区段，深度50~100 m范围存在4处低阻异常区域，区域呈现低阻特征，由浅部向深部扩张。图 7为工作面波场成像结果，在图中0~25，255~310，320~365，375~420，440~460，475~490 m区段波场曲线震动较为明显，将对应区段深度为35~42 m范围震动最为明显的波场曲线位置用红色虚线圈定，其中在255~490 m区段波场曲线振幅强且连续性较好，结合地质资料与视电阻率断面图结果分析，推测该区域存在废弃巷道，该电性界面为采空巷道的富水边界，深部区域波场曲线振幅相对较强的位置为含水采动裂隙发育影响。经过钻探验证，在开口切眼区段铁器较多钻孔未发现出水，在270~410 m区段多处钻孔发现采空区且钻孔出水，与综合分析结果吻合较好。

5 结论

1) 由瞬变电磁扩散场与地震波动场满足的数学关系式入手，结合全空间响应校正方法，通过波场变换反演可对井下瞬变电磁资料中的波场信息进行提取，在波场数据处理过程将井下小回线等效为小孔径天线，对不同测点的转换波场信号进行相关合成，理论模型计算结果表明低阻采空区异常纵向电性界面反映明显, 异常体侧边缘波场曲线震动减弱，可分辨采空区边界；

2) 工程应用试验效果表明常用的超前与工作面电阻率断面图可以确定采空区积水总体范围，与常规数据处理相比，波场变换成像技术能反映异常体的界面特征，从而降低体积效应对煤层气富集矿井含水采空区及采动裂隙异常解释的不利影响，综合视电阻率探测结果分析，验证了该方法对煤层气富集矿井采空区与采动裂隙的探测效果，可以应用于煤层气富集矿井采空区的精细探测工作中，在探明采空区的基础上，结合采空区岩体孔隙、积水量、含水饱和度等模型计算方法，可为矿井采空区煤层气资源量的评价工作提供帮助。