梁小星，謝先明，孫玉錚

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

0 引言

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

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

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

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

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

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

(1)

(2)

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

1.1 ML算法

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

(3)

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

(4)

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

1.2 MAP算法

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

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

(5)

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

(6)

(7)

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

(8)

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

1.3 TV+MAP算法

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

(9)

(10)

與高斯馬爾可夫隨機場模型不同的是，TV模型為非局部模型，選擇非局部先驗?zāi)芰勘苊饬藢植砍瑓?shù)矢量的估計，從而簡化模型，提升了算法的運算效率。但算法對地形的適應(yīng)性也因為沒有局部先驗而降低。

2 本文方法

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

2.1 基于Sobel算子的邊緣檢測

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

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

圖1 邊緣檢測示意圖

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

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

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

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

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

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

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

3 實驗驗證及分析

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

(13)

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

3.1 仿真實驗1

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

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

圖3 模擬城市地形高程反演實驗結(jié)果對比

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

表2 算法性能對比

3.2 仿真實驗2

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

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

圖4 Peaks地形高程反演實驗結(jié)果對比

表3 算法性能對比

3.3 仿真實驗3

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

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

圖5 Long’s Peak地形高程反演實驗結(jié)果對比

表4 算法性能對比

4 結(jié)束語

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