代高興，謝先明

桂林電子科技大學(xué)信息與通信學(xué)院，廣西 桂林 541004

基于修正嵌入式容積卡爾曼濾波的相位展開算法

代高興，謝先明

桂林電子科技大學(xué)信息與通信學(xué)院，廣西 桂林 541004

針對干涉圖展開問題，將Levenberg-Marquardt方法修正后的嵌入式容積卡爾曼濾波器應(yīng)用于纏繞相位圖像的展開中，結(jié)合基于修正矩陣束模型的局部相位梯度估計算法以及量化跟蹤策略，提出一種基于修正嵌入式容積卡爾曼濾波的相位展開算法。該算法不僅可以精確和快速地展開纏繞像元，還可以在展開纏繞像元的同時抑制相位噪聲，降低前置預(yù)濾波器的難度與復(fù)雜度，甚至可以在處理受噪聲污染不嚴(yán)重的干涉圖時免去前置預(yù)濾波處理步驟。試驗結(jié)果表明本文算法具有較高的效率和良好的穩(wěn)鍵性。

相位展開；嵌入式容積卡爾曼濾波器；Levenberg-Marquardt方法；量化跟蹤策略

干涉合成孔徑雷達(dá)(interferometric synthetic aperture radar,InSAR)可以高精度、高可靠性地獲取地表三維信息和高程變化信息，被廣泛應(yīng)用于大地測繪、海洋監(jiān)測、火山監(jiān)測和地震檢測等領(lǐng)域。干涉相位展開是lnSAR技術(shù)應(yīng)用中尤為關(guān)鍵的環(huán)節(jié)，其相位展開精度直接影響著InSAR系統(tǒng)高程測量的精度，故干涉圖展開算法的研究至關(guān)重要。

干涉圖展開算法研究至今大致可分為路徑跟蹤法[1-11]、最小范數(shù)法[12-13]、最優(yōu)估計法[14-23]三大類。其中經(jīng)典相位展開算法，包括枝切法[1]、質(zhì)量圖引導(dǎo)法[3-5]、網(wǎng)絡(luò)流法[7-9]、最小二乘法[12]等。為了得到較為理想的展開結(jié)果，通常要在相位展開之前使用前置預(yù)濾波器抑制干涉圖中的相位噪聲，即使用濾波器，如均值濾波器、Goldstein濾波器等，對干涉圖展開前濾波，但在盡可能地抑制相位噪聲的同時不免會模糊掉干涉條紋的邊緣特性?；跀?shù)據(jù)融合的最優(yōu)估計相位展開算法，包括擴(kuò)展卡爾曼濾波相位展開算法(EKFPU)[14-15]、不敏卡爾曼濾波相位展開算法(UKFPU)[16-20]、粒子濾波相位展開算法[21-22]等，利用非線性濾波器自身具有的噪聲抑制能力，可以在展開纏繞像元的同時抑制相位噪聲，在一定程度上彌補(bǔ)了經(jīng)典相位展開算法的不足，拓展了InSAR技術(shù)應(yīng)用的范圍。

為進(jìn)一步促進(jìn)基于數(shù)據(jù)融合的相位展開技術(shù)發(fā)展，本文將經(jīng)Levenberg-Marquardt方法[24-25]修正后的嵌入式容積卡爾曼濾波器[26-27](MECKF)、基于修正矩陣束模型(AMPM)的局部相位梯度估計算法[20]以及量化跟蹤策略[28]結(jié)合起來，提出一種新的基于數(shù)據(jù)融合的最優(yōu)估計相位展開算法——修正嵌入式容積卡爾曼濾波的相位展開(modified embedded cubature Kalman filter phase unwrapping,MECKFPU)。MECKFPU算法利用MECKF同時執(zhí)行相位展開和噪聲抑制，利用AMPM局部相位梯度估計算法快速和精確地從復(fù)干涉圖中獲取相位梯度信息，以及利用量化跟蹤策略從高質(zhì)量像元到低質(zhì)量像元的路徑快速地搜索最佳待展開像元，從而保證了MECKFPU精確和快速地展開纏繞像元，具有比經(jīng)典相位展開算法更穩(wěn)健的相位展開能力，且相比于同類型的EKFPU、UKFPU等算法，MECKFPU擁有更佳的相位展開效率和精度，能夠較好地解決干涉圖展開問題。

