基于交替方向乘子法的Capon層析SAR成像方法

劉 慧，郭馨宇，郭子夜，程碧輝

（北京建筑大學電氣與信息工程學院，北京 102616）

0 引 言

合成孔徑雷達層析成像技術(shù)（TomoSAR）是一種多基線干涉測量技術(shù)，通過多次航過目標，在不同高度形成基線，構(gòu)成第二合成孔徑，反演垂直于視線方向即斜距垂向的目標后向散射系數(shù)［1-2］。通過譜估計的方法可以得到目標散射體斜距垂向的位置，同時還可以得到后向散射體的散射強度信息，實現(xiàn)斜距垂向的合成孔徑成像［3］。

傅里葉變換法是最早被應用于層析SAR 三維成像的方法［4］，通過對同一距離-方位單元的數(shù)據(jù)序列進行傅里葉逆變換來獲取高度向上的后向散射系數(shù)，實現(xiàn)三維成像，但此算法受限于奈奎斯特采樣定理對均勻采樣的要求。根據(jù)瑞利準則［5］，斜距垂向上的分辨率與合成孔徑大小成反比，這就需要保證較多且較長基線，即足夠多的航過次數(shù)。但是多次航過成本高且難以實現(xiàn)，實際中航過間隔甚至超過了奈奎斯特采樣定理的限制，且非均勻基線會導致斜距垂向的低分辨率、高旁瓣問題。不同于冰原、森林等自然地貌，城市場景的三維成像需要高分辨率二維影像的支撐［6］。如何從多次航過且不滿足奈奎斯特采樣定理的高分辨率SAR影像中反演出斜距垂向的后向散射系數(shù)，實現(xiàn)城市場景的三維成像是層析SAR 成像算法研究的關(guān)鍵。

目前，有兩種層析SAR 成像的主流算法。一種是壓縮感知［7-11］（Compressive Sensing,CS）算法，另一種是譜估計算法［12-13］。壓縮感知能夠突破Nyquist 采樣定理的限制，使用稀疏的測量信號獲取信號離散采樣值，實現(xiàn)稀疏信號的重構(gòu)［14］。在層析SAR 成像問題中，待恢復的后向散射系數(shù)在斜距垂向上呈稀疏分布，因此可采用壓縮感知算法實現(xiàn)層析SAR 的三維成像。層析SAR 三維成像的壓縮感知算法主要有凸優(yōu)化算法［8］以及迭代貪婪算法［11］等。CS算法根據(jù)可解析的散射體數(shù)量進行三維重構(gòu)，但是受散射體數(shù)量的限制［7］。

而譜估計算法不考慮散射體的個數(shù)，直接根據(jù)散射體的功率譜強弱判別為目標或噪聲。其中基于波束形成理論的一類算法，從噪聲抑制的角度出發(fā)，通過加權(quán)處理，使信號輸出噪聲方差最小［15］的約束，計算獲得最佳加權(quán)，進而實現(xiàn)對斜距垂向信號的估計。文獻［16-17］采用波束形成的思想，在二維影像較少的情況下，通過譜估計，恢復斜距垂向上的后向散射系數(shù)，從而實現(xiàn)TomoSAR三維成像。此類基于譜估計的非參數(shù)化波束形成方法也可以看作是一種光譜分析信號估計器，非常適合處理分布式散射體的成像問題，它的分辨率很大程度上取決于真實孔徑，即基線的總跨度，能夠改善傅里葉變換成像的低分辨、高旁瓣、成像效果受系統(tǒng)誤差影響大等缺點，具有超分辨能力［18］。

