周東, 劉毛毛, 劉宗輝, 劉保東

(1.廣西大學(xué) 土木建筑工程學(xué)院，廣西 南寧 530004； 2.廣西大學(xué) 廣西防災(zāi)減災(zāi)與工程安全重點實驗室，廣西 南寧 530004； 3.南寧城建管廊建設(shè)投資有限公司，廣西 南寧 530219)

0 引言

探地雷達(dá)作為一種快速無損的地下目標(biāo)檢測技術(shù)，近年來已經(jīng)在路面病害定位和層厚檢測中得到廣泛應(yīng)用[1-4]。然而，目前探地雷達(dá)路面層位提取主要依靠解釋人員的經(jīng)驗或相關(guān)算法，人工或半自動拾取雷達(dá)剖面中強(qiáng)振幅同相軸連續(xù)的層位信息，存在主觀性強(qiáng)、解釋周期長、工作量大和每次僅能追蹤一個層位等缺點。因此，有必要研究一種路面多層位自動追蹤方法。

目前探地雷達(dá)層位識別和標(biāo)定主要是借鑒發(fā)展較為成熟的地震層位追蹤方法，基于波形相似特征[5]、人工神經(jīng)網(wǎng)絡(luò)[6]或圖像特征[7-8]自動拾取層位信息。由于地震波與電磁波在理論和數(shù)據(jù)處理中存在諸多相似性，已有許多學(xué)者將地震層位追蹤方法引入至探地雷達(dá)領(lǐng)域[9-12]。在路面檢測方面，Lahouar等[13]基于振幅閾值來檢測探地雷達(dá)強(qiáng)反射信號，并使用合成的子波信號與反射信號進(jìn)行相似計算來識別層反射數(shù)據(jù)。Loizos等[14]采用不同的介電常數(shù)估算方法對比分析了路面瀝青層厚度估算的準(zhǔn)確性。周輝林等[15]通過層界面檢測和模式識別等技術(shù)提取了探地雷達(dá)路基層厚特征，并在此基礎(chǔ)上研究了基于SVM的路基病害自動檢測算法。Le等[16]同樣基于SVM研究了探地雷達(dá)道路基層厚度的估計方法。Zhao等[17]研究了基于正則化反褶積技術(shù)的探地雷達(dá)瀝青層厚度預(yù)測算法。雖然已有許多學(xué)者采用不同方法對路面厚度和病害檢測展開了研究，并取得了一定進(jìn)展，但目前方法主要集中在單個層反射界面識別，且大多未考慮地下層界面同相軸反射的極性信息。

基于此，本文提出了一種快速、準(zhǔn)確的公路路面探地雷達(dá)多層位自動檢測方法。首先利用復(fù)信號分析獲取了探地雷達(dá)數(shù)據(jù)的瞬時相位余弦，并將其與合成的子波余弦矩陣數(shù)據(jù)進(jìn)行相似度分析后再計算瞬時相位余弦，提高相位同相軸的橫向連續(xù)性。根據(jù)深度方向相鄰層位線的振幅均方根平方值來確定層位數(shù)據(jù)，并通過設(shè)置層位線閾值和幅值閾值來過濾無關(guān)或不重要的層位數(shù)據(jù)。最后通過室內(nèi)數(shù)值模擬和現(xiàn)場探地雷達(dá)數(shù)據(jù)對本文方法的有效性和適應(yīng)性進(jìn)行了驗證。

1 方法理論

1.1 瞬時相位余弦

探地雷達(dá)采集得到的實信號可以表示如下：

y(t)=A(t)cos[ω0t+φ(t)] ,

(1)

式中：A(t)是關(guān)于時間t的函數(shù)，為瞬時振幅，主要與地下介質(zhì)電磁衰減特征和增益方法等有關(guān)；cos[ω0t+φ(t)]為瞬時相位余弦；[ω0t+φ(t)]為瞬時相位；ω0為中心頻率；φ(t)為關(guān)于時間t的相位函數(shù)。

