郝華東，劉國林，陳賢雷，曹振坦

(1.舟山市質(zhì)量技術(shù)監(jiān)督檢測研究院計(jì)量檢定測試中心，舟山 316021;2.山東科技大學(xué)測繪科學(xué)與工程學(xué)院，青島 266510)

0 引言

InSAR是一種極具潛力的微波遙感技術(shù)，在地形測繪、地震形變、火山運(yùn)動(dòng)、地面沉降以及滑坡監(jiān)測等方面都具有廣泛應(yīng)用。然而，在星載或機(jī)載InSAR遙感圖像上，由于雷達(dá)熱噪聲、陰影、幾何重疊、時(shí)間去相關(guān)和相位失真等因素引起的噪聲和偽信號往往導(dǎo)致干涉相位的不連續(xù)，進(jìn)而引起局部相位誤差[1]，這使得相位解纏成為國內(nèi)外相關(guān)學(xué)者的研究熱點(diǎn)。相位解纏算法主要有路徑跟蹤算法[2-4]、最小范數(shù)算法[5-7]、基于最優(yōu)估計(jì)方法[8-13]以及基于區(qū)域分割和編碼的相位解纏算法[14]等?？柭鼮V波相位解纏算法作為一種基于最優(yōu)估計(jì)的方法，通過將相位的解纏問題轉(zhuǎn)化為相位的狀態(tài)估計(jì)問題，運(yùn)用卡爾曼濾波技術(shù)，最終實(shí)現(xiàn)相位解纏最優(yōu)解的求取。其優(yōu)點(diǎn)是直接對SAR復(fù)干涉圖像進(jìn)行卡爾曼濾波處理，克服了傳統(tǒng)相位解纏算法在解纏之前需要進(jìn)行預(yù)濾波處理的弊端，能夠同時(shí)實(shí)現(xiàn)濾波和解纏，簡化了數(shù)據(jù)處理過程。

相位解纏的卡爾曼濾波狀態(tài)模型一般假設(shè)相位的狀態(tài)變化是一個(gè)緩慢的過程。在地形平坦或起伏不大時(shí)，常規(guī)卡爾曼濾波相位解纏算法通常能夠得到可靠的解纏結(jié)果，但在地形陡峭或起伏較大時(shí)，則因相位的狀態(tài)變化較快，導(dǎo)致解纏結(jié)果存在誤差傳遞[12]。因此，地形起伏地區(qū)的 SAR干涉圖解纏需要充分考慮地形因素的影響，通過引入相關(guān)的地形信息來提高解纏的精度。文獻(xiàn)[12]提出了一種顧及地形因素的卡爾曼濾波相位解纏算法，該算法采用線性調(diào)頻Z變換(chirp z transform，CZT)進(jìn)行局部頻率估計(jì)來實(shí)現(xiàn)對地形的自適應(yīng)模擬，但該方法的精度不高且計(jì)算效率偏低。為此，本文提出了一種基于SRTM DEM的卡爾曼濾波相位解纏方法，并通過采用InSAR實(shí)際圖像數(shù)據(jù)對該算法的可行性和有效性進(jìn)行了驗(yàn)證。

1 基于地形的卡爾曼濾波狀態(tài)空間模型

卡爾曼濾波的觀測方程以及相關(guān)參數(shù)的定義詳見文獻(xiàn)[12]。這里只介紹基于地形的卡爾曼濾波狀態(tài)空間模型。

當(dāng)SAR干涉相位在離散時(shí)間情況下，考慮地形影響的卡爾曼濾波相位解纏狀態(tài)空間模型為

式中:A為給定狀態(tài)空間模型的系統(tǒng)矩陣;B為輸入控制矩陣;η為與地形有關(guān)的輸入控制變量;w(k)為狀態(tài)噪聲;E{}表示數(shù)學(xué)期望;Q(k)是狀態(tài)噪聲方差陣;x(k)為k點(diǎn)處的真實(shí)相位。這里，將SRTM DEM引入與地形有關(guān)的輸入控制變量η。在引入之前，需要對SRTM DEM進(jìn)行一系列處理，首先填充SRTM DEM上的空白區(qū)，然后與SAR主圖像配準(zhǔn)，再模擬相位，最終將其結(jié)果應(yīng)用于狀態(tài)空間模型中。

