陳 誠(chéng), 劉 濤, 宋曉驥, 粟 毅, 金光虎

(國(guó)防科技大學(xué)電子科學(xué)學(xué)院, 湖南 長(zhǎng)沙 410073)

0 引 言

表層穿透雷達(dá)(surface penetrating radar,SPR)可實(shí)現(xiàn)對(duì)各種非可視隱蔽目標(biāo)直觀的定位和探測(cè)[1-2],并憑借其分辨率高、高效、便捷等優(yōu)點(diǎn),受到了研究者的密切關(guān)注,在建筑工程、安全檢查以及排雷防爆等領(lǐng)域均有應(yīng)用[3-4]。

SPR通常工作于反射模式下,由發(fā)射天線向被測(cè)介質(zhì)和目標(biāo)發(fā)射電磁波,在介質(zhì)材料內(nèi)的不連續(xù)處(如具有不同介電性能的兩種材料之間的邊界或介質(zhì)內(nèi)部出現(xiàn)的斷裂、損傷等),部分能量被反射回SPR,而一部分通過介質(zhì)繼續(xù)傳播。反射信號(hào)由接收天線接收并做處理分析,探測(cè)回波的幅度及相位信息可用來定位介質(zhì)中存在的不連續(xù)點(diǎn)。起初,SPR僅運(yùn)用于探測(cè)介質(zhì)中目標(biāo)有無及定位。隨著成像技術(shù)的發(fā)展,研究者將回波數(shù)據(jù)轉(zhuǎn)換為直觀的目標(biāo)圖像,更加符合人們的視覺感官,更為重要的是,通過成像處理,提高了方位向分辨率以及目標(biāo)的更多幾何信息,成像處理也因此成為SPR的一個(gè)發(fā)展重點(diǎn)。在SPR發(fā)展過程中,不斷有學(xué)者提出各類型的成像算法,如反向投影(back projection,BP)算法[5]、衍射層析(diffraction tomography, DT)算法[6]、距離偏移(range migration, RM)算法[7]等。利用SPR成像技術(shù)進(jìn)行無損探測(cè),還存在著一些困難和挑戰(zhàn),傳統(tǒng)成像算法依賴于準(zhǔn)確的介質(zhì)波速等參數(shù),而參數(shù)設(shè)置的誤差可能導(dǎo)致成像結(jié)果發(fā)生散焦而無法辨認(rèn)目標(biāo)。

對(duì)于SPR系統(tǒng),提高成像分辨率可從硬件上優(yōu)化系統(tǒng)設(shè)計(jì),提高發(fā)射信號(hào)帶寬或增加合成孔徑長(zhǎng)度,但過大的帶寬又會(huì)增加系統(tǒng)的色散,因此人們更趨向于在既定的硬件基礎(chǔ)上,通過算法處理提高成像分辨率。目前,帶寬外推(bandwidth extrapolation,BWE)方法發(fā)展較為成熟,實(shí)現(xiàn)簡(jiǎn)單,但是對(duì)數(shù)據(jù)和模型匹配程度的要求較高[8],并且當(dāng)場(chǎng)景中噪聲干擾較為明顯時(shí),BWE 方法成像結(jié)果難以令人滿意。譜估計(jì)(spectral estimation, SE)方法能夠替換傳統(tǒng)成像算法中的傅里葉變換功能,利用譜估計(jì)實(shí)現(xiàn)高分辨率成像能夠有效抑制旁瓣,但通常計(jì)算過程較為復(fù)雜[9-10]。稀疏信號(hào)處理近年來被廣泛應(yīng)用于雷達(dá)成像[11-13],特別是壓縮感知(compressive sensing,CS)方法[14-15],能夠突破奈奎斯特采樣定理限制,利用信號(hào)的稀疏性、測(cè)量矩陣的隨機(jī)性以及非線性優(yōu)化算法來實(shí)現(xiàn)對(duì)信號(hào)的壓縮采樣和重構(gòu),在雷達(dá)成像質(zhì)量改善方面效果顯著[16-17]。典型的稀疏特征增強(qiáng)方法包括l1范數(shù)正則化(L1-based regularization,L1R)[18]、K-Support[19]和Elastic Net[20]等。此外,馮俊杰等人提出了一種基于多測(cè)量向量塊稀疏信號(hào)重構(gòu)的成像算法[21],通過構(gòu)建稀疏模型并對(duì)稀疏信號(hào)優(yōu)化求解,提高了成像質(zhì)量。盡管上述方法均在一定程度上改善了成像結(jié)果,但復(fù)雜的建模及求解過程勢(shì)必增加成像處理時(shí)間,因此在保持運(yùn)算速度的條件下提高成像質(zhì)量對(duì)于SPR無損探測(cè)應(yīng)用意義重大。