1 MECKFPU算法

1.1 MECKF相位展開算法原理

利用干涉圖中相鄰干涉相位之間的關(guān)系，將歸一化復(fù)干涉的同向分量和正交分量分別作為干涉相位的兩個觀測值，在沿某一確定路徑下，相位展開系統(tǒng)模型如下[20]

(1)

式中，x(m,n)和x(a,s)分別表示(m,n)像元和(a,s)像元的真實干涉相位；g(m,n)|(a,s)表示干涉圖中(m,n)像元與(a,s)像元之間的相位梯度估計值，本文通過AMPM局部相位梯度估計算法[20]獲得；ζ(m,n)|(a,s)表示干涉圖中(m,n)像元與(a,s)像元之間的相位梯度估計誤差；v1,(m,n)和v2,(m,n)分別表示附加在復(fù)干涉正交分量和同向分量上的噪聲。

MECKF相位展開算法是在質(zhì)量圖引導(dǎo)策略指導(dǎo)下，同時利用待展開像元相鄰8個像元中已展開像元信息，沿高質(zhì)量區(qū)域到低質(zhì)量區(qū)域的路徑完成對纏繞像元遞推估計。

嵌入式容積卡爾曼濾波器采用三階嵌入式容積準(zhǔn)則，利用附有權(quán)值的積分點實現(xiàn)對非線性系統(tǒng)后驗均值和方差的逼近，具有良好的狀態(tài)估計精度和數(shù)值穩(wěn)定性。相應(yīng)積分點和權(quán)值計算如下

(2)

式中，nx表示狀態(tài)變量維數(shù)，本文nx為1；[δ]i表示[δ]的第i列，其δ及[δ]的詳細(xì)取值請參見文獻(xiàn)[26—27]，本文δ由經(jīng)驗取值為0.5。

根據(jù)相位展開系統(tǒng)模型，MECKF相位展開算法預(yù)測過程如下

(3)

式中，(m,n)和(a,s)分別表示待展開像元以及待展開像元相鄰8個像元中的已展開像元；γ(a,s)表示干涉圖(a,s)像元的相干系數(shù)(或偽相干系數(shù))；ψ表示(m,n)像元相鄰8個像元的集合；P(m,n)|(a,s)表示由(a,s)像元的估計誤差方差P(a,s)通過權(quán)系數(shù)d(a,s)加權(quán)得到的(m,n)像元的預(yù)測估計誤差方差的預(yù)測值；γi,(m,n)表示(m,n)像元積分點的預(yù)測值；Pxx,(m,n)表示狀態(tài)預(yù)測值X(m,n)的預(yù)測估計誤差方差；Q(m,n)|(a,s)表示g(m,n)|(a,s)的相位梯度估計誤差方差[20]。

MECKF相位展開算法利用Levenberg-Marquardt方法優(yōu)化預(yù)測過程中的預(yù)測估計誤差方差，以提高算法收斂性，減小待展開像元的展開誤差

(4)

式中，μ表示優(yōu)化Pxx,(m,n)的調(diào)節(jié)參數(shù)[24-25]，本文由經(jīng)驗取值為0.3；I表示nx維單位矩陣。

MECKF相位展開算法更新過程如下

(5)

式中，y(m,n)和Y(m,n)分別表示(m,n)像元的觀測值和觀測預(yù)測值；R(m,n)表示(m,n)像元的觀測預(yù)測協(xié)方差[20]；κ(m,n)表示(m,n)像元的增益矩陣；x(m,n)和P(m,n)分別表示(m,n)像元的狀態(tài)估計和狀態(tài)估計誤差方差。

1.2 結(jié)合量化跟蹤策略的二維MECKF相位展開算法(MECKFPU)

為減少M(fèi)ECKF相位展開算法搜索最佳待展開像元的時間消耗，引入文獻(xiàn)[28]提出的量化跟蹤策略來指導(dǎo)相位展開路徑。量化跟蹤策略通過歸一化和量化路徑引導(dǎo)圖，利用附有鏈表的優(yōu)先隊列快速排序算法實現(xiàn)最佳待展開像元的快速搜索工作。利用干涉圖偽相干系數(shù)圖和微分偏差圖構(gòu)造指導(dǎo)MECKF相位展開算法展開路徑的路徑引導(dǎo)圖[19],如下

