基于稀疏和低秩結構的層析SAR成像方法

趙 曜 許俊聰 全相印 崔 莉 張 柘

①(廣東工業(yè)大學 廣州 510006)

②(中國運載火箭技術研究院 北京 100076)

③(北京市遙感信息研究所 北京 100192)

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

⑤(蘇州市空天大數(shù)據(jù)智能應用技術重點實驗室 蘇州 215000)

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

1 引言

層析SAR成像(Tomographic SAR,TomoSAR)是一種先進的SAR干涉測量技術，該技術能夠監(jiān)測并獲取目標區(qū)域的高分辨率雷達影像，其中城市和森林區(qū)域是層析SAR成像技術的重要研究對象[1,2]。不同于傳統(tǒng)的合成孔徑雷達(Synthetic Aperture Radar,SAR)只能進行二維成像，層析SAR成像將合成孔徑雷達的合成孔徑原理擴展用于三維成像的高程方向，能夠重建散射體的三維信息并反演高程向剖面，其可以有效解決在二維SAR成像中處于同一散射單元內的目標散射點與雷達間斜距相等時存在的疊掩效應[3–5]。相比于干涉合成孔徑雷達(Interferometric SAR,InSAR)，層析SAR成像技術不僅可以獲得目標散射體的高程信息，同時還可以獲得散射體在高程向上的分布，能夠完全恢復真實三維場景[6]?；谝陨系膬?yōu)點，TomoSAR成為現(xiàn)在最受關注且最具潛力的三維成像方法之一。

目前國際上主流的層析SAR成像方法主要包括多航過成像與陣列SAR兩種，其體制不同但原理和信號處理方法是類似的。但多航過成像隨著航過數(shù)的增加，成像的時間會大大增加，時效性較差；而陣列SAR系統(tǒng)的成本、重量與復雜度也會隨著通道數(shù)的增加而急劇增加。因此，目前國內外的層析SAR研發(fā)方向都是著力于減少通道/航過數(shù)。德國宇航局(Deutsches zentrum für Luf-tund Raumfahrt,DLR)已經(jīng)做到了3～5軌的層析SAR成像[6]，還在繼續(xù)向下減少通道/航過數(shù)。我國提出的“微波視覺”[4]等新理論也在努力減少通道數(shù)。

層析SAR成像的信號處理方法研究始于20世紀80年代，在數(shù)十年的研究和發(fā)展中，其成像效果也在不斷改善。譜估計方法是層析SAR成像的經(jīng)典方法，包括CAPON算法[7,8]、MUSIC算法[9]等，但譜估計方法CAPON和MUSIC需要估計空間的協(xié)方差矩陣，而這個過程會降低方位-距離單元高程向的分辨率[10,11]。稀疏信號處理方法，尤其是壓縮感知(Compressed Sensing,CS)[12–14]作為近年來的研究熱點，也被應用到層析SAR成像中，該方法在相同基線數(shù)量下可以有效提高高程向的分辨率。文獻[15]針對目標區(qū)域的小波域中的稀疏表征，同時結合相應多基線測量的二階統(tǒng)計量，實現(xiàn)目標三維成像?；趬嚎s感知的層析SAR成像方法要求場景是點目標并且具有一定的稀疏性，壓縮感知的重構算法主要包括貪婪算法和凸優(yōu)化算法兩類，貪婪算法一般計算速度較快但是在低信噪比和低采樣數(shù)條件下重建精度較低，而凸優(yōu)化算法的重建精度較高卻計算效率較低。

城市和森林區(qū)域是層析SAR的重要應用對象。森林區(qū)域具有復雜的電磁散射特性，森林在高程向不具有典型的稀疏特性。同時城市區(qū)域中的建筑分布在高程向上屋頂部分通常具有稀疏特性，但地面上由于存在大量人工設施，例如花壇等，導致城市區(qū)域的地面不具有典型的稀疏特性。這為將稀疏信號處理應用于城市和森林區(qū)域提出了一定的挑戰(zhàn)。

KL變換(Karhunen Loeve Transform,KLT)是一種數(shù)據(jù)衍生正交變換，它充分利用數(shù)據(jù)的相關性得到目標的自適應表征形式，在MRI等場景中已經(jīng)得到了充分應用[16]。應用KL變換處理一個信號，其離散展開式為

其中，X(t) 為 輸入信號，為展開式的系數(shù)，en(t)為信號的特征函數(shù)，特征函數(shù)需要根據(jù)輸入信號來確定函數(shù)的類型。特征函數(shù)和信號自相關函數(shù)E{XkXl}有以下關系：

