曹蕓茜 吳仁彪 何煒琨

(1. 中國(guó)民航大學(xué)智能信號(hào)與圖像處理天津市重點(diǎn)實(shí)驗(yàn)室，天津 300300; 2. 中國(guó)民航大學(xué)基礎(chǔ)實(shí)驗(yàn)中心，天津 300300)

1 引言

機(jī)場(chǎng)作為城市和航空運(yùn)輸?shù)幕A(chǔ)設(shè)施，是交通運(yùn)輸?shù)闹匾M成部分。機(jī)場(chǎng)道面是機(jī)場(chǎng)建設(shè)的核心，跑道質(zhì)量的好壞將會(huì)直接影響機(jī)場(chǎng)與飛機(jī)運(yùn)行的安全，其中場(chǎng)道內(nèi)脫空和裂縫等隱性災(zāi)害的檢測(cè)和維護(hù)至關(guān)重要[1]。在不影響機(jī)場(chǎng)正常運(yùn)行且不破壞道面結(jié)構(gòu)的基礎(chǔ)上，對(duì)道面進(jìn)行無(wú)損檢測(cè)是必然趨勢(shì)。探地雷達(dá)(Ground Penetrating Radar, GPR)是一種有效的無(wú)損檢測(cè)工具，其具有探測(cè)速度快、操作方便、分辨率高等特點(diǎn)[2-3]。在機(jī)場(chǎng)場(chǎng)道的建設(shè)施工中，為了加強(qiáng)其承載能力，需要鋪設(shè)鋼筋進(jìn)行加固[4]。鋼筋的強(qiáng)反射回波甚至多次反射后的回波都會(huì)淹沒(méi)災(zāi)害目標(biāo)的回波，所以有效的抑制鋼筋強(qiáng)反射回波對(duì)災(zāi)害目標(biāo)的檢測(cè)是十分關(guān)鍵和必要的。

本文首先利用非一致性檢測(cè)對(duì)探地雷達(dá)接收數(shù)據(jù)進(jìn)行預(yù)處理，利用時(shí)變?cè)鲆婧瘮?shù)對(duì)鋼筋反射一次回波進(jìn)行粗抑制，然后對(duì)該數(shù)據(jù)在時(shí)間維進(jìn)行分段均衡，以濾除鋼筋在機(jī)場(chǎng)場(chǎng)道下各層分界面的多次反射回波，再利用時(shí)頻—時(shí)空域二重S變換方法對(duì)災(zāi)害目標(biāo)進(jìn)行檢測(cè)，并同時(shí)顯示其在場(chǎng)道下的空間位置。該方法能在災(zāi)害目標(biāo)回波信號(hào)被鋼筋多次強(qiáng)反射回波信號(hào)淹沒(méi)的情況下，有效抑制鋼筋對(duì)災(zāi)害目標(biāo)的影響，達(dá)到檢測(cè)道面下微小災(zāi)害目標(biāo)空間位置的目的，為機(jī)場(chǎng)養(yǎng)護(hù)施工提供直觀依據(jù)。

2 信號(hào)模型

探地雷達(dá)接收信號(hào)模型可以表示為：

w(t)=o(t)+d(t)+s(t)

(1)

如圖1所示，o(t)表示收發(fā)天線直接耦合波，在收發(fā)天線間距確定后直接耦合波分量一般相對(duì)固定。d(t)表示各層媒質(zhì)表面反射回波，由于探測(cè)環(huán)境往往是變化的，所以d(t)的變化性和不可預(yù)測(cè)性較大。s(t)表示媒質(zhì)內(nèi)目標(biāo)反射回波，也就是我們需要的回波分量。

圖1 GPR信號(hào)模型

我們將目標(biāo)反射回波s(t)進(jìn)行采樣，那么二維B-Scan數(shù)據(jù)中，雷達(dá)測(cè)試線l=1,2,…,L的第i=1,2,…,T個(gè)采樣點(diǎn)的目標(biāo)反射回波信號(hào)可以表示為sl(i)。

3 鋼筋回波抑制