2 基于SRTM DEM的卡爾曼濾波相位解纏算法實(shí)施過程

基于SRTM DEM的二維卡爾曼濾波相位解纏算法的基本流程如圖1所示。

圖1 基于SRTM DEM的卡爾曼濾波相位解纏算法流程Fig.1 Flow chart of Kalman filter phase unwrapping based on SRTM DEM

具體計(jì)算步驟[13]如下:

式中:矩陣A和B已在式(1)中定義;η(k，k)為局部頻率估計(jì);Q(k，k)為狀態(tài)噪聲方差陣。這里，矩陣A和B均取單位矩陣。設(shè)和 P(1，1)分別為初始相位以及對應(yīng)方差陣(根據(jù)經(jīng)驗(yàn)設(shè)定)。

式中:I為單位矩陣;J(k+1)為卡爾曼濾波的增益矩陣;r(k+1，k+1)為殘差;C(k+1，k)為線性化的觀測矩陣。且

這里，R(k+1，k+1)為觀測噪聲方差陣，具體計(jì)算詳見文獻(xiàn)[12]。

若式(8)中計(jì)算出的|r(k+1，k+1)|＞π，則

在對SAR干涉圖進(jìn)行相位解纏時(shí)，由于任一估計(jì)相位都是通過2個(gè)相鄰區(qū)域計(jì)算得到的，所以在卡爾曼濾波相位解纏算法中采用式(11)(12)進(jìn)行估計(jì)，即

式中:Mwm=[0.5，0];Mwn=[0，0.5];ηm(m-1，n)和ηn(m，n-1)分別表示方位向和距離向上的局部頻率估計(jì)，這里采用的是SRTM DEM模擬相位在方位向和距離向上的差分代替局部頻率估計(jì)。

3 實(shí)驗(yàn)結(jié)果分析

以2003年6月和2004年1月伊朗Bam地區(qū)2景ENVISAT SAR圖像為數(shù)據(jù)源，從圖像上選取地形起伏較大、細(xì)節(jié)較豐富的150像元×150像元大小的區(qū)域(圖2(a))，其相干圖如圖2(b)所示(圖中矩形框區(qū)域中的黑色部分噪聲污染嚴(yán)重，相干系數(shù)在0.1左右)。經(jīng)過與SAR主圖像配準(zhǔn)的SRTM DEM模擬得到的相位圖如圖2(c)所示。

圖2 真實(shí)數(shù)據(jù)相位圖像Fig.2 Real data phase images

運(yùn)用原卡爾曼濾波相位解纏方法(簡稱原方法)[11]、CZT顧及地形因素卡爾曼濾波相位解纏方法(簡稱CZT方法)[12]、本文基于SRTM DEM 的卡爾曼濾波相位解纏方法(簡稱本文方法)和目前較為流行的Snaphu網(wǎng)絡(luò)流相位解纏算法[8]分別進(jìn)行解纏，解纏結(jié)果如圖3(a)—(h)所示。

圖3 真實(shí)數(shù)據(jù)解纏相位圖像(rad)Fig.3 Real data unwrapped phase images(rad)

圖4 和圖5分別表示上述前3種方法在方位向 和距離向的頻率估計(jì)結(jié)果。

圖4 真實(shí)數(shù)據(jù)在方位向的頻率估計(jì)結(jié)果Fig.4 Results of real data frequency estimation in azimuth direction

圖5 真實(shí)數(shù)據(jù)在距離向的頻率估計(jì)結(jié)果Fig.5 Results of real data frequency estimation in range direction

3.1 頻率估計(jì)結(jié)果分析

