0.   引　　言

煤岩界面识别是实现井工煤矿智能无人化开采的关键^[1-2]。受工作面高粉尘、低照度及复杂地质条件等因素制约，诸多煤岩界面识别方法存在一定的应用局限性^[3]。探地雷达方法具有实时性高、适用范围广、非接触性等优点^[4]，可实现煤岩界面高精度随采动态探测^[5]。

煤岩界面探地雷达波组自动拾取技术是煤岩层位在线智能解译的基础^[6]。煤岩界面通常表现为非固定水平面分界，且具有一定粗糙度。根据瑞利准则，当煤岩界面均方根高度大于入射波波长的1/8时，探地雷达信号受显著影响，被认为是粗糙表面^[7]。利用高频探地雷达识别粗糙煤岩界面位置时，B-Scan图像主要特征包括由煤岩界面产生的连续曲线特征和局部粗糙界面产生的双曲线特征。此外，电磁波撞击粗糙煤岩界面时会产生多径效应，导致后向散射干扰波与目标界面回波在局部位置发生干涉，煤岩界面连续曲线特征被干扰双曲线破坏，为煤岩界面波组自动拾取带来更大挑战。

传统煤岩界面直接拾取方法适用于简单或理想条件下煤岩层位的追踪^[5]。许献磊等^[5]针对“串层”问题，提出了三级窗口算子，通过分析参考道与被追踪道之间数据序列的相关性，逐级调整窗口范围，以减少“串层”现象的发生。该方法可以有效提高煤岩层位追踪的准确性和稳定性，但对粗糙煤岩界面的适用性、参考道及窗口自动选择等方面仍存在优化的空间。

关于探地雷达目标界面波组自动拾取方法的研究多以双曲线特征检测为主。Radon变换^[8]和Hough变换^[9]应用较为广泛，将图像空间像素点映射到Radon空间或参数空间，利用累加器矩阵检测B-Scan图像的双曲线特征。为了减小参数空间的计算复杂度，通常引入边缘检测算法进行特征提取的预处理，如Sobel算子^[10]和Canny算子^[11]等。此外，模板匹配、统计分析、机器学习和深度学习等方法也相继提出并用于双曲线特征的检测。PASOLLI等^[12]通过预先设计标准模板，将其与分割后的GPR图像进行匹配，并使用汉明距离作为相似度测度来评估匹配效果。GADER等^[13]使用马尔可夫模型对GPR图像中的地雷目标信号进行建模，通过计算信号序列的概率来确定目标的位置。人工神经网络方法被应用于地雷探测^[14]和管线检测^[15]中的目标分类^[16]。深度卷积神经网络可用于复杂场景下的双曲线特征检测，如沥青路面破坏损伤检测^[17]、地下空间目标识别^[18]等。

利用边缘检测算法提取潜在目标界面的响应特征，并针对干扰双曲线进行检测可以间接拾取煤岩界面波组。但基于B-Scan图像的双曲线检测方法通常涉及了大量的计算资源和时间成本，并且不适用于粗糙煤岩界面导致的强电磁干涉干扰下的复杂场景。

笔者提出一种粗糙煤岩界面探地雷达波组自动拾取方法。以探地雷达A-Scan序列数据中波峰或波谷点为锚点^[19-20]，引入局部统计量自适应调整的振幅阈值和锚点间距阈值，以提高对较大埋深煤岩界面地质模型的鲁棒性。考虑粗糙煤岩界面曲线锚点分布特性，提出锚点邻域内就近生长方法，利用振幅变异系数作为锚点终止生长的判断依据，以解决电磁波干涉效应导致的多锚点或虚假曲线锚点问题。采用RANSAC迭代拟合算法和波形特征匹配方法进行锚点拟合和分类，以准确识别干扰双曲线特征和煤岩界面曲线特征，并利用拟合曲线顶点位置排除部分干扰双曲线。针对特殊地质条件下煤岩界面走向表现为双曲线特征问题，提出基于电磁波几何扩散效应的曲线锚点振幅补偿方法，并利用补偿后锚点振幅变化规律进行干扰双曲线的虚假检测控制。

