綜合孔徑輻射計偏微分方程近場圖像反演算法

姚現(xiàn)勛，尚曉舟，苗俊剛，李志平

(北京航空航天大學(xué) 電子信息工程學(xué)院，北京100191)

微波輻射計是一種微波無源遙感器，通過接收目標(biāo)的微波輻射信號來獲取目標(biāo)的亮溫信息，具有很好的安全性和隱蔽性，在反恐探測［1］、人體安檢［2］等領(lǐng)域具有廣泛的應(yīng)用前景.相比于X-ray，微波成像系統(tǒng)不僅可以檢測出隱藏在織物下的金屬物品，還可以檢測出陶瓷刀具、炸藥等危險品，獲得更加詳盡、準(zhǔn)確的信息.為了實時地獲取高分辨率微波圖像，綜合孔徑成像技術(shù)成為了人體安檢領(lǐng)域的研究熱點［3-4］.

綜合孔徑輻射計成像是將稀疏分布的小孔徑天線之間的干涉測量通過數(shù)字波束合成的辦法將其綜合成一個大的等效孔徑，并通過干涉測量獲得的可視度函數(shù)反演得到視場范圍內(nèi)的亮溫分布.理想情況下，在天線遠(yuǎn)場區(qū)，可視度函數(shù)與目標(biāo)亮溫分布之間滿足傅里葉變換關(guān)系.然而，對于人體安檢等領(lǐng)域的應(yīng)用，探測目標(biāo)處于天線陣的近場區(qū)域，此時可視度函數(shù)與目標(biāo)的亮溫分布不再滿足傅里葉變換關(guān)系［5］.因此，近場成像方法成為了綜合孔徑輻射計應(yīng)用于近距離探測領(lǐng)域的一個關(guān)鍵技術(shù).文獻(xiàn)［6］采用修正相位的點聚焦算法以重建相位修正后的可視度函數(shù)與目標(biāo)亮溫分布的傅里葉變換關(guān)系.然而，此算法對大視場范圍的擴(kuò)展目標(biāo)無法進(jìn)行準(zhǔn)確的近場修正，成像誤差較大.對于一維成像系統(tǒng)，文獻(xiàn)［7］借助數(shù)值計算的方法，利用天線陣列形式及目標(biāo)分布建立近場可視度函數(shù)矩陣，通過數(shù)值求逆的方式得到目標(biāo)的亮溫分布，消除近場誤差.然而，對于二維成像系統(tǒng)，隨著通道數(shù)目的增加，矩陣方程的穩(wěn)定性急劇惡化，從而使得圖像反演成為一種病態(tài)問題［8］.基于偏微分方程的正則化方法是求解病態(tài)問題的一種有效方法［9］.根據(jù)尺度空間理論，基于L2范數(shù)的各向同性擴(kuò)散模型有利于去噪，但在保留圖像邊緣和細(xì)節(jié)方面存在不足;基于L1范數(shù)的各向異性擴(kuò)散可以保持圖像邊緣細(xì)節(jié)，但會引入階梯效應(yīng)［10］.

本文將局部自適應(yīng)偏微分方程的思想引入到綜合孔徑輻射計近場成像領(lǐng)域.利用目標(biāo)分布的先驗信息，將整個視場范圍分為背景區(qū)域和目標(biāo)區(qū)域.在圖像背景區(qū)域利用L2范數(shù)進(jìn)行噪聲抑制，在目標(biāo)區(qū)域利用自適應(yīng)全變分模型來保留圖像邊緣.在濾除測量噪聲的同時更好地保留圖像邊緣細(xì)節(jié)信息.為了驗證該反演算法的有效性，采用一套48通道的二維綜合孔徑輻射計BHU-2DU進(jìn)行了仿真、實驗驗證.結(jié)果表明，此算法比傳統(tǒng)的偏微分方程方法更有效地減小了可視度測量噪聲對亮溫圖像反演的影響.

1 綜合孔徑輻射近場成像模型

