基于正交頻譜重構的微波三維成像方法

張 研, 王保平, 方 陽, 宋祖勛

(1.西北工業大學電子信息學院, 陜西 西安 710072;2.西北工業大學無人機特種技術重點實驗室, 陜西 西安 710065)

0 引 言

微波三維成像技術是一種工作在300 MHz～300 GHz波段的三維合成孔徑雷達(synthetic aperture radar, SAR)成像的技術[1-2],通過發射寬帶電磁波信號獲得距離向高分辨率,在方位向和高度向以合成孔徑的方式對目標觀測獲得方位向和高度向高分辨率。因其具有高分辨、低輻射、非接觸、三維立體成像等優點,在民用和軍用領域開展了廣泛的研究。在民用領域主要可應用于安檢、醫學診斷、災難救援、無損檢測等領域;在軍用和反恐領域主要可應用于戰場隱匿物體探測、重要軍事目標勘察、地雷探測等領域[2-14]。

微波成像的分辨率取決于發射信號帶寬、中心頻率和合成孔徑大小。受香農-奈奎斯特采樣定律的限制,傳統高分辨成像方法要求成像系統具有非常高采樣率,這對雷達系統的硬件處理帶來很高的難度。此外,三維成像需要在3個方向采集數據,導致回波數據量非常大,對數據處理也帶來一定挑戰[2,15-18]。將壓縮感知(compressed sensing,CS)技術[19-20]引入微波三維成像可以大大降低成像系統的采樣率和回波數據量,以欠采樣數據恢復出精確的目標圖像[21-23]。將稀疏度等先驗知識作為約束項的正則化CS成像方法進一步提高了成像精度[24-25]。總差分(total variation, TV)正則化CS成像方法以散射點之間的能量差進行稀疏表示,可以提高成像的抗噪性能,但會降低成像精度[26-27]。為了同時提高成像精度和抗噪性能,一種基于TV和l1范數的混合約束正則化CS成像方法被提出,該方法通過迭代循環調整TV約束和l1范數約束的權值系數,實現最優的成像結果[28-30]。目前的CS成像方法主要基于目標場景的稀疏表示,由于需要計算整個目標場景最優解,導致字典矩陣的內存占用巨大,對硬件處理帶來困難。此外,基于目標場景稀疏的CS成像結果是由離散的強散射點組成,對連續分布目標的細節特征體現較差,不利于目標識別。

本文首先分析了目標回波的頻域空間分布特性和回波信號模型,討論了目標回波頻譜對成像的影響。然后,根據所構建的目標回波波陣面的正交投影模型實現了目標回波頻譜的最優重構。最后,以成像結果的最小圖像熵為判別準則對重構的正則化參數進行了最優估計,進一步提高了頻譜重構精度和成像質量。計算機仿真實驗和微波暗室實測實驗表明,本文所提方法可以處理稀疏采樣數據,內存占用較CS方法大大降低,且具有較高的成像精度,尤其對目標細節特征的恢復較為準確。

1 成像模型分析

本文所討論的是平面掃描微波三維成像系統,采用直角坐標系表示成像幾何關系,如圖1所示。以目標場景中心o為坐標原點建立坐標系T(o,x,y,z),以掃描平面中心o′為坐標原點建立坐標系T′(o′,x′,y′,z),其中x′o′y′平面與平面xoy平行。雷達掃描平面的長和高分別為Lx和Ly,此即為成像系統在x′方向和y′方向的孔徑長度,目標成像場景中心o與掃描平面的法向距離為z0。成像系統采用單發單收方式進行信號采集,以步進頻方式發射寬頻帶電磁波信號,發射波瞬時頻率為f=f0+Δf,f0為初始頻率,Δf為頻率步進間隔。

圖1 成像幾何模型

設rt∈R3為目標場景中心到某一散射點P的距離矢量,ra∈R3為雷達天線到目標場景中心的距離矢量,則散射點與雷達天線之間的距離矢量為

