基于FDTD的空间步进频雷达煤层异质区探测正演分析

贺 顺，罗彩祯，杨志伟，贾远航

（1.西安科技大学 通信与信息工程学院，陕西 西安 710054；2.西安电子科技大学 雷达信号处理重点实验室，陕西 西安 710071）

0 引 言

现如今由煤层异质区如采空区、富水区、瓦斯突出区等引发的矿井事故多有发生，了解掘进面前方煤体的结构对煤矿安全生产尤为重要。目前，多种地球物理勘探技术已被广泛应用于煤矿井下地质预测中，如矿井直流电法、瞬变电磁法、地震勘探法等。矿井直流法受地形的影响较大，需要花费大量的时间进行设备的铺设；瞬变电磁法存在浅部探测盲区且容易受到铁磁介质的影响；在多种勘探技术中，探地雷达有抗干扰能力强、受地形的影响较小、操作方便、不需要进行庞杂的设备布置等优点。许多学者对探地雷达应用于煤矿井下探测进行了研究。文献[2⁃3]使用地质雷达探测井下瓦斯富集区并对其特征进行研究；文献[4]使用GprMax和Matlab数值软件对井下二维空间磁场模型进行正演模拟，研究不同激励源与中心频率对模拟结果的影响。文献[5]使用超宽带探地雷达开展了煤层厚度的探测试验。井下探测使用的探地雷达多为单发单收雷达系统，这种系统虽然操作简单，探测方便，但需要多次移动进行多点测量，可获得的信息量有限。

多输入多输出探地雷达探测一次可获得多组数据，相比单站雷达能获得更多的信息量，更加充分地利用了回波信号的空时信息，并提高了系统的空间分辨能力。文献[6⁃7]先后对阵列式探地雷达目标成像问题以及阵列式探地雷达信号极化场特征进行了研究；文献[8]对车载阵列探地雷达实时检测地雷进行了研究。目前，多输入多输出探地雷达多应用于检测地下藏雷、管线结构、路面地下空洞等方面，在煤矿井下探测的应用较少。

探地雷达通过增大信号带宽和天线合成孔径来提高雷达探测的距离分辨力和方位分辨力。硬件发展的水平限制了超宽带信号的产生。多输入多输出体制下的空间步进频雷达通过在相邻的阵元上依次发射频率步进的信号，在接收端合成一个宽带信号，解决了探地雷达系统中目标的距离分辨力与雷达信号带宽之间的矛盾，同时相比单站雷达方位分辨力又有所提高。本文主要研究多输入多输出雷达体制下的空间步进频信号在煤矿井下环境的传播，建立二维煤矿井下异质区环境模型，使用时域有限差分方法（Finite⁃difference Time⁃domain，FDTD）进行空间步进频信号探测异质区的数值模拟，分析不同异质区模型下的回波结果，得到对应模型下的传播规律。

1 空间步进频信号模型

空间步进频雷达对煤矿井下异质区进行探测时，将雷达天线阵列放置在井下巷道中。天线阵列同时平行于巷道地面与煤矿开采掘进面，如图1所示。发射天线向掘进面发射电磁波信号，信号在电性参数不同的物质分界面上发生反射和透射现象，接收天线接收反射信号，通过分析接收信号的时间、相位等信息对异质区的位置等特征进行解释。

图1 空间步进频雷达井下放置示意图

电磁波信号在介质中传播的快慢，影响回波响应出现的时间，根据相邻回波响应之间的时间差，可计算电磁波传播的距离。雷达发射信号在介质中的传播速度可根据式（1）求得：

式中：为电磁波在自由空间传播的速度；′为介质的相对介电常数。可知介电常数越大，信号在介质中的传播速度越慢。

空间步进频雷达实现空间上的频率步进信号，使用等间隔的均匀线阵组成线阵的每个阵元发射频率步进的信号。空间步进频雷达发射阵列结构如图2所示，发射天线为个阵元的均匀线阵，其阵元间距为。第个阵元的发射信号x ()为：

图2 空间步进频雷达发射阵列结构

（1）空间位置管理需要采用矢量数据结构。在地理位置表述时，矢量数据具有精确度高、数据量小的特性，以方便全面系统的阐述空间位置关系，对于数据图形、属性及数据恢复均能在短时间内实现，其具体格式如图4。

2 FDTD时域有限差分原理