此外在KL變換中，其展開式的系數(shù)Zn的范數(shù)為特征值λn，表示相對應的特征函數(shù)en(t)的能量值，且秩為非零特征值λn的個數(shù)。同時，KL變換在均方誤差準則下與其他正交變換相比效果更佳，是失真最小的一種變換，因此本文引入KL變換對目標的高程向進行處理，以獲得更佳的成像結果?；贙L變換的重構問題可以轉換為譜正則化的矩陣恢復問題，通過低秩矩陣的模型進行求解。由于目標區(qū)域在相鄰方位-距離單元的高程向分布具有較強相關性，本文針對低航過數(shù)和低通道數(shù)的情況提出了一種基于稀疏和低秩結構的層析SAR成像方法，利用目標區(qū)域的低秩結構特性，引入KL變換對目標的方位-距離單元的高程向進行處理，表征其低秩結構特性，以此獲得低航過數(shù)和低通道數(shù)情況下高分辨層析SAR成像結果。本文提出的基于KL變換的層析SAR成像方法不僅可以處理稀疏點目標場景，也可以處理分布式目標場景，其優(yōu)化模型比傳統(tǒng)壓縮感知方法更為復雜，傳統(tǒng)的壓縮感知重構算法不再適用，本文采用交替方向乘子法(Alternating Direction Method of Multipliers,ADMM)優(yōu)化算法，將復雜的原優(yōu)化問題分解為若干相對簡單的子問題，通過優(yōu)化變量交替投影的方式進行算法求解，其計算復雜度與傳統(tǒng)壓縮感知在同一量級。

本文結構如下：第2節(jié)給出層析SAR成像的信號模型；第3節(jié)介紹基于稀疏和低秩結構的層析SAR成像方法；第4節(jié)利用仿真和實測數(shù)據(jù)的實驗結果說明本文算法能夠更好地重建三維目標場景；第5節(jié)對全文進行總結。

2 層析SAR成像信號模型

層析SAR成像沿每一個軌道飛行獲得一副聚焦的二維SAR復數(shù)圖像，設gn是第n次獲取的二維SAR圖像，它可以理解為高程向散射單元的傅里葉變換的非均勻采樣結果[17]：

其中，x代表方位-距離二維平面位置，s代表高程向高度，γ(x,s)代 表高程向后向散射系數(shù)，λ代表載波波長，l代表觀測點到場景中心的斜距，bn代表第n條基線的高程向高度，ε代表噪聲。通過對式(3)進行離散化處理，可以將式(3)記為g=R(γ)+ε，其中R表示層析SAR成像的觀測算子。

通過將場景的后向散射系數(shù)γ(x,s)進行重新排列，可以分析目標在相鄰方位-距離單元的高程向分布的相關性，構造矩陣

其中，矩陣的列表示給定高度的方位-距離平面成像結果，也就是方位-距離平面二維成像結果按照方位向展開為一維向量，而矩陣的行表示給定方位-距離單元的高程向成像結果，矩陣的空間意義為三維網(wǎng)格每個格點的復散射信息。

以城市中的建筑為例，一個建筑物的外墻高度一致，所以一個高度的建筑物可以看作一個秩為1的成像結果矩陣，目標區(qū)域內建筑物的高度差異越復雜，則成像結果矩陣的秩越大。而對于森林區(qū)域，一種類型的樹木其高程向分布基本一致，不同類型的樹木其高程向分布則有一定差異，由此樹木種類越豐富，則成像結果中的秩越大。在常見的城市、森林場景中，典型目標如建筑物、樹木等，其種類較為一致，所以成像結果矩陣可預期具有先驗低秩性。因為矩陣的行線性相關，所以矩陣Γ的秩r ≤min(N,L) 。由矩陣的相關性，對矩陣Γ進行KL分解，可表示為以下3個矩陣的乘積

由此γ(x,s)可表示為r個高程基函數(shù)的加權線性組合

其中，ρi(x)是 矩陣UΣ的行，高程基函數(shù)νi(s)是矩陣V的列。

3 基于ADMM的層析SAR成像方法

本文在稀疏信號處理的基礎上利用相鄰方位-距離單元的高程向的相關性，引入KL變換來表征其低秩結構特性，提高目標區(qū)域層析SAR成像能力。根據(jù)低秩矩陣的重構理論，矩陣Γ ∈RN×L的秩小于等于 min(N,L)，可通過其量測值g=R(Γ)利用正則化方法進行高精度重建。由此形成如下稀疏與低秩結構相結合的層析SAR成像模型[18,19]：

