苗 翠，原 達(dá)，王冬雨，李文生

山東省高校智能信息處理重點(diǎn)實(shí)驗(yàn)室（山東工商學(xué)院），山東 煙臺(tái) 264005

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

對(duì)于探地雷達(dá)信號(hào)的雜波處理問(wèn)題，基于主成分分析（principal component analysis，PCA）[6]及魯棒主成分分析（robust principal component analysis，RPCA）[7]和奇異值分解（singular value decomposition，SVD）[8]等子空間技術(shù)往往表現(xiàn)出較好的性能，該類(lèi)方法將GPR圖像分解為對(duì)應(yīng)于雜波和目標(biāo)的分量，將主分量作為雜波，其余成分作為目標(biāo)信號(hào)來(lái)進(jìn)行處理，但該類(lèi)方法在GPR圖像含有多條雜波或多個(gè)目標(biāo)信號(hào)的情況下雜波不能完全去除。近兩年提出的形態(tài)學(xué)方法（morphological component analysis，MCA）[9]將GPR 圖像的雜波和目標(biāo)分量分別用曲波curvelet 和非采樣離散小波變換（unsampled discrete wavelet transform，UDWT）字典稀疏表示來(lái)進(jìn)行雜波抑制，雖然視覺(jué)效果表現(xiàn)不錯(cuò)，但算法復(fù)雜度較高，處理速度不夠快?；诙喾直媛实姆椒ㄔ贕PR 去雜波的處理中也取得了較好的效果，如多尺度雙邊濾波（multiscale directional bilateral filter，MDBF）[10]借助于雜波和目標(biāo)在反射波形中呈現(xiàn)出的幾何差異進(jìn)行多方向多尺度的分解。集成經(jīng)驗(yàn)?zāi)Ｊ椒纸猓╡nsemble empirical mode decomposition，EEMD）將探地雷達(dá)信號(hào)分解為一系列模函數(shù)（intrinsic mode functions，IMFs）并計(jì)算其置換熵（permutation entropy，PE），通過(guò)設(shè)置相關(guān)閾值來(lái)進(jìn)行噪聲和目標(biāo)的區(qū)分，能有效提高目標(biāo)的分辨率[11]，不過(guò)圖像質(zhì)量難以保證?；趶较蚧瘮?shù)（radial basis function，RBF）神經(jīng)網(wǎng)絡(luò)的方法也同樣被用于GPR 圖像的雜波濾除，使用零偏移格林函數(shù)作為訓(xùn)練數(shù)據(jù)的期望輸出，結(jié)果能提高探地雷達(dá)的垂直分辨率[12]，但對(duì)于目標(biāo)成分的保留不夠清晰。

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

1 基于概率模型的非負(fù)矩陣算法

1.1 構(gòu)建概率NMF模型

在NMF處理GPR圖像時(shí)，像素強(qiáng)度矩陣F看作是目標(biāo)成分Ftarget和雜波成分Fclutter的組合結(jié)構(gòu)：

將矩陣F進(jìn)行非負(fù)矩陣分解得到：

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

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

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

1.2 基于變分貝葉斯分解非負(fù)矩陣

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

噪聲方差的先驗(yàn)設(shè)為逆伽馬分布：

在PNMF 中，參數(shù)的后驗(yàn)分布使用變分貝葉斯（variational Bayes，VB）進(jìn)行推理，通過(guò)求解一個(gè)近似的后驗(yàn)分布去逼近真實(shí)的概率分布。在復(fù)雜的多變量參數(shù)θ情況下，VB將其分解為一組相互獨(dú)立的變量θi，表示如下：

其中，q(θ)對(duì)應(yīng)于參數(shù)θ的分布。根據(jù)文獻(xiàn)[20]，得到后驗(yàn)分布：

式（11）中，G(?)為伽馬分布，密度函數(shù)表達(dá)式見(jiàn)式（9），TN(?)為截?cái)喔咚狗植?，概率密度函?shù)表達(dá)式如下：

