基于MUSIC算法的Pol-InSAR相位估計(jì)方法

張曉玲，陳 欽，韋順軍

(電子科技大學(xué)電子工程學(xué)院 成都 610054)

極化干涉合成孔徑雷達(dá)(Pol-InSAR)技術(shù)能夠全天時(shí)、全天候獲取大面積數(shù)字高程圖像(DEM)，是目前SAR遙感技術(shù)發(fā)展的重要方向。極化干涉合成孔徑雷達(dá)不僅具有干涉合成孔徑雷達(dá)(InSAR)對散射體的空間分布和高度的敏感性，而且還具有極化合成孔徑雷達(dá)(PolSAR)對散射體的形狀和方向的敏感性。Pol-InSAR通過極化和干涉信息的有效組合，可以同時(shí)提取觀測對象的空間三維結(jié)構(gòu)特征信息和散射信息，為微波定量遙感、高精度數(shù)字高程信息和觀察對象細(xì)微形變信息的提取提供了可能性[1]。

傳統(tǒng)的極化干涉SAR干涉相位估計(jì)方法是直接選取全極化數(shù)據(jù)的一個極化通道(如HH)，利用InSAR圖像配準(zhǔn)的方法進(jìn)行干涉相位估計(jì)，得到配準(zhǔn)偏移量后應(yīng)用于所有極化通道進(jìn)行插值重采樣。由于極化組合方式不一樣，各極化通道的配準(zhǔn)偏移量存在微小的差別，直接以一個極化通道數(shù)據(jù)的偏移量為標(biāo)準(zhǔn)，將影響其他極化通道干涉相位估計(jì)的準(zhǔn)確性。文獻(xiàn)[2]提出了基于加權(quán)聯(lián)合單像素模型的干涉相位估計(jì)方法，是目前InSAR干涉相位估計(jì)精度很高的方法之一，但在估計(jì)干涉相位時(shí)需要首先確定配準(zhǔn)誤差的大小和方向，并且確定最優(yōu)權(quán)值時(shí)需要搜索，因此計(jì)算量較大。

本文提出一種基于MUSIC算法的多極化干涉相位估計(jì)方法，其基本思想是在Pol-InSAR圖像粗配準(zhǔn)的條件下，利用維納濾波器[3]獲得最佳濾波器參量(最優(yōu)權(quán)值)；再利用該參量構(gòu)造多極化聯(lián)合最佳加權(quán)矢量，由該矢量獲得協(xié)方差矩陣和相干系數(shù)矩陣；然后結(jié)合MUSIC算法[4]，特征分解協(xié)方差矩陣得到信號子空間和噪聲子空間，根據(jù)兩者的正交性構(gòu)造空間譜函數(shù)，通過譜峰搜索準(zhǔn)確地估計(jì)各極化通道的干涉相位。

1 多干涉相位估計(jì)模型

1.1 信號模型

在極化干涉合成孔徑雷達(dá)系統(tǒng)中，單次雙天線或重復(fù)軌道單天線測量得到的兩個信號可以分別表示為：

在上面的接收信號模型中，為了準(zhǔn)確地估計(jì)出各極化通道的干涉相位，各極化通道兩幅圖像的同一像素點(diǎn)i必須對應(yīng)地面的同一散射體。而實(shí)際接收到的各極化通道SAR圖像存在像素偏移，其相位差不能反映地面高度起伏的情況[5]。因此，首先進(jìn)行各極化通道SAR圖像配準(zhǔn)過程。

1.2 多極化SAR圖像精配準(zhǔn)

干涉復(fù)圖像對的配準(zhǔn)基本上分為兩個步驟：1) 進(jìn)行粗配準(zhǔn)[6-7]得到像素級的精度，對應(yīng)的配準(zhǔn)誤差為一個分辨單元；2) 通過精配準(zhǔn)得到亞像素級的精度，一般要求達(dá)到0.1～0.01分辨單元[8-10]。

