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

姚現勛，尚曉舟，苗俊剛，李志平

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

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

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

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

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

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

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

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

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

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

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

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

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

2.1 偏微分方程正則化方法

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

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

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

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

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

其對應的Euler-Lagrange方程為

2.2 局部自適應偏微分方程方法

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

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

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

其中

其中，Ωb為圖像的背景區域;Ωt為圖像的目標區域;Gσ為高斯濾波器，用來對噪聲圖像進行預處理.式(11)對應的Euler-Lagrange方程為

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

3 近場成像仿真

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

表1 BHU-2D-U系統主要參數Table1 Main specifications of BHU-2D-U

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

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

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

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

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

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

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

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

4 實驗驗證

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

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

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

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

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

5 結論

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

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

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

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

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］ 鄒謀炎.反卷積和信號復原［M］.北京:國防工業出版社，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］ 胡岸勇，苗俊剛.一種擴大8 mm波段綜合孔徑輻射計成像視場的方法［J］.紅外與毫米波學報，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，Figueiredo 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.