基于近似觀測和最小熵約束的SAR稀疏自聚焦方法

熊旭穎, 李 根, 馬彥恒

(陸軍工程大學石家莊校區(qū)無人機工程系, 河北 石家莊 050003)

0 引 言

壓縮感知和稀疏重構理論表明,當成像場景稀疏或可以被稀疏表示時,可以利用少量的雷達回波數(shù)據(jù)獲取高分辨和高質量的合成孔徑雷達(synthetic aperture radar,SAR)圖像,稀疏SAR成像技術能夠有效降低雷達接收機的數(shù)據(jù)采集與存儲的壓力,相關的成像方法受到了廣泛關注[1-4]。

傳統(tǒng)的稀疏SAR成像方法需要根據(jù)平臺的運動軌跡構建測量矩陣。對于機載SAR成像,由于慣導設備精度有限,平臺的運動軌跡通常難以精確測量,運動誤差的存在將導致成像理論模型和測量矩陣失配,這種失配將產生嚴重的散焦現(xiàn)象。因此,基于回波數(shù)據(jù)的自聚焦是稀疏SAR成像過程中必不可少的環(huán)節(jié)。Onhon等人在2012年提出了一種稀疏驅動的自聚焦方法[5],該方法將誤差相位引入測量矩陣，構建了一種稀疏自聚焦模型,并將模型的迭代求解分為SAR場景重建和誤差相位估計兩步,能夠校正一維和二維的相位誤差,該方法也被稱為基于兩步優(yōu)化的稀疏自聚焦方法?；趦刹絻?yōu)化的自聚焦思想,Yang等人提出了能夠估計觀測位置誤差的壓縮感知SAR成像方法[6]。Chen等人提出了一種基于參數(shù)化稀疏表示的SAR運動補償方法[7],進一步提高了成像質量。以上這些方法的自聚焦模型均是在時域構建測量矩陣,這種測量矩陣又被稱為精確觀測算子。當成像場景較大時,精確觀測算子將產生巨大的運算量和內存占用量[8-9],限制了稀疏SAR成像方法的應用。近年來,中科院電子所吳一戎院士團隊研究了基于近似觀測算子的稀疏微波成像方法[7，10],成功地將傳統(tǒng)的匹配濾波類SAR成像算法和壓縮感知理論相結合,極大地減小了測量矩陣的尺寸,使稀疏SAR成像方法能夠用于大場景成像,相關的成像算法受到了廣泛的關注[11-13]。對于運動誤差條件下稀疏SAR成像,Li等人將近似觀測算子引入稀疏自聚焦模型中,用于替代精確觀測算子,極大地降低了稀疏自聚焦方法的運算量[14]。Pu等人基于近似觀測算子構建了一種基于幅相誤差聯(lián)合校正的稀疏自聚焦模型,用于校正空變的相位誤差,在全采樣數(shù)據(jù)下有效估計了三維的運動誤差偏離[15]。

以上這些基于兩步優(yōu)化的稀疏自聚焦方法均采用最大似然估計器(maximum likelihood estimator,MLE)估計相位誤差,在該模型中,誤差相位在二維時域(回波數(shù)據(jù)域)引入,難以獲取誤差相位的初始解。在基于兩步優(yōu)化的稀疏自聚焦方法中,合適的誤差相位初始解對算法性能至關重要,當回波數(shù)據(jù)采樣率較低且重建散射點數(shù)量較多時,MLE收斂緩慢并且容易陷入誤差較大的局部最優(yōu)解,導致自聚焦失敗。本文基于近似觀測算子構建了一種新的稀疏自聚焦模型,在該模型中,誤差相位在聚焦圖像的傅里葉變換域引入,因此相位梯度自聚焦(phase gradient autofocus,PGA)[16]等經(jīng)典匹配濾波類自聚焦方法可用于提供誤差相位的有效初始解。所提模型基于兩步迭代的思想采用坐標下降法進行求解,為加快算法收斂速度,快速迭代閾值收縮(fast iterative shrinkage thresholding,FIST)[17]被用于重建場景散射系數(shù),同時引入了場景散射系數(shù)的最小熵約束,這一附加的先驗信息在加快算法迭代收斂速度的同時,使迭代結果更接近全局最優(yōu)解。

