杜 邦 仇曉蘭 張 柘 雷 斌 丁赤飚

①(中國科學(xué)院空間信息處理與應(yīng)用系統(tǒng)技術(shù)重點(diǎn)實(shí)驗(yàn)室 北京 100190)

②(中國科學(xué)院空天信息創(chuàng)新研究院 北京 100190)

③(中國科學(xué)院大學(xué) 北京 100049)

④(蘇州空天信息研究院 蘇州 215123)

1 引言

層析合成孔徑雷達(dá)(Synthetic Aperture Radar Tomography,TomoSAR)作為干涉合成孔徑雷達(dá)技術(shù)(Synthetic Aperture Radar Interfemetry,InSAR)的擴(kuò)展，通過獲取斜高向上不同位置的觀測數(shù)據(jù)形成斜高向的合成孔徑，具有高程向的分辨能力[1,2]。TomoSAR能夠解決由于SAR成像幾何帶來的疊掩問題，實(shí)現(xiàn)真正意義上的三維分辨。近年來，隨著星載SAR、機(jī)載SAR系統(tǒng)技術(shù)水平的提升，分辨率、信噪比等指標(biāo)不斷提高，通過星載或機(jī)載SAR層析觀測，已可實(shí)現(xiàn)觀測區(qū)域高精度三維點(diǎn)云的生成，甚至獲得建筑形變等更高維度的信息[3–5]。

TomoSAR在斜高向的分辨，本質(zhì)上是信號處理領(lǐng)域常見的譜估計(jì)問題，近年來已經(jīng)有許多相關(guān)研究成果。受到系統(tǒng)復(fù)雜度或觀測周期的制約，TomoSAR的基線數(shù)目一般是比較有限的，并且基線分布大多是非均勻的，這些問題限制了傳統(tǒng)譜估計(jì)算法的應(yīng)用。由于城區(qū)人造建筑設(shè)施通常具有較強(qiáng)的散射強(qiáng)度，并且在斜高向上的分布往往是稀疏的，因此，壓縮感知(Compressive Sensing,CS)方法[6,7]成為TomoSAR第三維重建時經(jīng)常采用的有效方法，它能夠突破瑞利分辨率的限制，有效降低對觀測次數(shù)的需求。其中，基于L1范數(shù)最小化的方法[8–14]具有很高的分辨能力和魯棒性，已成為Tomo-SAR三維成像的經(jīng)典方法。

然而，為了應(yīng)用壓縮感知理論，需要對連續(xù)的高程向進(jìn)行離散化處理，同時假定目標(biāo)正好位于離散化后的網(wǎng)格上。實(shí)際上，無論劃分多精細(xì)的網(wǎng)格，都無法保證目標(biāo)正好位于所劃分的網(wǎng)格上，也即必然存在離網(wǎng)格(Off-grid)的效應(yīng)，從而引起壓縮感知中的基失配(Basis Mismatch)的問題。文獻(xiàn)[15]考慮了觀測矩陣為傅里葉矩陣的情況，分析了基失配對壓縮感知算法恢復(fù)信號的影響，結(jié)果表明Off-grid目標(biāo)的能量會泄露到其鄰近網(wǎng)格上，導(dǎo)致重構(gòu)性能的下降。針對這個問題，一種常見的思路是劃分足夠密集的網(wǎng)格，但這一方面會導(dǎo)致計(jì)算復(fù)雜度的增加；另一方面，過于密集的網(wǎng)格[16–18]也會使得觀測矩陣中包含大量相干性較高的字典集，嚴(yán)重情況下會導(dǎo)致壓縮感知算法的失效。

近年來，也有文獻(xiàn)[19–21]針對壓縮感知理論的缺陷，提出了最小原子范數(shù)優(yōu)化(Atomic Norm Minimization,ANM)方法重構(gòu)連續(xù)域的值，徹底避免了網(wǎng)格效應(yīng)，這一類方法也被稱為無網(wǎng)格壓縮感知(Gridless Compressive Sensing,GCS)。然而ANM采用了全新的信號框架，存在運(yùn)算量大的挑戰(zhàn)，目前在TomoSAR實(shí)際數(shù)據(jù)處理中還少有應(yīng)用。

