質量圖引導的干涉合成孔徑雷達圖像相位解纏算法

郝紅星， 于榮歡， 張喜濤

(1. 裝備學院 復雜電子系統仿真實驗室, 北京 101416； 2. 裝備學院 研究生管理大隊, 北京 101416)

質量圖引導的干涉合成孔徑雷達圖像相位解纏算法

郝紅星1， 于榮歡1， 張喜濤2

(1. 裝備學院 復雜電子系統仿真實驗室, 北京 101416； 2. 裝備學院 研究生管理大隊, 北京 101416)

干涉合成孔徑雷達測量在數字高程數據獲取等遙感領域得到廣泛應用，相位解纏是干涉測量中的關鍵步驟。運用基于質量圖的馬爾科夫隨機場建模方法將相位解纏問題轉化為能量最小化問題，通過序列樹權值信息傳遞算法對所建立的馬爾科夫隨機場模型進行求解。論文提出通過最小方向線性插值方法對相位圖像中質量較差的區域進行修復，即根據質量較差的像素點所處的位置，利用最短的方向對質量差的像素點進行插值修復，以獲得完整的解纏結果。實驗結果表明：相對于不采用質量圖的方法以及Goldstein枝切法等解纏算法，所提出的算法能夠完成復雜真實地形對應的相位圖像解纏。

遙感信息；干涉合成孔徑雷達；相位解纏；數字高程模型；馬爾科夫隨機場

干涉合成孔徑雷達(InSAR)測量,已經廣泛應用于地形高程數據快速測量、沉降分析、海洋監測等。在軍事應用中，干涉合成孔徑雷達獲得的相位圖像相對于其他測量方式，能夠比較快速完成地形高程數據的測量。美國于21世紀初利用干涉測量方法開始了全球數字高程數據測量，“奮進號”航天飛機上的SRTM(Shuttle Radar Topography Mission)系統在234 h的全球性作業中獲得了地球北緯60°～南緯56°的干涉相位圖像，面積超過1.19×109km2[1]。近些年，各國相繼發射星載干涉合成孔徑雷達用于陸地測圖、區域性觀測和資源調查等，主要包括德國的TerraSAR-X/TanDEM-X、歐空局的Sentinel-1 A/B、意大利的Cosmo-Skymed星座項目、日本的ALOS/ALOS-2，以及我國的遙感衛星十五號、高分三號。

由于測量機制的局限性，干涉測量獲取的相位圖像并不能直接使用。通過合成孔徑雷達所獲取的干涉相位圖像由于地表反射率差別，天氣狀況、電子設備噪聲以及前序傳輸過程的影響，存在大量噪聲。這些噪聲導致在相位圖像中存在殘差點(相位不連續的點)，而殘差點影響最終的數字高程數據的正確性。另外，獲取的相位圖像是真實相位(也稱為絕對相位)的模2π值，并不能直接應用于高度信息的判讀。由于測量數據的這種周期特性，使得相位圖像的處理成為非常困難的逆問題。

為了將相位圖像處理逆問題轉化為可以解決的形式，進而獲取能夠使用的高程信息，大多數算法將絕對相位(與真實地形高程成比例的相位值)的估計分成2個階段：第一個階段，進行干涉相位的估計，即將獲得的相位圖像進行降噪，獲得降噪后的纏繞相位值；第二個階段，將第一個階段獲取的降噪后的纏繞相位進行絕對相位估計，這一階段又被稱為相位解纏或相位展開，本論文重點關注該階段的算法研究。

相位解纏是干涉相位圖像數據處理過程中的關鍵步驟。最初的相位解纏算法為多項式曲面建模的方法，例如文獻[2]中提出了運用低次多項式曲面近似進行相位解纏的方法。但是在實際的應用中，一個多項式曲面并不能精確地描述整個相位曲面，因此需要將相位平面分割為不同的部分，然后對每一部分運用不同的參數曲面模型進行建模。這種建模方式復雜度較高且解纏效果也不理想。目前，相位解纏算法主要可分為2大類：路徑跟隨算法和基于最優化的算法。路徑跟隨算法主要包括Goldstein枝切算法[3-4]、基于置信圖的枝切算法[5]、Flynn的最小不連續性算法[6]、路線積分法[7]等。該類算法的主要缺點是容易受到噪聲的影響，并且需要滿足Itoh條件，即任意2個像素點之間的相位差小于π，但是實際情況中具有非連續區域的相位圖像由于條紋相干性比較差，并不能很好地滿足該條件?；谧顑灮姆椒ㄖ饕ㄗ钚范數相位解纏[8]、最小二乘法算法[9]25和集群智能算法[10]5502等。其中，最小p范數算法并不能很好地保留相位圖像中的不連續區域，Guo等[9]25提出通過改進最小二乘法進行相位解纏，能夠有效改善解纏精度和解纏速度，但是其解纏過程中忽略了非連續相位導致的相干性較差區域的解纏，因此只能適用于連續平滑相位的解纏過程；Maciel等[10]5504提出基于集群智能的相位解纏算法，但算法實施過程中需要對智能策略進行選擇。