1 傳統(tǒng)的稀疏自聚焦模型

假設成像區(qū)域中共含有PQ個散射點,在一個合成孔徑時間內,SAR接收機采集到的離散二維回波數(shù)據(jù)可表示為

exp[-j4πR0(k,ηm)/λ]，

m=1,2，…，M;n=1,2，…，N

(1)

式中：ηm和τn分別表示離散的方位和距離采樣點；M和N分別表示距離和方位向上的采樣點數(shù)；g(k)表示第k個散射點的復散射系數(shù)；R0(k,ηm)表示成像區(qū)域中的第k個散射點在方位采樣時刻ηm的瞬時斜距；c和λ分別表示光速和發(fā)射信號波長；p(τ)表示發(fā)射的基帶線性調頻信號,具體為

p(τ)=rect(τ/Tr)exp(jπKrτ2)

(2)

式中：Tr表示發(fā)射信號的脈沖持續(xù)時間；Kr表示距離調頻率。

在考慮噪聲的情況下,式(1)可矩陣化表示為

s=Ξg+n0

(3)

式中：s∈CMN×1表示由2維回波組成的1維回波向量；Ξ=HA∈CMN×PQ表示精確觀測矩陣,其中H∈CMN×MN為降采樣矩陣,是對角元素僅含有0或1的對角矩陣,A∈CMN×PQ為測量矩陣,其具體形式見文獻[1]；g∈CPQ×1為散射點組成的1維向量；n0為1維的加性高斯白噪聲向量。

式(3)是基于精確觀測的壓縮感知模型,當成像場景稀疏時,散射系數(shù)向量g可以通過求解如下所示的L1范數(shù)正則化問題獲得：

(4)

式中：γ表示正則化參數(shù),由成像場景的稀疏度決定。

當SAR回波數(shù)據(jù)存在非空變的相位誤差時,常規(guī)的稀疏自聚焦模型在二維時域引入誤差相位[5],表示為

(5)

其中，Ψ為相位誤差矩陣,具體表示為

(6)

式中:φ(m)表示第m個方位采樣點的誤差相位。

文獻[5]將式(5)的求解轉化為兩步優(yōu)化問題,并采用坐標下降法進行求解,具體的求解過程如算法1所示。

算法 1 基于精確觀測的稀疏自聚焦算法初始化: i=0, g0=(Ψ0Ξ)Hs且Ψ0=Ψ0|?(m)=0步驟 1 gi+1=argming12s-ΨiΞg22+γg1{}步驟 2 Ψi+1=argminΨ12s-ΨΞgi+122+γgi+11{}步驟 3 i=i+1,然后回到步驟1。當gi+1-gi22gi22小于預設的門限值時,算法停止。

算法1中的步驟1在已知相位誤差的基礎上,對場景散射系數(shù)進行重建,是一個L1范數(shù)正則化問題,可采用迭代軟閾值(iterative soft thresholding,IST)[18]等方法進行求解;步驟2在已知場景散射系數(shù)的基礎上估計誤差相位,可以通過最小化均方誤差,利用MLE進行求解,則第m個方位采樣點的誤差相位為

(7)

其中，[Ξ]m和[s]m分別表示第m個方位采樣點的觀測子矩陣和回波數(shù)據(jù)。

算法1中的稀疏自聚焦模型需要將2維的散射場景和回波數(shù)據(jù)拉成1維向量,精確觀測矩陣Ξ的維度為MN×PQ,這將產生巨大的極大內存占用和運算量,限制了其在大場景成像中的應用。近似觀測算子的提出極大地降低了觀測矩陣的維度,提高了稀疏重構效率[11]。接下來,以線性調頻變標算法為例,簡要說明近似觀測算子的構造過程。

圖1給出了線性調頻變標算法的成像流程圖,其中，F(xiàn)FT表示快速傅里葉變換，IFFT表示逆快速傅里葉變換。用Γ(·)算子表示該成像過程,則Γ(·)可以矩陣化表示為

圖1 線性調頻變標算法成像流程圖Fig.1 The flow chart of chirp scaling algorithm for the linear frequency modulation

