工程科学与技术   2020, Vol. 52 Issue (6): 120-130
基于频谱能量分析的地质雷达探测图像判读
温世儒1, 杨晓华2, 郭元术3     
1. 江西理工大学 建筑与测绘工程学院,江西 赣州 341000;
2. 长安大学 公路学院,陕西 西安 710064;
3. 长安大学 信息工程学院,陕西 西安 710064
基金项目: 国家自然科学基金面上项目(41272285);江西省教育厅科学技术研究项目(GJJ170564)
摘要: 通过人工肉眼对地质雷达探测图像进行判读的方法易受判读人员主观性、经验性影响。为了规避这一不足,提出一种基于Counterlet等高变换和K-Means++均值聚类分析的频谱能量判读方法。以实际公路隧道为依托,经现场探测获取不良地质体的原始探测数据;采用IDSP(interactive detection-surveying prediction)探测数据分析软件生成原始图像,实施背景去除、滤波等时域、频域预处理以提高信噪比;基于子带分布系数采用Counterlet对预处理后的图像进行分解和重构,并采用K-Means++算法将重构后图像中的频率信息转化为颜色特征;利用MATLAB对颜色特征进行提取,并据此建立不良地质体颜色特征样本库,将原始探测图像与样本库进行匹配对比以实现自动判读。结果表明:采用Counterlet等高变换对多方向、多分辨率、多尺度的地质雷达图像进行分解与重构是可行的,曲线边缘逼近效果良好,重构后的图像无信息丢失;K-Means++算法能实现地质雷达灰度图像中能量—频率—色彩的转化,转化后的图像色彩突出、直观;频谱能量的匹配对比能较准确快速地实现自动判读及较好地规避个体主观性。
关键词: 公路隧道    地质雷达    预报    频谱能量    判读    
Interpretation Based on Frequency Spectrum Energy Analysis of Ground Penetrating Radar Detection Image
WEN Shiru1, YANG Xiaohua2, GUO Yuanshu3     
1. School of Architectural and Surveying & Mapping Eng., Jiangxi Univ. of Sci. and Technol., Ganzhou 341000, China;
2. School of Highway, Chang’an Univ., Xi’an 710064, China;
3. School of Info. Eng., Chang’an Univ., Xi’an 710064, China
Abstract: The interpretation method of GPR (ground penetrating radar) detection image by human eye is easily affected by the subjectivity and experience of interpreters. In order to avoid this shortcoming, a spectral energy interpretation method based on Counterlet contour transform and K-Means++ clustering analysis was proposed. Based on the actual highway tunnels, the original detection data of unfavorable geological bodies were obtained through on-site detection. IDSP (interactive detection-surveying prediction) detection data analysis software was used to generate the original image, and time domain and frequency domain preprocessing such as background removal and filtering were implemented to improve the signal-to-noise ratio. Based on subband distribution coefficient, Counterlet was used to decompose and reconstruct the preprocessed image and the K-Means++ algorithm was used to convert the frequency information in the reconstructed image into color features. The color features were extracted by MATLAB, and a sample library of color features of unfavorable geological bodies was established accordingly. The original detection image was matched and compared with the sample library to realize automatic interpretation. Facts and figures showed that it was feasible to decompose and reconstruct multi-directional, multi-resolution and multi-scale GPR images by using Counterlet contour transformation, the curve edge approximation effect was good and the reconstructed images had no information loss. K-Means++ algorithm can realize the conversion of energy to frequency to color in the gray-scale image of GPR, and the converted image has prominent and intuitive color. The matching and comparison of spectral energy can realize automatic interpretation more accurately and quickly and avoid individual subjectivity.
Key words: highway tunnel    ground penetrating radar    prediction    frequency spectrum energy    interpretation    

地质雷达后期判读是整个预报工作的重要组成部分,探测图像的判读特征是对不良地质体进行判读的关键依据,能否正确地对判读特征进行分析在很大程度上决定了能否准确地对不良地质体进行识别。当前,广泛采用的判读分析方法主要是通过对探测图像中的定性特征进行识别,进而作出相关的分析解答,这种判读分析方法长期以来为工程实践提供了较好的指导。

然而,这种判读分析方法对判读人员的主观实践经验具有很大的依赖性,不同的判读分析人员由于其实践经验和专业知识储备的不同,即使面对相同的地质雷达探测图像,其判读结果也容易出现偏差,甚至得出完全不一致的判读分析结果,这一缺陷显然限制了研究成果的有效推广与应用。因此,为了规避这一不足,寻求不以判读人员主观而改变的具有定量化特点的图像判读特征成为了地质雷达判读工作中需要迫切完成的重要工作。

根据地质雷达的探测原理可知,雷达天线接收到的电磁回波将通过数字–信号转换器转变为数字信号,最终以彩色或者灰度图像的形式进行显示[1-2]。由图像学可知,一副图像是由M行、N列个有限元素组成的,每个元素所在的行列都对应有定量化的颜色信息、纹理信息和形状信息,故一副图像的外在视觉特征可以最终被定量地描述为颜色特征、纹理特征和形状特征等数字特征,也说明反射回波的振幅、频率等波形特征具有与之相对应的数字特征。