現有的相位解纏算法針對連續相位能夠取得較好的解纏效果，但是對于真實數據，其解纏效果較差，尤其對于相干性較差區域的干涉相位圖像不能正確地進行解纏。本文研究了一種基于質量圖和馬爾科夫隨機場建模的相位解纏算法，有效結合了路徑跟隨和最優化算法的優勢。文獻[11]中也研究了基于質量圖的馬爾科夫隨機場模型應用于InSAR圖像的處理，但是該文獻中將對應模型應用于相位圖像的降噪問題，與本論文中解決的相位解纏問題是有區別的。

1 相位解纏問題及其馬爾科夫隨機場建模

在實際應用中，干涉合成孔徑雷達通過獲取測量信號的相位差得到測量的相位圖像，由于測量信號的周期特性，測量的相位圖像并不與實際測量曲面幾何對應，而是實際信號的模2π值，即在干涉圖像I中，干涉測量相位φ2π(p)與實際相位φ(p)之間的關系滿足

φ(p)=φ2π(p)+2k(p)π,?p∈I

(1)

式中，k(p)為像素點p處的解纏標簽值。相位解纏的過程為確定k(p),使其解纏后的相位圖像能夠滿足局部平滑性的特點。對于相位圖像φ2π中的相鄰像素點p和q，定義平滑性準則Epq(φ(p)-φ(q))為相鄰像素點差值絕對值的連續單調增函數，忽略噪聲影響，可將相位解纏問題建模為

(2)

式中，(p,q)∈e,標志p和q為相位圖像中的相鄰像素點，本文中e(馬爾科夫隨機場圖中的邊)的確定與質量圖有關。公式(2)是馬爾科夫隨機場的簡化形式，根據公式(1)對其進行進一步簡化，得到

2π(k(p)-k(q)))

(3)

式中:Z為整數的集合。相位解纏的過程為針對固定像素點p，確定整數值標簽k(p),使得式(3)取最小值的過程。

通常情況下，相位解纏問題中的e為相鄰像素點，這種方式導致相位圖像中的非連續區域在求解最小化過程中被平滑而獲得錯誤的解纏結果。本文將質量圖引入馬爾科夫隨機場模型，在解纏過程中能有效防止相位圖像中質量較差的區域誤差積累，從而獲得正確的絕對相位。e的確定過程如圖1所示。

1) 根據質量圖Q及設定的閾值η,將質量圖進行二值化得到Qb,即

(4)

2) 根據二值化質量圖確定馬爾科夫隨機場模型中邊的集合e。若像素p的四鄰域對應的Qb都為高質量區域，如圖1a)所示，則在馬爾科夫隨機場模型中考慮其4條邊；若像素p的四鄰域中有某一方向為質量較差的區域，如圖1b)所示，則在馬爾科夫隨機場中只考慮其3條邊；類似，其他2種情況邊的建模如圖1c)～圖1d)所示。

3) 若質量圖中出現孤立點，如圖1e)所示，則將質量圖中該點的質量置為0。

下面對相位解纏對應的馬爾科夫隨機場最優化問題進行求解。

圖1 待建模像素周圍的二值化質量圖情況(陰影像素為待建模像素)

2 馬爾科夫隨機場問題求解

本文運用序列樹權值信息傳遞算法求解相位解纏問題對應的馬爾科夫隨機場問題，前文已運用質量圖和平滑度相結合的方式對相位解纏問題進行了概述。通過分析可得，所建立的模型是只含有二元項的馬爾科夫隨機場問題的簡化形式，并提出運用信息傳遞算法求解，其算法步驟為：

第1步：初始化每個頂點的標簽值為k(p)=0,傳遞信息值為Mpq=0,目標函數下界最大值為Elowbound=0；

第2步：進行正向最優化，對于相位圖像中的每一個像素點(i,j)，計算該像素點右邊和下邊2條邊的信息傳遞，即Epq包括像素點(i,j)與(i+1,j)和像素點(i,j)與(i,j+1)2條邊對應的二元值。

