基于特征增強的非均勻采樣SAR三維稀疏成像

孫 豆, 邢世其，*, 高海峰, 龐 礴, 李永禎, 王雪松

(1.國防科技大學電子科學學院電子信息系統復雜電磁環境效應國家重點實驗室, 湖南 長沙 410073;2.陜西省地質調查院, 陜西 西安 710054)

0 引 言

在合成孔徑雷達(synthetic aperture radar,SAR)技術的發展中,很多研究工作致力于獲取飛機、車輛等人造目標高分辨率的三維成像結果[1-2]。通常情況下,通過層析SAR[3-4]或全息SAR[5-6]獲取三維波數域空間中的均勻密集采樣可以得到高分辨率的三維成像結果。然而,這種數據采集方式耗時長且成本較高[7-8]。此外,在軍事監察等多種領域中,獲取均勻且密集的采樣是不切實際的。作為一種新興技術,非均勻采樣SAR為獲取高分辨率的三維成像結果帶來了極大的方便[9-11]。通過多架搭載SAR系統的小型無人機(unmanned aerial vehicle,UAV)協同的非線性飛行[12-14]獲取SAR數據是一種收集三維非均勻采樣的可行方法。由于采集的數據稀疏且無人機平臺具有良好的機動性,這種非均勻采樣方式不僅高效并且大大節省了成本,適用于軍事偵察等實際應用中。

非均勻采樣的三維成像與均勻且密集采樣的三維成像有很大的不同,傳統的基于線性濾波的三維成像方法[5,15-16]并不適用。具體地說,非均勻采樣使得俯仰向和方位向之間高度耦合,因此不能先二維成像,再估計高度信息,三維分步成像[4-5,17]不再適用。此外,由于非均勻采樣是稀疏的,使用基于傅里葉變換的成像方法得到的成像結果會存在很強旁瓣,成像質量差[18-19]。作為一種模型匹配類方法,稀疏重構通過正則化約束可以得到稀疏的結果[20-21]。因此,結合稀疏重構進行三維成像有望得到三維高分辨成像結果。然而,三維稀疏成像面臨巨大的計算壓力,使重構具有很大挑戰性和難度[22-23]。Austin提出了一種基于數據插值和快速傅里葉變換相結合的稀疏成像方法[7],解決了計算規模大的困難,使得非均勻采樣的稀疏成像成為可能。但該方法中使用局部信息進行數據插值會產生較大的數據網格誤差,進而影響成像結果的精度。此外,基于理想點散射模型的稀疏成像雖然可以得到稀疏的成像結果[24-25],但對于分布式目標,斷裂開的稀疏的點并不符合其散射連續的本質,分布式目標的成像結果稀疏不利于后續的圖像解譯和目標識別。基于屬性散射模型的稀疏重建[26-27]同時估計目標的尺寸、姿態、位置和散射幅度等多個參數,可以得到符合分布式目標散射本質的重建結果。然而,這類方法的計算復雜度很高,并且由于待估參數多,數據耦合嚴重,較難得到全局最優解。

為了充分利用非均勻采樣數據,得到高質量的成像結果,本文提出一種基于特征增強的非均勻采樣SAR三維稀疏成像方法。首先,基于無數據插值的非均勻數據,將成像問題直接建模為三維空間中聯合的稀疏重構問題。在此模型基礎上,選取候選散射中心進行字典降維和模型降維。之后,通過建立散射中心之間的三維聯系,在模型中增加三維特征增強約束項,得到最終的三維稀疏重建模型。最后,結合高斯迭代法和優化的信號處理技巧,提出了一種高效的模型求解算法。通過仿真實驗驗證了該成像方法的有效性。與其他成像方法相比,本文所提方法具有以下優勢:避免了局部插值帶來的數據誤差,因此成像精度高;增加了特征增強約束項,因此保證了分布式目標成像結果的連續性;通過字典降維和信號處理技巧減小了計算規模,因此在提高成像質量的同時保證了計算效率和現有方法相當。

1 信號模型

相對于場景中心,雷達的方位角為φ,俯仰角為θ,并發射寬帶信號。假設雷達距離場景足夠遠,使用平面波模型,接收的信號可以表示為

(1)

式中,t表示時間;c是光速;s(t)是中心頻率為fc且帶寬為BW的已知寬帶信號;*表示卷積。場景的反射率函數由g(x,y,z;φ,θ)表示,x,y,z為目標在場景中的位置。

