基于聚類分析多通道InSAR聯(lián)合相位解纏算法?

薛永嬌,禹衛(wèi)東

(1.中國科學(xué)院電子學(xué)研究所,北京100190;2.中國科學(xué)院大學(xué),北京100039)

0 引言

干涉合成孔徑雷達(Interferometric Synthetic Aperture Radar,InSAR)測量技術(shù)是通過從兩幅SAR復(fù)數(shù)據(jù)圖像中提取的干涉相位差,利用地表高程信息和干涉相位差之間的關(guān)系,獲取地表三維信息的重要遙感技術(shù)[1-2]。傳統(tǒng)的單通道InSAR技術(shù)相位解纏算法通常要求相鄰像素點的相位變化不超過π,即要求高程重建區(qū)域具有空間連續(xù)性,在地形起伏較大的地區(qū),例如高山、城市等區(qū)域,傳統(tǒng)單通道InSAR無法有效解決復(fù)雜地形高程反演問題[2],并且單通道相位解纏方法運算較為復(fù)雜,受噪聲影響比較大,存在誤差傳遞現(xiàn)象,因此很難滿足實際反演精度的要求。

隨著InSAR技術(shù)的發(fā)展,多通道InSAR測量技術(shù)(包括多頻率InSAR測量技術(shù)[3]和多基線In-SAR測量技術(shù)[4])成為目前InSAR技術(shù)的研究熱點。多頻率InSAR測量技術(shù)通過使用不同的工作頻率或者對SAR信號頻譜進行子帶分割得到[5],通過載頻不同,擴大干涉處理高度模糊數(shù)。多基線InSAR測量技術(shù)通過在單航過飛行器上設(shè)置多個天線或者對同一區(qū)域多次航過,得到同一地區(qū)不同角度或者不同時間的SAR圖像,通過基線不同,擴大干涉處理高度模糊數(shù)[6]。模糊高度數(shù)的擴大有效解決了高程層疊問題和相位欠采樣問題,實現(xiàn)對復(fù)雜地形進行相位解纏,得到精確度更高的解纏相位信息[7]。

隨著多通道干涉SAR系統(tǒng)的發(fā)展,多通道高度重建算法也不斷地涌現(xiàn)出來。自從多通道InSAR高程重建算法提出以來,各類方法中最具代表性的算法有:1)基于多通道數(shù)據(jù)的相位概率統(tǒng)計特性直接進行高程重建,代表算法有多通道最大似然估計(ML)算法[8]和基于Markov隨機模型的MRF-MAP算法[9]。ML算法反演精度與工作頻率、相干系數(shù)以及干涉圖數(shù)目有很大關(guān)系,MRF-MAP算法能夠通過利用相鄰像素信息提高反演精度,但運算時間長、內(nèi)存壓力大。2)利用多通道InSAR數(shù)據(jù)的多樣性將單通道相位展開算法拓展為多通道相位展開算法,最具代表性的是最小平方法[10]。3)利用不同通道的纏繞相位和對應(yīng)高度模糊數(shù)之間的線性關(guān)系進行高程反演。最早在1994年由Xu提出,同時提出中國余數(shù)定理(CRT)的多通道相位解纏方法[11],但是該方法存在魯棒性較差的問題,對噪聲非常敏感,所以在實際應(yīng)用中比較差。2011年,于瀚雯在CRT算法的基礎(chǔ)上提出了一種快速聚類的高程重建(CA)方法[6],將具有相同模糊數(shù)向量的像素聚為一類,對某一類像素進行相位解纏繞,從而達到降噪和快速相位解纏的目的,然而數(shù)據(jù)量過大時,由于噪聲的積累,像素聚類數(shù)目選擇對高程反演精度產(chǎn)生很大影響。在CA算法基礎(chǔ)上,劉會濤等在2015年提出了密度聚類的改進聚類分析高度重建(CANOPUS)算法[12],利用多通道InSAR纏繞相位的線性組合,作為待聚類數(shù)據(jù),充分利用待聚類數(shù)據(jù)的空變特性,使得聚類精度有所提高、速度加快。

