基于稀疏貝葉斯學習的太赫茲電磁渦旋三維成像方法

蔣彥雯 范紅旗 李雙勛

(國防科技大學電子科學學院 長沙 410073)

1 引言

近年來，當在傳統電磁波上加載軌道角動量調制時，可在空間中疊加形成渦旋電磁波，利用不同形式天線可產生不同特性的軌道角動量電磁波，從而適用于不同的應用需求[1-3]。特別地，在陣列雷達成像技術中，利用渦旋電磁波照射，同一距離單元內不同方位向的目標所在位置的輻射場相位波前存在一定的差異性，從而使雷達接收到的目標散射回波攜帶有更多的目標信息，進而通過方位向的信息解耦最終可實現目標凝視成像[4-6]。對于渦旋電磁波照射下的目標成像，利用照射波束內波前的差異性，并且基于不同本征模態的正交性，通過多模態的照射可有望提升對目標的分辨性能[7]。在電磁渦旋成像技術的基礎上，結合合成孔徑雷達(Synthetic Aperture Radar,SAR)成像或逆合成孔徑雷達(Inverse Synthetic Aperture Radar,ISAR)成像特性，可進一步實現對目標或場景的三維(Three-Dimensional,3D)成像[8,9]。

在太赫茲頻段，由于雷達信號載頻高，極易發射大帶寬信號，實現較高的距離分辨率[10,11]，另外，相比于微波雷達，由于太赫茲波的波長短，在合成孔徑成像或者逆合成孔徑成像模式下可以獲得更高的橫向分辨率[12-14]。文獻[9]將太赫茲波的優勢與渦旋電磁波相結合，建立了基于太赫茲電磁渦旋ISAR的新體制成像模型，實現了對目標的高分辨三維成像，然而在已有的成像算法中，俯仰向信息需要通過極坐標下三角函數和方位向的插值得到，這是一個間接求解的過程，俯仰維成像分辨率較低。

為進一步提高目標俯仰維的成像分辨率，本文在前期工作的基礎上[9]，提出了一種基于稀疏貝葉斯學習(Sparse Bayesian Learning,SBL)的圖像重建方法，可直接在直角坐標系下求解目標的三維散射分布函數。利用稀疏貝葉斯學習方法進行圖像重建時，通過稀疏表示模型建立觀測矢量與待重建信號之間的關系，進而直接進行求解，求解過程能夠自動學習重建模型中的未知參數，且人工參數設置對重建結果的影響較少[15]。文中首先推導建立了電磁渦旋三維成像的回波模型；其次，構建了目標三維成像的稀疏表示模型，對SBL圖像重建過程進行了詳細介紹；最后，通過設置不同的成像實驗場景，對不同成像方法的成像性能進行了對比分析，并對不同信噪比(Signal-to-Noise Ratio,SNR)條件下SBL方法的重構性能開展了仿真實驗。

2 電磁渦旋三維成像模型與方法

2.1 電磁渦旋三維成像模型

在傳統ISAR成像模型下，利用帶寬信號和雷達與目標之間的相對運動可實現距離-方位二維成像。另外，在渦旋電磁波的照射下，結合帶寬信號可獲得距離-方位角的二維分辨能力。因此，本文將二者結合，建立了基于電磁渦旋ISAR的三維成像模型，如圖1所示，圖中均勻圓陣表示雷達，通過在線性調頻發射信號上加載軌道角動量調制而產生渦旋電磁波，假設目標上一理想散射點P的坐標可以表示為(x,y,z)，且目標繞Y軸旋轉，旋轉方向如圖中箭頭所示，轉角為θ∈[-Δθ/2,Δθ/2],Δθ表示目標旋轉的最大角度，圖中方位角φ為理想散射點P在XOY平面內的投影與X軸的夾角。

圖1 基于電磁渦旋ISAR的三維成像幾何Fig.1 Sketch map of the 3D imaging geometry based on electromagnetic vortex ISAR

一般來說，當目標由大量散射點構成，雷達接收到的目標總回波可表示為式(1)積分形式[9]

其中，f(x,y,z)表示目標三維散射分布函數，k為波數，l為軌道角動量模式數。

當雷達和目標之間轉動角度較小時，對式(1)所示雷達回波分別在頻率維k和方位維θ進行二維傅里葉逆變換，再在拓撲荷域l進行一維傅里葉變換，即可得到目標的散射分布函數

其中，f(x,z,φ)為方位-距離-方位角上的目標三維散射分布函數。因此，在已知目標方位x、距離z、方位角φ信息的情況下，要獲得直角坐標系下的目標三維散射分布函數，即f(x,y,z)，首先需要根據方位x、方位角φ的值計算得到俯仰向Y的分布，計算表達式為y=x·tanφ，然后在Y軸上將計算得到的非均勻網格插值為均勻網格，最終可得f(x,y,z)，實現目標三維成像[9]。

2.2 基于稀疏貝葉斯學習的成像方法

根據2.1節中建立的成像模型，式(2)中的成像過程可改寫為