本節(jié)研究鋼筋的強(qiáng)反射回波及多次反射后回波的抑制方法。首先利用基于廣義內(nèi)積的非一致性檢測(cè)對(duì)探地雷達(dá)接收數(shù)據(jù)進(jìn)行預(yù)處理，然后利用時(shí)變?cè)鲆婧瘮?shù)對(duì)鋼筋的一次強(qiáng)反射回波進(jìn)行粗抑制，最后利用時(shí)間域的分段均衡算法抑制鋼筋的多次反射回波。

3.1 預(yù)處理

我們將探地雷達(dá)接收到的回波信號(hào)w(t)進(jìn)行預(yù)處理，去除直接耦合波o(t)和各層媒質(zhì)表面反射回波d(t)，以提取出鋼筋及災(zāi)害目標(biāo)的反射回波s(t)。

首先去除探地雷達(dá)接收到的回波信號(hào)中收發(fā)天線間的直接耦合波o(t)。我們可以將雷達(dá)天線對(duì)準(zhǔn)開(kāi)闊無(wú)物體的空間，或?qū)⑻炀€置于暗室內(nèi)，錄取對(duì)空信號(hào)，此信號(hào)可以近似為直接耦合信號(hào)，再將雷達(dá)接收回波減去此直接耦合信號(hào)即有效去除了直接耦合波。

然后去除機(jī)場(chǎng)場(chǎng)道下各層媒質(zhì)表面反射回波d(t)。由于探測(cè)環(huán)境的變化，d(t)的變化性和不可預(yù)知性較大，所以，我們需要根據(jù)探測(cè)區(qū)域的數(shù)據(jù)對(duì)d(t)進(jìn)行有效的去除。這里，我們采用基于廣義內(nèi)積的非一致性檢測(cè)方法合理去除d(t)。

基于廣義內(nèi)積的非一致性檢測(cè)算法，其輸出定義為：

(2)

其中l(wèi)=1,2,…,L為雷達(dá)測(cè)試線數(shù)，xl為去除直接耦合波后的雷達(dá)采樣數(shù)據(jù)的矩陣表示，Rl為自相關(guān)矩陣，其可以表示為：

(3)

計(jì)算出雷達(dá)所有探測(cè)點(diǎn)的gl(GIP)后，給出閾值，大于此閾值則說(shuō)明與周?chē)町愝^大，則可認(rèn)為此部分探測(cè)區(qū)域存在目標(biāo)反射信號(hào)，我們將此類(lèi)區(qū)域數(shù)據(jù)濾除，剩下的小于閾值部分的回波即為各層媒質(zhì)表面強(qiáng)反射回波，選取該回波信號(hào)作為背景回波，將雷達(dá)接收回波減去背景回波即可以有效去除各層媒質(zhì)表面反射回波而有效的保留了目標(biāo)反射回波。

3.2 鋼筋回波抑制

在機(jī)場(chǎng)跑道使用過(guò)程中，初期常見(jiàn)的隱性災(zāi)害主要有脫空和裂縫[5- 6]。脫空主要出現(xiàn)在面層和上基層的交界處，裂縫主要出現(xiàn)在面層，而加固的鋼筋主要位于面層的上部，如圖2所示，所以鋼筋的強(qiáng)反射回波會(huì)嚴(yán)重影響甚至淹沒(méi)脫空和裂縫的反射回波，并且鋼筋與場(chǎng)道各層介面的多次反射回波也會(huì)嚴(yán)重影響場(chǎng)道災(zāi)害目標(biāo)的檢測(cè)，所以有效地抑制鋼筋的強(qiáng)反射回波對(duì)災(zāi)害目標(biāo)進(jìn)行檢測(cè)是十分關(guān)鍵和必要的。

圖2 鋼筋及災(zāi)害目標(biāo)分布模型

由于鋼筋的一次反射回波遠(yuǎn)遠(yuǎn)強(qiáng)于多次反射回波，且多次反射回波往往會(huì)與災(zāi)害目標(biāo)反射回波混疊，不易處理，所以我們將鋼筋回波抑制分為兩步，第一步利用時(shí)變?cè)鲆婧瘮?shù)對(duì)鋼筋一次反射回波進(jìn)行粗抑制，第二步利用時(shí)域分段均衡方法對(duì)鋼筋的多次反射回波進(jìn)行細(xì)抑制，以便提取出有效的災(zāi)害目標(biāo)反射回波。

