梁小星，謝先明，孫玉錚

(桂林電子科技大學(xué)，廣西 桂林 541004)

0 引言

干涉合成孔徑雷達(dá)(synthetic aperture radar interferometry，InSAR)可以高精度、高可靠性地獲取地表三維信息和高程變化信息，被廣泛應(yīng)用于海洋監(jiān)控、火山監(jiān)測(cè)、地震檢測(cè)和數(shù)字高程重建等領(lǐng)域。傳統(tǒng)的單通道InSAR高程反演技術(shù)受限于Itoh相位連續(xù)性假設(shè)，通常難以有效獲取不連續(xù)地形干涉圖的真實(shí)高程信息[1]。鑒于此，多種不受Itoh相位連續(xù)性條件限制的多基線InSAR高程反演技術(shù)被相繼提出。文獻(xiàn)[2-3]提出了基于中國(guó)余數(shù)定理(chinese remainder theorem，CRT)的多基線InSAR高程反演技術(shù)，用多幅基線滿(mǎn)足兩兩互質(zhì)的干涉圖構(gòu)建同余方程組，利用CRT求得唯一解，能有效克服相位模糊的問(wèn)題，可從理想狀態(tài)的無(wú)噪聲干涉圖中獲取精確的地形高程信息。Yuan等[4]提出了閉環(huán)魯棒的CRT多基線高程反演技術(shù)，提高了CRT多基線高程反演技術(shù)的噪聲魯棒性，可有效地從噪聲較小的多幅干涉圖獲取目標(biāo)高程信息。文獻(xiàn)[5-7]提出了多基線最大似然函數(shù)(maximum likelihood，ML)法，該方法利用真實(shí)相位與纏繞相位符合圓高斯分布的性質(zhì)，使用貝葉斯條件概率對(duì)真實(shí)相位進(jìn)行估計(jì)，在低噪聲情況下能得到較好的解纏結(jié)果，但在噪聲較大時(shí)，解纏結(jié)果出現(xiàn)大量“毛刺”點(diǎn)，算法魯棒性較差。袁志輝等[8]在ML的基礎(chǔ)上引入了壞點(diǎn)判斷和加權(quán)均值濾波，提高了算法精度。Ferraiuolo等[9]則提出了基于高斯馬爾可夫隨機(jī)場(chǎng)(Markov random fields，MRFs)的最大后驗(yàn)概率(maximum a posteriori，MAP)估計(jì)方法，利用馬爾可夫隨機(jī)場(chǎng)統(tǒng)計(jì)分布模型來(lái)描述高程的先驗(yàn)分布，將高程反演轉(zhuǎn)化為能量函數(shù)最小化的優(yōu)化問(wèn)題，再使用迭代條件模型進(jìn)行求解。基于馬爾可夫隨機(jī)場(chǎng)的MAP算法，需要估計(jì)復(fù)雜的超參數(shù)，耗時(shí)較長(zhǎng)。Ferraioli等[10]提出一種全變分(total variation，TV)模型的MAP算法(即TV+MAP算法)，避免了復(fù)雜的超參數(shù)估計(jì)，大大縮短了運(yùn)行時(shí)間。謝先明等[11]提出了一種融合殘差點(diǎn)信息的TV+MAP算法，通過(guò)殘差信息對(duì)高程圖中的壞點(diǎn)進(jìn)行適當(dāng)加權(quán)更新獲得最終估計(jì)高程，提高了TV+MAP算法的精度。針對(duì)能量函數(shù)最小化的研究，多種基于置信傳播[12]、順序樹(shù)重加權(quán)消息傳遞[13-14]，以及多標(biāo)簽圖割[15-17]的算法相繼被提出，并在文獻(xiàn)[18]中進(jìn)行了詳細(xì)的分析和比較。