rP=ra-rt

(1)

假設電磁波傳播介質為線性、各向同性、和均勻的,則在天線掃描位置(x′,y′)處接受到的散射點P(x,y,z)的反射波的電場強度標量可以表示為

E(x′,y′,x,y,z,k)=δ(x,y,z)e-jk|rP|

(2)

式中：k=4πf/c表示波數,c為光速;δ(·)表示某一位置上散射點的反射率(假設目標后向散射為各向同性的)。對于三維成像而言,可以將整個成像場景的回波近似為成像場景內所有散射點回波的線性疊加。因此,對于天線掃描位置(x′,y′),接收到的目標場景回波信號可以表示多個散射點在(x′,y′)處電場強度的積分,其表達式如下:

s(x′,y′,k)=

(3)

式中：R(x′,y′,x,y,z)表示散射點與雷達天線之間的瞬時斜距,其表達式如下:

R(x′,y′,x,y,z)=|rP|=

(4)

2 成像處理流程

2.1 頻譜重構

對式(3)在x′方向和y′方向進行傅里葉變換,可以將回波信號變換到三維頻域,其表達式如下:

s(kx,ky,kz)=

(5)

一種理想的處理方式是對上式分別在kx,ky和kz方向上進行傅里葉逆變換,即可得到目標場景的反射率分布,其表達式如下:

