可變先驗貝葉斯學習稀疏SAR成像

沈笑云, 廖仙華, 孫衛天, 夏亞波, 楊 磊

(1. 中國民航大學電子信息與自動化學院, 天津 300300;2. 中國民航大學天津市智能信號與圖像處理重點實驗室, 天津 300300)

0 引 言

合成孔徑雷達(synthetic aperture radar,SAR)作為一種軍事和遙感中廣泛使用的微波探測手段,具有全天時、全天候優點。在SAR成像中,經典的距離多普勒(range Doppler,RD)算法簡單高效,但抑噪能力差,無法對特定目標進行特征學習成像[1]。而如何實現高分辨SAR特征成像一直是業界研究熱點。近年來,壓縮感知(compressed sensing, CS)[2]理論的提出,為SAR特征學習成像提供了理論和技術支持[2-3]。CS充分利用雷達成像目標數據特征,只需少量采樣,就能實現原信號的無失真稀疏重建。

在CS框架下,目前常用的兩種高分辨SAR特征學習成像算法分別為凸優化類算法和貝葉斯學習類算法[4]。凸優化類算法通過正則先驗對特征進行建模[5]。其典型算法包括基追蹤算法、稀疏梯度投影算法、快速閾值迭代算法、交替方向乘子法[6-8]等。這類算法主要優點是運算效率高。但隨著成像場景復雜性和特征多樣性的增加,傳統單正則項已然無法滿足高分辨成像要求,基于多特征的多正則項協同增強成像已成為發展趨勢。然而隨著正則項的增多,正則項系數地選取成為一個高維空間最優參數搜索問題,因此導致算法效率極大降低。

貝葉斯學習類成像算法則是近幾年來隨著統計采樣技術發展起來的高分辨SAR成像算法[9-10]。該類算法采用貝葉斯統計學習框架,摒棄凸優化類算法繁瑣的正則參數調整步驟,并且表現出更強大的建模靈活性與參數自學習特性。目前,貝葉斯學習成像算法主要包括變分貝葉斯-期望最大化(variational Bayes expectation maximization, VB-EM)算法、稀疏貝葉斯學習(sparse Bayesian learning,SBL)算法等[11-13]。2016年,Yang等人將貝葉斯統計學習用于SAR動目標成像,并使用VB-EM算法進行稀疏信號求解,獲得較好的成像效果[14]。文獻[15]基于目標散射點的連續性特征進行稀疏逆SAR(inverse SAR, ISAR)成像恢復,獲得了較高精度的ISAR成像結果。但傳統貝葉斯學習類算法目前主要基于目標稀疏特征或連續性結構特征進行先驗建模。先驗分布通常為拉普拉斯分布或高斯分布,先驗建模固化,靈活性較差,難以適應復雜動態變化的目標特征,且采樣求解過程面臨高維矩陣求逆和隨機游走采樣效率低下的問題。

針對當前貝葉斯學習SAR成像算法中存在的問題,提出一種可變先驗貝葉斯(variable imaging prior Bayes,VIP-Bayes)學習高分辨成像算法。該算法首先將目標特征建模為一種VIP先驗,即廣義高斯分布(generalized Gaussian distribution,GGD)先驗[16-17]。該先驗通過其形狀參數的靈活可變性,可模擬適應不同的SAR目標場景成像特征,因此能很好地解決先驗固化問題。其形狀參數的靈活多變性賦予該先驗對不同SAR目標特征進行表征的可能性[18]。然后,在先驗建模基礎上,進行分層貝葉斯模型,構建似然函數與GGD先驗參數間的聯系,進而推導成像目標后驗分布。因所得后驗分布十分復雜,含有lp范數指數項,無法使用常規馬爾可夫鏈蒙特卡羅(Markov chain Monte Carlo,MCMC)采樣算法如吉布斯[19]采樣進行求解。本文引入一種哈密頓蒙特卡羅(Hamiltonian Monte Carlo,HMC)采樣算法[20]進行求解。HMC算法基于哈密頓動力學勢能與動能守恒思想,根據梯度信息進行變量更新迭代求解,取代了傳統MCMC算法的隨機游走,也無需像VB-EM算法對高維矩陣求逆過程,因此計算更加便捷[21-23]。但考慮到HMC算法對非平滑成像后驗分布無法求解,因此本文在HMC算法的基礎上,引入近端算子[24-26],提出近端-HMC(proximal HMC,P-HMC)采樣求解算法。P-HMC算法能有效解決SAR成像貝葉斯非平滑后驗采樣問題,因而能實現可變先驗貝葉斯學習高分辨稀疏成像。