式(1)可以理解為投影到x維的場景反射率函數的傅里葉變換。根據投影切片定理,反射率函數滿足

(2)

式中,G(kx,ky,kz)表示由接收信號r(t;φ,θ)和寬帶信號s(t)得到的波數域觀測信號。

波數域觀測的每個支撐集是波數域采樣點(kx,ky,kz)上的一個線段。以4πfc/c rad/m為中心,以雷達觀測角度(θ,φ)為方向的波數域采樣點表示為

(3)

式中,頻率f的采樣為f→fj,觀測角度的采樣為(θ,φ)→(θq,φq)。

本文的任務是根據式(2)中給出的已知波數據觀測和反射率函數之間的關系,對反射率函數進行估計。圖1給出了5架無人機協同飛行的雷達掃描軌跡,由一組隨機曲線組成。可以看出這種非均勻采樣方式形成了復雜的基線。圖2給出了圖1軌跡對應的波數域采樣情況,可以看出由觀測角度(即掃描軌跡)確定的波數域采樣是非均勻且稀疏的,這就使得成像難度變大。為了獲得非均勻且稀疏的波數域采樣的三維高分辨成像結果,下面將稀疏重構與非均勻采樣的數據特性相結合,進行稀疏成像建模。

圖1 UAV協同飛行的掃描軌跡

圖2 UAV協同飛行的波數據采樣

2 稀疏成像建模

一般情況下,SAR圖像的場景是稀疏的,特別是對于人造目標而言,其由少量的強散射體組成[17]。考慮到場景的稀疏性,將成像問題建模為稀疏重構問題有望提高成像結果的質量。為了高效率地得到高分辨率高質量的成像結果,下面分3步進行稀疏成像建模。

首先建立直接三維稀疏成像模型,避免將非均勻數據插值為均勻數據帶來的誤差。之后,進行模型降維處理,解決直接三維稀疏成像模型求解中面臨的計算量大、存儲壓力大的問題。為了保證分布式目標成像結果連續,最后在模型中增加三維特征增強約束項。

2.1 直接三維稀疏成像模型

在圖像重構空間中,我們定義一組N個候選散射中心的位置:

(4)

通常,這些位置從均勻的矩形網格中選擇。基于這些位置和波數域采樣,M×N維的字典矩陣表示為

A=[e-j(kx,mxn+ky,myn+kz,mzn)]m,n

(5)

式中,m=jq表示波數域中M個觀測的索引;n表示C中N個候選散射中心的索引。

為了避免數據插值帶來的誤差,這里直接利用式(5)中的字典矩陣建立直接三維稀疏成像模型。

結合式(5),將式(2)寫成矩陣形式,得到

b=Aβ

(6)

式中,向量β是要重建的三維圖像;向量b是M維的波數域觀測。

三維稀疏成像等價于解決下面的稀疏重構問題:

(7)

式(7)中的數學模型不是一個凸優化問題,而是一個NP難問題。由于lp范數具有產生稀疏和更精確的解的能力,且非常類似于l0范數[28],這里使用lp范數將式(7)中的模型松弛為下面的優化問題:

(8)

2.2 模型降維

在觀測場景中,我們感興趣的只是目標區域。并且通常情況下,目標區域以外的大部分區域的成像結果強度都很小甚至為零,特別是由強孤立散射體組成的目標。因此,只需對目標區域進行成像,進而獲取目標的成像結果。基于此,下面通過選取候選散射中心進行字典降維,同時降低稀疏成像模型的規模。

(9)

圖3 散射中心選取示意圖

(10)

候選散射中心選取后,字典矩陣更新為

(11)

由于選取后的候選散射中心數目為N1,字典矩陣A′的維數降為M×N1。同時,稀疏成像模型更新為

(12)

經過散射中心選取,式(12)中模型的規模相比于式(8)縮小了很多,計算量和存儲壓力也隨之變小。

2.3 三維特征增強

稀疏成像使用式(2)所示的點散射模型,成像結果通常是由幾個稀疏的點組成的。實際上,分布式目標的電磁散射模型并不是簡單的點散射模型。因此,在對分布式目標稀疏成像時,其成像結果會出現斷裂的情況,由一些斷裂開的點組成。這樣的成像結果不利于后續的目標識別和圖像解譯。