與此同時(shí)，Yu等[19]提出了基于聚類(lèi)分析的多基線InSAR高程反演算法，能在信噪比較高的干涉圖中獲得較好的結(jié)果，且該算法效率較高，不利之處在于該算法聚類(lèi)效果易受噪聲影響；Liu等[20]在此基礎(chǔ)上提出了一種基于密度聚類(lèi)的相位解纏算法，提高了聚類(lèi)算法的抗噪性能。Li等[21]提出了子空間聯(lián)合處理算法，通過(guò)子空間投影將真實(shí)相位與噪聲分離；隨后，索志勇等[22]在此基礎(chǔ)上對(duì)其信號(hào)空間維度估計(jì)的準(zhǔn)確性和算法復(fù)雜度進(jìn)行改進(jìn)，有效提高了算法性能。文獻(xiàn)[23-24]成功地將深度學(xué)習(xí)框架應(yīng)用到InSAR高程重建中，將高程重建轉(zhuǎn)化成相位模糊數(shù)的多分類(lèi)問(wèn)題，利用卷積神經(jīng)網(wǎng)絡(luò)進(jìn)行高程重建。

在上述算法中，ML、MAP和TV+MAP算法是目前多基線InSAR高程反演算法中最常用的算法。ML算法沒(méi)有考慮鄰域間的約束關(guān)系，算法魯棒性較差；MAP算法雖然引入了鄰域約束關(guān)系，但需要迭代估計(jì)超參數(shù)，運(yùn)行效率低，且不連續(xù)區(qū)域鄰域約束的引入將造成嚴(yán)重的誤差傳遞效應(yīng)；TV+MAP算法能有效避免超參數(shù)估計(jì)，但同樣未考慮不連續(xù)邊界鄰域約束導(dǎo)致的誤差傳遞問(wèn)題。為了解決上述問(wèn)題，本文提出一種基于邊緣檢測(cè)與路徑跟蹤策略的多基線InSAR高程反演方法。該算法分為2個(gè)步驟：第1步先用最大似然估計(jì)獲取粗略的地形高程，再用Sobel算子對(duì)濾波后的粗略地形高程進(jìn)行邊緣檢測(cè)，獲得地形的不連續(xù)邊界；第2步首先建立優(yōu)化的多基線InSAR高程反演模型，隨后利用單通道InSAR相位解纏技術(shù)中的路徑跟蹤策略引導(dǎo)構(gòu)建的多基線InSAR高程反演模型，沿高質(zhì)量像元到低質(zhì)量像元的路徑重新反演目標(biāo)高程，并在高程反演過(guò)程中根據(jù)邊緣檢測(cè)結(jié)果判斷是否引入鄰域約束，既有利于提高算法在連續(xù)區(qū)域的抗噪性，又可避免鄰域約束在不連續(xù)區(qū)域引起的誤差傳遞現(xiàn)象，從而達(dá)到增強(qiáng)算法魯棒性的目標(biāo)。

1 多基線InSAR高程反演及相關(guān)算法

在多基線InSAR高程反演模型中，地形高程、真實(shí)相位與纏繞相位之間存在的關(guān)系如式(1)、式(2)所示。

φi=φi(x)+2n(x)·π

(1)

(2)

式中：φi代表第i幅干涉圖的真實(shí)相位；φi代表第i幅干涉圖的纏繞相位(-π≤φi<π)；n代表整數(shù)倍的纏繞數(shù)；x代表對(duì)應(yīng)像元的坐標(biāo)位置；h為對(duì)應(yīng)的地表高程值；λ為雷達(dá)工作波長(zhǎng)；R為主天線到目標(biāo)場(chǎng)景中心的距離；α為水平基線角；Bi為第i條基線的長(zhǎng)度；κi為第i條基線對(duì)應(yīng)的高度靈敏度。

1.1 ML算法

ML算法假設(shè)真實(shí)相位與纏繞相位符合圓高斯分布[25]，則纏繞相位的概率密度函數(shù)表達(dá)如式(3)所示。

(3)

式中：bi=cos (κih-φi)；γi為第i幅干涉圖的相干系數(shù)。ML算法利用L幅單視干涉圖的聯(lián)合概率密度似然函數(shù)直接重建高程(式(4))。

(4)

ML算法利用了多基線InSAR干涉相位相互獨(dú)立的特性，有效減少了似然函數(shù)最大概率值的周期出現(xiàn)。

1.2 MAP算法

