張斌，胡慶榮，李爽，韋立登

(1.北京無線電測量研究所，北京 100854;2.中國航天科工集團(tuán)有限公司 第二研究院，北京 100854)

0 引言

InSAR技術(shù)通過利用2幅或多幅SAR圖像的相位信息可以獲取大范圍、高精度的地表三維信息。其中相位解纏是InSAR處理的關(guān)鍵步驟，其處理精度對DEM(digital elevation model)產(chǎn)生重要影響[1-2]。傳統(tǒng)的單基線相位解纏，大都假設(shè)Itoh條件滿足的前提下進(jìn)行的，其在相位噪聲大和陡峭地形情況下很難滿足高精度解纏要求，甚至失效，無法準(zhǔn)確還原地形[3]。而多基線相位解纏算法能打破Itoh條件限制，獲取更多的觀測區(qū)的信息，融合多幅干涉圖，減小干涉相位的模糊區(qū)間，降低解纏難度[4]。

近年來，多基線InSAR相位解纏得到了長足發(fā)展，Ghiglia[5]等提出多基線相位解纏繞方法的最小二乘法和最大似然函數(shù)法，F(xiàn)erraiuolo[6]等提出基于馬爾可夫隨機場的最大后驗概率估計方法，夏香根[7]等提出一種基于中國余數(shù)定理方法相位解纏，于瀚雯[8]提出基于聚類分析的多基線相位解纏繞算法。而在多基線最大后驗概率估計方法中，F(xiàn)erraioli[9]等提出了基于圖割的多基線相位解纏，其直接由多幅干涉條紋圖直接計算高程值，從而避開了相位解纏過程。但該類方法由于采用Ishikawa圖割[10]，其所建立的高度量化值越多，其相位解纏所用的內(nèi)存和處理所花的時間越大，因此A.shabou[11]等提出了一種基于二分法的多基線相位解纏方法，其所占用的內(nèi)存減少，但其每一個循環(huán)都要訪問所有的高度量化值，相應(yīng)的計算時間卻增加。相比直接估計圖像的高程值，用Ishikawa方法估計相位解纏所需的纏繞數(shù)(2π的整數(shù)倍)所花的內(nèi)存和處理時間無疑更占優(yōu)勢，而為了進(jìn)一步節(jié)省內(nèi)存和減少處理時間，本文據(jù)此提出了基于迭代的多基線圖割方法，在不影響全局最優(yōu)的基礎(chǔ)上通過對圖像每個點所需的纏繞數(shù)進(jìn)行分段，通過優(yōu)化迭代，使其最終達(dá)到最優(yōu)的纏繞數(shù)。

1 圖割模型

干涉復(fù)圖像可以表示為

φ=Aejφ+n，

(1)

式中：φ為真實干涉相位；n為零均值的復(fù)圓周高斯噪聲；φ為帶有噪聲的干涉相位。

根據(jù)貝葉斯準(zhǔn)則，基于最大后驗概率的相位解纏模型可表示為

(2)

式中：P(φ|φ)為在實際干涉相位下真實相位的概率;P(φ)為實際相位的概率，此處為一常數(shù);P(φ|φ)為假設(shè)已知真實相位的情況下實際相位的概率，即數(shù)據(jù)的似然函數(shù);g(φ)為真實干涉相位的先驗概率密度函數(shù)。

假設(shè)所有的像素點之間是相互獨立的，可知

(3)

式中：i為相應(yīng)的通道(基線)，此處共有C個基線;P為每個通道內(nèi)的元素，假設(shè)數(shù)據(jù)為M×N個元素；γi,p為第i個通道內(nèi)第p個像素的相干系數(shù)值，φi,p為對應(yīng)的實際相位值。

在一階馬爾可夫(MRF)模型中，相位的先驗概率密度函數(shù)可表示為

(4)

對式(2)取對數(shù)后，基于MRF模型的圖割相位解纏算法可表示為

(5)

式中：

2 基于迭代的圖割相位解纏算法

2.1 算法流程

基于式(5)的圖割相位解纏算法主要分為2種：第1種為Boykov等[12]提出的α-β移動迭代接近最優(yōu)值，算法處理速度快，但在噪聲比較大情況下或者欠采樣(地勢陡峭地區(qū)等)其整體最優(yōu)的魯棒性下降，局部解纏性能不佳;第2種方法為全局算法[9，13-14]，該算法將每個點可能所需的纏繞數(shù)作為標(biāo)簽集合，然后運用圖割算法在標(biāo)簽集合里尋找到整體最優(yōu)值，從而進(jìn)行相位解纏，此方法為整體最優(yōu)算法，魯棒性較高，但其占用內(nèi)存大，計算時間也很大，而A.Shabou[11]針對此問題提出了二分法的相位解纏來解決此問題，此算法雖然占用內(nèi)存變小但其運算時間大大增加。本文針對Ishikawa算法占用內(nèi)存大與計算時間長的缺點，提出一種基于迭代的多基線圖割算法，算法在不影響整體最優(yōu)的情況下，將每個點可能所需的纏繞數(shù)進(jìn)行初始分段，然后利用圖割算法，不斷逼近整體最優(yōu)值。