首先，對雷達回波采用卷積逆投影(Convolution Back-Projection,CBP)的方法進行成像，令G(k,θ,l)=k·Sr(k,θ,l)，且g(w,θ,l)為G(k,θ,l)的一維傅里葉逆變換。因此，式(3)中的第1步積分可表示為

其中，h(x,z,l)為不同軌道角動量模式數時目標二維散射分布函數，在不同θ值計算式(4)中對應的-xsinθ+zcosθ，然后將g(w,θ,l)沿不同的θ進行疊加，即可計算得到h(x,z,l)，這就是CBP成像算法的基本步驟[16]。

接下來，將式(4)代入式(3)，可以得到

式(5)完全符合傅里葉變換形式，因此，利用傅里葉變換與逆變換相對應的性質可知

將φ=arctan(y/x)代入式(6)等號右側的積分表達式中，可得

下面，再對X方向和Z方向的二維成像網格進行離散化，將相應方向的成像區域分別劃分為M,N個網格，在每一個離散值xm,zn對應的目標二維散射分布函數h(x,z,l)處，采用SBL方法進行一維重構，成像幾何如圖2所示，將式(7)進一步改寫為

圖2 基于SBL方法的三維成像幾何Fig.2 3D imaging geometry based on the SBL method

根據稀疏貝葉斯恢復的基本原理，在式(8)的基礎上構建如下信號模型

其中，H表示觀測矢量，A為對應的測量矩陣，y為 待求解的俯仰向Y的散射系數矢量，n表示噪聲。假設雷達發射信號的軌道角動量調制模式數的采樣點數為L，俯仰向的離散網格數為P，將式(9)改寫為

式(10)中忽略了xm,zn這兩個重復項，測量矩陣A可進一步表示為

在每一個xm,zn處，重復進行式(10)中的稀疏貝葉斯重構，最終可直接得到直角坐標系下的目標三維散射分布函數f(x,y,z)。

在上述重構過程中，xm,zn各有總共M,N個取值，這樣就需要完成M×N次貝葉斯重構，當M,N取值較大時，重構過程計算復雜度高且算法運行時間較長。在實際成像場景中，目標散射點通常是成稀疏分布的，并不是所有俯仰向都存在目標散射點，因此，為降低計算量，只需要在有目標散射點的距離-方位切片上進行一維重構。當l=0時，雷達發射信號僅為線性調頻信號，此時，h(x,z,0)為傳統ISAR成像結果，即目標在XOZ平面上的投影，當h(xm,zn,0)的幅度很小時，本文認為在xm,zn所在的距離-方位切片上不存在目標散射點。因此，本文對實際成像過程作進一步優化以節約計算成本，與文獻[17]中不同的是，為避免距離-方位切片上強散射點的影響，本文提出分區域幅度閾值設置方法，首先在距離-方位切片上尋找局部最大值點h(xm,zn,0)，以各局部最大值點h(xm,zn,0)為中心劃分W ×W個分辨單元為選定區域 Ω,W通常選擇為1/2距離-方位向點擴散函數主瓣寬度所占的分辨單元個數，設定幅度閾值為η·h(xm,zn,0)。當該區域內時，記為h(xmi,zni,0)，該距離-方位切片包含目標散射點，將h(xmi,zni,l)作為觀測矢量，構建測量矩陣，求解f(xmi,zni,y)；當時，不包含目標散射點，該距離-方位分辨單元內俯仰向散射強度均設置為0。最終，優化后的SBL成像處理流程如下。

Step 1：設置目標散射點，根據式(1)生成三維采樣下的雷達回波Sr(k,θ,l)；

Step 2：在不同軌道角動量模式數l=l1,l2,...,lN，分別對雷達回波進行二維CBP成像，得到目標二維散射分布函數h(x,z,l)；

Step 3：在距離-方位切片h(x,z,0)內，尋找所有局部最大值h(xm,zn,0)；

Step 5：采用稀疏貝葉斯學習的方法進行一維重構得到f(xmi,zni,y)；

Step 6：計算得到的所有離散值f(xm,zn,y)，最終得到目標三維成像結果。

3 仿真結果與分析

3.1 實驗場景設置

仿真實驗中，假設雷達發射信號中心頻率為330 GHz，信號帶寬為20 GHz，頻率采樣間隔為0.1 GHz。圓形陣列的陣元半徑設置為500λ，最大軌道角動量模式數取lmax=30，相應的軌道角動量模式數的變化范圍為[-30,30]，成像時設置方位向最大轉角為 Δθ=5°，對應轉角θ∈[-2.5o,2.5o]，仿真時設置方位向成像場景寬度為0.6 m，對應滿足不混疊距離的角度采樣點數可設置為126，設置各散射點的散射強度均為1，其幾何位置關系如圖3所示。

3.2 電磁渦旋三維成像結果

成像仿真中，本文的稀疏貝葉斯求解方法采用變分貝葉斯推斷(Variational Sparse Bayesian Inference,VSBI)[18]的方法，首先對圖3所示的目標散射點根據式(1)生成雷達回波，然后再采用本文提出的SBL方法進行三維成像，成像過程中，設置X,Y,Z3個方向的離散網格數均為M=N=P=601。

