基于稀疏重構的全極化SAR聯合多維重建

孫 豆 路東偉 邢世其*楊 瀟 李永禎 王雪松

①(國防科技大學電子信息系統復雜電磁環境效應國家重點實驗室 長沙 410073)

②(國防科技大學電子科學系 長沙 410073)

1 引言

合成孔徑雷達(Synthetic Aperture Radar,SAR)具有全天時全天候對目標場景進行觀測的優勢，日益成為目標識別的關鍵傳感器[1,2]。極化是電磁波的一個重要本質屬性，深刻反映目標形狀、結構和姿態等物理特征，對于提升目標識別能力具有重要作用[3,4]。近年來，SAR越來越朝著全極化和三維的方向發展[5–8]。三維信息使得雷達對目標的描述越來越精細，極化使得雷達對目標散射特性的描述更加完整，將全極化信息和三維散射信息相結合可以得到人造目標更為全面的重建結果，有利于提高目標分類識別精度。

梳理現有的文獻資料，將極化三維重建相關的研究工作進行分類。按照成像方法的不同，相關研究可以分為兩類。一類是基于傅里葉變換的成像方法[9,10]，這類方法是傳統二維SAR成像方法的三維擴展，具有快捷高效的優勢，但是這類方法只能單通道獨立成像，無法進行多通道聯合同時成像。另一類是基于稀疏重構的成像方法[11–13]，由于稀疏成像本質是一種基于模型匹配的方法，因此這類方法的成像結果分辨率高，且對旁瓣和噪聲的抑制效果好，此外，稀疏重構類的方法可以進行多通道聯合處理。按照成像方式的不同，相關研究可以分為兩類。一類是先進行二維成像再獲取高度維信息[9,11,14]，這類方法可以使用現有的二維成像方法，只需進一步對高度信息進行估計，其優勢在于計算復雜度低，效率高，但是成像結果受二維圖像配準質量的影響。另一類是直接三維成像[10,13,15]，即同時得到目標的三維信息，這類方法的優勢在于可以保證數據之間的關聯性，不需要進行額外的圖像配準，但是直接三維成像會面臨巨大的計算壓力。按照極化信息的提取方式的不同，相關研究也可以分為兩類。一類是各極化通道獨立處理[16,17]，即先每個通道分別成像，再提取極化散射矩陣。這種方式的優勢在于快捷，適用性廣，但是由于忽略了各個極化通道數據之間的關聯性，不能保證提取的極化散射矩陣準確。另一類是所有極化通道聯合處理[11,18]，即所有通道聯合成像，同時得到極化散射矩陣，這種聯合的方式保證了極化信息提取的準確度，可以有效克服獨立處理中散射中心失配的問題，但是聯合處理的計算復雜度較大。

現有文獻所采用的方法都是對上述3種分類的不同組合，但是還沒有方法將稀疏成像、直接三維成像和極化聯合處理這三者相結合。實際上，相比于基于傅里葉變換的成像，稀疏成像在分辨率、抗噪性等方面有著更多的優勢[19,20]。相比于直接三維成像，先二維后一維的成像方法很可能破壞散射的一致性，特別是每個二維成像結果中散射點的數目和位置估計可能不同，這種不一致性會使得高度信息的提取不準確。此外，由于每個極化通道觀測相同的區域，各通道數據高度相關，且之間存在著豐富的互補信息，三維信息和極化信息同時聯合獲取不破壞數據之間的關聯性，才能充分利用這些互補信息來提高估計精度和分辨率。因此，考慮到稀疏成像、直接三維成像和極化聯合處理這三種方式各自的優勢，將它們結合進行極化多維重建可以得到質量、精度更高的重建結果。