為了驗證所提VIP-Bayes算法的實用性與優越性,通過仿真與實測數據實驗,并與多種成像算法對比,驗證了所提算法的相對優越性。并利用相變熱力圖方法定量分析了所提算法在不同信噪比(signal to noise ratio, SNR)和降采樣率條件下的成像恢復性能。

1 SAR回波信號模型

根據SAR成像工作原理,推導出SAR成像信號模型。圖1為SAR成像幾何模型圖。本文以機載SAR成像為研究對象,假設雷達工作于聚束模式[27],如圖1所示,建立O-XYZ三維直角坐標系。

圖1 SAR 成像幾何模型Fig.1 Geometry model of SAR imaging

SAR載機平臺以一恒定速度v沿X軸正方向(方位向)預定航線飛行并不斷以正側視方向向地面目標場景發射線性調頻信號。地面場景如圖1所示,Q為場景中心。G點位于雷達正下方,第i個靜止目標散射點為Pi。q0(t)為參考天線相位中心。設地面目標散射點Pi到雷達天線相位中心q0(t)的斜距為Ri(t)=|R0+ri-q0(t)|,其中ri為Pi相對于場景中心Q的偏移矢量,R0為q0(t)到Q的參考斜距矢量,且R0=R0er,q0(t)=vtea,其中er、ea分別為距離向和方位向單位矢量。將Ri(t)在|R0-q0(t)|泰勒級數展開可得

(1)

(2)

載機雷達飛過地面場景后,地面靜止目標回波可表示為多散射點回波疊加形式:

S0(k,t)=

(3)

S0(k,t)=

(4)

式中：ka=-vt′/R0[14]；t′表示PFA插值后方位向時間變量。經距離壓縮后,可得距離壓縮域數據表達式為

(5)

式中：sinc(·)為回波的距離向包絡函數;λ為發射信號波長,為確定的雷達自身參數;t′為極坐標插值后的慢時間變量。

由式(5)可知,SAR回波可表示為距離向包絡與方位向線性相位乘積累計求和形式。因此,SAR回波可進一步表示為以下矩陣形式的線性方程:

Y=AX+CN

(6)

式中：Y∈CK×N為觀測回波；X∈CM×N為待成像的目標矩陣;M和N為方位向回波數和距離向采樣點數;CN∈CK×N為雜波或噪聲。在SAR模式下A∈CK×M為方位向傅里葉字典,構建如下:

(7)

式中：fd(·)為雷達運動在方位向產生的多普勒頻率。通常情況下,觀測回波可能為非完全數據,因此設K≤M。當K=M時,觀測回波為完全數據,即Y∈CM×N,與X維度相同;當K

2 SAR成像可變先驗統計建模

回波信號模型構建后,SAR高分辨成像問題可以轉化為求解式(6)線性方程中變量X的問題,因式(6)通常為非齊次欠定方程,因此不能直接使用X=A-1Y進行求解。SAR回波信號Y在距離壓縮域通常具有稀疏性,因此采用稀疏信號恢復方法。首先對噪聲進行建模,可建模為圓對稱復高斯分布,即似然函數為

(8)

2.1 廣義高斯先驗建模

傳統貝葉斯學習SAR稀疏成像先驗特征表征單一,先驗固化,本文提出可變成像先驗建模,即引入GGD先驗進行先驗表征。若以x為隨機變量,則標準的GGD先驗定義如下:

(9)

圖2 GGD概率密度分布Fig.2 Probability density function of GGD

(10)

式中:p=[p1,p2,…，pN]、α=[α1,α2,…，αN]為不同距離向單元的形狀參數和尺度參數;pn和αn分別為第n個距離向單元的形狀參數和尺度參數；‖·‖pn為lpn范數。考慮到目標先驗式(10)與似然分布式(8)不共軛,不能直接求解成像矩陣X的后驗分布。因此引入貝葉斯分層模型進行求解。

2.2 貝葉斯分層建模

(11)

(12)

式中:I[0,2]為定義在[0,2]的指示函數,即當pn∈[0,2]時,I[0,2](pn)=1,否則為0。因而GGD先驗可在[0,2]區間通過隨機采樣自學習獲得合適的形狀參數。同理,對于先驗分布的尺度參數α,由于對其先驗信息知之甚少,因此假設α服從Jeffreys無信息先驗,分布如下:

(13)