1.   粗糙煤岩界面建模与正演模拟

1.1   粗糙表面生成方法

粗糙煤岩界面可以简化为一维高斯粗糙表面模型^[21]。均方根高度和相关长度是描述其粗糙程度的重要统计参量，前者量化了表面高度的波动程度，后者描述了表面形貌的横向尺度。

均方根高度定义为

式中：$ {\delta }_{\mathrm{h}} $为均方根高度，m；$ h $为粗糙表面任一点相对于参考平面的高度，m；$ <\cdot >$为集合平均。

对于特定分布的高斯粗糙表面，均方根高度无法准确描述其粗糙特性。自相关函数表明了粗糙表面任意2点间的关联程度，高斯粗糙表面自相关函数定义为

式中：$ R\left({x}_{\mathrm{d}}\right) $为自相关函数；$ {x}_{\mathrm{d}} $为一维粗糙表面上任意2点之间的水平距离，m；$ l $为相关长度，m。相关长度定义为$ R\left({x}_{\mathrm{d}}\right) $衰减到初始值1/e对应的$ {x}_{\mathrm{d}} $值。

一维高斯粗糙表面可以利用蒙特卡罗方法^[22]模拟生成。其数学表达式为^[21]

式中：$ L $为粗糙表面长度，m；$ {x}_{n} $为粗糙表面第n个离散点；$ {K}_{j} $为离散波数，$ {K}_{j}=\text{2}\mathrm{\pi }j/l $，m⁻¹；$ F\left({K}_{j}\right) $为$ f\left({x}_{n}\right) $的空间傅里叶变换，定义为

式中：$ \Delta K $为谱域相邻的谐波样本的空间波数差；$ N\left(\mathrm{0,1}\right) $为均值为0、方差为1的高斯随机变量；$ W\left({K}_{j}\right) $为粗糙表面的功率谱密度，定义为

图1为利用蒙特卡罗方法生成的不同粗糙程度的高斯粗糙表面，均方根高度和相关长度全面描述了表面粗糙特性。对比图1a—图1c，随着均方根高度增大，表面起伏深度增加，随着相关长度减小，表面起伏频率增大。图1d和图1e中，高频分量对应快速起伏，低频分量对应慢速起伏，随着相关长度减小，高频分量增加，即表面起伏频率增大。

图  1  不同粗糙程度的高斯表面及起伏频率频谱

Figure  1.  Gaussian surfaces with different roughness levels

下载: 全尺寸图片 幻灯片

1.2   粗糙煤岩界面模型及正演模拟

建立如图2b所示粗糙煤岩界面地质模型。模型尺寸为1.5 m×1.2 m；空气层平均厚度为0.3 m；岩层的复杂走向采用样条插值方法定义，插值点为（0,1）、（0.3,0.9）、（0.6,1）、（0.9,0.95）、（1.2,0.85）、（1.5,0.5）；空气层、煤层、岩层的相对介电常数和电导率^[23]分别为$ \left[\text{1,0}\right] $、$ \left[\text{3,2}\text{×}{\text{10}}^{{-}\text{5}}\right] $、$ \left[\text{6,5}\text{×}{\text{10}}^{{-}\text{5}}\right] $；在空气‒煤界面、煤岩界面分别添加高斯粗糙表面，其均方根高度和相关长度分别为$ \left[\text{0.1,0.01}\right] $、$ \left[\text{0.165,0.01}\right] $。

图  2  粗糙煤岩界面地质模型及正演模拟

Figure  2.  Geological model of coal-rock interface and forward modeling results

下载: 全尺寸图片 幻灯片

采用时域有限差分方法进行探地雷达正演模拟。以中心频率为900 MHz的布莱克曼-哈里斯脉冲为激励源，其时域波形如图2a所示。激励源位于空气层上方，自起始位置0处沿测线方向移动，以自激励和自接收模式运行。电磁波沿探测方向传播，引入Mur二阶吸收边界条件，边界层厚度为0.1 m。模拟区域网格大小为0.01 m，时间步长为0.01 ns，运行1500个时间步，总模拟时间为 15 ns。空间步长为0.01 m，共记录50道数据。正演模拟B-Scan图像如图2c所示。