在截?cái)喔咚狗植贾?，x<0 時(shí)為密度為0 的高斯分布，其中φ(?)是服從N(0,1)分布的累積分布函數(shù)。得到W和H以及相關(guān)參數(shù)的后驗(yàn)分布后，將前K項(xiàng)作為高頻雜波對(duì)應(yīng)的分量并進(jìn)行如下分量表示：

其中K值代表矩陣的秩，通過(guò)該值可以得到雜波對(duì)應(yīng)的低秩矩陣。對(duì)于高頻雜波，前K項(xiàng)成分分量對(duì)雜波的貢獻(xiàn)最大，K值的選擇在實(shí)驗(yàn)部分進(jìn)行了具體論證。

2 實(shí)驗(yàn)

2.1 實(shí)驗(yàn)數(shù)據(jù)背景

本文實(shí)驗(yàn)分別采用GPR模擬數(shù)據(jù)和應(yīng)用環(huán)境采集的真實(shí)數(shù)據(jù)進(jìn)行方法驗(yàn)證分析，并通過(guò)不同算法的信噪比及視覺(jué)效果來(lái)驗(yàn)證本文算法的有效性及魯棒性。

2.2 模擬GPR數(shù)據(jù)

模擬數(shù)據(jù)通過(guò)gprMax軟件[21]的時(shí)域有限差分方法（FDTD）進(jìn)行模擬獲得。在模擬過(guò)程中，將鋁盤(pán)作為埋藏物，然后分別放置于3 種不同特性（干性土壤dry_sand、潮濕土壤damp_sand、濕潤(rùn)土壤wet_sand）的土壤環(huán)境下，埋藏物和3 種土壤的介質(zhì)參數(shù)來(lái)自于文獻(xiàn)[22]，其介電常數(shù)和導(dǎo)電率詳見(jiàn)表1。

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

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

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

峰值信噪比（PSNR）是圖像處理中最為普遍的一種客觀評(píng)價(jià)指標(biāo)，基于純凈參考圖像I與處理后圖像I?之間對(duì)應(yīng)像素點(diǎn)(i,j)的誤差來(lái)判斷圖像處理的質(zhì)量好壞。通常通過(guò)均方差（MSE）進(jìn)行定義：

其中M×N為圖像I的大小，在模擬GPR數(shù)據(jù)中，無(wú)雜波的參考圖像I通過(guò)含雜波的目標(biāo)反射圖像和只含雜波的圖像作差來(lái)得到。在使用概率模型進(jìn)行非負(fù)矩陣分解時(shí)，分別使用圖像像素強(qiáng)度和水平梯度作為輸入，通過(guò)PSNR 值的計(jì)算得到一組對(duì)比柱狀圖，如圖2 所示?？梢灾庇^看到在使用梯度進(jìn)行PNMF 探地雷達(dá)圖像雜波處理時(shí)，得到的PSNR值更高。

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

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

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

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

2.3 真實(shí)GPR數(shù)據(jù)

子空間技術(shù)中GPR圖像雜波去除的經(jīng)驗(yàn)[18]，選取第一主分量作為雜波成分，即K=1。

實(shí)驗(yàn)的真實(shí)數(shù)據(jù)是從含有地下管道的表面進(jìn)行采集的，地質(zhì)松軟平坦，多含細(xì)沙碎石。使用GSSI公司提供的探地雷達(dá)設(shè)備，其中天線的中心頻率為400 MHz，探地深度為1.5 m，數(shù)據(jù)采集模式使用距離采集，在地面做好測(cè)距標(biāo)識(shí)之后拖動(dòng)裝置進(jìn)行掃描，得到圖3 的A1，并依次在地下淺層處埋藏空箱和地雷模型得到圖3 的A2 和圖3 的A3。周?chē)幸桩a(chǎn)生電磁干擾的電線桿，會(huì)對(duì)采集的圖像造成干擾雜波。

如圖3的A1～A3所示。其雜波大體呈水平狀、數(shù)量多且與目標(biāo)成分部分重疊的特點(diǎn)。實(shí)驗(yàn)時(shí)選擇了信噪比和視覺(jué)質(zhì)量這兩項(xiàng)評(píng)價(jià)指標(biāo)來(lái)驗(yàn)證PNMF 算法的有效性和魯棒性。

