本文内容來源于《測繪學報》2024年第4期（審圖号GS京(2024)0714号）
InSAR相位解纏最大流/最小割權值改進算法
高延東
,1,2, 賈義琨,1,2, 李世金1,2, 陳宇1,2, 李懷展1,2, 鄭南山1,2, 張書畢1,2
1.中國礦業大學自然資源部國土環境與災害監測重點實驗室，江蘇 徐州 221116
2.中國礦業大學環境與測繪學院，江蘇 徐州 221116
摘要:InSAR已經在高精度DEM反演中得到了廣泛應用，相位解纏技術是影響DEM擷取精度的關鍵步驟之一，然而大梯度變化區域一直是影響解纏結果精度的核心問題。針對該問題，本文提出一種基于勢函數改進權值的最大流/最小割相位解纏算法。首先，針對PUMA模型權值設定不合理的問題進行研究，利用外部已有DEM擷取相位梯度變化的先驗資訊，将視窗最大絕對相位梯度值代入相應的勢函數公式中擷取權值；然後，通過對勢函數權值的設定進行門檻值調節，解決由于勢函數權值設定不合理導緻PUMA勢函數無法發揮作用引起的解纏錯誤問題，進而提高大梯度變化區域相位解纏精度；最後，通過仿真資料和真實TanDEM-X InSAR資料對本文算法進行驗證，并與已有方法進行對比。結果表明，仿真資料中本文算法可以提高至少44.93%的解纏精度，真實資料中在大梯度變化區域，本文算法較已有算法可以獲得更大範圍的有效解纏結果。關鍵詞： InSAR ; 勢函數 ; 最大流/最小割 ; 相位解纏
基金資助:國家自然科學基金(42001409)
作者簡介:高延東（1988—），男，博士，副教授，研究方向為單/多基線相位解纏、InSAR高精度DEM反演。E-mail：[email protected]
本文引用格式
高延東, 賈義琨, 李世金, 陳宇, 李懷展, 鄭南山, 張書畢. InSAR相位解纏最大流/最小割權值改進算法. 測繪學報[J], 2024, 53(4): 644-652 doi:10.11947/j.AGCS.2024.20220633
GAO Yandong. The improved max-flow/min-cut weight algorithm for InSAR phase unwrapping. Acta Geodaetica et Cartographica Sinica[J], 2024, 53(4): 644-652 doi:10.11947/j.AGCS.2024.20220633
閱讀全文
InSAR已經在DEM擷取、礦山安全監測及滑坡災害預警等領域中得到了廣泛應用[1-3]。相位解纏是InSAR資料處理的關鍵步驟之一，其結果精度将直接影響最終擷取産品的精度[4]。相位解纏旨在從纏繞在(－π，π]之間的相位中擷取真實相位[5-6]。然而大梯度變化區域相位解纏精度并不理想，這是影響InSAR擷取最終産品精度主要因素之一[7]。
相位解纏一直是InSAR資料處理領域的技術難點及研究熱點，主要分為兩類[8]。一類是以枝切法[9]、品質圖法[10]和最小不連續法[11]為代表的路徑跟蹤法，該類方法通過計算殘差點、設定枝切線、優化解纏積分路徑等手段，可以快速、精準地擷取解纏結果。然而，該類方法在低品質區域會産生誤差傳遞，甚至解纏空洞現象，嚴重影響解纏結果精度。另一類是以Lp範數法[12]、最小費用流(minimum-cost flow，MCF)方法[13]及統計費用流(statistical cost network flow，SNAPHU)[14]為代表的梯度估計優化法，該類方法旨在最小化相位梯度估計值與真實相位梯度之間的差異，進而擷取最終解纏結果，該類方法具有較好的穩定性。然而，該類方法容易産生全局性的誤差傳遞，影響最終的解纏結果[15]。此外，相位連續性假設是解纏方法常用的準則，該假設認為無噪聲的情況下，真實相位梯度等于纏繞相位梯度[16]，即相位梯度被限制在(－π，π]之間，然而，在真實情況中尤其是在大梯度變化區域，該假設是不适用的，是以大多數解纏方法在大梯度變化區域無法獲得理想結果[17]。針對以上問題，多基線相位解纏方法被提出。其中具有代表性的多基線相位解纏方法主要包括中國餘數定理[18]及兩階段規劃方法[19]，該類方法通過聯合不同基線長度的幹涉相位，可以擺脫相位連續性假設的限制，在大梯度變化區域也可以擷取較為理想的解纏結果。然而多基線相位解纏方法存在噪聲穩健性差的問題，不僅如此，現階段不是所有區域都有多基線方法所需的幹涉相位資料，這也是限制此類方法應用的主要因素之一。是以，如何通過單幅幹涉相位擷取大梯度變化區域更高精度的解纏結果就顯得尤為重要。
近年來，針對單幅幹涉相位大梯度變化區域相位解纏研究得到了廣泛關注。文獻[20]提出了一種基于礦區沉降預測模型輔助的相位解纏政策，該方法采用礦區沉降預測模型，結合礦區沉陷參數對礦區形變進行預測，然後進行疊代處了解決大梯度變化區域的解纏問題。文獻[21]提出了一種結合先驗資訊輔助的相位解纏方法，該方法可以通過結合先驗資訊解決大梯度變化區域的相位解纏問題。文獻[22]提出了一種基于外部DEM輔助的切片相位解纏算法，該方法通過外部DEM對單基線可解相位進行仿真，以此對待解纏相位進行切片處理，進而降低相位梯度變化，解決大梯度變化區域相位解纏問題。基于馬爾可夫随機場(Markov random field，MRF)，文獻[23]首次提出了最大流/最小割相位解纏算法(phase unwrapping max-flow/min-cut algorithm，PUMA)，并已證明凸勢函數對于處理連續相位具有很好的效果，非凸勢函數更适用于處理大梯度變化相位，是以PUMA針對大梯度變化區域解纏具有明顯的優勢。文獻[24]提出了線性過程MRF(line process-MRF，LP-MRF)方法，結合非凸勢函數，采用8鄰域線性視窗檢測相位連續性的方式設定權值，提高了PUMA在大梯度變化區域的相位解纏能力。文獻[25]提出了一種改進的PUMA(improved-PUMA，IM-PUMA)方法，采用拟凸勢函數自适應識别大梯度相位，并借助外部DEM擷取最大梯度先驗資訊輔助設定權值，提高了大梯度變化區域解纏結果的精度。IM-PUMA雖然可以獲得較為理想的結果，然而，其模型權值的設定在大梯度變化區域存在兩極分化現象，在地形平坦區域權值設定沒有達到理想的較大值，且在極端大梯度變化的情況下，會因為權值設定過小導緻一定的解纏誤差出現，是以IM-PUMA模型的權設定仍需進一步改進。
本文基于IM-PUMA方法，提出一種基于勢函數改進權值的PUMA。該方法首先利用外部輔助DEM擷取先驗相位梯度資訊，将視窗最大相位梯度代入相應的勢函數公式，并進行歸一化處理，進而擷取權值；采用對權值設定門檻值的方式，解決在極端大梯度變化區域由于權值設定過小勢函數無法發揮作用引起的解纏錯誤問題，提高整體相位解纏精度，進而擷取高精度DEM産品。經過仿真資料和真實TanDEM-X InSAR資料試驗，結果證明本文算法大梯度變化區域可以獲得較已有算法更高精度的結果。
1 相位解纏理論
1.1 相位解纏的基本原理
InSAR技術實際測量擷取的相位資訊是真實相位的主值，即纏繞相位，其值在(－π，π)之間。纏繞相位與真實相位之間相差2π的整數倍，真實相位與纏繞相位的關系為
(1)
式中，s為像素點；ψ(s)為真實相位；φ(s)為纏繞相位；k(s)為梯度模糊系數。對于單基線相位解纏，相位梯度被相位連續性假設限制在(－π，π)内，相位連續性假設可表示為
(2)
式中，Δψ(s)為像素s與像素s－1之間的真實相位梯度。該假設是現階段單基線相位解纏算法的必要前提，然而大梯度變化區域真實相位梯度變化往往大于纏繞相位梯度變化，這也是相位解纏誤差産生的主要因素之一。
1.2 IM-PUMA方法
PUMA是一種基于貝葉斯理論的相位解纏方法，PUMA解纏結果精度依賴于勢函數和權值設定兩個先驗資訊，傳統PUMA通過相幹系數設定勢函數所需的權值，然而相幹系數不足以表示相鄰像素之間的相位連續性，是以會對勢函數造成一定的幹擾，導緻解纏誤差增大。故此，針對傳統PUMA的模型缺陷，文獻[25]提出了IM-PUMA方法，該方法采用了更為合理的拟凸勢函數自适應識别大梯度相位，并利用外部DEM的梯度資訊來輔助設定權值，大大增加了勢函數權值的可靠性。IM-PUMA方法以PUMA為基礎，将幹涉圖抽象為MRF，其中MRF的節點對應幹涉圖像素，節點之間的邊對應幹涉圖中每對像素之間的纏繞相位梯度Δφ。以MRF為基礎，建構IM-PUMA方法的能量函數模型E(k|ψ)為
(3)
式中，p是幹涉圖像素點；P是幹涉圖像素點集合；Np是像素p的鄰域系統；Δψpq是像素p與鄰像素q之間的真實相位梯度；Δkpq是像素p與鄰像素q之間的梯度模糊系數內插補點；kp和kq分别是像素p與鄰像素q的梯度模糊系數；wpq代表權值，是表示像素p與鄰像素q之間相位連續性的名額。根據外部輔助DEM最大梯度先驗資訊計算獲得，像素p和鄰像素q組成的像素對稱作團Cpq，V(·)是勢函數
(4)
式中，a是調節參數。根據IM-PUMA能量函數模型，定義第t次疊代的能量函數值為
(5)
定義δ∈{0，1}，最小化第t次疊代的能量函數，則第t＋1次疊代的能量函數值表示為
(6)
式中，
和
分别表示像素p在第t和t＋1次疊代的相位模糊系數估計值。當
，
和
分别表示鄰像素q在第t和t＋1次疊代的相位模糊系數估計值；當
，在初始時刻
，
和
分别表示第t和t＋1次疊代的相位梯度估計值。由式(5)、式(6)可知，能量最小化問題可以映射為一個二進制優化問題，式(6)中團Cpq的能量函數值
可以表示成
、
、
、
這4種能量表達式，将這4種能量表達式根據文獻[26]的方法映射至一個賦有權值的“有向基本圖”中，将所有團的有向基本圖組合得到一個賦有權值的“有向圖”，對完整的有向圖進行最大流/最小割操作即可得到第t＋1次疊代的能量函數值[27]，直至能量函數值不再減少，即可得到最終的解纏相位。雖然IM-PUMA在大梯度變化區域較傳統PUMA可以獲得較為理想的結果，然而IM-PUMA的勢函數權值僅使用梯度值獲得，在大梯度變化區域權值設定上存在兩極分化問題，尤其是在較為極端的大梯度變化區域，勢函數在能量最小化過程中會産生無法正常發揮作用的現象，導緻IM-PUMA在大梯度變化區域仍舊存在一些缺陷。
2 基于勢函數改進權值的PUMA算法
針對上述問題，本文基于IM-PUMA算法提出了一種基于勢函數改進權值的PUMA算法。采用類似IM-PUMA算法中最大梯度的定義，令C表示最大權值，G定義為視窗最大絕對相位梯度，G(i，j)在以像素(i，j)為中心大小為2n＋1的滑動視窗内計算得到，以3×3視窗為例，如圖1所示。
圖1
圖1 3×3視窗示例
Fig. 1 3×3 window example
視窗内水準、垂直、45°方向和135°方向上的視窗絕對相位梯度分别記作Gh、Gv、Gd1和Gd2。4個方向的相位梯度可表示為
(7)
根據式(7)可以得到像素(i，j)處的視窗最大絕對相位梯度為
(8)
值得注意的是，本文提出的視窗最大絕對相位梯度與IM-PUMA中的最大梯度的差別在于，IM-PUMA中的視窗最大梯度可以是負值，而本文提出的視窗最大絕對相位梯度是根據4個方向梯度值的絕對值得到，即存在限制
(9)
由于式(9)限制的存在，可以完全避免IM-PUMA将權值歸一化之後在大梯度變化區域權值設定出現的兩極分化現象，無論是在大梯度上升或下降的區域，權值始終較小，符合大梯度變化處勢函數權值應較小的配置設定原則。
在IM-PUMA中，僅将所有像素的最大梯度(存在負值)進行歸一化得到權值，這與IM-PUMA能量函數模型中的拟凸勢函數模型并不能達到理想的結合。同樣，若僅将視窗最大絕對相位梯度G(i，j)進行歸一化得到權值用于能量最小化，也隻能夠消除大梯度變化區域兩極分化現象造成的影響，粗略地将相位連續性資訊表示出來。是以，為了能夠更好地利用相位梯度将相位連續性資訊合理地表示，本文提出将視窗最大絕對相位梯度G(i，j)結合拟凸勢函數公式進一步優化權值，像素(i，j)處的最大權值可以表示為
(10)
對所有像素的最大權值C作歸一化處理，使C的值在區間[0，1]内
(11)
式中，Cmin為C中的最小值；Cmax為C中的最大值。此時歸一化的最大權值C′并不是最終的勢函數權值。C′(i，j)接近1表示像素(i，j)處存在大梯度，應賦予對應勢函數較小的權值；C′(i，j)接近0表示像素(i，j)處地形平緩，應賦予對應勢函數較大的權值。是以像素(i，j)處的能量函數模型的勢函數權值w(i，j)由式(12)給出
(12)
幹涉圖能量函數模型的勢函數權值表達式為
(13)
值得注意的是，式(10)中對視窗最大絕對相位梯度G(i，j)賦予了平方項，這一操作使得原本較小的C′(i，j)變得更小，結合式(12)可知，w(i，j)的值會有所增加，即對視窗最大絕對梯度G(i，j)賦予平方項的操作增大了地形平坦區域的權值，使得地形平坦區域的權值有接近1的可能，符合地形平坦區域的相位應當具有較高可信度的原則。此外采用拟凸勢函數本身的公式計算得到的權值，與拟凸勢函數模型可以更加合理地比對，進一步提高了大梯度區域相位連續性資訊的可信度。
在較為極端大梯度變化區域，G(i，j)相較于其他區域會顯得非常大，即這種情況下的權值w(i，j)會非常接近0，這會導緻勢函數模型無法正常發揮作用。這是因為該權值認為此處是幾乎垂直的斷層，應當極力抹去勢函數在斷層處的作用。然而實際問題中較為極端的大梯度變化區域也存在一定梯度，是以，在這種情況下權值w(i，j)幾乎接近0顯然是不合理的。針對上述問題，本文提出一種通過對權值設定門檻值的方式來解決權值過小問題。
假設權值門檻值為M，通過門檻值對求得的勢函數權值w加以限制，限制後勢函數權值表達式W為
(14)
若w(i，j)＝M，則稱像素(i，j)處的視窗最大絕對相位梯度為臨界梯度GM，臨界梯度可根據外部輔助DEM資料來确定，根據臨界梯度GM和式(12)可以計算得到權值門檻值。臨界梯度GM的具體擷取方法為：分别計算外部DEM方位向和距離向的相位梯度，統計所有相位梯度分布規律，根據相位梯度分布情況确定GM，一般GM對應于所有相位梯度值的前10%~30%區間内，超過GM值的相位梯度均可認為是幹涉圖中較為極端的大梯度變化，需要對最大絕對梯度值超過GM的大梯度區域的權值進行補償，進而解決權值過小導緻勢函數無法發揮作用的問題。權值門檻值M取值範圍在[0.1，0.3]區間内。
3 試驗結果與分析
本文使用仿真資料和真實TanDEM-X InSAR資料對本文算法的有效性進行驗證分析，并與MCF、SNAPHU及IM-PUMA方法進行對比。本文采用均方根誤差作為相位解纏精度的定量評估名額
(15)
式中，N表示總的像素數；
是解纏相位；
是參考解纏相位。
3.1 仿真資料試驗
圖2為本文采用的仿真資料，資料大小為458×157。圖2（a）為仿真輔助解纏相位；圖2（b）為對應的纏繞幹涉相位；圖2（c）為IM-PUMA算法根據仿真解纏相位計算得到的權值；圖2（d）為本文算法根據仿真解纏相位計算得到的權值，本文仿真資料的權值門檻值為0.3。由圖2（a）和圖2（b）可以看出，仿真資料存在明顯的大梯度變化區域，這為驗證不同相位解纏算法在大梯度變化區域的性能創造了良好條件。為了可以更加精細地識别大梯度區域，本文将計算仿真資料權值的視窗大小設定為3×3。
圖2
圖2 仿真資料試驗資料集
Fig. 2 Simulation data experimental dataset
如圖3所示，圖3（a）—（d）分别為MCF、SNAPHU、IM-PUMA及本文算法的解纏結果，由圖可以看出大梯度變化區域4種方法均産生不同程度的解纏誤差；圖3（e）—（h）分别為不同方法解纏結果與圖2（a）的誤差，可以看出，MCF及SNAPHU産生了明顯的解纏誤差現象，IM-PUMA由于權值設定不合理的原因，導緻其在較極端大梯度變化區域也出現了明顯的解纏斑塊。由圖3（h）可以看出，本文算法雖然也産生了部分解纏誤差，但是已經消除大部分的解纏斑塊。本文用RMSE對不同方法的解纏誤差進行定量評估，見表1，由于較極端大梯度變化的影響，MCF、SNAPHU和IM-PUMA解纏結果的RMSE結果均不理想，而本文算法可以較好地解決大梯度變化區域的解纏問題，是以其RMSE是幾種解纏方法中精度最高的，這也驗證了本文算法在較極端的大梯度變化區域可以獲得較已有方法更好的解纏結果。
圖3
圖3 仿真資料試驗不同相位解纏方法的解纏結果和誤差
Fig. 3 Unwrapping results and error plots for different phase unwrapping methods with simulated data
表1 模拟資料的評估結果Tab. 1 Evaluation results of simulated data
| 相位解纏方法 | RMSE |
|---|---|
| MCF | 1.633 4 |
| SNAPHU | 1.385 3 |
| IM-PUMA | 1.999 5 |
| 本文算法 | 0.762 9 |
新視窗打開| 下載下傳CSV
3.2 真實資料試驗
為了進一步驗證本文算法的有效性，選取渭南地區真實TanDEM-X資料進行處理，如圖4所示。圖4（a）為真實資料的Goolgle圖像，圖4（b）為由外部輔助DEM模拟的解纏相位，圖4（c）為TanDEM-X真實幹涉相位，圖4（d）為IM-PUMA的權值，圖4（e）為本文算法的權值，本文真實資料的權值門檻值為0.3。由圖4（a）可以看出，研究區域存在大面積的大梯度變化區域，這會嚴重影響相位解纏精度。根據真實資料研究區域的情況，本文将計算真實資料權值的視窗大小設定為3×3。
圖4
圖4 真實資料試驗資料集
Fig. 4 Real data experimental dataset
圖5為不同相位解纏方法的解纏結果及誤差圖，圖5（a）和（e）分别為MCF方法的解纏結果和相應的解纏誤差，MCF方法的解纏結果存在明顯的解纏跳變誤差，這也說明了在大梯度變化區域MCF方法極易産生較大的解纏誤差。圖5（b）和（f）分别為SNAPHU方法的解纏結果和誤差，雖然SNAPHU方法可以獲得較MCF方法更好的解纏結果，但是仍然存在較大面積的解纏跳變誤差。圖5（c）和（g）分别為IM-PUMA的解纏結果和誤差，雖然仍存在一定面積的解纏跳變誤差，但相比于前兩種方法，IM-PUMA的解纏結果有明顯的提升。然而由于權值設定的問題，仍然産生了部分的解纏誤差。圖5（d）和（h）分别為本文算法的解纏結果和誤差，相比于前3種方法，本文算法取得了最為理想的解纏結果，存在最少面積的解纏跳變誤差。為進一步論證提出解纏方法的性能，本文利用RMSE對不同方法進行對比分析(表2)，由于大梯度變化區域解纏誤差的影響，MCF和SNAPHU方法的RMSE均較大；IM-PUMA的RMSE較前兩種方法精度有所提高，但仍然不是很理想；本文算法結果的RMSE在4種方法中最好的，這也進一步證明了本文算法在大梯度變化區域能夠取得較已有方法更好的結果。
圖5
圖5 真實資料試驗不同相位解纏方法的解纏結果和誤差
Fig. 5 Unwrapping results and error plots for different phase unwrapping methods with real data
表2 真實TanDEM-X資料的評估結果Tab. 2 Evaluation results of real TanDEM-X data
| 相位解纏方法 | RMSE |
|---|---|
| MCF | 5.387 0 |
| SNAPHU | 4.952 0 |
| IM-PUMA | 2.866 3 |
| 本文算法 | 2.442 8 |
新視窗打開| 下載下傳CSV
4 總結與展望
大梯度變化區域相位解纏是影響高精度DEM反演的關鍵問題之一。本文基于IM-PUMA方法，主要針對其算法模型權值設定不合理的問題進行研究，提出了一種基于勢函數改進權值的最大流/最小割相位解纏算法。該算法結合外部已有DEM的先驗相位梯度資訊，提出了一種基于視窗最大絕對相位梯度估計的權值設定方法，并采用拟凸勢函數公式獲得權值，使權值的設定與勢函數比對更為合理；此外，針對較極端大梯度變化區域的解纏誤差問題，本文采用對權值設定門檻值的方式，解決了極端大梯度變化區域權值設定過小勢函數無法發揮作用引起的解纏錯誤問題，進而提高大梯度變化區域的相位解纏精度。利用仿真資料和真實TanDEM-X InSAR資料對算法進行比較分析，試驗結果表明，在仿真資料中，與已有的相位解纏方法相比，本文算法在大梯度變化區域可以獲得更高精度的結果，解纏結果精度可以提高44.93%，較好地解決了IM-PUMA權值設定不合理引起的解纏誤差問題；真實TanDEM-X InSAR資料進行中，由于大梯度變化問題，已有的解纏算法均無法獲得理想的解纏結果，尤其是MCF産生了大範圍的解纏誤差，本文算法采用更為科學的權值設定方法，有效解決了大梯度變化區域的解纏誤差問題，試驗結果表明，本文算法可以獲得較已有方法更好的解纏結果。
初審：張豔玲複審：宋啟凡
終審：金 君
資訊