根據構建的層次貝葉斯模型,可推導獲得目標成像X的聯合后驗概率密度函數:

f(X,p,a,s2|Y)∝

(14)

根據貝葉斯分層模型,可繪制各參數或變量的有向無環圖,以對各變量相互關系進行表示,如圖3所示。

圖3 貝葉斯分層模型有向無環圖Fig.3 Directed acyclic graph of Bayesian hierarchical model

圖3中,Y為觀測數據,是顯隨機變量。X為待求解目標成像矩陣。p、α、σ2分別為GGD先驗分布的形狀參數、尺度參數、噪聲方差,均為隱隨機變量。a和b為逆伽馬分布形狀和尺度參數,是常量值。貝葉斯分層模型構建后,根據貝葉斯學習框架,需進一步對目標后驗分布進行采樣求解。

3 VIP-Bayes學習SAR成像算法

VIP-Bayes學習SAR成像算法中,目標成像建模后,后驗分布采樣求解為成像算法的關鍵步驟。統計采樣技術經過近幾十年的發展,目前已發展出M-H(Metropolis-Hasting)采樣、Gibbs采樣等經典采樣算法,近年來也出現HMC等新型采樣算法[21]。其中,Gibbs采樣算法適合于常規統計分布采樣。而M-H采樣適用范圍廣,但效率較Gibbs采樣算法低。新型HMC算法依據哈密頓動力學思想,可實現更復雜的凸分布函數采樣。考慮到VIP先驗下所得目標成像矩陣聯合后驗分布過于復雜,無法直接采樣求解。因此,本文將HMC采樣算法引入貝葉斯學習后驗采樣關鍵步驟中,以實現成像目標矩陣X的求解。根據式(14),首先分別推導各隨機變量的邊緣后驗分布。

3.1 超參數后驗分布計算

由式(8)和式(11),針對噪聲方差進行邊緣后驗推導,可得

f(s2|Y,X,p,α)∝

(15)

式(19)為常規逆伽馬分布,可直接根據該后驗分布式進行Gibbs采樣。對于可變GGD先驗的形狀參數和尺度參數,由式(10)、式(12)和式(13)可得邊緣后驗分布如下:

f(p|Y,X,α,s2)∝

(16)

(17)

式(16)中,指數項中含有lpn范數,為非常規概率分布,因此Gibbs采樣算法無法進行采樣求解。對于此類非常規分布,可使用M-H算法進行采樣求解,使形狀參數適應目標特征。針對式(17)尺度參數條件后驗分布,因其為逆伽馬分布,可直接根據式(17)進行Gibbs采樣。

3.2 成像結果后驗分布計算

對于待求解的目標成像矩陣X,根據式(8)、式(10)進行推導,獲得其邊緣后驗分布:

f(X|Y,p,α,s2)∝

(18)

式中:指數項中涉及l2范數及lpn范數求和。對于此類復雜非常規分布,Gibbs算法無法采樣求解,而M-H算法通過無方向性隨機游走采樣求解,效率較低。針對此問題,本文引入HMC采樣貝葉斯學習算法進行求解。

HMC采樣算法是一種基于哈密頓動力學思想提出的算法,2006年,Bishop將其引入統計機器學習中[21]。HMC算法以哈密頓方程為基礎,基于動能和勢能守恒原理提出,可實現如下分布采樣:

f(x)∝exp[-U(x)]

(19)

式中:U(x)為勢能函數,x為位移矢量。動能函數為K(q)=qTq/2σ2,q為物體動量矢量,HMC算法根據梯度信息進行蛙跳步驟[22]更新變量。

針對本文SAR高分辨成像問題中復雜分布式(24),根據HMC采樣原理,勢能函數為

(20)

考慮到勢能函數中‖X:n‖pnpn項可能出現l1范數,使得U(X)在某些點出現不可微情況,則HMC算法無法進行求解。因此筆者在HMC采樣貝葉斯學習算法基礎上,引入近端算子進行梯度近似,提出P-HMC采樣算法。接下來將對P-HMC算法對目標成像結果采樣策略進行詳細闡述。

3.3 近端-HMC采樣策略

近端方法是一種針對于非平滑函數進行的梯度近似策略[22]。設G(v)為待求解的非平滑目標函數,通常函數G(v)在x處的近端算子可表示為

(21)

式中:λ為近端尺度參數。通過近端算子的引入,HMC算法成為P-HMC算法,可解決非平滑勢能函數U(X)在迭代中無法求解梯度的問題。本文SAR成像中勢能函數U(X)可分解為