可以看到,图像的颜色特征、纹理特征及形状特征具有定量化特点,不因判读人员的主观而改变。因此,若能区别于常规定性的判读对地质雷达探测图像的上述定量数字特征进行提取识别,那么对于规避常规的定性判读所存在的主观性缺陷是大有益处的。

为此,温世儒等[3]以隧道衬砌钢拱架探测图像中的典型双曲线特征为例,提出一种基于改进BP神经网络的判识方法对图像中的双曲线特征进行自动判识,经网络构建与训练,成功地识别了实例图像,并利用MATLAB子函数对灰度图像的颜色直方图特征进行了提取。刘宗辉等[4]基于属性分析,以电磁波中心频率这一定量参数的变化规律为控制因素,通过S变换和子波模拟对溶洞、破碎带实施了预报判读分析。李尧等[5]以相位、振幅和频率为指标,针对溶洞、断层、裂隙等不良地质体提出了复信号分析技术,通过GPRMAX正演模拟的方法对不良地质体进行了预报分析。高永涛等[6]通过隧道典型不良地质体信号与正常地质体信号在时域、频域及时频域3个维度的对比分析,确定了信号最大振幅幅值、最大振幅位置、信号能量、频谱熵、时频谱熵及IMF1分量最大振幅为不良地质体信号辨识的6个典型波形特征,再通过二分类模型实施了自动辨识。Liu等[7]运用改进的双正交小波构造了运用于地质雷达信号定量分析的最优双正交小波基,提出基于该小波基的地质雷达典型图像信号定量分析法(QAGBW法),并成功将该法应用于模拟信号和典型空洞实测信号的定量分析中。考虑到地质雷达回波易受噪声等外界干扰,Tzanis[8]采用曲线变换(CT)方法对反射回波实施增强处理,然后采用多尺度计算模型提取识别了典型裂缝、断层及不同层理的反射信号。

可以看到,诸如上述研究均提出了基于定量特征对不良地质体进行预报分析。BP神经网络具有良好的收敛性和较低的计算冗余度,对于诸如双曲线等典型特征具有良好的识别性[3],但是该网络对于边缘奇异段的逼近性和方向性较差,易导致边缘信息丢失,不适合非典型图像的识别。文献[4-6]虽然采用了定量特征分析,但是局限于反射回波的波形特征,并未涉及探测图像本身的定量特征。文献[7-8]均从回波信号处理入手,设计了相应的算法和提取模型对典型图像的信号特征进行提取,但提取得到的依然是回波的波形特征。由此可见,目前还存在的不足之处有两点:一是当前多数识别研究主要集中在典型图像,对于大多数非典型图像自动识别方面的研究还不足;二是多数研究局限于波形特征,对探测图像本身特征开展的研究依然不足。因此,寻找一种适合于地质雷达探测非典型图像,具有良好容错性和曲线逼近性的自动识别方法是十分必要的。

对此,以广西六寨—河池高速公路瑶寨隧道和江西赣县—兴国高速公路郑枫隧道为依托开展现场探测,获取不良地质体的原始图像,提出一种基于Counterlet等高变换和K-Means++的频谱能量特征分析方法对探测图像进行自动识别研究。经工程实例验证,该方法具有较好的信息保真度和较快的计算速度,能为地质雷达非典型探测图像的自动识别提供相关参考。

1 工程概况与现场探测 1.1 工程概况

以广西六寨—河池高速公路瑶寨隧道和江西赣县—兴国高速公路郑枫隧道为依托开展现场探测。瑶寨隧道已经建成,是分离式岩溶长隧道[9]。郑枫隧道是一座上、下分离式四车道高速公路中的隧道,全长778 m,左、右线洞门均采用端墙式洞门。

郑枫隧道隧址区地处构造剥蚀低山丘陵地貌,山势总体呈南北走向,支脉和冲沟呈树杈状,地形起伏大,地表植被发育。进口位于一条南北走向的冲沟中段山腰,出口位于一条“Y”字型的冲沟中段边坡中下部。进、出口山体坡度陡峭,沟底见有基岩出露,常年流水,水质纯净。根据实地调查及钻孔揭露,地层结构较简单,表层均为第四系全新统残坡积(Q4el+dl)粉质黏土,下伏基岩为震旦系(Z)千枚岩。

根据区域地质资料结合现场地质,隧址区地形稍有起伏,无明显隆起或塌陷,未发现大型地质构造痕迹,钻探也未揭露断层、皱褶等地质构造。地下水有松散孔隙水、基岩裂隙水两种类型,孔隙水主要赋存于上覆土层,水量贫乏,基岩裂隙水主要分布于基岩节理裂隙中。隧址区未发现明显不良地质及特殊性岩土,但沿线围岩的整体性和稳定性较差,以Ⅳ、Ⅴ级围岩为主,开挖掘进时需防止塌方、大变形破坏及加强地质雷达超前预报。