另一種降低Off-grid對重構(gòu)性能影響的思路是將目標(biāo)偏離網(wǎng)格的誤差納入待估計(jì)的參數(shù)中[22,23]。例如文獻(xiàn)[24]和文獻(xiàn)[25]將網(wǎng)格偏離誤差建模為一個加性擾動。其中文獻(xiàn)[24]方法采用貝葉斯壓縮感知的框架，同樣面臨運(yùn)算量較大，在TomoSAR實(shí)際數(shù)據(jù)處理應(yīng)用較少的問題。而文獻(xiàn)[25]中提出的方法可以與TomoSAR中常采用的L1范數(shù)最小化進(jìn)行結(jié)合[24]，文章在單站雷達(dá)一維聚焦和SAR二維成像上驗(yàn)證了所提出算法的有效性。然而，TomoSAR三維成像觀測矩陣與上述一維和二維成像的理想傅里葉觀測矩陣存在較大差別，不能直接應(yīng)用。

本文首先考慮了在TomoSAR的基失配問題，首先推導(dǎo)了TomoSAR中針對Off-grid目標(biāo)的1階擾動模型，并進(jìn)一步將結(jié)合L1范數(shù)最小化的求解方法應(yīng)用到求解Off-grid TomoSAR模型中，最后用仿真數(shù)據(jù)和實(shí)測機(jī)載數(shù)據(jù)對所提方法進(jìn)行了驗(yàn)證，證明了方法的正確性和有效性。

本文后續(xù)內(nèi)容安排如下：第2節(jié)給出了傳統(tǒng)的TomoSAR成像模型及求解方法；第3節(jié)推導(dǎo)并分析了Off-grid目標(biāo)對三維成像的影響，給出了適用于Off-grid目標(biāo)的TomoSAR三維成像模型，以及結(jié)合L1最小范數(shù)優(yōu)化的求解方法；第4節(jié)進(jìn)行實(shí)驗(yàn)，通過點(diǎn)目標(biāo)仿真和中國科學(xué)院空天信息創(chuàng)新研究院的8通道機(jī)載陣列干涉SAR數(shù)據(jù)驗(yàn)證了所提方法的有效性；第5節(jié)進(jìn)行總結(jié)與展望。

2 TomoSAR信號模型

圖1給出了TomoSAR三維成像坐標(biāo)系下距離-斜高平面內(nèi)的成像幾何關(guān)系。假設(shè)對某個區(qū)域有N個相干的觀測，相當(dāng)于在斜高向上形成了等效合成孔徑，采用連續(xù)信號模型可以表示為

圖1 TomoSAR成像幾何關(guān)系Fig.1 The geometry of TomoSAR imaging

其中，gi表示第i個觀測采集到的復(fù)數(shù)觀測值，γ(s)是垂直斜距的高程向(也即斜高向)上散射體的散射系數(shù)分布，ξi=2bi/(λr)是 由等效基線bi，雷達(dá)波長λ和 斜距r決定的等效空間頻率。對上述連續(xù)信號模型進(jìn)行離散化后，表達(dá)為矩陣形式如式(2)。

其中，A是N×L大 小的觀測矩陣，L為s軸上劃分的網(wǎng)格數(shù)量，矩陣中第i行 第k列的元素為Aik=exp(-j2πξisk) ,sk是第k個預(yù)先劃分的網(wǎng)格位置，ε是噪聲向量。

在上述模型下，TomoSAR三維成像的目標(biāo)是通過求解式(2)獲得目標(biāo)的高程s及其后向散射系數(shù)γ(s)。L1范數(shù)最小化的方法能夠較好地處理該問題，其核心思想在于求解以下約束以獲得盡可能稀疏的估計(jì)。