各極化通道SAR圖像經(jīng)過粗配準(zhǔn)后，針對同一像素點(diǎn)i，副圖像較主圖像配準(zhǔn)誤差為一個分辨單元，并且誤差方向可能存在于副圖像像素i上下左右的任一方向。本文利用維納濾波[3]的相關(guān)理論，對各極化通道SAR圖像數(shù)據(jù)精配準(zhǔn)。在最小均方誤差準(zhǔn)則(MMSE)下，取副圖像中相應(yīng)像素點(diǎn)作為輸入信號，對主圖像像素i作最小均方估計(jì)。最后經(jīng)維納濾波得到的輸出信號就為對主圖像像素i精配準(zhǔn)后的副圖像相對應(yīng)的像素點(diǎn)。

因?yàn)檎`差單元存在于副圖像對應(yīng)像素i上下左右的任一方向，所以取副圖像中對應(yīng)像素i的一個3×3(也可以選取5×5)的矩陣，作為維納濾波器的輸入信號，建立如下的維納濾波器模型[3]：

圖1 維納濾波器模型

式中，S2(i)是由副圖像中對應(yīng)像素i的3×3矩陣所構(gòu)成的回波信號矢量；w(i)為對應(yīng)的維納濾波器參量；e(i)為主副圖像對應(yīng)像素點(diǎn)i的估計(jì)誤差。根據(jù)最小均方誤差準(zhǔn)則，最佳的濾波器參量wopt(i)應(yīng)使得性能函數(shù)——均方誤差為最?。?/p>

利用維納濾波理論對各極化通道SAR圖像精配準(zhǔn)過程中，因?yàn)檎`差單元就存在于副圖像對應(yīng)像素i的一個3×3矩陣?yán)铮恍枰x取該3×3矩陣作為維納濾波器的輸入信號，就可以得到精配準(zhǔn)后的副圖像相應(yīng)像素點(diǎn)。與文獻(xiàn)[2]相比，本文所提方法不需要首先估計(jì)出配準(zhǔn)誤差的大小和方向，降低了這一過程的復(fù)雜性。并且對副圖像每一像素點(diǎn)進(jìn)行同樣的操作，都可以得到主圖像對應(yīng)像素的最小均方估計(jì)，因此具有自適應(yīng)圖像配準(zhǔn)功能。

1.3 干涉相位估計(jì)

根據(jù)前面的信號模型式(1)和式(2)，每個極化通道主副圖像像素點(diǎn)i中含有相位噪聲，利用維納濾波實(shí)現(xiàn)各極化通道SAR圖像精配準(zhǔn)后，相位噪聲仍存在于主副圖像像素點(diǎn)i中。在圖像精配準(zhǔn)達(dá)到亞像素級精度的情況下，由于相位噪聲的存在，主副圖像直接共軛相乘估計(jì)干涉相位是不精確的。而精配準(zhǔn)后利用MUSIC算法特征分解出信號子空間和噪聲子空間，并且噪聲子空間正交于信號子空間，因此具有很強(qiáng)的相位噪聲抑制功能，可提高各極化通道干涉相位估計(jì)的精確度。

2 算法處理步驟

本文提出的基于MUSIC算法的多極化干涉相位估計(jì)方法，首先經(jīng)過維納濾波獲得最佳濾波器參量，實(shí)現(xiàn)多極化SAR圖像的精配準(zhǔn)，然后利用該參量構(gòu)造協(xié)方差矩陣，并結(jié)合MUSIC算法準(zhǔn)確地估計(jì)出信號子空間和噪聲子空間，再根據(jù)兩者的正交性構(gòu)造譜函數(shù)。圖2是利用該算法估計(jì)多極化干涉相位的處理流程。

圖2 本文算法處理流程

具體步驟如下：

1) 利用維納濾波圖像精確配準(zhǔn)。

首先利用傳統(tǒng)的相干系數(shù)法進(jìn)行粗配準(zhǔn)，使配準(zhǔn)誤差為一個分辨單元；接著通過維納濾波器實(shí)現(xiàn)各極化通道SAR圖像的精配準(zhǔn)，達(dá)到亞像素級的精度。維納濾波器中最佳維納濾波器參量w(i)的求解方法如下(以HH通道為例)：

