相位分塊與擬合法結合的InSAR相位解纏算法

馬靚婷，盧小平，余振寶

(河南理工大學 自然資源部礦山時空信息與生態修復重點實驗室，河南 焦作 454000)

0 引言

合成孔徑雷達干涉技術(interferometric synthetic aperture radar，InSAR)是在主動遙感基礎下發展的對地觀測技術，其利用同一目標區域的SAR復圖像對共軛相乘得到干涉圖，根據干涉圖的相位值，計算出2次成像過程中信號產生的路程差，從而獲取監測區域高精度的三維地形信息和微小的地形形變信息[1-2]。因全天時、全天候、方便快捷等優點，InSAR技術被廣泛用來監測地面沉降、地裂縫、火山活動、地震形變等[3]。

相位解纏是InSAR處理的重要環節，其結果的優劣直接影響到地形測量的精度。受限于合成孔徑雷達的成像和處理方式，直接利用影像獲得的干涉圖一般含有較大的噪聲，局部相位殘差點的增多會形成不可靠數據斑塊，使該區域相位解纏出現漏解或錯解，導致InSAR圖像恢復失敗，從而影響到形變提取的精度。因此，提高相位解纏精度是InSAR處理中提高形變精度的重要環節[4]。

相位解纏是通過在解纏路徑上進行積分從而還原真實目標信息。當干擾因素少、相位質量高時，能很好地還原真實相位信息；當干擾因素多時，誤差會通過積分進行積累與傳播，得到的相位數據與真實數據會有較大的差異。現有相位解纏方法主要分為路徑跟蹤解纏法[5-6]、最小范數解纏法[7]和網格規劃解纏法[8]3大類。路徑跟蹤法是通過設置合適的積分路徑，將誤差限制在一定區域內，防止相位誤差全局傳遞，其涵蓋了經典的Goldstein枝切法、質量圖引導法、最小范數法等。枝切法是利用殘差點的連接得到枝切線，最后沿著枝切線進行積分得到解纏結果，但枝切法易出現孤島現象[5]。質量圖引導法是在質量圖的引導下確定積分路徑，這種算法對干涉圖質量要求較高[9-10]。最小范數法是將相位解纏轉換成數學上的最小范數問題，其常用的是最小二乘法，但這種方法穿過殘差點會造成誤差的全局傳遞[11]。網格規劃法是將相位解纏問題轉化為求解費用流的網絡優化問題，主要有最小費用流和統計費用流等，但這種方法噪聲會沿著積分路徑傳遞，使得解纏結果不理想[12-13]。

針對當前相位解纏算法在相位圖像存在嚴重噪聲時解纏效果較差的問題，本文提出了基于相位分區與擬合法結合的InSAR干涉圖相位解纏算法。該算法：首先對纏繞相位進行分塊，將相位在相同區間的相鄰像素進行合并；塊內像素屬于同一包裹次數，通過線性擬合求取合適的補償系數K進行塊間解纏；最后利用曲面擬合的方法進行殘差塊解纏，合并所有解纏塊得到最終的解纏結果。

1 基本理論

1.1 相位解纏繞

在進行InSAR數據處理時，將主輔影像共軛相乘并取相位信息即可得到復干涉條紋圖，但通過共軛相乘得到的相位差，與真實相位差相差2Kπ，即真實相位與纏繞相位的關系如式(1)所示。

Φ=φ+2Kπ

(1)

式中：Φ表示真實相位；φ表示纏繞相位；K表示補償系數，其值為整數。

相位解纏是從纏繞相位φ確定補償系數K值，進而估計真實相位Φ的過程。相位解纏必須兼顧一致性和精確性。一致性是指解纏后的相位數據矩陣中任意兩點間的相位差與其路徑無關；精確性是指解纏后的相位數據能真實地恢復原始相位信息[14]。

1.2 區域生長算法

區域生長算法的基本思想是將具有相似性質的像素合并到一起。對每一個區域要先指定一個種子點作為生長的起點，然后將種子點周圍鄰域的像素點和種子點進行對比，將具有相似性質的像素合并起來繼續向外生長，直到沒有滿足條件的像素被包括進來為止[15]。

2 本文算法