圖3 目標散射點分布Fig.3 The distribution of point targets

圖4 l=0時距離-方位切片及不同幅度閾值劃分結果Fig.4 The range-azimuth image atl=0 and the results of different amplitude threshold setting method

圖4(a)為l=0時的歸一化的距離-方位切片h(x,z,0)，采用分區域幅度閾值設置方法時，根據點擴散函數的主瓣寬度，設置η=0.7，得到如圖4(b)所示的分區域幅度閾值劃分結果。區域劃分時，將l=0時的歸一化的距離-方位切片視為一幅圖像，采用圖像處理中尋找局部最大值的經典八連通方法進行處理[19]，首先會找到散射點P1所在分辨單元的局部最大值點，以該點幅度值的η倍為閾值，根據距離分辨率和方位分辨率與離散網格數計算選取W=40，對選定區域 Ω內的分辨單元作幅度閾值劃分，保證了后續稀疏貝葉斯重構方法的準確性。若采用文獻[17]中的幅度閾值方法，為所有h(x,z,0)設置單一幅度閾值，結合圖3和圖4(a)可以看出，除散射點P1之外，距離-方位平面上的其他散射點均由2個目標散射點投影形成，使得散射點P1的相對散射強度較弱，僅為0.35(遠低于1)，若設置閾值η=0.7，最終三維成像結果中將不包含散射點P1，而若將幅度閾值降低，設置為0.3，得到如圖4(c)所示的劃分結果。對比圖4(b)和圖4(c)可知，分區域幅度閾值設置方法得到的各散射點所在區域分辨單元個數是大致相同的，而圖4(c)中雖然能夠得到距離-方位切片上5個散射點所在區域，但P1點所在的分辨單元個數明顯少于其他散射點，這樣會造成P1點回波能量的大量損耗，進而影響后續重構結果。另外，對比圖4(b)和圖4(c)中選定區域分辨單元總數，圖4(c)中除P1外其他4個散射點由閾值0.3得到的區域面積過大，圖4(c)的分辨單元總數約超過圖4(b)中一倍，極大地增加了后續稀疏重構的計算量。

在分區域輻射閾值設置方法的基礎上，圖5為采用本文提出的VSBI方法重構的三維成像結果，為能直接獲得更清晰的觀測效果，將目標的三維成像結果分別投影到距離、方位和俯仰3個不同坐標平面內，得到不同的二維成像結果。在利用式(11)建立測量矩陣時，各矩陣元素的幅值根據散射點位置的不同而存在一定差異性，從而使得最終圖5中的重構得到的各散射點之間的幅度不盡相同。最后，從圖5中可以看出，包括散射點P1在內的所有散射點位置均得到了準確重構，驗證了本文所提的成像方法的有效性。

圖5 基于SBL方法的三維成像結果Fig.5 3D imaging results based on SBL method

3.3 不同成像方法的對比分析

本文提出的SBL方法對成像效果的提升主要體現在俯仰向Y的重構上，為進一步說明本文方法的優勢，下面將其與文獻[9]中所提的基于三維快速傅里葉變換(Fast Fourier Transform,FFT)和基于CBP和功率譜密度(Power Spectrum Density,PSD)估計兩種成像方法進行對比分析。仿真實驗中，在同一距離-方位分辨單元內設置不同俯仰向位置的4個點，散射點的三維坐標分別為Q1(0.15 m,0,0),Q2(0.15 m,0.01 m,0),Q3(0.15 m,0.1 m,0),Q4(0.15 m,0.115 m,0)，散射點分布如圖6(a)所示，各散射點X軸坐標相同，間隔分布在Y方向。

圖6 目標散射點分布及其成像結果對比Fig.6 The distribution of targets and comparison of imaging results

理想條件下，在x=0.15 m,z=0時的俯仰維剖面圖對比如圖6(b)所示，對比圖中成像結果可知，對于距離很近的Q1和Q2兩點、Q3和Q4兩點，文獻[9]中提出的基于三維FFT和基于CBP-PSD的兩種成像方法均無法在俯仰向實現分辨，而本文提出的SBL方法能夠準確地重構出Q1,Q2,Q3和Q44個散射點的位置，且成像分辨率優于0.01 m。

為進一步衡量SBL成像方法的性能，計算目標散射點的重構誤差，首先，定義散射系數的最小均方誤差(Mean Square Error,MSE)

圖7 重構散射系數MSE隨信噪比變化情況Fig.7 MSE of reconstructed scattering coefficient as a function of the SNR

4 結束語

本文將太赫茲ISAR與電磁渦旋相結合，建立了基于電磁渦旋ISAR的成像稀疏表示模型，并提出了基于稀疏貝葉斯學習的目標三維重構方法，通過SBL方法可直接由雷達回波重構目標在空間直角坐標系中的三維信息，極大地簡化了成像求解過程，提高了成像分辨率。后續研究將針對渦旋電磁波所特有的貝塞爾函數幅度特性并結合實際電磁渦旋雷達實驗系統進行深入分析，為新體制雷達三維成像技術的發展提供參考和借鑒。