以圖4、圖5中真實(shí)數(shù)據(jù)的方位向和距離向頻率估計(jì)結(jié)果對比來看，原方法中固定分析窗口獲得的結(jié)果雖然平滑，但是在噪聲污染嚴(yán)重的局部區(qū)域出現(xiàn)了錯(cuò)誤的估計(jì)結(jié)果。無論是方位向(圖4(a))還是距離向((圖5(a))的結(jié)果圖上都存在大量的殘差點(diǎn);CZT方法(圖4(b)和圖5(b))在平滑度上有了一定的提高，但是仍然存在一些殘差點(diǎn);而本文方法(圖4(c)和圖5(c))無論在平滑度和估計(jì)精度上都有明顯的優(yōu)勢，在噪聲污染嚴(yán)重的區(qū)域(圖2(b)中矩形框區(qū)域)能夠獲得相對平滑的估計(jì)結(jié)果。

3.2 相位解纏結(jié)果分析

3.2.1 定性分析

從圖3可以看出，原方法的解纏結(jié)果(圖3(a)(e))有明顯的誤差傳遞，且丟失部分邊緣信息，“拉線”現(xiàn)象較嚴(yán)重;CZT方法結(jié)果(圖3(b)(f))雖然對地形進(jìn)行了有效估計(jì)，但是仍然有誤差傳遞;而本文方法結(jié)果(圖3(c)(g))更加平滑，說明圖2(b)中所框的噪聲污染嚴(yán)重區(qū)域已經(jīng)成功解纏，幾乎沒有誤差傳遞，與Snaphu方法的解纏結(jié)果(圖3(d)(h))比較接近。從解纏重纏繞結(jié)果(圖6)上看，Snaphu方法(圖6(d))中噪聲仍然存在，故Snaphu方法不能實(shí)現(xiàn)對SAR干涉圖的濾波;但是卡爾曼濾波相位解纏算法能夠在對SAR干涉圖解纏的同時(shí)實(shí)施濾波，而本文方法(圖6(c))濾波效果也好于原方法(圖6(a))和CZT方法(圖6(b))，噪聲幾乎完全被去除。這充分說明本文方法在解纏和濾波方面的優(yōu)勢。

圖6 解纏重纏繞結(jié)果(rad)Fig.6 Rewrapped phase images of unwrapping result(rad )

3.2.2 定量分析

在對解纏結(jié)果進(jìn)行定性分析的基礎(chǔ)上，運(yùn)用定量分析的方法，分別從運(yùn)算時(shí)間、不連續(xù)點(diǎn)數(shù)目和解纏質(zhì)量ε值[15]3個(gè)方面評價(jià)本文基于SRTM DEM方法的性能:運(yùn)算時(shí)間越短，表明算法計(jì)算效率越高;不連續(xù)點(diǎn)數(shù)目越少，說明抗相位畸變的性能越好;ε值越小，則解纏質(zhì)量越高。

這里給出評價(jià)解纏結(jié)果質(zhì)量的ε值表達(dá)式[15]為

式中:M和N分別為圖像的列數(shù)和行數(shù);φ(i，j)是點(diǎn)(i，j)處的解纏相位;分別是與纏繞相位梯度相對應(yīng)的權(quán)重，此權(quán)重一般由相位質(zhì)量圖給出，其值一般在0～1之間;p的取值通常有0，1，2，這里p取1。本文實(shí)驗(yàn)中以相干圖作為相位質(zhì)量圖。

表1給出了原方法、CZT方法、本文方法和Snaphu方法的運(yùn)算時(shí)間、不連續(xù)點(diǎn)數(shù)目和ε值。

表1 4種相位解纏算法的運(yùn)算時(shí)間、不連續(xù)點(diǎn)數(shù)目和ε值Tab.1 Computation time、number of discontinuous point and ε value of 4 phase unwrapping methods

圖7給出了原始纏繞相位及原方法、CZT方法、本文方法和Snaphu方法的不連續(xù)點(diǎn)分布圖，其中原始纏繞相位的不連續(xù)點(diǎn)數(shù)目約為2 239個(gè)。

圖7 不連續(xù)點(diǎn)分布圖(白色為不連續(xù)點(diǎn))Fig.7 Discontinuity maps (white is discontinuous point)

通過表1和圖7可以看出，本文方法不僅運(yùn)算時(shí)間和ε值均小于原方法和CZT方法的，而且不連續(xù)點(diǎn)數(shù)目最少(僅剩2點(diǎn));在與Snaphu方法計(jì)算出的ε值相當(dāng)?shù)那闆r下，本文方法得到的不連續(xù)點(diǎn)數(shù)目相對更少，表明本文方法抗相位畸變性能強(qiáng)，解纏質(zhì)量相對較高。

4 結(jié)論