為了充分利用數據之間的關聯性帶來的額外信息，實現極化散射矩陣和目標三維信息的同時獲取，本文提出基于稀疏重構的全極化聯合多維重建方法。首先建立各極化通道的直接三維稀疏重建模型，然后對所有通道設置聯合稀疏約束，將全極化三維重建建模為多通道聯合稀疏重構問題。考慮到上述稀疏重建問題規模大、計算壓力大，本文通過數據插值對模型進行簡化，并結合三維快速傅里葉變換、共軛梯度法和牛頓迭代法給出了一種高效的模型求解方法。最后，對三維極化散射矩陣進行極化分解，得到包含目標三維位置信息以及散射類型的多維重建結果。本文提出的全極化聯合多維重建方法保證了所有極化通道中散射中心的數目和三維位置一致，有效地克服了散射中心失配的問題，確保了極化特征提取的準確性。此外，該方法大大降低了多維聯合處理面臨的計算壓力，能夠在不改變結果精度的情況下有效地獲取目標多維重建結果。基于仿真數據和電磁計算數據的實驗驗證了本文方法的有效性。

2 單通道信號模型

假設雷達距離場景足夠遠，位于相對于場景中心方位角?和俯仰角θ的位置，并發射中心頻率為fc，帶寬為BW的觀測信號，使用平面波模型，接收信號可以表示為

其中，G(k,0,0)為反射率函數的三維傅里葉變換在(k,0,0)處的值。

以角度(?,θ) 對和(k,0,0)旋轉，得到旋轉后的(x,y,z)和 (kx,ky,kz)。根據傅里葉變換的旋轉特性，式(3)變為

由式(4)可知，反射率函數是波數域采樣G(kx,ky,kz)的三維傅里葉逆變換。

3 聯合稀疏模型建立

對于每個極化通道，需要被重建的成像場景可以通過式(6)獲得

其中，gl(x,y,z;?,θ)為第l(l1,2,···,L)個極化通道的場景反射率函數，Gl(kx,ky,kz)是第l個極化通道的波數域采樣。

在三維重建空間中，定義N個位置為候選散射中心

通常，這些位置是從均勻的矩形網格中選擇的。基于這些位置，定義M ×N的字典矩陣

其中，m表示波數域采樣的索引，n表示C中N個位置的索引。

根據字典矩陣A，將式(6)寫成矩陣形式，可得

其中，N ×1的向量βl是要重建的第l個極化通道的三維場景，M ×1的向量bl表示第l個極化通道的波數域采樣。

每個極化通道單獨三維重建就是要解決下面的直接三維稀疏成像模型

對式(10)中的直接三維稀疏成像模型進行求解則可以同時得到目標的三維信息。相比于先二維稀疏成像再第三維稀疏重建的方法，直接三維稀疏成像避免了二維圖像稀疏支撐集的不一致。

對于每個極化通道，由于其觀測場景是相同的，且采樣方式不變，因此，各極化通道共享相同的字典矩陣A。基于此，建立聯合成像模型

其中β[β1···βl···βL]是聯合稀疏解矩陣，‖β‖0‖|β1|+···+|βl|+···+|βL|‖0表示稀疏解矩陣β中不為零的行數。

在聯合成像模型中，每個極化通道中βl的支撐集是相同的，只是不同極化通道的βl的幅度不同。聯合稀疏約束項‖β‖0‖|β1|+···+|βl|+···+|βL|‖0保證了不同極化通道的重建結果共享相同數量和位置的散射中心。聯合稀疏為散射系數的估計提供了額外的約束條件，使得在重建過程中可以保留更多有價值的信息。

4 全極化聯合多維重建方法

式(11)中的數學模型不是一個凸優化問題，而是一個NP難問題。為了用數值方法求解該模型，首先對目標函數進行放松，定義混合規范

其中βl(i)表示βl中的第i個元素。混合范數?2,p沿不同極化通道的方向計算?2范數，即計算β的張量，并沿散射位置的方向計算?p范數。

結合式(12)中的混合范數，式(11)中的模型被松弛為式(13)的優化問題

其中μ是控制稀疏性的正則化參數。聯合多維重建就是對式(13)中的模型進行求解。

為了求解式(13)，本文在Cetin等人[21]處理思想的基礎上，提出了一種快速迭代方法。定義代價函數

當代價函數J取最小值時，可以得到稀疏重建的結果。計算J對β的偏導數，得到

其中D(β)是N ×N的對角矩陣，b[b1···bl···bL]。

