基于擬合誤差消除的探地雷達(dá)圖像魯棒雙曲線識別模型

蘭 天 趙 毅 陳宏暢 龔俊波 王長軍 王 健 楊小鵬

（1.北京理工大學(xué)信息與電子學(xué)院雷達(dá)技術(shù)研究所，北京 100081；2.北京理工大學(xué)重慶創(chuàng)新中心，重慶 401120；3.北京中建建筑科學(xué)研究院有限公司，北京 100070）

1 引言

作為探測地下未知物體的無損工具，探地雷達(dá)（ground penetrating radar，GPR）已廣泛應(yīng)用于各個領(lǐng)域［1］，如軍事［2］、土木工程［3］、橋梁測量［4］、分層探測［5-6］等。脈沖探地雷達(dá)是一種常見的探地雷達(dá)系統(tǒng)，它向地下發(fā)射短電磁脈沖，并接收物體反射的信號，以探測地下區(qū)域。一般來說，探地雷達(dá)B-scan圖像中生成的雙曲線形特征［7］代表地下目標(biāo)。然而，探地雷達(dá)B-scan 圖像中的地下目標(biāo)通常被雜波所掩蓋，雜波主要由發(fā)射天線和接收天線之間的直接耦合、來自地面的反射以及來自地下不連續(xù)性的散射響應(yīng)引起［8］。因此，從探地雷達(dá)B-scan 圖像中自動提取雙曲線具有挑戰(zhàn)性。

由于霍夫變換技術(shù)［9］成功地用于變形形狀擬合，因此它被用于在GPR B-scan 圖像中提取雙曲線，以定位埋藏物體［10］。此外，在［11］中，根據(jù)未知參數(shù)與實(shí)驗(yàn)誤差的差異引入了加權(quán)因子，擴(kuò)展了霍夫變換。此外，霍夫變換技術(shù)還用于揭示物體位置、深度、半徑或速度［12］。由于基于霍夫變換的方法需要預(yù)先指定模型，計(jì)算量大，空間開銷大，因此這些方法在實(shí)際應(yīng)用中存在局限性。與霍夫變換不同，最小二乘法可以搜索和區(qū)分GPR B-scan 圖像中的二次曲線［13-14］。然而，如果沒有圖像分割步驟，基于霍夫變換的大多數(shù)方法只能識別圖像中的一條曲線。在文獻(xiàn)［15-20］中，基于機(jī)器學(xué)習(xí)的方法被用來提取B 掃描圖像中的雙曲線等目標(biāo)區(qū)域。其中，在文獻(xiàn)［15］中，提出了一個由反向傳播神經(jīng)網(wǎng)絡(luò)分類器、基于霍夫變換的模式識別方法、基于韋爾奇功率譜密度估計(jì)和邊緣檢測的圖像特征處理方法組成的處理系統(tǒng)模型。在文獻(xiàn)［16］中，Viola-Jones（VJ）算法被用來提取探地雷達(dá)數(shù)據(jù)中的目標(biāo)區(qū)域。文獻(xiàn)［17］中的方法在檢測時，采用一種基于神經(jīng)網(wǎng)絡(luò)的方法對噪聲信號進(jìn)行分類。在文獻(xiàn)［18］中，作者提出了包括數(shù)據(jù)過濾、檢測和目標(biāo)分類的雙曲線自動檢測和分類方法，在數(shù)據(jù)過濾部分利用了小波分析和Gabar 濾波器，在檢測分類部分提出了基于層次分析的目標(biāo)分類方法。在文獻(xiàn)［19］中，作者提出了一種基于相位對稱的手工特征，稱為定向矢量相位對稱直方圖（histogram of oriented vector phase symmetry，HOVPS），該特征作為分類器的特征輸入以實(shí)現(xiàn)對雙曲線的檢測。對于上述基于機(jī)器學(xué)習(xí)的方法，大多數(shù)特征需要專家來識別，并且訓(xùn)練數(shù)據(jù)并不足以滿足所有不同的檢測條件。此外，在文獻(xiàn)［20］中提出了一種結(jié)合深度學(xué)習(xí)網(wǎng)絡(luò)RetinaNet 的自動檢測雙曲線的方法，取得了0.58的平均精度。近年來，有一些研究實(shí)現(xiàn)了能夠自動檢測和擬合B-scan 圖像中雙曲線的集成模型［21-22］。在文獻(xiàn)［21］所提出的模型中，首先對Bscan 圖像進(jìn)行預(yù)處理，將其轉(zhuǎn)化為二值圖像。然后提出了列連接聚類算法（column-connection clustering，C3）來分離有交集的雙曲線。之后，采用基于神經(jīng)網(wǎng)絡(luò)的方法從C3 算法的輸出中識別雙曲線特征。雖然該模型包含了圖像處理、機(jī)器學(xué)習(xí)等方法，但該模型的C3算法沒有考慮雙曲線垂直向下開口這一重要特征。為了解決這一問題，文獻(xiàn)［22］提出了開放掃描聚類算法（open-scan-clustering algorithm，OSCA）。OSCA 通過從上到下掃描預(yù)處理后的二值圖像中向下的開口來進(jìn)行聚類。然后，基于拋物線擬合的判斷（fitting-based judgment，PFJ）方法從OSCA 的輸出中識別雙曲線特征。然而OSCA中的下開口特征是B-scan 圖像的像素級局部化特征。因此，OSCA 的輸出包含一些具有下開口特征的非雙曲點(diǎn)集。而且，PFJ 方法不能有效地消除這些非雙曲點(diǎn)集。為了去除這些非雙曲點(diǎn)簇，本文提出了一種魯棒集成模型。該模型由預(yù)處理方法、OSCA、基于代數(shù)距離的雙曲線擬合算法和基于擬合誤差的消除（fitting-errors-based eliminating，F(xiàn)EE）方法組成。其中，F(xiàn)EE 方法根據(jù)非雙曲點(diǎn)簇的擬合誤差遠(yuǎn)大于雙曲點(diǎn)簇的擬合誤差的特點(diǎn)消除了OSCA 生成的開口向下的非雙曲點(diǎn)簇，實(shí)現(xiàn)了GPR圖像中所有雙曲點(diǎn)簇的擬合和識別。

