基于改進馬爾可夫隨機場探地雷達有效信號提取方法

吳學禮,宋凱,史思遠,栗玉潔

(1.河北科技大學電氣工程學院,石家莊 050018; 2.河北省生產過程自動化工程技術研究中心,石家莊 050018)

探地雷達的探測原理是向地下發射高頻電磁波,經過地各類復雜地質的反射,接收機接收反射回波,形成地下情況資料。

在工程項目中,使用探地雷達對道路、涵洞、大壩進行工程質量探測時,因為工程環境的復雜,同時受到探地雷達自身系統結構的影響,回波信號中會包含各類噪聲,包括地表直達波和各類高頻雜波等。這些噪聲會將有效信號淹沒,導致成像效果較差,有效信息識別困難,可能產生誤判[1]。

對探地雷達數據進行雜波抑制和有效信號提取一直是中外學者研究熱點。空域方面主流方法包括均值濾波、中值濾波[2]、對稱近鄰平滑濾波(symmetric nearest neighbor,SNN)和K鄰近平滑濾波[3](Knearest neighbor,KNN)等降噪方法。均值濾波是一種線性濾波,該點像素值由周圍的像素均值代替,該方法簡單、快速,但數據圖像容易模糊,波形細節不清晰。中值濾波是一種非線性濾波,該點像素值由周圍像素的中間值代替,該方法濾波效果較好,但存在有效信號邊緣不準確的問題。對稱近鄰平滑濾波和K鄰近平滑濾波等降噪方法在降噪和邊緣細節保存方面效果較好,但對于噪聲程度較高的探地雷達數據情況效果一般。在頻域方面常見的方法有拉普拉斯變換、離散余弦變換(discrete cosine transform,DTC)和小波閾值濾波[4-5]等。前兩者主要通過時頻域的轉換,濾除高頻噪聲,然后進行逆變換。兩種方法對于高頻高幅噪聲濾除效果較好,但對低頻噪聲濾除效果不明顯。后者利用小波變換將信號分解,濾除噪聲小波系數,然后對圖像進行重構。該方法在有效信號提取和濾波方面效果均優,但是對于小波的選取有針對性,選取不同的小波效果差異較大。近期研究熱點機器學習和深度學習算法方面,有稀疏字典學習模型[6]以及卷積神經網絡算法[7]等方法,此類方法在降噪效果和信號邊緣保持方面效果均佳,但需要數量龐大的數據集進行訓練,但探地雷達實地數據較少,正演仿真速度較慢,導致學習效果不理想。

上述各種方法在雜波抑制效果和有效信號細節的保存方面各有優劣。面對過度濾波易造成有效信號細節保存不完整和過度關注細節導致雜波濾除效果差的這一主要矛盾點,首次引進的馬爾可夫隨機場模型進行探地雷達雜波抑制和有效信號提取。馬爾可夫隨機場數學理論由俄國數學家安德烈·馬爾可夫提出。后由英國統計學家Julian Besag將其應用在圖像分割和圖像降噪領域。原始馬爾可夫隨機場模型對探地雷達數據背景雜波抑制和有效信息的提取均有一定適用性,但未針對探地雷達數據特性構造能量函數,降噪后回波波形連續性不強,細節保存不完善。

針對上述問題,現改進馬爾可夫隨機場根據探地雷達數據特點對能量函數進行針對性構造,以增強抑噪和有效信號提取能力。算法使用大津算法結合改進馬爾可夫隨機場模型對探地雷達圖像進行降噪。大津算法計算出分割閾值,將噪聲和有效信號提取為二值圖像。使用改進馬爾可夫隨機場模型結合ICM迭代算法對二值化圖像進行降噪。然后根據探地雷達圖像特點,使用中值消去法濾除探地雷達地表直達波。最后使用仿真數據和現場剖面數據對算法的有效性和實用性進行驗證。

1 基于馬爾可夫隨機場雜波抑制

1.1 大津算法的閾值計算

大津算法(Oisu algorithm,OTSU)是一種對灰度圖像進行二值化處理的算法[8]。算法原理是通過最大化類間方差的方法求出其的劃分閾值,依次對閾值迭代找出圖像中方差最大的A、B兩個部分,最大限度地將圖像信息分為前景和背景,以此將探地雷達有效信號和雜波整體劃分出來,便于下一步濾波操作。

