基于變分模型的陣列三維SAR最優DEM重建方法

師 君張曉玲 韋順軍 向 高 楊建宇

(電子科技大學電子工程學院 成都 611731)

1 引言

陣列3維SAR[1－8]是指利用線型天線陣列代替傳統合成孔徑雷達的單天線結構，通過控制天線陣列的運動，合成2維等效陣列，并結合脈沖壓縮技術獲得3維空間分辨能力，以實現對觀測目標3維成像的技術。與傳統SAR相比，陣列3維SAR具有下視3維成像能力，在地形測繪、災害監測等領域具有廣泛的應用前景。

但是，由于陣列3維SAR需要在機翼上布設線陣天線以獲得該方向的分辨率，載機平臺尺寸限制使得其陣列方向分辨率遠遠低于距離向和航跡向。該缺陷嚴重制約了陣列3維SAR系統整體性能的提升。因此，需要針對陣列3維SAR數據特點，從成像/圖像處理的角度研究提高陣列向分辨率的方法。

由于3維SAR圖像顯著的稀疏特性[9]，稀疏重建技術，如L1正則化方法，OMP，基追蹤等技術目前已經被廣泛用于層析SAR、圓周SAR、軌道3維SAR的成像處理，在提高圖像分辨率、消除旁瓣等方面取得了一定的效果。但在深入研究后發現，直接將稀疏重建技術用于陣列3維SAR圖像增強問題仍存在一些問題。

首先，稀疏特性并不能完全描述陣列3維SAR圖像所具有的典型特征。在地形測繪任務中，目標散射點 3維空間中沿著某個特定的曲面(數字高程圖/DEM)分布，圖像增強應直接針對DEM圖進行處理。而在采用稀疏重建技術時，DEM 圖特有的“單值性”、“連續性”等特征被重述為“稀疏性”，實際上導致了約束條件的放松。單值性的缺失可能導致構造的高程圖存在一個x,y位置對應多個高度的不合理情況，最終導致DEM提取誤差。

其次，由于陣列3維SAR圖像為3維數據，采用稀疏重建技術可能需要構造上億維的線性方程組，并進行最優化求解，大大增加了算法的實現難度和可行性。而降維處理策略則會導致額外的模型誤差，進一步影響DEM增強效果。

針對上述問題，本文討論了基于變分模型的陣列3維SAR最優DEM重建方法，該方法直接將DEM 圖作為最優化目標，通過尋找最優化 DEM圖和對應的散射系數，實現最小二乘意義下的最優DEM重建。第2節簡要回顧了陣列3維SAR的原理和成像處理方法，第3節對陣列3維SAR最優DEM重建的變分模型進行了分析討論。第4節討論了變分模型的求解方法，并提出了一種基于變分模型的陣列3維SAR最優DEM重建方法。第5節對該方法的性能進行了分析比較，驗證了該方法的有效性。

2 陣列3維SAR原理

陣列3維SAR系統在機翼上安裝線陣天線，通過陣列天線的相對運動合成2維虛擬面陣，以獲得2維空間分辨率，并結合SAR雷達的距離分辨率實現對觀測區域的3維成像，如圖1所示。

圖1 陣列3維SAR工作幾何Fig.1 Geometry of LASAR

陣列3維SAR的高度向分辨率由脈沖壓縮技術獲得，根據脈沖壓縮理論，其分辨率為c/(2B)，其中，c為光速，B為發射信號帶寬。陣列3維SAR的沿航跡向分辨率通過平臺運動構成的合成孔徑獲得，根據合成孔徑理論，其分辨率為λ/(2θ)，其中，λ為波長，θ為合成孔徑角度。陣列 3維 SAR的切航跡向分辨率通過在機翼上布設陣列天線獲得，其分辨率為 λR/(κA)，其中，R為雷達到目標的距離，A 為陣列長度，κ為常數，1＜κ＜2。假設雷達工作于Ka波段，陣列長度10 m，飛行高度10 km，則可得到系統陣列向分辨率約為4～8 m，遠遠低于距離向和航跡向分辨率。