本文其余部分組織如下。預(yù)處理方法在第2節(jié)中介紹，然后在第3 節(jié)介紹OSCA，隨后在第4 節(jié)介紹雙曲線擬合算法，第5 節(jié)介紹FEE 方法，第6 節(jié)展示實(shí)驗(yàn)研究，最后在第7節(jié)得出結(jié)論。

2 預(yù)處理方法

在本節(jié)中，介紹了預(yù)處理過程，包括均值對消、自適應(yīng)閾值算法和開閉運(yùn)算。

2.1 均值對消

在應(yīng)用所提出的閾值分割算法之前，采用每個掃描道數(shù)據(jù)減去所有掃描道數(shù)據(jù)的平均值的操作來抑制雜波和噪聲。均值相減運(yùn)算前后的圖像如圖1所示。圖1（a）中的亮條被抑制，這表示發(fā)射和接收天線與地面反射的直接耦合。為了評估平均減法操作的有效性，利用信號雜波比（SCR）并定義為

圖1 均值相減操作Fig.1 Effect of mean subtraction operation

其中，BM×N表示B掃描圖像，TI×J表示B掃描圖像中的目標(biāo)區(qū)域。經(jīng)過均值相減運(yùn)算，SCR 得到改善，如表1所示。

表1 均值相減前后的SCRTab.1 SCR of images before and after mean subtraction operation

2.2 自適應(yīng)閾值算法

由于探地雷達(dá)圖像中垂直灰度值的變化代表了介質(zhì)不連續(xù)面的反射，因此對垂直灰度值變化點(diǎn)的灰度值進(jìn)行處理，計(jì)算閾值。梯度G(x，y)可由

其中V(x，y)是灰度值。由于均值對消運(yùn)算后存在殘余雜波和噪聲，若將閾值計(jì)算為

閾值過低，實(shí)驗(yàn)結(jié)果就會變差。因此，本文算法中的閾值定義為

