采用JVC算法生成枝切線的干涉相位解纏方法

王朝霞 劉永信 張 杰 范陳清

①(內蒙古大學計算機學院 呼和浩特 010021)

②(內蒙古大學電子信息工程學院 呼和浩特 010021)

③(自然資源部第一海洋研究所 青島 266061)

1 引言

3維成像雷達高度計[1，2]綜合了傳統高度計高精度的測高能力和合成孔徑雷達高分辨率的成像能力，使海面高程的面測量成為可能，是遙感技術的最新發展前沿。2016年9月15日，搭載世界首例3維成像雷達高度計的“天宮二號”空間實驗室成功發射，在國際上首次實驗驗證了采用小入射角和短基線干涉測量技術實現寬刈幅海平面高度厘米級精度測量的工作機理。因3維成像高度計測高采用了與干涉合成孔徑雷達(Interferometric Synthetic Aperture Radar， InSAR)基本相同的干涉測量技術[1]，故而在高程反演階段必須通過相位解纏這一環節將濾波后位于主值區間的纏繞相位轉換為絕對相位。相位解纏是干涉測量技術的關鍵環節，解纏的結果將直接影響最終的高程測量精度。在過去的數十年中，許多學者致力于相位解纏算法的研究和改進，但如何在強噪聲干擾環境中進一步提升干涉相位解纏的效率和速度，目前仍是一個富有挑戰的研究課題。

現有的微波干涉測量相位解纏方法主要有兩類，第1類是與展開路徑有關的局部法，典型的有枝切法、質量圖導引法[3]、最小斷點法等；第2類是與展開路徑無關的全局法，典型的有最小二乘法[4]、最小費用流法[5]、最小零范數法等。此外，也有綜合運用以上兩類方法或提出了基于混合數學模型的方法[6]，以在解決實際問題中尋求更好的實驗結果。其中，Goldstein等人[7]于1988年在研究In-SAR時提出的枝切法是一種非常經典的算法。該方法通過識別殘差點，并根據殘差原理設置正確的枝切線，從而選擇合適的積分路徑，實現相位解纏。但殘差原理只是不連續相位存在的充分非必要條件，而且枝切線的放置因方法不同往往結果也不相同，不合理的枝切線放置可能導致相位跳躍和無法解纏的“孤島問題”[8]。針對枝切法存在的缺點和問題，近年來國內外學者在理論研究和應用中進行了改進。蔣銳等人[9]提出了等效殘差點的概念，并基于等效殘差點設置枝切線，解決了積分路徑穿過殘差點密集區所引起的展開相位跳變問題。張妍等人[10]通過虛擬組合殘差點改進枝切線算法，在一定程度上解決了殘差點較密集區域在放置枝切線時易產生“孤島”的問題。De Souza等人[11]在獲取物體3維輪廓時，利用跳變相位解決了剩余殘差點的相位平衡問題，進而通過尋求最小化處理時間實現了動態全息技術中的干涉相位解纏。

本文在研究3維成像高度計高程反演時，提出了一種改進的枝切線干涉相位解纏方法，利用JVC(Jonker-Volgenant-Castanon)全局最優分配算法在正負殘差點之間放置枝切線，以進一步縮短枝切線總長度，避免“孤島”產生，提升干涉相位圖中像素的解纏率，降低高程反演誤差。

2 相位解纏與Goldstein枝切線算法

由于復數相位的周期性，對配準后兩幅相干復圖像作復共軛相乘得到的干涉相位圖中，相位值是被周期折疊后的相位主值，位于(-π，π]之間，與真實相位值之間相差2 π的整數倍。

其中，Φ(i，j)為真實相位值，Ψ(i，j)為纏繞相位值，k(i，j)為整數。相位解纏的目標就是恢復被模糊掉的周期分量2k(i，j)π。

一般情況下，干涉相位圖的采樣率均滿足Nyquist原理，因此干涉相位圖中相鄰像素點之差不超過半個周期，也就是差的絕對值不大于 π。在理想的無噪聲環境中，相位解纏只需根據干涉相位圖中像素值在距離向和方位向的偏導數進行簡單的積分即可。

然而在實際環境中，由于受噪聲等諸多因素的影響，干涉相位中存在不連續及不同方向上積分結果不一致的殘差點[12](相鄰4個像素點組成的正方形回路的相位差積分不為0，而是-2π(負殘差)或2π(正殘差))。根據相位解纏的基本原則，采用任何一種包含單個或多個正負殘差點數不等的積分路徑進行解纏，都會產生不一致的結果。