其中，μ是拉格朗日系數(shù)[26]，其與觀測數(shù)量及噪聲水平相關(guān)。式(3)可以進(jìn)一步轉(zhuǎn)化為凸優(yōu)化的問題，利用現(xiàn)有的工具箱如CVX,L1-LS,SPG-L1等進(jìn)行求解。

3 Off-gird TomoSAR

3.1 一階差分?jǐn)_動模型

上述離散網(wǎng)格模型下，如目標(biāo)恰好均位于離散網(wǎng)格上，則可以得到很好的求解結(jié)果，但如果目標(biāo)不恰好位于網(wǎng)格上，則會存在誤差，下面我們分析該誤差。傳統(tǒng)處理通常將網(wǎng)格均勻劃分，網(wǎng)格間距可以用Δs=smax/L進(jìn) 行表示，其中smax表示斜高向上的最大范圍，L代表被劃分的網(wǎng)格數(shù)?，F(xiàn)假設(shè)目標(biāo)與其最近的第k個網(wǎng)格存在偏差 Δk，也即目標(biāo)的真實(shí)位置為sk=(k+Δk)Δs，使用泰勒展開不難得到

其中，o(Δk2)包含了高階項(xiàng)。整理式(2)和式(4)，本文用組稀疏的形式表示Off-grid目標(biāo)的TomoSAR模型。

其中，A′可以看成是網(wǎng)格偏差引起的擾動矩陣，矩陣A′和 向量γ′的元素可分別表示為

新的噪聲向量ε′中除了原始噪聲ε之外，還包含了忽略的高階展開項(xiàng)。通過求解式(5)，不難進(jìn)一步通過式(7)求解出k。

3.2 局部閾值優(yōu)化

式(5)通過引入擾動項(xiàng)實(shí)現(xiàn)了比式(2)更加準(zhǔn)確的TomoSAR成像建模，然而，由于網(wǎng)格依然存在，如文獻(xiàn)[14]中分析指出，對于Off-grid的目標(biāo)，求解式(5)恢復(fù)的信號與真實(shí)信號相比，能量仍會以包絡(luò)的形式泄露到鄰近的網(wǎng)格之上，距離目標(biāo)真實(shí)位置近的網(wǎng)格獲得能量多，而距離目標(biāo)真實(shí)位置較遠(yuǎn)的網(wǎng)格獲得能量少。如圖2所示，實(shí)線表示真實(shí)位置上的目標(biāo)信號，虛線表示恢復(fù)出的信號，這會對后續(xù)的處理步驟產(chǎn)生影響。為了解決這個問題，本文引入局部閾值優(yōu)化(Local Optimizing Thresholding,LOT)[27]方法來選擇接近真實(shí)峰值的網(wǎng)格，判斷依據(jù)為

圖2 高程向目標(biāo)重建結(jié)果的簡單示意圖Fig.2 Sketch map of the reconstruct results in the elevation direction

通過式(8)可以得到由αt構(gòu)成的支撐集。

3.3 算法流程

基于上述模型構(gòu)建和支撐集確定方法，本文以基追蹤降噪(Basis Pursuit De-Nosing,BPDN)進(jìn)行L1范數(shù)最小化求解為例，結(jié)合局部閾值優(yōu)化方法，給出整個算法(Basis Pursuit Local Optimizing Thresholding,BPLOT)的流程，如表1所示。

表1 BPLOT的流程Tab.1 The process of BPLOT

4 實(shí)驗(yàn)結(jié)果與分析

為了證明本文提出的Off-grid TomoSAR模型和所提求解方法的有效性，下面分別進(jìn)行點(diǎn)目標(biāo)仿真實(shí)驗(yàn)和實(shí)際數(shù)據(jù)三維成像實(shí)驗(yàn)。文中所采用的實(shí)際數(shù)據(jù)為中國科學(xué)院空天信息創(chuàng)新研究院機(jī)載陣列干涉SAR實(shí)驗(yàn)系統(tǒng)獲得的山西運(yùn)城地區(qū)的8通道陣列干涉SAR數(shù)據(jù)[28]。