3.2.1粗抑制—時(shí)變?cè)鲆?/p>

機(jī)場(chǎng)場(chǎng)道在進(jìn)行鋼筋加固時(shí)有嚴(yán)格的規(guī)范并且有先驗(yàn)的鋼筋位置信息，我們可以根據(jù)鋼筋位置的先驗(yàn)知識(shí)，構(gòu)造合適的時(shí)變?cè)鲆婧瘮?shù)Pl(t)，其中l(wèi)=1,2,…,L為雷達(dá)測(cè)試線，將其采樣后記為Pl(i)，那么時(shí)變?cè)鲆婧蟮睦走_(dá)數(shù)據(jù)為

ul(i)=Pl(i)·sl(i),i=1,2,…,T,l=1,2,…,L

(4)

時(shí)變?cè)鲆婧瘮?shù)可以抑制鋼筋一次強(qiáng)反射回波時(shí)間窗內(nèi)數(shù)據(jù)，增強(qiáng)需要觀察的機(jī)場(chǎng)場(chǎng)道面層及基層時(shí)間窗內(nèi)數(shù)據(jù)，抑制不需要觀察的場(chǎng)道底部數(shù)據(jù)，即選擇需要觀察的數(shù)據(jù)予以增強(qiáng)，不需要觀察的數(shù)據(jù)予以抑制，對(duì)鋼筋的一次反射回波進(jìn)行粗抑制，此方法可以有效抑制鋼筋一次反射回波，而增強(qiáng)了災(zāi)害目標(biāo)反射回波，但是鋼筋的多次反射回波也同災(zāi)害目標(biāo)反射回波一同被增強(qiáng)，所以我們還要對(duì)鋼筋的多次反射回波進(jìn)一步抑制。

3.2.2細(xì)抑制—分段均衡

如果某一深度存在目標(biāo)，則此深度目標(biāo)所在位置必與同一深度其他位置的回波有差異且鋼筋對(duì)各層表面多次反射回波呈現(xiàn)較強(qiáng)的規(guī)律性，所以本文提出分段均衡的方法，對(duì)B-Scan數(shù)據(jù)進(jìn)行時(shí)間維分段，再利用數(shù)據(jù)信息求出分段矢量，對(duì)每一段數(shù)據(jù)進(jìn)行均衡，這樣可以對(duì)相對(duì)規(guī)律的數(shù)據(jù)進(jìn)行均衡濾波，對(duì)存在局部差異的數(shù)據(jù)起到增強(qiáng)的效果，有利于細(xì)節(jié)的辨認(rèn)。本文利用分段均衡的方法來(lái)抑制鋼筋的多次反射回波。

首先計(jì)算鋼筋粗抑制后的數(shù)據(jù)中所有雷達(dá)測(cè)試線l=1,2,…,L在第i個(gè)采樣點(diǎn)的平均值：

(5)

根據(jù)災(zāi)害目標(biāo)的厚度合理選取分段數(shù)據(jù)點(diǎn)數(shù)m=T/K，其中K為分段數(shù)目，T為總時(shí)間采樣點(diǎn)數(shù)。局部搜索第k=1,2,…,K段的最大值maxk及最小值mink，定義差值矢量為

d=[max1-min1max2-min2…maxK-minK]

(6)

為了方便將上式記為

d=[d1d2…dk]

(7)

設(shè)dmax為d中最大值，定義均衡放大矢量為

a=[d1/dmaxd2/dmax…dK/dmax]

(8)

記為a=[a1a2…ak]。

鋼筋粗抑制后的B-Scan接收數(shù)據(jù)用矩陣可以表示為u。根據(jù)上述分段方法，按照相同原理對(duì)u在時(shí)間維進(jìn)行同樣的分段。設(shè)uk為對(duì)u進(jìn)行分段的第k段數(shù)據(jù)，則uk可以表示為：

uk={ul(i),i=(k-1)m+1,(k-1)m+2,…,km;

l=1,2,…,L;k=1,2,…,K}