灰色圖像色域為[0,255],探地雷達圖像總像素是N個,灰度值是i的像素有ni個,迭代閾值記作k,圖像中灰度是i的概率記作pi,圖像平均灰度值記作mG(k),則

(1)

(2)

像素在A部分的概率記作pA(k),A部分的平均灰度記作mA(k),則

(3)

(4)

像素在B部分的概率記作pB(k),B部分的平均灰度是mB(k),則

pB(k)=1-pA(k)

(5)

(6)

根據數學關系推導可得

pA(k)mA(k)+pB(k)mB(k)=mG

(7)

pA(k)+pB(k)=1

(8)

求得A部分、B部分類間方差為

σ(k)2=pA(k)[mA(k)-mG(k)]2+

pB(k)[mB(k)-mG(k)]2

(9)

將式(7)和式(8)代入式(9),可得大津算法中類間方差公式為

(10)

根據式(10),依次進行迭代,當σ(k)值最大時,獲取最佳閾值,可以最大限度分割圖像。

k=argmax{σ(k)}

(11)

1.2 馬爾可夫隨機場模型搭建

無向圖模型是使用無向圖表示變量結點之間相關關系的一種模型。如果該模型的聯合概率分布滿足以下3個條件,則稱該聯合概率分布為馬爾可夫隨機場。

第一是成對馬爾可夫性,其計算公式為

P(Yu,Yv|Yo)=P(Yu|Yo)P(Yv|Yo)

(12)

式(12)中:u、v為兩個無邊直接連接的結點;Yu、Yv分別為其隨機條件;o為兩結點所有連接路徑上結點;Yo為其隨機條件;P(Yu,Yv|Yo)為在確定Yo下Yu、Yv同時發生的概率;P(Yu|Yo)和P(Yv|Yo)為在確定Yo下Yu、Yv分別發生的概率。

第二是局部馬爾可夫性,其計算公式為

P(Yv,Yo|Yw)=P(Yv|Yw)P(Yo|Yw)

(13)

式(13)中:w、v為所有直接連接的結點;Yw、Yv分別是其隨機條件;o為除w、v以外的所有結點;Yo為其隨機條件;P(Yv,Yo|Yw)為在確定Yw下Yv、Yo同時發生的概率;P(Yv|Yw)和P(Yo|Yw)為在確定Yo下Yu、Yv分別發生的概率。

第三是全局馬爾可夫性,其計算公式為

P(YA,YB|YC)=P(YA|YC)P(YB|YC)

(14)

式(14)中:A、B為兩個無邊直接連接的結點集合;YA、YB分別為其隨機條件;C為兩結點集合所有連接路徑上結點集合;YC為其隨機條件P(YA,YB|YC)是在確定YC下YA、YB同時發生的概率;P(YA|YC)和P(YB|YC)分別為在確定YC下YA、YB分別發生的概率。

使用大津法處理過的探地雷達二值化圖像包含有效信號和噪聲信息。根據探地雷達數據特點,有效信號較為明顯,噪聲含量相較有效信號不高,所以對于圖像像素來說,含噪圖像與降噪圖像有較強的相關性。同時,探地雷達噪聲特點是離散、分散的,有效信號的特點是相對集中、連續,所以推斷出像素點與自身周圍的像素點同樣有較強的相關性。根據上述特點,探地雷達數據符合馬爾可夫隨機場模型條件。

根據探地雷達圖像特點構建馬爾可夫隨機場模型進行噪聲濾除。設輸出降噪圖像是X,根據馬爾可夫隨機場性質,像素點只與周圍有直接邊連接的像素點相關,與其余的像素點條件獨立,故每個像素點xi,j與相鄰的同行和同列像素點xi±1,j和xi,j±1有強相關性,可以組成一個團{xi,j,xi±1,j,xi,j±1}。設可觀測圖像即帶噪的二值化圖像是Y,其圖像的每個像素點yi,j與xi,j同樣有較強的關聯性,可組成一個團{xi,j,yi,j}。同時,輸出像素點xi,j與其所在九宮格四角像素xi±1,j±1存在相對較弱關聯性,可組成團{xi,j,xi±1,j±1}。將各個團組合在一起可組成最大團,以此構建馬爾可夫隨機場模型,模型如圖1所示。

