郝紅星， 于榮歡， 張喜濤

(1. 裝備學(xué)院 復(fù)雜電子系統(tǒng)仿真實(shí)驗(yàn)室, 北京 101416； 2. 裝備學(xué)院 研究生管理大隊(duì), 北京 101416)

質(zhì)量圖引導(dǎo)的干涉合成孔徑雷達(dá)圖像相位解纏算法

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

(1. 裝備學(xué)院 復(fù)雜電子系統(tǒng)仿真實(shí)驗(yàn)室, 北京 101416； 2. 裝備學(xué)院 研究生管理大隊(duì), 北京 101416)

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

遙感信息；干涉合成孔徑雷達(dá)；相位解纏；數(shù)字高程模型；馬爾科夫隨機(jī)場

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

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

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

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

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

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

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

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

(1)

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

(2)

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

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

(3)

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

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

1) 根據(jù)質(zhì)量圖Q及設(shè)定的閾值η,將質(zhì)量圖進(jìn)行二值化得到Qb,即

(4)

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

3) 若質(zhì)量圖中出現(xiàn)孤立點(diǎn)，如圖1e)所示，則將質(zhì)量圖中該點(diǎn)的質(zhì)量置為0。

下面對相位解纏對應(yīng)的馬爾科夫隨機(jī)場最優(yōu)化問題進(jìn)行求解。

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

2 馬爾科夫隨機(jī)場問題求解

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

第1步：初始化每個(gè)頂點(diǎn)的標(biāo)簽值為k(p)=0,傳遞信息值為Mpq=0,目標(biāo)函數(shù)下界最大值為Elowbound=0；

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

(5)

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

第3步：進(jìn)行逆向最優(yōu)化，對于相位圖像中的每一個(gè)像素點(diǎn)(i,j)，計(jì)算該像素點(diǎn)左邊和上邊2條邊的信息傳遞，即Epq包括像素點(diǎn)(i,j)與(i-1,j)和像素點(diǎn)(i,j)與(i,j-1) 2條邊對應(yīng)的二元值,并計(jì)算目標(biāo)函數(shù)下界最大值

(6)

第4步：從圖像左上角像素點(diǎn)開始，計(jì)算每個(gè)像素點(diǎn)的標(biāo)簽值k(p),并計(jì)算當(dāng)前目標(biāo)函數(shù)值Eupbound。

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

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

3 低質(zhì)量區(qū)域修復(fù)算法

圖2 質(zhì)量低區(qū)域修復(fù)示意圖(陰影區(qū)域質(zhì)量圖為0)

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

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

2) 若像素點(diǎn)所在區(qū)域，橫向與縱向質(zhì)量圖值為0的像素點(diǎn)個(gè)數(shù)相等且沒有到達(dá)圖像邊緣，則采用雙線性插值方法進(jìn)行修復(fù)。如圖2中像素點(diǎn)q，其橫向和縱向待修復(fù)點(diǎn)的個(gè)數(shù)都為3，則通過雙線性插值進(jìn)行修復(fù)，即

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

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

4 實(shí)驗(yàn)結(jié)果

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

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

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

本文定義均方誤差來對解纏結(jié)果進(jìn)行定量分析，解纏結(jié)果的均方誤差定義為

圖3 真實(shí)地形解纏結(jié)果比較圖(區(qū)域1)

圖4 真實(shí)地形解纏結(jié)果比較圖(區(qū)域2)

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

綜上實(shí)驗(yàn)，本文的算法能夠?qū)?fù)雜相位圖像得到正確的解纏結(jié)果，并且對原始圖像中的非連續(xù)區(qū)域具有一定的保留作用。

5 結(jié) 束 語

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

)

[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].激光技術(shù)，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]王軍飛，彭軍還，楊紅磊，等.改進(jìn)金分發(fā)的InSAR相位解纏算法[J].測繪科學(xué)，2016，41(12):85-88.

[7]于宏，王成恩，于嘉鵬，等.基于粒子群算法的復(fù)雜產(chǎn)品裝配序列規(guī)劃[J].東北大學(xué)學(xué)報(bào),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]李泓宇,宋紅軍,王輝,等.質(zhì)量圖指導(dǎo)的改進(jìn)的基于 CMRF 模型的自適應(yīng) InSAR 相位濾波[J].科學(xué)技術(shù)與工程,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

部委級(jí)資助項(xiàng)目

郝紅星(1987—)，男，助理研究員，博士，主要研究方向?yàn)槎嗝襟w與虛擬現(xiàn)實(shí)、雷達(dá)信號(hào)處理。hongxinghao87@126.com

TN957； P237

2095-3828(2017)03-0008-06

A DOI 10.3783/j.issn.2095-3828.2017.03.002