4.1 點(diǎn)目標(biāo)仿真

首先本文進(jìn)行理想點(diǎn)目標(biāo)仿真實(shí)驗(yàn)，仿真參數(shù)參考了實(shí)測機(jī)載陣列干涉SAR飛行實(shí)驗(yàn)參數(shù)，見表2。在仿真中，本文將斜高向均勻劃分為128個網(wǎng)格，同時隨機(jī)生成目標(biāo)位置，目標(biāo)散射系數(shù)幅度均歸一化為1，相位設(shè)置為相互正交，為了方便顯示，將高程向s歸一化至( 0,1]之間。此處以城市地區(qū)常見的疊掩目標(biāo)稀疏度K=2和K=3的情況為例[29,30]。圖3展示了采用本文方法和傳統(tǒng)方法的點(diǎn)目標(biāo)仿真實(shí)驗(yàn)結(jié)果，藍(lán)色表示真值，紅色表示采用L1范數(shù)最小化求解式(3)的結(jié)果，黃色表示采用本文方法的結(jié)果。表3給出了點(diǎn)目標(biāo)位置的真值以及位置估計(jì)結(jié)果，其中L1范數(shù)最小化方法選取了離真值最近的估計(jì)結(jié)果(本文后續(xù)均采用該評價方法)。由仿真結(jié)果可見，采用本文方法得到了更加準(zhǔn)確的目標(biāo)位置估計(jì)結(jié)果，并且目標(biāo)的幅度也得到了更好的重建。

表3 圖3對應(yīng)的位置計(jì)算結(jié)果Tab.3 Point calculation results corresponding to Fig.3

圖3 點(diǎn)目標(biāo)仿真結(jié)果Fig.3 Point target simulation results

表2 仿真參數(shù)Tab.2 Simulation parameters

隨后，本文通過蒙特卡洛實(shí)驗(yàn)評估BPLOT在Off-grid TomoSAR中的統(tǒng)計(jì)性能。在每次隨機(jī)實(shí)驗(yàn)中，首先采用估計(jì)位置與真實(shí)位置的曼哈頓距離來比較本文方法與傳統(tǒng)方法對目標(biāo)位置的估計(jì)精度，曼哈頓距離的標(biāo)準(zhǔn)計(jì)算公式為

其中，sk表 示第k個散射點(diǎn)的真實(shí)位置，表示該散射點(diǎn)的估計(jì)位置。當(dāng)時，本文將該次隨機(jī)實(shí)驗(yàn)作為本文方法優(yōu)于傳統(tǒng)方法的一次抽樣，最后通過大量獨(dú)立重復(fù)實(shí)驗(yàn)，統(tǒng)計(jì)在不同信噪比下本文方法優(yōu)于傳統(tǒng)方法的概率。并且，為了更加嚴(yán)格地進(jìn)行比較，當(dāng)K＞1時，本文還分別計(jì)算和比較了每個散射點(diǎn)的位置估計(jì)精度，統(tǒng)計(jì)了本文方法每個散射點(diǎn)位置估計(jì)精度均優(yōu)于傳統(tǒng)方法的概率。