i、j分別是圖像像素點對應的行和列,作為像素點的索引

根據Hammersley-Clifford定理[9],聯合概率分布P(X,Y)等于表示構建模型上最大團上勢函數[10]的乘積,即

(15)

式(15)中:Q為構建的無向圖中最大團;XC、YC為團C分別對應的隨件變量;Z為歸一化常量,確保P(X,Y)不大1,處于概率分布的范疇,表達式為

(16)

P(X,Y)是概率分布,值一定不小于0,所以勢函數ΨC(XC,YC)是嚴格的正函數,故使用指數形式為

ΨC(XC,YC)=e-E(XC,YC)

(17)

式(17)中:E(XC,YC)中為一個定義在變量(XC,YC)上的實值函數,也稱為能量函數。

根據探地雷達圖像像素之間特點構造能量函數形式為

(18)

(19)

(20)

E(X,Y)=A(X,Y)+B(X,Y)+C(X,Y)

(21)

式中:A(X,Y)、B(X,Y)、C(X,Y)分別為各團構造的能量函數;M、N分別為原圖像X的行像素數和列像素數;α、β和η為各團能量函數的學習參數,均是正常數;α參數是無向圖中輸出圖像素點與四周直接相鄰像素點之間的相互關系的偏好量;β參數是無向圖中輸出圖像像素點與觀測圖像像素點之間的相互關系偏好量;η參數是無向圖中輸出圖像像素點與其所在九宮格四角像素點之間的相互關系偏好量。

1.3 ICM算法濾波

輸出降噪圖像為X,令初始X=Y,代入模型中,聯合概率分布最大值時求得每個像素點的值。

xi,j=arg max{P(X,Y)}

(22)

聯合概率分布P(X,Y)是勢函數ΨC(XC,YC)的乘積,從式(16)～式(18)可知,將E(XC,YC)能量函數相加獲得。E(XC,YC)能量函數值越小,聯合概率分布P(X,Y)值越大。所以可以將式(19)轉化能量函數最小值求每個像素點值,即

xi,j=arg min{E(XC,YC)}

(23)

圖像X是二值化圖像,當像素點xi,j被認定為噪聲時,其值會被設定為1,被濾入白色背景中。當被認定為有效信號時,其值被設定為-1,作為黑色像素點進行保留。然后對構建的馬爾可夫隨機場模型進行ICM算法[11]迭代,初始輸出值令X=Y。對xi,j分別賦值-1和1,對其分別計算出的能量函數的值進行比較,取較小值者。然后對下一個像素點進行相同計算取值,圖像X所有像素點至少迭代一次。ICM算法流程如圖2所示。

圖2 ICM算法

1.4 結合二值圖像對圖像降噪

通過帶噪二值圖像與降噪后的二值化圖像對比,對探地雷達圖像進行中值濾波。

(24)

式(24)中:pi,j為輸出濾波后圖像的像素點;gi,j為原始帶噪探地雷達圖像像素點;xi,j、yi,j分別為經過馬爾科夫隨機場模型濾波前后的像素點。對于當前位置對應的像素點,當該位置濾波前后的二值化圖像像素數值相同時,pi,j取值等于gi,j的值。當該位置濾波前后的二值化圖像像素數值不相同時,pi,j取值是該像素點九宮格內其余個點的中間值。通過該方式進行濾波,可以較大程度保留圖像的有效信號波形,只對噪聲進行抑制。

1.5 算法流程

原始圖像經大津算法分割,使用改進馬爾科夫隨機場模型生成去噪二值圖像,然后結合二值圖像使用中值定理對雷達數據進行降噪,最后利用均值法去除地表直達波,完成對原始數據的噪聲去除和有效信號提取工作,算法流程如圖3所示。

圖3 算法流程圖

本文算法在MATLAB2020b軟件中運行仿真,電腦CPU是i5,主頻2.66 GHz,程序大小為9 KB,運行時間為1.7 s。改進馬爾可夫隨機場算法是本算法中是時間復雜度最高的部分,復雜度是O(nlogn),適用于在探地雷達硬件中部署。