鑒于上述問題,本文提出一種基于稀疏信號(hào)處理的高效穿透成像增強(qiáng)算法,該方法在基于L1R稀疏成像方法的基礎(chǔ)上[22],結(jié)合雷達(dá)回波特性,利用交替方向乘子法(alternating direction method of multipliers, ADMM)[23],建立帶有增廣拉格朗日項(xiàng)的優(yōu)化方程,求解目標(biāo)函數(shù)并通過對(duì)偶迭代運(yùn)算得到稀疏特征增強(qiáng)后的復(fù)圖像數(shù)據(jù),以提高成像精度和效率。本文實(shí)驗(yàn)應(yīng)用仿真與實(shí)測(cè)數(shù)據(jù),通過與多種方法進(jìn)行對(duì)比,定量和定性的分析了本文方法的性能優(yōu)勢(shì),表明所提方法能夠在提升成像處理速度的同時(shí),有效提高成像結(jié)果信雜比。

1 基于L1R的稀疏成像方法

基于壓縮感知理論,將稀疏信號(hào)處理運(yùn)用于雷達(dá)成像中,可在低于奈奎斯特采樣率的條件下通過求解l1范數(shù)約束下的凸優(yōu)化問題實(shí)現(xiàn)對(duì)信號(hào)的精確重構(gòu),能夠有效抑制旁瓣和噪聲,從而提高成像質(zhì)量[11-22]。對(duì)于SPR系統(tǒng),令ΩR表示原始成像數(shù)據(jù)矩陣,ΩS表示基于稀疏成像方法重構(gòu)的圖像數(shù)據(jù)矩陣,兩者之間的差值N為雜波、噪聲等,即

N=ΩR-ΩS

(1)

基于稀疏信號(hào)處理理論,可通過基于L1R的方法對(duì)稀疏圖像進(jìn)行重構(gòu):

(2)

通過迭代方法求解式(2),得到迭代方程為

(3)

(4)

通過閾值收縮方法[24],即設(shè)定某個(gè)閾值,當(dāng)函數(shù)值大于某個(gè)閾值時(shí)對(duì)函數(shù)進(jìn)行線性縮小,而當(dāng)函數(shù)值小于閾值時(shí)歸零處理,由此估算出第i+1步迭代過程中的場(chǎng)景稀疏解ΔΩ(i+1):

(5)

式中：Sμη(·)為實(shí)數(shù)軟閾值算子；sgn(·)表示符號(hào)函數(shù)；η(i)為正則化參數(shù),影響著算法成像精度及場(chǎng)景稀疏度,其定義為

(6)

(7)

值得注意的是,迭代參數(shù)μ為0到1之間的一個(gè)小數(shù),控制著迭代算法的收斂速度。當(dāng)μ由0逐漸增大時(shí),算法的收斂速度變快,但同時(shí)稀疏重構(gòu)結(jié)果的精度會(huì)隨之降低。因此,在表層穿透雷達(dá)實(shí)際應(yīng)用中,應(yīng)綜合算法精度和計(jì)算速度兩者以選擇合適的μ值。

2 基于稀疏信號(hào)處理的穿透成像增強(qiáng)算法

根據(jù)ADMM算法思想,將基于L1R的稀疏成像方法分解為兩個(gè)子問題,即

則式(2)中的重構(gòu)問題轉(zhuǎn)換為

s.t.X=Y