(5)

式中：st為傳遞信息值的索引；Mst;l為標簽值l對應的傳遞信息值修正，其定義為Mpq;l-cMn;l(c為修正參數，Mn;l為像素點(i,j)的鄰域傳遞信息值之和)。

第3步：進行逆向最優化，對于相位圖像中的每一個像素點(i,j)，計算該像素點左邊和上邊2條邊的信息傳遞，即Epq包括像素點(i,j)與(i-1,j)和像素點(i,j)與(i,j-1) 2條邊對應的二元值,并計算目標函數下界最大值

(6)

第4步：從圖像左上角像素點開始，計算每個像素點的標簽值k(p),并計算當前目標函數值Eupbound。

第5步：若Eupbound-Elowbound<ε,則終止算法;否則，令Elowbound=0,并轉第2步。

若馬爾科夫隨機場中二元項為凸函數，則上述算法能夠收斂于最優解。為了防止一般情況下二元項為非凸時的算法無限迭代，可以設置最大迭代次數Imax。但在相位解纏問題中，由于本文算法通過質量圖來避免連續區域的平滑，因此可以設置二元項為凸函數。

3 低質量區域修復算法

圖2 質量低區域修復示意圖(陰影區域質量圖為0)

通過第2部分算法，對于質量圖為1的像素點能夠獲得正確的標簽值，因此可以獲得正確的解纏值，從而獲得絕對相位。但對于質量圖中值0對應的像素點，由于在建模過程中忽略了對應邊，因此需要單獨考慮。另外，由于質量圖中值為0的區域不可用信息較少，因此本文通過最短距離插值方法對該區域進行修復獲得絕對相位。根據待修復像素點所處的位置可將修復過程具體分3種情況：

1) 若像素點所在的區域，橫向與縱向質量圖值為0的像素點個數存在最小值且沒有到達圖像邊緣，則采用在該方向進行線性插值來進行修復。如圖2中像素點r，其橫向待修復點個數為4，縱向待修復點個數為3，因此其通過縱向進行修復，即φ(r)=L(φ(a),φ(e)),其中L(·)為線性插值函數。

2) 若像素點所在區域，橫向與縱向質量圖值為0的像素點個數相等且沒有到達圖像邊緣，則采用雙線性插值方法進行修復。如圖2中像素點q，其橫向和縱向待修復點的個數都為3，則通過雙線性插值進行修復，即

φ(q)=(L(φ(a),φ(e))+L(φ(d),φ(h)))/2

3) 若像素點所在區域，某一方向質量圖為0的像素點到達圖像邊緣，如圖2中像素點m，雖然其縱向待修復像素點個數為3，但由于其到達圖像邊緣，因此選擇橫向進行修復，即φ(m)=L(φ(c),φ(i))。若2個方向的待修復像素點都到達邊緣，因可利用信息太少，不對該區域進行修復。

4 實驗結果

為驗證所提出的算法，采用文獻[12]474所提供數據進行試驗驗證。該數據為美國Long’sPeak(Colorado)地區附近山地高程數據對應的干涉測量相位圖像，存在大量的懸崖等陡峭地形，相位圖像的相干性比較差，如圖3所示。算法中所使用的質量圖通過文獻中提供的相關系數圖可以獲得。在實驗中，采用的平滑函數為

實驗中，將本文提出算法與不采用質量圖引導的算法、Goldstein枝切法[3]和最小二乘算法[12]248進行比較。從圖3e)可以得出：本文算法能夠收斂到滿意解。圖3d)為本文算法的解纏結果，相對于不采用質量圖的方法，本文方法能對非連續區域取得較好的降噪結果。圖3g)和圖3h)分別為采用Goldstein枝切法和最小二乘算法的解纏結果。

圖4為對另一塊復雜地形相位圖像的解纏結果比較，對于圖4b)中矩形框中的非連續區域，本文算法能夠較好地解纏，如圖4c)所示；但不采用質量圖的方法不能對該區域正確解纏，如圖4d)所示。

本文定義均方誤差來對解纏結果進行定量分析，解纏結果的均方誤差定義為

圖3 真實地形解纏結果比較圖(區域1)

圖4 真實地形解纏結果比較圖(區域2)

算法圖3a)相位圖像圖4a)相位圖像Goldstein枝切法0.04700.0112最小二乘算法0.09820.0266非質量圖引導的算法0.08580.0158本文算法0.00700.0097