(8)

從式(8)中可以看出Γ(·)是僅含矩陣乘法和矩陣點乘運算的線性可逆算子,原始回波S可通過其逆運算獲得,即

(9)

其中，*表示矩陣的共軛。

在式(9)中，回波數(shù)據(jù)通過匹配濾波成像的逆過程獲得,Γ-1(·)被稱為近似觀測算子。計算Γ(·)和Γ-1(·)所需的維度與原始回波數(shù)據(jù)一致,極大地降低了觀測算子的內存占用和運算量。

為降低稀疏自聚焦模型中觀測矩陣Ξ的維度,可通過引入近似觀測算子構造如下式所示的稀疏自聚焦模型：

(10)

式中：L是一個僅包含0和1的二進制隨機矩陣,表示回波數(shù)據(jù)的欠采樣矩陣。設ξa表示方位向上的回波數(shù)據(jù)抽取率,ξr表示距離向上的回波數(shù)據(jù)抽取率,則L共有Nξa個非零列,每一列有Mξr個1。Φ表示相位誤差矩陣,具體形式為

Φ=diag[ejφ(1),ejφ(2),…,ejφ(M)]

(11)

式(10)中基于近似觀測的稀疏自聚焦模型同樣可以采用算法1進行求解,盡管該模型利用近似觀測算子降低了內存占用量,但誤差相位依舊在回波數(shù)據(jù)域引入,難以采用常規(guī)的自聚焦方法獲取有效的誤差相位初始解。誤差相位需要從0開始采用MLE進行迭代更新,當采樣率較低且重建的散射點數(shù)量較多時,該模型的迭代過程容易陷入誤差較大的局部最優(yōu)解。

2 基于近似觀測和最小熵約束的稀疏自聚焦

2.1 新的稀疏自聚焦模型

(12)

其中

(13)

在式(12)的稀疏自聚焦模型中,相位誤差Φ在聚焦圖像的傅里葉變換域引入,常規(guī)的匹配濾波類自聚焦方法可以用于提供誤差相位的初始解。

2.2 基于兩步優(yōu)化的模型求解方法

基于算法1中兩步優(yōu)化的思想,式(12)中所提的稀疏自聚焦模型可以采用坐標下降法進行求解,求解過程的偽代碼如算法2所示,主要內容包括誤差相位初始解Φ0的獲取,以及G和Φ的重建方法,具體的推導過程如下。

算法 2 基于近似觀測和最小熵約束的稀疏自聚焦算法初始化: i=0, Imax,V0=0, t0=1,Φ0=PGA(Γ(L☉S))步驟 1 基于FIST的場景散射系數(shù)重建Gi+1=argminGJ(G, Φi)步驟 1.1 Vi+1←Ψλ,δ(Vi-δFη(Φ?iΩ(L☉[Ω-1(Φi·F-1(Gi))-S])))步驟 1.2 ti+1=1+1+4t2i2步驟 1.3 Gi+1=Vi+1+ti-1ti+1(Vi+1-Vi)步驟 2 相位誤差估計Φi+1=argminΦJ(Gi+1, Φ)步驟 2.1 ?i(m)=angle([Ω(S)]m·[F-1η(Gi+1☉(1+log|Gi+1|2))]Hm)步驟 2.2 Φi=diag[ej?i(1),ej?i(2),…,ej?i(M)]步驟 2.3 i=i+1,然后重復步驟1～步驟5,當i=Imax時,算法停止。

2.2.1 誤差相位初始解的獲取

采用坐標下降法交替求解G和Φ時,僅能保證迭代收斂到局部最優(yōu)解。因此,合適的Φ0能夠減少算法的迭代次數(shù)并使收斂結果更接近全局最優(yōu)解。在所構建的稀疏自聚集模型中,PGA和最小熵等方法均可以用于提供Φ0。但在稀疏采樣的情況下,點目標聚焦后的點散布函數(shù)會出現(xiàn)較高的副瓣電平,大量高副瓣電平的存在會使聚焦后的SAR圖像存在大量的欠采樣噪聲,導致圖像熵值增加,因此最小熵方法不適用于稀疏采樣條件下Φ0的估計。稀疏采樣對強散射點的主瓣影響較小,在PGA方法中,隨著迭代次數(shù)的增加,窗長度逐漸減小,欠采樣產生的高副瓣電平可以被濾除,PGA方法更適用于稀疏采樣條件下Φ0的估計。