其中，R代表層析SAR成像的觀測算子，ψ(Γ)=,Φ為小波變換矩陣[15]，ψ(Γ)表示矩陣中行向量小波系數(shù)的l1范數(shù)，用以表征矩陣Γ的稀疏特性，表示矩陣Γ的Schatten-P范數(shù)，Γ=UΣV *表 示矩陣Γ的KL分解，Σ=diag([σ0,σ1,...,σr-1]) ,φ(Γ)用以表征模矩陣低秩結構特性，λ1和λ2表示模型的正則化參數(shù)。本文方法和CS方法都需要利用目標場景的先驗信息，對于CS而言，對于城市的建筑一般假設在高程向具有2～3個散射點，而森林樹冠、城市地面等由于不具有典型稀疏性，其散射點個數(shù)往往較多，需要設置合適的矩陣秩。在理論上，根據(jù)稀疏信號處理理論與陣列信號處理理論，三維重建所需通道數(shù)/軌數(shù)最少也應大于目標的自由度，自由度在本文中體現(xiàn)為成像結果矩陣的秩[6]。

由于式(7)同時使用到了稀疏和低秩結構，本文采用ADMM方法對其進行求解，式(7)可以轉化為

其中，K,S都是輔助變量，β1,β2為正則化參數(shù)。當β1,β2→∞時，式(8)等價于式(7)的解。對于式(8)，本文使用如下的3步交替最小化方法進行求解

其中，Q(Γ)=Γ Φ。

第2個子問題與標準核范數(shù)最小化問題類似，可將求解核范數(shù)最小化的迭代奇異閾值方法應用到式(10)中

其中，ui,vi,σi分 別是ΓL+1的 奇異向量和奇異值，為奇異值閾值函數(shù)，閾值函數(shù)的定義如下所示

第3個子問題需要用到多維閾值的方法，具體計算如下：

通過不斷迭代求解上述3個子問題，最終得到式(8)的解。算法流程如圖1所示。

圖1 算法流程圖Fig.1 Algorithm flowchart

4 實驗及結果分析

4.1 實驗數(shù)據(jù)

為了說明基于稀疏和低秩結構的層析SAR三維成像算法的有效性，本文使用了仿真數(shù)據(jù)和機載實測數(shù)據(jù)進行實驗，機載實驗包含兩組實測數(shù)據(jù)，一組為由德國宇航局(DLR)的機載F-SAR系統(tǒng)提供的HV極化通道的L波段數(shù)據(jù)[20]，另一組為中科院空天院機載陣列干涉SAR系統(tǒng)獲取的HH極化通道的Ku波段數(shù)據(jù)[21]。其中L波段數(shù)據(jù)總航過數(shù)為9,Ku波段數(shù)據(jù)總通道數(shù)為12。

近年，寧夏農村飲水工作取得了很大進步，農村水利基礎條件有了很大改善，自來水受益率達到了62%，農村的飲水基本得到了保障。隨著農村飲水和鄉(xiāng)鎮(zhèn)供水的持續(xù)發(fā)展，農村、鄉(xiāng)鎮(zhèn)居民生活用水量不斷增加，但小城鎮(zhèn)和廣大農村的居住區(qū)缺乏排水設施，更談不上污水的處理和利用，影響了可持續(xù)發(fā)展，部分地區(qū)群眾期望農村供排水設施同步發(fā)展。此外，一些工程管理仍在沿用計劃經(jīng)濟體制下的管理模式，管理意識淡薄，管理方式和管理手段落后，服務跟不上，與廣大群眾的期望仍有差距。

首先，仿真場景模擬了高程向0 m到90 m范圍的森林區(qū)域，在該場景中，本文對不同航過數(shù)、稀疏度進行了模擬。圖2(a)模擬了兩個散射體，一個表示地面，另一個表示樹冠散射中心，這兩個散射體分別位于高程向9 m處和60 m處，對應的后向散射系數(shù)相等。圖2(b)模擬了3個散射體，其中9 m處的表示地面，50 m和79 m處表示樹冠散射中心，3個散射體的后向散射系數(shù)一致。仿真場景的相位服從 (0,2π)的瑞利分布。仿真和L波段的F-SAR系統(tǒng)的實測數(shù)據(jù)所采用的雷達系統(tǒng)參數(shù)一致。雷達中心頻率為1.325 GHz，方位向分辨率為0.4 m，距離向分辨率為1.5 m，實驗分別采用9次航過數(shù)據(jù)和6次航過數(shù)據(jù)進行層析SAR成像。文獻[20]給出了觀測數(shù)據(jù)中的殘余相位補償文件和由基線間隔和中心斜距所計算的波數(shù)結果。

圖2 高程向歸一化能量分布Fig.2 Normalized energy distribution in the elevation direction