MAP算法是在ML算法的基礎(chǔ)上，引入馬爾可夫隨機(jī)場(chǎng)，約束相鄰像素點(diǎn)的高程范圍，是一種精度較高的多基線高程重建算法。MAP算法模型表達(dá)如式(5)所示。

hMAP=argmax lnFML(φ|h)·g(h，σ)

(5)

根據(jù)Hammersley-Clifford理論，任何類(lèi)型的馬爾可夫隨機(jī)場(chǎng)都可以用Gibbs分布表達(dá),如式(6)所示。

(6)

(7)

式中：Nx表示像素點(diǎn)x鄰域的像元集合，一般取8連通域；σxy為超參數(shù)集合σ的元素，代表高程的局部變化特征，一般可采用期望最大化算法[26](expectation maximization，EM)對(duì)其進(jìn)行估計(jì)。結(jié)合式(5)，MAP算法經(jīng)化簡(jiǎn)后的表達(dá)如式(8)所示。

(8)

對(duì)于式(8)的求解，本文實(shí)驗(yàn)中均使用Pascazio提出的條件迭代模型算法進(jìn)行10次迭代求解。相對(duì)于ML算法，MAP算法重建高程的精度更高。由于式(7)中復(fù)雜的局部約束關(guān)系，使得MAP算法適用于各類(lèi)地形，但該方法需要重復(fù)迭代估計(jì)超參數(shù)，算法耗時(shí)較長(zhǎng)。

1.3 TV+MAP算法

TV模型由于其對(duì)不同背景信息強(qiáng)大的適應(yīng)能力，成為圖像處理中最有用的先驗(yàn)?zāi)Ｐ?。在SAR應(yīng)用中，TV模型主要用來(lái)做圖像復(fù)原。Ferraioli等將TV模型成功應(yīng)用到InSAR高程重建，提出了TV+MAP算法，其TV能量函數(shù)表達(dá)如式(9)所示。

(9)

(10)

與高斯馬爾可夫隨機(jī)場(chǎng)模型不同的是，TV模型為非局部模型，選擇非局部先驗(yàn)?zāi)芰勘苊饬藢?duì)局部超參數(shù)矢量的估計(jì)，從而簡(jiǎn)化模型，提升了算法的運(yùn)算效率。但算法對(duì)地形的適應(yīng)性也因?yàn)闆](méi)有局部先驗(yàn)而降低。

2 本文方法

本文方法分為2個(gè)步驟：第1步直接利用多基線最大似然估計(jì)算法從多幅不同基線的干涉相位圖中獲取粗略的地形高程，隨后利用窗口為3×3的中值濾波與均值濾波算法對(duì)ML算法獲取的地形高程進(jìn)行濾波，接著用Sobel算子對(duì)濾波后的地形高程圖進(jìn)行邊緣檢測(cè)，獲得地形的不連續(xù)邊界；第2步首先建立優(yōu)化的多基線InSAR高程反演模型，隨后利用傳統(tǒng)單通道InSAR相位解纏技術(shù)中的路徑跟蹤策略引導(dǎo)構(gòu)建的多基線InSAR高程反演模型沿高質(zhì)量像元到低質(zhì)量像元的路徑反演目標(biāo)高程，并在高程反演過(guò)程中根據(jù)邊緣檢測(cè)結(jié)果判斷是否引入鄰域約束，既有利于提高算法在連續(xù)區(qū)域的抗噪性，又可避免鄰域約束在不連續(xù)區(qū)域引起的誤差傳遞現(xiàn)象，從而達(dá)到增強(qiáng)算法魯棒性的目標(biāo)。

2.1 基于Sobel算子的邊緣檢測(cè)