因此,在算法2中,Φ0由PGA方法對稀疏采樣的匹配濾波成像結果Γ(L⊙S)進行自聚焦得到。與全采樣條件下的自聚焦過程相同,應用PGA方法時,僅需選取約10%的平均能量較大的距離6單元進行自聚焦,最大迭代次數(shù)可固定為10次,具體的算法實現(xiàn)過程可參考文獻[17]。

2.2.2G和Φ的交替重建

在算法2的步驟1中,利用給定的相位誤差Φ來重建圖像,則式(12)的優(yōu)化問題可以表示為

(14)

式(14)是一個L1范數(shù)正則化問題,相較于IST方法,FIST方法[16]具有更快的收斂速度,因此算法1結合FIST方法進行求解。在步驟1中,Ψγ,δ(·)表示軟閾值算子,δ為搜索步長,通常取δ=1來保證算法收斂,參數(shù)γ的選擇取決于成像場景的稀疏度K0。本文取γ=|Vi|(K0),其中,|Vi|(K0)是矩陣Vi所有元素幅度排序(由大到小)第K0的幅度值。

在步驟2中,利用給定的場景散射系數(shù)Gi來估計相位誤差Φi,優(yōu)化問題可以表示為

(15)

在該優(yōu)化問題中,Gi是常數(shù),則有

(16)

基于最小方差和MLE求得第m個觀測位置的誤差相位為

(17)

式中：SΩ=Ω(S)表示相位校正的回波數(shù)據(jù);[X]m表示矩陣X的第m行。

為進一步加快求解過程中的收斂速度,本文利用成像場景稀疏的先驗信息,在重建相位誤差Φ時采用了見成像場景的最小熵約束,如步驟2.1所示,具體推導過程見第2.3節(jié)。

2.3 成像場景的最小熵約束

在稀疏SAR成像中,成像場景中散射強度較弱的背景信息可以看作是噪聲,因此SAR成像場景可以看作是稀疏的。在相同稀疏度的條件下,運動誤差將導致圖像散焦,使稀疏SAR圖像的熵增加。本節(jié)將通過對重建的稀疏SAR場景進行最小熵約束來進一步提高誤差相位的重建精度,從而減少迭代次數(shù),提高算法穩(wěn)定性。

采用匹配濾波原理成像時,相位補償后的聚焦SAR圖像可以表示為

G=Fη(ΦSΩ)

(18)

根據(jù)式(18)可得,二維散射矩陣G的第m行第n列元素G(m,n)可以表示為

(19)

用于估計誤差相位的圖像熵表示為

(20)

(21)

EG對相位誤差φ(l)的偏導數(shù)表示為

(22)

由于圖像熵EG是相位誤差的函數(shù),熵值最小時,熵對每一點相位誤差的偏導數(shù)為0,即

(23)

由于|G(m,n)|2=G(m,n)*G(m,n),所以

(24)

根據(jù)式(19)可知,G(m,n)相對于φ(l)的偏導數(shù)為

(25)

將式(25)代入式(24),再將式(24)代入式(22)可得

(26)

令式(26)等于0,可得誤差相位φ(l)的最大似然估計為

(27)

式(27)是基于圖像最小熵的誤差相位最大似然估計方法,對比式(27)和式(17)可以發(fā)現(xiàn),在稀疏自聚焦框架中,求解誤差相位時,對重建的場景散射系數(shù)Gi+1乘以(1+ln|Gi+1|2)的權系數(shù)即可實現(xiàn)最小熵約束,即步驟2.1。

對于任意的SAR成像模式,在獲取稀疏采樣的原始數(shù)據(jù)L⊙S和頻域成像算子Γ(·)后,采用算法2即可實現(xiàn)稀疏自聚焦成像。

3 實測數(shù)據(jù)驗證