图2b模拟了煤岩界面中常见的起伏、倾斜走向，及粗糙空气‒煤界面、粗糙煤岩界面。图2c中，空气‒煤界面及煤岩界面回波时延的变化规律与界面走向基本一致；同时，受粗糙界面影响产生了部分干扰双曲线杂波，与目标界面回波发生了干涉现象，在局部位置弱化甚至湮没了目标界面的响应特征。

2.   粗糙煤岩界面波组自动拾取方法

2.1   锚点提取及筛选

EDLine算法定义图像中梯度强度较大且变化显著的像素点为锚点，锚点的提出有效减少了线段绘制中的边缘点，是EDLines在计算效率上取得重要突破的原因之一^[19-20]。借鉴这一思想，将探地雷达数据中具有显著特征的关键点定义为锚点，即A-Scan数据中的波峰点或波谷点。

电磁波在煤层中传播并撞击介电常数更大的岩层时，电磁波会发生反相，导致入射波的波谷点在后向散射回波中变为波峰点，而以单一波峰点为锚点描述潜在目标界面特征时，会出现2层锚点描述同一煤岩界面的问题，增加了后续计算的复杂度。因此，考虑电磁波反相效应，针对布莱克曼-哈里斯脉冲，以界面回波波谷点为锚点。

局部极值方法通过遍历A-Scan序列数据，比较相邻数据点振幅的大小，可以准确提取极值锚点。但受非均质煤质影响，电磁波在其内部传播时会产生弱后向散射；同时，粗糙煤岩界面产生的电磁波波干涉效应也会发生波形部分重叠现象，在A-Scan图中表现为密集的波谷点。因此，有必要在局部极值方法中引入锚点振幅阈值及锚点间距阈值，以对锚点进行筛选。

设置锚点邻域与锚点间距阈值一致。探地雷达天线垂直分辨率指出了能够区分不同层介质的最小距离，以垂直分辨率为锚点间距阈值，可以减少电磁波干涉效应产生的无意义锚点。公式（6）提供了探地雷达天线垂直分辨率的经验公式：

式中：$ \lambda $为电磁波的波长，m。

采用局部极值方法提取锚点是基于A-Scan序列数据，而锚点位置的描述通常是以采样点为依据。为此，推导了利用采样点描述的锚点间距阈值计算公式：

式中：$ {R}_{\mathrm{T}} $为锚点间距阈值；$ {\varepsilon }_{\mathrm{r}} $为介质的相对介电常数；$ N $为采样点数，$ f $为天线中心频率，Hz；$ T $为时窗长度，ns。电磁波在煤岩层中传播$ {R}_{\mathrm{T}} $小于其在空气中传播，为充分保留潜在目标界面锚点，设置式（7）中$ {\varepsilon }_{\mathrm{r}} $为岩层介电常数$ {\varepsilon }_{\mathrm{r}\mathrm{R}} $。

计算候选锚点一维邻域内波点振幅的均值及标准差，并利用公式（8）计算锚点振幅阈值：

式中：$ {A}_{\mathrm{T}} $为锚点振幅阈值，V；$ {A}_{\mathrm{\mu }} $为邻域振幅均值，V；$ {A}_{\mathrm{\sigma }} $为邻域振幅标准差，V。

为解决引入$ {R}_{\mathrm{T}} $后目标界面锚点被舍弃的问题，当$ {R}_{\mathrm{T}} $内存在多个锚点时，保留振幅绝对值最大的锚点。此外，考虑煤岩界面锚点振幅要明显大于部分干扰杂波，为进一步减小无意义锚点数量，利用公式（9）保留目标锚点：

式中：$ {A}_{\mathrm{a}} $为目标锚点振幅，V；$ {A}_{\mathrm{a}\mathrm{\mu }} $为全域锚点振幅均值，V; $ k $为变量，取决于不同地层的煤岩电性参数和目标物埋深，此处取0.4。