在（4）中，Ve(x，y)代表大于最大灰度值的一定百分比的像素點(diǎn)灰度值，Vy(x，y)代表在G(x，y) ≠0處的像素點(diǎn)灰度值集合，MaxVy(x，y)代表集合中的最大灰度值，ρ是一個分?jǐn)?shù)（0 ＜ρ＜1）。為了分析ρ的影響，對所提出的閾值分割算法進(jìn)行了不同ρ值的測試。如圖2所示，ρ=0.5的閾值太低，只有一些非目標(biāo)區(qū)域被去除；ρ=0.7的閾值太高，雙曲特征變得不完整。因此，本文提出的閾值算法采用ρ=0.6。該閾值算法基于［22］中的自適應(yīng)閾值法。為了比較性能，這兩種方法在真實(shí)的B掃描圖像上進(jìn)行了評估，如圖3所示。［22］中提出的方法得到的閾值過低，只能去除一些暗區(qū)。如表2所示，圖3（c）的SCR高于圖3（b），因此本文提出的閾值算法性能優(yōu)于［22］中的閾值算法。

表2 不同閾值處理后的圖像SCR值Tab.2 SCR of images processed by different thresholding methods

圖2 所提出的方法不同ρ值下閾值處理后的圖像Fig.2 Thresholded image by the proposed method with different values of ρ

圖3 不同閾值方法處理后的圖像Fig.3 Thresholded image of different methods

2.3 開閉運(yùn)算

打開和關(guān)閉操作是兩個形態(tài)操作［23］。這兩種操作都是兩種基本形態(tài)操作的組合：擴(kuò)張操作和侵蝕操作。膨脹操作用于連接相似特征的區(qū)域，侵蝕操作用于去除明亮的噪聲點(diǎn)。二值圖像I通過結(jié)構(gòu)元素SE的打開操作定義為：

這表明I受到SE1 的侵蝕，接著是SE2 的膨脹。同理，通過SE對I的關(guān)閉操作定義為：

這表明了SE1 對I的膨脹以及SE2 對結(jié)果的侵蝕。

在所提出的處理方法中，首先應(yīng)用閉操作來融合窄裂縫和細(xì)長裂縫，消除小孔，填充輪廓中的空隙。然后使用開操作來平滑物體的輪廓，打破狹窄的峽部，消除小的突出。為了最大限度地保留圖像，在閉操作中使用半徑為3的圓形結(jié)構(gòu)元素，在開操作中使用半徑為2的元素。開閉操作后的處理結(jié)果如圖4（c）所示。

圖4 所提出的預(yù)處理方法后的圖像Fig.4 Images after the proposed preprocessing method

3 開口掃描聚類算法

經(jīng)過OSCA 處理后的B 掃描圖像，可以去除大部分非雙曲線點(diǎn)簇。OSCA 基于一些概念，包括“點(diǎn)段”、“下/上開口”和“點(diǎn)段參數(shù)存儲陣列”?！包c(diǎn)段”指的是連續(xù)的兩個以上的點(diǎn)。一個點(diǎn)段與另兩個不連接的點(diǎn)段水平重疊，它們之間在下一行中有間隙，上面的點(diǎn)段定義為開始向下的開口。同樣，向上開口也可以定義為位于連續(xù)兩行的三個連通點(diǎn)段，其中一個點(diǎn)段與前一行的另外兩個不連通點(diǎn)段水平重疊，且它們之間有間隙，較低的點(diǎn)段定義為開始一個向上開口?！包c(diǎn)段”和“下/上”開口如圖5 所示?！包c(diǎn)段參數(shù)存儲數(shù)組”是存儲點(diǎn)段參數(shù)的數(shù)組，在（7）中定義，（7）中的符號在表3中解釋。

表3 公式（7）中符號的含義Tab.3 Meaning of symbols in Eq.（7）

圖5 開口像素表示Fig.5 Pixel representation of opening