基于波束形成的譜估計算法主要有標準Capon 波束形成算法（Standard Capon Beamformer,SCB）［18］、雙約束魯棒Capon 波束形成算法（Doubly Constrained Robust Capon Beamforming,DCRCB）［19］。波束形成的譜估計算法性能優(yōu)劣主要可以從旁瓣高度、主瓣寬度及穩(wěn)健性三個方面來評價:低旁瓣可以有效抑制干擾，較窄的主瓣寬度可以提高目標的分辨能力，高穩(wěn)健性可以減小各種失配對信號反演結(jié)果造成的影響。文獻［12］采用SCB 算法解決了層析SAR 成像的問題，能夠得到高分辨率的同時還擁有更強的抗干擾能力，但是穩(wěn)健性較差；但在實際場景中，受導向向量和噪聲協(xié)方差矩陣誤差的影響，該方法魯棒性變差，信號處理性能嚴重下降。文獻［19］提出DCRCB 算法，在SCB 算法的基礎(chǔ)上，提高了魯棒性。

本文提出一種改進的波束形成優(yōu)化算法，在DCRCB 的基礎(chǔ)上，結(jié)合L1 范數(shù)的約束函數(shù)，構(gòu)建交替方向乘子法［20］（Alternating Direction Method of Multipliers,ADMM）的代價函數(shù)，將DCRCB 恢復的后向散射系數(shù)進行進一步稀疏優(yōu)化，實現(xiàn)層析SAR 的三維成像。ADMM 算法以增廣拉格朗日算法為基礎(chǔ)，將較為復雜的全局求解問題轉(zhuǎn)換為兩個或多個更易求解的簡單局部子問題。ADMM 算法在迭代中，各子問題可分別完成稀疏重構(gòu)和降噪運算，被分離的局部子問題代數(shù)式都較為簡單，均能較容易地求出確定的解，且不必對其進行收斂運算與約束操作。因此，ADMM 算法具有重建精度高的優(yōu)勢。

本文在第1 節(jié)中構(gòu)建了層析SAR 模型及反演層析SAR成像的波束形成理論；第2節(jié)介紹基于壓縮感知和譜估計方法的層析SAR 成像方法，并提出改進的層析SAR成像算法；第3節(jié)利用仿真參數(shù)和山西運城地區(qū)的8 通道機載陣列干涉SAR 真實數(shù)據(jù)進行實驗驗證，通過實驗對比，驗證了本文提出算法的有效性和優(yōu)勢；第4節(jié)對本文進行總結(jié)。

1 TomoSAR成像原理

1.1 TomoSAR成像模型

層析SAR 成像幾何模型如圖1 所示。圖1 中，x為SAR 平臺的航過方向，即方位向；r為距離向；s為斜距垂向；y為地距向；z垂直于x-y平面。θ為主影像的觀測視角。SAR 平臺經(jīng)過L次航過，得到L幅單視復數(shù)（Single Look Complex,SLC）圖像，經(jīng)過復圖像配準、去斜［21］以及相位誤差補償?shù)忍幚砗螅跋竦耐痪嚯x-方位單元構(gòu)成一個長度為L的序列g(shù)={}g0,g2,…,gL-1，其中g(shù)0表示參考影像，gl(l=1,2,…,L-1)表示輔影像。根據(jù)文獻［22］，沿斜距垂向上的散射點和影像gl(l=0,1,…,L-1)之間的關(guān)系表示為

圖1 層析SAR成像幾何模型

式中，γ(s)為觀測目標沿斜距垂向的連續(xù)域后向散射，Δs為斜距垂向的后向散射信號的分布范圍，為斜距垂向的空間觀測頻率，b⊥l為SAR 平臺第l次航過時與主影像構(gòu)成的垂直基線，λ為入射波長，r為斜距。層析SAR 成像就是已知觀測序列g(shù)={}g0,g2,…,gL-1，反演后向散射系數(shù)γ(s)，從而將斜距垂向疊掩的各散射體分離并獲知其散射強度以及位置的技術(shù)，實現(xiàn)SAR三維成像。

如果沿斜距垂向?qū)⑹剑?）進行離散化，則式（1）可以近似寫為如式（2）所示矩陣形式:

式中，向量gL×1=[g0,g2,…,gL-1]T為對目標的L次觀測數(shù)據(jù)，M是沿斜距垂向的離散化個數(shù)，A是一個L×M階的測量矩陣，n=[n1,n2,…,nL]T為噪聲向量，γM×1=[γ1,γ2,…,γM]T為待估計的后向散射系數(shù)向量。

