基于非負矩陣分解的GPR高頻雜波抑制

苗 翠，原 達，王冬雨，李文生

山東省高校智能信息處理重點實驗室（山東工商學院），山東 煙臺 264005

探地雷達（GPR）是一種常用的探地無損檢測技術[1]，對于地下埋藏物及地下結構的探索具有廣泛的應用[2-4]。在采集到的探地雷達信號中，目標的反射頻率往往是較強的，雜波可以定義為那些與目標散射特性無關但發生在同一采樣時間窗中且具有與目標波長相似的光譜特性的信號[5]。當部分雜波反射更高強度能量并掩蓋目標反射波時，稱這類雜波為高頻雜波，在圖像中表現為較高的像素強度。其來源可能是由發射天線和接收天線之間的穿透以及天線和地面之間的多次反射引起的，地面特性阻抗的局部變化以及材料中包含的一組小反射源也會引起高頻雜波。尤其對于淺層目標的檢測，當目標埋藏物越淺時，使用的探測頻率越大，這時所引起的高頻干擾也更為復雜。如何從復雜的GPR反射波中去除高頻雜波成分保留目標信號，是需要研究的重要課題。

對于探地雷達信號的雜波處理問題，基于主成分分析（principal component analysis，PCA）[6]及魯棒主成分分析（robust principal component analysis，RPCA）[7]和奇異值分解（singular value decomposition，SVD）[8]等子空間技術往往表現出較好的性能，該類方法將GPR圖像分解為對應于雜波和目標的分量，將主分量作為雜波，其余成分作為目標信號來進行處理，但該類方法在GPR圖像含有多條雜波或多個目標信號的情況下雜波不能完全去除。近兩年提出的形態學方法（morphological component analysis，MCA）[9]將GPR 圖像的雜波和目標分量分別用曲波curvelet 和非采樣離散小波變換（unsampled discrete wavelet transform，UDWT）字典稀疏表示來進行雜波抑制，雖然視覺效果表現不錯，但算法復雜度較高，處理速度不夠快。基于多分辨率的方法在GPR 去雜波的處理中也取得了較好的效果，如多尺度雙邊濾波（multiscale directional bilateral filter，MDBF）[10]借助于雜波和目標在反射波形中呈現出的幾何差異進行多方向多尺度的分解。集成經驗模式分解（ensemble empirical mode decomposition，EEMD）將探地雷達信號分解為一系列模函數（intrinsic mode functions，IMFs）并計算其置換熵（permutation entropy，PE），通過設置相關閾值來進行噪聲和目標的區分，能有效提高目標的分辨率[11]，不過圖像質量難以保證。基于徑向基函數（radial basis function，RBF）神經網絡的方法也同樣被用于GPR 圖像的雜波濾除，使用零偏移格林函數作為訓練數據的期望輸出，結果能提高探地雷達的垂直分辨率[12]，但對于目標成分的保留不夠清晰。

非負矩陣分解（non-negative matrix factorization，NMF）由Lee和Seung[13]于1999年在自然雜志上提出，該算法與與PCA、ICA（independent component analysis）、SVD 等子空間分解技術不同的是，NMF 中對分解矩陣的非負約束能夠使數據進行稀疏表示。分解的結果中不含有負值，使數據更具有可解釋性和可使用性。并且由于其簡便高效的特點，已經被廣泛應用于數據降維、特征選擇、圖像融合、文本聚類、語音增強等方面[14-17]。對于NMF 在GPR 圖像雜波處理的應用中[18]，將GPR 圖像自身的像素強度值來進行NMF 分解，去掉第一主成分來進行雜波抑制。而后又提出了魯棒非負矩陣分解（robust nonnegative matrix factorization，RNMF）[19]，在NMF 基礎上加入誤差矩陣進行分解，具有更好的實時性。在本文中，基于NMF 提出了概率非負矩陣分解（PNMF）應用于GPR雜波抑制，選取圖像的水平梯度作為非負矩陣分解的輸入，對于參數后驗分布的求取，使用變分貝葉斯方法進行推理。變分貝葉斯是貝葉斯推理中的一種近似推理算法，對于概率模型中參數的求取具有速度快、穩定性好的特點，與傳統的NMF 相比，具有更高的信噪比和更好的魯棒性。

1 基于概率模型的非負矩陣算法

1.1 構建概率NMF模型