本文分別對城市地區(qū)常見的稀疏度K=1～3的情況進(jìn)行了測試，在仿真中添加信噪比為5 dB至20 dB的加性噪聲，進(jìn)行了1000次的獨(dú)立重復(fù)實(shí)驗(yàn)，得到的結(jié)果如圖4和圖5所示。圖4給出了每個散射點(diǎn)估計(jì)位置與真值位置的偏差的統(tǒng)計(jì)平均值隨信噪比的變化，圖中虛線為本文方法結(jié)果，實(shí)線為傳統(tǒng)L1范數(shù)最小化的結(jié)果，可見在不同稀疏度和不同信噪比下，本文方法對每個散射點(diǎn)都能得到更加準(zhǔn)確的位置估計(jì)結(jié)果。當(dāng)信噪比大于10 dB時，本文方法對目標(biāo)位置估計(jì)精度能達(dá)到0.2個網(wǎng)格間距以內(nèi)。圖5給出了本文方法優(yōu)于傳統(tǒng)方法的概率隨信噪比的變化曲線，可以看出：如以總曼哈頓距離來比較，當(dāng)信噪比為5 dB時，本文方法有約70%的概率優(yōu)于傳統(tǒng)方法；如當(dāng)每個散射點(diǎn)的估計(jì)精度均優(yōu)于傳統(tǒng)方法時才統(tǒng)計(jì)為一次更精確的估計(jì)，則當(dāng)信噪比大于13 dB時，本文方法有50%以上的概率能夠更精確地估計(jì)出每個目標(biāo)的位置。上述結(jié)果驗(yàn)證了本文方法相對于傳統(tǒng)On-grid假設(shè)L1范數(shù)最小化方法的優(yōu)越性。

圖4 散射點(diǎn)位置估計(jì)值與真值距離的統(tǒng)計(jì)平均值Fig.4 The average distance between the estimated scattering point position and the true value

圖5 本文方法優(yōu)于傳統(tǒng)方法的概率與SNR的關(guān)系Fig.5 Probability of the proposed algorithm providing better estimation versus SNR

需要說明的是，受到觀測數(shù)量只有8個的限制，本文方法信噪比較低且稀疏度較大(如K=3)時，估計(jì)精度也會受到較大影響。為了考察本文方法在TomoSAR三維成像中的成功概率與觀測數(shù)量的關(guān)系，本文固定SNR=10 dB，在觀測數(shù)量為8～20的條件下分別進(jìn)行蒙特卡洛實(shí)驗(yàn)，仿真時基線總長度不變，基線間隔隨著觀測數(shù)量的增多而均勻變小。本文將每個散射點(diǎn)DBPLOT＜DOn-gridL1且DBPLOT小于1/8網(wǎng)格間距的隨機(jī)實(shí)驗(yàn)作為一次“成功”且本文方法獲得更高精度的概率抽樣，從而統(tǒng)計(jì)本文方法更精確重建的概率。實(shí)驗(yàn)結(jié)果如圖6所示，可見在信噪比為10 dB時，對于K=2的情況，觀測次數(shù)為8時本文方法就有50%以上的概率更精確地實(shí)現(xiàn)“成功”重建；K=3的情況，觀測次數(shù)大于11時，本文方法有50%以上的概率更精確地實(shí)現(xiàn)“成功”重建。上述仿真結(jié)果為本文方法的工程應(yīng)用提供了參考。

圖6 獲得更精確估計(jì)的概率與觀測數(shù)量的關(guān)系Fig.6 Probability of the proposed algorithm providing better estimation versus number of acquisitions

為進(jìn)一步評估本方法在重建目標(biāo)幅度和相位方面的優(yōu)越性，本文比較了BPLOT和傳統(tǒng)L1范數(shù)最小化得到的散射系數(shù)與仿真設(shè)置目標(biāo)散射系數(shù)的均方誤差(Root Mean Square Error,RMSE)：

其中，Nt代表了進(jìn)行獨(dú)立仿真實(shí)驗(yàn)的次數(shù)。在該仿真實(shí)驗(yàn)中，我們固定觀測數(shù)量為8，并且只統(tǒng)計(jì)本方法每個散射點(diǎn)的位置估計(jì)精度均優(yōu)于傳統(tǒng)方法時的散射系數(shù)估計(jì)結(jié)果，從而與圖4和圖5結(jié)果相對應(yīng)。實(shí)驗(yàn)結(jié)果如圖7所示，方形標(biāo)記代表K=1，菱形標(biāo)記代表K=2，叉號標(biāo)記代表K=3的情況。仿真結(jié)果可見，本文方法對于Off-grid目標(biāo)可以更精確地恢復(fù)目標(biāo)散射系數(shù)。