1.2 基于波束形成估計理論的TomoSAR信號模型

多基線層析SAR 系統(tǒng)的各條航跡沿斜距垂向排列，各個航跡的天線中心可以看作是沿著斜距垂向構(gòu)成的一個相控陣天線陣列，因此上述成像問題可以采用波達方向（Direction-of-Arrival,DOA）估計理論［23］來解決。

由式（2）所表示的線性方程組構(gòu)成的SAR 層析成像信號模型，矢量γ由M個樣本組成，這些樣本是沿斜距垂向位置分別為的散射體的未知連續(xù)復反射率。矢量n為加性噪聲。測量矩陣A由導向矢量組成，每個矢量維度為L，其中包含由位置給定的干涉相位信息，對于每一已知的sm，有

式中，kl=2πξl(l=1,2,…,L-1)，kl(l=1,2,…,L-1)還可以寫為

式中，kl(l=1,2,…,L-1)代表參考衛(wèi)星和第l個輔星之間的雙向垂直波數(shù)。其中，dl(l=1,2,…,L-1)為z方向上參考衛(wèi)星和第l個輔星之間的基線，r是目標到參考衛(wèi)星的斜距，θ代表入射角度，λ是載波波長。

由于維度L通常遠低于樣本數(shù)M，式（2）中的線性方程有無窮多個解，這意味著未知數(shù)個數(shù)M比式（2）中的方程個數(shù)L多。此外，加性噪聲和乘性噪聲的存在增加了求解的不確定性。因此，為了確保對未知矢量γ求解唯一性，必須對式（2）中的測量矩陣A施加一些約束。

通常認為，后向散射矢量γ、噪聲n和觀測向量g為均值為零的復高斯隨機過程，其自相關(guān)矩陣［17］可表示為

式中（·）+代表埃爾米特共軛，{E}是期望算子，因子N0是白噪聲功率的功率譜密度。

在層析SAR 中，通常后向散射矢量γ是不相關(guān)的，因此相關(guān)矩陣Rγ為對角矩陣，則其主對角線上的元素所構(gòu)成的矢量b為后向散射矢量γ的二階統(tǒng)計量，即功率譜，該功率譜可用來表征層析SAR 三維成像斜距垂向上的后向散射。為求解出后向散射功率譜矢量b，可基于DOA估計理論的波束形成算法，利用導向矩陣A以及信號和噪聲統(tǒng)計的先驗信息［24］解決該問題。

2 TomoSAR三維成像方法

2.1 壓縮感知層析SAR成像方法

基于壓縮感知的層析SAR 成像主要解決的數(shù)學問題如式（8）所示:

式中AL×M為測量矩陣，由基線、斜距、波長等已知參數(shù)確定，觀測所得的SAR 影像為gL×1，γM×1為待重建的稀疏信號，即斜距垂向的后向散射。

本文采用壓縮感知中的貪婪算法開展層析SAR 三維成像的研究。正交匹配追蹤算法［11］（Orthogonal Matching Pursuit,OMP）是貪婪算法中的經(jīng)典算法。該算法采用迭代的方法得到待恢復的稀疏信號。算法初始化γ(0)=0，在第n次迭代時，用測量矩陣AL×M與殘差r(n)=g-A·γ(n)進行施密特正交，然后找出前K（K表示后向散射的稀疏度）個最大元素對應的索引值添加到索引集Λ(n)。設(shè)第n次迭代滿足結(jié)束條件，則恢復稀疏向量如式（9）所示:

OMP算法的偽代碼迭代流程如表1所示。

表1 正交匹配追蹤算法（OMP）偽代碼

2.2 基于Capon波束形成譜估計層析SAR成像方法

