基于同心方形網格插值處理的柱面SAR成像算法

丁 麗 何華港 王 韜 儲得苗

①(上海理工大學光電信息與計算機工程學院 上海 200093)

②(上海理工大學生物醫學工程研究所 上海 200093)

1 引言

柱面毫米波合成孔徑雷達(Cylindrical millimeter-wave Synthetic Aperture Radar, CSAR)采用圓跡掃描實現了目標場景360°無盲區檢測[1]，相比于傳統直線軌跡SAR, CSAR具有高空間覆蓋和空間分辨率、高信噪比等優勢，在安檢成像、目標檢測等領域得到了廣泛的研究[2-6]。目前，CSAR成像算法主要基于傅里葉變換理論，需要通過2維插值解決柱面掃描引入的方位維和距離維向波數域非均勻填充的問題，但是這兩個維度的強耦合性，采用2維逐點遍歷插值的算法復雜度高，導致整個成像算法實時性差。

文獻[7]首次提出基于波前重建的CSAR成像算法，為后續柱面成像算法的發展奠定了理論基礎。文獻[8]提出的CSAR成像算法利用球面波分解理論，通過傅里葉變換、相位補償、2維插值等頻域處理方式得到了目標3維圖像。文獻[9]提出了改進的3維CSAR成像算法，通過對高度維進行運動補償，消除了陣列天線在高度維運動造成的運動誤差，提高了成像性能。文獻[10]提出了一種降維時域相關法用于CSAR，在距離域和方位域采用相干累加，使3維時域相關降維至2維時域相關，但是該算法需要對空間目標逐點進行相干累加成像，成像效率低，無法滿足實時成像的需求。文獻[11]在聚束SAR成像處理中利用Chirp-Z變換(Chirp-Z Transform, CZT)實現方位維和距離維重采樣替代2維插值，提高了算法效率。但是CSAR由于采用圓跡掃描使得徑向波數域是非均勻采樣，而CZT變換具有從均勻數據重采樣到均勻數據點的特點，使得無法使用CZT處理CSAR中的2維插值。因此，對于CSAR中波數域方位維與距離維的2維非均勻性和強耦合性造成成像實現效率低、難以開展并行處理的問題仍需進一步研究。

為解決方位維和距離維的強耦合性，本文提出一種基于同心方形網格處理的插值方法?；谥娉上裰胁〝涤驍祿植继攸c的分析，首先，在每一個高度維下對方位維和距離維波數域數據在極坐標系下進行徑向擴充補零，得到膨脹的虛擬同心圓環帶，實現波數域頻帶利用率最大化。接著，將補0后分布在極坐標系下的數據沿著徑向1維插值到同心方形網格中。在進行徑向1維插值之后，原始波數域方位維和距離維的2維強耦合已經可以在直角坐標系下進行分解，因此沿著波數域方位維×距離維橫截面方形的對角線進行區間劃分，把數據劃分為兩個正交不交疊的基本垂直點數據區間和基本水平點數據區間，至此實現對方位維和距離維的解耦。再分別進行方位維和距離維的獨立1維插值，使數據重采樣到均勻網格中，避免了傳統2維插值在遍歷2維網格時所耗的時間復雜度和在待插值點查找局部鄰域所需的計算復雜度等難題，大大減小插值的整體運算量。最后，將兩部分數據進行拼接得到直角坐標系下均勻分布的波數域數據，利用快速傅里葉逆變換重構目標圖像。實驗證明，改進算法插值處理速度較傳統算法提高了7倍，符合算法復雜度理論推導，證明所提算法能大幅度提高成像效率。

2 柱面SAR成像算法

圖1(a)為柱面3維成像系統模型。以場景中心O為原點建立直角坐標系，X軸、Y軸和Z軸分別表示為方位維、距離維和高度維。天線陣列的高度為Lz，以R為半徑進行圓周掃描，θ表示某時刻天線旋轉方位角，且θ ∈[π,-π]。設天線陣列任意采樣的位置坐標為(x′,y′,z′)=(Rcosθ,Rsinθ,z′)，某一散射點目標的位置坐標為(x,y,z)=(rcosφ,rsinφ,z)，對應的后向散射系數為σ(x,y,z)。