在成像處理算法方面，由于陣列3維SAR的等效2維天線陣列結構復雜，RD算法、ωK算法適用性較差，無法滿足任意模式3維SAR成像處理要求。因此，陣列3維SAR成像處理主要采用后向投影算法等時域處理技術。

采用后向投影算法進行成像處理后得到的3維SAR圖像如圖2所示，其中，x方向為陣列方向，y方向為平臺運動方向，z方向為高度向。可以看出，其能反映出目標的典型幾何特征。但是，觀察其單個高度切片的圖像可以看出，其圖像在陣列方向分辨率較低，且由于稀疏陣列的采用，導致圖像像素間的串擾較大。

3 最優DEM重建的變分模型

3.1 變分模型

在地形測繪任務中，陣列3維SAR主要觀測目標為起伏地形，可寫成散射系數與高程相對于水平面坐標系的函數的形式，即：

其中，x,y表示 2維水平面的變量，z(x,y)反映了x,y處地形的高度值，σ(x,y)反映了x,y處的復值散射系數。

假設陣列 3維 SAR系統的 3維模糊函數為χ(x,y,z)，則采用后向投影算法得到的成像結果可由高程函數z(x,y)和散射系數函數σ(x,y)表示為：

圖2 陣列3維SAR成像處理結果Fig.2 Imaging results of LASAR

其中，u,v表示2維水平面的變量，I(x,y,z)表示利用系統模糊函數和散射系數估計得到的 3維 SAR圖像。

利用式(2)和實際系統成像結果，可建立式(3)所示最優化問題模型：

其中，e(x,y,z)表示實際 3維圖像與估計圖像間的誤差，表示實際系統得到的3維SAR圖像。

式(3)表明，陣列3維SAR圖像增強問題的可描述為針對z(x,y)和σ(x,y)的最優化問題：尋找高程函數z(x,y)和散射系數函數σ(x,y)，使得該高程函數和散射系數函數對應的估計結果與實際測量結果的誤差最小。由于該最優化問題的優化變量為關于自由變量x和y的函數，其在數學上屬于典型的變分問題。

為了便于分析與算法設計，需要進一步將3維模糊函數 χ(x,y,z)寫成3個1維模糊函數乘積的形式，即：

其中，χR(?),χL(?)和χA(?)分別為距離向、陣列方向和航跡向的模糊函數。一般情況下，距離向和航跡向的系統分辨率優于或等于3維圖像像素間隔，χR(?)和χA(?)可近似看作沖激函數，則式(2)可近似表示為：

此時，式(3)的最優化問題可沿航跡向劃分為一組相互獨立的最優化問題：

其中，“/y”表示“給定航跡向距離單元y”，χ(x)=χL(x)。式(6)將二元變分優化問題轉化為一組一元變分問題，大大降低了算法實現難度。

觀察式(6)可以發現，一方面，由于最優化變量z(x)與模糊函數復合，該最優化問題屬于典型的非線性變分，采用歐拉-拉格朗日方程法求解較為困難。另一方面，在z(x)已知條件下，關于散射系數σ(x,y)的最優化問題屬于典型的最小二乘問題，可以采用偽逆、稀疏重建等方法進行求解。

因此，針對式(6)的最優化問題可分成兩個階段：首先，給定z(x)，求解最小二乘問題，得到對應的最優化散射系數和最優化代價函數值；其次，搜索所有可能z(x)，比較最優化代價函數值，選擇最小的一個作為最優化z(x)其對應的最優化散射系數作為期望的散射系數。

3.2 局部化近似

為了降低求解該變分問題的難度，首先需要對其進行局部化近似處理，其核心是忽略陣列向模糊函數的旁瓣(在超分辨問題中，像素間隔遠遠小于陣列向分辨率，故散射點間的旁瓣串擾影響遠遠小于臨近區域散射點間的主瓣串擾)。在此近似條件下，某個x位置處高度z和其對應的散射系數的改變只對其臨近區域的數據產生影響，如圖3所示。