在NMF處理GPR圖像時，像素強度矩陣F看作是目標成分Ftarget和雜波成分Fclutter的組合結構：

將矩陣F進行非負矩陣分解得到：

其中圖像大小為M×N,F∈RM×N,W∈RM×K為基矩陣，H∈RK×N為系數矩陣，且W和H都具有非負約束條件，即s.t.W,H≥0。高頻雜波的強度很大，一般選取前K個分量對應于雜波成分，K

其中D(?)表示距離度量函數，用來度量F分解為WH后與F的距離。為了衡量矩陣分解前后的相似性，使用KL 散度（Kullback-Leibler divergence）作為度量指標進行計算，此時的參數W、H的更新規則如下：

在基矩陣W和系數矩陣H的更新中分母1表示的是單位矩陣。上述使用拉格朗日乘子法進行非負矩陣求解時通常具有簡單快速的特點，這種非概率的方法容易導致過擬合，不具有穩定性。而貝葉斯方法通過定義滿足非負約束的先驗分布，再結合實際觀察的數據進行后驗推理，能大大減少過擬合的現象，提高非負矩陣分解算法的性能。在概率模型下的NMF 中，將待分解的潛在基矩陣W和系數矩陣H看作是兩個隨機變量，數據F中的值來自于W和H的乘積，以及一些高斯噪聲E，如下表達式：

1.2 基于變分貝葉斯分解非負矩陣

在進行變分貝葉斯推理之前，對變量的先驗分布P(θ)進行假設，由于PNMF 的非負性約束，設W和H的先驗服從指數分布：其中，ε(x;λ)=λexp(-λx)u(x)是指數分布的密度函數，λk>0，u(x)為單位階躍函數。

噪聲方差的先驗設為逆伽馬分布：

在PNMF 中，參數的后驗分布使用變分貝葉斯（variational Bayes，VB）進行推理，通過求解一個近似的后驗分布去逼近真實的概率分布。在復雜的多變量參數θ情況下，VB將其分解為一組相互獨立的變量θi，表示如下：

其中，q(θ)對應于參數θ的分布。根據文獻[20]，得到后驗分布：

式（11）中，G(?)為伽馬分布，密度函數表達式見式（9），TN(?)為截斷高斯分布，概率密度函數表達式如下：

在截斷高斯分布中，x<0 時為密度為0 的高斯分布，其中φ(?)是服從N(0,1)分布的累積分布函數。得到W和H以及相關參數的后驗分布后，將前K項作為高頻雜波對應的分量并進行如下分量表示：

其中K值代表矩陣的秩，通過該值可以得到雜波對應的低秩矩陣。對于高頻雜波，前K項成分分量對雜波的貢獻最大，K值的選擇在實驗部分進行了具體論證。

2 實驗

2.1 實驗數據背景

本文實驗分別采用GPR模擬數據和應用環境采集的真實數據進行方法驗證分析，并通過不同算法的信噪比及視覺效果來驗證本文算法的有效性及魯棒性。

2.2 模擬GPR數據

模擬數據通過gprMax軟件[21]的時域有限差分方法（FDTD）進行模擬獲得。在模擬過程中，將鋁盤作為埋藏物，然后分別放置于3 種不同特性（干性土壤dry_sand、潮濕土壤damp_sand、濕潤土壤wet_sand）的土壤環境下，埋藏物和3 種土壤的介質參數來自于文獻[22]，其介電常數和導電率詳見表1。

表1 材料的電磁特性Table 1 Electromagnetic properties of materials

獲得的GPR 圖像如圖1 中的A1～A3 所示，對應于鋁盤分別埋藏于干性土壤dry_sand、潮濕土壤damp_sand和濕潤土壤wet_sand環境下的反射結果，可以從圖中清楚的看到，它們的目標反射波形隨著環境的改變而發生位置高度和向下開口大小的變化，但是都擁有相同特性的水平雜波。在使用本文算法PNMF進行實驗時，首先使用圖像自身的像素強度代入計算，得到圖1 中B1～B3的處理結果，可以看到該水平雜波并沒有完全去除，并且B1、B2 中的雙曲波目標頂點部分也被削弱了。當通過圖像像素值的水平梯度進行PNMF計算時，結果對應于圖1 中的C1～C3，可以看出，處理結果不僅可以很好的保留雙曲波信號，還對雜波進行了明顯的抑制。