本節(jié)通過機載SAR的實測數(shù)據(jù)來驗證所提算法的有效性,選用文獻[14]中基于近似觀測的稀疏自聚焦方法作為參考算法,采用圖像熵來評估SAR圖像的聚焦質量,在相同的稀疏度下,熵越小圖像質量越高。除此之外,取全采樣條件下的PGA自聚焦結果作為全局最優(yōu)解,用于衡量稀疏采樣條件下算法的自聚焦效果。與參考算法相比,在每一步的迭代中,所提算法2需要額外進行一次用于最小熵約束的矩陣點乘運算和步驟2.1中用于加快收斂速度的矩陣加法運算,相對于迭代過程中近似觀測算子的運算,這些附加的運算量可以忽略。因此,所提方法和參考算法單次迭代的運算量是近似一致的,算法效率可通過收斂所需的迭代次數(shù)來體現(xiàn)。

為便于直觀比較算法的收斂速度和精度,本文算法和對比算法中的迭代次數(shù)均固定為30,為便于比較圖像的聚焦效果,所有圖片的最大亮度一致。實測數(shù)據(jù)對應的機載SAR參數(shù)如表1所示,采用線性調頻變標算子作為近似觀測算子。本文的算法模型能夠用于距離和方位二維稀疏采樣回波數(shù)據(jù)的自聚焦成像,為便于分析回波采樣率對算法收斂性的影響,實驗過程中僅在方位向上對實測數(shù)據(jù)進行抽取。方位向上的回波數(shù)據(jù)抽取率為ξa,ξa=1表示全采樣回波數(shù)據(jù)。在表1對應的實測數(shù)據(jù)中,分辨率約為0.35 m,全采樣回波數(shù)據(jù)的方位向過采樣率為1.2(即PRF與奈奎斯特采樣率的比值,其中奈奎斯特采樣率為回波信號的方位向帶寬)。因此,在稀疏采樣的條件下,回波信號相對奈奎斯特采樣的欠采樣率為1.2ξa,即當ξa=0.5時,欠采樣率為0.6。成像場景的相對稀疏度定義為α=K0/PQ,表示重建的散射點占場景所有散射點的比例。

表1 機載SAR參數(shù)Table 1 Parameters of airborne SAR

圖2給出了全采樣和方位稀疏采樣(ξa=0.5)的實測數(shù)據(jù)成像結果,其中圖2(a)為成像場景的谷歌地圖。對比圖2(b)和圖2(c)可以看出,在全采樣數(shù)據(jù)下,PGA可以實現(xiàn)理想的運動補償。而從圖2(d)中可以看出,回波數(shù)據(jù)的方位向稀疏采樣使PGA聚焦后的圖像存在大量欠采樣噪聲,嚴重降低圖像質量。

圖2 實測SAR數(shù)據(jù)成像結果Fig.2 Imging results of real SAR data

為同時補償運動誤差和抑制欠采樣噪聲,分別采用本文和文獻[14]中的稀疏自聚焦方法對該回波數(shù)據(jù)進行處理,為比較算法對場景散射點的重構能力,固定ξa為0.5,相對稀疏度α分別設置為0.1和0.5,自聚焦結果分別如圖3和圖4所示。從圖3(a)和圖3(b)中可以看出,當重建的散射點數(shù)量較少時(α=0.1),兩種算法都可以較好地估計和補償運動誤差并抑制欠采樣噪聲,圖3(c)表明,兩種算法在30次迭代內都收斂到了全局最優(yōu)解附近,而所提方法具有更快的收斂速度,僅需約7～8次迭代即可收斂。對于多數(shù)SAR圖像,僅重建少數(shù)強散射點將導致圖像信息損失。在欠采樣的條件下,先驗信息的利用程度將極大的影響散射點數(shù)量的重建能力,從圖4(a)和圖4(b)中可以看出,參考方法僅利用了成像場景稀疏的先驗信息,當重建的散射點數(shù)量大量增加時,場景的稀疏性將變弱,在欠采樣的條件下重建的成像場景將有較大誤差,進而導致參考方法不能準確地估計和補償運動誤差,致使迭代結果依舊存在明顯的方位散焦。由于利用了場景最小熵的附加先驗信息,且提供了更精確的初始誤差相位,所提方法在欠采樣條件下能夠有效重建更多的散射點。觀察圖4(c)可以發(fā)現(xiàn),隨著迭代次數(shù)的增加,參考方法逐漸收斂到誤差較大的局部最優(yōu)解,而所提方法則可以快速收斂到全局最優(yōu)解附近,具有良好的自聚焦效果,極大地減少了圖像信息損失。