(8)

進(jìn)一步,結(jié)合增廣的拉格朗日多乘子方法建立優(yōu)化方程,得到

(9)

式中:U表示對(duì)偶變量矩陣;(·)H表示對(duì)矩陣進(jìn)行共軛轉(zhuǎn)置;λ為拉格朗日乘子系數(shù),將約束條件通過對(duì)偶變量與原本的最小化問題耦合成拉格朗日函數(shù),并引入拉格朗日乘子項(xiàng),將目標(biāo)函數(shù)的解約束于可行域內(nèi),使得算法更加穩(wěn)健。

式(9)求解的過程可以分解為3個(gè)子問題,即分別固定3個(gè)變量中的兩個(gè)變量,更新其中一個(gè)變量,具體步驟如下：

步驟 1固定變量Y和U,估計(jì)X,即

步驟 2固定變量X和U,估計(jì)Y,即

步驟 3固定變量X和Y,估計(jì)U,即

Uk+1=Uk+λ(Xk+1-Yk+1)

其中,上標(biāo)k表示圖像矩陣第k次的迭代值。

步驟1中的問題可以理解為凸優(yōu)化問題,通過對(duì)回歸系數(shù)增加懲罰項(xiàng)來解決最小二乘法的思想,可以求解出

Xk+1={(λ+1)I}-1{ΩR+λ(Yk-Uk)}

(10)

求解步驟2中變量Y時(shí),同樣采用類似于式(5)中的閾值收縮方法,值得注意的是,通常雷達(dá)回波為復(fù)數(shù)據(jù),如果采用實(shí)數(shù)軟閾值算子會(huì)損失雷達(dá)回波的復(fù)數(shù)特性,因此本文在求解變量Y時(shí),通過將變量函數(shù)與其行列式的比值替代符號(hào)函數(shù)得到復(fù)數(shù)軟閾值算子,即

Yk+1=(Xk+1+Uk)/

|Xk+1+Uk|max{|Xk+1+Uk|,μ/λ}

(11)

然后,根據(jù)步驟3中的求解公式,對(duì)變量U進(jìn)行迭代,并計(jì)算出殘差:

Rk=‖Xk-Zk‖F(xiàn)

(12)

3 實(shí)驗(yàn)結(jié)果和分析

這里通過仿真實(shí)驗(yàn)和實(shí)測(cè)實(shí)驗(yàn)兩種方式對(duì)算法性能進(jìn)行驗(yàn)證,并通過與L1R、K-Support和Elastic Net 這3種稀疏方法處理后成像結(jié)果進(jìn)行對(duì)比,證明算法在成像質(zhì)量和速度上的優(yōu)勢(shì)。

3.1 仿真實(shí)驗(yàn)設(shè)計(jì)

仿真通過CST電磁仿真軟件獲取場(chǎng)景數(shù)據(jù),實(shí)驗(yàn)介質(zhì)為表面略有起伏(不超過10 mm)的干燥壤土,相對(duì)介電常數(shù)為2.44,壤土內(nèi)部埋藏r=35 mm,高h(yuǎn)=45 mm的塑料圓柱,相對(duì)介電常數(shù)為2.8,距離壤土表面d=20 mm,用于模擬埋藏在土壤中的地雷。發(fā)射信號(hào)頻率為5 GHz。天線掃描面距離介質(zhì)表面50 mm,掃描面大小為200 mm×200 mm,掃描步進(jìn)間隔為10 mm。仿真實(shí)驗(yàn)場(chǎng)景如圖1所示。

圖1 仿真實(shí)驗(yàn)場(chǎng)景圖

3.2 實(shí)測(cè)實(shí)驗(yàn)設(shè)計(jì)

實(shí)驗(yàn)1使用的系統(tǒng)包括一臺(tái)矢量網(wǎng)絡(luò)分析儀(安捷倫 N5230C)及一套可以精確定位的二維移動(dòng)掃描平臺(tái),接收天線與發(fā)射天線安裝在平臺(tái)上,通過計(jì)算機(jī)控制平臺(tái)的移動(dòng)與矢量網(wǎng)絡(luò)分析儀數(shù)據(jù)的接收和處理,實(shí)驗(yàn)場(chǎng)景如圖2所示。

