徐安林, 張 毓, 周 峰

(1. 中國人民解放軍63921部隊, 北京 100094; 2. 西安電子科技大學(xué)電子工程學(xué)院, 陜西 西安 710071)

0 引 言

逆合成孔徑雷達(inverse synthetic aperture radar,ISAR)能夠全天時、全天候地提供飛機和衛(wèi)星等非合作目標(biāo)的高分辨二維圖像,并從中獲取其尺寸、結(jié)構(gòu)等形狀特征及運動特征[1]。因此,高分辨ISAR成像在空間目標(biāo)監(jiān)視與彈道目標(biāo)防御中發(fā)揮著重要作用[2]。

在ISAR觀測中,通過發(fā)射大時寬-帶寬積信號及脈沖壓縮技術(shù)即可獲得較高的距離分辨率,而方位高分辨則需要通過方位向的長時間相干積累實現(xiàn)。然而在實際情況中,ISAR往往面臨著非常復(fù)雜的環(huán)境,導(dǎo)致無法對目標(biāo)進行長時間連續(xù)觀測,從而產(chǎn)生方位缺損。同時,受到強噪聲或干擾的影響,目標(biāo)的回波通常存在低信噪比的問題。此外,由于傳統(tǒng)平動補償中的自聚焦方法精度不足,并且信號傳輸過程中會受到大氣影響,目標(biāo)回波中往往存在著隨機相位誤差[3]。面對上述情況,傳統(tǒng)成像方法的性能會迅速下降甚至失效。因此,迫切需要研究適用于復(fù)雜觀測環(huán)境的高分辨ISAR成像方法。

在光學(xué)區(qū),目標(biāo)滿足點散射模型,其ISAR像在距離域和方位域具有稀疏性。因此,可以利用稀疏信號重構(gòu)方法實現(xiàn)復(fù)雜觀測環(huán)境下的高分辨ISAR成像[4-7]。這類方法首先構(gòu)建稀疏觀測模型,進而利用數(shù)值優(yōu)化或貝葉斯方法從低維回波信號中重構(gòu)出目標(biāo)散射點分布?，F(xiàn)有基于數(shù)值優(yōu)化的求解方法雖然計算量較小,但在低信噪比情況下成像結(jié)果存在散焦[8-12]。傳統(tǒng)參數(shù)化貝葉斯方法[13-17]往往根據(jù)現(xiàn)有經(jīng)驗來選取先驗分布,導(dǎo)致其數(shù)據(jù)描述的靈活性弱,重構(gòu)誤差較大。而非參數(shù)貝葉斯方法將隨機過程作為模型的先驗分布,使得其參數(shù)空間能夠隨觀測數(shù)據(jù)變化,從而獲得比參數(shù)化貝葉斯方法更高的靈活性,為解決稀疏信號的準(zhǔn)確重構(gòu)問題提供了新思路。

本文提出一種基于非參數(shù)貝葉斯的ISAR高分辨成像方法。該方法首先引入Beta過程非參數(shù)先驗構(gòu)建概率圖模型,并通過Gibbs采樣求解模型的后驗分布,同時利用最大似然方法估計相位誤差,進而實現(xiàn)低信噪比、稀疏孔徑等復(fù)雜觀測環(huán)境下的高分辨聚焦成像。實驗結(jié)果表明,該方法能夠在復(fù)雜觀測環(huán)境下獲得優(yōu)于現(xiàn)有典型數(shù)值優(yōu)化方法及參數(shù)化稀疏貝葉斯重構(gòu)方法的成像結(jié)果。

1 信號模型

1.1 信號建模

對于平穩(wěn)運動的小轉(zhuǎn)角目標(biāo)[18],其散射點的瞬時斜距可以表示為ΔRk(tm)≈xkωrottm+yk[19],其中ωrot表示目標(biāo)旋轉(zhuǎn)角速度,xk和yk分別表示第k個散射點的橫坐標(biāo)和縱坐標(biāo),tm表示慢時間。

令距離脈壓后的回波矩陣為Y∈CNa×Nr,其中Na表示方位單元數(shù),Nr表示距離單元數(shù)。同時,假設(shè)每個距離單元的回波相互獨立,考慮到隨機相位誤差的距離不變性,則第n個距離單元回波s(m)滿足:

(1)

式中:n∈[1,Nr];m為方位單元序號;k∈[1,K]表示散射點的序號;Ak表示散射點k的幅度,且Aka=Akexp(-j4πyk/λ),ar(yk)表示與yk有關(guān)的sinc波形函數(shù);exp(jφm)表示隨慢時間變化的隨機相位誤差;ωk=-2xkωrot/(λ·PRF)表示第k個散射點的多普勒,PRF為脈沖重復(fù)頻率,λ為波長;ε(tm)為加性噪聲。

當(dāng)回波存在缺損時,令可觀測回波序號為ma,且可觀測回波數(shù)為Ma。對于第n個距離單元,其稀疏觀測模型可表示為