干涉式綜合孔徑輻射計測量的是不同天線單元接收信號的復(fù)相關(guān)值，即可視度函數(shù).二元干涉儀是構(gòu)成綜合孔徑輻射計的基本單元，其干涉測量的幾何關(guān)系如圖1所示.天線陣列位于z=0平面上，輻射面源位于z=h平面上.

圖1 近場干涉測量幾何關(guān)系圖Fig.1 Near field geometric diagram of interferometry

將面輻射源離散成P個小面源.理想情況下，對任意兩天線單元的接收信號進(jìn)行復(fù)相關(guān)運算，測得目標(biāo)的可視度函數(shù)為［5］

其中，T(θkj，i，φi)為第i個小面源的輻射亮溫;rk，i和rj，i分別為該小面源到兩接收天線之間的距離;Ω為單元天線的固體角;Fn為歸一化的天線電壓方向圖;為考慮了空間去相關(guān)效應(yīng)的條紋洗滌函數(shù);對于理想窄帶系統(tǒng)，去相關(guān)效應(yīng)≈1，可以忽略.綜合孔徑輻射計的遠(yuǎn)場條件通常定義［11］:

其中，D為天線陣的孔徑尺寸;λ為系統(tǒng)工作頻率所對應(yīng)的波長.此時，式(1)描述的可視度函數(shù)變?yōu)槟繕?biāo)修正亮溫的傅里葉變換，直接對測得的可視度函數(shù)進(jìn)行逆傅里葉變換就可獲得目標(biāo)的亮溫分布.對于近場情況，測量的可視度函數(shù)不僅與天線陣列的分布有關(guān)，且與目標(biāo)場景的空間位置有關(guān).因此，綜合孔徑近場成像的復(fù)雜性使得其不存在精確的解析解.考慮到實際系統(tǒng)中只能測得可視度函數(shù)離散采樣點的值，式(1)可轉(zhuǎn)化成離散矩陣方程的形式:

其中，G為系統(tǒng)響應(yīng)矩陣;M為可視度函數(shù)采樣點個數(shù).理想情況下，可直接通過矩陣求逆獲得精確的近場亮溫圖像.然而，為了保證足夠的離散精度，通常情況下有P≥3M［8］.因此，此方程組為欠定方程組，解不唯一.同時考慮到實際系統(tǒng)不可避免的測量噪聲，其近場可視度模型式(2)轉(zhuǎn)化為

其中n為可視度采樣測量噪聲.上述矩陣方程的穩(wěn)定性主要是由系統(tǒng)響應(yīng)G矩陣的條件數(shù)決定，并隨著去相關(guān)效應(yīng)的加劇、接收通道數(shù)目的增加而增大.實際中為了達(dá)到所需的視場范圍和空間分辨率，成像系統(tǒng)通常需要較多的接收通道，從而使得綜合孔徑輻射計近場圖像反演成為一個病態(tài)求逆問題，即微小的測量噪聲也可能導(dǎo)致反演圖像的完全失真.正則化是將病態(tài)問題良性化的一個有效方法，其基本思想是利用解的先驗知識，構(gòu)造附加約束或者改變求解策略，使得逆問題的解變得穩(wěn)定和確定.

2 基于偏微分方程的成像算法

2.1 偏微分方程正則化方法

圖像反演中，正則化方法是通過極小化下面的約束誤差方程來求得目標(biāo)的亮溫分布［9］:

其中，第1項表示估計值與實際數(shù)據(jù)的誤差;第2項為解的邊界約束;Ω為圖像的支持域;μ＞0為正則化參數(shù).對于基于L2范數(shù)的偏微方程成像模型，邊界約束為

代入式(5)，可得其Euler-Lagrange方程為

從上式可看出，擴(kuò)散項 ·(T)的擴(kuò)散系數(shù)為1，即在各個方向上具有相同的擴(kuò)散能力.因此，該方法能夠濾除噪聲，但會損失圖像的細(xì)節(jié)和邊緣信息.

為了保持圖像邊緣，基于L1范數(shù)的全變分模型將式(5)中的邊界約束改寫為［10］