(6)

式中，D(m,n)表示像元(m,n)的微分偏差；β表示調(diào)整參數(shù)，詳細(xì)請參見文獻(xiàn)[19]，本文由經(jīng)驗取值為0.2；normalization[·]表示歸一化處理；round{·}表示四舍五入取整操作；M(m,n)表示(m,n)像元經(jīng)歸一化和量化處理后的量化路徑引導(dǎo)值；N表示量化等級，本文取500。MECKFPU算法流程如圖1。

圖1 MECKFPU算法流程Fig.1 Flow chart of MECKFPU

2 試驗結(jié)果與分析

2.1 模擬數(shù)據(jù)試驗

2.1.1 結(jié)合預(yù)濾波處理的金字塔地形干涉圖展開試驗

模擬干涉圖主要參數(shù)：軌道高度為590 km，基線長度為610 m，下視角為45°，波長為0.04 m，基線傾角為10°，地面分辨率為6×6 m，場景為380 m的金字塔地形(256×256)，干涉圖見圖2。圖3為MECKFPU展開經(jīng)均值濾波器(3×3)處理信噪比為3.01 dB圖2(c)后的結(jié)果?？梢钥闯觯琈ECKFPU的展開誤差絕大分集中在[-0.5,0.5]。表1列出了枝切法、質(zhì)量圖引導(dǎo)法、迭代最小二乘法、網(wǎng)絡(luò)流法、EKFPU、UKFPU和MECKFPU處理不同信噪比干涉圖時的均方根誤差(MSRE)。可以直觀地看出，MECKFPU的展開精度最高。表2列出了相同運(yùn)算環(huán)境條件下上述各算法展開圖2(c)的運(yùn)行時間，顯然，MECKFPU時間消耗遠(yuǎn)遠(yuǎn)小于質(zhì)量圖引導(dǎo)法、EKFPU以及UKFPU算法，即使與枝切法、迭代最小二乘法以及網(wǎng)絡(luò)流法等算法相比，MECKFPU時間消耗也是可以接受的。

2.1.2 無預(yù)濾波處理的金字塔地形干涉圖展開試驗

為了進(jìn)一步分析MECKFPU的性能，直接用MECKFPU展開圖2(b)，結(jié)果見圖4?？梢奙ECKFPU對于不經(jīng)預(yù)濾波處理的干涉圖也可以得到理想的展開結(jié)果。表3列出了質(zhì)量圖引導(dǎo)法、迭代最小二乘法、網(wǎng)絡(luò)流法、EKFPU、UKFPU和MECKFPU處理不同信噪比(無預(yù)濾波處理的)干涉圖時的MSRE?？梢?，相比于表1，質(zhì)量圖引導(dǎo)法和EKFPU的展開精度明顯降低，且干涉圖信噪比越低其展開精度越低；迭代最小二乘法和網(wǎng)絡(luò)流法不僅展開精度下降明顯，而且算法較不穩(wěn)定；UKFPU和MECKFPU的展開精度稍有變化，但無明顯下降，且在干涉圖信噪比較低的情況下，仍有較為理想的展開精度，遠(yuǎn)遠(yuǎn)高于其他算法的展開精度，其中MECKFPU的展開精度更高。由此可得，即使對于較低信噪比金字塔地形干涉圖，MECKFPU仍有較高的展開精度。

2.1.3 復(fù)雜多山地形干涉圖展開試驗

用多山地形場景取代上一節(jié)的金字塔地形場景，即可得復(fù)雜多山地形干涉圖，見圖5。圖6為MECKFPU對信噪比為1.42 dB圖5(b)的處理結(jié)果?？梢?，MECKFPU的展開誤差大部分集中在0附近。選取金字塔地形干涉圖展開試驗中表現(xiàn)出較好能力的UKFPU和MECKFPU做進(jìn)一步比較，直接展開不經(jīng)預(yù)濾波處理的不同信噪比復(fù)雜多山地形干涉圖，其MSRE在表4中列出，可看出，MECKFPU展開精度更高。綜上可得，MECKFPU處理低信噪比復(fù)雜多山地形干涉圖仍有較高的展開精度。

圖2 金字塔地形干涉圖Fig.2 Synthetic interferogram over pyramid