采用上述方法提取潜在粗糙煤岩界面曲线的目标锚点，如图3所示。其中， $ N $、$ \varepsilon $、$ f $、$ T $为正演模拟预设参数，分别取值1500、6 F/m、$ \text{9×}{\text{10}}^{\text{8}} $ Hz、$ \text{15×}{\text{10}}^{\text{‒}\text{9}} $ s。 对比图3和图2c可以看出，提出的锚点自动提取方法基本保留了潜在目标界面曲线上的锚点，并弱化了其他干扰波形的响应特征。

图  3  锚点自动提取与筛选

Figure  3.  Automatic extraction and filtering results of anchor points

下载: 全尺寸图片 幻灯片

2.2   锚点区域生长

图4所示为电磁波在光滑及粗糙空气‒煤界面处的响应特征，可以看出电磁波在粗糙界面的后向散射特性较为复杂，同一界面锚点对应的采样点呈现不确定性分布。

图  4  平直/粗糙空气‒煤界面对后向散射波时延的影响

Figure  4.  Influence of rough interfaces on time delay of backscattered waves

下载: 全尺寸图片 幻灯片

为准确拟合煤岩界面曲线形态，设置以锚点为左边界中心点的右侧矩形区域，为锚点区域生长邻域。同时，考虑煤岩界面曲线走向较为平缓且具有一定不确定性，设置其邻域不存在特定方向偏好，即允许锚点在领域内均匀分布。锚点区域生长邻域边界由公式（10）确定：

式中：$ d $为邻域道数，通常取1。

粗糙煤岩界面起伏高度大于天线垂直分辨率时，会产生杂波，与目标界面回波反生干涉，导致锚点分布的连续性被破坏，即可能出现邻域内多锚点或虚假曲线锚点问题。杂波与目标回波发生干涉后，只在局部区域内产生波形的重叠，相位差的存在导致了目标回波振幅的结构性增强或破坏性减弱，但在干涉区域外，其振幅的变化保持原有规律，即具有一定的连续性。同时，相对于电磁波在粗糙表面的散射，电磁波干涉效应引起的锚点振幅会发生突跃式变化。因此，可以利用锚点振幅的波动特征作为锚点终止生长的判据。

变异系数常用于衡量数据的波动性，其数学表达式为

式中：$ \sigma $是数据集振幅标准差；$ \mu $是振幅均值。通常$ CV{\leqslant 0.2} $表示数据点变化较平稳和连续。

锚点区域生长采用邻域内就近生长方法，当邻域内无锚点时停止生长。利用式（11）处理多锚点问题，并对已生长的曲线进行锚点验证，分别选择曲线中包括任一锚点在内前3个锚点计算其振幅变异系数，如果$ CV\text{}\text{ > 0.2} $，则停止生长。锚点区域生长结果如图5所示。

图  5  锚点区域生长

Figure  5.  Schematic of anchor point region growth results

下载: 全尺寸图片 幻灯片

对比图5和图2b可以看出，锚点区域生长方法可以准确检测潜在目标界面曲线的锚点，同时振幅变异系数的引入可以有效克服电磁波干涉效应导致的虚假曲线锚点问题。

2.3   曲线分类识别

利用探地雷达识别煤岩界面时，B-Scan图像中干扰杂波通常由点目标散射导致，在B-scan表现为开口向下双曲线特征。但考虑粗糙界面对锚点分布特性的影响，可以利用开口向下的抛物线进行曲线锚点的拟合。

RANSAC算法是常用于处理非线性回归问题的迭代拟合方法，通过随机采样和一致性检验，可以有效寻找最优模型拟合参数^[24]。针对锚点区域生长后的共曲线锚点子集，当其锚点数量大于3时，任选3个锚点进行RANSAC迭代拟合，以决定系数$ {R}^{2} $为评价指标，进行曲线分类。同时，针对电磁波干涉效应导致完整曲线被切分问题，再次利用RANSAC迭代算法进行同类曲线的合并，其中，参与迭代的为共曲线上的一组锚点。拟合结果如图6所示。

