網(wǎng)絡(luò)流與曲面擬合結(jié)合的相位解纏方法

， ，

( 1.廣西無(wú)線寬帶通信與信號(hào)處理重點(diǎn)實(shí)驗(yàn)室， 廣西桂林 541004； 2.桂林電子科技大學(xué)計(jì)算機(jī)與信息安全學(xué)院， 廣西桂林 541004； 3.桂林電子科技大學(xué)信息與通信學(xué)院， 廣西桂林 541004)

0 引言

相位解纏是干涉合成孔徑雷達(dá)(Interference Synthetic Aperture Radar，InSAR)干涉測(cè)量的關(guān)鍵流程[1-3]，它的準(zhǔn)確程度直接影響到InSAR生成的數(shù)字高程模型的精確性[3]。干涉圖的相位展開方法有多種，主要分為4類：第一類是路徑跟蹤法[4-5]，第二類是最小范數(shù)法[6-7]，第三類是網(wǎng)絡(luò)流法[8]，第四類是其他方法。

路徑跟蹤法最常用的是枝切法。利用枝切法進(jìn)行相位解纏是Goldstein等首先提出的[4]。枝切法在相干性差、殘差點(diǎn)較多的區(qū)域難以選擇準(zhǔn)確的枝切線，容易形成一個(gè)個(gè)被枝切線包圍的孤立區(qū)間，即“解纏孤島”，積分路徑都無(wú)法到達(dá)從而無(wú)法解纏。“解纏孤島”由于沒(méi)有明顯的聯(lián)系,單獨(dú)解纏時(shí)會(huì)造成很大誤差。最小范數(shù)法最典型的就是最小二乘法，它是一種全局算法。解纏結(jié)果中不存在不能進(jìn)行相位展開的區(qū)域，但由于此算法穿過(guò)而非繞過(guò)殘差點(diǎn)，很容易使相位誤差傳遞到整個(gè)區(qū)域。這就造成了整體上結(jié)果為連續(xù)的，但實(shí)際上沒(méi)有一點(diǎn)相位的值是準(zhǔn)確的[9]。網(wǎng)絡(luò)流法是通過(guò)全局范圍搜索最優(yōu)路徑和最短枝切線來(lái)求得最小化問(wèn)題的最優(yōu)解。在網(wǎng)絡(luò)流成熟的條件下網(wǎng)絡(luò)流算法可大大提高運(yùn)算的效率[10]，與此同時(shí)保證一定的精度，在誤差的限制方面有著良好的效果，可以防止誤差傳遞到整個(gè)區(qū)域內(nèi)，對(duì)于噪聲較為嚴(yán)重的干涉圖也可以進(jìn)行解纏并取得較為精確的結(jié)果[11-13]。但是在噪聲較大的環(huán)境下網(wǎng)絡(luò)流法的解纏結(jié)果會(huì)將噪聲不可避免地沿著積分路徑傳遞，造成解纏結(jié)果的嚴(yán)重“尖峰毛刺”現(xiàn)象，使得解纏結(jié)果不理想。

結(jié)合殘差點(diǎn)的思想、曲面擬合法以及網(wǎng)絡(luò)流法提出了一種新的相位解纏算法。新算法根據(jù)相關(guān)度找出解纏結(jié)果的噪聲點(diǎn)，然后運(yùn)用最小曲面擬合法重構(gòu)值代替噪聲點(diǎn)的數(shù)值，最后沿著積分路徑更正有誤差傳遞的結(jié)果。逐一處理所有噪聲點(diǎn)，最終的解纏結(jié)果可以看作是由多個(gè)最小曲面多次貼合的結(jié)果。本文首先介紹相關(guān)理論，然后介紹新算法的具體步驟和推導(dǎo)，最后通過(guò)實(shí)驗(yàn)驗(yàn)證可行性。

1 新算法描述

網(wǎng)絡(luò)流算法最早見(jiàn)于Costantini等提出的基于網(wǎng)絡(luò)流的相位解纏方法[8],從根本上來(lái)講此種方法仍屬于路徑跟蹤法。這種方法是將相位解纏問(wèn)題轉(zhuǎn)化為最小化問(wèn)題，通過(guò)在全局范圍內(nèi)搜索路徑和最短枝切線來(lái)求得最小化問(wèn)題的最優(yōu)解。在一個(gè)M×N大小的方格網(wǎng)內(nèi),設(shè)φ和φ分別表示解纏和未解纏的相位,則有

φ(i,j) =φ(i,j) + 2πn

(1)