需要特别说明的是,瑶寨隧道属于岩溶隧道,围岩地质条件较差,共穿越了3条断层。Ⅳ、Ⅴ级围岩的占比达到了76%,在开挖过程中已揭露大量的空腔型溶洞、填充型溶洞、富水溶洞、破碎带、岩溶管道,时常发生塌方、落石、涌水、突泥灾害,是典型的烂洞子工程。为了降低施工风险、保障施工安全,引入地质雷达探测技术实施围岩地质预报。

1.2 现场探测

在现场地质调查的基础上,采用GSSI—3000型低功耗地质雷达配备收发一体屏蔽式天线开展现场连续线测和点测。掌子面上的测线、测点按照“网格型”布置,其中:测点沿测线布置,间距为50~100 cm;测线在纵横向的间距均为1~2 m,其数量视掌子面具体情况(如:是否有局部塌方、是否有金属体干扰、掌子面是否严重凹凸不平等)而定,但在纵横向均不少于2条,如图1所示。考虑到地质雷达容易受到外界噪声、锚杆等金属体的干扰,实际探测时将台车往洞口方向搬离掌子面5~10 m,同时清除锚杆、钻机、钢筋网等金属干扰物体。

图1 测线测点布置示意图 Fig. 1 Sketch of scan line and point

根据奈奎斯特定律,为了保证采样的有效性,连续线测时,天线在掌子面上的移动速度小于10 cm/s;点测时,各个测点之间的间距根据掌子面的平整度进行调整,范围为1.0~1.5 m[10]。部分探测参数如表1所示,其他参数根据系统默认取值。

表1 探测参数 Tab. 1 Detection parameters

根据《铁路工程物理勘探规程》(TB10013—2004),以及参考《公路隧道地质雷达检测技术规程》(DB35/T 957—2009)的规定,每座隧道的标定应不少于1处,每次实测不少于3点。取算术平均值作为最终的相对介电常数。为了提高探测质量,在每次实际探测时,均实施了现场标定:1)采用简易钻杆在掌子面上钻孔(掏槽眼、辅助眼和周边眼位置各钻取1孔,钻孔高度以便于操作为准),孔径为50 mm、孔深为40~80 cm不等(规程建议不少于15 cm),视围岩钻入难易程度而定;2)在钻孔位置进行连续触发探测,并根据探测剖面获取与钻孔深度相同位置处的电磁波双程走时;3)根据式(1)计算相对介电常数;4)取3个相对介电常数的算术平均值作为输入系统中的最终取值。

${\varepsilon _{\rm{r}}} = {\bigg(\frac{{0.3t}}{{2d}}\bigg)^2}$ (1)

式中:t为电磁波双程走时,ns;d为目标体的厚度,即钻孔深度,m。

2 方案与技术路线

地质雷达属于电磁勘探技术,探测时需按照主机电脑的设定向目标地层发射具有一定频率的入射电磁波。电磁波本质上属于交替变化的电场和磁场在空间中的传播,当遇到能使电磁场发生改变的介质时,电磁波将产生反射和折射[11]。一般而言,地层本身是由固相、液相和气相组成的三相体,三相组成分子均属带电粒子,能与电磁波产生褶积,从而改变电场和磁场,并据此产生反射与折射,最终被雷达系统所接收。然而,目标地层周围的金属体(台车、锚杆、钢筋网等)、管状电缆、外加电磁场等也会与电磁波产生褶积,由此导致雷达系统接收到的电磁回波不但来自于真实地层,而且还来自于外加干扰体。

干扰体造成的电磁回波会对判读分析造成不利影响,因而首先需要对原始探测图像进行滤波、背景去除、道间平衡、反褶积等时间域和频率域处理,以提高信噪比。在此基础上,采用LP滤波器对处理后的图像进行分解和重构以得到多方向、多尺度和多分辨率信息;然后,采用K-Means++算法对重构后图像中的子带分布系数进行聚类处理,以便将频率信息转化为颜色特征;其次,利用MATLAB对颜色特征进行提取,并据此建立不良地质体颜色特征样本库;最后,将原始探测图像与样本库进行匹配对比以实现自动判读。上述方案和总体技术路线如图2所示。

图2 技术路线 Fig. 2 Technical scheme

近年来,神经网络、小波变换、S变换等人工深度学习算法在地质雷达探测与检测后期解译判读分析方面得到了较为成熟的应用,实现了由常规定性判读向定量判读的转变。目前,国内外多数研究主要集中在波形与频谱特征,如幅值、回波主频、波速等的定量提取,较少对探测图像本身的图形定量特征进行研究。事实上,探测图像本身正是波形与频谱特征的图像化,作者正是基于这一点提出有别于波形与频谱特征的图像特征判读分析方法。

3 原始图像预处理

原始图像预处理的目的:一是为了消除或者弱化由于外界干扰造成的“假信号”以提高目标土层真实反射信号;二是为了消除相同目标体多次反射形成的多次波叠加,最终目的是为了提高信噪比。

3.1 滤波去噪