為獲取探地雷達(dá)數(shù)據(jù)的瞬時相位余弦信息，需要將探地雷達(dá)數(shù)據(jù)進(jìn)行Hilbert變換后構(gòu)造復(fù)信號，在復(fù)信號分析中利用瞬時振幅和實信號來獲取瞬時相位余弦。實信號經(jīng)Hilbert變換后可表示為

(2)

由式(1)、(2)可以得到瞬時振幅：

(3)

瞬時相位余弦：

(4)

式中，實信號瞬時振幅不能為0，雖然在實際信號分析中這種情況很少，但若存在這種情況時瞬時相位余弦值應(yīng)置為0。

圖1顯示了一個歸一化的雷達(dá)子波波形及其瞬時相位余弦波形。圖中可以看出瞬時相位余弦波形的波峰、波谷和零值點位置與子波具有較好的對應(yīng)關(guān)系。由于該子波最大振幅絕對值在波峰處，因此在本文中認(rèn)為該子波的相位極性為正相位。

圖1 雷達(dá)子波及其瞬時相位余弦

1.2 瞬時相位余弦橫向增強(qiáng)處理方法

探地雷達(dá)在實際工程應(yīng)用中，受現(xiàn)場探測環(huán)境和天線自身耦合干擾，采集得到的數(shù)據(jù)中往往會包含有衍射干擾。衍射干擾會使得同相軸局部不連續(xù)，給基于相位特征的層位自動追蹤造成困難。由于公路路面的探地雷達(dá)數(shù)據(jù)主要是橫向?qū)臃植?，因此在進(jìn)行層位識別前，增強(qiáng)相位橫向連續(xù)性是有必要的。

地下層分布介質(zhì)的雷達(dá)反射波形在測線方向上具有連續(xù)性、漸變性和相似性等特點，而且反射波波形會保持對稱的Ricker子波[13]。由于瞬時相位余弦相當(dāng)于將雷達(dá)信號在波峰和波谷處進(jìn)行歸一化，降低了層反射數(shù)據(jù)的波形偏差，使得層反射數(shù)據(jù)與Ricker子波的瞬時相位余弦具有較高的相似度。若將子波瞬時相位余弦數(shù)據(jù)在測線方向上進(jìn)行拓展，并將其與雷達(dá)反射數(shù)據(jù)進(jìn)行相關(guān)性分析，則橫向?qū)臃植嫉姆瓷鋽?shù)據(jù)會具有較高的相似性。若在此基礎(chǔ)上重新計算相關(guān)分析后的瞬時相位余弦，則能夠進(jìn)一步降低測線方向上的波形偏差，增強(qiáng)橫向同相軸的連續(xù)性。兩個相同大小數(shù)據(jù)的相關(guān)系數(shù)可通過

(5)

為增強(qiáng)余弦剖面同相軸的橫向連續(xù)性，本文提出了一種瞬時相位余弦橫向增強(qiáng)處理方法。該方法在獲取探地雷達(dá)數(shù)據(jù)的瞬時相位余弦后，主要處理步驟如下：

1)構(gòu)建標(biāo)準(zhǔn)核矩陣。設(shè)置一個大小為m×n的標(biāo)準(zhǔn)子波瞬時相位余弦矩陣作為標(biāo)準(zhǔn)核矩陣，其中m為波長大??；n為測線方向上拓展的道數(shù)，值越大，橫向增強(qiáng)能力越強(qiáng)，但垂直分辨率越低。

2)計算相關(guān)系數(shù)矩陣。對于瞬時相位余弦數(shù)據(jù)中某樣點，選取以該樣點為中心、大小為m×n的窗口數(shù)據(jù)，計算該窗口數(shù)據(jù)與標(biāo)準(zhǔn)核矩陣的相關(guān)系數(shù)值并代替該樣點。

3)再次計算瞬時相位余弦。計算相關(guān)系數(shù)矩陣數(shù)據(jù)的瞬時相位余弦，即可得到橫向增強(qiáng)處理后的結(jié)果。

1.3 層位自動識別方法