圖3 相對稀疏度α=0.1時的稀疏自聚焦結果Fig.3 Sparse autofocuing results with relative sparsity α=0.1

圖4 相對稀疏度α=0.5時的稀疏自聚焦成像結果Fig.4 Sparse autofocuing results with relative sparsity α=0.5

在本文方法中,構建了一種新的稀疏自聚焦模型,使得PGA方法能夠被用于提供相位誤差的初始解,最小熵(minimum entropy, ME)約束能被用于為相位誤差的最大似然估計提供更多的先驗信息,同時所提方法采用FIST方法重構場景散射信息,進一步加快了收斂速度。為驗證PGA提供初始解和ME約束在自聚焦算法中的重要性,在本文的稀疏自聚焦模型的基礎上分別構建了不含PGA和ME約束的FIST+MLE算法,含有PGA提供初始解的PGA+FIST+MLE算法,以及含有PGA和最小熵約束的PGA+FIST+ME+MLE算法。其中FIST+MLE方法的成像效果可近似看作為文獻[14]中的算法,但事實上,受益于FIST方法的高效率,FIST+MLE的收斂速度會好于文獻[14]中的IST+MLE方法。圖5和圖6分別給出了方位采樣率和相對稀疏度對FIST+MLE,PGA+FIST+MLE和PGA+FIST+ME+MLE 3種算法收斂速度影響的實測數(shù)據(jù)分析結果?？梢钥闯?本文方法在不同的采樣率和相對稀疏度下都有很好的穩(wěn)定性。從圖5(b)和圖5(c)中可以看出,含有PGA提供初始解的兩種算法的初始圖像熵較小,這表明在相對較高的采樣率下,PGA方法可以提供有效的初始解,能夠減少算法所需的迭代次數(shù)。同時可以看出,在PGA提供初始解的基礎上,隨著方位采樣率的升高,自聚焦所需的迭代次數(shù)迅速減小。從圖5(a)中可以看出,當采樣率較低時,欠采樣噪聲將降低PGA方法對初始解的估計精度,此時自聚焦算法所需的迭代次數(shù)較多,ME約束能夠明顯提高收斂速度,并避免算法陷入誤差較大的局部最優(yōu)解。從圖6中可以看出,PGA提供初始解以及ME約束對提高算法收斂速度避免陷入誤差較大的局部最優(yōu)解的作用非常顯著,隨著場景重建散射點數(shù)量的增多,PGA+FIST+MLE方法將逐漸收斂到誤差較大的局部最優(yōu)解,而PGA+FIST+ME+MLE方法的收斂結果更接近全局最優(yōu)解。

圖5 α=0.4時,采樣率對收斂速度的影響Fig.5 Influence of the sampling rate on convergence rate with α=0.4

圖6 ξa=0.5時,相對稀疏度對收斂速度的影響Fig.6 Influence of the relative sparsity on convergence rate with ξa=0.4

4 結束語

基于近似觀測的稀疏SAR成像方法將匹配濾波類成像算法和稀疏成像方法結合在一起,極大提高了稀疏成像方法的效率,在此基礎上本文提出了一種基于近似觀測和最小熵約束的稀疏自聚焦方法,該方法將匹配濾波類成像算法中常用的PGA和最小熵等自聚焦方法結合到SAR的稀疏自聚焦成像中。實測數(shù)據(jù)成像結果表明,與常規(guī)的稀疏自聚焦方法相比,所提方法進一步利用了成像場景的最小熵約束,且為誤差相位提供了更精確的初始解,能夠在稀疏采樣的條件下,有效增大散射點的重建數(shù)量,減小迭代次數(shù),使迭代結果更接近全局最優(yōu)解。本文所構建的稀疏自聚焦模型僅能校正非空變的相位誤差,下一步將進一步研究空變相位誤差下的稀疏自聚焦方法。