高效UKF相位解纏算法

劉媛媛，謝先明，田憲輝，李春，曾慶寧

(1.桂林電子科技大學(xué) 信息與通信學(xué)院，廣西 桂林 541004；2.廣西科技大學(xué) 電氣與信息工程學(xué)院，廣西 柳州 545006)

0 引言

相位解纏是干涉合成孔徑雷達(dá)測(cè)量(interferometric synthetic aperture radar，InSAR)、光學(xué)干涉測(cè)量、合成孔徑聲吶(interferometric synthetic aperture sonar，InSAS)以及核磁共振成像等數(shù)據(jù)處理過(guò)程中的重要步驟之一[1]。由于三角函數(shù)的周期性，從干涉測(cè)量中獲取的表征目標(biāo)參數(shù)的干涉相位被限制在其相位主值(-π，π]區(qū)間，俗稱纏繞相位；從纏繞相位中恢復(fù)反應(yīng)目標(biāo)參數(shù)的真實(shí)干涉相位，即為所謂的相位解纏[2]。在理想情況下，通過(guò)對(duì)纏繞相位加上2π的整數(shù)倍即可實(shí)現(xiàn)相位解纏。然而，干涉測(cè)量過(guò)程中的相位跳變、相位噪聲、雷達(dá)陰影等干擾因素均會(huì)顯著增加相位解纏的難度，為此研究者們相繼提出了各種各樣方法，來(lái)解決干涉圖的相位解纏問(wèn)題。