當?βJ0時，則得到了散射系數的估計。考慮到式(15)中的2AHA+μpD(β)可近似看作Hessian矩陣，因此可以采用近似高斯迭代法求解?βJ0，得到迭代式：

其中?n+1是在迭代中的步長，為了使算法快速收斂，在迭代過程中采用變步長?n+1(?n)0.9。當小于預設閾值時，算法退出，此時就得到了所有極化通道的三維重建結果。為了進一步提取目標極化特性，還可對式(16)得到的極化散射矩陣進行Cameron分解[22]，最終就得到了包含極化信息和目標三維信息的多維重建結果。

實際上，式(16)的迭代計算還存在一些問題。由于本文方法進行三維直接重建，同時獲取目標的全極化三維信息，不可避免地面臨計算量巨大的問題，即式(16)中的字典矩陣維度很大，無法直接存儲計算。具體地，字典矩陣A是M ×N的，N與場景的大小和分辨率有關,M與波數域的采樣點數有關。假設要重建的場景大小為5 m×5 m×5 m，每個維度的分辨率是0.05 m，那么三維重建結果的維度是100×100×100，則N1×106。假設頻率的采樣點數為100，方位向的采樣點數為100，俯仰向的采樣點數為100，則M1×106。在這種情況下，由于M ×N1×1012，字典矩陣A需要的存儲空間大小約為3GB。除了要存儲巨大的字典矩陣之外，每次迭代都不可避免地要進行字典矩陣的乘法運算，如計算Aβl，需要M ×N1×1012次乘法和加法計算，因此直接對式(16)進行迭代的計算復雜度很高。為了解決上述計算規模大的問題，本文借鑒文獻[12]中的思路，結合三維非均勻快速傅里葉變換(3-D Non-Uniform Fast Fourier Transform,3-D NUFFT)和共軛梯度法求解式(16)中的部分。

由式(5)可知，kx,ky,kz與f,θ,?之間的關系不是線性的，因此kx,ky,kz處于非均勻網格上。對于字典矩陣雖然x,y,z處于均勻的網格上，但由于kx,ky,kz是非均勻的，因此字典矩陣不可以被看作三維傅里葉變換矩陣。為了可以使用三維快速傅里葉變換(3-D Fast Fourier Transform,3-D FFT)代替字典矩陣，就需要在稀疏重構前對每個極化通道的波數域采樣進行預處理，使得波數域采樣處于均勻的網格上。

這里使用3-D NUFFT對非均勻網格kx,ky,kz上的每個極化通道的波數域采樣Gl(kx,ky,kz)進行插值，插值后的波數域采樣處于均勻網格上。對應于插值后的波數域采樣，字典矩陣則更新為式(16)更新為

表1總結了本文所提出的全極化聯合多維重建方法的具體步驟。

該聯合多維重建方法實現了極化散射矩陣和目標三維信息的同時獲取，充分利用了數據之間的關聯性帶來的額外信息。實際上，本文方法的關鍵在于反映了聯合稀疏性的矩陣包含了所有極化通道和全部三維的散射系數。聯合稀疏性為散射系數的估計提供了額外的約束條件，保證了不同極化通道、不同維度的稀疏支撐集是一致的。不管散射中心是不是在所有極化通道都有散射信息，聯合稀疏對4個極化通道進行聯合，同時估計散射中心在各個極化通道的散射信息，有效地克服了獨立處理和三維分步成像中散射中心失配的問題，進而保證了獲取的極化散射矩陣的準確性。此外，同時獲取極化和三維信息不可避免地要面對計算規模大的問題。本文方法結合3-D FFT和共軛梯度法對模型進行簡化，避免了巨大字典矩陣的存儲，大大降低了計算壓力，使聯合多維重建可以在不改變解的精度的前提下高效地進行。

表1 全極化聯合多維重建方法的步驟Tab.1 Steps of full polarization joint multi-dimensional reconstruction method

5 實驗結果與分析

本節基于仿真數據和電磁計算數據開展實驗，驗證本文提出方法的有效性。實驗在一臺配備Intel Core I5-6500 CPU和12 GB RAM的計算機上的MATLAB R2016b中進行。