本文針對CA算法噪聲魯棒性差、反演精度較低的問題提出了基于聚類分析的多通道InSAR聯(lián)合高程重建算法(CARE)。第1節(jié)首先介紹基于聚類分析的多通道高程重建算法原理以及該算法存在的問題,第2節(jié)詳細介紹聚類分析和區(qū)域擴展的聯(lián)合相位解纏算法原理,第3節(jié)給出聯(lián)合高程重建算法在仿真和實測數(shù)據(jù)集上的重建結(jié)果,第4節(jié)給出分析和結(jié)論。

1 基于聚類分析的InSAR相位解纏算法原理

在多通道的干涉合成孔徑雷達系統(tǒng)中,φl為第l通道的纏繞相位值(-π≤φl＜π),?l為第l通道的實際相位值,纏繞相位和實際相位之間滿足關(guān)系:

以雙頻InSAR系統(tǒng)為例,實際的地形高度h和絕對相位?l之間滿足關(guān)系式:

式中,λl為第l通道的波長,B⊥為垂直基線長度,r為主天線到目標的斜距,θ為主天線到目標的斜視角。對于同一觀測區(qū)域,由式(1)和式(2),可得到雙頻InSAR系統(tǒng)中纏繞相位和頻率的約束關(guān)系:

由此,可以得到不同通道波長、模糊數(shù)向量和纏繞相位之間的關(guān)系為

式中,k,λ和φ分別表示相應(yīng)通道對應(yīng)像素模糊數(shù)向量、波長和纏繞相位,下標為通道數(shù)。由式(4)可知,干涉圖中某一像素的信息可以表示成線性方程,其斜率和截距由λ1,λ2,φ1和φ2共同確定。該像素的模糊數(shù)向量[k1,k2]是該直線上的整數(shù)點。因此,具有相同模糊數(shù)向量的所有像素對應(yīng)同一直線。在低噪情況下,截距存在明顯的聚類現(xiàn)象,如圖1所示,圖1(a)是待聚類數(shù)據(jù)。

聚類分析算法利用相同模糊數(shù)向量的像素具有式(4)中相同截距的特性,將具有相同模糊數(shù)向量的像素聚為一類。其算法流程為:

1)獲得截距直方圖和直方圖包絡(luò);

2)獲得直方圖包絡(luò)峰值,將截距峰值內(nèi)所有像素歸為一類;

3)直方圖中峰值點的橫坐標作為該類中截距,通過搜索方法獲得滿足式(4)的[k1,k2],這組[k1,k2]即為該峰值截距對應(yīng)的模糊數(shù)向量,通過式(1)即可獲得解纏相位值。

在InSAR圖像處理中,數(shù)據(jù)量大,很容易因為噪聲的積累淹沒某個聚類中心,導(dǎo)致該類消失,噪聲越高,聚類的準確性越不可靠。為了減少噪聲積累,本文提出只針對突變地形進行基于聚類分析的多通道相位解纏。

圖1 截距直方圖包絡(luò)

2 聯(lián)合相位解纏算法基本原理

在InSAR相位圖中,在地形連續(xù)的高相干區(qū)域,相鄰像素之間的相位變化較小,在隨機散射嚴重或者地形起伏較大的區(qū)域,相鄰像素之間的相位變化較大,相位產(chǎn)生突變[15]。因為多通道InSAR系統(tǒng)對突變地形解纏精度高,所以本算法通過邊緣提取,獲得突變地形位置信息。

2.1 邊緣提取

干涉相位圖中相位跳變明顯的區(qū)域通過圖像邊緣檢測可以得到十分明顯的邊緣曲線[14]。在一幅M×N的干涉圖中,利用式(4)中截距的空變特性信息,將截距矩陣設(shè)為邊緣提取矩陣:

為了找到邊緣,可以在矩陣ψ中利用一階或者二階導(dǎo)數(shù)來獲取相位變化,從而提取干涉相位圖的邊緣曲線:

式中,ψm,n∈N(ψi,j),N(ψi,j)為像素點ψi,j的鄰域(如4鄰域,8鄰域)。在InSAR圖像中,因為地形復(fù)雜導(dǎo)致相位變化不規(guī)則,因此,本文利用prewitt算子進行X方向和Y方向邊緣提取,利用sobel算子進行對角線方向的邊緣提取。prewitt表達式和sobel表達式如下所示:

邊緣提取后得到的曲線圖E′=[e i,j],i∈[0,M-1],j∈[0,N-1],該曲線圖對研究突變相位關(guān)系時,并不能提供有效信息,因此需要對邊緣曲線圖進行形態(tài)骨架提取。處理后可表示為