綜上實驗，本文的算法能夠對復雜相位圖像得到正確的解纏結果，并且對原始圖像中的非連續區域具有一定的保留作用。

5 結 束 語

本文運用質量圖與馬爾科夫隨機場建模相結合的方式對復雜相位解纏問題進行求解，將基于馬爾科夫隨機場建模方法應用于相位解纏問題中，所研究的建模和求解算法對于真實的復雜相位圖像能夠取得較好的解纏效果。特別是文中提出的通過最小方向線性插值方法對相位圖像中質量較差的區域進行修復的算法，能有效修補由于質量圖較差導致解纏后絕對相位的不連續性。但是，算法的效果和效率與質量圖的選取有一定關系，下一步的研究應側重于如何生成或者選擇質量圖以取得更加精確的解纏結果。

)

[1]SEAL D,ROGEZ F.SRTM As-flown mission timeline[R].Pasadena,California：JPL NASA,2000：1.

[2]FRIEDLANDER B,FRANCOS J M.Model based phase unwrapping of 2-d signals [J].IEEE Transactions on Signal Processsing，1996(44):2999-3007.

[3]曾凡光，吳光敏，MAI J D，等.Goldstein枝切法對存在間斷相位缺陷的解纏研究[J].激光技術，2014，38(3)：335-341.

[4]HUANG Q,ZHOU H,DONG S,et al.Parallel branch-cut algorithm based on simulated annealing for large-scale phase unwrapping [J].IEEE Transactions on Geoscience and Remote Sensing,2015，53(7):3833-3846.

[5]JIAN G,Reliability-map-guided phase unwrapping method [J].IEEE Geoscience and Remote Sensing Letters,2016，13(5):716-720.

[6]王軍飛，彭軍還，楊紅磊，等.改進金分發的InSAR相位解纏算法[J].測繪科學，2016，41(12):85-88.

[7]于宏，王成恩，于嘉鵬，等.基于粒子群算法的復雜產品裝配序列規劃[J].東北大學學報,2010,31(2):261-264.

[8]LU Y G,WANG X Z,ZHANG X P.Weighted least-squares phase unwrapping algorithm based on derivative variance correlation map[J].Optik,2007(118):62-66.

[9]GUO Y,CHEN X T,ZHANG T.Robust phase unwrapping algorithm based on least squares[J].Optics and Lasers in Engineering，2014，63:25-29.

[10]MACIEL L D S，ALBERTAZZI A G.Swarm-based algo-rithm for phase unwrapping[J].Applied Optics,2014,53(24):5502-5509.

[11]李泓宇,宋紅軍,王輝,等.質量圖指導的改進的基于 CMRF 模型的自適應 InSAR 相位濾波[J].科學技術與工程,2015 (25):55-60.

[12]GHIGLIA D ,PRITT M.Two-dimensional phase unwrapping:theory,algorithms,and software[M].Hoboken,NJ:Wiley,1998:248-474.

(編輯：李江濤)

A Phase Unwrapping Method of Interferometric SAR Phase Image Based on Quality Map Guided Model

HAO Hongxing1， YU Ronghuan1， ZHANG Xitao2

(1. Science and Technology on Complex Electronic System Simulation Laboratory, Equipment Academy, Beijing 101416,China； 2. Department of Graduate Management, Equipment Academy, Beijing 101416, China)

Interferometric synthetic aperture radar (SAR) has been widely used in remote sensing field, such as digital elevation data acquisition, and phase unwrapping is one of the key steps in interferometric measurement. The paper transfers phase unwrapping to the minimization of an energy function based on the quality map guided Markov Random Field model, and solves the Markov Random Field problem with the algorithm of sequence tree weight message transfer. The paper also proposes to restore the pixels with phase image of poor quality by minimum direction linear interpolation, that is, based on the position of the pixels of poor quality, the interpolation repair is made to them in the shortest direction, to achieve complete unwrapping. Experiments show that the proposed algorithm can unwrap complex phase images generated by the varied real terrain, compared with the unwrapping algorithm like Goldstein branch cut method and methods without the quality map.

remote sensing information; interferometric synthetic aperture radar (InSAR); phase unwrapping; digital elevation model; Markov Random Field

2017-04-17

部委級資助項目

郝紅星(1987—)，男，助理研究員，博士，主要研究方向為多媒體與虛擬現實、雷達信號處理。hongxinghao87@126.com

TN957； P237

2095-3828(2017)03-0008-06

A DOI 10.3783/j.issn.2095-3828.2017.03.002