圖7 RMSE的計(jì)算結(jié)果Fig.7 RMSE results

4.2 機(jī)載飛行數(shù)據(jù)

為驗(yàn)證本方法的實(shí)用性，本文將其應(yīng)用于實(shí)測機(jī)載數(shù)據(jù)的三維成像處理中。圖8(a)給出了整個實(shí)驗(yàn)場景的谷歌地球影像，圖8(b)為機(jī)載SAR獲取的該區(qū)域的二維圖像，場景中的主要目標(biāo)為人造建筑。

在整個觀測場景中，本文選擇圖8中紅框標(biāo)記的最后一排建筑物的一個局部(圖9紅框標(biāo)記所示)作為實(shí)驗(yàn)對象，采用本文方法和L1范數(shù)最小化方法進(jìn)行三維成像處理。圖10顯示了本文方法對所選建筑物的 3-D 點(diǎn)云重建結(jié)果，實(shí)地考察該樓的高度在53 m左右[31]，重建結(jié)果與真值比較接近。為了清晰顯示本文方法和傳統(tǒng)方法的結(jié)果差異，我們截取一個方位向的距離-斜高向切片，如圖11、圖12所示，藍(lán)色點(diǎn)為直接使用L1范數(shù)最小化方法的結(jié)果，紅色點(diǎn)為采用本文方法的結(jié)果?？梢姡苯硬捎肔1范數(shù)最小化時，很難保證目標(biāo)位于預(yù)設(shè)網(wǎng)格之上，導(dǎo)致原本只有1個目標(biāo)的周圍會出現(xiàn)錯誤的多個估計(jì)；當(dāng)采用本文方法后，這些“虛假”目標(biāo)消失了，并且散射點(diǎn)的位置分布在視覺上更加自然，驗(yàn)證了本文方法的有效性。

圖8 從谷歌地球獲取的實(shí)驗(yàn)地區(qū)于2015年的光學(xué)影像和對應(yīng)的SAR圖像Fig.8 Optical image of the experiment area obtained from Google Earth in 2015 and SAR image of the corresponding area

圖9 圖8(a)中最后一排建筑光學(xué)影像Fig.9 Optical image of the last row of buildings Fig.8(a)

圖10 圖9中紅框建筑的三維重建結(jié)果Fig.10 TomoSAR results of the building marked in Fig.9

圖11 本文方法與傳統(tǒng)L1范數(shù)最小化的方法對圖8中紅線標(biāo)記的方位向切片的成像結(jié)果比較Fig.11 Comparison of our method and the conventional process of L1 minimization in the slant range plane for the azimuth-bin marked in Fig.8

圖12 圖11中區(qū)域的放大圖Fig.12 Zoom in of the circled area in Fig.11

5 結(jié)束語

本文探討了現(xiàn)有基于壓縮感知的TomoSAR處理算法中很少考慮的由Off-grid目標(biāo)所引起的網(wǎng)格失配問題，提出了基于擾動項(xiàng)的Off-gird TomoSAR成像模型，并采用結(jié)合局部閾值的L1最小范數(shù)優(yōu)化方法求解所提出的模型。最后通過點(diǎn)目標(biāo)仿真，以及在實(shí)測機(jī)載陣列干涉SAR數(shù)據(jù)，驗(yàn)證了所提模型的正確性及所提方法的有效性，通過與傳統(tǒng)L1范數(shù)最小化的處理結(jié)果進(jìn)行對比，證明本文方法能夠更精確地估計(jì)目標(biāo)位置，并得到更準(zhǔn)確的目標(biāo)散射系數(shù)估計(jì)值。