算法的流程如下：

(1) 估計整幅圖的纏繞數(shù)，選取合適的分段(分段間隔為2π)。

(2) 在分段內(nèi)，更新相位值，構(gòu)建改進(jìn)的Ishikawa圖網(wǎng)絡(luò)模型。

2.2 圖網(wǎng)絡(luò)的建立

圖網(wǎng)絡(luò)是圖割算法的基礎(chǔ)，其將每一條邊賦予相應(yīng)權(quán)值并進(jìn)行計算，利用網(wǎng)絡(luò)圖的最大流最小割來尋求當(dāng)前迭代中的最優(yōu)值。

圖1單向箭頭表示數(shù)據(jù)邊或約束邊，雙向箭頭表示交互邊，p,q,r表示數(shù)據(jù)的元素(p,q為相鄰的元素)，割集如圖中Cut粗體線段所示。

圖1中的頂點集合為干涉圖中像素點的數(shù)目再加上S,T2個節(jié)點，即

i=1,2,…,L}∪{s}∪{t}，

(6)

式中：L為算法中所選取的分段的數(shù)目;M為干涉圖的數(shù)；N為干涉圖的列。

圖1的有向邊集可分為4個部分，分別為數(shù)據(jù)邊εD，約束邊εC，交互邊εI，補充邊εA。

數(shù)據(jù)邊如圖1中向上的箭頭表示，其代表數(shù)據(jù)的似然函數(shù)項，即

(7)

(8)

數(shù)據(jù)邊的權(quán)值為

(9)

(10)

約束邊如圖1中向下箭頭表示是為了在進(jìn)行最大流最小割時，每個像素有且僅有一個最優(yōu)值，此時有

(11)

(12)

約束邊的權(quán)值表示為

(13)

交互邊如圖1中雙向箭頭表示，體現(xiàn)能量函數(shù)的先驗信息，此時有

(14)

(15)

在交互邊中為了達(dá)到公式(5)的平衡，需要加入補充邊，此時有

(16)

(17)

在此可以將εA對應(yīng)的權(quán)值cp,q(i)和cp,q(j)加到εD對應(yīng)邊的權(quán)值去。

通過構(gòu)建網(wǎng)絡(luò)圖模型將能量最小問題轉(zhuǎn)化為求解網(wǎng)絡(luò)圖的最小割，如圖1中粗體Cut所示，得到圖的最大流最小割的同時也就得到了當(dāng)前迭代中所有像素點的纏繞數(shù)，再依據(jù)算法流程向下執(zhí)行。

3 實驗校驗

本文提出的方法在占用內(nèi)存和計算時間方面均有一定程度的改進(jìn)，為了對上述算法進(jìn)行驗證，本節(jié)通過仿真中進(jìn)行了說明。仿真實驗在Matlab上運行，其中圖割算法中的最大流最小割算法用C++實現(xiàn)[15]。

用三基線進(jìn)行仿真實驗，其中長基線原始高斯面如圖2所示,像素大小為64×64,對原始相位圖圖2加入相干系數(shù)為0.9的高斯白噪聲。圖3為加入高斯白噪聲后對應(yīng)的纏繞圖，從纏繞圖中可以看出該高斯面出現(xiàn)了較大面積的欠采樣情況，即打破了Itoh條件，此時用傳統(tǒng)的單基線相位解纏方法難以得到較好的結(jié)果，因此適合采用多基線相位解纏方法。相應(yīng)的短基線纏繞相位圖如圖4所示，長基線與短基線的比例為7：3，本次仿真實驗與文獻(xiàn)[6]中的Ishikawa方法進(jìn)行了對比，實驗結(jié)果如圖5，6所示。

對相位解纏前后的誤差(與原始相位之差)進(jìn)行分析比較，則有error迭代=5.922 8 rad,errorIshikawa=-232.477 9 rad，迭代算法的誤差較小，更為接近原始相位，這與迭代算法分段后能更好地兼顧局部最優(yōu)有關(guān)系。從仿真可以看出迭代算法能打破Itoh條件限制，較好地處理相位變化較快的區(qū)域，且運算時間相比Ishikawa算法有所減少。