圖3 MECKFPU處理經(jīng)預(yù)濾波的金字塔地形干涉圖結(jié)果Fig.3 The MECKFPU solution of pre-filtered interferogram over pyramid

表2 各算法運(yùn)行時間

圖4 MECKFPU處理無預(yù)濾波的金字塔地形干涉圖結(jié)果Fig.4 The MECKFPU solution of without-filtered interferogram over pyramid

圖5 復(fù)雜多山地形干涉圖Fig.5 Synthetic interferogram over peaks

Tab.4Phaseunwrappingaccuracyofdifferentmethodswithoutpre-filteringprocedure

rad

2.2 實測數(shù)據(jù)試驗

在實測數(shù)據(jù)試驗中，用MECKFPU算法展開經(jīng)均值濾波器(3×3)處理后的局部Etna火山干涉圖(256×256)，見圖7?？梢姡归_相位平滑連續(xù)，重纏繞相位圖干涉條紋與圖7(a)中干涉條紋保持一致，且干涉條紋清晰，幾乎不存在相位噪聲。這表明MECKFPU算法有效地抑制干涉圖中的相位噪聲，并獲得了較好展開結(jié)果。綜上所述，對于實際地形干涉圖，MECKFPU仍有理想的展開效果和展開效率，具有良好的魯棒性。

3 結(jié) 論

本文提出的MECKFPU是將經(jīng)Levenberg-Marquardt方法修正后的嵌入式容積卡爾曼濾波器、AMPM局部相位梯度估計算法以及量化跟蹤策略相結(jié)合的結(jié)果，并在模擬與實測數(shù)據(jù)處理試驗中驗證其有效性。本文算法與部分常用經(jīng)典算法包括枝切法、質(zhì)量圖引導(dǎo)法、迭代最小二乘法、網(wǎng)絡(luò)流法、EKFPU和UKFPU等算法在相關(guān)數(shù)據(jù)處理試驗中進(jìn)行了比較。結(jié)果表明，MECKFPU可以在快速和精確地展開干涉圖的同時抑制相位噪聲，具有較高的效率和良好的穩(wěn)健性，能夠有效地解決干涉圖展開問題。

圖7 MECKFPU展開結(jié)果Fig.7 The MECKFPU solution

[1] GOLDSTEIN R M，ZEBKER H A，WERNER C L.Satellite Radar Interferometry：Two-dimensional Phase Unwrapping[J].Radio Science，1988，23(4)：713-720.

[2] ASUNDI A，ZHOU Wensen.Fast Phase-unwrapping Algorithm Based on a Gray-scale Mask and Flood Fill[J].Applied Optics，1998，37(23)：5416-5420.

[3] ZHAO Ming，HUANG Lei，ZHANG Qican，et al.Quality-guided Phase Unwrapping Technique：Comparison of Quality Maps and Guiding Strategies[J].Applied Optics，2011，50(33)：6214-6224.

[4] FANG Suping，MENG Lei，WANG Leijie，et al.Quality-guided Phase Unwrapping Algorithm Based on Reliability Evaluation[J].Applied Optics，2011，50(28)：5446-5452.

[5] ZHONG Heping，TANG Jinsong，ZHANG Sen，et al.An Improved Quality-guided Phase-unwrapping Algorithm Based on Priority Queue[J].IEEE Geoscience and Remote Sensing Letters，2011，8(2)：364-368.

[6] GHIGLIA D C，PRITT M D.Two-dimensional Phase Unwrapping：Theory，Algorithms，and Software[M].New York：John Wiley and Sons Inc，1998.

[7] CHEN C W.Statistical-cost Network-flow Approaches to Two-dimensional Phase Unwrapping for Radar Interferometry[D].Stanford，California：Stanford University，2001.

[8] COSTANTINI M.A Novel Phase Unwrapping Method Based on Network Programming[J].IEEE Transactions on Geoscience and Remote Sensing，1998，36(3)：813-821.

[9] CHEN C W，ZEBKER H A.Network Approaches to Two-dimensional Phase Unwrapping：Intractability and Two New Algorithms[J].Journal of the Optical Society of America A，2000，17(3)：401-414.

[10] GAO Dapeng，YIN Fuliang.Mask Cut Optimization in Two-dimensional Phase Unwrapping[J].IEEE Geoscience and Remote Sensing Letters，2012，9(3)：338-342.