為了同時保持稀疏成像對于旁瓣的有效抑制和對點散射目標的高分辨成像,并增強分布式目標,得到分布式目標連續且平滑的成像結果,我們給稀疏成像模型中增加三維特征增強約束項,即通過建立相鄰散射中心散射強度之間的關聯,約束能量分布趨于均勻。模型改進為

(13)

在三維空間中,目標散射中心的相鄰散射中心有13個,分別位于直角坐標系的3個方向上和斜對角線的10個方向上。通過仿真發現,選取直角坐標系的3個方向上的相鄰散射中心足夠實現區域的平滑,因此三維特征增強只對這3個方向進行約束,D矩陣表示為

(14)

三維特征增強約束項表示為

(15)

式中,Dx|β′|,Dy|β′|,Dz|β′|表示每個目標散射中心與圖4所示的x,y,z3個方向上的相鄰散射中心之間的強度變化情況。對于選取后的候選散射中心β′內某一序號為q的目標散射中心,其在x,y,z3個方向上的前向相鄰散射中心的序號分別為qx,qy,qz,如圖4所示。基于此,則Dx,Dy,Dz為

圖4 三維特征增強示意圖

Dx=Λ+X,Dy=Λ+Y,Dz=Λ+Z

(16)

式中,

Λ=diag(-1)

(17)

(18)

三維特征增強項建立了相鄰散射中心幅度變化的約束。通過調整正則化參數λ2,可以約束相鄰散射中心之間幅度變化程度,達到對分布式目標成像結果進行增強的目的。

3 稀疏成像算法

經過上面的推導,最終要求解的模型為式(13)。根據式(13),定義代價函數J(β′):

(19)

當代價函數J(β′)取最小值時,就得到了稀疏重構的結果。受Cetin的處理思想[29]啟發,首先計算J(β′)對β′的偏導:

(20)

式中,Λ1(β′)=diag(|(β′)i|p-2)是N1維的對角矩陣,H(β′)為

H(β′)=ΦH(β′)DxHΛx(β′)DxΦ(β′)+

ΦH(β′)DyHΛy(β′)DyΦ(β′)+ΦH(β′)DzHΛz(β′)DzΦ(β′)

(21)

式中,

Λx(β′)=diag(|(Dx|β′|)i|p-2)

(22)

Λy(β′)=diag(|(Dy|β′|)i|p-2)

(23)

Λz(β′)=diag(|(Dz|β′|)i|p-2)

(24)

Φ(β′)=diag(exp(-jφ[(β′)i]))

(25)

式中,φ[(β′)i]表示(β′)i的相位。

(26)

(27)

(28)

實際上,式(27)的計算規模很大。為了降低計算量和存儲壓力,進而高效地進行迭代計算,這里采用一種優化的信號處理技巧。

算法1總結了所提出的稀疏成像算法的具體步驟。

算法1 基于特征增強的非均勻采樣SAR三維稀疏成像算法1.使用3D-NUFFT進行成像,得到初始成像結果β⌒;2.根據式(9)和初始成像結果β⌒,確定目標區域C′;3.根據目標區域C′以及序號為q的目標散射中心在x,y,z3個方向上的前向相鄰散射中心的序號qx,qy,qz,確定Dx,Dy,Dz;4.根據目標區域C′、波數域觀測G(kx,ky,kz)和波數域采樣kj,qx,kj,qy,kj,qz,計算并存儲A′Hb和A′HA′;5.根據式(27)迭代計算β～′k+1,當‖β～′k+1-β～′k‖22/‖β～′k‖22<τ時,得到解β～′=β～′k+1;6.根據式(28),由β～′拼接得到完整的三維成像結果β～。

4 仿真實驗與分析

為了評價本文提出成像方法的性能,本節基于電磁仿真數據開展仿真實驗,并將本文方法和3D-NUFFT,Austin方法進行對比。實驗在一臺配備Intel Core I5-6 500 CPU和12 GB RAM的計算機上的Matlab R2016b上進行。

4.1 成像性能分析

選取單個三面角作為仿真目標,從抗噪性、采樣分布依賴性、成像結果分辨率、計算時間這4個方面分析本文成像方法的性能。回波數據由電磁仿真軟件生成,雷達中心頻率為10 GHz,帶寬為2 GHz,三面角放置在場景中心。