其對應(yīng)的Euler-Lagrange方程為

2.2 局部自適應(yīng)偏微分方程方法

為了實現(xiàn)較高的空間分辨率，需要較大的綜合孔徑尺寸，在這種情況下，為了減少接收機(jī)通道數(shù)量降低運算復(fù)雜度，系統(tǒng)的天線間隔需要增加.然而，較大的天線間隔將造成完全混疊的視場范圍，因此，需要采用背景對消［11-12］的辦法獲得一定的視場范圍.雖然這種方法不能獲得目標(biāo)的絕對亮溫，但是，在安檢等應(yīng)用中，相對亮溫圖像已經(jīng)能夠滿足應(yīng)用需求.在背景對消情況下，系統(tǒng)獲得的差分可視度可表示為

其中，Vt+b0表示目標(biāo)場景Tt+b0的可視度函數(shù);Vb+b0表示背景場景Tb+b0的可視度函數(shù);下標(biāo)t表示視場范圍中存在目標(biāo)時的目標(biāo)區(qū)域;下標(biāo)b表示視場范圍中不存在目標(biāo)時的目標(biāo)區(qū)域;下標(biāo)b0表示視場范圍中除目標(biāo)區(qū)域外的背景區(qū)域;T'表示待求解的差分亮溫圖像，為目標(biāo)亮溫與背景亮溫之差.

經(jīng)背景對消處理后，待求解的差分亮溫圖像一般可以分為兩個特征區(qū)域:背景區(qū)和目標(biāo)區(qū)，且具有明顯的分界線.根據(jù)此先驗信息，提出一種局部自適應(yīng)偏微分方程的綜合孔徑輻射計近場成像方法，在不同的圖像區(qū)域采用不同的約束，其模型為

其中

其中，Ωb為圖像的背景區(qū)域;Ωt為圖像的目標(biāo)區(qū)域;Gσ為高斯濾波器，用來對噪聲圖像進(jìn)行預(yù)處理.式(11)對應(yīng)的Euler-Lagrange方程為

上式表明，該模型在圖像背景區(qū)域采用基于L2范數(shù)的偏微分方程模型，以便盡可能地抑制噪聲;在目標(biāo)區(qū)域采用自適應(yīng)全變分模型，邊緣處較大，p≈1保持圖像邊緣;平坦區(qū)較小，p≈2抑制噪聲.根據(jù)目標(biāo)區(qū)域的圖像特征自適應(yīng)的選擇擴(kuò)散方向，從而在保持邊緣的同時有效地抑制了階梯效應(yīng).實際中，首先采用點聚焦FFT［6］獲得質(zhì)量較差的場景亮溫圖像，從中提取出目標(biāo)的輪廓，再采用本文的方法對目標(biāo)圖像進(jìn)行精確反演.對方程(13)，可采用兩步迭代方法(TwIST)進(jìn)行快速求解［13］.

3 近場成像仿真

為了分析上述圖像反演算法的正確性及可行性，根據(jù)現(xiàn)有系統(tǒng)BHU-2D-U的技術(shù)參數(shù)進(jìn)行了成像仿真，其主要參數(shù)如表1所示.

表1 BHU-2D-U系統(tǒng)主要參數(shù)Table1 Main specifications of BHU-2D-U