通過生成該陣列，OSCA 的操作思路可以從二值圖像轉(zhuǎn)移到“點(diǎn)段參數(shù)存儲陣列”。根據(jù)上述定義，通過掃描每一行中的點(diǎn)段和相鄰行中的點(diǎn)段，可以找到向下/向上的開口。在第一次掃描中，該算法的目標(biāo)是找到向上的開口。找到向上開口后，標(biāo)記向上開口的點(diǎn)段，并連續(xù)標(biāo)記其下方的相鄰點(diǎn)段。當(dāng)標(biāo)記出向上開口的底部點(diǎn)段時，算法會檢查該點(diǎn)段是否開始向下開口。如果出現(xiàn)這種情況，則會出現(xiàn)x 形相交的情況，在第二次掃描時對這種情況進(jìn)行操作。否則，這個向上開放簇的所有點(diǎn)段將被刪除。在第二次掃描中，當(dāng)發(fā)現(xiàn)向下的開口［見圖6（a）］時，算法將檢查該點(diǎn)段是否開始向上的開口。如果沒有，則在點(diǎn)段參數(shù)存儲陣列中記錄開放范圍[S，E]，見圖6（b），并構(gòu)建新的點(diǎn)集。向下開口點(diǎn)簇的頂部區(qū)域可以向上尋找小于或等于下面開口點(diǎn)簇的點(diǎn)段，直到?jīng)]有點(diǎn)段滿足這一要求［見圖6（c）］。通過掃描左右兩側(cè)向下開口下方重疊的點(diǎn)段可以得到尾部［見圖6（c）］。頂部和尾部的所有點(diǎn)將被添加到構(gòu)造的開放點(diǎn)集群中。如果發(fā)生這種情況，X 形交叉上的處理見文獻(xiàn)［22］。OSCA 算法的偽代碼見算法1。通過使用OSCA，大多數(shù)沒有下開口的非雙曲點(diǎn)集在B掃描圖像中被消除。

圖6 OSCA中概念的像素表示Fig.6 Pixel representation of concepts in OSCA

公式（7）中各參數(shù)的含義如表3所示：

4 基于代數(shù)距離的雙曲線擬合算法

在雙曲線擬合算法之前，介紹了一種擬合點(diǎn)提取方法。通過在每個向下開放點(diǎn)簇的每列中自動提取一組均位點(diǎn)，得到擬合點(diǎn)集。首先，逐列掃描每個向下開口點(diǎn)簇，找到每列的上端點(diǎn)PU=(x0，y1)和下端點(diǎn)PL=(x0，y2)；然后得到各列的擬合點(diǎn)PF=(x0，(y1+y2)/2)。如圖7所示，擬合點(diǎn)用黃色標(biāo)記。

圖7 擬合點(diǎn)（黃色標(biāo)記）Fig.7 Fitting points（tag yellow）

擬合點(diǎn)集定義為PF={(xi，yi)∣1 ≤i≤n}，擬合問題建模為約束最小二乘問題。在本文擬合算法中，雙曲線方程可以表示為

該雙曲線標(biāo)準(zhǔn)方程可以轉(zhuǎn)化為一般方程

式（9）可以簡單地表示為

從一點(diǎn)到曲線的代數(shù)距離定義為

曲線到點(diǎn)簇的代數(shù)距離的平方和為

由于所有擬合點(diǎn)都需要在擬合雙曲線的中心(x0，y0)以下，且擬合雙曲線的中值需要位于雙曲線點(diǎn)簇的開放范圍[S，E]內(nèi)，因此需要增加兩個約束條件，表示為

根據(jù)式（10）、式（12）、式（14）、式（15），可以將雙曲線擬合算法表示為凸優(yōu)化問題：

其中≤表示兩個向量之間的元素小于或等于

該凸優(yōu)化問題也是一個約束線性最小二乘問題，可以用內(nèi)點(diǎn)凸四元規(guī)劃算法［24］來求解。利用該擬合算法，可以初步擬合OSCA得到的下開口點(diǎn)集。

5 基于擬合誤差的剔除方法

采用雙曲線擬合算法對所有下開口點(diǎn)集進(jìn)行擬合后，擬合點(diǎn)簇(1 ≤i≤N)（N＞2）的擬合誤差定義為