圖3 DEM增強問題局部化Fig.3 Localization of the DEM enhancement problem

假設主瓣寬度為 2L+1，當前狀態為 x，則改變當前位置的值將導致x?L到x+L段數據值(包含所有高度)發生改變。相應地，x?L到x+L數據段值的改變則會對x?2L到x+2L段散射系數產生影響。基于此特征，滿足最優化問題式(6)的x處最優高度z，同樣應滿足式(7)所示最優化問題：

式(7)的物理意義為：尋找以x為中心的4L+1點最優路徑使得以x為中心的2L+1點數據對應的均方誤差最小。通過進一步分析可以發現，在無誤差條件下，式(7)局部化問題的最優路徑同樣也是最優化問題式(6)的最優解。

為了便于計算最小二乘問題，需要通過數據重排技術，進一步將式(7)表示為矩陣形式，該過程可以通過簡單的代數知識解決，此處不再贅述。當L=2時，給定高度h其對應的矩陣公式為：

需要注意的是，此處的散射系數σ(u,h)為 x-z切片內的散射系數，其為x和高度h的函數，而非式(1)中x與y的函數。此時，σ(u,h)具有稀疏性，即 ：σ(x,h)=0,h ≠ z(x)。 一 般 情 況 下 ，φ 為(2L+1)×(4L+1)矩陣，σ[h]為 4L+1 點列向量，d[h]為2L+1 點列向量。

對于整個2維相關數據段區域，其對應的矩陣表示為：

其中，Nz為3維圖像沿z方向的像素點數。需要注意的是，根據其物理意義，σ只包含 4L+1 個非零分量。另外，如果事先已知σ的某些分量為0，則可以去除式(9b)中的對應行，以簡化最小二乘問題。

4 基于OMP算法的變分模型求解

為了求解式(6)的全局最優化問題，需要消減其可行解集以達到降低運算量的目的。本文方法的基本策略是將全局最優化問題按照陣列方向分解為一系列的局部最優化問題，每個階段只決定一個x位置對應的最優高度，并通過滑窗迭代的結構實現對整個x-z切片的最優化搜索。

通過觀察可以發現，由于滑窗結構的采用，當優化x位置時，其先前位置的最優化高程已經獲得，而無需進行重復搜索。另一方面，x位置的后續位置只用于估計當前局部最優化問題的最小均方誤差，其各位置的高程將在后續迭代中計算，因此可以放松后續位置的“單值映射條件”。

基于上述考慮，式(7)可轉化為改進的稀疏重建問題，即：

其中，第1個約束條件表示位置x?2L到x?1對應的高程為先前獲得的最優化路徑； Θ[x]為一個較大的復數，表示x狀態處存在一個“顯著的”散射系數；Sp(σ)表示σ的稀疏度，第3個約束條件表示整個優化問題的稀疏度為4L+1。

該最優化問題可以通過對OMP算法進行改進求解。第1個約束條件表明在進行OMP算法迭代初始，需要將x?2L 到x?1位置處已知高度對應于矩陣Ψ的指標集作為初始支撐集傳遞給OMP算法；第3個約束條件表明，OMP算法的迭代次數應為2L次(前 2L個高程已知，x處對應的高程需要遍歷，在單次OMP調用時也相當于已知)。由于Θ[x]未知，且其取值隨著x的變化而變化，第2個約束條件處理起來較為困難，本文采用式(11)估計式(10)的最優代價函數(詳見附錄)。

至此，基于OMP算法的變分模型求解算法如下：

目標：求解式(7)定義的局部最優化問題。

輸入：給定矩陣Ψ，先前高程z(u)，觀測數據d,L和當前狀態x。

步驟1 高程遍歷循環：遍歷當前位置的所有可能高程，并進行如下操作：

(1) 選擇z(x)，先前高程以及所有后續位置所對應的矩陣Ψ的行向量，獲得Ψ的子矩陣；