现场实际探测所得原始图像,由于杂波、多次波的影响,往往无法分辨真实的波形(图3,经开挖验证为一空腔型溶洞的探测图像)。此时,需要进行去噪处理。滤波去噪的方法包含两类:一类是背景去除;另一类是平滑和垂直滤波。通过去噪处理,可以良好地提高信噪比,利于判读分析。采用IDSP分析软件进行去噪时,主要的去噪参数如表2所示,叠加次数等其他参数取系统默认值。

表2 滤波去噪参数 Tab. 2 Filtering parameters

图3 未处理原始图像 Fig. 3 Raw detection image

表2可知,高通、低通参数是被保留电磁回波的频率上下限,只有介于该频率范围内的电磁回波才会被保留,该频率范围外的回波则会被滤除,在处理后的图像上不再显示,从而提高真实信号的反映度。

3.2 反褶积

反褶积处理的目的在于对电磁反射子波进行压缩,以及衰减多次波,从而在通道上仅留下反射系数以提高垂向分辨率[12]图4为对图3实施背景去除、滤波去噪和反褶积处理后得到的图像,可见杂波和多次波被滤除,真实反射信号变得清晰。图4中,白框所示为滤波处理后的预判位置,黑框所示为开挖揭露后的位置,可见二者基本一致。

图4 预处理后图像 Fig. 4 Detection image after preprocessing

在判读中之所以要进行预处理,其原因在于地质雷达接收系统属于宽频带接收装置,不但会接收真实围岩的反射信号,也会接收非真实围岩的反射信号,导致原始图像中包含了大量的非真实围岩反射信号,如多次波反射信号、金属体反射信号、手机干扰信号等。这些假信号掩盖了真实的反射信号,如果不滤除此类假信号,那么在后续K-means++聚类处理时,提取得到的是真假信号混合在一起的无用信号,这将导致后续的样本库无法被建立,甚至其匹配与自动识别也会无效。

滤波和反褶积预处理的优势在于消除假信号而留取真实信号,以保证判读分析不被假信号所干扰,提高判读分析的准确性。

4 图像变换 4.1 算法选择

一直以来,Fourier变换是图像信号处理的首选工具,在图像处理领域被广泛采用。该变换属于全信号处理技术,会对全时间域内的信号进行变换处理。但Fourier变换在局部信号复杂、分辨率和图像均匀性要求高的情况下,对于局部信号的时频分析,分辨率无法自动与相应信号的频率进行调整,其适应性和适用性较差[13]。为了弥补这一不足,小波变换被提出并相继得到了广泛应用。小波变换的显著优点是,在全局信号范围内能根据信号频率的高低自动进行分辨率调整,从而使得变换前后的图像分辨率更加均匀、连续。

对于地质雷达探测图像而言,其图像存在大量的局部化信号,从信号分辨率单因素来看,小波变换亦能满足变换前后的分辨率要求。然而,小波变换采用的基函数是缺乏方向性的正方形支撑域,这容易使得边缘逼近性差[14-15]。地质雷达图像的边缘恰恰不是连续的,而是带有方向的细小曲线,由此导致小波变换无法对地质雷达图像的边缘信息进行“稀疏”表示。

因此,对于地质雷达图像而言,为了保证全局信号及局部边缘信号均能得到连续逼近,需寻求一种具有良好多尺度、多分辨率和多方向信息的图像变换算法。

自然条件下的图像,其频率、分辨率和曲线连续性都是多方向、多尺度的。为了更好地对此类图像进行变换处理,Ridgelet变换和Curvelet变换被相继提出[16]。Ridgelet和Curvelet变换在方向性、尺度性和分辨率方面较小波变换有了极大地改善,但对于曲线边缘的逼近效果依然不佳。为此,Do M N和Martin Vetterli于2002年提出了Counterlet等高变换。

Counterlet等高变换采用长方形基结构。该结构的方向不是固定的,能随着图像曲线的方向不断进行调整,即方向自适应性,由此使得该结构在图像处理中具有良好的多方向性,适用于对多频率和多尺度图像信息进行处理变换,尤其适用于对带有方向的细小曲线或线段进行处理。Counterlet变换的逼近曲线详见文献[17]。

因此,选择Counterlet变换对经过预处理的地质雷达图像进行分解与重构处理。

4.2 图像分解与重构

除直线或者类直线段外,地质雷达图像中有许多细小的曲线段,具有典型的多方向、多尺度和多分辨率特征。采用Counterlet变换对地质雷达图像进行分解与重构的目的,就是为了得到探测图像的多方向、多尺度和多分辨率信息。

4.2.1 LP滤波器

拉普拉斯滤波器(LP滤波器)是Counterlet变换中所采用的一种重要滤波器,其基本作用是将图像中的曲线样条按照频率的高低进行分离。

地质雷达探测时,入射电磁波经与地层介质发生电磁褶积后,其反射能量(电磁强度)会产生改变,具有不同反射能量的子波对应于不同的频率。地层介质的破碎性、富水性、填充性等物理特征对反射回波的能量具有直接影响,相同或者相近的物理特征对应于相似的回波能量及其频率。因此,若能将地质雷达探测图像中的曲线样条按照频率的高低进行分离,就能对相近物理特征土层介质的反射回波进行提取分类,从而为实现自动识别奠定基础。因此,选用具有此功能的LP滤波器进行地质雷达图像的分解与重构。