式中,skel[]為骨架提取算子,E″為處理后的邊緣曲線矩陣。

實際提取得到的并非是某一區(qū)域的邊緣,而是邊緣的輪廓,并不是有效的突變位置信息,所以對檢測到的邊緣進行二次邊緣檢測,得到相鄰邊緣E。E是一個二值矩陣,有效記錄了相位突變的位置信息。二次邊緣檢測也保證了待解纏像素24鄰域中有效已解纏像素的個數(shù),提高鄰域像素對待解纏像素的影響力,加快區(qū)域擴展速度。

為了保證聚類的準確性,本文提出對截距矩陣進行均勻采樣,提取的突變數(shù)據(jù)作為初步待聚類數(shù)據(jù)。為了進一步降低噪聲對聚類準確性的影響,本文將待聚類數(shù)據(jù)中低相干區(qū)域數(shù)據(jù)剔除,作為最終的待聚類數(shù)據(jù)集D。然后根據(jù)多通道CA相位解纏算法進行相位解纏,得到區(qū)域擴展種子像素集S。這樣在保證數(shù)據(jù)多樣性的前提下降低數(shù)據(jù)量,從而降低噪聲積累,提高聚類準確性。

2.2 區(qū)域擴展相位解纏

區(qū)域擴展處理是從種子點開始將其鄰域中滿足條件的像素合并到同一區(qū)域,直到處理完所有相鄰像素。InSAR干涉圖像中的相位解纏繞過程實際是將像素從“纏繞相位”擴展為“解纏相位”的過程[15]。在CA算法中,即使不會因噪聲積累淹沒聚類中心,含噪點的解纏相位仍然誤差較大,因此本文提出利用數(shù)據(jù)集D的解纏相位S作為區(qū)域擴展種子點。利用區(qū)域擴展進行其余像素的相位解纏,按照鄰域已經(jīng)解纏像素的“影響力”大小,選擇和估算待解纏像素的相位。

鄰域已經(jīng)解纏像素的“影響力”分為兩種:像素距離“影響力”和相位連續(xù)性“影響力”。將待解纏像素的24鄰域看作5階矩陣P5×5,8鄰域視為3階矩陣P3×3。P中已解纏像素所在位置值為1,待解纏像素所在位置值為0。設(shè)相位梯度矩陣為G3×3,以目標待解纏像素為中心,兩個像素的半徑中,若某一方向兩個像素均解纏,相位梯度矩陣G3×3在該方向的值為1,否則為0。圖2中,黑色像素為目標像素,白色像素為未解纏像素,網(wǎng)格和斜線分別表示8鄰域已解纏像素和16鄰域已解纏像素。將目標像素受到24鄰域中已經(jīng)解纏像素的距離“影響力”和相位連續(xù)“影響力”如下:

式中,ω1,ω2,ω3為歸一化影響權(quán)重系數(shù),且ω1＜ω2＜ω3。當(dāng)目標像素鄰域影響力大于某一閾值,即可對該像素解纏相位進行初步估計。設(shè)目標像素24鄰域中16鄰域和8鄰域已經(jīng)解纏像素的相位矩陣為H5×5和H3×3,則相位的初步估算值為

式中,N1,N2,N3分別為矩陣P5×5,P3×3和G3×3中非零值個數(shù)。得到的估算解纏相位需要通過可靠性檢測。設(shè)估算解纏相位預(yù)測可靠性參數(shù)為d p,是初步估算解纏相位和鄰域中解纏相位差值的加權(quán)平均值,其表達式如下:

當(dāng)d p＜π/2時,初步估算模糊數(shù)為

估算解纏相位:

這樣,對于一些含噪聲較大的像素點,估算相位不能通過可靠性檢測,可以忽略該像素,繼續(xù)進行區(qū)域擴展,待其周圍解纏相位個數(shù)增多,估算相位可靠性增大,同時,適當(dāng)放寬估算可靠性檢測條件為d p＜π或者d p＜2π,直到所有像素解纏完成。

圖2 待解纏像素鄰域的“影響力”模型

3 實驗結(jié)果

