基于数据子空间的探地雷达图像岩体结构分析

左 杰，欧加利，刘金莲

(1. 中交第三航务工程局有限公司江苏分公司，江苏 连云港 222042；2. 长沙理工大学 交通运输学院， 湖南 长沙 410076)

0 引 言

探地雷达探测是一种有效的地球物理探测方法，已经广泛地应用于岩土工程、环境和地质调查中，如地质灾害、隧道地质超前预报探测，以及地下岩体结构特征探测等[1]。

正在建设的深茂铁路江门至茂名段新台隧道位于广东省江门市新会区罗坑镇境内，隧道全长6 258 m，隧道拱顶埋深0～260 m，总体走向为北东-南西向，隧道通过低山丘陵区。新台隧道洞身穿越多条南北向断层，隧道附近上方分布多个人工水库，隧址区地下水主要为孔隙水、基岩裂隙水和构造裂隙水，区域断层破碎带构造裂隙水发育，为富水区。当隧道开挖至富水带后易产生小规模地质灾害，对施工造成影响。为防止隧道施工中的地质灾害，保证施工安全和施工质量，在隧道开挖施工中采取了地质雷达方法地质超前预报，以探测、了解隧道施工掌子面前方地质情况。然而，由于隧道地质条件的复杂性和探测环境中的各种干扰，有用的探地雷达探测数据信息被各种噪声、杂波淹埋，使得对探地雷达数据图像结果分析、解译困难[2]。

研究探地雷达回波记录数据的子空间数据处理分析方法，提取强背景噪声和杂波干扰下岩体结构信息。在数据处理中，将探地雷达记录的原始回波数据分解成不同的数据子空间；在这些子空间中选取包含有多次波、噪声、杂波等无用的干扰信息构成的子集，用一种简单的信号近似方法，作为新的信号分量组合。关于子空间中建模方法，近年来已在信号和图像处理分析与研究中得到成功的应用，如用奇异值分解方法进行离散有限长时间信号建模、分析处理数据，增强有用信号信噪比提取特征信息[3-4]。

奇异值分解数值分析与处理方法近年来已经在探地雷达信号处理中开始应用，如使用探地雷达探测地雷的信号处理分析中，针对地雷浅埋目标体特点，研究基于奇异值分解的子空间分解方法处理探地雷达剖面图像数据，压制无用的直达波，提取浅层掩埋体地雷特征信息，得到较好的结果[4]。还有一些其它的信号子空间分析方法，如基本量分析法(PCA)、独立分量分析法(ICA)、经验模态法(EMD)等，都已应用于探地雷达数据处理与分析研究中，这些方法压制杂波和背景散射干扰波，取得了较好的研究成果。对于探地雷达记录的剖面二维信号，研究表明可以较为有效地用曲波分析方法表示，这是因为探地雷达记录剖面中有用的岩体结构信息往往以曲线方式表达，而在曲波空间中，这类特征信息可以用较少的曲波系数来表达。在曲波和小波变换分析中，无用的噪声信号分量和有用的岩体结构信号分量分别有不同的统计特征，这样噪声和杂波可以在这些变换域中通过设置适当的阀值或门限除去。

多信号分类的超分辨率谱估计方法(MUSIC方法)是近年研究的一种子空间方法[5-6]，该方法是基于信号空间理论的准谱估计方法，是在基于高信噪比条件下相关矩阵特征矢量分解的多信号高分辨率成像方法。M.SUM等人在探地雷达路面结构评价应用中，开展了用MUSIC数值方法分析处理超薄路面材料结构层测量研究，模拟、实验研究了复杂与重叠雷达反射波时间识别与估计提取的数值计算，较好的获取了薄层路面结构信息[7]。用该方法处理分析探地雷达数据剖面信号，提取复杂环境强干扰信号下有用岩体结构信息，是一种新途径。

一般地，可以将探地雷达探测到的目标回波数据作为子空间数据向量线性组合[8]，并设定用子空间基函数来构建这个低维空间。用这种子空间方法模拟探测目标的有效性还取决于子空间的选取。研究目的在于，应用基于子空间方法,更好地获取在强干扰环境下探地雷达信号中相对较弱的岩体结构信息。

1 信号子空间模型及方法

1.1 信号子空间模型

在探地雷达地质探测中，通常事先在探测区域设定测线，雷达探测天线沿设定的测线移动，雷达天线发射电磁波并记录回波信号，如图1(a)。