(22)

式中:

其中,L(X:n)為可微函數,可直接求解梯度;G(X:n)為不可微函數,需借助近端算子求解:

(23)

式中:V:n為X:n的替代變量,若pn=1,即‖X:n‖pnpn表示X:n的l1范數,此時proxλ G(X:n)為復數軟閾值[24]。若pn≠1且G(X:n)非平滑,可使用前后向算法進行迭代近似求解[25]。

對于平滑函數項L(X:n),求解其梯度如下:

(24)

根據近端算子與梯度的關系[22],本文中SAR目標成像矩陣X中任意第n個距離向采樣元素Y:n使用P-HMC算法進行變量更新,第r次迭代時更新如下:

(25)

(26)

(27)

接受概率表示為

(28)

式中:

綜上,根據超參數與目標成像結果后驗計算與采樣求解算法,可總結出VIP-Bayes成像算法流程。

4 VIP-Bayes成像算法流程

SAR目標成像可以理解為求解成像結果矩陣X的過程。圖4為VIP-Bayes高分辨SAR成像算法流程。主要包括參數與變量初始化、超參數采樣求解與成像結果求解3個步驟。具體如下。首先,對各變量及參數進行初始化,設置算法總迭代次數為S,Buin-in期門限為T。Buin-in期指貝葉斯學習采樣求解過程中樣本收斂至目標分布前的階段,該階段樣本不穩定,應舍去。

圖4 VIP-Bayes成像算法流程圖Fig.4 Solution flowchart of VIP-Bayes imaging algorithm

接著,進行3個超參數采樣求解,如圖4中藍色虛線框所示。即根據式(15)對噪聲方差進行Gibbs采樣,根據式(16)和式(17)對可變成像先驗的形狀參數、尺度參數分別使用M-H算法、Gibbs采樣算法進行采樣求解。當3個參數更新的迭代次數大于Buin-in門限T之后,得到的樣本則為可用的收斂后樣本值。最后,對于成像結果矩陣X,可對距離向第n個單元成像結果Xn依次進行采樣求解,最終合成目標成像結果。對成像結果矩陣X的采樣求解使用前文所提的P-HMC算法,如圖4中紅色虛線框所示。該采樣步驟中,首先進行動量變量q的初始采樣,使其服從高斯分布。然后求得第r次迭代時q和X:n,并依據式(25)～式(27)進行變量更新迭代,再通過式(28)中“接受-拒絕”步驟篩選獲得的目標成像結果樣本。最后,通過對所得樣本求期望,即可獲得高分辨SAR成像結果。

5 實驗驗證

為驗證算法的有效性,本文首先選取仿真數據進行實驗。通過與經典RD算法、SBL算法、l2范數正則約束下ADMM凸優化算法(以下簡稱ADMM算法)進行成像對比實驗,驗證了所提VIP-Bayes成像算法的優越性。最后,本文利用相變熱力圖分析法,定量衡量各算法在不同信噪比與降采樣率下的恢復性能趨勢,從而驗證所提算法的優越恢復性能。

5.1 仿真數據實驗

首先使用一組仿真數據進行實驗,仿真參考恢復圖像如圖5所示,該數據由57個散射點構成,有序分布于中間十字區域內。實驗中仿真實驗參數如表1所示。

表1 仿真實驗參數設置表

圖5 仿真數據參考恢復圖像Fig.5 Simulation data reference restored image

根據上述仿真參數,在仿真數據中加入噪聲并進行降采樣處理,分別使用RD算法、SBL算法、ADMM算法及VIP-Bayes算法進行成像恢復。實驗中設置SNR=5 dB,降采樣率為0.8。恢復后的結果如圖6所示。其中圖6(a)為RD算法恢復結果,圖6(b)～圖6(d)分別為SBL、ADMM算法和VIP-Bayes算法恢復的結果,對比可以看出,RD算法恢復結果中仍存在大量背景噪聲。而SBL和ADMM算法恢復結果中,雖然噪聲被大部分消除,但十字場景中部分散射點細節信息丟失。VIP-Bayes算法成像結果最佳,消除噪聲的同時很好地保留了目標特征的完整性。

圖6 仿真數據成像結果Fig.6 Imaging data results of simulation