根據(jù)上述參數(shù)，為滿足遠(yuǎn)場成像條件，BHU-2D-U的天線陣與探測目標(biāo)的距離需大于39.1m.然而，為了獲得1 m×2 m的視場范圍，其成像距離限制約為3 m，屬于近場區(qū)域，無法直接采用傅里葉變換進(jìn)行成像.仿真的原始場景如圖2(a)所示，人體亮溫設(shè)為280 K，周圍環(huán)境設(shè)為0 K，人體所攜帶的金屬槍亮溫設(shè)為150K，人體腿部攜帶的刀具亮溫設(shè)為210 K.考慮到實際系統(tǒng)中，由于有限積分時間引入的可視度測量噪聲.當(dāng)接收機(jī)的頻率響應(yīng)近似為矩形時，測量噪聲的方差可近似為［14］其中，B為系統(tǒng)帶寬;τ為有效的積分時間;TA為天線溫度;TR為接收機(jī)等效溫度;Λ(x)為三角函數(shù).系統(tǒng)為單邊帶結(jié)構(gòu)，因此Λ(2Δf/B)=0.忽略接收機(jī)不理想性，當(dāng)積分時間為0.5 s時，可視度函數(shù)的信噪比約為33.4 dB，分別采用L2范數(shù)、全變分及本文所述方法應(yīng)進(jìn)行圖像反演.仿真結(jié)果分別如圖2(b)、圖2(c)、圖2(d)所示.

圖2 不同反演算法的仿真結(jié)果Fig.2 Simulation results of different reconstruction algorithms

對比結(jié)果可以看出，由于測量噪聲的影響，采用L2范數(shù)反演圖像(圖2(b))的噪聲抑制并不充分，且人體攜帶的刀具幾乎無法識別.這是因為該算法在整幅圖像都采用同一種約束，很難在保持整幅圖像邊緣細(xì)節(jié)和抑制噪聲之間折衷.全變分方法雖然在保持圖像邊緣細(xì)節(jié)的同時具有一定的去噪能力，但由于噪聲抑制不充分，在目標(biāo)區(qū)域及周圍背景區(qū)域產(chǎn)生了階梯效應(yīng).局部自適應(yīng)的方法將待反演圖像分區(qū)，克服了整幅圖像都采用同一種約束的缺點，使得其反演獲得的圖像具有更低的背景噪聲及更突出的邊緣信息，同時有效地抑制了階梯效應(yīng).上述反演圖像中的波紋振蕩是由于有限孔徑截取造成的Gibbs效應(yīng)，可以通過加窗的方式進(jìn)行改善，但會降低空間分辨率.為了更加客觀地描述不同反演算法的性能，采用反演圖像與理想圖像的均方根誤差(TRMSE)來衡量反演圖像的質(zhì)量:

其中，T0為理想圖像;T1為反演圖像.利用上式，分別計算出采用L2范數(shù)、全變分及局部自適應(yīng)方法獲得的亮溫圖像的 TRMSE分別為 10.139 9，9.3593和7.2055.因此，采用本文提出的局部自適應(yīng)方法獲得的亮溫圖像的質(zhì)量均優(yōu)于其他成像算法.仿真中，不同算法的迭代停止準(zhǔn)則均設(shè)置為前后兩次迭代圖像的均方根誤差小于0.1.為了驗證不同測量噪聲水平情況下該反演算法的性能，分別選取積分時間 t為0.1，1，10 s時進(jìn)行仿真，反演結(jié)果的相對均方誤差如表2所示.

表2 不同反演算法的均方根誤差TRMSETable2 TRMSEof different reconstruction algorithms

由表2可以看出，不同測量噪聲水平下，局部自適應(yīng)算法均具有較好的適用性，一定程度上減小了可視度測量噪聲對亮溫圖像的影響.仿真中，正則化參數(shù)μ需根據(jù)不同可視度測量噪聲水平進(jìn)行調(diào)整.另一方面，由于實際中無法準(zhǔn)確獲知目標(biāo)與天線陣的距離信息，只能進(jìn)行近似估計.下面對存在距離誤差的情況下進(jìn)行仿真，成像距離設(shè)為3 m，估計誤差范圍為3±1 m.同樣，在積分時間為0.5 s時仿真結(jié)果如圖3所示.

圖3 存在距離誤差時各反演算法的均方根誤差TRMSEFig.3 TRMSEof different reconstruction algorithms caused by imaging distance error