5.1 重建性能分析

選取三面角，偶極子，30°二面角和45°二面角作為仿真目標，進行全極化聯合三維重建，并從目標類型依賴性、極化散射矩陣估計精度、噪聲敏感度這3個方面分析本文方法的性能。表2給出了4個仿真目標的信息，其中三面角，偶極子和45°二面角的稀疏支撐集不一致，30°二面角的稀疏支撐集一致。回波數據由MATLAB生成，使用理想點散射模型，表3記錄了仿真實驗時的各項參數。

表2 仿真目標信息Tab.2 Information of simulated targets

表3 目標的仿真參數Tab.3 Simulation parameters of simulated targets

圖1給出了仿真目標的4個極化通道的成像結果，其中左列是HH,HV,VH和VV極化通道下的二維投影結果，右列是HH,HV,VH和VV極化通道下的三維結果。以HH極化通道成像結果的峰值為0 dB，圖1中所有子圖顯示成像結果的幅度閾值為–20 dB。觀察HH通道的成像結果有3個點目標，對應三面角、偶極子、30°二面角；HV通道和VH通道的成像結果有兩個點目標，對應30°二面角和45°二面角；VV通道的成像結果有兩個點目標，對應三面角和30°二面角。各個極化通道的成像結果和目標的極化散射矩陣可以對應上，目標極化散射矩陣分量是0的通道的成像結果沒有該目標。

圖1 各個極化通道的仿真目標成像結果Fig.1 Each polarization channel’s imaging results of simulated targets

圖2給出了仿真目標的多維重建結果，其中圖2(a)為全極化三維成像結果，圖2(b)為Cameron分解結果。對比圖2(a)的全極化三維成像結果和圖1中的各個極化通道的三維成像結果，可以看出圖2(a)融合了所有極化通道的結果，包含了全部的四個目標，對應三面角、偶極子、30°二面角和45°二面角。圖2(b)中不同的圖標表示由Cameron分解得到的不同的散射類別。對比圖2(b)中的結果和表2中的目標信息可知，各散射類型判別結果與實際情況相符合，驗證了本文方法對極化散射矩陣提取的有效性。

為了進一步分析本文提出的聯合多維重建方法對極化散射矩陣的估計精度，表4給出了仿真目標的極化散射矩陣估計結果，并同時列出了各個極化通道分別獨立多維重建的估計結果進行對比。變型前的結果是多維重建得到的極化散射矩陣原始值。為了和目標真實的極化散射矩陣進行對比分析，我們對估計得到的每個散射矩陣進行歸一化處理，將幅度和相位提到矩陣外，得到變形后的極化散射矩陣估計值。對于變形后的極化散射矩陣估計值，藍色部分表示各個目標的幅度估計結果，紅色部分表示各個目標的極化散射矩陣估計結果，指數部分表示與當前場景有關的相對相位，不影響結果的分析。從表4中聯合多維重建方法得到的變型后的結果可以看出，各個目標的散射幅度相近，這與表2中目標真實的幅度情況一致。由于NUFFT插值和3-D FFT會帶來增益的變化，因此幅度的估計結果不是1。與表2中目標真實的極化散射矩陣對比，可以看出對于三面角、偶極子和45°二面角，聯合多維重建方法得到的極化散射矩陣的估計值和真實值完全一致，對于30°二面角，聯合多維重建方法得到的極化散射矩陣估計結果基本上和真實值一致，交叉極化通道的估計結果略微大于真實值。觀察表4中獨立多維重建法得到的變型后的結果，可以看出極化散射矩陣的估計結果并不準確。具體地，30°二面角的散射幅度和其他目標的散射幅度差異很大，這與表2中目標真實的幅度情況不一致。此外，對于30°二面角，獨立多維重建方法得到的極化散射矩陣估計結果與真實值差異很大，交叉極化通道的估計結果遠大于真實值。因此，各個極化通道分別獨立多維重建不能保證極化散射矩陣提取的準確性，而聯合多維重建方法則可以得到較準確的極化散射矩陣的估計結果。