譜估計算法中的波束形成技術(shù)獲取層析SAR成像所需的后向散射系數(shù)，從噪聲抑制的角度出發(fā)，通過加權(quán)處理，使信號輸出噪聲方差最小E[σ2]=h+mRghm的約束，計算獲得最佳加權(quán)，進而實現(xiàn)對斜距垂向信號的估計。其設(shè)計原理是讓濾波器對感興趣方位的信號無失真地輸出。本文基于式（7）研究了Capon 波束形成譜估計的算法，其最佳加權(quán)的約束方程為

該約束方程可由拉格朗日因子法求解，有

式（11）中的自相關(guān)矩陣Rg可由樣本協(xié)方差矩陣代替，J代表觀測數(shù)目。濾波器的輸出功率構(gòu)成的對角線矩陣為后向散射γ的自相關(guān)矩陣。由此可得出式（5）的Capon波束形成譜估計結(jié)果，具體的可寫為

由于Capon濾波器要求數(shù)據(jù)協(xié)方差矩陣可逆，對角加載法是針對式（7）所示數(shù)據(jù)協(xié)方差矩陣Rg的秩不足時，提高Capon魯棒性而普遍采用的一種方法［13］。一般來說，理想模型化的數(shù)據(jù)自相關(guān)矩陣Rg是不等于樣本數(shù)據(jù)協(xié)方差矩陣G的。由于理想信號的波達方向和真實信號的波達方向之間存在的差異，當預期響應和實際響應之間存在差異時，Rγ和對角加載矩陣可以解釋為轉(zhuǎn)向矢量誤差造成的差異［24］。

式中U和Γ通過特征分解G=UΓU+得到，U的列包含G的特征向量，對角矩陣Γ的主對角線聚集了各向量的特征值，因子對應于從下式獲得的拉格朗日常數(shù):

2.3 基于DCRCB譜估計的交替方向乘子法

對于如下的分布式凸優(yōu)化問題［20］:

式中，f1(x)為損失函數(shù)，f2(y)為正則項，C，D，E為常系數(shù)矢量。對于對偶變量(x,y)，由于其由各自獨立的函數(shù)約束，因此可以采用分解協(xié)調(diào)的形式，將小局部子問題的解通過協(xié)調(diào)以找到大局或全局問題的解。交替乘子法正是基于這樣的思想而求解上述凸優(yōu)化問題的。交替乘子法以增廣的拉格朗日函數(shù)作為基礎(chǔ)，融合雙重分解和增廣拉格朗日方法的優(yōu)點，用于約束優(yōu)化。具體的對于上述凸優(yōu)化問題，交替乘子法構(gòu)建如下的代價函數(shù):

式中，μ為拉格朗日乘子向量，向量維度取決于線性方程組Cx+Dy=E的維度，c為交替乘子法的正則化因子。

論文引入交替方向乘子法解決Capon 譜估計的層析SAR 成像最優(yōu)化問題。層析SAR 成像的分布式凸優(yōu)化問題構(gòu)建如下:

式中，ρ是正則化因子，通常取值為0.5～1.0，引入Capon 譜估計的層析SAR 成像結(jié)果作為中間變量，=(m=1,2,…,M)，由式（18）可以構(gòu)建對于層析SAR 成像問題的交替方向乘子法（ADMM）優(yōu)化函數(shù):

通過迭代的方法即可估計出最優(yōu)的成像結(jié)果。

2.4 交替方向乘子法的處理流程

交替方向乘子法通過交替迭代運算實現(xiàn)對式（18）所表示的優(yōu)化問題的求解。具體如下:

式中，γk+1表示第k+1 次迭代時對γ估計結(jié)果，表示第k+1 次迭代時對重新估計的結(jié)果，重新估計可以有效抑制噪聲，μk+1表示第k+1 次迭代時拉格朗日乘子向量更新，c為常數(shù)，通常可在0.03～0.05 范圍內(nèi)取值。γk+1與可通過對式（18）求偏導并將偏導置零來求解得到上述運算的解析表達式如下:

式中，I表示單位矩陣，S(·)(·)表示軟門限函數(shù)，對于任意向量e與實標量門限ζ，有

式中，ui=Sζ(ei)，其中ei和ui分別表示向量e和u的第i個元素。