干涉圖相位解纏算法研究至今大致可分為以下三類：路徑跟蹤法、最小范數(shù)法、基于非線性濾波的狀態(tài)估計(jì)法。路徑跟蹤法重要代表有枝切法(branch-cut,BUT)[3]、質(zhì)量引導(dǎo)法(quality guidance phase unwrapping，QGPU)[4]、網(wǎng)絡(luò)流方法[5]等算法。枝切法通過(guò)識(shí)別殘差點(diǎn)以及建立相應(yīng)的枝切線，積分時(shí)避開枝切線來(lái)實(shí)現(xiàn)相位解纏，在干涉圖相位噪聲較小、殘差點(diǎn)較少的情況下解纏速度快、精度高。然而，當(dāng)干涉圖噪聲過(guò)大、殘差點(diǎn)較多且分布較密集時(shí)，會(huì)產(chǎn)生“孤島”現(xiàn)象。質(zhì)量引導(dǎo)法則是在相位質(zhì)量圖的指導(dǎo)下沿高質(zhì)量到低質(zhì)量像元的積分路徑進(jìn)行解纏，在相位噪聲較小時(shí)可獲得較為穩(wěn)健的解纏結(jié)果；而當(dāng)干涉圖相位噪聲較大時(shí)，相位解纏過(guò)程中的誤差將難以避免地沿著解纏路徑傳遞，從而降低相位解纏精度。網(wǎng)絡(luò)流算法的典型代表為最小費(fèi)用流算法，該算法將相位解纏問(wèn)題轉(zhuǎn)化為網(wǎng)絡(luò)費(fèi)用最小化問(wèn)題，通常能從干涉相位噪聲較小的干涉圖中獲得較為穩(wěn)健的相位解纏結(jié)果，但當(dāng)干涉圖信噪比較低或相位殘差點(diǎn)較多時(shí)，該方法相位解纏精度與效率下降嚴(yán)重。最小范數(shù)法中的典型代表有基于快速傅里葉變換的無(wú)權(quán)重最小二乘法[6]、加權(quán)最小二乘法[7]等算法。這類方法將相位解纏問(wèn)題轉(zhuǎn)換成了數(shù)學(xué)上的最小范數(shù)問(wèn)題，通過(guò)在纏繞相位微分與解纏相位微分之間建立合適代價(jià)函數(shù)來(lái)求解解纏相位的估計(jì)值，這類方法通常能連續(xù)與平滑地解纏相位，但易降低解纏相位的動(dòng)態(tài)范圍，導(dǎo)致干涉圖條紋細(xì)節(jié)信息失真嚴(yán)重。基于非線性濾波的狀態(tài)估計(jì)法因其特有的抗相位噪聲能力而引起人們?cè)絹?lái)越多的關(guān)注，這類方法降低了干涉圖前置濾波對(duì)相位解纏的限制，包括擴(kuò)展卡爾曼濾波相位解纏算法[8-9]、無(wú)味卡爾曼濾波相位解纏算法(unscented Kalman filter phase unwrapping，UKFPU)[10-11]，以及粒子濾波類相位解纏算法[12-13]。其中，無(wú)味卡爾曼濾波相位解纏算法使用基于無(wú)味變換的sigma點(diǎn)來(lái)捕獲狀態(tài)變量的均值和方差，能夠減少非線性系統(tǒng)模型線性化造成的相位信息丟失，且利用路徑跟蹤策略來(lái)指導(dǎo)相位解纏程序沿高質(zhì)量區(qū)域到低質(zhì)量區(qū)域解纏干涉圖，在許多干涉圖相位解纏實(shí)例中都獲得了較好的解纏結(jié)果。然而該相位解纏算法需要在相位解纏過(guò)程不斷地搜索最佳待解纏像元，其時(shí)間耗費(fèi)較大，不利于實(shí)時(shí)性要求較高的一些應(yīng)用場(chǎng)景，因此如何保持相位解纏精度(亦或是提高相位解纏精度)的同時(shí)提高解纏效率成為亟待解決的問(wèn)題。Xie等[14]提出一種改進(jìn)的無(wú)味卡爾曼濾波相位解纏算法，是將量化路徑引導(dǎo)策略、基于相位質(zhì)量信息的像素分類策略與無(wú)味卡爾曼濾波相結(jié)合的結(jié)果，多種干涉圖相位解纏實(shí)例表明該方法在保持較好解纏精度的同時(shí)解纏效率有所提升；Xie[15]提出一種迭代無(wú)味卡爾曼濾波相位解纏算法，是將迭代無(wú)味卡爾曼濾波器、基于修正矩陣束的相位梯度估計(jì)器和基于堆排序的高效質(zhì)量引導(dǎo)策略相結(jié)合的結(jié)果，干涉圖相位解纏實(shí)例表明該方法相位解纏精度以及解纏效率都有所提升。除上述相位解纏方法外，最近有一些基于可靠性掩模的加權(quán)最小二乘相位解纏算法[16]、相位解纏的CKF(cubature Kalman filter，CKF)局部多項(xiàng)式系數(shù)遞推估計(jì)法[17]、深度學(xué)習(xí)算法[18-19]被提出對(duì)干涉圖進(jìn)行相位解纏，這些方法通常是有效的，在一些干涉圖相位解纏實(shí)例中獲得了較好的解纏結(jié)果，這些方法有利于拓展相位解纏技術(shù)的應(yīng)用領(lǐng)域。

把可靠性掩模圖、像元擴(kuò)散策略以及無(wú)味卡爾曼濾波算法結(jié)合起來(lái)，本文提出一種高效UKF相位解纏算法。首先，根據(jù)干涉圖生成相應(yīng)的枝切線分布圖；其次，利用二階差分函數(shù)計(jì)算干涉圖中各像元的質(zhì)量權(quán)值，并設(shè)置閾值將質(zhì)量權(quán)值二值化為“0”和“1”；再次，由枝切線分布圖以及二值化質(zhì)量權(quán)值矩陣生成可靠性掩模圖，該可靠性掩模圖將干涉圖分為權(quán)值為“1”的高質(zhì)量像元(非枝切線上的二值化質(zhì)量權(quán)值為“1”的像元)、權(quán)值為“0”的低質(zhì)量像元(枝切線上的像元、二值化質(zhì)量權(quán)值為“0”的像元以及枝切線圍成的閉環(huán)區(qū)域中的像元)兩部分；然后，由可靠性掩模圖確定解纏路徑，先利用UKF相位解纏程序按照像元擴(kuò)散策略解纏高質(zhì)量像元，余下未解纏像元根據(jù)已解纏像元信息，利用UKF相位解纏程序按照行(或列)的方式進(jìn)行解纏。實(shí)驗(yàn)結(jié)果表明，高效UKF相位解纏算法能夠高效與穩(wěn)健地處理干涉圖的相位解纏問(wèn)題。