鄰域約束是MAP以及TV+MAP算法的核心，有利于在干涉圖連續(xù)區(qū)域得到高精度的相位解纏結(jié)果。然而，在不連續(xù)區(qū)域鄰域約束的引入將造成嚴(yán)重的誤差傳遞效應(yīng)，導(dǎo)致MAP和TV+MAP算法魯棒性降低。如何準(zhǔn)確地在連續(xù)區(qū)域引入鄰域約束且在不連續(xù)區(qū)域阻斷鄰域約束是提升相位解纏算法精度的關(guān)鍵。為解決這個(gè)問(wèn)題，本文先用ML算法獲取粗略的地形高程圖，隨后利用窗口為3×3的中值濾波算法與均值濾波算法對(duì)ML算法獲取的地形高程圖進(jìn)行濾波，接著再用Sobel算子對(duì)濾波后的地形高程圖進(jìn)行邊緣檢測(cè)，能夠有效地檢測(cè)出地形的不連續(xù)邊界并得到標(biāo)記邊界像元的二值圖像，記為edge。

關(guān)于邊緣檢測(cè)中的閾值設(shè)置問(wèn)題，以圖1(a)中的美國(guó)Long’s Peak公園真實(shí)地形高程數(shù)據(jù)為例，在傳統(tǒng)單基線InSAR高程重建算法中，當(dāng)干涉圖相鄰像元之間高程超過(guò)其基線所對(duì)應(yīng)高度模糊數(shù)的二分之一時(shí)，則無(wú)法滿(mǎn)足相位連續(xù)性假設(shè)。以長(zhǎng)基線干涉圖對(duì)應(yīng)的高度模糊數(shù)的二分之一作為邊緣檢測(cè)閾值，可以檢測(cè)出更多邊緣細(xì)節(jié)，但這些細(xì)節(jié)邊緣并不是真正意義上的不連續(xù)邊界，如圖1(a)和圖1(b)白色矩形框所示。使用短基線所對(duì)應(yīng)高度模糊數(shù)的二分之一作為邊緣檢測(cè)閾值，可以更有效檢測(cè)出真正的不連續(xù)邊界，如圖1(c)所示。

圖1 邊緣檢測(cè)示意圖

2.2 優(yōu)化的高程反演模型

優(yōu)化的高程反演模型如式(11)、式(12)所示。

圖2以不連續(xù)像元可能出現(xiàn)的一種分布情況進(jìn)行一般性分析，通過(guò)在不連續(xù)邊界周?chē)赃m應(yīng)地調(diào)整鄰域約束強(qiáng)弱，能夠有效避免誤差傳播并提高算法的魯棒性。

圖2 中心像元鄰域不連續(xù)邊界示意圖

算法流程及步驟如下所示。

算法開(kāi)始步驟 1 (邊緣提取)步驟1.1:利用多基線ML算法從多幅不同基線的干涉相位圖中獲取粗略的地形高程。步驟1.2:利用窗口為3×3的中值濾波與均值濾波算法對(duì)多基線ML算法獲取的地形高程圖進(jìn)行濾波。步驟1.3:利用Sobel算子對(duì)濾波后的地形高程圖進(jìn)行邊緣檢測(cè),得到標(biāo)記了不連續(xù)邊界點(diǎn)的二值圖像edge。步驟2 (高程重建)步驟2.1:在多基線InSAR干涉圖相位質(zhì)量圖中選擇5～10個(gè)質(zhì)量較好的像元,利用ML算法獲取其高程估計(jì)值。步驟2.2:在已反演高程的像元中,將它們4鄰域中未解纏的像元按相位質(zhì)量從大到小排列,作為鄰接列表。步驟2.3:取出鄰接列表中質(zhì)量最好的一個(gè)像元x,根據(jù)邊緣檢測(cè)得到的二值圖像進(jìn)行判斷。若該像元為不連續(xù)邊界像元,則利用式(11)獲取高程估計(jì)值;若該像元為連續(xù)像元,則計(jì)算該像元8鄰域中非邊界已解纏像元的高程平均值hyunwrap,并利用式(12)獲取高程估計(jì)值,將像元x的4鄰域中未高程反演的像元嵌入鄰接列表。步驟2.4:判斷鄰接列表是否為空,非空則執(zhí)行步驟2.3;為空,則算法結(jié)束,得到估計(jì)高程。算法結(jié)束