4 結(jié)束語

本文提出了一種基于迭代的多基線InSAR圖割相位解纏算法。該算法針對Ishikawa圖割方法處理相位解纏時占用內(nèi)存大、計算時間長的問題，通過對纏繞數(shù)進(jìn)行分段，采取迭代的方式逐步逼近最優(yōu)值。算法通過構(gòu)建能量函數(shù)模型，利用一階馬爾可夫模型的先驗信息，構(gòu)建改進(jìn)的Ishikawa圖網(wǎng)絡(luò)模型，通過最大流最小割算法及最小化能量函數(shù)準(zhǔn)則逐步優(yōu)化相位值，從而完成相位解纏。仿真實驗驗證該迭代算法能夠減少相位解纏所需的內(nèi)存和處理時間。實驗是在低噪情況下進(jìn)行，對于噪聲比較大的區(qū)域，其解纏的魯棒性有待進(jìn)一步驗證。此外，對實測數(shù)據(jù)的處理也將在下一步展開。

參考文獻(xiàn)：

[1] BAMLER R,HARTL P. Synthetic Aperture Radar Interferometry[J].Inverse Problems,1998,14(4):1-54.

[2] 曲長文,侯海平,周強.合成孔徑雷達(dá)三維成像技術(shù)及特點評述[J].現(xiàn)代防御技術(shù),2010,38(4):110-115.

QU Chang-wen,HOU Hai-ping,ZHOU Qiang.An Overview of Three-Dimensional Imaging Technology and Technical Characteristic of SAR[J].Modern Defence Technology,2010,38(4):110-115.

[3] LOTH K.Analysis of the Phase Unwrapping Problem[J].Applied Optics,1982 ,21(14):2470.

[4] FERRAIUOLO G,PASCAZIO V,SCHIRINZI G.Maximum a Posteriori Estimation of Height Profiles in InSAR Imaging[J].IEEE Geoscience and Remote Sensing Letters,2004,1(2):66-70.

[5] GHIGLIA D C,WAHL D E.Interferometric Synthetic Aperture Radar Terrain Elevation Mapping from Multiple Observations[C]∥Proceeding of IEEE Digital Signal Processing Workshop IEEE Press,1994:33-36.

[6] FERRAIUOLO G,MEGLIO F,vito Pascazio,et al.DEM Reconstruction Accuracy in Multichannel SAR Interferometry[J].IEEE Transaction on Geoscience and Remote Sensing,2009,47(1):191-201.

[7] XIA Xiang-gen,WANG Gen-yuan.Phase Unwrapping and a Robust Chinese Remainder Theorem[J].IEEE Signal Processing Letters,2007,14(4):247-250.

[8] YU Han-wen,LI Zhen-fang,BAO Zheng.A Cluster-Analysis-Based Efficient Multibaseline Phase Unwrapping Algorithm[J].IEEE Transaction Geoscience and Remote Sensing,2011,49(1):478-487.

[9] FERRAIOLI G,SHABOU A,TUPIN F,et al.Multichannel Phase Unwrapping with Graph Cuts[J].IEEE Geoscience and Remote Sensing Letters,2009,6(3):562-566.

[10] ISHIKAWA H.Exact Optimization for Markov Random Fields with Convex Priors[J].IEEE Trans Pattern Anal Mach Intell,2003,25(10):1333-1336.

[11] SHABOU A,DARBON J,TUPIN F.Multilabel Partition Moves for MRF Optimization[J].Image and Vision Computing,2013,31(1):14-30.

[12] BOYKOV Y,VEKSLER O,ZABIH R.Fast Approximate Energy Minimization via Graph Cuts[J].IEEE Trans.Pattern Anal Mach Intell,2001,23(11):1222-1239.

[13] JOSé M,Dias Bioucas,Val?dao Gon?alo.Phase Image :Unwrapping and Denoising with Diversity and Multi-Resolution[J].IEEE Transaction on Geoscience and Remote Sensing,2001,33(3):579-589.

[14] Jerome Darbon,Marc Sigelle.Image Restoration with Discrete Constrained Total Variation Part II:Levelable Functions,Convex Priors and Non-Convex Cases[J].Journal of Mathematical Imaging and Vision,2006:277-291.

[15] BOYKOV Y,Vladimir Kolmogorov.An Emperiment Comparison of Min-Cut/Max-Flow Algorithms for Energy Minimization in Vision[J].IEEE Trans.Pattern Anal Mach Intell,2004,26(9):1124-1137.