圖2 實(shí)驗(yàn)1場(chǎng)景圖

實(shí)驗(yàn)1采用的目標(biāo)為一塊厚度為44 mm的混凝土試件,介電常數(shù)為6.5,在其背面布設(shè)兩根平行的金屬棒用于模擬混凝土墻內(nèi)的鋼筋,兩根金屬棒之間間隔為5 cm,其中一根金屬棒內(nèi)部存在1 cm的空隙用于模擬鋼筋斷裂,介質(zhì)與目標(biāo)的詳細(xì)位置關(guān)系如圖3所示。

圖3 實(shí)驗(yàn)1目標(biāo)與介質(zhì)位置關(guān)系

掃描頻率范圍設(shè)置為4.5～6.5 GHz,頻段內(nèi)共計(jì)201個(gè)頻點(diǎn),在頻率范圍內(nèi)均勻分布。天線距離墻體0.01 m,掃描區(qū)域范圍為0.25 m×0.25 m,采樣間隔為5 mm。將系統(tǒng)接收到的回波數(shù)據(jù)進(jìn)行成像,在此基礎(chǔ)上利用本文方法與其他幾種稀疏算法進(jìn)行處理,得到特征增強(qiáng)結(jié)果。

實(shí)驗(yàn)2為埋設(shè)于龍骨墻內(nèi)的“NUDT”金屬字樣,龍骨墻由石膏墻面和膠合板組成(整體厚度約3 cm),如圖4所示。使用本實(shí)驗(yàn)室研發(fā)的SPR系統(tǒng)對(duì)墻體內(nèi)部進(jìn)行掃描獲取回波數(shù)據(jù),并按實(shí)驗(yàn)1同樣處理流程進(jìn)行處理。

圖4 實(shí)驗(yàn)2場(chǎng)景圖

3.3 實(shí)驗(yàn)結(jié)果與分析

圖5(a)和圖5(b)分別給出了對(duì)仿真場(chǎng)景使用算法處理前后的成像結(jié)果。

圖5 仿真實(shí)驗(yàn)結(jié)果圖

從仿真結(jié)果可以看出,由于壤土表層的非均勻散射等原因,未進(jìn)行特征增強(qiáng)時(shí),成像結(jié)果中雜波和噪聲明顯,影響圖像質(zhì)量。通過使用本文算法進(jìn)行處理后,能夠在保留目標(biāo)完整性的前提下,對(duì)成像結(jié)果中雜波和噪聲進(jìn)行抑制。兩個(gè)實(shí)測(cè)實(shí)驗(yàn)場(chǎng)景的實(shí)驗(yàn)結(jié)果分別如圖6和圖7所示。

圖6 實(shí)驗(yàn)1結(jié)果圖

圖7 實(shí)驗(yàn)2結(jié)果圖