图  6  拟合曲线分类识别

Figure  6.  Schematic of fitted curve classification and identification results

下载: 全尺寸图片 幻灯片

由图6可以看出，粗糙煤岩界面处产生了多条干扰双曲线。双曲线顶点位于引起电磁干涉效应的粗糙点位置，通常与煤岩界面曲线局部重合或相近。因此，针对开口向下且$ {\text{R}}^{\text{2}}{\geqslant 0.9} $的拟合抛物线，计算其顶点坐标，并以公式（10）锚点邻域边界为顶点邻域大小，考虑曲线拟合误差，$ d $取值为3。当邻域内存在其他拟合曲线时，认为该顶点所在曲线为干扰双曲线，可以直接去除。处理后结果如图7所示。

图  7  煤岩界面处干扰双曲线滤除结果

Figure  7.  Filtering results of interference hyperbolas at coal-rock interface

下载: 全尺寸图片 幻灯片

图7中，曲线分为潜在目标界面曲线及干扰双曲线，但结合图2b可以看出，当煤岩界面走向较为平缓时，部分潜在煤岩界面曲线被分类为干扰双曲线。利用波形特征匹配方法进行检测存在一定局限性，尤其是当煤岩界面曲线符合双曲线特征时，有必要针对干扰双曲线进行二次检测。

2.4   虚假检测控制

电磁波在介质中传播，其振幅变化主要由电磁波入射角度、色散、几何扩散、吸收和散射等效应所决定^[25]。当电磁波天线中心频率、煤岩层电性参数确定后，其影响变量为电磁波传播距离、入射角度和界面粗糙度，此时电磁波振幅表示为^[26]

式中：$ {A}_{0} $为初始振幅，V；$ z $为传播距离，m；$ \alpha $为介质吸收系数；$ \text{R}\left(\text{θ}\right) $为与入射角度相关的反射系数；$ \text{S}\left({\delta }_{\mathrm{h}}\right) $为与介质界面粗糙度相关的散射系数。

推导利用采样点描述的传播距离表达式：

式中：$ c $为电磁波在真空中传播速度，$ \text{3×}{\text{10}}^{\text{8}} $ m/s；$ T $为时窗长度，ns；$ N $为采样点数；$ \Delta N $为当前锚点与参考锚点的采样间隔；$ {\varepsilon }_{\mathrm{r}\mathrm{c}} $为煤质的相对介电常数。

电磁波在煤介质中传播的吸收系数$ \alpha $可以通过其电性参数来确定，通常与电导率、介电常数和磁导率等因素有关。对于煤等一般绝缘介质，吸收系数可以表示为

式中：$ \omega $为电磁波的角频率，rad/s，$ \omega =\text{2}\mathrm{\pi }f $；$ \mu $为煤介质的磁导率，通常与真空磁导率相近，1.257$ \times $10^‒6 H/m；$ \sigma $为煤介质的电导率，S/m。

对于垂直极化的电磁波，反射系数$ R\left(\theta \right) $可以表示为

式中：$ \theta $为电磁波入射角度。

Kirchhoff近似下，散射系数$ S\left({\delta }_{\mathrm{h}}\right) $可以表示为^[27]

式中：$ k $为电磁波波数，m^‒1，$ k=\text{2}\mathrm{\pi }/\lambda $。

对于连续分布的粗糙煤岩界面，小尺度范围内可认为电磁波垂直撞击分界面，此时其振幅简化为：

式中：$ {A}_{\mathrm{Z}\mathrm{\theta }} $为电磁波垂直入射煤岩界面振幅，V。

干扰双曲线由电磁波在粗糙煤岩界面处发生后向散射形成，可认为其对应同一均方根高度，即双曲线顶点位置处，其振幅为

式中：$ {A}_{\mathrm{Z}{\delta }_{\mathrm{h}}} $为粗糙界面同一均方根高度下振幅，V。

公式（17）和（18）中，$ {A}_{\mathrm{Z}\mathrm{\theta }} $和$ {A}_{\mathrm{Z}{\delta }_{\mathrm{h}}} $值还与电磁波的传播距离相关。利用式（12）—式（14），计算考虑电磁波几何扩散效应的锚点补偿振幅：