图1 探地雷达探测模型及数据分析Fig. 1 Model with fractures and voids under subsurface surveyed by a GPR and data analysis

探测中，探地雷达采集的记录数据表示了事先设定的某一测线的扫描结果B-scan，如图1(b)。而在该测线中某一特定探测点记录的是探地雷达的1条扫描线回波信号A-scan，如图1(c)。该扫描线可以反映该点地下地质情况，这里用函数A(xi,t)表示。考虑地下目标结构反射体、背景散射和噪声等，A(xi,t)可以表示为

A(xi,t)=si(t)+bi(t)+ni(t)

(1)

式中：si(t)为目标反射信号；bi(t)表示直达波；ni(t)表示噪声。1组A-scan 构成探测测线的雷达图剖面B-scan，用方程(2)表示：

B=[A(x1,t),A(x2,t),…A(xN,t)]

(2)

式中：B表示沿测线探测的记录剖面，N为总的扫描线记录数。从图1(b)可以看出,探测中记录的岩体结构裂隙等特征信号被淹没在较强的杂波和噪声中。从图1 (c)的一维扫描线记录可以看出，探测的岩体裂隙反射信号被强噪声、多次波和杂波等所干扰、掩埋，图1中的第2层裂隙结构特征在雷达图中很难看到。由此看来，对这类复杂剖面信号雷达图的地质分析与解译比较困难。

1.2 方 法

为简单起见，设V=[v1,v2,v3,…,vm]为信号子空间。这里vi为子空间向量，且各向量之间是线性独立的，m是设定的子空间数。这样，探测雷达回波信号可以表示为加权子空间基向量和，如式(3)。

(3)

考虑一维雷达信号A(xi,t)，并将其分解为有用目标体、杂波和噪声3个子空间。子空间分解意在保存信号子空间，而摒弃噪声与杂波。与目标反射信号相比，杂波能量一般相对较小，但在实际探测中，如在多层岩体结构特征环境中，这类雷达记录图像信噪比较低，杂波信号会屏蔽或掩埋较深的或复杂的目标体信号。研究噪声信号的平稳性统计特点，进行用傅立叶子空间去噪处理，提高雷达信号图像信噪比。在去杂波干扰处理中，使用空间平滑处理杂波相干信号，由此来压制多次波干扰，提高弱特征信号提取能力。

模拟实验中，设定探地雷达某条扫描线的反射系数为L(1,2,3,…L)，L

函数(Pmusic)峰值点位置来估计,如方程(4)。

(4)

式中：α(τ)目标反射信号延时向量，表示地下目标特征体的位置深度。EN为L×(L-k) 残余矩阵，包含有杂波和噪声(L-k)个特征向量 ，如图2是一维模拟结果。通过这种子空间分解方法，岩体结构、裂隙等特征信息可以被提取，即使对于较深埋的岩体结构特征弱回波信号，也有较好的结果。

从图2(a)中可以看出，含有2层结构特征的接收信号雷达回波图中包络有多次波、杂波和噪声，通过文中的子空间MUSIC方法，将雷达接收的各种信号分别投影到不同的子空间之中，这样抑制无用的干扰信号，而有效地提取有用的信息。在文中模拟实验中，L为4，模拟的第2层岩层结构回波信号被有效的提取，于此同时杂波和多次波几乎全部被除去。

图2 雷达图一维MUSIC分析Fig. 2 MUSIC spectral estimation of 1D radargram

2 现场实验研究

正在建设的深茂铁路江门至茂名段新台隧道位于广东省江门市新会区罗坑镇境内，隧道通过的地层主要是第四系残破积层(Q4el+dl)和燕山期花岗岩(γ53)。该隧道区花岗岩节理、裂隙整体不发育、局部地段较发育，隧道附近分布多个人工水库，如青翠水库、百步梯水库、龙门水库、长坑水库等，部分裂隙与水库相连通，受水库水补给，富水性较强。在文中研究的探地雷达数据处理方法实验数据采集中，选取新台隧道DK148+644断面掌子面作探地雷达超前地质预报数据采集试验场地，使用GSSI GPR设备，100 MHz天线。如图3是正在施工的新台隧道,图4是实际探测数据采集过程。

图3 建设中的新台隧道Fig. 3 Xintai tunnel in construction

图4 现场探地雷达探测Fig. 4 GPR data acquisition at site

3 实验数据及分析