(9)

利用每一段數(shù)據(jù)的均衡放大系數(shù)對(duì)相應(yīng)的該段數(shù)據(jù)進(jìn)行加權(quán)均衡，得到第k段數(shù)據(jù)的均衡處理結(jié)果：

yk=akuk

(10)

其中，k=1,2,…,K。那么，B-Scan接收數(shù)據(jù)的均衡處理結(jié)果可以表示為：

Y=[y1y2…yK]T

(11)

4 災(zāi)害目標(biāo)檢測(cè)

鋼筋回波抑制后的數(shù)據(jù)在災(zāi)害目標(biāo)附近可能還會(huì)存在少許干擾回波，這也會(huì)干擾體積或厚度較小的災(zāi)害目標(biāo)，譬如小裂縫、脫空薄層的檢測(cè)。本文提出了時(shí)頻—時(shí)空域二重S變換災(zāi)害目標(biāo)檢測(cè)方法，可以有效的檢測(cè)體積及厚度較小的災(zāi)害目標(biāo)。

S變換[7-9]是基于短時(shí)窗傅里葉變換和小波變換提出來(lái)的一種時(shí)頻分析方法，該方法吸收了兩者的優(yōu)點(diǎn)，同時(shí)克服了它們的缺點(diǎn)。

設(shè)函數(shù)h(t)∈L2(R)(L2(R)表示能量有限的函數(shù)空間)的S變換的定義如下：

(12)

S變換的基本小波定義為：

(13)

(14)

由式(14)可以看出，S變換采用寬度隨頻率呈反比變化的高斯窗函數(shù)，低頻段的時(shí)窗較寬，可以獲得較高的頻率分辨率，高頻段的時(shí)窗較窄，可獲得較高的時(shí)間分辨率[10]。利用高頻段時(shí)間分辨率較高的特性可以有效的檢測(cè)出體積及厚度較小的災(zāi)害目標(biāo)。

首先選取均衡后數(shù)據(jù)Y中的A-Scan數(shù)據(jù)yl(t)，雷達(dá)測(cè)試線l∈[1,L]。對(duì)yl(t)進(jìn)行時(shí)間—頻率維S變換：

(15)

根據(jù)S變換結(jié)果選取時(shí)間分辨率較高的頻率段f0…fN，固定該頻率段中的一個(gè)頻率點(diǎn)fn(n∈[0,N])，在該頻率點(diǎn)下對(duì)yl(t)再次進(jìn)行S變換：

(16)

對(duì)B-Scan數(shù)據(jù)Y中的所有測(cè)試線l=1,2,…,L的數(shù)據(jù)yl(t)分別進(jìn)行S變換，即可得到一個(gè)固定頻點(diǎn)fn的時(shí)間—距離維S變換結(jié)果：

Sfn(τ,l)={Sl(τ,fn),l=1,…,L}

(17)

將選取的頻率段f0…fN中所有頻點(diǎn)的時(shí)間—距離維S變換數(shù)據(jù)疊加，進(jìn)行頻域均衡濾波：

(18)

根據(jù)疊加結(jié)果不但可以有效檢測(cè)出體積及厚度較小的災(zāi)害目標(biāo)，還可以得到災(zāi)害目標(biāo)的空間位置。

綜上所述，機(jī)場(chǎng)場(chǎng)道鋼筋加固區(qū)的災(zāi)害目標(biāo)檢測(cè)方法的流程如圖3所示。

圖3 災(zāi)害目標(biāo)檢測(cè)方法流程

5 仿真結(jié)果與分析

為了驗(yàn)證本文方法的有效性，實(shí)驗(yàn)中所有數(shù)據(jù)均由GPRMAX2.0[11]高保真仿真軟件產(chǎn)生，采用雙基地探地雷達(dá)，選取天線間距為0.1 m，距地面高度為0.04 m，發(fā)射信號(hào)采用ricker波脈沖形式，中心頻率為900 MHz，采集數(shù)據(jù)時(shí)間窗為25 ns，采集測(cè)試線間距為0.02 m，共采集118道數(shù)據(jù)。仿真模型為一般機(jī)場(chǎng)場(chǎng)道基本四層模型，層厚從上到下依次為0.34 m,0.20 m,0.20 m,0.26 m，介電常數(shù)依次為9,12,15,22。

