劉媛媛,謝先明,田憲輝,李春,曾慶寧
(1.桂林電子科技大學 信息與通信學院,廣西 桂林 541004;2.廣西科技大學 電氣與信息工程學院,廣西 柳州 545006)
相位解纏是干涉合成孔徑雷達測量(interferometric synthetic aperture radar,InSAR)、光學干涉測量、合成孔徑聲吶(interferometric synthetic aperture sonar,InSAS)以及核磁共振成像等數據處理過程中的重要步驟之一[1]。由于三角函數的周期性,從干涉測量中獲取的表征目標參數的干涉相位被限制在其相位主值(-π,π]區間,俗稱纏繞相位;從纏繞相位中恢復反應目標參數的真實干涉相位,即為所謂的相位解纏[2]。在理想情況下,通過對纏繞相位加上2π的整數倍即可實現相位解纏。然而,干涉測量過程中的相位跳變、相位噪聲、雷達陰影等干擾因素均會顯著增加相位解纏的難度,為此研究者們相繼提出了各種各樣方法,來解決干涉圖的相位解纏問題。
干涉圖相位解纏算法研究至今大致可分為以下三類:路徑跟蹤法、最小范數法、基于非線性濾波的狀態估計法。路徑跟蹤法重要代表有枝切法(branch-cut,BUT)[3]、質量引導法(quality guidance phase unwrapping,QGPU)[4]、網絡流方法[5]等算法。枝切法通過識別殘差點以及建立相應的枝切線,積分時避開枝切線來實現相位解纏,在干涉圖相位噪聲較小、殘差點較少的情況下解纏速度快、精度高。然而,當干涉圖噪聲過大、殘差點較多且分布較密集時,會產生“孤島”現象。質量引導法則是在相位質量圖的指導下沿高質量到低質量像元的積分路徑進行解纏,在相位噪聲較小時可獲得較為穩健的解纏結果;而當干涉圖相位噪聲較大時,相位解纏過程中的誤差將難以避免地沿著解纏路徑傳遞,從而降低相位解纏精度。網絡流算法的典型代表為最小費用流算法,該算法將相位解纏問題轉化為網絡費用最小化問題,通常能從干涉相位噪聲較小的干涉圖中獲得較為穩健的相位解纏結果,但當干涉圖信噪比較低或相位殘差點較多時,該方法相位解纏精度與效率下降嚴重。最小范數法中的典型代表有基于快速傅里葉變換的無權重最小二乘法[6]、加權最小二乘法[7]等算法。這類方法將相位解纏問題轉換成了數學上的最小范數問題,通過在纏繞相位微分與解纏相位微分之間建立合適代價函數來求解解纏相位的估計值,這類方法通常能連續與平滑地解纏相位,但易降低解纏相位的動態范圍,導致干涉圖條紋細節信息失真嚴重。基于非線性濾波的狀態估計法因其特有的抗相位噪聲能力而引起人們越來越多的關注,這類方法降低了干涉圖前置濾波對相位解纏的限制,包括擴展卡爾曼濾波相位解纏算法[8-9]、無味卡爾曼濾波相位解纏算法(unscented Kalman filter phase unwrapping,UKFPU)[10-11],以及粒子濾波類相位解纏算法[12-13]。其中,無味卡爾曼濾波相位解纏算法使用基于無味變換的sigma點來捕獲狀態變量的均值和方差,能夠減少非線性系統模型線性化造成的相位信息丟失,且利用路徑跟蹤策略來指導相位解纏程序沿高質量區域到低質量區域解纏干涉圖,在許多干涉圖相位解纏實例中都獲得了較好的解纏結果。然而該相位解纏算法需要在相位解纏過程不斷地搜索最佳待解纏像元,其時間耗費較大,不利于實時性要求較高的一些應用場景,因此如何保持相位解纏精度(亦或是提高相位解纏精度)的同時提高解纏效率成為亟待解決的問題。Xie等[14]提出一種改進的無味卡爾曼濾波相位解纏算法,是將量化路徑引導策略、基于相位質量信息的像素分類策略與無味卡爾曼濾波相結合的結果,多種干涉圖相位解纏實例表明該方法在保持較好解纏精度的同時解纏效率有所提升;Xie[15]提出一種迭代無味卡爾曼濾波相位解纏算法,是將迭代無味卡爾曼濾波器、基于修正矩陣束的相位梯度估計器和基于堆排序的高效質量引導策略相結合的結果,干涉圖相位解纏實例表明該方法相位解纏精度以及解纏效率都有所提升。除上述相位解纏方法外,最近有一些基于可靠性掩模的加權最小二乘相位解纏算法[16]、相位解纏的CKF(cubature Kalman filter,CKF)局部多項式系數遞推估計法[17]、深度學習算法[18-19]被提出對干涉圖進行相位解纏,這些方法通常是有效的,在一些干涉圖相位解纏實例中獲得了較好的解纏結果,這些方法有利于拓展相位解纏技術的應用領域。
把可靠性掩模圖、像元擴散策略以及無味卡爾曼濾波算法結合起來,本文提出一種高效UKF相位解纏算法。首先,根據干涉圖生成相應的枝切線分布圖;其次,利用二階差分函數計算干涉圖中各像元的質量權值,并設置閾值將質量權值二值化為“0”和“1”;再次,由枝切線分布圖以及二值化質量權值矩陣生成可靠性掩模圖,該可靠性掩模圖將干涉圖分為權值為“1”的高質量像元(非枝切線上的二值化質量權值為“1”的像元)、權值為“0”的低質量像元(枝切線上的像元、二值化質量權值為“0”的像元以及枝切線圍成的閉環區域中的像元)兩部分;然后,由可靠性掩模圖確定解纏路徑,先利用UKF相位解纏程序按照像元擴散策略解纏高質量像元,余下未解纏像元根據已解纏像元信息,利用UKF相位解纏程序按照行(或列)的方式進行解纏。實驗結果表明,高效UKF相位解纏算法能夠高效與穩健地處理干涉圖的相位解纏問題。
利用干涉圖相鄰像元干涉相位之間的關系,干涉圖相位解纏系統方程可參見文獻[11]。文獻[11]提出的無味卡爾曼濾波相位解纏算法是UKF相位解纏程序與路徑跟蹤策略相結合的結果。無味卡爾曼濾波相位解纏算法利用路徑跟蹤策略來指導UKF相位解纏程序沿高質量區域到低質量區域解纏干涉圖,具有較高相位解纏精度的同時時間耗費較大。
把可靠性掩模圖、像元擴散策略以及無味卡爾曼濾波算法相結合,本文提出一種高效UKF相位解纏算法。該方法主要包括以下4個步驟。
步驟1:根據枝切線生成原理生成干涉圖的枝切線分布圖。
1)識別干涉圖中的殘差點,利用最小回路積分進行判斷并標記位置,生成殘差點矩陣。
2)按照從左至右、從上至下的原則,以搜索到的第一個殘差點為中心。
3)首先在3×3的窗口中搜索其他殘差點,如果搜索到極性相反的點,連接兩點生成枝切線,并標記該枝切線極性為0。
4)如果搜索到極性相同的點,連接兩點生成枝切線,并將極性值標記為該枝切線極性值,并以當前的殘差點為中心繼續搜索。
5)若在鄰域內沒有搜索到殘差點,則擴大搜索范圍為5×5、7×7、……,直到圖像邊界,連接最近的邊界點,生成枝切線。
6)重復以上過程,直至枝切線連接并平衡所有殘差點為止。值得注意的是,在殘差點連接過程中,已達到“電荷”平衡的殘差點不用重復連接,可減少枝切線數量。
步驟2:利用二階差分函數計算干涉圖中各像元的權值[16],并設置閾值將該權值二值化為“0”和“1”,0表示低質量像元,1表示高質量像元。
(1)
(2)
式中:W為纏繞算子,將相位差值限制于(-π,π]之間;φi,j表示干涉圖(i,j)像元纏繞相位;Ri,j表示(i,j)像元的可信度。二值化質量權值表示如式(3)所示。
(3)
式中:θ表示可信度閾值;qi,j表示(i,j)像元二值化質量權值。
步驟3:由枝切線分布圖以及干涉圖二值化質量權值矩陣生成可靠性掩模圖,非枝切線上的二值化質量權值為“1”的像元為權值為“1”的高質量像元,枝切線上的像元、二值化質量權值為“0”的像元以及枝切線圍成的閉環區域中的像元等像元為權值為“0”的低質量像元。
步驟4:把可靠性掩模圖、像元擴散策略以及UKF相位解纏程序結合起來,按以下步驟對干涉圖進行解纏。
1)首先選取非枝切線上可信度最高的像元為起點,其纏繞相位作為解纏相位估計值,其估計誤差方差預設為0.6;利用UKF相位解纏程序解纏4個鄰接像元中的高質量像元,隨后將已解纏像元鄰接像元中的未解纏高質量像元依次存儲在“鄰接列”中。
2)按順序取出“鄰接列”中的待解纏像元,利用UKF相位解纏程序對該像元進行解纏,隨后將解纏像元鄰接像元中的未解纏高質量像元依次存儲在“鄰接列”中。
3)若“鄰接列”不為空,回到步驟2),若為空,則轉步驟4)。
4)根據已解纏像元信息,利用UKF相位解纏程序按行(或列)的方式逐一解纏余下的未解纏像元。
圖1進一步給出了本文算法解纏步驟示意圖。圖1(a)中“黃色”像元(非枝切線上可信度最高的像元)為起始像元,“綠色”像元為高質量像元,“紅色”像元為低質量像元,先利用UKF相位解纏程序進行解纏“綠色”像元;“紅色”像元將在步驟4)中進行解纏。圖1(b)給出了執行步驟1)之后的結果,其中,“黑色”像元為已解纏像元,“藍色”像元為存儲在“鄰接列”的待解纏像元,其余像元為未處理的像元。需要注意的是,若干涉圖條紋十分密集,無法利用UKF相位解纏程序從選擇的起點進行解纏時,可以選擇可信度次高的像元為起點,按照圖1(b)方式進行解纏。