[11] FLYNN T J.Two-dimensional Phase Unwrapping with Minimum Weighted Discontinuity[J].Journal of the Optical Society of America A，1997，14(10)：2692-2701.

[12] GHIGLIA D C，ROMERO L A.Robust Two-dimensional Weighted and Unweighted Phase Unwrapping that Uses Fast Transforms and Iterative Methods[J].Journal of the Optical Society of America A，1994，11(1)：107-117.

[13] 陳強(qiáng)，楊瑩輝，劉國祥，等.基于邊界探測的InSAR最小二乘整周相位解纏方法[J].測繪學(xué)報，2012，41(3)：441-448.

CHEN Qiang，YANG Yinghui，LIU Guoxiang，et al.InSAR Phase Unwrapping Using Least Squares Method with Integer Ambiguity Resolution and Edge Detection[J].Acta Geodaetica et Cartographica Sinica，2012，41(3)：441-448.

[14] LOFFELD O，NIES H，KNEDLIK S，et al.Phase Unwrapping for SAR Interferometry: A Data Fusion Approach by Kalman Filtering[J].IEEE Transactions on Geoscience and Remote Sensing，2008，46(1)：47-58.

[15] NIES H，LOFFELD O，WANG R.Phase Unwrapping Using 2D-Kalman Filter: Potential and Limitations[C]∥Proceedings of 2008 IEEE International Geoscience and Remote Sensing Symposium.Boston，MA：IEEE，2008：IV-1213-IV-1216.

[16] WAGHMARE R G，MISHRA D，SAI SUBRAHMANYAM G R，et al.Signal Tracking Approach for Phase Estimation in Digital Holographic Interferometry[J].Applied Optics，2014，53(19)：4150-4157.

[17] 謝先明.結(jié)合濾波算法的不敏卡爾曼濾波器相位解纏方法[J].測繪學(xué)報，2014，43(7)：739-745.DOI：10.13485/j.cnki.11-2089.2014.0102.

XIE Xianming.An UKF Phase Unwrapping Algorithm with a Pre-filtering Procedure[J].Acta Geodaetica et Cartographica Sinica，2014，43(7)：739-745.DOI：10.13485/j.cnki.11-2089.2014.0102.

[18] CHENG Zhongtao，LIU Dong，YANG Yongying，et al.Practical Phase Unwrapping of Interferometric Fringes Based on Unscented Kalman Filter Technique[J].Optics Express，2015，23(25)：32337-32349.

[19] XIE Xianming，ZENG Qingning.Efficient and Robust Phase Unwrapping Algorithm Based on Unscented Kalman Filter，the Strategy of Quantizing Paths-guided Map，and Pixel Classification Strategy[J].Applied Optics，2015，54(31)：9294-9307.

[20] XIE Xianming，LI Yinghui.Enhanced Phase Unwrapping Algorithm Based on Unscented Kalman Filter，Enhanced Phase Gradient Estimator，and Path-Following Strategy[J].Applied Optics，2014，53(18)：4049-4060.

[21] MARTINEZ-ESPLA J J，MARTINEZ-MARIN T，LOPEZ-SANCHEZ J M.A Particle Filter Approach for InSAR Phase Filtering and Unwrapping[J].IEEE Transactions on Geoscience and Remote Sensing，2009，47(4)：1197-1211.

[22] WAGHMARE R G，SUKUMAR P R，SUBRAHMANYAM G R K S，et al.Particle-filter-based Phase Estimation in Digital Holographic Interferometry[J].Journal of the Optical Society of America A，2016，33(3)：326-332.

[23] KULKARNI R，RASTOGI P.Simultaneous Estimation of Phase Derivative and Phase Using Parallel Kalman Filter Implementation[J]. Measurement Science and Technology，2016，27(6)：065203.

[24] MORé J J.The Levenberg-Marquardt Algorithm：Implementation and Theory[M]∥WATSON G A.Numerical Analysis.Berlin：Springer，1978：105-116.

[25] 侯代文，殷福亮.基于迭代中心差分卡爾曼濾波的說話人跟蹤方法[J].電子與信息學(xué)報，2008，30(7)：1684-1689.