其次，L波段F-SAR實測數(shù)據(jù)的采集區(qū)域為德國Trockenborn的部分區(qū)域，圖3所示為該區(qū)域的光學圖像，主要由森林和平地構成，非常適合驗證本文方法在森林區(qū)域層析SAR成像的成像效果。Ku波段的實測數(shù)據(jù)為峨眉區(qū)域的城市數(shù)據(jù)，圖4(a)為該區(qū)域的光學圖像，圖4(b)為該區(qū)域的SAR圖像，該雷達系統(tǒng)的參數(shù)如表1所示。

表1 Ku波段雷達系統(tǒng)參數(shù)Tab.1 Ku-band radar system parameters

圖3 F-SAR觀測場景的光學圖像Fig.3 Optical image of F-SAR measured scene

圖4 峨眉區(qū)域的圖像Fig.4 Image of Emei area

4.2 仿真實驗結果

仿真實驗中，本文使用3種不同類型的方法對仿真數(shù)據(jù)進行重建，這3種方法分別是經(jīng)典譜估計算法CAPON[7,8]、壓縮感知[12–14]方法和基于稀疏和低秩結構的層析SAR成像方法。通過這些方法的仿真實驗結果對比，以體現(xiàn)本文方法在層析SAR成像中的有效性。

首先在全航過(航過數(shù)為9)的情況下，分別采樣以上方法對仿真數(shù)據(jù)進行重建。圖5給出了航過數(shù)為9的情況下這3種方法的重建結果，可以看出CS和本文方法均能夠準確恢復散射體的高程向位置，對歸一化能量的值估計也大致準確，而CAPON算法對歸一化能量值的估計存在誤差，而且分辨率比另外兩個方法要低。

接下來將航過數(shù)減少為6，同樣采樣以上方法對仿真數(shù)據(jù)進行重建。圖6給出了6條航過數(shù)的情況下這3種方法的重建結果，可以看出本文方法在減少航過數(shù)的情況下依然能夠獲得與航過數(shù)為9的情況下幾乎一致的結果；而CS雖然可以恢復散射體的高程向位置，但歸一化能量的估計出現(xiàn)誤差，并且出現(xiàn)了虛假目標；CAPON算法恢復的散射體位置發(fā)生了偏移，無法準確估計目標的高程向位置。

通過對圖5和圖6的4個圖進行分析，可以發(fā)現(xiàn)本文方法在不同航過數(shù)和不同稀疏度條件下的結果均能夠將多個散射體分離并進行重建，且對歸一化能量的估計較好。

圖5 航過數(shù)為9時各算法獲得的高程向后向散射能量分布與仿真結果的對比Fig.5 Comparison of the simulation result of backscattered energy distribution of the elevation obtained by each algorithm with 9 interferograms

圖6 航過數(shù)為6時各算法獲得的高程向后向散射能量分布與仿真結果的對比Fig.6 Comparison of the simulation result of backscattered energy distribution of the elevation obtained by each algorithm with 6 interferograms

最后，本文針對兩個散射體的場景分析噪聲對航過數(shù)為9和航過數(shù)為6的仿真結果的影響，考察不同信噪比各種成像方法性能的對比。本文通過比較先驗后向散射系數(shù)與重建的后向散射系數(shù)的均方誤差(Mean Squared Error,MSE)說明算法的有效性，其具體計算公式如下：

其中，yi表示先驗的后向散射系數(shù)真值，為各方法重建后的后向散射系數(shù)。

實驗結果如圖7和圖8所示。通過圖7可知航過數(shù)為9的情況下，CAPON和其他方法相比誤差較大，且不夠穩(wěn)定，而CS和本文方法隨著信噪比的變化均有逐步減小的趨勢，但本文方法擁有更低的均方誤差，成像結果更加穩(wěn)定。當航過數(shù)減少為6時，從圖8可知本文方法的重建結果在不同信噪比下的均方誤差均優(yōu)于CAPON方法和CS方法。

圖7 航過數(shù)為9時不同信噪比情況下各算法的均方誤差Fig.7 Mean square error of each algorithm at different SNR values with 9 interferograms

圖8 航過數(shù)為6時不同信噪比情況下各算法的均方誤差Fig.8 Mean square error of each algorithm at different SNR values with 6 interferograms

4.3 機載實驗結果

在機載實驗中，本文仍然采用CAPON,CS和基于稀疏和低秩結構的層析SAR成像方法。首先對F-SAR系統(tǒng)森林數(shù)據(jù)進行層析成像處理。通過比較3種不同方法對實測森林區(qū)域數(shù)據(jù)的成像結果，由此驗證本文方法在層析SAR的森林區(qū)域成像的性能優(yōu)勢。