圖1 像元擴散法解纏示意圖
為了驗證本文算法性能,利用不同算法包括BUT、非加權迭代最小二乘法(unweighted iterative least squares,ILS)、QGPU、UKFPU,以及本文算法,在同一MATLAB軟件環境(Intel i7-6700U@3.40 GB CPU+8 GB RAM)下解纏不同的模擬干涉圖和實測干涉圖,并對各算法解纏結果進行比較分析。
圖2為兩幅不同模擬干涉圖(256像素×256像素)。圖3(a)~圖3(b)分別為圖2(a)~圖2(b)含噪聲的纏繞相位圖,其信噪比依次為3.01 dB、2.18 dB。依次用BUT、ILS、QGPU、UKFPU以及本文算法對上述兩幅干涉圖進行解纏,解纏結果分別如圖4~圖8、圖9~圖13所示。

圖2 模擬干涉圖

圖3 含噪纏繞相位圖

圖4 BUT算法解纏相關處理結果

圖5 ILS算法解纏結果

圖6 QGPU算法解纏結果

圖7 UKFPU算法解纏結果

圖8 本文算法解纏結果

圖9 BUT算法解纏相關處理結果

圖10 ILS算法解纏結果

圖11 QGPU算法解纏結果

圖12 UKFPU算法解纏結果

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

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