LP滤波器的分解过程,可以简单地描述为:按照设定的迭代次数,重复不断地进行下采样和低通滤波处理,直至达到需要的分解尺度[18]图5为LP分解结构图(图中各符号仅表示含义,不是某个量值,没有单位),其分解过程可以描述为:先输入一个原始信号x;然后,通过低通滤波器H进行低通滤波处理;接着,采用采样矩阵M对低通滤波处理后的信号进行下采样,并由此得到低通信号c,此时再采用采样矩阵M对c进行上采样,并使用合成滤波器G进行合成滤波处理,进而得到信号P;最后,对Px进行差运算得到输出信号d

图5 LP分解结构图 Fig. 5 Decomposition structure chart of LP

对于地质雷达探测图像而言,x为经预处理后得到的图像,对该图像不断重复下采样和低通滤波,即可完成高频、低频的分离。

原始图像被分离成高频成分和低频成分,得到的是各方向、各尺度上相对应的高频信息子带和低频信息子带,为了对与子带相对应的轮廓信息进行提取,还必须得到各子带上的Counterlet系数[19-20]。为此,还须对分离图像进行LP重构运算。

LP重构是分解的逆过程,利用与分解算子相对偶的算子来实现,其算法在相关教材中可见,此处不再赘述,重构的基本程序如图6所示。

图6 LP重构结构图 Fig. 6 Reconstruction structure chart of LP

图6可见,LP重构过程为:对高通输入信号d做和运算得到信号P;然后,利用滤波器H对P进行低通滤波;接着,利用采样矩阵M对滤波后的信号进行下采样运算,将低通信号c与下采样得到的信号做和运算,并利用采样矩阵M进行上采样运算;最后,利用合成滤波器G进行合成,并与d做和运算而得到信号x

4.2.2 探测图像分解与重构

以瑶寨隧道YK46+980掌子面现场探测图像为实例进行LP分解与重构。图7为掌子面现场拍摄照片。

图7 YK46+980掌子面 Fig. 7 Working face of YK46+980

首先开展去噪预处理以抑制干扰波,并提高信噪比,继而进行分解与重构,处理过程如下:

1)图像预处理

图8为去噪预处理后得到的掌子面前方30 m范围内的灰度图像(纵向为探测深度,但是对于图像分解与重构而言,该坐标意义不大,因而省略,下同)。

图8 去噪灰度图像 Fig. 8 Filtering gray image

2)LP分解与重构

图5可知,理论上图像的分解级数(用nlevel表示)可以是无限大或小的,直至图像被无限缩小,频率信息完全丢失。但为了保证频率信息的完整,实际分解级数取值通常小于10[21-22]图9(a)(c)为nlevel=2、3、4时的分解图像,可见随着级数的提高,图像逐渐缩小,当nlevel=4时,部分高频子带信息开始丢失,说明nlevel的取值不能大于4。

图9 LP实例分解与重构 Fig. 9 Decomposition and reconstruction example of LP

为此,选择nlevel=3进行分解运算较为合理。图9(d)为重构后得到的图像,由图9(d)可见,LP重构后的图像与原图能保持高度一致,信息不会丢失。

5 K-Means++聚类与样本库建立 5.1 K-Means++聚类

LP分解与重构将地质雷达探测图像分离成高频和低频图像,分解与重构过程中各方向尺度上的Counterlet特征系数是反射回波频率及其能量的数字表达,其数据量较大,单纯地依靠特征系数在实际工作中不利于判读分析。一个可行的方法是:通过K-Means++聚类分析,对特征系数进行提取及实施聚类采集,并以类为单位赋予相应的颜色,进而得到直观、清晰的探测图像彩色特征图,不同的颜色代表了不同的反射频谱能量。

K-Means++聚类分析首先需确定聚类中心的个数[23]。根据瑶寨隧道及郑枫隧道现场探测与验证,地质体的类型共划分为5类,即:溶洞、富水围岩、Ⅲ级围岩、小破碎围岩、断层破碎带,其中,溶洞包括空腔型溶洞和填充型溶洞,因此聚类中心K的取值选5。然后,按照该聚类中心对特征系数进行聚类及赋色运算,从而把探测图像转化为彩色图像,实现“反射回波频率—特征系数—颜色”的转变,颜色的深浅表示反射能量的高低。

K-Means++聚类是图像处理领域的成熟算法,已被广泛采用,详细算法见文献[24]。

5.2 样本库建立

基于第5.1节分析,可以得到与围岩物理特征相对应的反射能量频谱,也就是彩色频谱图。据此,对实际探测所得,经工程验证后的5类地质体的探测图像进行上述去噪预处理、LP变换及K-Means++聚类运算获取彩色频谱图后,利用MATLAB建立样本库,基本步骤如下:首先,建立一个“*.mdl”格式的新窗口,并调用simulink函数建立目标库“mylib”,同时调用“setpath”函数设置存储路径;然后,打开“slblocks”文件模板,并利用“browser”函数显示将目标库命名为“DIZHILEIDAYANGBENKU”;点击命令“open”打开样本库,并用函数“show”进行显示。