實(shí)驗(yàn)1和實(shí)驗(yàn)2的結(jié)果表明,幾種算法均能對(duì)成像結(jié)果實(shí)現(xiàn)特征增強(qiáng)并有效抑制旁瓣和噪聲,提高圖像分辨率,但是在改善效果及處理時(shí)間上有所差異。從成像效果改善效果來看,本文算法效果最好,在實(shí)驗(yàn)1成像結(jié)果中,能夠清晰分辨金屬棒的斷裂處,K-support方法次之,而L1R和Elastic Net方法對(duì)雜波和噪聲的抑制效果表現(xiàn)相對(duì)較差,特別是在實(shí)驗(yàn)1成像結(jié)果中,還殘留著較為明顯的干涉條紋,且難以判斷金屬棒是否斷裂;表1給出了實(shí)驗(yàn)1和實(shí)驗(yàn)2中數(shù)據(jù)通過各種方法處理的時(shí)間,從處理時(shí)間上來看,本文算法處理數(shù)據(jù)速度略優(yōu)于L1R和Elastic Net方法(約3～8倍),明顯優(yōu)于K-support方法(約7～24倍),這是因?yàn)長(zhǎng)1R方法通過最小化l1范數(shù)來實(shí)現(xiàn)稀疏化,能夠避免求解l0范數(shù)NP難問題或是貪婪算法求解效率低的問題,但是L1R方法在追求稀疏化而犧牲了對(duì)數(shù)據(jù)局部結(jié)構(gòu)特征的保留能力,不能完全保證特征增強(qiáng)后圖像接近稀疏表示ΩS,Elastic Net同時(shí)引入l1范數(shù)和l2范數(shù),通過正則化參數(shù)來調(diào)節(jié)l1范數(shù)和l2范數(shù)的權(quán)重,以平衡參數(shù)模型稀疏與光滑性來使求解更加準(zhǔn)確,但是正則化參數(shù)的選擇會(huì)對(duì)成像結(jié)果產(chǎn)生影響,K-Support方法將矩陣稀疏性的凸松弛與l2范數(shù)相結(jié)合,能夠在實(shí)現(xiàn)模型稀疏性的同時(shí)保證求解精度,但計(jì)算效率低于其他幾種方法。這幾種方法均局限在實(shí)數(shù)域處理圖像,而本文算法在基于L1R的稀疏成像方法的基礎(chǔ)上,結(jié)合雷達(dá)回波復(fù)數(shù)特性,利用ADMM算法,建立帶有增廣拉格朗日項(xiàng)的優(yōu)化方程,能夠保證得到的稀疏解更為精確,且在提高成像效果的同時(shí)通過對(duì)偶迭代運(yùn)算保證了成像速度。

表1 不同方法處理時(shí)間對(duì)比

為了定量分析圖像質(zhì)量的改善程度,本文采用目標(biāo)信雜比(signal to clutter ratio,SCR)來對(duì)成像效果進(jìn)行評(píng)估,即圖像目標(biāo)區(qū)域平均功率與雜波區(qū)域平均功率的比值,其定義為

(13)

式中:Pt表示目標(biāo)區(qū)域平均功率;Pc表示雜波區(qū)域平均功率；Rt和Rc分別表示目標(biāo)和雜波區(qū)域；Nt和Nc分別表示目標(biāo)和雜波區(qū)域內(nèi)像素點(diǎn)個(gè)數(shù)；I(p)和I(q)分別表示目標(biāo)和雜波區(qū)域各點(diǎn)的像素值大小。

表2給出了不同方法處理后成像結(jié)果的SCR。通過對(duì)比4種方法在各個(gè)實(shí)驗(yàn)場(chǎng)景的SCR,可以看出本文提出的稀疏特征增強(qiáng)方法在改善SCR方面均優(yōu)于其他幾種方法,K-Support方法次之,Elastic Net方法和L1R方法效果改善效果相對(duì)較弱,表明本文方法能夠有效地增強(qiáng)圖像質(zhì)量,抑制雜波和噪聲。

表2 不同方法處理后成像結(jié)果SCR參數(shù)

4 結(jié) 論

本文針對(duì)表層穿透雷達(dá)應(yīng)用于無損檢測(cè)成像時(shí)難以兼顧精度與效率的問題,提出一種基于稀疏信號(hào)處理的穿透成像增強(qiáng)算法。該方法在基于L1R的稀疏成像方法的基礎(chǔ)上,結(jié)合交替方向乘子法分解優(yōu)化問題通過建立帶有增廣拉格朗日項(xiàng)的優(yōu)化方程,在優(yōu)化最小化框架下求解最優(yōu)化函數(shù),最后根據(jù)目標(biāo)函數(shù)通過對(duì)偶迭代運(yùn)算得到稀疏特征增強(qiáng)后的復(fù)圖像數(shù)據(jù),在保持運(yùn)算速度的同時(shí)實(shí)現(xiàn)SPR高質(zhì)量成像。但本文方法在圖像中存在強(qiáng)雜波時(shí)對(duì)成像效果改善有限,因此在使用本文方法時(shí),通常需要先對(duì)原始數(shù)據(jù)進(jìn)行雜波抑制處理。