InSAR干涉圖雖然經過了濾波處理，但干涉圖中仍存在噪聲，噪聲使得纏繞相位出現殘差點和低相干區域，使干涉圖解纏結果不全面或解纏誤差較大。基于分塊與合并策略的相位解纏算法是一類高效的方法[16-18]，該類方法把整個纏繞相位圖分成若干塊并執行塊間的相位解纏，合并所有塊可得到最終的解纏結果。現有的解纏算法無法將誤差和孤島現象同時避免，因此本文提出針對InSAR干涉圖的分區與擬合法結合的相位解纏算法。

2.1 相位分區

本文將相位區域擴張法(prelude)中的相位分區用于InSAR干涉圖相位解纏中，當分區間隔較大時相位塊內會出現纏繞的現象，分區間隔較小時相位塊較為零碎，不利于殘差塊的識別，影響算法的效率和精度。因此，本文選擇以π/3為區間長度，對InSAR干涉圖纏繞的相位矩陣實行分區，將相位信息在同一個區域的相鄰像素點合并為像素塊。纏繞相位的相位區間都在(-π，π]間，根據選定的區間長度，將相位分為(-π，-2π/3]，(-2π/3，-π/3]，(-π/3，0]，(0，π/3]，(π/3，2π/3]，(2π/3，π]6個區間。遍歷整個相位矩陣，對其中相位屬于各個分區的像素點進行標記，如圖1(a)所示，將屬于同一分區且相連的像素合并為塊，得到若干相位處于同一區間的像素塊，如圖1(b)所示。合并像素采用四鄰域法，即像素周圍的4個與其相鄰的同區像素可以進行合并，如圖1(c)所示，塊內像素擁有同一包裹數K。設置閾值為50，像素塊大于等于50為正常塊，像素塊小于50為殘差塊。

圖1 相位解纏過程示意圖

2.2 擬合法解纏

在理想狀態下干涉圖不存在噪聲，相位梯度小于π，選擇一個點作為起始點，可直接進行積分解纏，但現實中噪聲、地形起伏和相干性較低等現象給相位解纏帶來了困難。利用線性擬合求取補償系數K值的方法對噪聲不敏感，在區域生長解纏繞中利用臨近的已解纏繞像素的相位信息求取K值可快速得到相鄰像素的真實相位信息，且能夠適應較大的噪聲水平和相位的快速變化等情況[19]。

擬合法相位解纏是基于相鄰像素真實相位變換不大于π，使用相位分區方法對相位圖像進行分塊，塊內相鄰像素的相位變化小于給定相位間隔，即相鄰像素之間的相位差小于π/3且塊內像素擁有同一包裹數K，將塊間相位解纏問題轉化為線性擬合法求K值。選取塊間相鄰的2列像素，如圖1(d)所示，其2列像素的相位關系為式(2)所示。

Y=X+2Kπ

(2)

式中：X表示已解纏像素塊的相位值；Y表示未解纏像素塊的相位值；K為補償系數，其值為整數。

根據式(2)可知k為

(3)

式中：k為線性擬合系數，認為其最接近的整數值為整數補償系數K；Xi表示已解纏像素的相位值；Yi表示未解纏像素的相位值。

纏繞相位塊之間的相位解纏采用區域生長方式進行。把塊作為一個解纏處理的基本單元，塊內像素擁有相同的包裹次數K，選擇第一個塊為起始塊，距離起始塊中心位置最近的塊為生長塊。生長塊的解纏繞即是找到一個最佳的相位包裹次數K，選擇起始塊與生長塊相臨近的2組像素，進行線性擬合，找到最優的補償系數進行相位解纏，并合并解纏后的2塊，之后進行下一個塊的解纏，直至所有正常的塊完成解纏與合并。

2.3 殘差塊曲面擬合

當所有的正常塊完成相位解纏后，利用曲面擬合的方法對殘差塊進行解纏。根據解纏后相位是連續曲面的原理，通過殘差塊周圍像素的真實相位對殘差塊進行解纏，將殘差塊內像素點10×10窗口內所有經過解纏的像素進行最小二乘曲面擬合，將其結果作為殘差像素的真實相位值，選擇第一個像素為起點利用區域生長算法將殘差塊內的像素解纏完畢。

2.4 本文算法描述

本文所提相位解纏算法主要分3個部分，即進行相位分區，將像素點分為若干個塊；進行塊間的相位解纏；殘差像素塊解纏。其算法的具體實現有以下4步。

1)將輸入的InSAR干涉圖纏繞的相位矩陣實行相位分區，將相位信息在同一個區間且相鄰的像素點合并為像素塊，將像素小于50的像素塊標記為殘差塊。