式中，n為整數(shù)，φ(i,j)∈[-π，π]，相位解纏就是從φ(i,j)到φ(i,j)的過(guò)程。根據(jù)二維行為解纏原理，定義相鄰像素點(diǎn)間的差分估計(jì)[14]：

Δφ1(i,j)=φ(i+1,j)-φ(i,j)+2πn1

(2)

Δφ2(i,j)=φ(i,j+1)-φ(i,j)+2πn2

(3)

式中，n1(i,j)，n2(i,j) 為基于先驗(yàn)知識(shí)選取，使Δφl(shuí)(i,j)∈[-π,π)(l=1,2) 成立的整數(shù)值。由于積分路徑的不同，Δφl(shuí)(i,j)∈[-π,π)(l=1,2) 并不能和相鄰點(diǎn)的差分保持一致，因而定義以下差分殘差：

由于k1(i,j)，k2(i,j)都是非常小的數(shù)，可以用如下的最小化問(wèn)題來(lái)估計(jì)殘差：

如圖1所示，根據(jù)網(wǎng)絡(luò)流理論，式(5)這個(gè)最小化問(wèn)題可以轉(zhuǎn)化為最小費(fèi)用流來(lái)解決。最小費(fèi)用流問(wèn)題的輸入為各個(gè)節(jié)點(diǎn)的度，即各個(gè)殘差的值和每條流的費(fèi)用，式中c1(i,j)，c2(i,j)可以根據(jù)像素點(diǎn)的相關(guān)度得出，而該問(wèn)題的輸出為各條流的流量。式(5)的最小化問(wèn)題可轉(zhuǎn)化為如下公式：

圖1 網(wǎng)絡(luò)圖

(6)

min(0,k1(i,j))

(7)

min(0,k2(i,j))

(8)

在計(jì)算出k1(i,j)和k2(i,j)之后，再計(jì)算相位梯度，最終有

(9)

由后面的實(shí)驗(yàn)結(jié)果可以看出，通過(guò)網(wǎng)絡(luò)流方法的解纏結(jié)果和枝切法、最小二乘法相比有很大的優(yōu)勢(shì)，沒(méi)有“解纏孤島”的現(xiàn)象，并且解纏精度要比其他兩種方法高。但是在噪聲較大的環(huán)境下網(wǎng)絡(luò)流法的解纏結(jié)果會(huì)將噪聲不可避免地沿著積分路徑傳遞，造成解纏結(jié)果的嚴(yán)重“尖峰毛刺”現(xiàn)象，造成解纏結(jié)果的失真。很明顯這些隨機(jī)出現(xiàn)的“毛刺”點(diǎn)是由于噪聲的干擾引起的。可以設(shè)計(jì)方法識(shí)別出這些有誤差的點(diǎn)，然后根據(jù)這個(gè)點(diǎn)與周圍其余點(diǎn)的聯(lián)系重新構(gòu)造出更準(zhǔn)確的值。通過(guò)誤差點(diǎn)識(shí)別方法找出這些由于噪聲出現(xiàn)的“毛刺”點(diǎn)的集合，假設(shè)這些點(diǎn)(i,j)∈S1(i=0,1,…，n-1，j=0,1,…，m-1)，設(shè)集合S1的殘差值k1(i,j)，k2(i,j)：

i=0,1,…,n-1,j=0,1,…,m-1

(10)

i=0,1,…,n-1,j=0,1,…,m-1

(11)

誤差點(diǎn)的重構(gòu)值是根據(jù)其與周圍點(diǎn)相關(guān)程度和曲面趨勢(shì)進(jìn)行的重構(gòu)，經(jīng)過(guò)重構(gòu)后的相位值在一定程度上降低了噪聲的干擾，減少了噪聲的疊加,所以可以得出

從而得出經(jīng)過(guò)重構(gòu)后的網(wǎng)絡(luò)圖的目標(biāo)函數(shù)：

結(jié)合式(5)、式(12) 和式(13)可以得出

目標(biāo)函數(shù)最小化是相位解纏的本質(zhì)，以上的推導(dǎo)可以證明通過(guò)對(duì)網(wǎng)絡(luò)流算法解纏結(jié)果中的噪聲點(diǎn)進(jìn)行精準(zhǔn)識(shí)別、定位與重構(gòu)后，可以有效改善“毛刺”現(xiàn)象，并且能夠進(jìn)一步最小化目標(biāo)函數(shù)，從而提高解纏結(jié)果的準(zhǔn)確度。考慮到解纏算法的時(shí)間成本以及內(nèi)存要求，對(duì)于誤差點(diǎn)的識(shí)別和利用曲面擬合重構(gòu)數(shù)值分別設(shè)計(jì)了以下方法實(shí)現(xiàn)。