表2 各算法解纏時間 s
圖14為經中值濾波處理過后截取的部分Enta火山干涉圖,各算法解纏結果如圖15所示,其解纏相位重纏繞結果如圖16所示。

圖14 局部Enta火山干涉圖

圖15 各算法解纏結果

圖16 各算法解纏相位重纏繞結果
圖17為三峽實測干涉圖,各算法解纏結果如圖18所示,解纏相位重纏繞結果如圖19所示。
由圖15(a)以及圖18(a)可以看出,BUT算法解纏相位圖存在部分無法解纏以及明顯不一致的區域。由圖15(b)以及圖18(b)可見,ILS算法解纏結果較為光滑,但根據圖16(a)以及圖19(a)可以發現,ILS算法重纏繞結果相較于圖14以及圖17原始干涉圖條紋差別較大,大量干涉條紋細節信息丟

圖17 三峽實測干涉圖

圖18 各算法解纏結果

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

表3 各算法解纏時間 s
本文提出一種高效UKF相位解纏算法,該算法是可靠性掩模圖、像元擴散策略以及無味卡爾曼濾波算法相結合的結果。模擬數據實驗以及實測數據實驗結果證明了高效UKF相位解纏算法有效性,相對于傳統算法如BUT、ILS以及QGPU,具有更穩健的相位解纏能力,且在時間耗費較少的同時其解纏精度與UKFPU算法相當。然而,所提出的方法處理復雜地形仍有改進空間,且相位解纏時間消耗還有優化的空間。努力尋找高效穩健的相位解纏技術是我們未來努力的目標。