通过上述步骤,对原始图像进行滤波、去噪预处理,然后采用LP滤波器进行分解与重构,以获取多方向、多尺度、多分辨率信息,接着利用K-means++进行聚类处理以转化、提取其频谱能量图像,最后只需将5类地质体的彩色频谱能量图像导入即可建立探测图像的彩色样本库,简单易于操作。对于后续探测所得图像,经去噪预处理—LP变换—K-Means++聚类后,调用imread颜色读取函数分别读取样本彩色图像与目标彩色图像的颜色分量,就可以完成图片的匹配比对,实现围岩物理特征的自动判读。

6 工程实例

以瑶寨隧道右线YK47+213及左线ZK47+190掌子面为实例段进行探测图像的频谱能量自动识别。

6.1 实例段地质概况与现场探测

右线YK47+213掌子面地质概况:埋深约52 m,岩性为白云质灰岩;围岩的整体完整性较好,掌子面右上侧局部较破碎,整体呈块状结构;岩面及周边围岩干燥无水,中层状构造,右上侧破碎部可见方解石及次生矿物填充;整体稳定性较好;根据地质勘察报告,该段围岩级别初步划分为Ⅳ级。左线ZK47+190掌子面地质概况:埋深约16 m,岩性为白云质灰岩;围岩完整性较差,掌子面较破碎,节理、裂隙发育,呈块状或块碎状结构;岩面及周边围岩干燥无水,薄或中层状构造;结构面由方解石及潮湿的次生矿物填充;整体稳定性较差;根据地质勘察报告,该段围岩级别初步划分为Ⅴ级。

为了保障施工安全,避免施工风险,按照现场探测的要求,施工时采用GSSI—3000型地质雷达开展超前地质预报以探明前方围岩的地质情况。图10为掌子面的现场探测灰度图像(已去噪处理)。

图10 现场探测图像 Fig. 10 Field detection image

6.2 LP变换与K-Means聚类

对现场探测图像进行LP变换处理,通过分解与重构对图像多尺度、多方向和多频率特征进行提取。图11为经变换处理后得到的特征提取图像。

图11 特征提取图像 Fig. 11 Image of feature extraction

经LP变换处理后,对特征提取图像进行K-Means++聚类处理运算,以便将图像中的频率、能量转化为颜色,从而实现灰度图像到直观彩色图像的转变。图12图11经聚类处理后得到的相应彩色图像。

图12 聚类处理图像 Fig. 12 Cluster image

对比图1112可见:1)灰度图像经K-Means++聚类处理运算后变成了直观的彩色图像,灰度图中不同反射强度的区域被赋予了不同的颜色;2)反射越强(反射能量越大)的区域,在彩色图中的颜色就越深,彩色图中的橙色区域对应于灰度图中的强反射区域,蓝色区域对应于灰度图中的弱反射区域;3)反射强度(反射能量大小)相近的区域,聚类处理时被赋予了相同的颜色。

6.3 自动识别与验证

图12(目标图像)导入至已建MATLAB样本库,调用MATLAB自带的imread颜色读取函数对目标图像和样本图像的颜色分量进行读取匹配,以便对掌子面前方围岩的地质情况进行预报。经匹配比对表明:1)图12与样本库中的12张Ⅲ级围岩的实测图像吻合,表明YK47+213前方围岩很可能是Ⅲ级围岩;2)图12与样本库中的17张Ⅳ级破碎围岩的实测图像相吻合,表明ZK47+190前方围岩很可能是破碎围岩。图13图12(b)自动识别后的窗口,可见系统自动给出了相应的围岩级别及安全提示。

图13 自动识别窗口 Fig. 13 Automatic identification window

现场探测及上一轮锚喷支护施工完毕并稳定后,施工单位随即进行了光面爆破作业,经爆破开挖发现:YK47+213前方围岩完整性较好,呈块状结构,中、厚层状构造,掌子面及周边围岩干燥,属于Ⅲ级稳定性围岩;ZK47+190前方围岩整体较破碎,结构面发育,尤其是掌子面左下侧部位岩体破碎,呈块状或块碎状结构;岩面及周边围岩较干燥,中、薄层状构造;结构面由方解石填充;整体稳定性较差,属于Ⅳ级或者Ⅳ级偏弱围岩。图14为爆破开挖后的掌子面,此时的里程分别为YK47+218.5、ZK47+194,可见基于频谱能量的样本库自动识别比对较为成功。在图14(b)中,可以看到台车,这似乎会对视觉效果造成影响,但该图片无法用其他图片替换,因为当时围岩较破碎,整体稳定性差,为了防止坍塌掉块造成人员损害,施工方及时将台车放置在掌子面前方,因此不论从哪个方向进行拍照,都无法避免台车。

图14 开挖后掌子面 Fig. 14 Working face of excavation