本文提出的基于聚類分析的多頻率聯(lián)合高程重建算法,通過對仿真數(shù)據(jù)(128像素×128像素)、基于美國Isolation Peak國家公園的仿真數(shù)據(jù)(458像素×157像素)和巴塞羅那部分城區(qū)多基線InSAR實測數(shù)據(jù)(96像素×131像素)進行驗證。

3.1 仿真數(shù)據(jù)實驗結(jié)果

仿真數(shù)據(jù)建立了一個128行×128列的模擬突變模型。如圖3所示,數(shù)據(jù)表征最大高度為450 m,最大高度差為450 m。高度模糊數(shù)分別為42.22 m和57.47 m,相鄰兩點的最大干涉相位差分別為21.32π和15.66π。表1是仿真數(shù)據(jù)的主要參數(shù)。

表1 基于模擬地形的仿真系統(tǒng)參數(shù)

圖3(a)為無噪聲仿真數(shù)據(jù)的三維立體圖,圖3(b)和圖3(c)為通道1和通道2對應(yīng)的含噪聲去平地干涉相位圖,圖像的信噪比為20 dB,圖3(d)為截距圖,即為CA算法的待聚類數(shù)據(jù),圖3(e)為通過對待聚類數(shù)據(jù)進行邊緣檢測得到的突變地形的位置信息,圖3(f)為通過采樣和二次邊緣檢測后得到的本文提出算法的待聚類數(shù)據(jù)。

圖3 仿真數(shù)據(jù)

下面,通過仿真數(shù)據(jù)來比較最大似然估計(MLE)相位解纏算法、CA算法及本文算法的性能,如圖4所示。

圖4 仿真數(shù)據(jù)高程重建結(jié)果

從圖4(a)可以看出,最大似然估計(ML)相位解纏結(jié)果中存在大量噪點。由圖4(b)可以看出,傳統(tǒng)CA算法解纏結(jié)果較好,但噪點數(shù)目仍然不少。圖4(c)本文算法相位解纏結(jié)果,圖4(d)為本文算法得到的高程反演誤差,可以看出該算法不僅解纏精度高,而且噪點數(shù)目少。算法性能對比如表2所示。

表2 基于模擬地形的算法性能對比

第二個仿真實驗采用的是美國Isolation Peak國家公園DEM數(shù)據(jù)。基于地形的多通道仿真系統(tǒng)頻率分別為4.9 GHz和3.6 GHz,干涉數(shù)據(jù)的相干系數(shù)約為0.85,其他參數(shù)如表3所示。圖5(a)和圖5(b)分別是高頻和低頻數(shù)據(jù)獲得的干涉相位圖,圖5(c)是截距圖,圖5(d)是由截距圖進行邊緣檢測得到的突變地形位置數(shù)據(jù),圖5(e)是經(jīng)過數(shù)據(jù)均勻采樣和突變檢測后得到的待聚類數(shù)據(jù)。

表3 基于真實地形的仿真系統(tǒng)參數(shù)

圖6(a)顯示了Isolation Peak地形的真實DEM值,圖6(b)是最大似然估計算法得到的DEM數(shù)據(jù),圖6(c)是聚類分析算法得到的DEM數(shù)據(jù),圖6(d)是本文算法得到的DEM數(shù)據(jù),圖6(e)是本文算法得到的高程反演誤差值。對比圖6(c)和圖6(d)可以發(fā)現(xiàn),本文算法解纏結(jié)果比CA算法解纏結(jié)果中噪點明顯減少,并且本文算法的誤差值要明顯小于CA算法。算法性能對比如表4所示。

表4 基于真實地形的算法性能對比

由表2和表4中不同算法相位解纏精度和運算時間對比發(fā)現(xiàn),本文提出算法較CA算法雖然運行時間相比于CA算法要慢,但是解纏精度提高10%,本文算法運行時間和最大似然估計算法時間相當(dāng),在提高解纏精度的同時,一定程度上保證了運行效率。

圖5 Isolation Peak仿真數(shù)據(jù)

圖6 Isolation Peak數(shù)據(jù)集高程重建結(jié)果

3.2 實測數(shù)據(jù)實驗