為分析本文方法對噪聲的敏感度，仿真中給波數域采樣G(kx,ky,kz)增加噪聲，在圖像域信噪比(Signal to Noise Ratio,SNR)為13 dB,18 dB,23 dB時分別進行三維聯合重建，并記錄極化散射矩陣估計結果如表5所示。從表5中變型后的結果可以看出，當SNR為18 dB和23 dB時，藍色部分表示的各個目標的散射幅度相近，這與表2中目標真實的幅度情況一致。此外，紅色部分表示的各個目標的極化散射矩陣估計結果也與表2中目標真實的極化散射矩陣基本一致。當SNR=13 dB時，可以看出各個目標的散射幅度估計結果差異略大，特別是30°二面角和45°二面角的散射幅度估計結果相差0.09。此外，與表2中目標真實的極化散射矩陣對比，可以看出對于三面角、偶極子和45°二面角，其極化散射矩陣的估計值和真實值基本一致，但對于30°二面角，受噪聲影響，其極化散射矩陣的估計值和真實值略有差異。

圖2 仿真目標全極化聯合多維重建結果Fig.2 Full polarization joint multi-dimensional reconstruction results of simulated targets

表4 仿真目標的極化散射矩陣估計結果Tab.4 Polarization scattering matrix estimation results of simulated targets

表5 不同SNR下仿真目標的極化散射矩陣估計結果Tab.5 Polarization scattering matrix estimation results of simulated targets under different SNR

根據上述仿真分析可得，本文方法的性能不受目標類型影響，不論目標的稀疏支撐集是否一致，本文方法都可以較準確地估計出極化散射矩陣。相比較于各個極化通道分別獨立多維重建不能保證極化散射矩陣提取的準確性，本文的聯合多維重建方法對極化散射矩陣的提取準確且估計精度高。此外，本文方法具有一定的抗噪性，在SNR=13 dB時，極化散射矩陣的估計結果開始變差。

5.2 典型目標重建結果

使用Slicy模型的電磁計算數據進行仿真實驗。該模型由三面角、二面角、圓柱等組成，圖3給出了其三維CAD模型及尺寸信息。回波數據由CST微波工作室生成，采用彈跳射線法。表6記錄了Slicy仿真實驗中的各項參數。

圖3 Slicy的CAD模型Fig.3 CAD model of Slicy

表6 Slicy的仿真參數Tab.6 Simulation parameters of Slicy

圖4給出了Slicy模型的4個極化通道的成像結果，其中左列是HH,HV,VH和VV極化通道下的二維投影結果，右列是HH,HV,VH和VV極化通道下的三維結果。以HH極化通道成像結果的峰值為0 dB，圖4中所有子圖顯示成像結果的幅度閾值為–30 dB，可以看出成像結果稀疏，旁瓣基本都被去除掉，驗證了本文基于稀疏重構的重建方法在旁瓣抑制上的有效性。觀察兩個主極化通道，其成像結果相似，都可清晰地看到Slicy模型的外形輪廓，Slicy模型中的兩個帽型結構，一個三面角和兩個二面角都被重建出來，成像結果的位置和尺寸與圖3中的模型吻合。觀察兩個交叉極化通道，目標的散射相對較弱，輪廓不太清晰，有些部分沒有被重建出來，且兩個交叉極化通道的結果略有不同，Slicy模型的一個帽型結構的強度有些差異。

圖4 各個極化通道的Slicy成像結果Fig.4 Each polarization channel’s imaging results of Slicy

圖5給出了Slicy模型的多維重建結果，其中圖5(a)為全極化三維成像結果，圖5(b)圖為Cameron分解結果。對比圖5(a)的全極化三維成像結果和圖4中的各個極化通道的三維成像結果，可以看出圖5(a)和圖4中的主極化通道的結果基本一致，不存在散射中心失配的問題，這也就驗證了本文的聯合重建方法的有效性，即保證了不同極化通道、不同維度的稀疏支撐集一致。圖5(b)中不同的圖標表示由Cameron分解得到的不同的散射類別。對比圖5(b)中的結果和圖3中的模型可知，各散射類型判別結果基本與實際情況相符合，驗證了本文方法對極化散射矩陣提取的有效性。但由于Cameron分解自身的局限性，散射類型判別存在一定的不確定性，Slicy模型中的二面角和帽型結構沒法分辨。