2)實現相鄰塊間的相位解纏。利用線性擬合法求取補償系數K進行相鄰塊間的相位解纏，根據區域生長算法將正常塊解纏完畢。

3)當正常像素塊全部合并完成后，對標記的殘差像素塊進行曲面擬合解纏。

4)合并所有解纏結果，完成整個區域的相位解纏工作。

3 實驗結果與分析

本文對模擬數據進行實驗并用實測數據驗證。針對InSAR相位的特性，從解纏的準確度和算法運行的時間對解纏結果進行評價，除主觀視覺評價外，用均方根誤差和算法運行時間對實驗進行客觀評價。

3.1 模擬實驗

本實驗為MATLAB仿真的地形圖，用基于雷達傳感器參數和軌道數據的方法模擬InSAR干涉圖[20]。首先模擬無噪相位圖，如圖2(a)所示，加入噪聲后進行相位纏繞形成纏繞相位，如圖2(c)所示,并將其中50像素×50像素大小的范圍加重噪聲，以驗證算法在高強度噪聲下的效果。

圖2 模擬數據示意圖

將本文算法與Goldstein枝切法、質量圖引導的路徑追蹤法、四向剪切最小二乘法解纏結果進行比較與分析。圖3(a)～圖3(l)為不同解纏方法結果。從目視效果看，Goldstein枝切法解纏結果在高強度噪聲區域出現了孤島現象，無論是曲面圖和俯視圖都存在明顯的錯誤。質量圖引導的路徑追蹤法解纏結果在噪聲較弱、相位質量好的區域解纏效果較好，在噪聲密集區域解纏結果不連續。四向剪切最小二乘法解纏結果相位曲面圖較為連續，但據其誤差直方圖可知解纏結果誤差較大。本文算法解纏結果與原曲面圖相位較為一致，在噪聲密集區域依然有好的效果且誤差值較小。

圖3 模擬數據實驗結果

表1是對模擬數據解纏結果的定量比較。從中可看出，四向剪切最小二乘法運算速度最快，但誤差全局傳遞造成其精度最低。質量圖引導法均方根誤差小于四向剪切最小二乘法，但其運行時間大于四向剪切最小二乘法。本文算法的解纏結果均方根誤差最小，相較于四向剪切最小二乘法減少了44%，相較于質量圖引導法減少了38%，且運行時間與質量圖引導法相當。

表1 不同方法的模擬干涉圖解纏結果定量比較

3.2 實測數據驗證

為驗證本文算法的有效性，利用山西平朔地區2011年12月17日和2012年2月27日的RadarSat-2實測數據進行驗證。對其數據進行配準、干涉等處理得到其真實相位干涉數據，截取其中500像素×500像素大小進行實驗處理，如圖4(a)所示。

圖4(b)～圖4(h)顯示了不同方法的解纏結果。從目視效果看，Goldstein枝切法出現了孤島現象，存在明顯的錯誤。質量圖引導法在低質量區域解纏效果不佳，其解纏結果俯視圖上出現多個解纏效果不連續的區域。四向剪切最小二乘法解纏結果較為連續，但誤差較大。本文算法結果較為光滑避免了很多尖峰毛刺現象。

圖4 實測數據實驗結果

表2是對實測數據解纏結果的定量比較。可以看出，Goldstein枝切法運算時間最長且解纏失敗，四向剪切最小二乘法運算時間最短但其均方根誤差最大。質量圖引導法運算時間較長，但其精度高于四向剪切最小二乘法。本文算法解纏結果精度最高且運算時間適中。

表2 不同方法的實測干涉圖解纏結果定量比較

4 結束語

本文提出了一種相位分區與擬合法結合的相位解纏算法，該算法首先將獲得的相位圖像分為多個塊，塊內相位都在給定的相位區間內，把像素個數小于給定閾值的塊歸類為殘余像素；然后利用區域生長的擬合方法，依次進行塊與塊之間相位解纏繞和殘余像素相位解纏繞，合并后得到最終的解纏結果。通過對模擬數據進行實驗并用實測數據驗證，實驗結果表明，該算法無論是目視效果還是定量指標分析均優于其他算法，減小了誤差傳播的范圍，提高了相位解纏的精度，但該算法受相位連續性約束面對迭掩區域和復雜地形相位跳變問題時仍存在較大的改進空間。