圖5給出了3種不同采樣方式沿俯仰角和方位角的掃描軌跡。三面角的散射回波數據依據這3種采樣方式獲取。雖然3種掃描軌跡不同,但這3個軌跡的方位角和俯仰角范圍分別都是[86°,96°]和[18°,42°]。因此,比較它們的成像結果是合理的。這3種不同采樣方式的波數據采樣情況如圖6所示。可以看出每種采樣方式的波數域采樣都是非均勻且稀疏的,但波數域采樣的分布情況各不相同,方式1是最稀疏的,方式3最密集。

圖5 3種采樣方式的掃描軌跡

圖6 3種采樣方式的波數據采樣

依據圖5中3種采樣方式下的回波數據,使用本文方法、3D-NUFFT和Austin方法對三面角進行成像得到如圖7所示的結果。圖7中,設置顯示成像結果的幅度閾值為-20 dB,每個子圖的中間青色部分是目標的三維成像結果,周圍灰色部分是X-Y、Y-Z和X-Z方向的二維投影結果。從圖7可以看出,不同采樣方式下的3D-NUFFT成像結果存在差異,且均不理想,旁瓣很高。而對于Austin方法和本文方法,每種采樣方式下的成像結果相似,都只有一個橢球,完全去除了多余的旁瓣,實現了單個三面角的稀疏三維成像,且這兩種方法不受采樣方式的影響,采樣分布依賴性低。為了更直觀地對比這3種方法,圖8給出了成像結果在距離向、方位向和俯仰向的切片圖。從圖中可以看出,每種采樣方式下,3D-NUFFT方法得到的結果的主瓣最寬,且明顯可以看見旁瓣。因此,該方法分辨率差,并且旁瓣抑制能力差。對于Austin方法和本文方法,其結果均具有顯著的稀疏特性,主瓣較窄,并且都抑制掉了旁瓣,兩種方法具有一定的旁瓣抑制能力。然而,相比較于Austin方法,本文方法結果的主瓣更窄,因此具有更好的分辨率。

圖7 三面角的成像結果

圖8 三面角成像結果的切片圖

圖9給出了不同軌跡和不同信噪比下本文方法和Austin方法的成像結果分辨率對比結果。這里參考文獻[30]中對分辨率的定義,取三維成像結果中大于-6 dB的橢球體積作為成像分辨率值。由于 3D-NUFFT的成像結果較差,這里沒有選取此方法進行對比分析。從-40 dB到10 dB,在不同的SNR下進行了50次試驗。從圖9可以看出,當信噪比大于-30 dB時,所有曲線趨于穩定,分辨率值基本保持不變,且本文方法的分辨率總是優于Austin方法。此外,不同采樣方式下本文方法的分辨率值基本相同,而Austin方法在3種采樣方式下的分辨率值略有不同。當信噪比小于-30 dB時,隨著信噪比的降低,所有曲線都開始上升,兩種方法的成像質量越來越差。

圖9 分辨率對比結果

根據以上結果可得,Austin方法和本文方法基本上具有相同的抗噪性能,且均不受采樣方式的影響,對采樣分布的依賴性低,但是在成像結果分辨率上,本文方法優于Austin方法。

表1給出了3D-NUFFT、Austin和本文方法計算時間的對比結果。可以看出,3種不同方法的計算時間差距不大,計算效率基本處于一個水平。3D-NUFFT所需的計算時間最少,而Austin方法由于包含了3D-NUFFT中的插值步驟,所需的時間比3D-NUFFT略多。本文方法直接對非均勻數據進行稀疏重構,因此相比于Austin方法,算法的收斂需要更多的時間。

表1 不同成像方法的計算時間比較

4.2 典型目標成像結果

為進一步驗證本文成像方法的有效性,下面選取兩個典型目標開展仿真實驗,并比較本文方法成像結果和其他成像方法的不同。

為了對比本文方法和其他方法在成像結果精度上的差異,首先對一個飛機目標進行仿真實驗。圖10給出了飛機的三維CAD模型,由13個三面角組成。

圖10 飛機的三維CAD模型

飛機的長度和寬度分別為0.7 m和0.6 m,在場景中的傾角為30°,因此成像場景的高度為0.35 m。回波數據由電磁仿真軟件生成,雷達中心頻率為10 GHz,帶寬為2 GHz,雷達沿俯仰角和方位角的掃描軌跡與圖5(a)所示的采樣方式1一致。

圖11給出了使用不同成像方法得到的飛機三維成像結果。

圖11 飛機的成像結果