從圖3可以看出，隨著成像距離誤差增大，圖像反演誤差逐漸增加，且前向距離誤差比后向距離誤差的影響更大.對比不同反演算法的結(jié)果可知，當(dāng)存在有距離誤差時，采用局部自適應(yīng)算法的反演結(jié)果的相對均方誤差均小于其他兩種算法.同時表明，局部自適應(yīng)算法對成像距離誤差具有更好的適應(yīng)性.上述仿真的實驗平臺為一臺CPU型號為3.30GHz Intel(R)Core(TM)i3-2120、內(nèi)存為8 GB的計算機(jī)，仿真軟件為Matlab 2013.在此平臺上，該算法的平均運行時間約為1.6 s.

4 實驗驗證

為了進(jìn)一步驗證本文提出的近場圖像反演算法的有效性，利用現(xiàn)有的一套8mm波段二維綜合孔徑輻射計系統(tǒng)BHU-2D-U(如圖4所示)，對人體進(jìn)行近距離的成像實驗.在進(jìn)行圖像反演之前，需先測得BHU-2D-U的系統(tǒng)響應(yīng)G矩陣.本實驗中，利用一個安裝在機(jī)械掃描架上的外部單點源逐點掃描的方式獲得系統(tǒng)響應(yīng)G矩陣［15-16］，此處系統(tǒng)視場范圍與G矩陣測量范圍一致.然而，由于機(jī)械掃描架尺寸的限制，本實驗中BHU-2D-U的視場范圍被限制為70 cm×70 cm.

圖4 被動毫米波成像系統(tǒng)BHU-2D-UFig.4 Passive millimeter-wave imaging system of BHU-2D-U

實驗場景如圖5(a)所示，被測人員位于BHU-2D-U系統(tǒng)正前方約2 m處，右手在胸前持有一個V型的金屬架.為了形成一個均勻的冷背景，在人體背后放置了一塊與地面大約成45°夾角的鋁板，用來反射天空亮溫.實驗中，積分時間設(shè)置為0.5 s.圖5(b)給出了采用本文提出的局部自適應(yīng)方法獲得的近場毫米波圖像.實驗中，首先利用近場點聚焦FFT算法獲得模糊的亮溫圖像，以劃分待反演亮溫圖像的背景區(qū)域與目標(biāo)區(qū)域.可以看出，雖然在光學(xué)圖像中金屬架的一部分被隱藏在衣服下面，但在毫米波圖像中可以清晰地識別出該金屬架的輪廓.為了對比不同成像方法的性能，分別采用L2范數(shù)、全變分方法進(jìn)行毫米波亮溫圖像反演，其結(jié)果分別如圖5(c)、圖5(d)所示.實驗中，考慮到實際系統(tǒng)誤差的存在，將算法的迭代停止準(zhǔn)則均設(shè)置為前后兩次迭代圖像的均方根誤差小于0.3.

圖5 不同反演算法的實驗結(jié)果Fig.5 Experimental results of different reconstruction algorithms

實驗結(jié)果表明，采用L2范數(shù)的反演圖像由于平滑特性導(dǎo)致圖像分辨率較低，金屬架的輪廓并不清晰.全變分法提高了分辨率，但由于較高的噪聲引入的階梯效應(yīng)非常明顯，金屬架輪廓部分已出現(xiàn)不連續(xù)情況，影響對人體隱匿物品的識別.局部自適應(yīng)方法既去除了背景區(qū)域噪聲又抑制了目標(biāo)區(qū)域的階梯效應(yīng)，具有較好的圖像反演效果.

5 結(jié)論

本文根據(jù)毫米波人體安檢圖像特點，提出了一種基于局部自適應(yīng)偏微分方程的綜合孔徑輻射計近場圖像反演算法.此算法根據(jù)待反演亮溫圖像的分布特性，對不同的區(qū)域采用不同的約束，以達(dá)到去除噪聲且保持圖像細(xì)節(jié)邊緣的目的.對此算法進(jìn)行了仿真和實驗驗證，結(jié)果表明:

1)該算法對近場綜合孔徑輻射計圖像反演具有明顯的降噪和保持邊緣細(xì)節(jié)效果，為安檢應(yīng)用中進(jìn)一步的違禁目標(biāo)分割、識別提供良好基礎(chǔ).