2) 估計(jì)信號子空間和噪聲子空間。

對各極化通道每一像素i執(zhí)行上面的處理過程，就可以準(zhǔn)確地估計(jì)出多極化通道干涉相位。

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

本文選用的極化干涉數(shù)據(jù)是PiSAR全極化(HH,HV, VV)干涉實(shí)測數(shù)據(jù)，實(shí)驗(yàn)中天線1的SAR圖像(主圖像)像素大小為550×450，天線2(副圖像)像素大小為450×450。主圖像和副圖像的總體幅度圖分別如圖3a和圖3b所示，圖中對應(yīng)像素存在明顯的偏移，在進(jìn)行干涉處理前需要對兩幅SAR圖像進(jìn)行配準(zhǔn)。

圖3 HH通道幅度圖

各極化通道SAR圖像經(jīng)過粗配準(zhǔn)后，各極化通道的干涉相位圖如圖4所示。

圖4 粗配準(zhǔn)的極化干涉相位

圖5 多極化干涉相位圖

采用本文基于MUSIC算法的多極化干涉相位估計(jì)方法，得到的HH、HV和VV多極化干涉相位圖如圖5所示。對比圖4和圖5的相位圖可以看出，本文干涉相位估計(jì)方法相對于粗配準(zhǔn)算法可以獲得更加清晰的干涉相位條紋圖，并且極大地抑制了相位噪聲，提高了干涉相位圖的質(zhì)量。

干涉相位圖殘差點(diǎn)數(shù)[13]是衡量干涉相位圖質(zhì)量的主要指標(biāo)之一，殘差點(diǎn)數(shù)過多容易導(dǎo)致后續(xù)相位展開失敗而無法獲取地形高程。經(jīng)過圖像配準(zhǔn)后，殘差點(diǎn)數(shù)越少則證明配準(zhǔn)精度越高。為了比較圖像配準(zhǔn)精度，對相干系數(shù)粗配準(zhǔn)法、文獻(xiàn)[2]中方法及本文方法獲取的干涉相位圖進(jìn)行殘差點(diǎn)統(tǒng)計(jì)，得到的不同極化干涉相位中的殘差點(diǎn)數(shù)如表1所示。

表1 各極化通道干涉相位圖殘差點(diǎn)數(shù)統(tǒng)計(jì)

另外，文獻(xiàn)[2]中的最優(yōu)權(quán)值w的搜索過程還大大增加了計(jì)算量。假設(shè)SAR圖像大小為M×N，從相鄰像素中獲取獨(dú)立同分布的樣本數(shù)為L，則在圖像精配準(zhǔn)過程中，文獻(xiàn)[2]方法所需加法次數(shù)為12(L?1)MN，乘法次數(shù)為12(L+1)MN；本文方法所需加法次數(shù)為3(L?1)MN，乘法次數(shù)為(3(L+1)+1) ×MN。該方法由于避免了w的搜索過程，大大降低了計(jì)算量。同時(shí)，由于利用了MUSIC算法，極大地抑制了相位噪聲，提高了干涉圖的質(zhì)量。

傳統(tǒng)方法以一個極化通道SAR圖像的配準(zhǔn)誤差對其他極化通道圖像作配準(zhǔn)處理是不準(zhǔn)確的，而本文方法估計(jì)了不同極化通道的最佳濾波器參量，即各極化通道的配準(zhǔn)誤差不一樣，從而可準(zhǔn)確得到各極化通道的干涉相位圖。通過該方法準(zhǔn)確地估計(jì)出了多極化干涉相位，提高了極化干涉SAR不同極化組合干涉相位的精度，為極化干涉SAR復(fù)雜場景地形地表參數(shù)精確反演提供了基礎(chǔ)。

4 結(jié) 論