為避免相位解纏后出現不一致的現象，Goldstein枝切線算法首先查找并確定干涉相位圖中所有的正負殘差點，然后在數量相等的相鄰正負殘差點之間放置“平衡”枝切線，如果正負殘差點數量不等，則把剩余的殘差點連接到圖像最近的邊界上，這樣也認為是“平衡”的。進而控制積分路徑，使其繞過這些枝切線，采用式(3)進行相位解纏，即可得到一致解纏結果。因所有積分路徑包含的正負殘差點值總和都為0，所以不會形成全局誤差。但同時在枝切線另一側的像素點將會有相位不連續情況出現。為使不連續像素點數量最小，Goldstein算法在生成枝切線時要求滿足總長度最短的原則。此外，在生成枝切線時，在殘差點相對密集部分有時形成多條枝切線環繞的像素區域(即“孤島現象”)，造成被環繞區域像素無法解纏的情況，這些都是改進枝切線算法時需要重點考慮的問題。

3 基于JVC算法生成枝切線的解纏方法

JVC算法[13]是由Jonker， Volgenant和Castanon 3人聯合提出的一種全局最優線性分配算法，常用于兩個數據點集之間的一對一分配。JVC算法將線性分配問題描述為：在約束條件

成立的前提下，求最小值

其中，c[i，j]為代價，x[i，j]為權值。

(2) 采用上述基于JVC算法生成枝切線的方法，將殘差點分為兩類，分別在殘差點與其最近邊界點之間或正負殘差點之間放置枝切線；最后遍歷所有的正負殘差點，若還有未連接枝切線的，則在其與最近邊界點之間放置枝切線，以此確保枝切線的極性平衡。

(3) 在枝切線相對較為稀疏的區域，選取一個非枝切線上的像素點作為起點，繞過所有枝切線采用式(3)的方法進行積分，即可得到解纏相位。

(4) 對于枝切線上未能解纏的像素點，根據解纏后相位圖應具有連續性的原理，將枝切線上點的像素值用7×7鄰域內已解纏像素點的均值替換，完成全部像素點的相位解纏。

4 方法驗證

4.1 3維成像高度計海面仿真干涉相位圖驗證

為驗證算法的可行性和有效性，本文根據“一發雙收”天線3維成像高度計的成像原理[1]仿真了Ku波段的海面干涉相位圖。首先采用P-M海浪譜和雙尺度模型模擬海洋表面，求得海面上每一點的高程數據。然后利用Delaunay三角剖分法將海面劃分為三角形小面元，采用物理光學(Physical Optics，PO)模型及其Kirchhoff近似解計算仿真海面區域的后向散射系數。最后根據3維成像雷達高度計的成像原理使用后向投影(Back Projection， BP)算法進行仿真成像，得到兩幅仿真3維成像高度計海面相干復圖像。經圖像配準、去平地效應及濾波后得到如圖1(a)所示的干涉相位。

在3維成像高度計海面仿真干涉相位圖中提取殘差點，共得到正殘差點578個，負殘差點577個。分別采用本文算法和Goldstein枝切線算法生成枝切線，結果如圖1(b)和圖1(c)所示。Goldstein枝切線算法生成的枝切線總長度為3313像素，本文算法生成的枝切線總長度為1723像素。

為驗證本文算法的可行性和有效性，分別采用Goldstein枝切線算法、四向加權最小二乘法[14]、快速傅里葉變換(Fast Fourier Transform， FFT)算法[15]和本文算法對3維成像高度計海面仿真干涉相位圖進行解纏，結果如圖2所示。Goldstein枝切線算法的像素解纏率為76.99%，本文算法的像素解纏率為88.03%，采用7×7鄰域內的解纏像素均值替換枝切線上點的像素值后，解纏率可達100%。觀察解纏后的4幅干涉相位圖，可見四向加權最小二乘法解纏后相鄰相素之間的連續性在4種算法之中最好，階躍式跳變最少，本文算法與FFT算法次之，傳統Goldstein枝切線算法相對最差。

為定量比較4種算法的解纏結果，利用式(9)分別將解纏相位轉換為圖3(a)—圖3(d)所示的相對數字高程。

其中，h 為相對數字高程， H為3維成像高度計雷達平臺高度(本文取393 km)， Rm為雷達主天線到海面高程點的斜距， λ為電磁波波長，φ 為解纏后的相位值， B為基線長度(本文中取10 m)。

統計圖3(a)—圖3(d)所示數字高程與圖3(e)所示的原始模擬數字高程的均方根誤差(Root Mean Square Error， RMSE)，并在同樣的軟硬件環境(Intel Core i5-4590 4核，16 GB DDR3 1666 MHz內存，Win7 x64旗艦版操作系統)下分析4種算法的時間復雜度及實測運行時間，詳細對比結果如表1所示。表中，分析時間復雜度時， M， N分別為干涉相位圖的行、列數，m ， n分別為正、負殘差點的個數，s=max(m，n)， W為四向加權最小二乘法控制加權窗口的尺寸， I為FFT算法的迭代次數。

圖1 3維成像高度計海面仿真干涉相位及其枝切線示意圖

圖2 3維成像高度計海面仿真干涉相位解纏結果