需要说明的是,样本库的匹配精度与样本库中围岩地质条件的类型数量及其容量有关,类型越多、容量越大,则精度越高。为了进一步提高匹配精度,后续还需进一步细化样本库中围岩地质条件的类型并增加其容量,这也是作者后续还需继续完成的重要工作。此外,准确获取(如:合理布置测线测点、合理设置探测参数、有效避免干扰信号等)原始探测图像是判读分析的前提[25],但获取原始探测图像属于现场探测的范畴,本文旨在分析一种能更好地规避判读个体主观性、更加快捷的自动判读方法,因而不再赘述现场探测的其他有关内容。

7 结 语

为了规避地质雷达判读分析时的个体主观性和经验性不利影响,作者基于理论分析、Counterlet变换和K-Means++聚类分析,以实际工程为依托分析探讨了一种利用频谱能量实施自动识别的判读分析方法,结论如下:

1)地质雷达的智能化判读能良好地规避判读个体的主观性,今后在地质雷达实践中应加以重视;

2)采用Counterlet等高变换对地质雷达图像进行分解与重构是可行的,曲线边缘逼近效果良好,能获得多方向、多分辨率、多尺度信息,重构后的图像无信息丢失,保真度良好;

3)K-Means++算法能将地质雷达灰度图像转化为彩色频谱能量图,转化后的图像色彩突出、直观,能量越大,颜色越深;

4)基于MATLAB样本库的匹配对比能较准确、快速地实现自动判读,可自动给出围岩级别和施工建议。