為減少解釋人員主觀判斷的影響，提高公路路面探地雷達(dá)層位識別精度，文中提出了一種基于瞬時相位余弦的探地雷達(dá)層位自動識別方法。對于一個大小為M×N的雷達(dá)數(shù)據(jù)，首先計算其瞬時相位余弦，獲取各道瞬時相位余弦的極值點及空間位置信息。將獲取的相位極值點采用二值化的方法保存在矩陣C中(大小為M×N)，0表示該點瞬時相位余弦為非極值點，1表示極值點(波峰或波谷)。

為方便后續(xù)層位追蹤和極性識別，在相位極值數(shù)據(jù)中將具有相同極性的點進(jìn)行連線。具體方法是：以相位極值點為中心，設(shè)置深度方向搜索窗口大小，追蹤測線方向相鄰道具有相同極性的層位極值點，并將追蹤得到的相位極值點連接為層位線。為防止層位線追蹤出現(xiàn)串層和偏離，僅追蹤搜索窗口內(nèi)只有單個極值點時的數(shù)據(jù)。根據(jù)層反射單個波長內(nèi)相鄰相位極值點距離約為1/4波長的特點，搜索窗口大小一般設(shè)置為1/8～1/4波長。

在本文中，基于反射波的波形特點，認(rèn)為層位點在反射波的絕對振幅最大值處，并且該層位點的極性取決于絕對振幅最大值處的瞬時相位余弦值。實際層位線提取需要設(shè)置一個大小為1～1.5倍波長的深度搜索窗口，若某一層位線在深度搜索窗口內(nèi)存在上一條和下一條層位線，并且該層位線上的振幅均方根平方值大于另外兩條層位線上的振幅均方根平方值，則保留該層位線數(shù)據(jù)。深度搜索窗口設(shè)置為1～1.5倍波長主要是為了防止出現(xiàn)層位追蹤串層，并且盡量包含完整的子波反射特征。

為減少一些無關(guān)(不重要)的同相軸反射數(shù)據(jù)，可以通過設(shè)置層位線閾值Lmin過濾一些長度較短的層位線數(shù)據(jù)。對于層分布較廣或密集采樣的雷達(dá)數(shù)據(jù)，層位線閾值應(yīng)適當(dāng)增大。此外，在實際工程應(yīng)用中，由于層反射分界面兩側(cè)存在介電差異，層反射信號幅值相對較大，而背景或噪聲數(shù)據(jù)的信號幅值相對較小。因此，可以通過設(shè)置振幅閾值A(chǔ)min來進(jìn)一步去除幅值較小的層位線。經(jīng)過上述步驟，能夠追蹤得到強(qiáng)反射區(qū)域?qū)游痪€數(shù)據(jù)及其極性。

對于輸入的探地雷達(dá)數(shù)據(jù)，本文層位自動追蹤方法主要步驟如下：

1)獲取瞬時相位余弦信息。計算雷達(dá)數(shù)據(jù)的瞬時相位余弦，并進(jìn)行橫向增強(qiáng)處理，記錄瞬時相位余弦極值點對應(yīng)的空間位置、振幅和極性信息。

2)層位追蹤。設(shè)置搜索窗口大小，追蹤測線方向上相鄰道相同極性的數(shù)據(jù)，并連接為層位線。

3)層位識別。設(shè)置深度搜索窗口大小，若某層位線在深度窗口內(nèi)存在上一條和下一條層位線，并且該層位線上的振幅均方根平方值大于相鄰的另外兩條層位線，則認(rèn)為該層位線為有效層位線，保留該層位線。

4)層位線過濾。根據(jù)設(shè)置的層位線閾值和振幅閾值去除較短和幅值較弱的層位線。

5)極性識別。層位線對應(yīng)的瞬時相位余弦值大于0為正相位，反之為負(fù)相位。

2 數(shù)值模擬