2 實驗與分析

探地雷達可以對大壩、道路等工程項目質量進行檢測,查看其中是否有漏洞、空洞以及鋼筋是否錯位等影響工程質量的缺陷[12]。使用GPRMAX軟件進行建模以及正演計算可展示探地雷達在工程項目中的使用效果。工程項目中所使用的材料介質不同,其電磁學性質有很大差異。根據麥克斯韋方程,雷達波穿越不同的介質,其反射速度和衰減率不同,導致回波的波形有不同的表現。這是探地雷達探明地下情況的工作原理。主要工程介質的電磁參數如表1[13]所示。

表1 主要工程介質的電磁參數

設置一個寬度2 500 mm,厚度是450 mm的帶筋混凝土板。板內有3根直徑是40 mm的鋼筋,形心間距均是1 000 mm,埋深均是225 mm。板內隨機位置設置一處形狀為三角形,介質為空氣的缺陷,用于模擬滲漏點,建模圖像如圖4所示。中心天線使用ricker源進行模擬,頻率是900 MHz,位于混凝土板上方500 mm處,正演計算仿真圖像如圖5所示。

圖4 模型圖

圖5 正演計算仿真圖

對正演計算仿真圖添加椒鹽噪點以模擬工程現場實際情況,含噪圖像如圖6所示。

圖6 含噪圖

使用改進馬爾可夫隨機場濾波法進行降噪,首先通過大津算法生成二值圖,劃分出前景和背景。前景包括有效信號、噪點、地表直達波,帶噪二值化圖像如圖7所示。

圖7 帶噪二值圖

使用改進馬爾可夫隨機場濾波法結合ICM迭代濾波算法,對帶噪二值圖像進行降噪。根據地表直達波的特點將其濾除,得到去噪二值圖像,如圖8所示。

將去噪二值圖提取的有效信號與原始帶噪圖像進行對比結合濾波,得到降噪后的圖像。同時使用均值濾波法、對稱近鄰平滑濾波法、K鄰近平滑濾波法、中值濾波法、馬爾可夫隨機場濾波法等方法對比降噪效果。降噪對比仿真圖像如圖9所示。

圖9 不同方法濾波效果對比

通過觀察圖像對比不同方法濾波效果可以發現:①6種方法均對雜波的抑制有一定的效果;②均值濾波法、稱近鄰平滑濾波法和K鄰近平滑濾波法濾波后的圖像中可以觀察到較為明顯的噪點;③中值濾波圖像無噪點,但因無提取有效信號過程,波形細節保存不完整。馬爾可夫隨機場濾波法對雜波的抑制較好,背景雜波幾乎均被濾除,但有效信號在馬爾可夫隨機場濾波中信息保存不足,導致波形邊緣準確度較差;④改進馬爾可夫隨機場濾波法因在濾波過程中提取出有效信號完整,故與原始數據信號偏差更低,更好地保留了有效信號的信息。

隨機選取原始正演圖像、中值濾波圖像、均值濾波圖像、對稱近鄰平滑濾波法、K鄰近平滑濾波法、馬爾可夫隨機場濾波法和改進馬爾可夫隨機場濾波圖像中同一單道波,對波形進行分析,對比降噪效果。單道波圖對比如圖10所示。

圖10 單道波形對比圖

由圖10(a)單道波形對比可以看出,均值濾波法波形不平滑,存在突變,前者突變角度多為銳角,后者突變因濾去高頻雜波的原因,突變角度更大。兩者與原圖像波形相似度最低說明濾波效果最差。從圖10(b)單道波形細節對比可以看出,對稱近鄰平滑濾波法、K鄰近平滑濾波法波形與原圖像數值差別不大,但是準確度不高,導致圖像擬合度較低。改進馬爾可夫隨機場濾波法相較于中值濾波法、馬爾可夫隨機場濾波法,與原始圖像單道波重合度更高,說明對雜波的抑制更好,對有效信號的細節保留更準確,提取有效信號更完整。

