0 引言
目前，矿井突水已成为仅次于煤与瓦斯突出的第二大煤矿灾害。矿井富水区超前探测的主要方法为瞬变电磁法。国内外学者针对半空间的瞬变电磁法做了大量研究[1-4]，但对矿井全空间瞬变电磁法的研究很少，主要是将全空间分为上下半空间来计算层状介质的瞬变电磁响应特征[5-6]；针对矿井瞬变电磁法的一维正演、反演计算做了大量的基础理论和实际应用研究，对三维正演有一定的研究，但对三维反演和成像的研究比较少，尚未形成完整的理论体系。时域有限差分(Finite-Difference Time-Domain，FDTD)方法是目前常用的矿井富水区超前探测瞬变电磁法，其时间步长的选取受限于Courant-Friedrich-Lewy(CFL)稳定性条件，即受限于空间步长。为了确保成像精度，须将成像探测区域剖分成足够密的网格。但网格剖分密度越大，FDTD方法的计算效率越低。Crank-Nicolson时域有限差分(CN-FDTD) 方法能够解决传统FDTD方法受限于CFL稳定性条件的问题，其时间步长选择不受空间步长选择的限制，在相同网格剖分密度和成像精度条件下，可明显提高计算效率，节约计算空间[7-11]。
本文首次将逆时偏移成像算法和CN-FDTD方法引入矿井富水区陷落柱成像研究，建立了矿井富水区陷落柱三维空间模型，分析了激励源线圈频率和角度对成像分辨率的影响，并对富水区陷落柱进行了成像试验。
1 逆时偏移成像算法
逆时偏移成像包括正向延拓、逆时延拓和成像3个步骤：正向延拓是波场沿着时间轴正方向传播到最大时刻；逆时延拓是从接收波场的最终时刻开始，向时间轴的负方向延拓至零时刻；成像为应用成像条件得到矿井富水区异常体的真实构造信息[12-14]。根据文献[15-17]，通过互相关成像条件可得探测点(x,y,z)在矿井富水区的偏移成像结果：
(1)
式中：S(x,y,z,t)为t时刻(x,y,z)处的时间域正向延拓波场；R(x,y,z,t)为t时刻(x,y,z)处的时间域逆向延拓波场。
S(x,y,z,t)R(x,y,z,t)为t时刻对整个矿井三维空间的波场做一次成像运算并存储，然后将所有成像结果进行求和叠加而获得的最终偏移成像结果。考虑到矿井全空间三维条件下的计算效率和存储空间问题，采用激发时间成像条件进行富水区三维空间成像，可节约大量物理空间和读取时间，具有较好的保幅性[18-19]。井下富水区激发时间成像条件：
(2)
(3)
式中:(x′,y′,z′)为富水区异常体在逆推过程中某空间位置坐标;td(x′,y′,z′)为激励源到富水区异常体的初至时间。
2 CN-FDTD方法
逆时偏移成像算法的关键在于正演问题，本文首次将CN-FDTD方法作为逆时偏移成像的正推算法引入矿井富水区探测。CN-FDTD方法用于数值求解热方程和类似的偏微分方程，是一种时间域上的二阶方法。该方法隐含在时间上，可写为隐式Runge-Kutta方法，并且在数值上是稳定的[20]。
以Maxwell方程组三维FDTD离散形式的电场强度计算为例说明CN-FDTD方法原理。
(4)
式中：Ex为矿井全空间介质x轴方向的电场强度；εx为矿井全空间介质x轴方向的介电常数；Hy，Hz分别为井下全空间介质y轴、z轴方向的磁场强度。
将式(4)展开可得
(5)
式中：n为时间步长指标；Hy(i，j，k)，Hz(i，j，k)分别为矿井全空间介质在点(i，j，k)处的y轴、z轴方向磁场强度。
将式(5)后一时刻形式中前向差分替换为后向差分，且由于
对其取2次中心差分可得
(6)
式中：Δt为时间步长；Δx,Δy,Δz为矿井探测区域沿x轴、y轴、z轴剖分的网格尺寸。
式(6)可改写为
(7)
矿井全空间介质y轴、z轴方向的电场强度计算与此类似。
3 矿井富水区陷落柱成像试验
3.1 陷落柱三维空间模型建立
在矿井掘进过程中，富水区陷落柱已造成多次突水事故。陷落柱一般呈圆柱状，在导通水源之后易充水，逐步形成富水区。结合矿井地质结构，建立陷落柱三维空间地球物理模型，采用CN-FDTD逆时偏移成像算法探测富水区陷落柱的具体形状、位置和范围。
掘进工作面前方陷落柱三维空间模型如图1所示。模型在y轴正方向24 m处分层。分界面左侧为巷道，其介质为空气，相对介电常数为1;右侧取相对介电常数为4.09的褐煤作为煤质。激励源线圈和接收线圈采用偶极装置，其中激励源线圈为边长1 m的正方形，位于分界面左侧5 m处，且平行于xoz平面，激励源函数采用Ricker；接收线圈为边长2 m的正方形，位于分界面上，且平行于xoz平面。陷落柱走向与巷道掘进方向垂直。
图1 陷落柱三维空间模型
Fig.1 Three-dimensional model of karst collapse pillar
3.2 激励源线圈频率选择
通过试验研究激励源线圈频率对矿井富水区陷落柱成像分辨率的影响，并根据试验结果选取合适的激励源线圈频率。分别取峰值频率为25，35，45，55，65 MHz的Ricker作为激励源。陷落柱模型尺寸为40 m×80 m×40 m(x×y×z)，x，y，z方向的网格剖分密度为200，300，200。陷落柱半径设置为10 m，其中心轴线与分界面相距28 m。模型其他结构均与图1相符。
矿井富水区陷落柱在xoy平面的偏移成像如图2所示。可看出分界面的位置大约位于y=24 m处，能量最强区域即为接收线圈所在区域，为y=19～24 m，与陷落柱三维空间模型相符。随着Ricker峰值频率的增加，激励源线圈左侧波纹的波峰与波谷间距逐渐减小。当激励源线圈发出的电磁波从空气传播至煤层时，其传播方向在分界面处发生改变，进入煤层后遇到富水区陷落柱产生强烈的电磁反射波。随着激励源线圈频率的增大，成像分辨率变高，但探测深度会变小，因此应当选取频率适合的Ricker作为激励源。图2中，当Ricker峰值频率为65 MHz时，陷落柱轮廓清晰可见，成像分辨率较高，因此本文选用峰值频率为65 MHz的Ricker作为激励源。
(a) 激励源线圈峰值频率为25 MHz
(b) 激励源线圈峰值频率为35 MHz
(c) 激励源线圈峰值频率为45 MHz
(d) 激励源线圈峰值频率为55 MHz
(e) 激励源线圈峰值频率为65 MHz
图2 不同激励源线圈峰值频率下陷落柱xoy平面偏移成像
Fig.2 Migration imaging of karst collapse pillar inxoyplane under different peak frequency of excitation source coil
3.3 激励源线圈角度选择
通过试验研究激励源线圈角度对矿井富水区陷落柱成像分辨率的影响，并根据试验结果选取合适的激励源线圈角度。选取方形激励源线圈，分别垂直于xoz平面、与xoz平面呈45°夹角、平行于xoz平面。陷落柱模型尺寸为40 m×130 m×40 m(x×y×z)，x，y，z方向的网格剖分密度为180，350，180，陷落柱半径设置为5 m，其中心轴线与分界面相距86 m。选取峰值频率为65 MHz的Riker作为激励源，模型其他结构与图1相符。
不同激励源线圈角度下富水区陷落柱在xoy平面的偏移成像如图3所示。可看出分界面位置大约在y=24 m处，能量最强区域即为接收线圈所在区域，为y=19～24 m，均与陷落柱三维空间模型相符。当激励源线圈发出的电磁波从空气传播至煤层时，其传播方向在分界面处发生改变。当方形激励源线圈与xoz平面呈45°夹角或平行于xoz平面时，电磁波遇到富水区陷落柱时产生强烈的电磁反射波。当激励源线圈垂直于xoz平面时，无法看到陷落柱xoy平面的偏移成像；当激励源线圈与xoz平面呈45°夹角时，陷落柱左侧轮廓较清晰，右侧轮廓较模糊；当激励源线圈与xoz平面平行时，陷落柱轮廓清晰，与陷落柱三维空间模型相符。因此，本文选择平行于xoz平面的方形激励源线圈。
(a) 方形激励源线圈垂直于xoz平面
(b) 方形激励源线圈与xoz平面呈45°夹角
(c) 方形激励源线圈平行于xoz平面
图3 不同激励源线圈角度下陷落柱xoy平面偏移成像
Fig.3 Migration imaging of karst collapse pillar inxoyplane under different angle of excitation source coil
3.4 富水区陷落柱成像
上述研究表明方形激励源线圈峰值频率为65 MHz且平行于xoz平面时成像效果最好，因此选取上述参数对富水区陷落柱进行成像试验。陷落柱模型尺寸为40 m×80 m×40 m(x×y×z)，x，y，z方向的网格剖分密度为200，300，200，陷落柱半径设置为10 m，其中心轴线与分界面相距28 m，模型其他结构均与图1相符。
矿井富水区陷落柱xoy平面和yoz平面模型如图4所示，偏移成像如图5所示。从图5可看出，分界面位置大约在y=24 m处。当激励源线圈发出的电磁波从空气传播至煤层时，其传播方向在分界面处发生改变。采用CN-FDTD逆时偏移成像算法得到的矿井富水区陷落柱成像在y=42 m处出现强烈的电磁反射波，与图4相符。陷落柱xoy平面和yoz平面偏移成像中，其轮廓均与图4相符。根据xoy平面和yoz平面偏移成像可描绘出富水区陷落柱的三维空间模型。
图4 陷落柱xoy平面和yoz平面模型
Fig.4 Karst collapse pillar models inxoyplane andyozplane
(a)xoy平面
(b)yoz平面
图5 陷落柱xoy平面和yoz平面偏移成像
Fig.5 Migration imaging of karst collapse pillar inxoyplane andyozplane
4 CN-FDTD方法计算效率分析
在相同的模型和网格剖分密度下，对基于CN-FDTD方法和传统FDTD方法的矿井富水区陷落柱成像的计算效率进行试验对比。陷落柱模型尺寸为40 m×80 m×40 m(x×y×z)，x，y，z方向的网格剖分密度为200，300，200，设置80个观测点。陷落柱半径设置为10 m，其中心轴线与分界面相距28 m。在分析计算效率时选择与图1相同的模型。
CN-FDTD方法和FDTD方法的数值计算特征曲线如图6所示。在数值计算过程中，时间步长Δt与空间步长Δx，Δy，Δz之间必须满足CFL稳定性条件。而CN-FDTD方法是无条件稳定的，时间步长Δt的选取不受空间步长Δx，Δy，Δz的限制，且在计算过程中采用多线程模式。CN-FDTD方法在时间步长取值远大于CFL稳定性条件时，其与FDTD方法的数值计算特征曲线拟合良好。
(a) 数值计算特征曲线
(b) 局部放大
图6 CN-FDTD方法与传统FDTD方法的数值计算特征曲线
Fig.6 Numerical calculation characteristic curves of CN-FDTD method and traditional FDTD method
CN-FDTD方法与传统FDTD方法的性能指标对比见表1。可看出CN-FDTD方法的运行时间为传统FDTD方法的1/3，内存消耗约为传统FDTD方法的87%，从而证明CN-FDTD方法计算效率高，内存占比小，用于矿井复杂环境时，CN-FDTD方法具有明显优势。
表1 CN-FDTD方法与FDTD方法的性能指标对比
Table 1 Performance index comparison between CN-FDTD method and traditional FDTD method
5 结语
将逆时偏移成像算法和CN-FDTD方法引入矿井富水区陷落柱三维成像研究中，解决了传统FDTD方法受限于CFL稳定性条件的问题，提高了矿井瞬变电磁法的计算效率，节省了内存成本。试验表明，CN-FDTD逆时偏移成像算法可对矿井富水区陷落柱进行高分辨率成像，精确确定其形状、位置和范围。该算法为矿井富水区陷落柱超前探测提供了一种新思路。
参考文献：
[1] 张鹏.中国煤炭矿井物探技术现状及展望[J].工矿自动化,2017,43(3):20-23.
ZHANG Peng.Status and prospect of coal mine geophysical exploration technology in China[J].Industry and Mine Automation,2017,43(3):20-23.
[2] 周金,程久龙,温来福.矿井瞬变电磁法反演方法研究进展与展望[J].煤矿安全,2017,48(4):180-183.
ZHOU Jin,CHENG Jiulong,WEN Laifu.Research progress and prospect on inversion methods of mine transient electromagnetic method[J].Safety in Coal Mines,2017,48(4):180-183.
[3] 孟强华,窦文武.瞬变电磁探测在煤矿采空区积水边界探测中的应用[J].煤炭技术,2018,37(2):123-125.
MENG Qianghua,DOU Wenwu.Application of transient electromagnetic detection in goaf water boundary detection in coal mine[J].Coal Technology,2018,37(2):123-125.
[4] 王鹏.厚煤层露天矿地下水瞬变电磁法探测[J].煤矿安全,2019,50(1):153-156.
WANG Peng.Detection of groundwater in thick coal seam open pit by TEM[J].Safety in Coal Mines,2019,50(1):153-156.
[5] ROY K K,VERMA S K,MALLICK K.Deep electromagnetic exploration[M].Berlin:Springer,1999.
[6] KRIVOCHIEVA S,CHOUTEAU M.Whole-space modeling of a layered earth in time-domain electromagnetic measurements[J].Journal of Applied Geophysics,2002,50(4):375-391.
[7] 常江浩.煤矿富水区矿井瞬变电磁响应三维数值模拟及应用[D].徐州:中国矿业大学,2017.
[8] 张广博.掘进工作面富水区瞬变电磁法多分量探测物理模拟及应用[D].徐州:中国矿业大学,2016.
[9] 许新刚.矿井含水构造瞬变电磁场响应及透视技术研究[D].徐州:中国矿业大学,2017.
[10] 于景邨,常江浩,苏本玉,等.老空水全空间瞬变电磁法探测三维数值模拟研究[J].煤炭科学技术,2015,43(1):95-99.
YU Jingcun,CHANG Jianghao,SU Benyu,et al.Study on whole space transient electromagnetic method prospect three dimensional numerical modeling of gob water[J].Coal Science and Technology,2015,43(1):95-99.
[11] 刘尧.基于有限差分的瞬变电磁三维数值模拟[D].北京:中国地质大学(北京),2015.
[12] 徐兴荣,王西文,王宇超,等.基于波场分离理论的逆时偏移成像条件研究及应用[J].地球物理学进展,2012,27(5):2084-2090.
XU Xingrong,WANG Xiwen,WANG Yuchao,et al.Study and application of imaging condition for reverse-time migration based on wave-fields separation[J].Progress in Geophysics,2012,27(5):2084-2090.
[13] 朱尉强,黄清华.探地雷达衰减补偿逆时偏移成像方法[J].地球物理学报,2016,59(10):3909-3916.
ZHU Weiqiang,HUANG Qinghua.Attenuation compensated reverse time migration method of ground penetrating radar signals[J].Chinese Journal of Geophysics,2016,59(10):3909-3916.
[14] 刘学建,刘伊克.表面多次波最小二乘逆时偏移成像[J].地球物理学报,2016,59(9):3354-3365.
LIU Xuejian,LIU Yike.Least-squares reverse-time migration of surface-related multiples[J].Chinese Journal of Geophysics,2016,59(9):3354-3365.
[15] 鲁兴林,钱荣毅.地质雷达有限差分逆时偏移方法研究[J].地球物理学进展,2017,32(2):885-890.
LU Xinglin,QIAN Rongyi.Ground-penetrating radar finite-difference reverse time migration[J].Progress in Geophysics,2017,32(2):885-890.
[16] LU Xinglin,SONG Ao,QIAN Rongyi,et al. Anisotropic reverse-time migration of ground-penetrating radar data collected on the sand dunes in the Badain Jaran Desert[J].IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing,2018,11(2):647-654.
[17] 程飞.多井井间三维波场数值模拟与逆时偏移成像研究[D].武汉:中国地质大学,2015.
[18] 王京.隧道空间有限元地震波场模拟与逆时偏移成像[D].武汉:中国地质大学,2017.
[19] YANG Y.The three-dimensional unconditionally stable FDTD algorithm based on Crank-Nicolson method[C]//IEEE Antennas and Propagation Society International Symposium,Albuquerque,2006:81-84.
[20] WANG Zhiyong,DING Hao,LU Guijin,et al. Reverse-time migration based optical imaging[J].IEEE Transactions om Medical Imaging,2016,35(1):273-281.