井下探地雷达探测异质区属于电磁类探测方法。FDTD十分适用于探地雷达建模，因为它能有效地模拟具有不同几何形状和介电特性的地下分层结构。FDTD方法将微分形式的时域麦克斯韦方程转化为差分方程，使用迭代的方法求解有限空间、时间范围内的电场和磁场分量。一般情况下，对于非均匀、各项同性、线性的媒质材料，其麦克斯韦电磁方程组微分形式如下所示：

式中：表示电场强度矢量，单位为V/m；表示磁场强度矢量，单位为A/m；表示电位移矢量，单位为C/m；表示磁感应强度，单位为T；表示电流密度矢量，单位为A/m；J 表示磁通量密度矢量，单位为V/m；表示电荷密度，单位为C/m。

同时有电磁场对材料介质作用的关系式，即本构方程：

式中：表示物质材料的介电常数，单位为F/m；表示物质材料的电导率，单位为S/m；表示物质材料的磁导率，单位为H/m；σ表示计算磁损耗的磁阻率，单位为Ω/m。

FDTD方法对计算空间进行网格的剖分，图3为著名的Yee元胞，元胞在坐标轴方向上的大小称为该方向上的空间步长，规定每一个（或）场分量周围有4个（或）场分量环绕，应用这种离散方式将含时间变量的麦克斯韦旋度方程转化为差分方程。

图3 FDTD离散中的元胞

对于二维问题，所有物理量与无关，∂∂=0。以TM波为例，其麦克斯韦旋度方程为式（8），式（9）~式（11）为推导所得任意节点处电场和磁场的FDTD差分公式，式中分别为（，+1 2），（+1 2，），（，）。

由于使用差分方法对麦克斯韦方程进行数值计算会产生色散现象，即波的传播速度会随波长而改变。这种色散会引起波形的畸变，且网格越粗色散越严重。为了使差分方程的计算结果能更近似于原方程的解，需要对时间步长和空间步长进行一定的限制，在二维问题中应该满足式（12）。

3 正演模拟

3.1 实验模型及参数设置

在煤矿开采的过程中，地质条件复杂多变，将存在于煤体中的非煤区域称为异质区。煤层中常见的异质区有陷落柱、断层、采空区、富水区、瓦斯富集区等。本文对采空区、富水区以及岩层三种主要煤矿异质区进行模拟。未考虑巷道中各种仪器、线缆等产生的干扰，并假设介质为均匀同性，建立掘进工作方向煤层内部的物理模型如图4所示。模型设置异质区上方煤层厚度为0.5 m，采空区厚度为0.5 m，富水区厚度为0.5 m，薄泥岩层厚度为0.5 m。

图4 异质区物理模型

本文使用收发分置的天线阵，收发天线均为5阵元的均匀线阵，阵元间隔为0.1 m。收发天线阵水平放置，与掘进方向煤层表面平行。收发天线位于同一水平面，发射天线阵放置距煤层表面0.8 m，接收天线阵放置距煤层表面0.6 m。发射天线发射时长为1 ns的正弦信号，选定线阵一端的阵元作为第一阵元，其所发射信号的中心频率为200 MHz，其余阵元的中心频率依次以200 MHz的频率增量进行步进，这样第5个阵元发射信号的中心频率为1 GHz。

模拟中的最高频率选择为阵元中最高中心频率的3倍，由最高频率确定的最小波长为0.1 m，计算得到空间步长不应超过0.01 m，本文实验中3个二维异质区模型在方向与方向上的空间网格步长均设定为0.003 m。由式（12）可算得时间步长为7.076 0×10s。本文使用完美匹配层（Perfectly Matched Layer，PML）作为吸收边界条件对空间进行截断。异质区模型中所用介质的电性参数如表1所示。

表1 模型电性参数

3.2 传播特性分析

不改变发射天线以及接收天线的位置，对上述三种异质区进行数值计算，得到接收天线5个阵元处随时间迭代次数变化的电场信号数值。将5个接收点结果进行绘图显示，得到三种异质区的仿真结果如图5~图7所示。

图5 采空区结果

图7 薄泥岩层结果

图5为采空区模拟结果，煤与空气界面表现为一条明显的轴线，出现在第1 157次时间迭代处，通过计算可知该点时间为8.187 ns。采空区上下表面回波对应的反射轴线分别出现在第2 337次与第2 789次时间迭代处，计算得到对应时间分别为16.537 ns和19.735 ns。某位置的时间等于该位置的时间迭代次数乘以时间步进间隔，时间步进间隔可由式（12）求得。