實測數(shù)據(jù)使用的是巴塞羅那部分城區(qū)干涉SAR圖像。這兩幅圖像由TerraSAR-X以重復(fù)軌道方式獲得。兩幅干涉圖對應(yīng)的基線長度分別為151.5 m和86.6 m,高度模糊數(shù)為36.64 m和64.12 m。圖7(a)、圖7(b)分別為基線151.5 m和86.6 m對應(yīng)干涉圖。圖8(a)為最大似然估計高程重建結(jié)果,圖8(b)為CA算法高程重建結(jié)果,圖8(c)為本文提出算法重建結(jié)果,圖8(d)為本文算法重建后的三維高程圖。

圖7 不同基線對應(yīng)Barcelona部分城市區(qū)域干涉圖

圖8 Barcelona實測數(shù)據(jù)集高程重建結(jié)果

4 結(jié)束語

本文提出一種聯(lián)合聚類分析和區(qū)域擴展的多通道干涉SAR高程重建算法,聚類分析算法對突變地形解纏精度高,區(qū)域擴展相位解纏算法對連續(xù)地形細節(jié)恢復(fù)效果好,本文算法結(jié)合兩種算法的優(yōu)點,與傳統(tǒng)CA算法相比,通過邊緣提取和對數(shù)據(jù)的均勻采樣降低待聚類數(shù)據(jù)量,減小數(shù)據(jù)噪聲積累。同時,均勻?qū)Ь垲悢?shù)據(jù)進行采樣,避免了區(qū)域擴展算法中噪聲的傳遞,噪聲魯棒性好、相位解纏精度高、計算量適中。該算法適合數(shù)據(jù)量大、解纏精度要求高的干涉SAR圖像相位解纏。但本文算法仍然存在不足,以后工作中,仍然需要對邊緣檢測和采樣獲得的待聚類數(shù)據(jù)進行優(yōu)化,減少區(qū)域擴展種子數(shù)據(jù)帶來的重建誤差。

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

[2]YU H,XING M,BAO Z.A Fast Phase Unwrapping Method for Large-Scale Interferograms[J].IEEE Trans on Geoscience and Remote Sensing,2013,51(7):4240-4248.

[3]PASCAZIO V,SCHIRINZI G.Estimation of Terrain Elevation by Multifrequency Interferometric Wide Band SAR Data[J].IEEE Signal Processing Letters,2001,8(1):7-9.

[4]FERRETTI A,PRATI C,ROCCA F.Multibaseline InSAR DEM Reconstruction:the Wavelet Approach[J].IEEE Trans on Geoscience and Remote Sensing,1999,37(2):705-715.

[5]王偉.多頻率InSAR高程信息重建算法研究[D].哈爾濱:哈爾濱工業(yè)大學(xué),2006:18-19.

[6]于瀚雯.單/多基線相位解纏繞技術(shù)研究[D].西安:西安電子科技大學(xué),2012:75-87.

[7]YU H,LI Z,BAO Z.A Cluster-Analysis-Based Efficient Multibaseline Phase-Unwrapping Algorithm[J].IEEE Trans on Geoscience and Remote Sensing,2011,49(1):478-487.

[8]PASCAZIO V,SCHIRINZI G.Multifrequency In-SAR Height Reconstruction Through Maximum Likelihood Estimation of Local Planes Parameters[J].IEEE Trans on Image Processing,2002,11(12):1478-1489.

[9]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.

[10]靳國旺,張紅敏,徐青,等.多波段InSAR差分濾波相位解纏方法[J].測繪學(xué)報,2012,41(3):434-440.

[11]YUAN Z,LI F,DENG Y,et al.Multichannel In-SAR DEM Reconstruction Through Closed-Form Robust Chinese Remainder Theorem[J].IEEE Geoscience and Remote Sensing Letters,2013,10(6):1314-1318.

[12]劉會濤,張歡,邢孟道,等.改進的多基線相位解纏繞CANOPUS算法[J].系統(tǒng)工程與電子技術(shù),2015,37(8):1767-1772.

[13]FANG S,MENG L,WANG L,et al.Quality-Guided Phase Unwrapping Algorithm Based on Reliability Evaluation[J].Applied Optics,2011,50(28):5446.

[14]郭春生.優(yōu)化的區(qū)域增長InSAR相位解纏算法[J].中國圖象圖形學(xué)報,2006,11(10):1380-1386.

[15]畢海霞,魏志強.基于區(qū)域識別和區(qū)域擴展的相位解纏算法[J].電波科學(xué)學(xué)報,2015,30(2):244-251.