圖1 柱面3維成像系統模型

在實際系統中，CSAR回波信號是體目標在3維空間內散射系數的體積分，信號的幅度隨著距離的衰減忽略不計，因為它對圖像的聚焦產生很小的影響，所以在(θ,w,z)域中接收到目標的回波信號為[12]

其中，k為波數，令Φ為相位

對式(1)兩邊做高度維和角度維的傅里葉變換得

其中，kθ和kz分別為波數域角度維分量和高度維分量。根據駐定相位法，對相位Φ沿z′方向求導取0點得

根據式(4)，令ε為天線到目標視線的俯仰角，則可以得到

因此利用正切與正弦三角函數關系解得駐相點z′

將式( 6 ) 代入式( 3 ) 中得Φ=-kr同理可求θ的駐相點，對相位Φ沿θ方向進行求導并令其為0可得

根據式(7)，將CSAR沿著高度維投影到XOY平面上，其投影后的空間幾何關系如圖1(b)所示，其中，天線到目標直線分別與天線到原點直線之間的夾角，和與目標到原點直線之間的夾角記為α和β，且有θ=π-α-β+φ。利用圖1(b)中所示三角形的幾何關系可得

根據三角函數變換求得駐相點θ

將式(9)代入到Φ中得

令Φ=Φ1+Φ2，其中

將式(11)代入式(2)中得

式(11)表明相位量Φ1與積分變量無關，故對式(12)兩邊乘以exp(-jΦ1)進行相位補償得

對式(13)兩邊進行kθ維傅里葉逆變換得

令式(14)中被積函數的相位為

根據駐定相位法，將式(15)在kθ維進行求導且令其為0，得到kθ=rkrsin(θ-φ)，代入式(14)

其中，kx是方位向波數，ky是距離維波數

式(17)表明波數域數據在直角坐標系下的方位維kx和距離維ky是非均勻分布的[14]，將直角坐標系下非均勻分布的S6(θ,w,kz)通過方位維和距離維的插值，得到均勻分布的頻率數據

其中，interpkx,ky{·}表示插值函數，和分別表示插值后均勻采樣的波數域方位維分量和距離維分量，對式(18)進行3維傅里葉逆變換可得到目標重建圖像

3 同心方形網格插值方法

根據式(16)，給定任一高度維波數，CSAR波數域是一個在極坐標系下表現出為kr非均勻分布、kθ為均勻分布的圓環。通過插值從非均勻極坐標系獲得直角坐標系下均勻分布的波數域值，獲得式(18)所需的呈現同心方形環帶的波數域填充樣式，如圖2所示。盡管如此，式(17)表明波數域數據在方位維和距離維滿足存在強耦合性，傳統插值需要通過2維逐點遍歷實現，即每一個待插值點在2維空間尋找鄰域點，運算復雜度高。因此，本文利用同心方形網格的分布對稱性與CSAR波數域在極坐標下的采樣特點，提出一種波數域數據的2維解耦方法，將2維逐點遍歷插值分解為兩個獨立的1維數據插值處理，大大減小插值所耗的運算量，完成2維插值的高效實現。改進后的插值步驟如圖3所示，分為徑向kr數據補0、徑向插值、平面數據分區、波數域方位向與距離向獨立插值、數據拼接5個步驟來實現。

圖2 柱面SAR成像數據分布圖

步驟1 徑向補0如圖3(a)-圖3(c)所示，得到膨脹的虛擬同心圓環帶。

所需的同心方形帶填充樣式如圖3(a)中實線框定的區域所示，能完全覆蓋原始帶寬引入的同心圓環波數域填充區域(實心環帶所覆蓋的區域)，這里對S5(θ,ω,kz)頻譜數據在徑向進行最內環和最外環的補0膨脹，找到能夠覆蓋同心方形區域的虛擬同心圓環帶。kr_min和kr_max表示虛擬同心圓環帶的內外圓半徑kr_min，由圖3(a)和圖3(b)可以得到kr_min=其中min(kr)和max(kr)為初始圓環的內外圓半徑。根據原始環帶波數域徑向分量kr的采樣規則進行補0，如圖3(c)所示。該步驟保證波數域頻帶利用率最大化，與傳統直接插值方法相比，所提方法考慮到徑向非均勻采樣的特點先通過補0操作確定覆蓋所需同心方形環帶的環帶填充區域可以減少區間的邏輯計算和索引查找等運算量。