在步驟2.1中，選取5～10個(gè)點(diǎn)利用多基線ML算法獲取其高程估計(jì)值作為初始參照點(diǎn)。理論上只需要選取一個(gè)質(zhì)量最好的點(diǎn)獲得其高程估計(jì)值作為初始參照點(diǎn)即可，但為了增加鄰接列中像素點(diǎn)數(shù)量，提高算法的穩(wěn)健性，建議選取5～10個(gè)點(diǎn)。

3 實(shí)驗(yàn)驗(yàn)證及分析

為驗(yàn)證算法的有效性，且更好地突出邊緣檢測(cè)和優(yōu)化TV高程模型結(jié)合路徑跟蹤策略在本文方法中的作用，本文針對(duì)封閉式線狀突變地形、連續(xù)地形以及開(kāi)放型線狀突變地形數(shù)據(jù)進(jìn)行實(shí)驗(yàn)，并在仿真實(shí)驗(yàn)中將本文方法與多基線InSAR高程重建技術(shù)中經(jīng)典的ML、MAP以及TV+MAP算法進(jìn)行對(duì)比，檢驗(yàn)本文方法性能。本文采用歸一化均方根誤差計(jì)算實(shí)驗(yàn)結(jié)果誤差，計(jì)算方法如式(13)所示。

(13)

式中：hunwrap為估計(jì)高程；h為真實(shí)高程值；M和N分別為方位向和距離向的像元個(gè)數(shù)；p和q為像元二維坐標(biāo)。

3.1 仿真實(shí)驗(yàn)1

在實(shí)驗(yàn)1中采用如圖3(a)所示大小150像素×150像素的模擬城市地形。該模擬城市地形圖中存在多處突變邊緣，且相鄰點(diǎn)之間的最大高程差為150 m。InSAR系統(tǒng)參數(shù)為表1所示，3條基線分別為202.5 m、179.4 m、143.6 m。

圖3(c)是信噪比為7.4 dB的長(zhǎng)基線干涉圖。圖3(d)、圖3(e)分別為ML估計(jì)高程圖、經(jīng)過(guò)中值濾波和均值濾波的ML估計(jì)高程圖。圖3(f)為Sobel算子對(duì)圖3(d)進(jìn)行邊緣檢測(cè)后得到的標(biāo)記了不連續(xù)邊界的二值圖像。圖3(g)至圖3(i)分別為MAP 10次迭代估計(jì)的結(jié)果、TV+MAP估計(jì)結(jié)果和本文方法估計(jì)結(jié)果。圖3(j)至圖3(l)分別為MAP、TV+MAP和本文方法對(duì)應(yīng)的估計(jì)高程誤差。從圖3(d)至圖3(f)可以看出，ML估計(jì)高程雖然存在大量的毛刺點(diǎn)，但其地形的邊緣特性依然清晰可見(jiàn)，經(jīng)過(guò)中值和均值濾波處理后，能得到保持良好邊緣特性的高程圖，使用Sobel算子進(jìn)行邊緣檢測(cè)能夠有效提取出地形的不連續(xù)邊界點(diǎn)。表2給出了各算法的高程重建的誤差，以及針對(duì)突變區(qū)域的高程重建誤差。結(jié)合表2以及圖3(g)至圖3(l)可以看出，本文方法針對(duì)該地形的高程重建效果與TV+MAP算法效果相當(dāng)，都能夠較好地估計(jì)出地形高程，且誤差主要集中在不連續(xù)的邊界處。由于該地形高程值較為單一，因此本文方法相較于TV+MAP算法精度提高并不明顯。

圖3 模擬城市地形高程反演實(shí)驗(yàn)結(jié)果對(duì)比

表1 多基線InSAR系統(tǒng)基本參數(shù)

表2 算法性能對(duì)比

3.2 仿真實(shí)驗(yàn)2

為了進(jìn)一步驗(yàn)證本文方法在連續(xù)區(qū)域的處理效果，在實(shí)驗(yàn)2中，采用如圖4(a)所示大小為256像素×256像素的連續(xù)地形Peaks進(jìn)行實(shí)驗(yàn)。圖4(b)為Peaks地形的二維高程圖，圖4(c)是含噪聲的長(zhǎng)基線干涉圖，模擬InSAR系統(tǒng)參數(shù)與實(shí)驗(yàn)1相同，3條基線分別為415.5 m、272.2 m、151.8 m。