圖11(a)中設置顯示成像結果的幅度閾值為-10 dB,圖11(b)和圖11(c)中設置顯示成像結果的幅度閾值為-20 dB。如圖11(a)所示,由3D-NUFFT得到的成像結果質量較差,具有高的旁瓣,成像結果中看不到飛機的輪廓。圖11(b)中Austin方法得到的成像結果可以看出飛機的大致輪廓,但是成像結果的質量一般。雖然Austin方法減少了大部分的旁瓣,但成像結果中仍有一些旁瓣沒有被去掉,對旁瓣的抑制能力不夠,圖11(b)中用橙色圈出了這些旁瓣。此外,Austin方法成像結果的精度也一般,沒有得到機翼上的兩個三面角的成像結果,圖11(b)中用紅色圈出了這一缺失的部分。實際上,由于Austin方法中包含基于局部數據信息的插值步驟,數據插值帶來的額外誤差影響了成像效果,造成了圖11(b)中成像結果精度不高。圖11(c)給出了本文方法的成像結果,可以看出成像結果準確度很高,與真實場景吻合,飛機的輪廓非常清晰,旁瓣被完全去除掉了,且13個三面角在圖11(c)中都被顯示出來了。因為本文方法不需要數據插值,保證了數據的真實性,因此相比其他成像方法,本文方法不僅對旁瓣抑制能力強,并且獲得的成像結果精度更高且更稀疏。

最后,為了測試本文方法對分布式目標的成像效果,這里使用Slicy模型進行仿真實驗,圖12給出了其三維CAD模型。Slicy模型由三面角、二面角、圓柱等組成,其中的二面角、圓柱是典型的分布式目標。回波數據由電磁仿真軟件生成,雷達中心頻率為10 GHz,帶寬為2 GHz,雷達沿俯仰角和方位角的掃描軌跡與圖5(a)所示的采樣方式1一致。

圖12 Slicy的三維CAD模型

圖13給出了使用不同成像方法得到的Slicy模型三維成像結果。為了驗證本文方法對分布式目標成像的有效性,圖13中還增加了本文方法在沒有三維特征增強約束項時的成像結果進行對比。圖13(a)中設置顯示成像結果的幅度閾值為-10 dB,圖13(b)～圖13(d)中設置顯示成像結果的幅度閾值為-20 dB。如圖13(a)所示,通過3D-NUFFT獲得的成像結果較差,并且具有很高的旁瓣。圖13(b)中Austin方法去除了大部分旁瓣,但成像結果中仍有一些旁瓣沒有被去掉,對旁瓣的抑制能力不夠,圖13(b)中用橙色圈出了這些旁瓣。此外,Austin方法的成像結果的精度較低,二面角的成像結果出現斷裂的情況,且有些部分沒有被重構出來,圖13(b)中用紅色圈出了這些部分。圖13(c)給出了本文方法在沒有三維特征增強約束項時的成像結果,該結果稀疏,完全去除掉了旁瓣,對旁瓣抑制能力強,并且Slicy模型每個部分都被重構出來了,但是圖13(c)中用紫色圈出的分布式目標的成像結果是斷裂的,由多個離散點組成,并不符合分布式目標連續散射的機理。本文方法的成像結果如圖13(d)所示,由于增加了特征增強約束項,該結果不僅保持了稀疏成像對旁瓣的抑制,完全去除了旁瓣,并且得到了連續的分布式目標成像結果,避免了稀疏的點組成分布式目標的情況。圖13(d)中Slicy模型的每個部分都被重構出來,和真實場景匹配且成像結果連續,這就驗證了本文方法對分布式目標成像的有效性。

圖13 Slicy的成像結果

5 結 論

本文提出了一種基于特征增強的非均勻采樣SAR三維稀疏成像方法。該方法直接對非均勻采樣數據進行稀疏成像建模,避免了局部插值對結果帶來的不良影響,保證了數據的真實性。由于模型中增加了特征增強約束項,因此該方法不僅保持了稀疏成像對旁瓣的有效抑制,且保證了分布式目標成像結果連續。字典降維和信號處理技巧減小了計算規模,因此在提高成像質量的同時,該方法保證了算法計算效率和其他現有方法相當。從實驗結果來看,本文方法在抗噪性、采樣分布依賴性、計算效率這3個方面和其他方法的性能相當。在旁瓣抑制能力、成像結果分辨率和精度方面,本文方法明顯優于其他方法,且本文方法可以得到分布式目標連續的稀疏成像結果。