結(jié)合空間譜估計(jì)理論，本文提出了基于MUSIC算法的Pol-InSAR多極化干涉相位估計(jì)方法。該方法通過維納濾波理論得到各極化通道的最佳濾波器參量，構(gòu)造多極化干涉信號模型，利用包含干涉相位的信號子空間與噪聲子空間的正交性，準(zhǔn)確地估計(jì)出多極化干涉相位。同時(shí)該方法具有自適應(yīng)圖像配準(zhǔn)功能和很強(qiáng)的噪聲抑制能力。與文獻(xiàn)[2]相比，本文方法得到的圖像配準(zhǔn)精度更高，并且降低了計(jì)算量。實(shí)測數(shù)據(jù)證明了該方法的正確性和有效性。

[1] 吳一戎, 洪文, 王彥平. 極化干涉SAR的研究現(xiàn)狀和啟示[J]. 電子信息學(xué)報(bào), 2007, 29(5): 1258-1262.WU Yi-rong, HONG Wen, WANG Yan-ping. The current status and implications of polarimetric SAR interferometry[J]. Journal of Electronics & Information Technology, 2007, 29(5): 1258-1262.

[2] 李海, 廖桂生. InSAR自適應(yīng)圖像配準(zhǔn)的干涉相位估計(jì)方法[J]. 電子學(xué)報(bào), 2007, 3(35): 420-425.LI Hai, LIAO Gui-sheng. An estimation method for InSAR interferometric phase with adaptive image coregisreation [J].ACTA Electronica Sinica, 2007, 3(35): 420-425.

[3] 龔耀寰. 自適應(yīng)濾波[M]. 2版. 北京: 電子工業(yè)出版社,2003.GONG Yao-huan. Adaptive-filtering[M]. 2nd ed. Beijing:Publishing House of Electronics Industry, 2003.

[4] 何子述, 夏威. 現(xiàn)代數(shù)字信號處理及其應(yīng)用[M]. 北京:清華大學(xué)出版社, 2009.HE Zi-shu, XIA Wei. Modern digital signal processing and using[M]. Beijing: Tsinghua University Press, 2009.

[5] 皮亦鳴, 楊建宇, 付毓生, 等. 合成孔徑雷達(dá)成像原理[M]. 成都: 電子科技大學(xué)出版社, 2007.PI Yi-ming, YANG Jian-yu, FU Yu-sheng, et al.The imaging theory of Synthetic Aperture Radar[M]. Chengdu:University of Electronic Science and Technology Press,2007.

[6] ROSEN P A, HENSLEY S, JOUGHIN I R, et al. Synthetic aperture radar interferometry[J]. Proceeding of the IEEE,2000, 88(3): 333-382.

[7] BAMLER R, HARTL P. Synthetic aperture radar interferometry[J]. Inverse Problem, 1998, 14: 1-54.

[8] EICHEL P H, GHIGLIA D C. Spotlight SAR interferometry for terrain elevation mapping and interferometric change detection[R]. [S.l.]: Sandia National Labs Tech Report,1993.

[9] LANARI R, FORNARO G. Generation of digital elevation models by using SIR-C/X-SAR multifrequency two-pass interferometry: the etna case study[J]. IEEE Trans on GRS,1996, 34(5): 1097-1114.

[10] LEE J S, PATHANASSIOU K P. A new technique for nosie filtering of SAR interferometric phase images[J]. IEEE Trans on GRS, 1998, 36(5): 1456-1465.

[11] LOMBARDINI F, MONTANARI M, GINI F. Reflectivity estimation for multibaseline interferometric radar imaging of layover extended sources[J]. IEEE Trans on Signal Processing, 2003, 51(6): 1508-1519.

[12] 王永良, 陳輝, 彭應(yīng)寧, 等. 空間譜估計(jì)理論與算法[M].北京: 清華大學(xué)出版社, 2004.WANG Yong-liang, CHEN Hui, PENG Ying-ning, et al.Space spectrum estimation theory and algorithm[M].Beijing: Tsinghua University Press, 2004.

[13] 王超, 張紅, 劉智. 星載合成孔徑雷達(dá)干涉測量[M]. 北京: 科學(xué)出版社, 2002.WANG Chao, ZHANG Hong, LIU Zhi. The satellite synthetic sperture radar interferometry[M]. Beijing:Science Press, 2002.