以同一参考采样点对双曲线锚点传播距离进行振幅补偿后，由公式（17）和（18）可以看出，粗糙煤岩界面锚点振幅变化与其均方根高度有关，而随机粗糙煤岩界面均方根高度通常呈现不确定变化；双曲线反射波振幅则与入射角度相关，随着探地雷达天线的移动，振幅先增大后减小，在顶点位置处达到最大值。利用式（19）对双曲线锚点振幅进行补偿，如图8所示。

图  8  考虑几何扩散效应的振幅补偿前后示意

Figure  8.  Compensated amplitude considering geometric diffusion effects

下载: 全尺寸图片 幻灯片

图8中，对比补偿前后干扰双曲线和粗糙煤岩界面曲线锚点振幅变化规律，干扰双曲线锚点振幅随着电磁波入射角度的增大而减小，粗糙煤岩界面锚点振幅随界面均方根高度的随机变化而无规律波动，其振幅变化规律与上述理论分析一致。

提取共抛物线锚点振幅值并采用滑动窗口方法计算每个窗口内的振幅斜率，窗口大小设置为2。当存在3个连续窗口振幅斜率同向时，定义该曲线为干扰双曲线，否则为煤岩界面曲线。同时，默认检测到的首个连续分布的界面为空气‒煤壁界面，可直接排除。煤岩界面曲线检测结果如图9所示。对比图2b和图9可以看出，该方法可以准确拾取具有复杂走向的粗糙煤岩界面位置。

图  9  随机走向粗糙煤岩界面曲线检测结果

Figure  9.  Detection results of randomly oriented rough coal-rock interface curves

下载: 全尺寸图片 幻灯片

3.   试验验证

利用GSSI SIR-20型号探地雷达（图10a）在山西某煤矿掘进工作面对顶板余煤厚度进行探测。试验地点为待掘进巷道开口处，煤层平均厚度为3.31 m，直接顶平均厚度为5.33 m的细砂岩/砂质泥岩。煤层倾向大致为北西向，煤岩层倾角平均为4°。收发一体天线中心频率为400 MHz，时窗长度为25 ns，采样点数为512，测线沿水平方向布置，共采集500道数据。

图  10  试验用探地雷达及B-Scan图像

Figure  10.  GPR and scanned B-Scan images for experiments

下载: 全尺寸图片 幻灯片

图10b为实测雷达数据，经降噪、去直流及滤波等处理后，煤岩界面反射波清晰可见，且存在多个干扰双曲线波形，与目标界面回波发生较强干涉。根据前序分析， 400 MHz天线下该煤岩界面为粗糙表面。

读取探地雷达数据，利用提出方法进行粗糙煤岩界面波组自动拾取。设置砂岩介电常数为6 F/m，其他参数采用雷达主机预设值和默认值。检测过程及结果如图11所示。

图  11  煤岩界面自动拾取过程及结果

Figure  11.  Coal-rock interface detection results

下载: 全尺寸图片 幻灯片

图11中，粗糙煤岩界面波组拾取位置与实际位置一致。该方法输入参数仅为探地雷达天线的中心频率、时窗长度、采样点数等探测参数及岩层的介电常数，可实现煤岩界面波组的自动拾取。

4.   结　　论

1）结合蒙特卡洛方法和样条插值方法构建了粗糙煤岩界面地质模型，并采用时域有限差分方法正演模拟了起伏走向和倾斜走向粗糙煤岩界面的电磁散射特性。

2）提出了粗糙煤岩界面探地雷达波组自动拾取方法，其变量参数仅为已知的探测参数和岩层介电常数，可实现煤岩界面的自动拾取。

3）提出了考虑电磁波几何扩散效应的振幅补偿方法，并分别构建了由粗糙表面均方根高度和电磁波入射角度主导的电磁波振幅衰减数学模型。

4）正演模拟及工程试验结果表明，该方法可准确检测粗糙煤岩界面位置。