步驟2 徑向插值，如圖3(d)-圖3(f)所示，確定所需同心方形帶內波數域徑向的均勻采樣。

為了最大化利用有效頻譜數據，需要將插值前數據即空心點插值到正方實心點上。如圖3(e)所示，由于kr是非均勻的，假設原始徑向kr的采樣個數為NR，先將kr均勻化

圖3(e)為某一原始數據點位置與待插數據點位置的幾何關系，可以得到krv(n)的表達式，徑向插值后數據如圖3(f)所示

步驟3kx-ky平面正交分區 ，如圖3(g)-圖3(h)所示。

經過步驟2，圖3(f)的同心方形帶填充呈現出分維度分塊的均勻性。根據正方形對角線進行2維的正交區間劃分，獲得兩組數據，每組數據包含兩個兩兩不交疊相對的數據塊，分別記為基本垂直點數據和基本水平點數據。圖3(g)為基本垂直點數據Sv，該區間數據在波數距離維ky是均勻分布的，但是波數方位維kx為非均勻；圖3(h)為基本水平點數據SH，該區間數據在波數方位維kx是均勻分布的，但是波數距離維ky為非均勻。根據波數域分布特點，可以對2維插值進行分區降維處理，實現對kx和ky維度的解耦。

步驟4 兩個維度獨立的1維插值，如圖3(i)-圖3(j)所示。

根據分區，將基本垂直點數據Sv沿著方位維利用3次樣條插值法對kx進行內插，基本水平點數據SH沿著距離維利用3次樣條插值法對ky進行內插。注意到在以正方形對角線劃分數據時，對于正方形上半部分數據，從上往下每行數據點的個數依次減少2個，下半部分則反之。所以對于基本垂直點數據Sv的上半部分進行內插時，從上往下每行插值點的個數依次減少2個，下半部分則相反，如圖3(i)所示?？招狞c為插值前的數據，實心點代表插值后的數據。對基本水平點數據SH亦是如此，基本水平點數據SH重采樣后數據如圖3(j)實心點所示。

步驟5 數據拼接，如圖3(k)所示。

將基本垂直點數據Sv和基本水平點數據SH重采樣后的數據點進行拼接得到完整的直角坐標系下分布的均勻數據，拼接后完整的數據如圖3(k)所示，其中實心點為插值后的有效數據，空心圓點為置0數據點。最后將處于直角坐標系下均勻分布的數據進行3維傅里葉逆變換即可得到重建圖像。

4 算法運算量分析

本節對比所提插值方法與逐點遍歷插值方法的運算量，以算法所需的浮點運算量(FLoating point Operations Per Second, FLOPS)來進行衡量[10]。令Nker是1維插值核長度，角度維度采樣點數為Nθ，距離維采樣點數Nf，高度維采樣點數為Nz，插值后均勻網格的方位維和距離維大小分別為Nx和Ny，即Nx=Ny。

2維逐點遍歷(細胞單元法)插值所需運算量為

同心方形網格插值所需運算量為

分析所提插值算法與傳統算法的運算量，定義運算量比γ=C2/C1，將式(22)和式(23)代入可得

Nker一般取8,16或32，這里取Nker=8，且由圖3(i)可知，Nf ＜Nx/2，則

CSAR場景下通常取Nθ ＜4Nx，所以γ ＜1/6?？芍岵逯堤幚矸椒ㄅc傳統插值處理方法運算量比小于1/6，說明所提插值方法運算量小，實際成像中耗時少于傳統插值處理方法，理論上提升的速度最少可達到6倍以上。

5 實驗驗證

本節通過數值仿真和實測數據對比分析本文所提插值方法與傳統2維逐點遍歷插值方法的成像時間與成像效果。實驗參數見表1。成像處理的軟件平臺為MATLAB R2018a，為了更好地對比分析兩種插值方法成像質量和成像時間，成像處理過程中不使用任何濾波處理和并行計算。