由圖3及表1可知，雖然四向加權最小二乘法解纏后相位值的連續性和算法運行時間優于其他3種算法，但解纏相位對應的數字高程誤差大于本文算法和FFT算法，這是因為四向加權最小二乘法在進行相位值擬合時導致了誤差的全局擴散。采用4次迭代FFT算法時，雖然運行時間優于本文算法，但轉換后的高程誤差約為本文算法的2倍。本文算法的執行時間優于Goldstein枝切線算法、略次于FFT算法，但解纏精度均高于其他3種算法，能夠滿足3維成像高度計厘米級測高誤差的要求。

就時間復雜度而言，本文算法主要取決于殘差點的稠密程度，越稀疏時間復雜度越低。當m ＜M， n ＜N ， 1+log2s ＜W2時，本文算法的時間復雜度將會低于其他3種算法；當nlog2s ＞MN時，本文算法的時間復雜度將高于其他3種算法。介于上述兩種情況中間時，本文算法的時間復雜度一般低于Goldstein枝切線法，高于四向加權最小二乘法，接近FFT算法。

圖3 海面仿真干涉相位圖解纏后對應的相對數字高程

表1 3維成像高度計海面仿真干涉相位圖4種算法解纏結果對比

4.2 實際Etna火山地區干涉相位圖驗證

為說明算法的通用性和有效性，利用本文算法和Goldstein枝切線算法、四向加權最小二乘法、FFT算法分別對ERS1/2兩次經過Etna火山地區時獲取的相干SLC(Single Look Complex)數據進行解纏。本文從經過多視處理的兩幅復圖像中分別選取500像素×500像素的區域進行實驗，經配準、去平地效應、濾波后的干涉相位圖如圖4(a)所示。提取殘差點，得到正殘差點3383個、負殘差點3389個。本文算法和Goldstein枝切線算法生成的枝切線如圖4(b)、圖4(c)所示，本文算法和Goldstein枝切線算法生成的枝切線長度分別為13546像素和15664像素。解纏后的結果如圖5所示。將解纏相位轉換為圖6(a)—圖6(d)所示的相對數字高程，并與圖6(e)所示對應經緯度的ASTER GDEMV2 30 m分辨率數字高程數據[16]進行比較，計算RMSE，4種算法的詳細對比結果如表2所示(統計算法執行時間的硬件環境與4.1節相同)。

圖4 Etna火山地區干涉相位及枝切線示意圖

圖5 Etna火山地區干涉相位圖解纏結果

由圖5、圖6及表2可知，雖然FFT算法和四向加權最小二乘法兩種全局性算法解纏后相位的連續性和執行時間優于本文算法，但解纏后相位所代表的高程走勢與原干涉相位圖不一致，且與GDEMV2之間的RMSE相對較大，這是由于兩種算法均引起了誤差的全局性擴散。傳統Goldstein枝切線算法解纏雖未改變高程走勢，但未解纏相位和相鄰像素值之間的階躍式跳變較多，解纏結果連續性較差，且與GDEMV2的RMSE較大。本文算法兼顧了解纏的精確性和相位連續性，消除了大部分相鄰相素之間階躍式跳變，避免了誤差的全局擴散。

分別對比圖2(b)、圖2(c)與圖4(b)、圖4(c)中兩種算法生成的枝切線，可以看出采用本文算法生成枝切線可以有效防止其相互連接形成“孤島”，使得干涉相位圖中除枝切線以外的像素點都可以通過積分進行解纏。這是因為，不論是與邊界點直接相連生成枝切線，還是通過JVC線性分配后在正負殘差點對之間放置枝切線，每個殘差點只存在于一條枝切線上，任何一條枝切線都是兩點之間的連線(忽略像素數字化時取整的因素后可近似為直線段)，是無法連成閉環的。

圖6 Etna火山地區干涉相位解纏后對應的相對數字高程

表2 Etna火山地區干涉相位圖4種算法解纏結果對比

5 結束語

本文為縮短枝切線相位解纏算法中生成枝切線的總長度，提升干涉相位圖的解纏率和準確性，提出了采用JVC全局最優線性分配算法在殘差點之間放置枝切線，以此對Goldstein枝切線算法進行改進。通過對3維成像高度計海面仿真干涉相位圖和Etna火山地區干涉相位圖進行解纏實驗，結果表明在同樣的軟硬件運行環境中，本文算法在生成枝切線總長度、像素解纏率和運行時間上均優于Goldstein枝切線法，解纏精度優于Goldstein枝切線法、四向加權最小二乘法和FFT算法，同時能夠較為有效地防止“孤島現象”產生，為干涉相位解纏提供了一種可行的方法。然而，在殘差點數量少且正負殘差點距離較近的干涉相位圖中，對Goldstein枝切線算法的改進并不明顯；在信噪比低、殘差點密集的干涉相位圖(或局部區域)中，算法運行時間相對其他算法較長，且無法完全避免“孤島現象”的產生，這是下一步需要重點考慮的優化方向。