ya=EaAaxa+εa

(2)

式中:ya=[s(ma1),s(ma2),…,s(maMa)]T表示可觀測回波向量;Ea=diag{exp(jφma1),exp(jφma2),…,exp(jφmaMa)}表示相位誤差矩陣;xa表示待重構(gòu)的方位像;εa為噪聲向量;Aa∈CMa×Da表示過冗余字典,且其第i行可以表示為

Ai·=exp(j2πωdamai),i∈[1,Ma];da∈[1,Da]

(3)

式中:Da表示構(gòu)建的多普勒網(wǎng)格數(shù)。當(dāng)Ma≤Da時可實現(xiàn)超分辨成像。

重新在實數(shù)域定義信號模型得到:

y=Φω+ε

(4)

1.2 概率建模

為了充分利用散射中心的稀疏特性,引入Beta過程非參數(shù)貝葉斯先驗[20]。若隨機測度H服從參數(shù)為a、b和H0的Beta過程,那么H可由下式表示:

(5)

并記作H～BP(αH0),其中K→∞,delta函數(shù)δBk(B)當(dāng)且僅當(dāng)B=Bk時為1;否則為0。同時,每一個獨立同分布地從概率測度H0中生成的原子Bk都有一個權(quán)值πk與之對應(yīng)。

若某個無限維的二值列向量z的第k個元素zk服從Bernoulli分布:

zk～Bernoulli(πk)

(6)

那么新的測度,X(B)=∑kzkδBk(B)則由參數(shù)為H的Bernoulli過程產(chǎn)生,記作X～BeP(H)。

在式(4)所示的稀疏信號重構(gòu)問題中,y∈RP表示實數(shù)域有效回波向量,Φ∈RP×K表示實數(shù)域多普勒字典,權(quán)向量ω∈RK即為實數(shù)域待重構(gòu)方位像。首先,將權(quán)向量ω分解為兩個向量元素相乘的形式,即

y=Φ(s°z)+ε

(7)

式中:°表示向量的哈達瑪積?？梢钥闯?ω的稀疏性由二值向量z決定,且其幅值由向量s決定。本文利用Beta過程獲取稀疏二值向量z。

令二值向量z的先驗分布滿足Bernoulli分布,即

zk～Bernoulli(πk)

(8)

且參數(shù)πk服從如下的Beta分布:

(9)

此外,給幅度向量s和噪聲向量ε引入高斯先驗:

(10)

(11)

為了使模型更加靈活,再對s和ε的精度參數(shù)γs和γε引入Gamma分布先驗:

γs～Gamma(c,d)

(12)

γε～Gamma(e,f)

(13)

根據(jù)上述概率建模,基于Beta過程線性回歸(Beta process linear regression, BPLR)的層級概率模型可以表示為

(14)

其對應(yīng)的概率圖模型如圖1所示。

2 模型求解及分析

由于各距離單元對應(yīng)的方位像及相位誤差矩陣均未知,因此本文利用循環(huán)迭代的方法進行優(yōu)化求解。首先通過Gibbs采樣算法求得待重構(gòu)的稀疏方位像,再通過最大似然方法估計相位誤差。交替進行上述兩個步驟,直到滿足收斂條件。

2.1 基于Gibbs采樣的方位像估計

分別對模型各參數(shù)進行采樣。

(1) 對z=[z1,z2,…,zK]采樣,令Φ=(φ1,φ2,…,φK),其中zk后驗分布為Bernoulli分布:

(15)

(16)

p0=1-πk

(17)

(2) 對sk=[s1,s2,…,sK]采樣,其中sk服從高斯分布:

(18)

方差Σsk和均值μsk可以表示為

(19)

(20)

(3) 對πk采樣,πk服從Beta分布:

(21)

(4) 對γs采樣,γs服從Gamma分布:

(22)

(5) 對γε采樣,γε同樣服從Gamma分布:

(23)

2.2 超參數(shù)設(shè)置

由前文公式的推導(dǎo)可以看出,Gibbs采樣的結(jié)果依賴于超參數(shù)a、b的設(shè)置。因此,分別定義超參數(shù)a、b:

(24)

(25)

2.3 相位誤差估計

在完成所有距離單元的方位像重構(gòu)后,即可獲得本次迭代的ISAR圖像X。在此基礎(chǔ)上,給有效回波矩陣Y0∈CMa×Nr的每一列引入相互獨立的多維復(fù)高斯分布,可得

(26)

式中:Y·r和X·r分別表示Y0和X的第r列。

精度α0服從如下的Gamma分布:

p(α0)=Gamma(α0|v1,v2)

(27)

Ea的似然函數(shù)為

(28)

通過最小化負對數(shù)似然函數(shù),可得Ea的估計值為

(29)

(30)

式中:Yi·和Ai·分別表示Y0和A的第i行。