圖5 Slicy全極化聯合多維重建結果Fig.5 Full polarization joint multi-dimensional reconstruction results of Slicy

使用衛星模型的電磁計算數據進行仿真實驗，圖6給出了其三維CAD模型，衛星中心位于坐標原點，太陽能兩帆板外側的距離為8.21 m，星體最長尺寸3.55 m。回波數據由電磁計算軟件[23]生成，采用積分方程快速計算方法。表7記錄了衛星仿真實驗中的各項參數。

圖7給出了衛星的4個極化通道的成像結果，其中左列是HH,HV,VH和VV極化通道下的二維投影結果，右列是HH,HV,VH和VV極化通道下的三維結果。以HH極化通道成像結果的峰值為0 dB，圖7中所有子圖顯示成像結果的幅度閾值為–30 dB，可以看出成像結果稀疏，衛星的輪廓很清晰，旁瓣基本都被去除掉了，進一步驗證了本文基于稀疏重構的重建方法在旁瓣抑制上的有效性。觀察兩個主極化通道，都可清晰地看到衛星模型的外形輪廓，衛星主體的柱狀結構和衛星兩側的太陽能帆板都被重建出來，與圖6中的模型吻合。不過兩個主極化通道的成像結果略有不同，HH極化通道中衛星兩個帆板的散射強度明顯大于VV極化通道。觀察兩個交叉極化通道，其成像結果相似，目標的散射相對較弱，輪廓不太清晰，衛星兩側的太陽能帆板沒有被重建出來，只有衛星主體柱狀結構被重建出來。

圖6 衛星的CAD模型Fig.6 CAD model of satellite

圖8給出了衛星的多維重建結果，其中圖8(a)為全極化三維成像結果，圖8(b)為Cameron分解結果。對比圖8(a)的全極化三維成像結果和圖7中的各個極化通道的三維成像結果，可以看出圖8(a)和圖7中的HH極化通道的結果基本一致，因此不存在散射中心失配的問題，進一步驗證了本文聯合重建方法的有效性。圖8(b)中不同的圖標表示由Cameron分解得到的不同的散射類別，可以看出衛星的兩個帆板上的散射類型為圓柱體，兩個帆板下沿以及帆板與衛星主體連接處的散射類型為平板/三面角，衛星主體上部的散射類型為圓柱體，衛星主體中部比較復雜，散射類型較多，衛星主體下部的散射類型為平板/三面角和帽型/二面角的組合。由于Cameron分解自身的局限性，二面角和帽型結構，三面角和平板沒法分辨。對比圖8(b)中的結果和圖6中的模型，各散射類型判別結果基本與實際情況相符合，進一步驗證了本文方法對極化散射矩陣提取的有效性。

表7 衛星的仿真參數Tab.7 Simulation parameters of satellite

圖7 各個極化通道的衛星成像結果Fig.7 Each polarization channel’s imaging results of satellite

圖8 衛星全極化聯合多維重建結果Fig.8 Full polarization joint multi-dimensional reconstruction results of satellite

6 結論

本文提出一種基于稀疏重構的全極化聯合多維重建方法。該方法對所有極化通道和所有維度的散射系數進行了聯合稀疏約束，不僅保證了不同極化通道、不同維度的稀疏支撐集一致，并且充分利用了數據之間的關聯性帶來的額外信息。此外，該方法解決了同時獲取極化和三維信息面臨的計算規模大的問題，使聯合多維重建可以在不改變解的精度的前提下高效地進行。實驗結果表明本文方法的性能不受目標類型影響，具有一定的抗噪性，可以實現極化散射矩陣和目標三維信息同時獲取，有效地克服了散射中心失配的問題，得到了準確的目標多維重建結果。此外，該方法保持了稀疏重構的優勢，有效地抑制了旁瓣，得到的重建結果分辨率高。