左 杰，歐加利，劉金蓮

(1. 中交第三航務工程局有限公司江蘇分公司，江蘇 連云港 222042；2. 長沙理工大學 交通運輸學院， 湖南 長沙 410076)

0 引 言

探地雷達探測是一種有效的地球物理探測方法，已經廣泛地應用于巖土工程、環(huán)境和地質調查中，如地質災害、隧道地質超前預報探測，以及地下巖體結構特征探測等[1]。

正在建設的深茂鐵路江門至茂名段新臺隧道位于廣東省江門市新會區(qū)羅坑鎮(zhèn)境內，隧道全長6 258 m，隧道拱頂埋深0～260 m，總體走向為北東-南西向，隧道通過低山丘陵區(qū)。新臺隧道洞身穿越多條南北向斷層，隧道附近上方分布多個人工水庫，隧址區(qū)地下水主要為孔隙水、基巖裂隙水和構造裂隙水，區(qū)域斷層破碎帶構造裂隙水發(fā)育，為富水區(qū)。當隧道開挖至富水帶后易產生小規(guī)模地質災害，對施工造成影響。為防止隧道施工中的地質災害，保證施工安全和施工質量，在隧道開挖施工中采取了地質雷達方法地質超前預報，以探測、了解隧道施工掌子面前方地質情況。然而，由于隧道地質條件的復雜性和探測環(huán)境中的各種干擾，有用的探地雷達探測數(shù)據信息被各種噪聲、雜波淹埋，使得對探地雷達數(shù)據圖像結果分析、解譯困難[2]。

研究探地雷達回波記錄數(shù)據的子空間數(shù)據處理分析方法，提取強背景噪聲和雜波干擾下巖體結構信息。在數(shù)據處理中，將探地雷達記錄的原始回波數(shù)據分解成不同的數(shù)據子空間；在這些子空間中選取包含有多次波、噪聲、雜波等無用的干擾信息構成的子集，用一種簡單的信號近似方法，作為新的信號分量組合。關于子空間中建模方法，近年來已在信號和圖像處理分析與研究中得到成功的應用，如用奇異值分解方法進行離散有限長時間信號建模、分析處理數(shù)據，增強有用信號信噪比提取特征信息[3-4]。

奇異值分解數(shù)值分析與處理方法近年來已經在探地雷達信號處理中開始應用，如使用探地雷達探測地雷的信號處理分析中，針對地雷淺埋目標體特點，研究基于奇異值分解的子空間分解方法處理探地雷達剖面圖像數(shù)據，壓制無用的直達波，提取淺層掩埋體地雷特征信息，得到較好的結果[4]。還有一些其它的信號子空間分析方法，如基本量分析法(PCA)、獨立分量分析法(ICA)、經驗模態(tài)法(EMD)等，都已應用于探地雷達數(shù)據處理與分析研究中，這些方法壓制雜波和背景散射干擾波，取得了較好的研究成果。對于探地雷達記錄的剖面二維信號，研究表明可以較為有效地用曲波分析方法表示，這是因為探地雷達記錄剖面中有用的巖體結構信息往往以曲線方式表達，而在曲波空間中，這類特征信息可以用較少的曲波系數(shù)來表達。在曲波和小波變換分析中，無用的噪聲信號分量和有用的巖體結構信號分量分別有不同的統(tǒng)計特征，這樣噪聲和雜波可以在這些變換域中通過設置適當?shù)拈y值或門限除去。

多信號分類的超分辨率譜估計方法(MUSIC方法)是近年研究的一種子空間方法[5-6]，該方法是基于信號空間理論的準譜估計方法，是在基于高信噪比條件下相關矩陣特征矢量分解的多信號高分辨率成像方法。M.SUM等人在探地雷達路面結構評價應用中，開展了用MUSIC數(shù)值方法分析處理超薄路面材料結構層測量研究，模擬、實驗研究了復雜與重疊雷達反射波時間識別與估計提取的數(shù)值計算，較好的獲取了薄層路面結構信息[7]。用該方法處理分析探地雷達數(shù)據剖面信號，提取復雜環(huán)境強干擾信號下有用巖體結構信息，是一種新途徑。

一般地，可以將探地雷達探測到的目標回波數(shù)據作為子空間數(shù)據向量線性組合[8]，并設定用子空間基函數(shù)來構建這個低維空間。用這種子空間方法模擬探測目標的有效性還取決于子空間的選取。研究目的在于，應用基于子空間方法,更好地獲取在強干擾環(huán)境下探地雷達信號中相對較弱的巖體結構信息。

1 信號子空間模型及方法

1.1 信號子空間模型

在探地雷達地質探測中，通常事先在探測區(qū)域設定測線，雷達探測天線沿設定的測線移動，雷達天線發(fā)射電磁波并記錄回波信號，如圖1(a)。

圖1 探地雷達探測模型及數(shù)據分析Fig. 1 Model with fractures and voids under subsurface surveyed by a GPR and data analysis

探測中，探地雷達采集的記錄數(shù)據表示了事先設定的某一測線的掃描結果B-scan，如圖1(b)。而在該測線中某一特定探測點記錄的是探地雷達的1條掃描線回波信號A-scan，如圖1(c)。該掃描線可以反映該點地下地質情況，這里用函數(shù)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為總的掃描線記錄數(shù)。從圖1(b)可以看出,探測中記錄的巖體結構裂隙等特征信號被淹沒在較強的雜波和噪聲中。從圖1 (c)的一維掃描線記錄可以看出，探測的巖體裂隙反射信號被強噪聲、多次波和雜波等所干擾、掩埋，圖1中的第2層裂隙結構特征在雷達圖中很難看到。由此看來，對這類復雜剖面信號雷達圖的地質分析與解譯比較困難。