交替方向乘子算法通過式（20）進行迭代更新，當滿足迭代結(jié)束條件

則輸出=γk+1，其中ε為極小數(shù)為稀疏優(yōu)化后的后向散射系數(shù)。具體的交替方向乘子算法流程如圖2 所示。圖2 中，首先進行初始化，然后通過式（20）進行迭代更新，最后輸出最優(yōu)化后向散射系數(shù)。

圖2 交替方向乘子算法流程圖

2.5 算法復雜度

壓縮感知OMP 算法包括選擇和更新兩個環(huán)節(jié)。首先是索引集Jn的選擇，即表1迭代過程的第一步，涉及排序算法，時間復雜度為O(KLM)；K表示稀疏度，L、M分別表示測量矩陣的維度。壓縮感知OMP 算法第n次更新迭代的環(huán)節(jié)涉及到最小二乘法求解問題，時間復雜度為O(nN)。OMP 最多在2K次迭代終止，因此，OMP 總運行復雜度為O(KLM)。

本文提出的基于DCRCB 譜估計的交替方向乘子算法，在迭代過程中的運算量主要集中在更新估計向量γ時的矩陣求逆過程。從式（20）可以看出，當A和c固定時，待求逆的矩陣為一常數(shù)矩陣，這使得求逆運算可從循環(huán)中分離出來且僅為1次，其運算量為O(LM2)次浮點操作。可見兩種算法的運算復雜度都與測量矩陣維度有關(guān)，在成像過程中，受基線和瑞利分辨率的限制，M的取值過大時，對成像結(jié)果的影響不大，因此兩種算法的計算復雜度可視為同一量級。表2 列出涉及兩種算法的時間復雜度。

表2 兩種算法的時間復雜度

3 實驗結(jié)果

為驗證提出算法的成像結(jié)果，下面分別使用仿真數(shù)據(jù)和真實數(shù)據(jù)進行實驗，其中真實數(shù)據(jù)采用中國科學院空天信息創(chuàng)新研究院在《雷達學報》官方網(wǎng)站公開的機載陣列干涉SAR 系統(tǒng)獲取的山西運城地區(qū)的8通道數(shù)據(jù)［25］，進行大范圍城區(qū)建筑物層析SAR 三維重建。實驗環(huán)境的參數(shù)如表3所示。

表3 機載陣列干涉SAR系統(tǒng)參數(shù)

3.1 層析SAR仿真實驗

利用表2所示的機載陣列干涉SAR系統(tǒng)參數(shù)，我們仿真了如圖3 所示的非均勻基線分布時4 種算法的成像，結(jié)果如圖4 所示。從圖4（a）、（b）和（c）中我們不難看出，壓縮感知算法和以Capon、DCRCB 為代表的波束形成算法在非均勻基線成像時對建筑物的頂、底以及墻體各面的還原效果都不佳，點云間斷且有丟失現(xiàn)象，兩種波束形成算法還有大量的噪點存在；而圖4（d）所示的本文算法對建筑物形狀具有較強的還原能力，點云形狀完整且噪聲抑制能力更優(yōu)，同時成像結(jié)果更加清晰，這說明本文算法的層析SAR 成像方法在真實數(shù)據(jù)的應用場景中更加具有優(yōu)勢。

圖3 非均勻基線分布

圖4 仿真實驗成像結(jié)果

3.2 層析SAR三維成像方法驗證

運城地區(qū)的實驗場景的Google Earth光學影像和SAR影像分別如圖5（a）和（b）所示。

圖5 實驗場景的SAR圖像及其對應的光學影像

為了驗證算法，論文針對圖5中紅線所示的方位向，即影像的第2 000 行進行了層析SAR 三維成像反演。圖6 給出了分別采用CS 算法、Capon 算法、DCRCB 算法和本文提出算法的成像結(jié)果，并對實驗結(jié)果對比分析。

圖6 運城地區(qū)真實數(shù)據(jù)SAR層析成像結(jié)果