為了定量分析仿真數據成像結果分辨率性能,根據上述仿真實驗參數,選取其中單個散射點進行點目標仿真成像獲得仿真結果如圖7所示,其中圖7(a)為單個點目標成像的三維結果圖,形狀為二維Sinc函數。圖7(b)和圖7(c)分別為距離向和方位向頻譜圖,均為未加窗函數抑制旁瓣結果。通過圖7(a)測量獲得距離向分辨率為1.1 m,方位向分辨率為0.87 m,與理論計算值接近。在分辨率分析中,峰值旁瓣比(pulse sidelobe ratio,PLSR)常用來評估成像性能。從圖7(b)和圖7(c)獲得距離向和方位向分別為-13.23 dB和-18.22 dB。根據上述點目標分辨率分析結果,該實驗充分證明了所提算法在仿真數據上的有效性。

圖7 點目標成像結果Fig.7 Results of point target imaging

5.2 實測數據實驗

SAR實測數據選取美國Sandias實驗室的機載SAR實測數據集。實驗選取成像場景為停靠的靜止飛機,降采樣率為0.6,即K/M=0.6,通過方位向的隨機降采樣后,分別使用RD、 SBL、 ADMM算法與本文所提算法進行成像對比。結果如圖8所示。

對比成像結果可看出,圖8(a)RD算法成像結果中噪聲較多,對目標成像不夠清晰。圖8(b)中已經無法辨識成像目標,只能觀察到極少數強散射點。圖8(c)使用ADMM算法進行成像,可以看出目標整體結構特征較為完整,但噪聲抑制能力不足,仍存在一定量的背景噪聲。

圖8 SAR實測數據成像結果Fig.8 SAR imaging results based on raw data

圖8(d)為使用VIP-Bayes算法進行成像獲得的結果。從圖8中可以看出,結果中保留并增強了目標飛機的強散射稀疏特征,又保留了飛機整體的全局結構特征,同時噪聲去除較為徹底,能很清晰地識別成像飛機地機翼,尾翼機身等各部分結構。該實驗結果很好地體現了可變成像先驗對復雜目標場景表征的動態靈活性與算法的有效性。

5.3 相變分析實驗

相變分析方法最早由Donoho提出,可用于評估各算法在不同參數下成像恢復性能。在SAR成像中,常用相變熱力圖(phase transition diagram,PTD)來衡量某一算法在不同SNR、降采樣參數下的恢復概率[30-31]。本實驗中,使用如圖5所示的仿真數據,分別應用SBL算法, ADMM算法和 VIP-Bayes算法進行成像恢復實驗,并將恢復后的成像結果與圖5參考圖像進行相關度計算,相關度越接近1則說明成像恢復性能最佳,即恢復概率越大。在本實驗中,恢復概率從大到小分別使用colorbar中由大到小的顏色值表示。實驗中設置的恢復閾值為0.85,即概率大于0.85的區域為可恢復區。相變圖實驗結果如圖9所示,其中橫坐標為降采樣率,從0到1變化。縱坐標為SNR,變化范圍為[-10,11]。紅色部分表示算法對目標場景恢復概率接近于1,藍色則表示恢復概率接近于0。對比3種算法相變熱力圖可看出,圖10(a)紅色面積最小,經計算可得可恢復概率為19.76%;圖9(b)次之,可恢復區對應概率為22.86%;圖9(c)可恢復區面積最大,可恢復概率為40.95%。因此,實驗表明VIP-Bayes成像算法恢復性能優于其他兩種算法。另外,在SNR較低的區域,如圖9(c)所示,當SNR=-2 dB時,圖中仍有一半區域為紅色,說明VIP-Bayes算法在SNR較低時仍有不錯的成像恢復性能。

圖9 不同算法PTD結果Fig.9 PTD results of different algorithms

6 結 論

本文針對高分辨SAR成像中,傳統貝葉斯學習成像算法先驗固化問題,提出一種可變成像先驗貝葉斯學習高分辨成像算法。該算法基于動態可變GGD先驗對目標特征進行靈活表示,并引入P-HMC采樣算法應用于復雜目標成像后驗求解。本文實驗部分對SAR實測數據進行成像實驗,證明了所提算法相比于其他傳統算法具有較好的成像恢復性能。另外,發現凸優化類算法和貝葉斯類算法具有許多共同特征,凸優化類算法正則項與貝葉斯學習算法先驗分布具有一定的等價性。凸優化類算法正則項多樣,而統計學習類算法已知可用統計先驗較少,但無需繁瑣參數調整過程。若能將凸優化類算法與高效統計采樣算法結合起來,充分發揮兩者優勢,這將會具有很廣闊的應用前景,也將是后續研究工作重點。