1.1 誤差點(diǎn)處理

受到枝切法中識(shí)別“殘差點(diǎn)”思想的啟發(fā)，按照如下方法對(duì)誤差點(diǎn)進(jìn)行識(shí)別并處理。根據(jù)二維解纏原理連續(xù)相位的假設(shè)，首先設(shè)置π為閾值，再依次檢驗(yàn)所有解纏后的像素點(diǎn)。在像素點(diǎn)周圍3×3的區(qū)域內(nèi)進(jìn)行檢測(cè)，如果像素點(diǎn)與其所有周圍的8個(gè)點(diǎn)的差值的絕對(duì)值都大于閾值，那么識(shí)別此像素點(diǎn)為誤差點(diǎn)。如圖2(a)所示，3×3區(qū)域的像素點(diǎn)中心點(diǎn)明顯比其周圍的8個(gè)像素點(diǎn)的值大很多。因?yàn)榻?jīng)過(guò)網(wǎng)絡(luò)流處理后的數(shù)據(jù)在大概率上可以代替真實(shí)相位值，那么基于連續(xù)相位的解纏假設(shè)，如果出現(xiàn)上述的像素點(diǎn)就可以認(rèn)定此點(diǎn)由于噪聲的干擾造成了數(shù)據(jù)的嚴(yán)重失真，可以將此點(diǎn)識(shí)別記錄，認(rèn)定為誤差點(diǎn)，稱之為E1誤差點(diǎn)。

圖2 E1,E2誤差點(diǎn)三維示意圖

式中，DI是像素點(diǎn)的相位值，而CI是根據(jù)像素點(diǎn)的相關(guān)度得出的權(quán)值。最終得到的D5就是利用最小曲面擬合構(gòu)造的新像素點(diǎn)的值。

如圖3(b)所示的12個(gè)像點(diǎn)，其中N5和N8為3×4區(qū)域識(shí)別的誤差點(diǎn)，而其余的10個(gè)點(diǎn)均是認(rèn)定的可靠數(shù)據(jù)。在這種情況下，將這個(gè)區(qū)域分割成兩個(gè)3×3的區(qū)域，如圖4所示。

圖3 E1,E2誤差點(diǎn)二維示意圖

圖4 E2誤差點(diǎn)分割示意圖

(16)

圖5 E3誤差點(diǎn)俯視圖

1.2 新算法步驟

根據(jù)以上描述的方法，實(shí)現(xiàn)的算法步驟如下：

步驟1：根據(jù)纏繞相位信息計(jì)算出像素點(diǎn)間相位的相關(guān)系數(shù)，利用原始相位信息以及各個(gè)像素點(diǎn)相關(guān)系數(shù)等信息，結(jié)合圖論中圖的概念構(gòu)造網(wǎng)絡(luò)圖；

步驟2：對(duì)初步構(gòu)造的圖，以圖論的最大流最小割定理(max flow/min cut theory)為基礎(chǔ)的增廣路徑算法求得最小費(fèi)用路徑；

步驟3：根據(jù)求得的路徑對(duì)相位進(jìn)行積分求得初步的解纏結(jié)果，并保存積分路徑；

步驟4：對(duì)初步解纏結(jié)果進(jìn)行第1輪誤差點(diǎn)識(shí)別，識(shí)別出所有E1點(diǎn)，使用最小曲面擬合方法進(jìn)行數(shù)據(jù)重構(gòu)，并利用重構(gòu)值沿著步驟3的積分路徑修正積分路徑的值；

步驟5：進(jìn)行第2輪誤差點(diǎn)識(shí)別，識(shí)別出所有E2點(diǎn)，使用最小曲面擬合方法進(jìn)行數(shù)據(jù)重構(gòu)，并利用重構(gòu)值沿著步驟3的積分路徑修正積分路徑的值；

步驟6：進(jìn)行第3輪誤差點(diǎn)識(shí)別，識(shí)別出所有E3點(diǎn)，使用最小曲面擬合方法進(jìn)行數(shù)據(jù)重構(gòu)，并利用重構(gòu)值沿著步驟3的積分路徑修正積分路徑的值；

步驟7：重復(fù)步驟6至少兩次，識(shí)別出由E4點(diǎn)轉(zhuǎn)而來(lái)的E3點(diǎn)，使用最小曲面擬合方法進(jìn)行數(shù)據(jù)重構(gòu)，并利用重構(gòu)值沿著步驟3的積分路徑修正積分路徑的值。

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