其中Di為第i個擬合點(diǎn)簇的式（13）矩陣，xi為第i個雙曲線擬合問題的解。OSCA 的輸出點(diǎn)集可以分為非雙曲點(diǎn)集和雙曲點(diǎn)集。為了區(qū)分這兩類點(diǎn)簇，定義了非雙曲點(diǎn)簇的判斷標(biāo)準(zhǔn)為

如圖8 所示，由于非雙曲點(diǎn)集不具有完全的雙曲線特征，非雙曲點(diǎn)集的擬合誤差會比雙曲點(diǎn)集的擬合誤差大得多，應(yīng)用這種消去方法可以消去大多數(shù)非雙曲點(diǎn)集。

圖8 三個非雙曲點(diǎn)集的像素表示Fig.8 Pixel representation of three non-hyperbolic point clusters

6 合成數(shù)據(jù)和實(shí)測數(shù)據(jù)實(shí)驗(yàn)結(jié)果

在本節(jié)中，展示和分析了在合成數(shù)據(jù)集和實(shí)測數(shù)據(jù)集上的實(shí)驗(yàn)結(jié)果。本節(jié)最后對本文提出的非雙曲線點(diǎn)集消除方法與［22］中提出的PFJ方法進(jìn)行了比較研究。

6.1 合成數(shù)據(jù)實(shí)驗(yàn)

利用電磁模擬器gprMax［25］生成模擬探地雷達(dá)B掃描圖像。該開源軟件允許用戶通過設(shè)置不同的相對介電常數(shù)和電導(dǎo)率來指定不同的介質(zhì)，并在不同的場景中設(shè)置不同的參數(shù)來指定不同的埋藏目標(biāo)的大小。基于單層鋼筋網(wǎng)格場景，采用gprMax 軟件生成綜合數(shù)據(jù)。如圖9所示，6根鋼筋埋在gprMax定義的土壤介質(zhì)中。發(fā)射天線和接收天線靠近分界面，并以步長Δx=0.8 m的方式從x=0.2 m移動到x=1.8 m。

圖9 仿真單層鋼筋網(wǎng)格場景的剖視圖Fig.9 Sectional view of the synthetic single layer rebar mesh scene

圖10（a）和圖10（b）分別為均值對消后的原始B 掃描圖像和預(yù)處理后的圖像。OSCA 和雙曲線擬合算法的結(jié)果如圖10（c）和12（d）所示。由于介質(zhì)的異質(zhì)性，圖10（a）原始B 掃描圖像中存在一些雜波。在沒有其他物理干擾的情況下，圖10（c）的OSCA 結(jié)果中沒有出現(xiàn)非雙曲點(diǎn)簇。在本次合成數(shù)據(jù)實(shí)驗(yàn)中，本文提出的B 掃描處理模型自動擬合并檢測出所有正確的雙曲點(diǎn)集。合成數(shù)據(jù)集實(shí)驗(yàn)的具體結(jié)果如表4所示。

表4 合成數(shù)據(jù)實(shí)驗(yàn)結(jié)果Tab.4 Obtained results on synthetic data sets

圖10 所提方法處理仿真數(shù)據(jù)結(jié)果Fig.10 Images of the synthetic data processed by the proposed model

6.2 實(shí)測數(shù)據(jù)實(shí)驗(yàn)

數(shù)據(jù)集采集由具有1500 MHz 天線的LTD-2600探地雷達(dá)設(shè)備進(jìn)行，用于檢測圖11中修復(fù)水泥路面下的鋼筋網(wǎng)。鋼筋網(wǎng)由半徑6 mm，相鄰距離40 cm的鋼筋組成，埋在地表下10 cm處。

圖11 GPR設(shè)備和數(shù)據(jù)收集Fig.11 GPR equipment and data collection