(2) 記錄z(x)和先前高程所對應的Ψ的子矩陣的行序號，作為OMP算法的初始支撐集；

(3) 調用OMP算法；

(4) 記錄OMP算法的殘留誤差和z(x)對應的散射系數σ(z)。

步驟 2 尋找σ(z)的最大值，并利用式(11)計算J～[x](z)。

步驟 3 尋找J～[x](z)的最小值對應的位置作為當前位置x的最優高程。

輸出：當前位置x的最優高程。

通過遍歷所有沿陣列方向的位置x，即可得到當前x-z切片內的最優高程函數。

根據上述分析，基于OMP算法的變分模型求解方法的運算量為：

5 性能分析

由于超分辨率條件下圖像增強問題屬于典型的病態問題，其性能與信噪比密切相關。當無噪聲條件下，很多算法都能獲得較好的增強效果；但隨著噪聲的增加，不同算法的性能差異很大。本節仿真實驗在固定超分辨倍數的條件下，分析不同算法的重建誤差，并對基于變分模型的DEM重建算法的噪聲性能進行分析比較。

仿真地形包括起伏山區地形和城市區域兩種，圖像尺寸為150像素×200像素。根據前面的分析可知，由于不同沿航跡單元內距離-陣列切片分的問題相互獨立。本仿真也可以看作在相同信噪比條件下進行的不同地形200次蒙特卡洛實驗。為了便于對比，標準的OMP算法[10－12]和LASSO算法[13－15]也被用于陣列3維SAR圖像增強，其基本過程是抽取3維SAR圖像中每個航跡向和距離向柵格的陣列向數據，調用稀疏重建方法得到處理結果，然后對選擇每個水平面柵格節點處最大值作為該處高程的估計。

假設L=4，散射系數幅度為1(隨機相位)，噪聲標準差為0.1，變分模型方法、OMP和LASSO算法對起伏地形的處理結果如圖4所示。對比可以看出，采用變分模型的處理方法可以在較大噪聲條件下實現起伏地形的精確 DEM 增強，其對應的DEM殘差圖如圖4(d)所示，其中只包含155個誤差點。與之相比，OMP算法和LASSO算法所得到的處理結果則存在較大的噪聲，其誤差點分別為12691個和5926個。

圖4 起伏地形的最優DEM重建結果(L = 4，噪聲標準差0.1)Fig.4 Enhancement results of a mountain area with L = 4 and the standard deviation 0.1

圖5為仿真城市區域的DEM增強結果。觀察圖5(a)所示的觀測地形可以發現，由于城市區域房頂等建筑結構與地面平行，其沿陣列方向抽取的數據表現出了“連續性”而非“稀疏性”。此時采用稀疏重建模型的算法處理將更為困難。圖5(b)為變分模型的處理結果，可以看出其能夠較好地重建城市區域地形。與之相反，圖5(e)和圖 5(f)的OMP算法、LASSO算法得到的結果則由于稀疏性的破壞導致處理效果嚴重下降。通過進一步統計發現，變分模型法、OMP和LASSO算法得到的DEM圖像的殘留誤差分別為119,24473和14940。

通過上述分析比較可以看出，由于利用了DEM圖的單值性特征，基于變分模型的DEM重建方法性能要遠高于基于稀疏信號模型的重建方法，且在處理起伏地形、城市區域時都能得到較好的處理結果。

6 結論

本文提出了基于變分模型的陣列3維SAR最優DEM重建方法，該方法直接將DEM圖作為最優化目標，通過尋找最優化DEM圖和對應的散射系數，實現最小二乘意義下的最優DEM重建。仿真分析表明，該方法可以實現各種地形(山區、城市)等的穩健DEM增強，其性能遠優于OMP算法和正則化方法。

圖5 城市區域的最優DEM重建結果(L = 4，噪聲標準差0.1)Fig.5 Enhancement results of an urban area with L = 4 and the standard deviation 0.1