圖4為本實(shí)驗(yàn)中所采用的機(jī)場(chǎng)道面仿真數(shù)據(jù)原始幾何模型圖，該模型中含有五根鋼筋，一條裂縫和一個(gè)脫空，鋼筋均勻分布在面層中距地表面17 cm，每條鋼筋間距為50 cm，裂縫長(zhǎng)18 cm，寬5 mm，脫空寬20 cm，厚7 mm。

圖4 原始數(shù)據(jù)幾何模型

圖5為探地雷達(dá)接收到的二維B-Scan回波數(shù)據(jù)。由于地表反射回波很強(qiáng)，淹沒(méi)了地下目標(biāo)的回波信息，所以首先必須有效的去除地表反射回波，提取需要的回波信息。圖6為利用非一致性檢測(cè)濾除背景回波后的目標(biāo)反射回波。背景的四層媒質(zhì)表面反射回波被有效濾除，凸顯出地下目標(biāo)的反射回波，但是由于場(chǎng)道淺層鋼筋反射回波很強(qiáng)，且存在鋼筋多次反射回波，嚴(yán)重影響甚至淹沒(méi)鋼筋下災(zāi)害目標(biāo)的回波，無(wú)法有效地檢測(cè)出災(zāi)害目標(biāo)，尤其是細(xì)小的裂縫。

圖6 目標(biāo)反射回波數(shù)據(jù)

根據(jù)圖6所示數(shù)據(jù)，選擇合適的時(shí)變?cè)鲆鏁r(shí)窗，對(duì)強(qiáng)金屬回波及不需要觀察的場(chǎng)道深層處數(shù)據(jù)進(jìn)行抑制，對(duì)需要觀察的數(shù)據(jù)段進(jìn)行增強(qiáng)，如圖7所示。圖8為鋼筋回波粗抑制后B-Scan數(shù)據(jù)，脫空的反射回波被增強(qiáng)，可以觀察到，但是由于裂縫上下表面比較窄，反射回波很弱，并且鋼筋與地下各層媒質(zhì)表面的多次反射回波也比較強(qiáng)，所以裂縫還是不能有效觀測(cè)。圖9為分段均衡后B-Scan數(shù)據(jù)，鋼筋在各層媒質(zhì)間的多次反射回波被有效均衡，突出了脫空的反射回波，但仍無(wú)法有效檢測(cè)出災(zāi)害目標(biāo)。

圖7 時(shí)變?cè)鲆?/p>

圖8 鋼筋回波粗抑制

圖9 鋼筋回波細(xì)抑制

圖10為均衡濾波后采集數(shù)據(jù)進(jìn)行時(shí)間—頻率維S變換的時(shí)頻分析結(jié)果，我們選取其中時(shí)間分辨率較高的高頻段2500～2900 MHz數(shù)據(jù)，先確定該頻段內(nèi)的一個(gè)頻點(diǎn)，對(duì)均衡后B-Scan數(shù)據(jù)進(jìn)行時(shí)間—距離維S變換，將所有頻段內(nèi)的S變換結(jié)果進(jìn)行頻域均衡濾波，得到如圖11所示結(jié)果，可以明顯檢測(cè)出脫空及裂縫的存在。圖12為不進(jìn)行均衡濾波而直接用相同的方法進(jìn)行S變換的結(jié)果，雖然能夠檢測(cè)到脫空，但是不能有效地檢測(cè)到裂縫的存在，且鋼筋的多次反射回波嚴(yán)重影響了災(zāi)害目標(biāo)的檢測(cè)。

圖10 S變換時(shí)頻分析

圖11 災(zāi)害目標(biāo)檢測(cè)結(jié)果

圖12 未分段均衡檢測(cè)結(jié)果

6 結(jié)論