参考文献
[1]
Bai Dawei,Du Bingrui,Zhang Penghui. Weak signal processing technology of low frequency GPR based on Hilbert-Huang transform and its application in permafrost area gas hydrate exploration[J]. Geophysical and Geochemical Exploration, 2017, 41(6): 1060-1067. [白大为,杜炳锐,张鹏辉. 基于希尔伯特–黄变换的低频探地雷达弱信号处理技术及其在天然气水合物勘探中的应用[J]. 物探与化探, 2017, 41(6): 1060-1067. DOI:10.11720/wtyht.2017.6.10]
[2]
Li Jinghe,He Zhanxiang,Yang Jun,et al. Scale and rotation statistic-based self-adaptive function for ground penetrating radar denoising in curvelet domain[J]. Acta Physica Sinica, 2019, 68(9): 74-83. [李静和,何展翔,杨俊,等. 曲波域统计量自适应阈值探地雷达数据去噪技术[J]. 物理学报, 2019, 68(9): 74-83.]
[3]
Wen Shiru,Yang Xiaohua,Wu Xia. Research on intelligent identification and extraction of GPR image property based on BP neural network[J]. Highway, 2018, 63(7): 312-317. [温世儒,杨晓华,吴霞. 基于BP神经网络的探地雷达图像特征判识与提取研究[J]. 公路, 2018, 63(7): 312-317.]
[4]
Liu Zonghui,Liu Maomao,Zhou Dong,et al. Recognition method of typical anomalies in karst tunnel construction based on ground penetrating radar attribute analysis[J]. Rock and Soil Mechanics, 2019, 40(8): 3282-3290. [刘宗辉,刘毛毛,周东,等. 基于探地雷达属性分析的典型岩溶不良地质识别方法[J]. 岩土力学, 2019, 40(8): 3282-3290.]
[5]
Li Yao,Li Shucai,Liu Bin,et al. Forward simulation and complex signal analysis of borehole radar detection for underground adverse geological bodies[J]. Rock and Soil Mechanics, 2017, 38(1): 300-308. [李尧,李术才,刘斌,等. 钻孔雷达探测地下不良地质体的正演模拟及其复信号分析[J]. 岩土力学, 2017, 38(1): 300-308.]
[6]
Gao Yongtao,Xu Jun,Wu Shunchuan,et al. An intelligent identification method to detect tunnel defects based on the multidimensional analysis of GPR reflections[J]. Chinese Journal of Engineering, 2018, 40(3): 293-301. [高永涛,徐俊,吴顺川,等. 基于GPR反射波信号多维分析的隧道病害智能辨识[J]. 工程科学学报, 2018, 40(3): 293-301.]
[7]
Liu H R,Ling T H,Li D Y,et al. A quantitative analysis method for GPR signals based on optimal biorthogonal wavelet[J]. Journal of Central South University, 2018, 25(4): 879-891. DOI:10.1007/s11771-018-3791-y
[8]
Tzanis A. The curvelet transform in the analysis of 2-D GPR data:Signal enhancement and extraction of orientation-and-scale-dependent information[J]. Journal of Applied Geophysics, 2015, 11(5): 145-170.
[9]
Wen Shiru.Study on advanced detection technology of highway tunnel based on ground penetrating radar in karst area[D].Xi’an:Chang’an University,2015.
温世儒.基于地质雷达的岩溶地区公路隧道超前探测技术研究[D].西安:长安大学,2015.
[10]
Wen Shiru,Yang Xiaohua,Wu Xia,et al. Study on detection method of GPR for karst cave on tunnel side wall ahead of excavation face[J]. Journal of Highway and Transportation Research and Development, 2014, 31(10): 93-96. [温世儒,杨晓华,吴霞,等. 掌子面前方边墙溶洞的地质雷达探测方法研究[J]. 公路交通科技, 2014, 31(10): 93-96. DOI:10.3969/j.issn.1002-0268.2014.10.015]
[11]
Liu Zonghui,Wu Heng,Zhou Dong. Application of spectrum inversion method in GPR signal processing for tunnel lining detection[J]. Chinese Journal of Geotechnical Engineering, 2015, 37(4): 711-717. [刘宗辉,吴恒,周东. 频谱反演法在探地雷达隧道衬砌检测中的应用研究[J]. 岩土工程学报, 2015, 37(4): 711-717.]
[12]
Abbass M Y,Kim H W,Abdelwahab S A,et al. Image deconvolution using homomorphic technique[J]. Signal,Image and Video Processing, 2018, 13(4): 703-709.
[13]
Uruma K,Konishi K,Takahashi T,et al. Colorization-based image coding using graph Fourier transform[J]. Signal Processing(Image Communication), 2019, 74: 266-279. DOI:10.1016/j.image.2018.12.011
[14]
Nistane V,Harsha S,Sun Y,et al. Performance evaluation of bearing degradation based on stationary wavelet decomposition and extra trees regression[J]. World Journal of Engineering, 2018, 15(5): 646-658. DOI:10.1108/WJE-12-2017-0403
[15]
Chen Jun,Ge Shuangcheng,Zhao Yonghui,et al. Wavelet processing and its explanation on the GPR image of hidden seawall hazard[J]. Chinese Journal of Underground Space and Engineering, 2015, 11(Supp1): 337-341. [陈军,葛双成,赵永辉,等. 海堤隐患雷达探测图像的小波处理及解释[J]. 地下空间与工程学报, 2015, 11(增刊1): 337-341.]
[16]
Zhao B,Ren Y,Gao D K,et al. Fuzzy ridgelet neural network prediction model trained by improved particle swarm algorithm for maintenance decision of polypropylene plant[J]. Quality and Reliability Engineering International, 2019, 35(4): 1231-1244. DOI:10.1002/qre.2456
[17]
王润秋.信号分析与处理[M].北京:石油工业出版社,2015.
[18]
Zhang Lijuan,Huang Jianwu,Xu Xuejun. Prediction of time-dependent one-dimension temperature field of pavement on the basis of laplace transform[J]. Journal of South China University of Technology (Natural Science Edition), 2017, 45(11): 10-16. [张丽娟,黄建武,许薛军. 基于拉普拉斯变换的路面一维时变温度场预测[J]. 华南理工大学学报(自然科学版), 2017, 45(11): 10-16.]
[19]
Li Jianing.The research and application of counterlet transform in tunnel measured ground penetrating radar images processing[D].Xi’an:Chang’an University,2014.
李佳宁.Contourlet变换在隧道实测探地雷达图像处理中的应用研究[D].西安:长安大学,2014.
[20]
Wang Xu,Wang Guozhong,Fan Tao. Block-based adaptive compressed sampling of depth image[J]. Application Research of Computers, 2016, 33(3): 903-906. [王旭,王国中,范涛. 深度图像的分块自适应压缩感知[J]. 计算机应用研究, 2016, 33(3): 903-906. DOI:10.3969/j.issn.1001-3695.2016.03.060]
[21]
Wu Yiquan,Meng Tianliang,Wu Shihua,et al.Remote sensing images segmentation of rivers based on two-dimensional reciprocal gray entropy[J].Journal of Huazhong University of Science & Technology (Natural Science Edition),2014,42(12):70–74.
吴一全,孟天亮,吴诗婳,等.基于二维倒数灰度熵的河流遥感图像分割[J].华中科技大学学报(自然科学版),2014,42(12):70–74.
[22]
Rafael C,Richard E.数字图像处理[M].2版.阮秋琦,阮宇智,译.北京:电子工业出版社.2017.
[23]
Tang T L,Chen S Y,Zhao M,et al. Very large-scale data classification based on K-means clustering and multi-kernel SVM[J]. Soft Computing, 2019, 23(11): 3793-3801. DOI:10.1007/s00500-018-3041-0
[24]
Wang Haiyan,Cui Wenchao,Xu Peidi,et al. An optimized k-means++ algorithm guided by local probability[J]. Journal of Jilin University (Science Edition), 2019, 57(6): 1431-1436. [王海燕,崔文超,许佩迪,等. 一种局部概率引导的优化K-means++算法[J]. 吉林大学学报(理学版), 2019, 57(6): 1431-1436.]
[25]
Shi Gang,Xie Yongli,Yang Xiaohua,et al. Effect of GPR signal processing based on least square deconvolution[J]. Journal of Chang’an University (Natural Science Edition), 2014, 34(4): 104-108. [石刚,谢永利,杨晓华,等. 探地雷达信号的最小平方反褶积处理效果[J]. 长安大学学报(自然科学版), 2014, 34(4): 104-108.]