表1 實驗參數

5.1 單點目標仿真驗證

仿真實驗中設置的點目標為距離場景圓心0.01 m、角度方向為20°、高度方向為0.01 m的單點目標，則點目標在直角坐標系下坐標為(0.009 4, 0.003 4,0.01)，其中X表示方位維，Y表示距離維，Z表示高度維。圖4給出了傳統2維逐點遍歷插值算法(傳統算法)與本文算法成像的剖面圖。

圖4 點目標成像剖面圖

從3個方向剖面圖可以看出，兩種算法成像結果位置位置基本一致。為了定量評估點目標成像性能，表2給出了點目標在3個維度的分辨率、峰值旁瓣比和積分旁瓣比，其運算量參數對應于表1。由表2可知，改進算法成像3個維度的分辨率與傳統算法成像的分辨率基本一致，各項測量參數非常接近，說明兩種算法點目標成像質量相差不大，并且本文算法插值處理時間只有傳統遍歷插值處理時間的14%，提升速度達到7倍，說明了本文算法大大提高了成像的速度，驗證了算法的高效性。

表2 點目標成像質量參數

5.2 多點目標仿真

以距離場景圓心0.1 m、高度維間隔為0.01 m構建柱形目標點陣進行仿真。圖5給出了本文算法及傳統算法成像3個維度的成像投影圖。從投影圖可以看出，成像位置與仿真設置的位置一致，且設置的點目標均能清晰分辨出，驗證了算法成像的精確性。

圖5 多點目標成像

5.3 實測數據驗證

實驗采用矢量網絡分析儀、喇叭天線、高精度旋轉步進電機和豎直向步進電機搭建實驗平臺如圖6(a)所示。圖6(b)為測試目標鋁箔。

圖6 柱面成像實驗成像場景

圖7給出了目標物品用兩種算法成像的3個維度投影圖，從投影圖可以看出成像結果整體均勻程度和聚焦質量良好，成像目標清晰。對于無參考圖像的質量評價，采用圖像熵和圖像信噪比來說明兩種算法的成像效果[15]。傳統算法3個維度投影圖的圖像熵分別為5.6, 2.2, 2.1，改進算法3個維度投影圖的圖像熵分別為3.6, 2.1, 1.1，由此看出改進算法的圖像熵更小，圖像聚焦效果較傳統算法好；傳統算法3維投影圖的圖像信噪比分別為21.6, 21.3, 19.3，改進算法3維投影圖的圖像信噪比為22.4, 23.8,24.2，通過信噪比可知改進算法的圖像信噪比更大，圖像所含噪聲較少。

圖7 實測數據成像圖

表3給出了兩種算法的運算量和成像時的插值處理時間，可以看出本文算法的運算量與傳統算法相比減少了一個數量級，而且本文算法的插值處理時間約為傳統2維逐點遍歷插值處理時間的14%，提升速度達到7倍，符合算法復雜度理論分析結果。雖然表3中運算量減小倍數10和實際處理時間提速比7不完全相等，是因為在計算算法所需的浮點運算量時，忽略邏輯運算、索引查找等運算。并且在實際實現時，所耗時間還和代碼的編寫風格以及成像處理的軟件平臺MATLAB R2018a的運行模式有關。綜合來說，本文改進的算法插值部分的運算量大幅度減少，插值處理時間短，在保證成像質量的基礎上，大幅提高了成像算法的效率。

表3 理論運算量與插值處理時間

6 結論

本文對CSAR成像算法插值處理進行研究，針對在大角度下方位向與距離向波數域耦合性強，傳統2維逐點遍歷插值處理運算量大和耗時長的問題，提出了基于同心方形網格插值處理的方法。所提方法實現了方位向與距離向波數域解耦合，只涉及兩個獨立的1維插值處理，易于后期硬件實現時的并行化處理。實驗結果表明，在保證成像質量的基礎上，本文算法的插值處理時間約為傳統2維逐點遍歷插值處理時間的14%，提升速度達到7倍。因此，本文提出的基于同心方形網格插值處理的方法能夠有效解決CSAR 3維成像在大角度下方位向與距離向波數域耦合性強、2維插值處理運算量大導致耗時長的問題。