用笔者的方法处理分析上述实际实验数据，图5(a)是记录雷达不清晰的原始图像数据，图5(b)是不清晰的傅里叶滤波预处理结果，图5(c)是图5(b)中第400道扫描线一维记录数据曲线。

图5 雷达原始数据预处理分析Fig. 5 Analysis and preprocessing of GPR raw data profile

分析图5雷达图像，可以估计道隧道掌子面前方10～15 m处岩体结构异常。但仅仅从这个图像解释，不能准确地确定其具体位置和规模。图5中原始雷达图像数据信噪比较低，这表明雷达电磁场信号在地下传播过程中能量衰减或散射吸收。对原始雷达图像数据进行预处理，信号的幅度得以放大增强，但同时噪声也得到某种程度的放大，有用的目标结构信息仍掩埋在较强的噪声之中。处理中采用频域子空间分解，进行带通滤波，提高雷达图像信号的信噪比，去除平稳的随机噪声信号，保存目标特征雷达信号。

图6是笔者的方法处理结果，从图中可以看出岩体结构裂隙可以清楚的提取，且这些异常区边界点位置也可以较准确的确定。

图6 MUSIC子空间方法分析结果Fig. 6 Processed result with MUSIC subspace

图7 现场实验布设Fig. 7 Layout of field experiment

为验证笔者方法的有效性，在实际的雷达探测采集区进行了钻孔及孔内成像分析，钻孔深度20 m，图7是钻孔布置和雷达测线图。

图8 孔内成像结果及对应的MUSIC曲线Fig. 8 Imaging borehole view and corresponding MUSIC curve

图8是第5号钻孔的孔内成像分析图及相对应的岩体结构异常位置分析MUSIC曲线峰值点。分析图8孔内成像结果及对应MUSIC子空间处理的一维雷达扫描线曲线。可以看出，MUSIC曲线在 5.6 m, 10.5～11.0 m和 15.5 m处分别有3个峰值估计点，可以判断为疑似岩体结构异常点，而这些峰值点对应的孔内成像图正确定了岩体结构裂隙位置。由此看出，对于复杂雷达图像数据，用MUSIC算法处理分析，能较为有效地提取岩体结构特征信信息。

4 结 语

试验研究了基于子空间的探地雷达图像数据岩体结构特征分析与处理方法。对于强干扰、复杂环境下探地雷达图像剖面数据，首先进行基于傅里叶子空间的预处理，提高图像的信噪比；此后，用基于MUSIC算法的超分辨率谱估计处理分析雷达图像数据，进而能准确地提取岩体结构特征位置和规模。模拟实验及现场实验数据处理分析表明，方法简单、有效，为探地雷达记录图像数据的地质分析与解译提供了新途径。

参考文献(References)：

[1] SHAO W, BOUZERDOUM A, PHUNG S L. Sparse representation of GPR traces with application to signal classification[J].IEEETransactionsonGeoscience&RemoteSensing, 2013, 51(7):3922-3930.

[2]VUKSANOVIC B, BOSTANUDIN C, HIDZIR H, et al. Discarding unwanted features from GPR images using 2DPCA and ICA techniques[J].InternationalJournalofInformationandElectronicsEngineering, 2013, 3(3) :317-323.

[3] TIVIVE F H C, BOUZERDOUM A, AMIN M G. An SVD-based approach for mitigating wall reflections in the wall radar imaging[J].RadarConference,2011,6(1):519-524.

[4] GUNATILAKA A H, BAERTLEIN B A. Subspace decomposition technique to improve GPR imaging of antipersonnel mines[J].ProceedlingsofSPIE-theInternationalSociety,2000,4038:24-28.

[5] SAKAMOTO T, TAKIAND H, SATO T. Experimental study on super-resolution techniques for high-speed UWB radar imaging of human bodies[J].PIERSOnline, 2009, 5(6): 596-600.

[6] SHRESTHA S. Signal processing of ground penetrating radar using spectral estimation techniques to estimate the position of buried targets[J].EURASIPJournalonAppliedSignalProcessing, 2003(12):1-12.

[7] SUN M, BASTARD C L, PINEL N, et al. Estimation of time delay and interface roughness by GPR using modified MUSIC[J].SignalProcessing, 2017, 132(C):272-283.

[8] HO K C, GADER P D, WILSON J N, et al. Subspace processing of GPR signals for vehicle-mounted landmine detection system[J].ProceedingsofSPIE-theInternationalSocietyforOpticalEngineering, 2008, 6953(1D):1-10.