為驗證本文方法的可靠性，利用GPRMAX 3.0[18]軟件模擬了一個地下多界面反射數(shù)據(jù)，幾何模型如圖2所示。模型大小為13.0 m×5.0 m，網(wǎng)格為0.005 m×0.005 m。幾何模型中主要包括3層反射介質(zhì)，在第二層中包含了2個具有相同層厚的傾斜介質(zhì)，各層介質(zhì)均為各向同性、不導(dǎo)電、非磁性的材料，從上到下各層介質(zhì)的介電常數(shù)分別為6、4、8、6。雷達(dá)天線從模型上方0.5 m的空氣層底部由左至右水平移動，發(fā)射和接收天線間隔0.4 m，激勵源為 200 MHz 的Ricker子波。采樣時窗為100 ns，每隔 0.1 m采集一道數(shù)據(jù)，共采集了130道數(shù)據(jù)。此外，為避免邊界效應(yīng)，實際模型區(qū)域在各邊界處向外延伸了2 m。

圖2 數(shù)值模擬幾何模型

圖3a為純凈數(shù)值模擬雷達(dá)剖面，可以看出圖中除了由反射面產(chǎn)生的同相軸反射外，還存在一些幅值較弱的偽影反射。為使模擬結(jié)果更接近實際情況，在數(shù)值模擬結(jié)果中添加了信噪比為15 dB的高斯白噪聲。在添加高斯噪聲后，偽影反射基本被遮蔽，但大致能看出幾何模型的邊界分布范圍(圖3b)。由于第二層右側(cè)介質(zhì)的介電常數(shù)高于模型其他介質(zhì)，電磁波在該介質(zhì)中傳播時速度較慢，因此在雷達(dá)剖面中該區(qū)域上下反射界面更遠(yuǎn)。

a—純凈無噪聲數(shù)據(jù);b—添加高斯噪聲

圖4為圖3b含噪數(shù)據(jù)的瞬時相位余弦及其橫向增強(qiáng)處理后的結(jié)果。由圖可見，噪聲的存在同樣使得瞬時相位余弦剖面中的同相軸反射大量被遮蔽，不利于基于同相軸的層位追蹤(圖4a)。而橫向增強(qiáng)處理后，瞬時相位余弦剖面中同相軸的連續(xù)性得到提高，同相軸反射特征得到顯著增強(qiáng)，有利于進(jìn)一步的層位追蹤(圖4b)。

a—橫向增強(qiáng)處理前;b—橫向增強(qiáng)處理后

利用本文方法追蹤了圖3b含噪數(shù)據(jù)的層位，結(jié)果如圖5所示。層位追蹤結(jié)果與模型分界面反射基本一致，其中黃、紅兩種顏色分別表示該層位為正、負(fù)相位極性。各界面處反射波的相位極性取決于分界面的反射系數(shù)，如第一個分界面ε1>ε2，此時反射系數(shù)為正，相位極性與入射波一致。而第二個分界面ε1<ε2，反射系數(shù)為負(fù)，相位極性與入射波相反。數(shù)值模擬試驗表明本文層位追蹤方法不僅能夠準(zhǔn)確自動識別出多個層位數(shù)據(jù)，還可以識別出層位的極性。

圖5 層位追蹤結(jié)果

3 現(xiàn)場案例

選取廣西北海市城市公路路面探地雷達(dá)數(shù)據(jù)進(jìn)行測試，該公路結(jié)構(gòu)包括瀝青面層、混凝土基層和水泥穩(wěn)定碎石層。探地雷達(dá)采集儀器選用意大利IDS公司生產(chǎn)的K2型探地雷達(dá)，天線中心頻率為900 MHz?，F(xiàn)場采集設(shè)置的時窗為60 ns，每道采集512個樣點，每米采集20道。為更好顯示路基探地雷達(dá)特征，本文后續(xù)僅展示前256個樣點(30 ns)數(shù)據(jù)。

圖6為原始探地雷達(dá)剖面及其測線第15 m處的單道波形。從圖中可以看出該區(qū)域主要存在2個相對平行的層反射，第一層和第二層反射波形狀均與Ricker子波類似，其中第一層反射在10～30樣點之間，第二層反射在60～80樣點之間。此外，單道波形圖中還可以看出上下兩層反射信號的最大絕對振幅分別在波谷和波峰處，可知上下兩層反射信號的極性相反。這是因為當(dāng)電磁波從瀝青面層(相對介電常數(shù)為3～5)傳播至混凝土基層(相對介電常數(shù)為6～9)時，反射系數(shù)為負(fù)，反射波的相位與入射波相反。而電磁波由混凝土基層傳播至水泥穩(wěn)定碎石層(相對介電常數(shù)為4～6)時，反射系數(shù)為正，此時反射波的相位與入射波相同，因此第一層與第二層的同相軸反射相位極性出現(xiàn)反轉(zhuǎn)。