1 UKF相位解纏程序方法和技術(shù)流程

利用干涉圖相鄰像元干涉相位之間的關(guān)系，干涉圖相位解纏系統(tǒng)方程可參見文獻(xiàn)[11]。文獻(xiàn)[11]提出的無(wú)味卡爾曼濾波相位解纏算法是UKF相位解纏程序與路徑跟蹤策略相結(jié)合的結(jié)果。無(wú)味卡爾曼濾波相位解纏算法利用路徑跟蹤策略來(lái)指導(dǎo)UKF相位解纏程序沿高質(zhì)量區(qū)域到低質(zhì)量區(qū)域解纏干涉圖，具有較高相位解纏精度的同時(shí)時(shí)間耗費(fèi)較大。

2 高效UKF相位解纏算法

把可靠性掩模圖、像元擴(kuò)散策略以及無(wú)味卡爾曼濾波算法相結(jié)合，本文提出一種高效UKF相位解纏算法。該方法主要包括以下4個(gè)步驟。

步驟1：根據(jù)枝切線生成原理生成干涉圖的枝切線分布圖。

1)識(shí)別干涉圖中的殘差點(diǎn)，利用最小回路積分進(jìn)行判斷并標(biāo)記位置，生成殘差點(diǎn)矩陣。

2)按照從左至右、從上至下的原則，以搜索到的第一個(gè)殘差點(diǎn)為中心。

3)首先在3×3的窗口中搜索其他殘差點(diǎn)，如果搜索到極性相反的點(diǎn)，連接兩點(diǎn)生成枝切線，并標(biāo)記該枝切線極性為0。

4)如果搜索到極性相同的點(diǎn)，連接兩點(diǎn)生成枝切線，并將極性值標(biāo)記為該枝切線極性值，并以當(dāng)前的殘差點(diǎn)為中心繼續(xù)搜索。

5)若在鄰域內(nèi)沒(méi)有搜索到殘差點(diǎn)，則擴(kuò)大搜索范圍為5×5、7×7、……，直到圖像邊界，連接最近的邊界點(diǎn)，生成枝切線。

6)重復(fù)以上過(guò)程，直至枝切線連接并平衡所有殘差點(diǎn)為止。值得注意的是，在殘差點(diǎn)連接過(guò)程中，已達(dá)到“電荷”平衡的殘差點(diǎn)不用重復(fù)連接，可減少枝切線數(shù)量。

步驟2：利用二階差分函數(shù)計(jì)算干涉圖中各像元的權(quán)值[16]，并設(shè)置閾值將該權(quán)值二值化為“0”和“1”，0表示低質(zhì)量像元，1表示高質(zhì)量像元。

(1)

(2)

式中：W為纏繞算子，將相位差值限制于(-π，π]之間；φi，j表示干涉圖(i，j)像元纏繞相位；Ri，j表示(i，j)像元的可信度。二值化質(zhì)量權(quán)值表示如式(3)所示。

(3)

式中：θ表示可信度閾值；qi，j表示(i，j)像元二值化質(zhì)量權(quán)值。

步驟3：由枝切線分布圖以及干涉圖二值化質(zhì)量權(quán)值矩陣生成可靠性掩模圖，非枝切線上的二值化質(zhì)量權(quán)值為“1”的像元為權(quán)值為“1”的高質(zhì)量像元，枝切線上的像元、二值化質(zhì)量權(quán)值為“0”的像元以及枝切線圍成的閉環(huán)區(qū)域中的像元等像元為權(quán)值為“0”的低質(zhì)量像元。

步驟4：把可靠性掩模圖、像元擴(kuò)散策略以及UKF相位解纏程序結(jié)合起來(lái)，按以下步驟對(duì)干涉圖進(jìn)行解纏。