圖1 不同PNMF對3幅模擬GPR圖像的噪聲抑制效果對比Fig.1 Comparison of noise suppression effects of different PNMF on three simulated GPR images

峰值信噪比（PSNR）是圖像處理中最為普遍的一種客觀評價指標，基于純凈參考圖像I與處理后圖像I?之間對應像素點(i,j)的誤差來判斷圖像處理的質量好壞。通常通過均方差（MSE）進行定義：

其中M×N為圖像I的大小，在模擬GPR數據中，無雜波的參考圖像I通過含雜波的目標反射圖像和只含雜波的圖像作差來得到。在使用概率模型進行非負矩陣分解時，分別使用圖像像素強度和水平梯度作為輸入，通過PSNR 值的計算得到一組對比柱狀圖，如圖2 所示。可以直觀看到在使用梯度進行PNMF 探地雷達圖像雜波處理時，得到的PSNR值更高。

圖2 像素PNMF和梯度PNMF的PSNR對比Fig.2 PSNR comparison of pixel PNMF and gradient PNMF

由視覺質量和PSNR 對比得到，在PNMF 算法對GPR圖像的去雜波過程中，選擇其水平梯度作為輸入矩陣進行處理效果更佳（后續的PNMF 指的都是梯度PNMF）。

在GPR 圖像雜波處理中，雜波成分一般存在于矩陣分解后的前幾個主分量中，由此分別選取K=1，2，3，4，5進行處理后PSNR值的比較，結果如表2所示。數據結果顯示，干土壤中鋁盤反射的GPR 圖像在K取2 的時候去雜波后的PSNR值較高，而在潮濕和濕潤的土壤環境反射下，K取1時PNMF的峰值信噪比較高。結合

表2 不同K 值下PNMF算法的PSNRTable 2 PSNR of PNMF algorithm with different K valuesdB

2.3 真實GPR數據

子空間技術中GPR圖像雜波去除的經驗[18]，選取第一主分量作為雜波成分，即K=1。

實驗的真實數據是從含有地下管道的表面進行采集的，地質松軟平坦，多含細沙碎石。使用GSSI公司提供的探地雷達設備，其中天線的中心頻率為400 MHz，探地深度為1.5 m，數據采集模式使用距離采集，在地面做好測距標識之后拖動裝置進行掃描，得到圖3 的A1，并依次在地下淺層處埋藏空箱和地雷模型得到圖3 的A2 和圖3 的A3。周圍有易產生電磁干擾的電線桿，會對采集的圖像造成干擾雜波。

如圖3的A1～A3所示。其雜波大體呈水平狀、數量多且與目標成分部分重疊的特點。實驗時選擇了信噪比和視覺質量這兩項評價指標來驗證PNMF 算法的有效性和魯棒性。

圖像的信噪比是用于比較處理后圖像與原圖像質量的參數，信噪比的數值越大，圖像的質量就越好。信噪比的計算往往需要借助于原始的純凈參考圖像，但由于實際測量的GPR圖像缺乏參考圖像，本實驗使用圖像的均值與方差之比來計算信噪比進行度量。真實的GPR數據見圖3中的A1～A3，分別使用NMF[18]、RNMF[19]和本文的PNMF 在GPR 圖像上進行處理，其SNR 結果以折線圖的形式呈現在圖4。可以看出在GPR 圖像的雜波處理中，RNMF 算法的信噪比高于NMF，而PNMF 算法得到的信噪比則比RNMF和NMF算法都要高。

圖3 不同算法對3幅真實GPR圖像的雜波抑制效果對比Fig.3 Comparison of clutter suppression effects of different algorithms on three real GPR images

圖4 NMF、RNMF和PNMF算法處理后的SNRFig.4 SNR value by NMF，RNMF and PNMF algorithms

為了進一步驗證本算法在GPR雜波處理中的魯棒性，分別在真實圖像中撒上密度為0.2、0.3和0.5的椒鹽噪聲，使用SVD、PCA、NMF、RNMF 和本文的PNMF 進行信噪比的計算，結果如表3所示。

表3 加入噪聲后不同算法的信噪比值Table 3 SNR of different algorithms after adding noisedB

從表3中的實驗數據可以看出，在實驗圖像中添加不同程度的噪聲時，本文算法的去噪效果均要優于其他經典的去雜波算法，信噪比SNR 值相比其他算法都有一定的提高，說明本文算法對不同的噪聲環境具有較強的魯棒性。