本文針對現(xiàn)有卡爾曼濾波相位解纏算法在地形起伏較大或噪聲高時(shí)無法準(zhǔn)確反演地表形變信息這一情況，通過在卡爾曼濾波狀態(tài)空間模型中引入SRTM DEM模擬相位來估計(jì)地形因素對解纏的影響，提出了一種基于SRTM DEM的卡爾曼濾波相位解纏算法。該算法通過利用SRTM DEM地形信息指導(dǎo)解纏，提高了運(yùn)算的速度和精度。實(shí)際數(shù)據(jù)解纏結(jié)果表明，該算法的運(yùn)算效率和解纏質(zhì)量有了明顯的提高，尤其能夠提高低相干區(qū)域的解纏精度。

[1]張 紅，王 超，劉 智.星載合成孔徑雷達(dá)干涉測量[M].北京:科學(xué)出版社，2002:100-107.Zhang H，Wang C，Liu Z.Spaceborne synthetic aperture Radar interferometry[M].Beijing:Science Press，2002:100-107.

[2]Goldstein R M，Zerber H A，Werner C L.Satellite Radar interferometry:Two-dimensional phase unwrapping[J].Radio Science，1988，23(4):713-720.

[3]Xu W，Cumming I.A region-growing algorithm for InSAR phase unwrapping[J].IEEE Trans Geosci Remote Sens，1999，37(1):124-134.

[4]Flynn T J.Two-dimensional phase unwrapping with minimum weighted discontinuity[J].J Opt Soc Am A，1997，14(10):2692-2701.

[5]Ghiglia D C.Romero L A.Robust two-dimensional weighted and unweighted phase unwrapping that uses fast transforms and iterative methods[J].J Opt Soc Am A，1994，11(1):107-117.

[6]Pritt M D，Shipman J S.Least-squares two-dimensional phase unwrapping using FFTs[J].IEEE Trans Geosci Remote Sens，1994，32(3):706-708.

[7]Pritt M D.Phase unwrapping by means of multigrid techniques for interferometric SAR[J].IEEE Trans Geosci Remote Sens，1996，34(3):728-738

[8]Chen C W.Statistical-cost network-flow approaches to two-dimensional phase unwrapping for radar interferometry[D].California，USA，Stanford:Stanford University，2001.

[9]Kr?mer R.Auf kalman-filtern basierende phasen-und parameter estimation zur l?sung der phasen vieldeutigkeits problematik bei der H?henmodellerstellung aus SAR-interferogrammen[D].Germany，Siegen:Univer-sit?t-GH Siegen，1989.

[10]Loffeld O，Nies H，Knedlik S，et al.Phase unwrapping for SAR interferometry—a data fusion approach by Kalman filtering[J].IEEE Trans Geosci Remote Sens，2008，46(1):47-58.

[11]劉國林，郝華東，陶秋香.卡爾曼濾波相位解纏及其與其他方法的對比分析[J].武漢大學(xué)學(xué)報(bào):信息科學(xué)版，2010，35(10):1174-1178.Liu G L，Hao H D，Tao Q X.Kalman filter phase unwrapping algorithm and comparison and analysis with other methods[J].Geomatics and Information Science of Wuhan University，2010，35(10):1174-1178.

[12]劉國林，郝華東，閆 滿，等.顧及地形因素的卡爾曼濾波相位解纏算法[J].測繪學(xué)報(bào)，2011，40(3):283-288.Liu G L，Hao H D，Yan M，et al.Phase unwrapping algorithm by using Kalman filter based on topographic factors[J].Acta Geodaetica et Cartographica Sinica，2011，40(3):283-288.

[13]郝華東.卡爾曼濾波在InSAR相位解纏中的應(yīng)用研究[D].青島:山東科技大學(xué)，2010.Hao H D.Study of Kalman filter applied in phase unwrapping of InSAR[D].Qingdao:Shandong University of Science and Technology，2010.

[14]詹總謙，錢 俊，舒 寧.基于區(qū)域分割和編碼的相位解纏方法[J].武漢大學(xué)學(xué)報(bào):信息科學(xué)版，2002，27(3):316-320.Zhan Z Q，Qian J，Shu N.A method of phase unwrapping based on region-segmentation and encoding[J].Geomatics and Information Science of Wuhan University，2002，27(3):316-320.

[15]Ghiglia D C，Pritt M D.Two-dimensional phase unwrapping:Theory，algorithms and software[M].New York:John Wiley ＆ Sons Inc，1998.