圖像的信噪比是用于比較處理后圖像與原圖像質(zhì)量的參數(shù)，信噪比的數(shù)值越大，圖像的質(zhì)量就越好。信噪比的計(jì)算往往需要借助于原始的純凈參考圖像，但由于實(shí)際測(cè)量的GPR圖像缺乏參考圖像，本實(shí)驗(yàn)使用圖像的均值與方差之比來(lái)計(jì)算信噪比進(jìn)行度量。真實(shí)的GPR數(shù)據(jù)見(jiàn)圖3中的A1～A3，分別使用NMF[18]、RNMF[19]和本文的PNMF 在GPR 圖像上進(jìn)行處理，其SNR 結(jié)果以折線圖的形式呈現(xiàn)在圖4?？梢钥闯鲈贕PR 圖像的雜波處理中，RNMF 算法的信噪比高于NMF，而PNMF 算法得到的信噪比則比RNMF和NMF算法都要高。

圖3 不同算法對(duì)3幅真實(shí)GPR圖像的雜波抑制效果對(duì)比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

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

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

從表3中的實(shí)驗(yàn)數(shù)據(jù)可以看出，在實(shí)驗(yàn)圖像中添加不同程度的噪聲時(shí)，本文算法的去噪效果均要優(yōu)于其他經(jīng)典的去雜波算法，信噪比SNR 值相比其他算法都有一定的提高，說(shuō)明本文算法對(duì)不同的噪聲環(huán)境具有較強(qiáng)的魯棒性。

同時(shí)，滿足雜波抑制效果的前提下，算法的時(shí)間復(fù)雜度也是影響性能的關(guān)鍵因素。因此本文對(duì)比了基于NMF 和RNMF 算法運(yùn)行所需的時(shí)間，運(yùn)行時(shí)所使用的計(jì)算機(jī)操作系統(tǒng)為Windows 7，處理器為Inter Core i7-6700 3.40 GHz，8 GB 的安裝內(nèi)存和64 位的操作系統(tǒng)。運(yùn)行結(jié)果的具體數(shù)值見(jiàn)表4。

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

表4 算法的運(yùn)行時(shí)間Table 4 Running time of algorithms

在圖像處理領(lǐng)域，視覺(jué)質(zhì)量是最直觀也最重要的評(píng)價(jià)指標(biāo)，效果圖如圖3 所示。針對(duì)3 種不同測(cè)量環(huán)境下的GPR 圖像（圖3（a）組）分別使用PCA（圖3（b）組）、SVD（圖3（c）組）傳統(tǒng)的雜波去除方法和已經(jīng)在GPR圖像中應(yīng)用過(guò)的NMF（圖3（d）組）、RNMF（圖3（e）組）方法進(jìn)行實(shí)驗(yàn)，并與本文的PNMF（圖3（f）組）進(jìn)行比較。其中圖3（a）組是真實(shí)的GPR數(shù)據(jù)，可以看出3張?zhí)降乩走_(dá)圖像的雜波形態(tài)不盡相同，但都屬于像素強(qiáng)度高于目標(biāo)反射的高頻雜波。如圖3的A1中有一條高頻水平狀雜波在經(jīng)過(guò)SVD、PCA、NMF和RNMF算法處理后水平雜波都得到了一定程度的削弱，但只有本文的PNMF方法在去除雜波的同時(shí)保留了區(qū)域1處的雙曲波，且區(qū)域2處的雜波處理也最干凈。背景相對(duì)純凈的同時(shí)很好地保留了目標(biāo)信號(hào)的雙曲波。

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

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

以上實(shí)驗(yàn)結(jié)果說(shuō)明了在GPR圖像的高頻雜波抑制過(guò)程中，使用本文的PNMF 算法能獲得較高的信噪比，且視覺(jué)上使GPR 圖像背景清晰，突出了目標(biāo)信息。從而驗(yàn)證了本文算法的有效性和魯棒性。

3 結(jié)論

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