為了驗(yàn)證算法的有效性和正確性，分別用枝切法、等權(quán)最小二乘法、網(wǎng)絡(luò)流法以及改進(jìn)算法進(jìn)行實(shí)驗(yàn)，并將結(jié)果進(jìn)行對(duì)比。實(shí)驗(yàn)分為仿真數(shù)據(jù)實(shí)驗(yàn)和實(shí)測(cè)數(shù)據(jù)實(shí)驗(yàn)。模擬數(shù)據(jù)中加入不同程度的隨機(jī)噪聲和椒鹽噪聲，構(gòu)成兩組不同程度的信噪比進(jìn)行仿真實(shí)驗(yàn)，用以驗(yàn)證不同信噪比情況下方法的正確性和精確度[15]。評(píng)定的參數(shù)用運(yùn)行時(shí)間和均方根誤差(REMS)兩個(gè)參數(shù)來(lái)衡量，均方根誤差計(jì)算公式如式(17)所示。實(shí)驗(yàn)環(huán)境為Matlab R2004a，計(jì)算機(jī)運(yùn)算條件為Intel(R) Xeon(R) CPU E3-1241 v3+Windows 7專業(yè)版64 bit+8 G內(nèi)存。

2.1 模擬數(shù)據(jù)仿真

如圖6所示，模擬數(shù)據(jù)是Matlab軟件中peaks函數(shù)產(chǎn)生模擬相位曲面，大小為512像素×512像素。對(duì)模擬曲面加入隨機(jī)噪聲和椒鹽噪聲實(shí)現(xiàn)信噪比(SNR)為6.11 dB和-0.51 dB的模擬數(shù)據(jù)。圖7(a)是信噪比為6.11 dB仿真相位的曲面圖，圖7(b)是纏繞相位的俯視圖。將數(shù)據(jù)進(jìn)行纏繞后分別運(yùn)用枝切法、等權(quán)最小二乘法、網(wǎng)絡(luò)流法和新方法進(jìn)行解纏仿真。解纏結(jié)果如圖8所示，運(yùn)行結(jié)果如表1所示。

圖6 模擬數(shù)據(jù)示意圖

圖7 SNR=6.11 dB纏繞相位圖

算法運(yùn)行時(shí)間/s均方根誤差/rad枝切法25.891.002最小二乘法2.82361.4995網(wǎng)絡(luò)流法24.590.795新算法30.190.5733

圖9(a)是信噪比為-0.51 dB模擬數(shù)據(jù)的纏繞圖，圖9(b)是圖9(a)的俯視圖。將數(shù)據(jù)分別運(yùn)用枝切法、最小二乘法、網(wǎng)絡(luò)流法和新算法進(jìn)行解纏仿真。解纏結(jié)果如圖10所示，仿真運(yùn)行結(jié)果匯總?cè)绫?所示。

圖8 SNR=6.11 dB模擬數(shù)據(jù)實(shí)驗(yàn)結(jié)果

圖9 SNR=-0.51 dB纏繞相位圖

算法運(yùn)行時(shí)間/s均方根誤差/rad枝切法42.27632.4258最小二乘法2.5721.6207網(wǎng)絡(luò)流法28.47341.3214新算法35.03620.5861

圖10 SNR=-0.51 dB模擬數(shù)據(jù)實(shí)驗(yàn)結(jié)果

2.2 真實(shí)數(shù)據(jù)實(shí)驗(yàn)

圖11 實(shí)測(cè)數(shù)據(jù)

算法運(yùn)行時(shí)間/sσ枝切法127.6868—最小二乘法2.57403.1157網(wǎng)絡(luò)流法48.92193.0344新算法59.05882.6157

圖12 實(shí)測(cè)數(shù)據(jù)實(shí)驗(yàn)結(jié)果

3 結(jié)束語(yǔ)

本文的方法借鑒枝切法思想，采用網(wǎng)絡(luò)流法和曲面擬合的方法，對(duì)網(wǎng)絡(luò)流算法增加了最小曲面擬合的步驟。以最小的時(shí)間代價(jià)，利用網(wǎng)絡(luò)流解纏算法過(guò)程中的中間參數(shù)和初步解纏結(jié)果，對(duì)噪聲引起的錯(cuò)誤進(jìn)行了可靠重構(gòu)，并沿著解纏積分路徑更正了積分過(guò)程中誤差的傳遞。無(wú)論從模擬數(shù)據(jù)實(shí)驗(yàn)還是真實(shí)數(shù)據(jù)實(shí)驗(yàn)來(lái)看，新方法都大大地增加了相位解纏精度，并證明了方法的可靠性，相比其他方法具有一定的優(yōu)越性。