在全航過數(shù)(航過數(shù)為9)的情況下，使用前文提到的3種方法對F-SAR的實測森林數(shù)據(jù)進行層析SAR成像，具體成像結果如圖9所示。通過對比3種方法的成像結果，可以看出3種方法均能夠將實測森林數(shù)據(jù)中的地面和森林區(qū)域分開，但CAPON方法與另外兩種方法相比分辨率較低。

圖9 航過數(shù)為9時各算法獲得的森林區(qū)域重建結果Fig.9 Reconstruction results of the forest area by each algorithm with 9 interferograms

隨機抽取6次航過數(shù)據(jù)，依然使用以上3種方法在降低航過數(shù)的情況下對實測森林數(shù)據(jù)進行層析SAR成像。圖10為航過數(shù)為6的情況下不同算法的成像結果。通過對比圖10中各方法的成像結果，可以看出本文方法在減少航過數(shù)的情況下依然可得到高精度的成像結果，CAPON方法則出現(xiàn)了比較明顯的分辨率下降，而CS對地面的成像效果出現(xiàn)地面不完整的現(xiàn)象。

圖10 航過數(shù)為6時各算法獲得的森林區(qū)域重建結果Fig.10 Reconstruction results of the forest area by each algorithm with 6 interferograms

其次，我們將上面所說的3種方法應用于中科院空天院陣列干涉SAR系統(tǒng)獲取的峨眉區(qū)域數(shù)據(jù)，對其進行層析成像處理。通過比較3種方法在峨眉區(qū)域的層析成像結果，說明本文方法在城市區(qū)域成像的有效性。

在全通道數(shù)(通道數(shù)為12)的情況下，各方法的成像結果如圖11所示，可以看出CAPON方法對峨眉區(qū)域數(shù)據(jù)的建筑物進行層析成像結果較差，難以形成建筑外形，而CS和本文方法均能夠區(qū)分出地面與建筑，且層析成像結果比較清晰。

圖11 通道數(shù)為12時各算法的峨眉區(qū)域重建結果Fig.11 Reconstruction results of Emei area by each algorithm with 12 channels

圖12為通道數(shù)為8的情況下不同算法的成像結果。通過與圖12對比可以看出CAPON方法對峨眉區(qū)域進行層析SAR成像的結果仍然較差，無法分辨建筑外形，CS雖然能夠區(qū)分地面與建筑，但其成像結果中建筑頂部出現(xiàn)較多偽影，而本文方法在減少通道數(shù)的情況下依然能夠區(qū)分峨眉區(qū)域的地面和建筑，并且基本不會出現(xiàn)偽影現(xiàn)象。

圖12 通道數(shù)為8時各算法的峨眉區(qū)域重建結果Fig.12 Reconstruction results of Emei area by each algorithm with 8 channels

本文對3種方法的成像結果進行距離向上的切片，并把3種方法的切片結果置于一張圖上進行比較，如圖13所示。

對圖13進行分析，可以發(fā)現(xiàn)CAPON的成像結果中無法識別出建筑物；雖然CS方法能夠對建筑和地面進行成像，但其在高度150 m附近出現(xiàn)了比較明顯的偽影現(xiàn)象，如圖中黑圈所示；而本文方法不僅對建筑和地面進行了高精度成像，并且基本沒有出現(xiàn)偽影現(xiàn)象。

圖13 各方法成像結果的距離向切片對比Fig.13 Comparison of range slices of imaging results of each method

表2是各方法對峨眉數(shù)據(jù)進行成像的時間對比，實驗運行平臺：CPU為i9-10900X，內存為64GB，使用MATLAB軟件計算。不難看出CAPON耗時非常少，本文方法耗時約為CS的1.67倍，基本在同一計算量級。

表2 各方法耗時對比Tab.2 Time comparison of each method

5 結束語

本文提出了一種基于稀疏和低秩結構的層析SAR三維成像方法。本文從目標區(qū)域的特征入手，引入KL變換對相鄰方位-距離單元高程向的低秩結構特性進行表征，構建稀疏和低秩結構相結合的層析SAR成像模型，采用ADMM方法推導層析SAR重構算法，獲得高質量的三維重建結果。仿真與實測數(shù)據(jù)的實驗結果表明，相比于傳統(tǒng)的譜估計方法和稀疏信號處理方法，本文方法可以有效實現(xiàn)地面、樹冠散射中心和建筑物成像，且在較少航過數(shù)和通道數(shù)時，保證目標區(qū)域三維重構精度。