图6为富水区模型的模拟结果，可以清楚地看到在煤与空气界面形成的响应后仅有一条清楚的响应轴线，该反射轴线表示煤层与富水区上表面的回波。由于水的介电常数大，电磁波在水中的衰减速度快。当信号穿过煤层进入富水区之后很快被吸收，因此仅能看到煤与水上表面形成的回波响应。

图6 富水区结果

图7为薄泥岩层模型的模拟结果，与采空区模型相似，同样可以看到煤与泥岩上下交界面处的两个回波，但与采空区相比较为微弱。

实验设置天线阵平行于掘进煤层表面，即发射信号为垂直入射，根据双程走时′和速度可计算煤体到异质区上表面厚度=′2，电磁波在煤中的传播速度利用式（1）求得。表2为根据实验数据对模型厚度进行计算的结果，表3为对模型内介质速度进行计算的结果。速度从大到小排序为空气、煤、泥岩、水。根据仿真数据得到的计算结果与理论值存在一定的误差，厚度的误差不超过0.05 m，速度的误差不超过8%。除仿真过程中产生的误差外，还有结果读取和计算时不精确造成的误差。

表2 模型厚度分析 m

表3 模型介质速度分析

异质区的物理模型不变，改变异质区上方煤层厚度为5 m，采空区、富水区、泥岩区三种异质区的厚度均为5 m，和方向的空间步长调整为0.01 m，其余参数不变，得到的结果如图8~图10所示。同样为采空区与泥岩区表现为两条响应轴线，泥岩区的响应较弱，富水区表现为一条响应轴线，与厚度为0.5 m模型下的模拟结果一致，仅到达时间不同。

图8 5 m厚采空区模拟结果

图9 5 m厚富水区模拟结果

图10 5 m厚泥岩区模拟结果

为了更好地观察空间步进频信号在整个仿真区域的传播情况，计算了三种模型空间中所有网格的电场强度，绘制频率为1 GHz情况下二维仿真模型平面上的电场分布图，如图11所示。

图11 三种模型的电场分布图

电场分布图中的横纵坐标分别为模型掘进方向与垂直掘进方向上的距离坐标。设置=0的界面为煤矿开采工作面，横坐标的正半轴表示掘进面上方的巷道；横坐标的负半轴表示掘进面下方的煤层内部。

右侧各发射阵元发射频率步进的信号向左侧煤层中进行传播。图11a）为采空区模型下的电场分布情况，电磁波信号由自由空间传入煤层中，采空区部分波纹与煤层上表面自由空间处波纹一致；图11b）为富水区的电场分布情况，由于能量衰减快在传播中消耗殆尽，电磁波基本不能穿过水层；图11c）为薄泥岩层的电场分布情况，电磁波从煤层到泥岩中波纹变密。电磁波信号传播过程中波纹有密有疏，波纹的稀疏情况随波长大小发生变化。在频率不变的情况下，波速越大波长就越大，波纹越稀疏，根据电场分布图得到模型中介质速度由大到小排序分别为空气、煤、泥岩、水，与理论一致。

4 结 语

本文基于时域有限差分方法对空间步进频信号在三种不同井下异质区模型中的传播情况进行正演模拟，模拟结果表明：通过使用空间步进频信号的时域波形进行数值计算，得到的煤层厚度与正演模型设置的参数基本吻合；在三种异质区正演结果中，采空区模型与薄泥岩层模型下均可看到信号在异质区上下表面产生的两条反射轴线，但采空区的回波更为明显，而富水区模型下仅可看到水煤上交界面的反射轴线；根据模拟的电场分布图，通过波纹疏密程度可以对介质的速度大小进行大致的排序；信号传播时衰减的快慢不同，在水中的衰减最大，泥岩次之，空气衰减最小。

由以上仿真结果可知：进行实际探测时，若探测数据出现明显的回波异常区域，则可以判断该处存在与煤体不同的异质区，并可根据异常回波出现的时间点对异质区域上方的煤层厚度进行估计；若信号在产生响应后急速衰减则可判断该区域存在水，验证了使用多发多收体制下的空间步进频雷达对煤体异质区探测的可行性。