圖6 原始雷達(dá)剖面和測線第15 m處單道波形

為驗證本文瞬時相位余弦橫向增強(qiáng)處理方法的有效性，分別追蹤了原始雷達(dá)數(shù)據(jù)瞬時相位余弦橫向增強(qiáng)處理前、后的層位，結(jié)果如圖7所示。圖中可看出未經(jīng)橫向增強(qiáng)處理的層位追蹤結(jié)果存在局部不連續(xù)現(xiàn)象，而經(jīng)橫向增強(qiáng)處理后第一層和第二層的層位線更加連續(xù)、完整，表明瞬時相位余弦橫向增強(qiáng)處理有助于提高層位線追蹤的連續(xù)性。此外，從圖中還可以看出第二層層位線起伏較大，在其下方還存在一條較短的正相位層位線。表明該區(qū)域混凝土基層存在局部沉降現(xiàn)象，根據(jù)相位極性可知下方較短的層反射界面下方介質(zhì)的介電常數(shù)小于水泥穩(wěn)定碎石層，并且該異常體反射信號幅值較強(qiáng)，因此推斷該異常反射體可能為局部脫空。經(jīng)室內(nèi)測試，本文方法層位點拾取窗口大小設(shè)置為5，層位線長度閾值為150，幅值閾值為原始雷達(dá)數(shù)據(jù)的絕對平均幅值。

圖7 橫向增強(qiáng)處理前(a)后(b)層位追蹤結(jié)果對比

為進(jìn)一步檢驗本文層位追蹤方法的有效性，將其與目前常用的基于波形相似特征的層位自動追蹤方法進(jìn)行對比，結(jié)果如圖8所示。圖中青色星號為傳統(tǒng)基于波形相似方法層位追蹤的種子點，第一層和第二層種子點分別為(22，800)、(68，800)。從圖中可以看出傳統(tǒng)方法追蹤第一層層位時，在剖面左側(cè)和右側(cè)均出現(xiàn)了層位追蹤偏離現(xiàn)象，而在追蹤第二層層位時存在層位追蹤串層現(xiàn)象。而本文方法不需要設(shè)置種子點，降低了人為干預(yù)成本，還能夠自動提取異常體強(qiáng)振幅同相軸反射的層位數(shù)據(jù)。現(xiàn)場案例分析結(jié)果驗證了本文方法的有效性和適應(yīng)性。

圖8 兩種方法層位追蹤結(jié)果對比

4 結(jié)論

針對傳統(tǒng)層位追蹤方法普遍存在主觀性強(qiáng)、耗時耗力和每次僅能追蹤1個層位等問題。本文提出了一種探地雷達(dá)公路路面多層位自動追蹤方法。通過復(fù)信號分析獲取了雷達(dá)數(shù)據(jù)的瞬時相位余弦，并通過構(gòu)建子波余弦矩陣對瞬時相位余弦進(jìn)行橫向增強(qiáng)處理，提高了相位數(shù)據(jù)同相軸反射的橫向連續(xù)性；基于反射信號波形類似Ricker子波的特點，利用雷達(dá)數(shù)據(jù)的波峰、波谷、幅值和極性等信息來確定層位數(shù)據(jù)；最后根據(jù)層位線閾值、振幅閾值過濾一些無關(guān)或不重要的層位線，實現(xiàn)有效強(qiáng)反射同相軸層位自動追蹤。數(shù)值模擬和現(xiàn)場案例分析結(jié)果驗證了本文方法的有效性和適應(yīng)性。本文方法不僅可以自動追蹤多個層位數(shù)據(jù)，還可以估計出層位的極性，研究成果可為探地雷達(dá)地層分析、層反射介質(zhì)反演、沉積體系域解釋提供技術(shù)支持。