2.4 算法流程

交替進行方位像估計和相位誤差估計,直到收斂。特別地,將收斂條件定義為

(31)

最終,基于Beta過程的高分辨ISAR成像流程圖如圖2所示,其中T表示采樣次數(shù)。在算法初始化的步驟中,設(shè)置相位誤差矩陣Ea為單位陣,二值向量z為全零向量;對s的元素進行隨機初始化,并將π設(shè)為較小的值。同時,將超參數(shù)c、d、e和f設(shè)為10-4,S和F的初始值滿足1

3 實驗及分析

本節(jié)首先進行蒙特卡羅實驗,對所提算法重構(gòu)的準(zhǔn)確性以及對于噪聲的穩(wěn)健性進行分析。然后利用YAK-42飛機實測數(shù)據(jù)的成像結(jié)果驗證所提算法的有效性

3.1 蒙特卡羅實驗

本節(jié)通過一維稀疏信號重構(gòu)實驗將所提算法的性能與正交匹配追蹤(orthogonal matching pursuit, OMP),L1范數(shù)優(yōu)化[21],基于變分推斷(variational inference, VI)的相關(guān)向量機(relevance vector machine, RVM),即RVM-VI進行比較。其中,前兩種方法為典型數(shù)值優(yōu)化方法,RVM為典型的參數(shù)化稀疏貝葉斯學(xué)習(xí)方法,該方法采用Gamma-Gaussian先驗進行稀疏建模。

(32)

由圖3(b)可以看出,隨著信噪比增大,4種算法重構(gòu)誤差逐漸減小,并且在所有信噪比下,所提方法均能獲得最小的重構(gòu)誤差,其次是RVM、L1范數(shù)優(yōu)化方法以及OMP方法。由于稀疏貝葉斯學(xué)習(xí)方法通過引入概率分布獲得了高階統(tǒng)計信息,因此其重構(gòu)誤差明顯小于OMP、L1范數(shù)優(yōu)化等典型的數(shù)值優(yōu)化方法。同時,由于Beta過程先驗?zāi)軌蚋鶕?jù)數(shù)據(jù)特性自適應(yīng)調(diào)整參數(shù)空間,從而獲得比RVM等參數(shù)化貝葉斯方法更高的靈活性[22],因此其重構(gòu)誤差最小。

3.2 實測數(shù)據(jù)成像結(jié)果

本節(jié)分別采用存在隨機相位誤差的YAK-42飛機稀疏孔徑及短孔徑實測數(shù)據(jù)驗證所提算法有效性。

稀疏孔徑數(shù)據(jù)如圖4(a)所示,回波距離和方位向采樣點數(shù)分別為256和512,信噪比為-3 dB,回波方位缺損率為50%。各算法成像結(jié)果如圖4(b)～圖4(e)所示。由圖可知,OMP、L1范數(shù)優(yōu)化等數(shù)值優(yōu)化方法無法有效抑制噪聲,雖然基于RVM的成像結(jié)果虛假點較少,但目標(biāo)輪廓不夠清晰。所提方法的成像結(jié)果背景干凈且聚焦良好。

表1給出了4種算法稀疏孔徑成像結(jié)果的圖像熵?？梢钥闯?所提方法的ISAR圖像熵最小,表示其聚焦效果最優(yōu)。

表1 稀疏孔徑不同算法的圖像熵

對于短孔徑數(shù)據(jù),目標(biāo)有效回波所占方位單元數(shù)為200,回波缺損率61%,信噪比為-3 dB。各算法成像結(jié)果如圖5(b)～圖5(e)所示。

由圖5可知,OMP、L1范數(shù)優(yōu)化等數(shù)值優(yōu)化方法對應(yīng)圖像包含較多虛假點,同時RVM的成像結(jié)果目標(biāo)輪廓不夠清晰。而所提方法成像結(jié)果背景干凈且聚焦良好。4種算法在短孔徑情況下的圖像熵如表2所示,對比可知所提方法的ISAR像熵最小,聚焦效果最優(yōu)。

表2 短孔徑不同算法的圖像熵

4 結(jié) 論

本文提出基于Beta過程的高分辨ISAR成像方法。該方法首先將高分辨成像問題轉(zhuǎn)化為稀疏信號重構(gòu)問題,進而通過引入Beta過程先驗構(gòu)建層級概率模型。在此基礎(chǔ)上,交替利用Gibbs采樣算法求得稀疏方位像并利用最大似然方法估計隨機相位誤差,最終得到聚焦良好的高分辨ISAR像。實驗結(jié)果證明在低信噪比、稀疏孔徑、短孔徑等復(fù)雜觀測條件下,所提算法性能優(yōu)于現(xiàn)有成像方法。

未來將構(gòu)建更靈活的概率模型,并尋找高效的模型求解方法。同時研究存在強干擾等復(fù)雜環(huán)境中機動目標(biāo)的高分辨成像方法。