同時，滿足雜波抑制效果的前提下，算法的時間復雜度也是影響性能的關鍵因素。因此本文對比了基于NMF 和RNMF 算法運行所需的時間，運行時所使用的計算機操作系統為Windows 7，處理器為Inter Core i7-6700 3.40 GHz，8 GB 的安裝內存和64 位的操作系統。運行結果的具體數值見表4。

從表4 中的數據可以看出，對于圖3 的A1 和A2，PNMF的運行時間比NMF和RNMF都縮短了4倍左右，圖3 的A3 的運行時間與RNMF 相差不大，但總體的運行速度還是有所提高，這是因為PNMF中的變分貝葉斯屬于一種近似推理，能夠在較低的時間復雜度下獲得原問題的近似解。并且面對不同雜波類型時的抑制效果良好。

表4 算法的運行時間Table 4 Running time of algorithms

在圖像處理領域，視覺質量是最直觀也最重要的評價指標，效果圖如圖3 所示。針對3 種不同測量環境下的GPR 圖像（圖3（a）組）分別使用PCA（圖3（b）組）、SVD（圖3（c）組）傳統的雜波去除方法和已經在GPR圖像中應用過的NMF（圖3（d）組）、RNMF（圖3（e）組）方法進行實驗，并與本文的PNMF（圖3（f）組）進行比較。其中圖3（a）組是真實的GPR數據，可以看出3張探地雷達圖像的雜波形態不盡相同，但都屬于像素強度高于目標反射的高頻雜波。如圖3的A1中有一條高頻水平狀雜波在經過SVD、PCA、NMF和RNMF算法處理后水平雜波都得到了一定程度的削弱，但只有本文的PNMF方法在去除雜波的同時保留了區域1處的雙曲波，且區域2處的雜波處理也最干凈。背景相對純凈的同時很好地保留了目標信號的雙曲波。

圖3的A2的GPR圖像上方含有多條較細的高頻水平狀雜波，雜波之間相隔緊密，并覆蓋雙曲波頂部部分信息。在PCA（B2）中雜波的像素強度有所減弱但目標部分也變暗了；SVD（C2）中將第一條最強的雜波去除了，但后半段區域2處的雜波沒有本文的PNMF算法效果純凈；在NMF 算法（D2）中，在PCA 的效果基礎上對下面幾條雜波也進行了部分減弱；該圖像中RNMF（見圖3 的E2）算法表現最差，沒有起到雜波抑制的效果。PNMF對于該圖中的幾條雜波雖沒有完全抑制，但相比之下的背景是較為干凈的，且還對目標的邊緣細節進行了增強，見圖3的F2。

圖3的A3中含有3條較粗的高頻水平狀雜波（方框區域），雜波之間相隔較遠且與目標雙曲波部分存在大量重疊。可以看到，PCA、SVD、NMF、RNMF對這3條雜波的去除效果都很好，見B3、C3、D3 和E3，但是覆蓋在雙曲波上的雜波干擾情況并未得到改善（如區域1和區域2）。PNMF（圖F3）能準確高效地將雜波與目標分離，從而保留目標的雙曲波信息，使處理結果清晰自然。

以上實驗結果說明了在GPR圖像的高頻雜波抑制過程中，使用本文的PNMF 算法能獲得較高的信噪比，且視覺上使GPR 圖像背景清晰，突出了目標信息。從而驗證了本文算法的有效性和魯棒性。

3 結論

本文提出了一種改進非負矩陣分解的GPR圖像雜波抑制方法PNMF，與傳統的NMF不同的是，本文采用概率方法進行非負矩陣分解，得到的是分解后矩陣參數的后驗分布，然后選取矩陣分解后的第一主成分作為雜波進行去除。對于真實的復雜環境，概率模型往往比非概率方法具有更強的穩定性。PNMF 使用變分貝葉斯作為非負矩陣分解的近似推理方法，追求速度的同時獲得了較好的魯棒性，并且在使用過程中，結合雜波特點提取了探地雷達圖像的水平梯度作為輸入而不是像素強度。最后的實驗對比分析表明，在GPR 圖像的雜波抑制問題上，本文算法具有良好的有效性和魯棒性。后續工作將加入深度學習的知識，結合非負矩陣分解與深度卷積網絡對圖像數據進行深層次的特征提取，以此提升GPR圖像背景去除的純凈度和目標識別的準確率。