? {δ(x,y,z)e-j[kxx+kyy+kz(z+z0)]·

(6)

直接通過傅里葉逆變換運算效率太慢,一般計算機處理傅里葉變換使用的是快速傅里葉逆變換(inverse fast Fourier transform,IFFT)。IFFT要求數據具有均勻采樣的格式。但在三維頻域,如式(6)所示的回波信號頻譜在kx-ky-kz坐標系下不具有均勻分布的性質,這是因為kz與k,kx和ky之間具有如下關系:

(7)

則式(5)變為

s(kx,ky,k)=

(8)

圖2 相位面投影示意圖

(9)

則回波頻譜變為

(10)

式(10)經過IFFT就可以得到目標反射率分布。

通過變量替換實現頻譜投影最常用的方法是插值,但是插值存在一定的誤差和近似,而插值的過程本身也不是一個優化過程。此外,如果雷達成像系統在掃描平面和發射信號均采用隨機稀疏采樣的方式進行數據采集,則存在兩個問題,第一個問題是數據不滿足香農-奈奎斯特采樣定律;第二個問題是數據為非均勻采樣。這兩個問題會導致原始回波數據s在x′方向和y′方向無法進行傅里葉變換,且插值精度會較差。

為了應對稀疏采樣以及通過優化重構的方法實現數據的投影變換,首先需要推導由回波信號到待重構頻譜之間的矩陣運算關系。設sP×Q×R為原始回波數據矩陣,P,Q和R分別為x′方向,y′方向和z方向采樣點數,回波信號的頻譜可以表示為以下矩陣運算：

(11)

(12)

(13)

(14)

(15)

(16)

(17)

式中：λ為正則化參數,用以調節重構精度和稀疏性之間的平衡。λ的取值與先驗知識相關,如待重構信號的稀疏度等。通過凸優化方法求解式(17)時,由于λ在實際應用中往往難以提前測量和估計,為了保證重構結果盡可能接近最優解,λ一般在0～1之間取較小的值。凸優化方法處理稀疏性不是很好的信號時仍能保持較高的精度,但其在處理大規模優化模型時效率太低,如三維SAR成像。因此，本文采用正交匹配追蹤(orthogonal matching pursuit, OMP)算法進行重構。

在正交匹配追蹤算法中,迭代終止條件依賴于對帶重構信號稀疏度的估計,這將在第2.2節具體討論。為了便于矩陣運算,將F拉伸成一維矩陣,設χλ為迭代的殘差矩陣,th為判別迭代是否終止的殘差閾值,t為迭代次數。具體重構流程如下所示。

算法1 目標波譜正交重構算法輸入:F,Φ初始化:χλ=F,t=0迭 代:步驟1計算最大投影向量:Φη=max(ΦηF),η=1,2,…步驟2計算正交投影系數:αη=ΦTηFiΦTηΦη步驟3將重構結果填充在矩陣中:F′η=αη步驟4計算殘差χλ=F-αηΦη步驟5判斷:如果|χλ|≤th,停止迭代;否則,返回步驟1。 輸出:F′

對F′作3D-IFFT就可以得到目標的反射率分布。

2.2 基于最小圖像熵的參數優化

通過對目標頻譜進行最優重構,然后再通過3D-IFFT得到目標的反射率分布(即目標圖像)雖然可以在一定程度上提高成像效果,并且可以處理稀疏采樣數據。但前提是選取合適的殘差閾值th來判別迭代是否終止。不同的目標,不同成像環境中的雜波干擾,參數誤差等都會影響最優迭代次數。因此，需要通過對殘差閾值th進行最優估計。

正交投影可以消除目標頻譜在kx,ky和kz方向上的耦合,使得目標精確聚焦。根據微波成像機理和信息熵理論,只有當成像參數為最優參數時,目標的聚焦性最好,此時目標圖像具有最小的熵。因此，將最小圖像熵作為判定準則,對正則化參數進行估計,可以求解出最優的正則化參數,并實現最優成像結果。圖像熵的定義如下所示：

(18)

式中：

(19)

式中：I(x,y)為二維成像結果。由于頻譜正交投影重構是逐波數采樣進行的,因此只能通過對x′方向和y′方向進行IFFT得到目標的二維成像結果。在實際處理流程中對每次正交重構的殘差閾值th進行調整,求出各th候選值對應成像結果的圖像熵,其最小值對應的th和頻譜F′就是我們所需要的最優解。為了提高計算速度,可以采用梯度下降算法或者牛頓法等進行求解。

定義imax為的采樣點數,γ為迭代循環次數,ε為梯度下降算法的誤差,Δth0為正交投影重構的誤差閾值的初始步進值。根據最小圖像熵進行參數優化的完整成像處理流程如圖3所示。

圖3 成像整體流程圖

3 仿真與實驗

本節根據平面掃描三維成像模型構建成像實驗系統,分別通過計算機生成的仿真實驗數據和在微波暗室測量獲取的真實回波數據對本文所提方法進行驗證。計算機仿真數據實驗主要對所提方法的關鍵步驟進行定性定量分析。微波暗室測試數據主要驗證本文所提方法對真實連續目標的成像效果。根據CS理論中重構字典的非相關條件,本實驗所有的稀疏采樣均為隨機稀疏采樣,采樣位置的隨機分布通過M序列實現。

3.1 點目標仿真分析

在計算機仿真實驗中,目標如圖4所示,為3個點目標,成像系統主要參數設置如表1所示。

圖4 點目標模型

表1 成像仿真參數

從圖5(a)可以看出25%隨機稀疏采樣條件下傳統頻域插值成像方法完全散焦,說明傳統頻域插值方法無法處理隨機稀疏采樣數據。從圖5(b)可以看出所提方法僅產生少量散焦情況,仍然較好地對目標的三維空間分布進行了表征。

圖5 25%采樣數據成像結果

頻譜變換的精度對成像質量至關重要,對回波信號分別在電磁波發射方向和掃描平面進行稀疏采樣,分析不同稀疏采樣率下的均方誤差(mean square error, MSE),如圖6所示。由于稀疏重構結果與完整頻譜相比會存在一些缺損,但少量缺損不影響成像,因此正交重構的MSE在計算誤差時只計算非零值點。從圖6中可以看出, sinc插值和非均勻快速傅里葉變換(non-uniform fast Fourier transform,NUFFT)的MSE在稀疏采樣率大于70%時增長較小,在稀疏采樣率小于70%時迅速增加。本文所提方法的MSE在采樣率大于30%時增長比較緩慢,在稀疏采樣率小于30%時迅速增加。

圖6 不同數據采樣率對頻譜重構的影響

3.2 微波暗室測試實驗分析

本研究搭建了微波暗室半物理仿真平臺進行了平面掃描3D成像測試系統,測試系統如圖7(a)所示。采用矢量網絡分析儀、功率放大器、發射探頭和接收探頭來模擬一臺雷達系統,模擬雷達系統放置在可升降掃描架上,掃描架被放置在一條導軌上,通過掃描架的升降和沿導軌移動實現平面掃描。測試采用逐列掃描和單發單收的方式,成像參數設置如表2所示。

表2 微波暗室成像參數

目標選用一把刀、2個礦泉水瓶、以及2個易拉罐組合成一個目標體,如圖7(b)所示。這些物體被用塑料膠帶固定在一個塑料泡沫支架上,刀柄和刀刃均為金屬;2個礦泉水瓶的大小相同,但所裝水的體積不同,左側礦泉水裝滿水至瓶口,右側大約裝了3/4瓶水;2個易拉罐沒有裝任何物體。

圖7 微波暗室測試方案圖

首先對本文所提方法在100%采樣數據下的有效性進行驗證,分別采用基于16點SINC插值的頻域插值方法,基于OMP重構的CS成像方法,基于l1-TV重構的CS成像方法和本文所提方法,成像結果如圖8所示。圖8(a)為頻域插值方法成像結果,由于易拉罐的反射率較低,為了顯現出易拉罐的成像結果,采用圖像能量峰值下降6 dB作為顯示門限。從圖中可以看出,刀和礦泉水瓶的成像結果存在一定變形,尤其是刀刃的尖端部分和左側礦泉水的瓶口部分;右側礦泉水瓶比左側礦泉水瓶高度略低,這是因為塑料的電磁波反射能力很弱,而水的電磁波反射能力較強,成像結果實際上是瓶中的水;由于采用較低的圖像顯示門限,雖然易拉罐的反射能量可以顯示出來,但是導致兩個礦泉水瓶中間位置出現了殘影。

在CS成像處理中,成像場景大小為0.6 m×0.4 m×0.4 m,成像場景離散網格化為31×21×21。分別采用了OMP算法和l1-TV重構算法。OMP成像結果如圖8(b)所示,刀、礦泉水瓶和易拉罐均可以顯示出來,沒有旁瓣對成像結果造成干擾。但由于成像結果為離散點組成,只能體現出目標的大致輪廓,對目標形狀的細節體現較差,不利于目標識別。l1-TV重構成像結果如圖8(c)所示,雖然相比OMP算法提高了成像結果的連續性,但是由于算法對目標反射率進行平滑處理,降低收斂速度,使得成像結果精度下降,出現了較多雜散點。

圖8 100%采樣數據成像結果

圖8(d)是本文所提方法成像結果。成像結果采用了圖像峰值下降2 dB作為顯示門限。本文所提方法成像結果與真實目標外形更加接近,可以清晰的看出刀刃尖端和左側礦泉水瓶口的外形特征;兩個易拉罐在圖中清晰可見,并且幾乎沒有雜散點對成像結果造成干擾。

為了驗證本文所提方法在稀疏采樣條件下的成像能力,分別用頻域插值方法、OMP重構方法、l1-TV重構方法和本文所提方法對25%稀疏采樣的數據進行成像處理,成像結果如圖9所示。從圖9(a)可以看出,頻域插值方法圖像完全散焦,說明頻域插值方法無法處理稀疏采樣數據;從圖9(b)可以看出,基于OMP重構的CS成像結果在25%稀疏采樣數據下基本可以聚焦出目標形狀,但發生了一定的變形,且存在少量雜散點;從圖9(c)可以看出,基于l1-TV重構的CS成像結果目標圖像有少量變化,但是產生了較多的雜散點;從圖9(d)可以看出,本文所提方法在25%稀疏采樣數據條件下,仍然可以較好的聚焦出目標形狀,且相比100%采樣數據條件圖像質量的下降小于其他3種方法。

圖9 25%采樣數據成像結果

為了對實驗結果進行定量分析,分別采用峰值信噪比(peak signal to noise ratio, PSNR)和結構相似度(structural similarity, SSIM)對成像質量進行分析,并比較不同算法的計算內存占用,數據采用100%采樣數據,如表3所示。比較4種方法的PSNR可以看出,基于OMP重構的CS成像結果具有最高的PSNR,這是因為OMP只提取最強的若干個散射點。本文所提方法的PSNR略低于OMP重構方法,l1-TV重構方法和頻域插值方法的PSNR較低,因為l1-TV重構方法和頻域插值方法均存在一些雜散點,而所提方法成像結果基本無雜散點。

表3 成像性能對比

對于SAR成像,SSIM更能符合人眼視覺對成像效果的要求。從表3中可以看出,本文所提方法具有最高的SSIM,其次是l1-TV重構方法和頻域插值方法,而OMP重構方法的SSIM最低。

分析這幾種方法的內存占用可以看出,頻域插值方法的占用內存最小,由于不需要構建字典矩陣,因此內存基本等于數據大小;本文所提方法占用內存略大于頻域插值方法,由于本文所提方法的字典矩陣只針對每一個波數面構造,因此字典矩陣較小;OMP重構方法和l1-TV重構方法的占用內存遠大于前兩種方法,因為這兩種方法都需要根據整個成像場景構造字典矩陣,且在三維成像模式下,字典矩陣大小會呈指數增加,嚴重增加硬件處理難度。

幾種成像方法的計算時間差別巨大,傳統頻域插值方法速度非常快;OMP重構方法具有快速收斂的優點,但由于字典矩陣巨大,因而處理時間比頻域插值法有大幅度增加;l1-TV重構方法收斂速度較慢,且有復雜的優化迭代過程,因此處理時間最長;所提方法通過逐次對頻譜列向量進行正交重構,大幅度降低了字典矩陣大小,雖然計算時間比傳統頻域插值方法長,但比OMP和l1-TV重構方法具有明顯減少。

圖10為不同數據采樣率下的成像結果SSIM曲線。從圖10中可以看出在100%采樣數據條件下,本文所提方法具有最高的SSIM值,之后依次是CS-l1-TV成像方法,頻域插值成像方法和CS-OMP成像方法。隨著數據采樣率的不斷下降,在100%數據采樣率至30%數據采樣率的區間內,本文所提方法和兩種CS成像方法均保持了較好的穩定性,SSIM值下降較為緩慢;而頻域插值方法的SSIM值下降非常明顯;當數據采樣率低于30%時,本文所提方法和兩種CS成像方法的SSIM值開始出現快速下降,當數據采樣率為20%時SSIM均在0.5以上,說明此時圖像相比于目標出現了較大變形,但仍能基本聚焦,當數據采樣率為10%時4種算法的SSIM值均小于0.4,說明此時圖像已經完全散焦。結合圖8,圖9和圖10說明本文所提方法相比于頻域插值方法,CS-l1-TV成像方法和CS-OMP成像方法,與真實目標最為相似,且具有較好的稀疏采樣成像能力。

圖10 稀疏采樣成像相似性分析

4 結 論

本文提出了一種基于正交頻譜重構的微波三維成像方法,該方法采用垂直投影模型對目標的頻譜進行表征,然后根據最小圖像熵準則對正則化參數進行了優化估計,通過正則化稀疏重構實現了目標波譜的最優重構,最后通過對重構后的目標場景頻譜做三維IFFT得到高質量的成像結果。計算機仿真實驗和微波暗室測試實驗表明,所提方法具有較好的稀疏采樣成像能力,對目標幾何結構重構較為準確,且具有較低的內存占用。