附錄

本附錄將推導式(11)，該問題可以重述為：

給定線性方程組Ax=y和其OMP算法解*x，計算該線性方程組在約束條件xi=c下的殘留誤差，其中，i屬于支撐集。

為了分析上述問題，可將 OMP算法的殘留誤差寫為：

另一方面，Ax=y且xi=c可以寫為：

其中，A↓i表示剔除第i列后A矩陣的子矩陣，x↓i表示剔除第i個分量后x向量的子向量，αi表示矩陣A的第i列。

由于i屬于支撐集，式(A-3)可表示為：

根據OMP算法的正交性，則有：

即式(11)。

[1]Matthias W and Markus G.Initial ARTINO radar experiments[C].20108th European Conference on Synthetic Aperture Radar (EUSAR),Aachen,Germany,2010: 1-4.

[2]Han Kuo-ye,Wang Yan-ping,Tan Wei-xian,et al..Efficient pseudopolarformat algorithm for down-looking linear-array SAR 3-D imaging[J].IEEE Geoscience and Remote Sensing Letters,2015,12(3): 572-576.

[3]Peng Xue-ming,Hong Wen,Wang Yan-ping,et al..Downward looking linear array 3D SAR sparse imaging with wave-front curvature compensation[C]. International Conference on Signal Processing,Communication and Computing (ICSPCC),Kunming,2013: 1-4.

[4]Zhang S,Zhu Y,and Kuang G.Imaging of downward-looking linear array three-dimensional SAR based on FFT-MUSIC[J].IEEE Geoscience and Remote Sensing Letters,2015,12(4):885-889.

[5]Peng Xue-ming,Hong Wen,Wang Yan-ping,et al..Polar format imaging algorithm with wave-front curvature phase error compensation for airborne DLSLA three-dimensional SAR[J].IEEE Geoscience and Remote Sensing Letters,2014,11(6): 1-5.

[6]Wei Shun-jun,Zhang Xiao-ling,and Shi Jun.Sparse autofocus via bayesian learning iterative maximum and applied for LASAR 3-D imaging[C].2014 IEEE Radar Conference,Cincinnat,USA,2014: 666-669.

[7]Zhang Ying-jie,Han Kuo-ye,Wang Yan-ping,et al..Study on motion compensation for airborne forward looking array SAR by time division multiplexing receiving[C].2013 Asia-Pacific Synthetic Aperture Radar (APSAR),Tsukuba,Japan,2013:392-395.

[8]Zhuge Xiao-dong and Yarovoy A G.A sparse aperture MIMO-SAR-Based UWB imaging system for concealed weapon detection[J].IEEE Transactions on Geoscience and Remote Sensing,2011,49(1): 509-518.

[9]Shi Jun,Zhang Xiao-ling,Xiang Gao,et al..Signal processing for microwave array imaging: TDC and sparse recovery[J].IEEE Transactions on Geoscience and Remote Sensing,2012,50(11): 4584-4598.

[10]Needell D and Tropp J A.CoSaMP: iterative signal recovery from incomplete and inaccurate samples[J].Applied and Computational Harmonic Analysis,2009,26(3): 301-321.

[11]Wang Jian,Kwon S,and Shim B.Generalized orthogonal matching pursuit[J].IEEE Transactions on Signal Processing,2012,60(12): 6202-6216.

[12]Wang Jian and Shim B.On the recovery limit of sparse signals using orthogonal matching pursuit[J].IEEE Transactions on Signal Processing,2012,60(9): 4973-4976.

[13]Tibshirani R.Regression shrinkage and selection via the lasso[J].Journal Of The Royal Statistical Society Series B-Methodological,1996,58(1): 267-288.

[14]EfronB,Hastie T,Johnstone I,et al..Least angle regression[J].Annals of Statistics,2004,32(2): 407-499.

[15]Rosset S and Zhu Ji.Piecewise linear regularized solution paths[J].Annals of Statistics,2007,35(3): 1012-1030.