圖4(d)、圖4(e)分別為ML估計(jì)高程圖、經(jīng)過(guò)中值濾波和均值濾波的ML估計(jì)高程圖。圖4(f)為Sobel算子對(duì)圖4(e)進(jìn)行邊緣檢測(cè)后得到的不連續(xù)邊界標(biāo)記圖。從圖中可以看出，ML估計(jì)高程布滿(mǎn)了密集的毛刺點(diǎn)，經(jīng)過(guò)濾波處理后仍然存在部分毛刺點(diǎn)，由于地形本身是連續(xù)的，這些濾波處理后的毛刺點(diǎn)在邊緣檢測(cè)中被成功過(guò)濾。從圖4(g)至圖4(i)和表3的結(jié)果可以看出，在連續(xù)區(qū)域中，本文方法中優(yōu)化的TV高程反演模型在路徑跟蹤策略的引導(dǎo)下進(jìn)行遞推估計(jì),對(duì)算法的高程重建精度有很大提升。

圖4 Peaks地形高程反演實(shí)驗(yàn)結(jié)果對(duì)比

表3 算法性能對(duì)比

3.3 仿真實(shí)驗(yàn)3

為了進(jìn)一步綜合性地驗(yàn)證本文算法在連續(xù)區(qū)域及不連續(xù)區(qū)域的性能，在實(shí)驗(yàn)3中，采用如圖5(a)所示，大小為458像素×152像素的美國(guó)國(guó)家公園Long’s Peak地形的真實(shí)數(shù)字高程，該地形右側(cè)有一條垂直斷裂帶將地形分割成2塊連續(xù)的區(qū)域，高程范圍為0～127 m，最大高程差為127 m。圖5(c)是長(zhǎng)基線干涉圖，模擬InSAR系統(tǒng)參數(shù)與實(shí)驗(yàn)1相同，3條基線分別為415.5 m、272.2 m、151.8 m。

圖5(d)、圖5(e)分別為ML估計(jì)高程圖、經(jīng)過(guò)中值濾波和均值濾波的ML估計(jì)高程圖。圖5(e)為Sobel算子對(duì)圖5(d)進(jìn)行邊緣檢測(cè)后得到的不連續(xù)邊界標(biāo)記圖。圖5(g)至圖5(i)分別為MAP 10次迭代估計(jì)的結(jié)果、TV+MAP估計(jì)結(jié)果和本文方法估計(jì)結(jié)果，從圖中可以看出，本文方法高程重建效果最好。圖5(j)至圖5(m)分別為ML估計(jì)、MAP 10次迭代估計(jì)、TV+MAP估計(jì)以及本文方法估計(jì)高程與真實(shí)高程對(duì)比的誤差圖，結(jié)合表4可以看出，本文方法有效提取出了地形不連續(xù)邊界，通過(guò)阻斷不連續(xù)邊界處的鄰域約束，避免誤差傳遞，在連續(xù)區(qū)域引入鄰域約束的同時(shí)，結(jié)合路徑跟蹤策略進(jìn)行遞推估計(jì)，有效地提高了算法的魯棒性，且算法的運(yùn)行效率也與TV+MAP模型相當(dāng)。

圖5 Long’s Peak地形高程反演實(shí)驗(yàn)結(jié)果對(duì)比

表4 算法性能對(duì)比

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

論文提出了一種基于邊緣檢測(cè)與路徑跟蹤策略的多基線InSAR高程反演算法。通過(guò)提取地形的不連續(xù)邊界，在連續(xù)區(qū)域引入鄰域約束并結(jié)合路徑跟蹤策略進(jìn)行高程重建，有效地提高了算法在連續(xù)區(qū)域的魯棒性；在不連續(xù)邊界處阻斷鄰域約束，避免了誤差傳遞現(xiàn)象。通過(guò)3組仿真實(shí)驗(yàn)證明本文方法對(duì)于不同的地形進(jìn)行高程重建時(shí)都具有較高的精度、效率和穩(wěn)健性。