從圖6 可以看出，CS 算法將目標看作一組稀疏向量，需要提前設(shè)定稀疏度進行成像，但在實際場景中，強噪聲可能也會替代目標，導致目標散射體丟失；也可能由于噪聲的影像，會導致散射體的錯位或偏移。Capon 算法、DCRCB 算法和本文算法是波束形成算法，從信號處理的角度出發(fā)，只單純對所有散射點的回波信號進行濾波，本質(zhì)上是從成像角度解決層析SAR 第三維的成像反演問題，并不會和CS 算法一樣依賴于稀疏度。Capon算法由于其測量矩陣A所具有的正交性，功率譜也呈現(xiàn)較尖銳的形狀，能夠?qū)⑴园暌种频煤艿停敯粜暂^差，在實際場景中極易收到干擾導致誤差和噪聲嚴重。DCRCB 算法在Capon 算法的基礎(chǔ)上增強了魯棒性，但同時就會帶來分辨率的損失。本文提出的算法則進一步提高分辨率的同時，并降低旁瓣，可提升層析SAR 成像結(jié)果的分辨率及抗干擾能力。

圖7 給出了第2 000 行第460 列的分辨單元的成像結(jié)果。可以看出，CS 算法和Capon 算法功率譜呈現(xiàn)較尖銳的形狀，魯棒性較差，但分辨率較高。而DCRCB 算法和本文算法增強了魯棒性，本文算法相對DCRCB 算法降低旁瓣的同時更好地提高了分辨率。

圖7 某一距離-方位單元成像結(jié)果（第2 000行第460列）功率強度

3.3 山西運城層析SAR三維成像結(jié)果

圖8 展示了基于上述4 種算法對山西運城某地區(qū)的層析SAR 三維重建結(jié)果，其中，圖8（a）和（b）所示的CS 算法和Capon 算法噪點過多，建筑立面分辨情況不佳，圖8（c）所示的DCRCB 算法分辨率有所提升，但噪點過多現(xiàn)象依然存在，最后本文提出的基于交替方向乘子法的Capon 譜估計算法如圖8（d）所示，在減少噪點的同時最大程度重建樓體輪廓，在不采用任何降噪手段處理時也能得到清晰成像結(jié)果。

圖8 運城地區(qū)實驗場景層析SAR三維重建結(jié)果

為了進一步解決噪聲問題，本文采用DBSCAN方法提取有效點云［26］，利用點云的密度信息，提取其中密度較大的點云聚類。圖9 給出了本文所提及的4 種算法經(jīng)過DBSCAN 提取點云之后對整個建筑群的三維重建結(jié)果。從圖9的處理結(jié)果來看，本文提出的基于交替方向乘子法的Capon 估計算法對后續(xù)的去噪方法沒有強依賴關(guān)系，而其他3種成像方法必須經(jīng)過去噪處理才能獲得建筑的清晰成像結(jié)果。

圖9 經(jīng)DBSCAN降噪后的運城地區(qū)實驗場景層析SAR三維重建結(jié)果

4 結(jié)束語

本文提出了一種基于交替方向乘子法的Capon 譜估計算法進行層析SAR 三維成像，并大規(guī)模重建了山西運城某城區(qū)的建筑物群，同時在噪點情況、分辨率、主瓣以及旁瓣等方面，與現(xiàn)有算法包括CS 算法、Capon 譜估計算法以及DCRCB 譜估計算法進行了SAR 層析成像結(jié)果的對比分析。通過實驗，表明了基于傳統(tǒng)CS 算法的層析SAR 成像在散射點選擇問題上的一些弊端，而本文中提及的3 種基于波束形成的層析SAR 成像算法則從信號處理的角度，對所有散射點進行全面的反演工作，但這類算法需要目標信號的導向向量及噪聲干擾的協(xié)方差矩陣，這一限制使其在真實場景中的性能大幅度降低，本文在對導向向量加以約束來提高其魯棒性的基礎(chǔ)上，通過L1范數(shù)，進一步提升此類算法的分辨率，減少噪聲，獲得更優(yōu)成像結(jié)果，目的是突出信號處理算法的強還原能力。