2)該算法在不同噪聲水平下均可將反演圖像性能提升約20%，且對近場成像距離的誤差具有較好的適應(yīng)性，進(jìn)一步增加了該算法的應(yīng)用價值.

然而，為了實現(xiàn)實時快速成像的目的，在以后的工作中需要對本文算法的運算效率進(jìn)行重點研究.

References)

［1］ Kolinko V G，Lin S，Shek A，et al.A passive millimeter-wave imaging system for concealed weapons and explosives detection［C］//Proceedings of SPIE-The International Society for Optical Engineering.Orlando:SPIE，2005，5781:85-92.

［2］ Wikner D A.Progress in millimeter-wave imaging［C］//Proceedings of SPIE-The International Society for Optical Engineering.Bellingham，WA:SPIE，2011，7936:1-9.

［3］ Nova E，Romeu J，Torres F，et al.Radiometric and spatial resolution constraints in millimeter-wave close-range passive screener systems［J］.IEEE Transactions on Geoscience and Remote Sensing，2013，51(4):2327-2336.

［4］ Zheng C，Yao X X，Miao J G，et al.Initial results of a passive millimeter-wave imager used for concealed weapon detection BHU-2D-U［J］.Process In Electromagnetics Research C，2013，43:151-163.

［5］ Peichel M，Suess H，Suess M.Microwave imaging of the brightness temperature distribution of extended areas in the near and far field using two-dimensional aperture synthesis with hight spatial resolution［J］.Radio Science，1998，33(3):781-801.

［6］ Tanner A B，Lambrigsten B H，Gaier T M，et al.Near field characterization of the GeoSTAR demonstrator［C］//IEEE Geoscience and Remote Sensing Symposium.Piscataway，NJ:IEEE，2006:2529-2532.

［7］ Zhang C，Wu J，Liu H，et al.Imaging algorithm for synthetic aperture interferometric radiometer in near field［J］.Science China Technological Sciences，2011，54(8):2224-2231.

［8］ Anterrieu E.A resolving matrix approach for synthetic aperture imaging raidometers［J］.IEEE Transactions on Geoscience and Remote Sensing，2004，42(8):1649-1656.

［9］ 鄒謀炎.反卷積和信號復(fù)原［M］.北京:國防工業(yè)出版社，2001.Zou M Y.Deconvolution and signal reconvery［M］.Beijing:National Defence Industry Press，2001(in Chinese).

［10］ Perona P，Malik J.Scale-space and edge detection using anisotropic diffusion［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，1990，12(7):629-639.

［11］ Laursen B，Skou N.Synthetic aperture radiometry evaluated by a two-channel demonstration model［J］.IEEE Transactions on Geoscience and Remote Sensing，1998，36(3):822-832.

［12］ 胡岸勇，苗俊剛.一種擴(kuò)大8 mm波段綜合孔徑輻射計成像視場的方法［J］.紅外與毫米波學(xué)報，2013，32(1):50-55.Hu A Y，Miao J G.A method to enlarge the FOV for 8 mm band synthetic aperture radiometer［J］.International Journal of Infrared and Millimeter Waves，2013，32(1):50-55(in Chinese).

［13］ Bioucas-Dias J M，F(xiàn)igueiredo A T.A new TwIST:two-step iterative shrinkage/thresholding algorithms for image restoration［J］.IEEE Transactions on Image Processing，2007，16(12):2992-3004.

［14］ Ruf C S，Swift C T，Tanner A B，et al.Interferometric synthetic aperture microwave radiometry for the remote sensing of the earth［J］.IEEE Transactions on Geoscience and Remote Sensing，1988，26(5):597-611.

［15］ Tanner A B，Swift C T.Calibration of synthetic aperture radiometer［J］.IEEE Transactions on Geoscience and Remote Sensing，1988，31(1):257-267.

［16］ Yao X X，Zheng C，Zhang J，et al.Near field image reconstruction algorithm for passive millimeter-wave imager BHU-2D-U［J］.Process in Electromagnetics Research C，2013，45:57-72.