HOU Daiwen，YIN Fuliang.Iterated Central Difference Kalman Filter Based Speaker Tracking[J].Journal of Electronics & Information Technology，2008，30(7)：1684-1689.

[26] 張鑫春，郭承軍.均方根嵌入式容積卡爾曼濾波[J].控制理論與應(yīng)用，2013，30(9)：1116-1121.

ZHANG Xinchun，GUO Chengjun.Square-root Imbedded Cubature Kalman Filtering[J].Control Theory & Applications，2013，30(9)：1116-1121.

[27] 張鑫春.INS/GNSS深組合導(dǎo)航系統(tǒng)的非線性研究[D].成都：電子科技大學(xué)，2014.

ZHANG Xinchun.Research on the Nonlinearity of INS/GNSS Deeply Integration[D].Chengdu：University of Electronic Science and Technology of China，2014.

[28] 鐘何平，唐勁松，張森，等.利用量化質(zhì)量圖和優(yōu)先隊列的快速相位解纏算法[J].武漢大學(xué)學(xué)報(信息科學(xué)版)，2011，36(3)：342-345.

ZHONG Heping，TANG Jinsong，ZHANG Sen，et al.A Fast Phase Unwrapping Algorithm Based on Quantized Quality Map and Priority Queue[J].Geomatics and Information Science of Wuhan University，2011，36(3)：342-354.

A Fresh Phase Unwrapping Method Based on Modified Embedded Cubature Kalman Filter

DAI Gaoxing,XIE Xianming

School of Information and Communication，Guilin University of Electronic Technology，Guilin 541004，China

A phase unwrapping algorithm is proposed that based on modified embedded cubature Kalman filter for interferometric fringes.The algorithm is the result through combining an embedded cubature Kalman filter modified by a Levenberg-Marquardt method，a robust phase gradient estimator based on amended matrix pencil model，and a quantization path-following strategy.This method can accurately estimate unambiguous unwrapped phase of interferometric fringes by applying the modified embedded cubature Kalman filter to perform phase unwrapping and noise suppress simultaneously along the path routed by the quantization path-following strategy，which is beneficial to simplify the complexity and the difficulty of pre-filtering procedure followed by phase unwrapping procedure，and even can remove the pre-filtering procedure.Results obtained with synthetic data and real data show more acceptable solutions with the proposed method，compared with some of the most used algorithms.

phase unwrapping； modified embedded cubature Kalman filter； Levenberg-Marquardt method；quantization path-following strategy

The National Natural Science Foundation of China(Nos. 41661092； 61461011)； The State Key Program of National Natural Science of Guangxi Zhuang Autonomous Region (No. 2016GXNSFDA380009)； The Natural Science Foundation of Guangxi Zhuang Autonomous Region (No. 2014GXNSFBA118273)； The Guangxi Wireless Broad Band Communication and Signal Processing Key Laboratory 2014/2015 Director Fund Project (No. GXKL061503)； 2017 Postgraduate Research and Innovation Project of Guilin University of Electronic Technology(No. 2017YJCX23)

DAI Gaoxing (1991—)，male，master of science，majors in InSAR phase unwrapping.

XIE Xianming

E-mail： xxmxgm@163.com

代高興，謝先明.基于修正嵌入式容積卡爾曼濾波的相位展開算法[J].測繪學(xué)報，2017,46(12)：1998-2005.

10.11947/j.AGCS.2017.20160665.

DAI Gaoxing,XIE Xianming.A Fresh Phase Unwrapping Method Based on Modified Embedded Cubature Kalman Filter[J]. Acta Geodaetica et Cartographica Sinica,2017,46(12):1998-2005. DOI:10.11947/j.AGCS.2017.20160665.

P236

A

1001-1595(2017)12-1998-08

國家自然科學(xué)基金(41661092；61461011)；廣西自然科學(xué)基金重點項目(2016GXNSFDA380009)；廣西自然科學(xué)基金(2014GXNSFBA118273)；廣西無線寬帶通信與信號處理重點實驗室2014/2015年主任基金(GXKL061503)；2017年桂林電子科技大學(xué)研究生科研創(chuàng)新項目(2017YJCX23)

叢樹平)

2016-12-26

2017-07-19

代高興(1991—)，男，碩士，研究方向為InSAR相位解纏。

E-mail： 1984393901@qq.com

謝先明