如圖12 所示，（a）和（f）分別為數(shù)據(jù)1 和數(shù)據(jù)2經(jīng)過均值對消之后的圖像結(jié)果。通過預(yù)處理方法，將B 掃描圖像轉(zhuǎn)化為包含下開口點(diǎn)集和一些不規(guī)則區(qū)域的二值圖像，如圖12（b）和（g）所示。然后通過OSCA 去除不規(guī)則區(qū)域，識別出所有下開口點(diǎn)集，但保留了一些具有下開口特征的非雙曲點(diǎn)集，如圖12（c）和（h）所示。隨后利用FEAA 對OSCA 結(jié)果進(jìn)行處理，剔除了非雙曲線點(diǎn)集，保留了所有目標(biāo)雙曲點(diǎn)集，如圖12（d）和（i）所示。最后進(jìn)行雙曲線擬合得到結(jié)果，如圖12（e）和（j）所示。在實(shí)測數(shù)據(jù)集上得到的結(jié)果如表5 所示。圖13 為在進(jìn)行FEAA方法處理時，各擬合點(diǎn)集的誤差與所有點(diǎn)集誤差均值的2 倍比較結(jié)果。在圖13（a）中，展示了圖12（c）中（數(shù)據(jù)1）中9 個點(diǎn)集的擬合誤差，從圖中可以看出，第一個點(diǎn)集的誤差大于誤差均值的2 倍，因此相比于圖12（c），在圖12（d）中消除了第一個非雙曲點(diǎn)集。在同樣的情況下，在圖13（b）中，展示了圖12（h）中（數(shù)據(jù)2）中6 個點(diǎn)集的擬合誤差，從圖中可以看出，第一和第二點(diǎn)集的誤差大于誤差均值的2 倍，因此相比于圖12（h），在圖12（i）中剔除了這兩個非雙曲點(diǎn)集。

表5 實(shí)測數(shù)據(jù)集實(shí)驗(yàn)結(jié)果Tab.5 Obtained results on real datasets

圖12 由所提方法處理實(shí)測數(shù)據(jù)集結(jié)果Fig.12 Images of real datasets processed by proposed model

圖13 各擬合點(diǎn)集誤差與所有點(diǎn)集誤差均值的2倍比較結(jié)果Fig.13 Comparison between the error of each fitting point cluster and double mean of errors of all point cluster

為了進(jìn)一步驗(yàn)證FEE 方法的有效性，進(jìn)行了PFJ 與FEE 方法的對比實(shí)驗(yàn)。圖14（a）為進(jìn)行拋物線擬合的結(jié)果，且式（21）為數(shù)據(jù)集1 中非雙曲點(diǎn)集的擬合拋物線方程。根據(jù)PFJ 方法，該非雙曲點(diǎn)集無法消除，如圖14（b）所示。而這種非雙曲點(diǎn)集可以通過FEE 方法消除，如圖12（c）和12（d）所示。FEE方法和PFJ得到的雙曲線點(diǎn)集數(shù)如表6所示。

表6 FEE和PFJ方法得到的雙曲線點(diǎn)集數(shù)Tab.6 Number of obtained hyperbolic clusters of FEE and PFJ

圖14 PFJ處理數(shù)據(jù)集1結(jié)果Fig.14 Images of dataset 1 processed by PFJ

y=-0.0027x2+0.1820x+710.3430 （21）

7 結(jié)論

本文提出了一種基于非雙曲點(diǎn)集消除方法的魯棒集成模型：基于擬合誤差的消除（FEE）方法，用于自動識別和擬合探地雷達(dá)（GPR）B 掃描圖像中的雙曲線。在我們的模型中，首先對B 掃描圖像進(jìn)行閾值分割開閉操作處理，然后采用OSCA 方法對下開口點(diǎn)集進(jìn)行分離，隨后采用基于代數(shù)距離的雙曲線擬合算法。最后，利用FEE 方法可以對所有的雙曲點(diǎn)集進(jìn)行檢測和擬合。綜合仿真和實(shí)測數(shù)據(jù)實(shí)驗(yàn)結(jié)果，驗(yàn)證了所提模型的有效性和魯棒性。對于所提出的方法，待改進(jìn)部分為雙曲線識別部分，可以結(jié)合深度學(xué)習(xí)領(lǐng)域的識別方法以提高方法的魯棒性。此外，可以采用先利用深度學(xué)習(xí)方法進(jìn)行檢測，后進(jìn)行擬合的方法思路。