1)首先選取非枝切線上可信度最高的像元為起點(diǎn)，其纏繞相位作為解纏相位估計(jì)值，其估計(jì)誤差方差預(yù)設(shè)為0.6；利用UKF相位解纏程序解纏4個(gè)鄰接像元中的高質(zhì)量像元，隨后將已解纏像元鄰接像元中的未解纏高質(zhì)量像元依次存儲(chǔ)在“鄰接列”中。

2)按順序取出“鄰接列”中的待解纏像元，利用UKF相位解纏程序?qū)υ撓裨M(jìn)行解纏，隨后將解纏像元鄰接像元中的未解纏高質(zhì)量像元依次存儲(chǔ)在“鄰接列”中。

3)若“鄰接列”不為空，回到步驟2)，若為空，則轉(zhuǎn)步驟4)。

4)根據(jù)已解纏像元信息，利用UKF相位解纏程序按行(或列)的方式逐一解纏余下的未解纏像元。

圖1進(jìn)一步給出了本文算法解纏步驟示意圖。圖1(a)中“黃色”像元(非枝切線上可信度最高的像元)為起始像元，“綠色”像元為高質(zhì)量像元，“紅色”像元為低質(zhì)量像元，先利用UKF相位解纏程序進(jìn)行解纏“綠色”像元；“紅色”像元將在步驟4)中進(jìn)行解纏。圖1(b)給出了執(zhí)行步驟1)之后的結(jié)果，其中，“黑色”像元為已解纏像元，“藍(lán)色”像元為存儲(chǔ)在“鄰接列”的待解纏像元，其余像元為未處理的像元。需要注意的是，若干涉圖條紋十分密集，無(wú)法利用UKF相位解纏程序從選擇的起點(diǎn)進(jìn)行解纏時(shí)，可以選擇可信度次高的像元為起點(diǎn)，按照?qǐng)D1(b)方式進(jìn)行解纏。

圖1 像元擴(kuò)散法解纏示意圖

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

為了驗(yàn)證本文算法性能，利用不同算法包括BUT、非加權(quán)迭代最小二乘法(unweighted iterative least squares，ILS)、QGPU、UKFPU，以及本文算法，在同一MATLAB軟件環(huán)境(Intel i7-6700U@3.40 GB CPU+8 GB RAM)下解纏不同的模擬干涉圖和實(shí)測(cè)干涉圖，并對(duì)各算法解纏結(jié)果進(jìn)行比較分析。

3.1 模擬數(shù)據(jù)實(shí)驗(yàn)

圖2為兩幅不同模擬干涉圖(256像素×256像素)。圖3(a)～圖3(b)分別為圖2(a)～圖2(b)含噪聲的纏繞相位圖，其信噪比依次為3.01 dB、2.18 dB。依次用BUT、ILS、QGPU、UKFPU以及本文算法對(duì)上述兩幅干涉圖進(jìn)行解纏，解纏結(jié)果分別如圖4～圖8、圖9～圖13所示。

圖2 模擬干涉圖

圖3 含噪纏繞相位圖

圖4 BUT算法解纏相關(guān)處理結(jié)果

圖5 ILS算法解纏結(jié)果

圖6 QGPU算法解纏結(jié)果

圖7 UKFPU算法解纏結(jié)果

圖8 本文算法解纏結(jié)果

圖9 BUT算法解纏相關(guān)處理結(jié)果

圖10 ILS算法解纏結(jié)果

圖11 QGPU算法解纏結(jié)果

圖12 UKFPU算法解纏結(jié)果

圖13 本文算法解纏結(jié)果