1.2 方 法

為簡單起見，設V=[v1,v2,v3,…,vm]為信號子空間。這里vi為子空間向量，且各向量之間是線性獨立的，m是設定的子空間數(shù)。這樣，探測雷達回波信號可以表示為加權子空間基向量和，如式(3)。

(3)

考慮一維雷達信號A(xi,t)，并將其分解為有用目標體、雜波和噪聲3個子空間。子空間分解意在保存信號子空間，而摒棄噪聲與雜波。與目標反射信號相比，雜波能量一般相對較小，但在實際探測中，如在多層巖體結構特征環(huán)境中，這類雷達記錄圖像信噪比較低，雜波信號會屏蔽或掩埋較深的或復雜的目標體信號。研究噪聲信號的平穩(wěn)性統(tǒng)計特點，進行用傅立葉子空間去噪處理，提高雷達信號圖像信噪比。在去雜波干擾處理中，使用空間平滑處理雜波相干信號，由此來壓制多次波干擾，提高弱特征信號提取能力。

模擬實驗中，設定探地雷達某條掃描線的反射系數(shù)為L(1,2,3,…L)，L

函數(shù)(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 現(xiàn)場實驗研究

正在建設的深茂鐵路江門至茂名段新臺隧道位于廣東省江門市新會區(qū)羅坑鎮(zhèn)境內，隧道通過的地層主要是第四系殘破積層(Q4el+dl)和燕山期花崗巖(γ53)。該隧道區(qū)花崗巖節(jié)理、裂隙整體不發(fā)育、局部地段較發(fā)育，隧道附近分布多個人工水庫，如青翠水庫、百步梯水庫、龍門水庫、長坑水庫等，部分裂隙與水庫相連通，受水庫水補給，富水性較強。在文中研究的探地雷達數(shù)據處理方法實驗數(shù)據采集中，選取新臺隧道DK148+644斷面掌子面作探地雷達超前地質預報數(shù)據采集試驗場地，使用GSSI GPR設備，100 MHz天線。如圖3是正在施工的新臺隧道,圖4是實際探測數(shù)據采集過程。

圖3 建設中的新臺隧道Fig. 3 Xintai tunnel in construction

圖4 現(xiàn)場探地雷達探測Fig. 4 GPR data acquisition at site

3 實驗數(shù)據及分析

用筆者的方法處理分析上述實際實驗數(shù)據，圖5(a)是記錄雷達不清晰的原始圖像數(shù)據，圖5(b)是不清晰的傅里葉濾波預處理結果，圖5(c)是圖5(b)中第400道掃描線一維記錄數(shù)據曲線。

圖5 雷達原始數(shù)據預處理分析Fig. 5 Analysis and preprocessing of GPR raw data profile

分析圖5雷達圖像，可以估計道隧道掌子面前方10～15 m處巖體結構異常。但僅僅從這個圖像解釋，不能準確地確定其具體位置和規(guī)模。圖5中原始雷達圖像數(shù)據信噪比較低，這表明雷達電磁場信號在地下傳播過程中能量衰減或散射吸收。對原始雷達圖像數(shù)據進行預處理，信號的幅度得以放大增強，但同時噪聲也得到某種程度的放大，有用的目標結構信息仍掩埋在較強的噪聲之中。處理中采用頻域子空間分解，進行帶通濾波，提高雷達圖像信號的信噪比，去除平穩(wěn)的隨機噪聲信號，保存目標特征雷達信號。

圖6是筆者的方法處理結果，從圖中可以看出巖體結構裂隙可以清楚的提取，且這些異常區(qū)邊界點位置也可以較準確的確定。

圖6 MUSIC子空間方法分析結果Fig. 6 Processed result with MUSIC subspace

圖7 現(xiàn)場實驗布設Fig. 7 Layout of field experiment

為驗證筆者方法的有效性，在實際的雷達探測采集區(qū)進行了鉆孔及孔內成像分析，鉆孔深度20 m，圖7是鉆孔布置和雷達測線圖。

圖8 孔內成像結果及對應的MUSIC曲線Fig. 8 Imaging borehole view and corresponding MUSIC curve

圖8是第5號鉆孔的孔內成像分析圖及相對應的巖體結構異常位置分析MUSIC曲線峰值點。分析圖8孔內成像結果及對應MUSIC子空間處理的一維雷達掃描線曲線?？梢钥闯觯琈USIC曲線在 5.6 m, 10.5～11.0 m和 15.5 m處分別有3個峰值估計點，可以判斷為疑似巖體結構異常點，而這些峰值點對應的孔內成像圖正確定了巖體結構裂隙位置。由此看出，對于復雜雷達圖像數(shù)據，用MUSIC算法處理分析，能較為有效地提取巖體結構特征信信息。

4 結 語

試驗研究了基于子空間的探地雷達圖像數(shù)據巖體結構特征分析與處理方法。對于強干擾、復雜環(huán)境下探地雷達圖像剖面數(shù)據，首先進行基于傅里葉子空間的預處理，提高圖像的信噪比；此后，用基于MUSIC算法的超分辨率譜估計處理分析雷達圖像數(shù)據，進而能準確地提取巖體結構特征位置和規(guī)模。模擬實驗及現(xiàn)場實驗數(shù)據處理分析表明，方法簡單、有效，為探地雷達記錄圖像數(shù)據的地質分析與解譯提供了新途徑。

參考文獻(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.