考慮到鋼筋的強(qiáng)反射回波對(duì)場(chǎng)道鋼筋加固區(qū)災(zāi)害目標(biāo)檢測(cè)的影響，本文首先提出時(shí)變?cè)鲆婧瘮?shù)和分段均衡的方法對(duì)鋼筋的多次反射回波進(jìn)行抑制，然后提出時(shí)頻—時(shí)空域二重S變換方法對(duì)災(zāi)害目標(biāo)進(jìn)行檢測(cè)。仿真結(jié)果表明本文算法可以有效去除鋼筋的多次反射回波對(duì)厚度及體積較小的常見(jiàn)災(zāi)害目標(biāo)回波的影響，對(duì)災(zāi)害目標(biāo)進(jìn)行檢測(cè)并顯示其在場(chǎng)道下的空間位置，為機(jī)場(chǎng)場(chǎng)道的養(yǎng)護(hù)施工提供直觀的依據(jù)。

[1] 翁興中,蔡良才.機(jī)場(chǎng)道面設(shè)計(jì)[M].北京:人民交通出版社,2007.

Weng Xingzhong, Cai Liangcai. Design of Airport Runway[M].Beijing: China Communications Press, 2007.(in Chinese)

[2] 何煒琨, 吳仁彪, 劉家學(xué). 基于探地雷達(dá)的機(jī)場(chǎng)場(chǎng)道脫空層檢測(cè)及厚度估計(jì)[J]. 信號(hào)處理, 2011, 27(10):1547-1551.

He Weikun, Wu Renbiao, Liu Jiaxue. Void-layer Detection and Depth Determination in Runways based on GPR[J]. Signal Processing, 2011, 27(10):1547-1551.(in Chinese)

[3] Hoarau Q, Ginolhac G, Atto A M, et al. Robust Adaptive Detection of Buried Pipes using GPR[C]∥24th European Signal Processing Conference, 2016:533-537.

[4] Harry M J. Ground penetrating radar theory and applications[M]. London: Elsevier, 2009.

[5] 何煒琨. 基于探地雷達(dá)的機(jī)場(chǎng)場(chǎng)道質(zhì)量監(jiān)測(cè)關(guān)鍵技術(shù)研究[D]. 天津：天津大學(xué), 2012.

He Weikun. Enabling Techniques for Runways Quality Surveillance via Ground Penetrating Radar[D]. Tianjin: The University of Tianjin, 2012.(in Chinese)

[6] 龐希斌, 徐進(jìn), 盧小賓, 等. 地質(zhì)雷達(dá)在機(jī)場(chǎng)跑道缺陷檢測(cè)中的應(yīng)用——以“5·12”汶川地震后九寨黃龍機(jī)場(chǎng)檢測(cè)為例[J]. 西南民族大學(xué)學(xué)報(bào), 2008, 34(6):1096-1100.

Pang Xibin, Xu Jin, Lu Xiaobin, et al. Application of the ground penetrating radar to the inspection of the airfield runway defect——Taking the detection of Jiuhuang Airport after the 5.12 Wenchuan earthquake for example[J]. Journal of Southwest University for Nationalities Natural Science Edition, 2008, 34(6):1096-1100.(in Chinese)

[7] Yan Jin, Liu Jie. Parameter estimation of frequency hopping signals based on the Robust S-transform algorithms in alpha stable noise environment[J]. International Journal of Electronics & Communications, 2016, 70(5): 611- 616.

[8] Kazemi K, Amirian M, Dehghani M J. The S-transform using a new window to improve frequency and time resolutions[J]. Signal, Image and Video Processing, 2014, 8(3): 533-541.

[9] 劉志成, 王殿偉. 非線性調(diào)頻信號(hào)的自適應(yīng)時(shí)頻濾波算法[J]. 信號(hào)處理, 2015, 31(3):356-363.

Liu Zhicheng, Wang Dianwei.Adaptive Time-Frequency Filter Method of Nonlinear Frequency Modulation Signal[J]. Journal of Signal Processing, 2015, 31(3):356-363.(in Chinese)

[10]Huang Nantian, Peng Hua, Cai Guowei, et al. Power Quality Disturbances Feature Selection and Recognition Using Optimal Multi-Resolution Fast S-Transform and CART Algorithm[J]. Energies, 2016, 9(11): 1-21.

[11]GprMax2.0[OL]. http:∥www. gprmax.org.