圖4、圖9第一行從左至右分別代表圖3(a)～圖3(b)的纏繞相位圖的殘差點(diǎn)、枝切線、枝切法解纏結(jié)果，圖5～圖8、圖10～圖13的左列、中間列和右列分別表示ILS、QGPU、UKFPU、本文算法解纏相位結(jié)果、解纏相位誤差和解纏相位誤差直方圖。表1為各算法解纏不同信噪比干涉圖的均方根誤差。表1中的“-”表示隨著噪聲增大，枝切法形成一些閉環(huán)區(qū)域，這些閉環(huán)區(qū)域無(wú)法解纏，因此無(wú)法統(tǒng)計(jì)均方根誤差。由圖4～圖8及圖9～圖13可以看出：由于干涉圖信噪比較低，BUT算法解纏相位中存在部分無(wú)法解纏的區(qū)域；ILS算法雖然解纏結(jié)果較為平滑，由其解纏相位誤差圖、解纏相位誤差直方圖以及表1可知其解纏誤差較大；QGPU算法解纏過(guò)程中誤差沿著積分路徑傳遞，導(dǎo)致部分區(qū)域誤差較大；而UKFPU以及本文算法相較于BUT、ILS、QGPU算法具有更穩(wěn)健的相位解纏能力，且UKFPU與本文算法解纏精度相當(dāng)。表2列出了上述算法解纏上述干涉圖的平均運(yùn)行時(shí)間。可以看出，本文算法解纏時(shí)間遠(yuǎn)小于QGPU、UKFPU算法。因此，本文算法在保持穩(wěn)健的相位解纏能力的同時(shí)，其時(shí)間消耗亦是可以接受的。

表1 不同信噪比下不同算法均方根誤差

表2 各算法解纏時(shí)間 s

3.2 實(shí)測(cè)數(shù)據(jù)實(shí)驗(yàn)

圖14為經(jīng)中值濾波處理過(guò)后截取的部分Enta火山干涉圖，各算法解纏結(jié)果如圖15所示，其解纏相位重纏繞結(jié)果如圖16所示。

圖14 局部Enta火山干涉圖

圖15 各算法解纏結(jié)果

圖16 各算法解纏相位重纏繞結(jié)果

圖17為三峽實(shí)測(cè)干涉圖，各算法解纏結(jié)果如圖18所示，解纏相位重纏繞結(jié)果如圖19所示。

由圖15(a)以及圖18(a)可以看出，BUT算法解纏相位圖存在部分無(wú)法解纏以及明顯不一致的區(qū)域。由圖15(b)以及圖18(b)可見，ILS算法解纏結(jié)果較為光滑，但根據(jù)圖16(a)以及圖19(a)可以發(fā)現(xiàn)，ILS算法重纏繞結(jié)果相較于圖14以及圖17原始干涉圖條紋差別較大，大量干涉條紋細(xì)節(jié)信息丟

圖17 三峽實(shí)測(cè)干涉圖

圖18 各算法解纏結(jié)果

圖19 各算法解纏相位重纏繞結(jié)果

失，故該方法解纏結(jié)果精度有限。而由圖15(c)～圖15(e)以及圖18(c)～圖18(e)可知，QGPU、UKFPU、本文算法解纏結(jié)果較好，且由圖16(b)～圖16(d)以及圖19(b)～圖19(d)可見，其重纏繞干涉圖條紋與原始干涉圖條紋基本一致。與此同時(shí)，不同于QGPU算法，UKFPU算法以及本文算法重纏繞相位圖中殘留散斑噪聲較少，這表明這兩種算法能在相位解纏的同時(shí)有效去除干涉圖中相位噪聲，有利于獲得更高的解纏精度。表3給出了上述算法的運(yùn)行時(shí)間，可知本文算法解纏時(shí)間遠(yuǎn)小于QGPU、UKFPU算法。

表3 各算法解纏時(shí)間 s

4 結(jié)束語(yǔ)

本文提出一種高效UKF相位解纏算法，該算法是可靠性掩模圖、像元擴(kuò)散策略以及無(wú)味卡爾曼濾波算法相結(jié)合的結(jié)果。模擬數(shù)據(jù)實(shí)驗(yàn)以及實(shí)測(cè)數(shù)據(jù)實(shí)驗(yàn)結(jié)果證明了高效UKF相位解纏算法有效性，相對(duì)于傳統(tǒng)算法如BUT、ILS以及QGPU，具有更穩(wěn)健的相位解纏能力，且在時(shí)間耗費(fèi)較少的同時(shí)其解纏精度與UKFPU算法相當(dāng)。然而，所提出的方法處理復(fù)雜地形仍有改進(jìn)空間，且相位解纏時(shí)間消耗還有優(yōu)化的空間。努力尋找高效穩(wěn)健的相位解纏技術(shù)是我們未來(lái)努力的目標(biāo)。