峰值信噪比[14](peak signal to noise ratio,PSNR)和結構相似度(structural similarity,SSIM)兩個評價標準,可以從整體空域角度對降噪效果進行評價。峰值信噪比是通過降噪后圖像與原始圖像像素間的均方誤差通過計算得到的結果,其結果越大,表示降噪后的圖像相較于初始圖像還原度越高,對雜波的抑制效果越好。計算公式為

(25)

式(25)中:Xi,j和Yi,j分別為探地雷達原始圖像和降噪后的圖像,圖像像素數量是MN個。

結構相似度是衡量兩個圖像相似程度評價方法。其包括亮度函數、對比度函數和結構函數3個評價指標。計算結果數值越大,兩個圖像相似度越高,表明圖像降噪效果較好,有效信號細節保存準確完整。公式為

(26)

從峰值信噪比、結構相似度2個方面衡量6種方法的降噪效果如表2所示。

表2 6種方法濾波效果對比

根據表2展示結果可以看出,與觀察結果一致,均值濾波、對稱近鄰平滑濾波法、K鄰近平滑濾波法的峰值信噪比和結構相似度數值均較低,說明降噪效果和有效信號細節保存結果不理想。中值濾波法效果結構相似度數值、峰值信噪數值較高,但低于馬爾可夫隨機場濾波,說明馬爾可夫隨機場提取有信號后濾波對波形細節完整度保存更好。本文提出濾波方法兩項評價指標較其他方法數值最高,峰值信噪比達52.527 1 dB,說明降噪效果在各類方法中最好。同時與原圖像結構相似度為0.998 1,與原始數據相似度最高,說明有效信號保存更完整。同時結合觀察法和單道波形對比,說明該方法改進效果明顯,在濾除背景雜波和有效信號提取的完整性有較好的效果,對馬爾可夫隨機場降噪改進效果明顯。

根據圖像特點,選用中值法濾除地表直達波,同時對圖像進行細節增強處理,使得有效信號表現更清晰,最終濾波圖像如圖11所示。

圖11 最終濾波圖像

3 工程實例

使用中國電波傳播研究所研制的LTD-60道路綜合檢測系統探地雷達實測數據進行算法有效性測試。圖像是在對某城市主干道進行安全檢測時得到實測剖面圖,如圖12所示。經過改進馬爾可夫隨機場濾波法進行背景降噪的二值圖像,如圖13所示,可以從二值圖中較為清晰地發現道路下邊的缺陷[15]。選取降噪效果較好的濾波圖像對實測圖像進行效果比對,如圖14所示。同中值濾波、馬爾可夫隨機場濾波法相比,使用本文提出的濾波方法圖像背景噪聲和地表直達波得到有效抑制,圖像細節表現相比其他方法更加清晰,有利于工作人員對道路基礎病害的發現和辨別,表明該方法對實測數據處理有較好的效果,在工程探測項目中有一定的實踐價值。

圖12 實測剖面圖

圖13 實測剖面二值圖

圖14 不同濾波方法實測剖面圖對比

4 結論

探地雷達圖像在工程項目實踐中,易受到大量噪聲污染,造成有效信號被淹沒,造成有效信號識別、提取困難。本文研究提出了一種基于改進馬爾可夫隨機場探地雷達有效信號處理方法。利用大津法將有效信號、噪聲和背景分離。運用改進馬爾可夫隨機場原理結合ICM迭代法對圖像進行降噪,后利用濾波后的二值化圖像與原始圖像結合、濾波,以達到降噪效果。根據實驗結果得出以下結論。

(1)本文方法與各類經典濾波法相比,在各項評價指標中評分最高。同馬爾可夫隨機場濾波法相比,在峰值信噪高出5 dB以上、結構相似度高出0.04以上。表明改進馬爾可夫隨機場濾波法抑噪效果更好,與原始圖像相似度更高,相較于馬爾可夫隨機場濾波,改進效果明顯。

(2)本文方法在對實測剖面圖進行雜波抑制和有效信號提取中被證明有效。在單道波形對比中與原始圖像有更好的擬合度,說明有效信號提取準確,濾波后的二值圖像,可以更好地輔助辨識的圖像波形,幫助專業人員更加準確地分析判斷工程缺陷。