宋 瑞，劉 廣，PERSKI Zbigniew，郭華東

一種改進(jìn)的SAR干涉圖綜合濾波算法研究

宋 瑞1，劉 廣2，PERSKI Zbigniew3，郭華東2

(1.北京航空航天大學(xué)儀器科學(xué)與光電工程學(xué)院，北京 100191;2.中國(guó)科學(xué)院對(duì)地觀測(cè)與數(shù)字地球科學(xué)中心，北京 100094;3.Polish Geological Institute National Research Institute，Carpathian Branch，Poland 31560)

SAR干涉圖作為相位信息的載體，其質(zhì)量直接影響對(duì)研究區(qū)域形變狀況的進(jìn)一步分析，采取有效的濾波算法能抑制干涉圖相位噪聲，提高干涉測(cè)量精度。在獲得的干涉相位圖中，由于礦區(qū)開(kāi)采而造成的地表沉降會(huì)體現(xiàn)出近環(huán)狀相位條紋的特征。針對(duì)這一特點(diǎn)，對(duì)傳統(tǒng)的基于梯度的濾波算法做出了改進(jìn)，并結(jié)合Goldstein頻域?yàn)V波和改進(jìn)的梯度自適應(yīng)濾波，提出了一種適用于礦區(qū)沉降形成的SAR干涉相位模式濾波方法。選取河北峰峰煤礦的PALSAR干涉相位圖作為實(shí)驗(yàn)數(shù)據(jù)，對(duì)該濾波方法做出了詳細(xì)的性能評(píng)價(jià)和對(duì)比。結(jié)果表明，采用本文提出的綜合濾波方法在顯著降低實(shí)驗(yàn)區(qū)SAR干涉圖相位噪聲的同時(shí)，也很好地保持了相位分辨率，使由于礦區(qū)沉降而造成的形變相位環(huán)的邊緣形態(tài)更加清晰。

SAR干涉圖;Goldstein濾波;改進(jìn)的梯度自適應(yīng)濾波;礦區(qū)地表沉降

0 引言

合成孔徑雷達(dá)干涉測(cè)量(InSAR)技術(shù)綜合了合成孔徑雷達(dá)成像原理和干涉測(cè)量技術(shù)，能夠全天時(shí)、全天候?qū)δ繕?biāo)進(jìn)行觀測(cè)，具有覆蓋范圍廣、空間分辨率高、可以實(shí)施大范圍內(nèi)連續(xù)地表監(jiān)測(cè)的特點(diǎn)［1］，在大地測(cè)量界得到廣泛關(guān)注。合成孔徑雷達(dá)差分干涉測(cè)量(DInSAR)技術(shù)是InSAR應(yīng)用的一個(gè)拓展，理論上可以得到非常精確的數(shù)字高程模型和毫米量級(jí)的地表形變信息［2］。

我國(guó)近幾年礦業(yè)迅速發(fā)展，但一些次生災(zāi)害如地表沉降問(wèn)題也隨之出現(xiàn)，它不僅影響到周邊的環(huán)境，也給國(guó)家造成了巨大的經(jīng)濟(jì)損失。目前，采礦工程和巖石工程領(lǐng)域基于大量數(shù)學(xué)方法已經(jīng)衍生出一些對(duì)礦區(qū)地表沉降預(yù)測(cè)的方法［3－4］，但都存在一定的局限性。相比傳統(tǒng)手段，將DInSAR技術(shù)應(yīng)用于礦區(qū)地表沉降的監(jiān)測(cè)具有很多優(yōu)勢(shì)［5－6］，不僅可以降低成本，而且可以準(zhǔn)時(shí)動(dòng)態(tài)監(jiān)測(cè)整個(gè)雷達(dá)圖像覆蓋范圍內(nèi)地表沉降偏移量。目前，國(guó)際上一些學(xué)者和研究團(tuán)體正致力于尋求地表沉降和地下開(kāi)采活動(dòng)的相關(guān)關(guān)系［7－9］。

通過(guò)DInSAR技術(shù)得到的干涉圖不可避免地存在相位噪聲，這直接影響了相位解纏的精度和計(jì)算效率，對(duì)干涉測(cè)量的結(jié)果分析也造成了干擾。經(jīng)過(guò)近幾十年的發(fā)展，針對(duì)InSAR干涉圖濾波降噪問(wèn)題已經(jīng)涌現(xiàn)出很多算法。經(jīng)典的 Lee 濾波［10］和Weiner［11－12］濾波是以最小均方誤差估計(jì)為理論基礎(chǔ)，通過(guò)建立線性噪聲模型實(shí)現(xiàn)降噪;常用的均值濾波［13］、中值濾波［14－15］和梯度濾波［16－17］是局部統(tǒng)計(jì)自適應(yīng)濾波的典型代表;頻域的Goldstein濾波［18］通過(guò)分析噪聲的高頻響應(yīng)，利用噪聲不同于一般信號(hào)的頻譜特性，對(duì)頻率域的系數(shù)進(jìn)行處理，實(shí)現(xiàn)降噪的目的。這些算法都具有一定局限性，在實(shí)際運(yùn)用中并不能夠達(dá)到理想的濾波效果。分析采礦原理與地表沉降之間的關(guān)系發(fā)現(xiàn)，在InSAR干涉圖中由于礦區(qū)開(kāi)采造成的地表沉降將會(huì)呈現(xiàn)出典型的近環(huán)狀干涉條紋。如何針對(duì)這種干涉圖進(jìn)行濾波，在減少噪聲的同時(shí)，更好地從中提取形變相位環(huán)的完整信息是我們關(guān)心的重點(diǎn)，這對(duì)進(jìn)一步分析礦區(qū)沉降深度、沉降面積以及精確的沉降地理方位都有重要的意義。

本文對(duì)傳統(tǒng)的基于梯度的濾波算法做出了改進(jìn)，改進(jìn)后的算法在保持環(huán)形邊緣的細(xì)節(jié)方面具有良好效果。采用頻域Goldstein濾波和改進(jìn)后的梯度自適應(yīng)濾波相結(jié)合的二級(jí)濾波算法對(duì)礦區(qū)沉降監(jiān)測(cè)的SAR干涉圖進(jìn)行處理，有效濾除了用于礦區(qū)沉降監(jiān)測(cè)的SAR干涉圖的相位噪聲。

1 干涉圖濾波算法

1.1 Goldstein干涉圖濾波

干涉圖中往往含有大量高頻噪聲，而干涉條紋的空間頻率則相對(duì)穩(wěn)定?；谠肼朁c(diǎn)和一般信息點(diǎn)頻率特性的不同，將濾波處理從常用的空間域轉(zhuǎn)換到頻率域，對(duì)噪聲點(diǎn)顯著的頻譜特性分析處理，能夠?qū)崿F(xiàn)對(duì)干涉圖的頻域?yàn)V波。Goldstein和Werner提出了一種頻域加權(quán)的干涉圖濾波方法［18］，該算法能較好地抑制InSAR干涉圖的相位噪聲，同時(shí)兼顧了保持相位邊緣信息的優(yōu)點(diǎn)。

首先，在干涉圖復(fù)圖像B(x，y)中選取合適的濾波窗口，通過(guò)二維傅里葉變換(FFT)獲得圖像的頻譜數(shù)據(jù)BW(u，v)(u和v表示空間頻率)進(jìn)一步得到干涉圖的頻譜幅值|BW(u，v)|，即

然后選取0～1之間合適的濾波參數(shù)，通過(guò)式對(duì)頻譜函數(shù)進(jìn)行加權(quán)處理，式中α為Goldstein濾波的關(guān)鍵參數(shù);

最后，對(duì)加權(quán)處理后的干涉圖頻譜進(jìn)行二維傅里葉反變換(IFFT)，得到空間域?yàn)V波后的干涉圖BG(x，y)，即

Goldstein的濾波參數(shù)α可以通過(guò)經(jīng)驗(yàn)人為設(shè)定。實(shí)驗(yàn)結(jié)果表明，隨著α的增大(不超過(guò)1)，濾波效果增強(qiáng)，但是帶來(lái)的負(fù)面影響則會(huì)使圖像的分辨率降低，圖像的一些區(qū)域?qū)⒆兊媚：?。所以，α的取值?yīng)該適當(dāng)，在達(dá)到較好的降噪性能的同時(shí)要保證圖像仍然具有較高的分辨率，使兩者之間達(dá)到平衡。

1.2 基于梯度的自適應(yīng)濾波

基于梯度的自適應(yīng)濾波是以每個(gè)點(diǎn)的梯度作為衡量標(biāo)準(zhǔn)，通過(guò)梯度計(jì)算模板中每個(gè)像素點(diǎn)對(duì)模板中心像素點(diǎn)的貢獻(xiàn)。對(duì)復(fù)數(shù)形式存儲(chǔ)的干涉圖數(shù)據(jù)，首先應(yīng)將原始數(shù)據(jù)的實(shí)部與虛部進(jìn)行分離，然后通過(guò)以下3個(gè)步驟，完成后將實(shí)部與虛部重組，生成濾波后的干涉圖。

1)以一階差分計(jì)算每個(gè)像素點(diǎn)在水平和垂直方向的梯度分量，即

式中:Gx(x，y)為像素點(diǎn)的水平梯度分量;Gy(x，y)為像素點(diǎn)的垂直梯度分量。計(jì)算圖像邊界像素點(diǎn)的梯度時(shí)，可以先采用最鄰近法對(duì)圖像的邊界進(jìn)行整體擴(kuò)充，擴(kuò)充后便可計(jì)算出原始邊界像素點(diǎn)的梯度。

2)結(jié)合每個(gè)點(diǎn)的方向梯度，采用

計(jì)算出相應(yīng)的權(quán)系數(shù)。式中k值代表該濾波算法的平滑尺度，決定了保留多大幅值的邊緣。當(dāng)某點(diǎn)的梯度大于k時(shí)，視該點(diǎn)為邊緣信息點(diǎn)，將其權(quán)系數(shù)置為0，該點(diǎn)不參與平滑處理;梯度小于k的點(diǎn)按式(7)計(jì)算權(quán)系數(shù)，導(dǎo)致非邊緣點(diǎn)權(quán)系數(shù)變大，較大程度參與了平滑，達(dá)到既去除噪聲又保留邊緣的效果。通常情況下，可以通過(guò)統(tǒng)計(jì)整張圖像中每個(gè)像素點(diǎn)的梯度后，以平均梯度的一定倍數(shù)作為k的大小，但具體取值要根據(jù)圖像特征而定。

3)得到權(quán)系數(shù)矩陣后，通過(guò)式(8)計(jì)算模板內(nèi)中心點(diǎn)濾波后的像素值，即式中:f(x，y)為干涉圖的實(shí)部或虛部;n為濾波算法的迭代次數(shù)，取值應(yīng)根據(jù)相位干涉圖受噪聲污染的程度而定，一般迭代3～5次即可達(dá)到較好效果。

1.3 改進(jìn)的梯度自適應(yīng)濾波

由于傳統(tǒng)的基于梯度濾波算法只考慮了水平和垂直兩個(gè)方向，不能兼顧其他方向的梯度分量，求得的方向往往偏向于水平和垂直方向，在礦區(qū)沉降的干涉相位圖中，不能很好地跟蹤近環(huán)狀邊緣的走向。為了兼顧其他方向上的精確性，對(duì)傳統(tǒng)的基于梯度的濾波算法做出了改進(jìn)，增加了兩個(gè)方向的梯度，采用 0°，90°，135°和 45°方向的梯度來(lái)計(jì)算每個(gè)像素點(diǎn)的權(quán)重值，則有

采用上述4個(gè)方向的梯度計(jì)算權(quán)系數(shù)值時(shí)，每個(gè)點(diǎn)包含的梯度信息更加豐富，尤其針對(duì)本文所研究的近環(huán)狀邊緣更加有效，因此，對(duì)應(yīng)的權(quán)系數(shù)公式應(yīng)變?yōu)?/p>

另外，使用3×3模板通過(guò)加權(quán)函數(shù)式(8)進(jìn)行平滑時(shí)，由于在一些邊緣信息點(diǎn)附近同時(shí)存在噪聲點(diǎn)，一個(gè)模板內(nèi)可能會(huì)有大部分像素點(diǎn)的權(quán)系數(shù)會(huì)被置零，因此模板內(nèi)真正參加平滑的像素點(diǎn)會(huì)大大減少，這會(huì)影響濾波結(jié)果的精度，甚至出現(xiàn)錯(cuò)誤值，這是該算法產(chǎn)生誤差的最大因素。為了減小這一誤差，使相位圖中近環(huán)狀邊緣的信息更加完整，對(duì)此算法做出第二點(diǎn)改進(jìn)，即通過(guò)可變矩形的濾波窗口實(shí)現(xiàn)完全自適應(yīng)的梯度濾波。首先，在模板內(nèi)計(jì)算中心點(diǎn)的像素值前，對(duì)權(quán)系數(shù)不為0的像素點(diǎn)個(gè)數(shù)做出判斷，如果實(shí)際參與式(8)計(jì)算的點(diǎn)數(shù)小于3×3模板內(nèi)像素點(diǎn)總數(shù)的一半(約為5個(gè))時(shí)，將濾波窗口擴(kuò)大為5×5;如果擴(kuò)大后的濾波窗口內(nèi)權(quán)系數(shù)不為零的點(diǎn)數(shù)仍然小于5，可以對(duì)窗口再進(jìn)行擴(kuò)充為9×9。Smax為最大濾波窗口，可以通過(guò)設(shè)定Smax限制濾波窗口的尺度，實(shí)驗(yàn)中定義Smax為9×9較為合適。在不同大小的濾波窗口下，原始濾波函數(shù)變換為

式中i代表濾波窗口邊長(zhǎng)。

1.4 綜合濾波算法

本文提出的Goldstein濾波和改進(jìn)的梯度自適應(yīng)濾波相結(jié)合的綜合濾波方法的整體流程為:首先，通過(guò)Goldstein濾波算法，利用FFT變換得到復(fù)圖像的頻譜數(shù)據(jù)，用加權(quán)函數(shù)對(duì)頻譜數(shù)據(jù)進(jìn)行處理后，由FFT反變換得到濾波后的復(fù)圖像;然后，將頻域?yàn)V波后的復(fù)數(shù)干涉圖分離成實(shí)部和虛部?jī)蓚€(gè)部分，分別通過(guò)平滑尺度k、迭代次數(shù)n和最大濾波窗口Smax的控制對(duì)數(shù)據(jù)集進(jìn)行改進(jìn)的局部梯度自適應(yīng)濾波;最后，將實(shí)部與虛部重組，生成濾波后干涉圖，通過(guò)干涉圖復(fù)數(shù)據(jù)計(jì)算干涉相位。

頻域Goldstein以噪聲的高頻特性為基礎(chǔ)，濾波過(guò)程中在邊緣保持方面效果較好;梯度自適應(yīng)濾波有較強(qiáng)的濾除孤立噪聲點(diǎn)的特性，但是在圖像邊緣與噪聲混疊的情況下，可能對(duì)圖像邊緣造成一定破壞，算法改進(jìn)后邊緣信息得到豐富，但是仍然存在離散分布的誤差區(qū)域。采用二級(jí)濾波手段，通過(guò)Goldstein頻域?yàn)V波相當(dāng)于先對(duì)圖像作了一次預(yù)處理，將邊緣處的噪聲密度大大降低，接著再進(jìn)行改進(jìn)的梯度自適應(yīng)濾波處理圖像。實(shí)驗(yàn)結(jié)果表明，綜合濾波方式使礦區(qū)沉降區(qū)域的干涉相位環(huán)更加清晰、邊緣保持效果更好，同時(shí)大大降低了相位殘差點(diǎn)個(gè)數(shù)，為礦區(qū)沉降信息的進(jìn)一步分析帶來(lái)了很大的方便。

2 煤礦開(kāi)采與礦區(qū)地表沉降的關(guān)系

由于地質(zhì)狀況的差異，不同礦區(qū)開(kāi)采后的沉降程度也有所不同。通常在采空區(qū)形成后的1～3個(gè)月內(nèi)沉降最為迅速，通常能達(dá)到每天幾個(gè)厘米，而后逐漸趨于穩(wěn)定，殘余的沉降將持續(xù)近一年［20］。通過(guò)SAR對(duì)地表大范圍成像，理想情況下礦區(qū)地表形變會(huì)在干涉圖上清楚地反映出來(lái)，從而實(shí)現(xiàn)對(duì)地表形變及地下開(kāi)采活動(dòng)的監(jiān)測(cè)。圖1描述了煤礦開(kāi)采與礦區(qū)地表沉降之間的關(guān)系。

圖1 煤礦開(kāi)采與礦區(qū)地表沉降的關(guān)系［19］Fig.1 Demonstration of mine extraction and subsidence［19］

從圖中可以看出，發(fā)生沉降的礦區(qū)地表形成了一個(gè)塌陷槽，其采空中心具有最大沉降距離。受到塌陷的影響，沉降深度由采空中心向四周逐漸衰減，形成了具有一定塌陷角的漏斗形地形結(jié)構(gòu)，其整體塌陷面積受到塌陷角大小的影響。因此，根據(jù)等高線原理，在通過(guò)DInSAR技術(shù)獲得的雷達(dá)干涉圖中，發(fā)生地表沉降的采礦區(qū)會(huì)表現(xiàn)為近環(huán)狀相位條紋，圓環(huán)中心具有較大的相位信息，隨著與沉降中心圓心距的逐漸增大，環(huán)狀結(jié)構(gòu)向外逐層擴(kuò)散，與其相對(duì)應(yīng)的相位值也不斷衰減，意味著地表沉降距離由中心向外部逐漸減少，在SAR干涉相位圖中整體表現(xiàn)為一個(gè)形變信息規(guī)律變化的相位環(huán)。根據(jù)以上對(duì)礦區(qū)沉降SAR干涉圖特征的分析，有針對(duì)性的采用最佳的噪聲濾波方法，使圖像中礦區(qū)沉降產(chǎn)生的形變相位環(huán)更加清晰、邊緣更加完整，將大大提高通過(guò)干涉相位圖分析獲得礦區(qū)沉降深度及沉降范圍的精度。

3 綜合濾波算法實(shí)驗(yàn)

選取河北峰峰煤礦2007年7月30日和2007年9月14日間隔46 d的ALOS的PALSAR合成孔徑雷達(dá)數(shù)據(jù)，通過(guò)DInSAR技術(shù)獲得其形變干涉圖，圖像覆蓋了礦區(qū)的幾個(gè)重要采礦點(diǎn)。

首先，分別采用改進(jìn)前后的梯度自適應(yīng)算法進(jìn)行濾波實(shí)驗(yàn)。圖2選取的干涉相位圖包含了3處明顯的近環(huán)狀形變條紋，用于對(duì)比分析改進(jìn)后的算法在邊緣細(xì)節(jié)保持方面體現(xiàn)出的優(yōu)勢(shì)。

圖2(a)采用原始的梯度濾波算法，此算法僅考慮水平和垂直方向的梯度，信息過(guò)于單一，導(dǎo)致部分相位環(huán)邊緣出現(xiàn)塊狀的誤差區(qū)域，采用改進(jìn)的梯度自適應(yīng)濾波算法后多數(shù)塊狀誤差區(qū)域的面積有所減少甚至消除(圖2(b))?？傮w上，通過(guò)增加45°，135°梯度信息的計(jì)算和做出可變?yōu)V波窗口的調(diào)整，圖像在濾波后整體體現(xiàn)的層級(jí)信息更加分明。

首先對(duì)形變干涉圖進(jìn)行Goldstein濾波，濾波窗口大小為32×32，為了兼顧濾波的降噪特性和保持圖像分辨率的要求，濾波參數(shù)α在0.4～0.6之間調(diào)整，不宜過(guò)大或過(guò)小，具體的取值對(duì)濾波結(jié)果的影響將在后文進(jìn)行分析。將Goldstein濾波后的干涉復(fù)圖像存為實(shí)部和虛部?jī)蓚€(gè)數(shù)據(jù)集，分別進(jìn)行改進(jìn)的梯度自適應(yīng)濾波。綜合濾波結(jié)束后，再通過(guò)實(shí)部和虛部求解相位圖，觀察濾波前后相位圖所發(fā)生的改變。實(shí)驗(yàn)結(jié)果將與單獨(dú)采用Goldstein濾波和改進(jìn)的梯度自適應(yīng)濾波的結(jié)果進(jìn)行對(duì)比，以一維剖面圖和相位殘差點(diǎn)［21］作為評(píng)價(jià)指標(biāo)進(jìn)行了定量分析，最重要的是通過(guò)目視解譯的方法［22］，觀察采用本文的濾波方法后，核心礦區(qū)在相位圖中所體現(xiàn)出的更加清晰的近環(huán)狀形變紋理。

圖3為原始干涉相位圖(圖3(a))與幾種濾波方法所得結(jié)果圖(圖3(b)－(d))的對(duì)比，可以對(duì)濾波效果做出目視評(píng)價(jià)。

圖2 改進(jìn)前后梯度自適應(yīng)濾波后的相位圖Fig.2 Interferogram using original and the improved gradient－ based adaptive filter

圖3 (1) 不同方法濾波的干涉相位圖Fig.3 (1)Interferogram using different filters

圖3(2) 不同方法濾波的干涉相位圖Fig.3 (2)Interferogram using different filters

圖3 (b)采用了Goldstein濾波，對(duì)圖像的邊緣細(xì)節(jié)保持良好，幾處明顯的沉降環(huán)在濾波結(jié)束后依然很完整，但是降噪性能有限，整體相位圖的噪聲密度仍然較高。圖3(c)采用了改進(jìn)的梯度自適應(yīng)濾波，圖像噪聲密度雖然相比圖3(b)降低了很多，但是圖像細(xì)節(jié)破壞嚴(yán)重，幾處沉降環(huán)內(nèi)出現(xiàn)了塊狀的誤差區(qū)域，降低了圖像質(zhì)量;并且，在圖像噪聲點(diǎn)密度較高的情況下，直接采用梯度自適應(yīng)濾波需要迭代多次才能達(dá)到精度要求，濾波效率也受到一定限制，因此不宜單獨(dú)使用梯度自適應(yīng)濾波。圖3(d)采用了本文提出的Goldstein濾波和改進(jìn)的梯度自適應(yīng)濾波相結(jié)合的綜合濾波方法，可以看出無(wú)論是邊緣、細(xì)節(jié)的保持，還是降噪方面，相比圖3(b)，(c)都得到了很大的提高。這是因?yàn)槭紫冉?jīng)過(guò)Goldstein頻域?yàn)V波預(yù)處理，使礦區(qū)沉降環(huán)的形態(tài)更加清晰明顯，噪聲密度也得到一定程度的下降;接著進(jìn)行的梯度自適應(yīng)濾波中，由于上一步的處理使圖像的邊緣噪聲得到了分離，此時(shí)進(jìn)行改進(jìn)的梯度自適應(yīng)濾波時(shí)不會(huì)出現(xiàn)圖3(c)中圓環(huán)形態(tài)受到破壞的結(jié)果;并且因?yàn)镚oldstein濾波降低了圖像其他范圍的噪聲密度，再進(jìn)行梯度自適應(yīng)濾波時(shí)孤立噪聲的分布已經(jīng)相比原始圖像更加分散，更適合基于梯度的濾波算法來(lái)處理;在算法效率方面，經(jīng)過(guò)Goldstein濾波的預(yù)處理后，改進(jìn)的梯度自適應(yīng)濾波的迭代次數(shù)減少，效率也得到提高。綜合分析圖3(d)，本文算法將兩種濾波方法的優(yōu)點(diǎn)進(jìn)行了有機(jī)的結(jié)合，得到了更加理想的濾波效果，在針對(duì)礦區(qū)沉降干涉相位圖濾波方面具有較大優(yōu)勢(shì)。

通過(guò)一維剖面圖可以對(duì)濾波結(jié)果進(jìn)行更直觀的分析。選擇第240列對(duì)濾波前后各像素的相位值進(jìn)行了統(tǒng)計(jì)(圖4)。

圖4 不同方法濾波后干涉圖的第240列剖面圖Fig.4 Profile chart of the 240th column of SAR interferogram using different filters

圖4中相位值突變起伏的點(diǎn)被認(rèn)為是噪聲點(diǎn)，而經(jīng)過(guò)濾波算法后，相位值變化斜率仍然較大的點(diǎn)被認(rèn)為是邊緣信息點(diǎn)。對(duì)統(tǒng)計(jì)圖進(jìn)行具體的分析可以得到，濾波之前相位圖中還有大量相位噪聲，因此圖4(a)中的相位變化曲線沒(méi)有表現(xiàn)出一定的規(guī)律性，噪聲點(diǎn)與邊緣點(diǎn)發(fā)生混疊，導(dǎo)致相位信息一直存在較大的起伏，無(wú)法明確地體現(xiàn)邊緣出現(xiàn)的位置。圖4(b)中，經(jīng)過(guò)Goldstein濾波后大量噪聲點(diǎn)被剔除，邊緣信息變得更加明顯，細(xì)節(jié)方面相比(a)圖也更加平滑，但依然存在一些毛刺。圖4(c)采用改進(jìn)后的梯度自適應(yīng)濾波，雖然相比圖4(b)毛刺更少，表面上有更好的降噪效果，但這也是采用梯度自適應(yīng)濾波后圖像細(xì)節(jié)受損的原因;由于圖像中梯度的區(qū)域分布特征，導(dǎo)致圖像濾波后一些噪聲點(diǎn)及正常信息點(diǎn)都成塊狀分布，使剖面統(tǒng)計(jì)圖顯得更加平滑，因此單獨(dú)采用梯度自適應(yīng)濾波會(huì)降低相位圖的質(zhì)量。圖4(d)采用本文的綜合濾波方法，取得了較好的濾波效果，既發(fā)揮了Goldstein濾波的平滑優(yōu)點(diǎn)，又達(dá)到了梯度自適應(yīng)濾波去除毛刺的效果，同時(shí)一定程度上保證了干涉相位圖較高的分辨率。

1998年，Goldstein將干涉相位圖中噪聲引起的相位不連續(xù)定義為殘差點(diǎn)，并利用殘差點(diǎn)的數(shù)目來(lái)表征相位圖質(zhì)量［21］。其原理是在干涉相位圖中統(tǒng)計(jì)由鄰近4個(gè)像素組成的閉合環(huán)纏繞相位差，若殘差和不為零，則說(shuō)明有殘差點(diǎn)存在。表1是通過(guò)統(tǒng)計(jì)相位殘差點(diǎn)的個(gè)數(shù)對(duì)濾波效果進(jìn)行的定量分析結(jié)果，其中采用本文綜合算法濾波后相位殘差點(diǎn)數(shù)最少，客觀地體現(xiàn)了這種綜合濾波方法的優(yōu)勢(shì)。

表1 不同濾波算法的定量評(píng)價(jià)Tab.1 Quantitative evaluation of different filtering algorithms

為了達(dá)到濾波效果并保持圖像分辨率，實(shí)驗(yàn)過(guò)程中將α取值設(shè)置在0.4～0.6區(qū)間內(nèi)，k的取值分別達(dá)到10，30，50倍平均方向梯度，并對(duì)采用以上組合參數(shù)進(jìn)行濾波的效果進(jìn)行了定量統(tǒng)計(jì)和對(duì)比分析，以求得到最佳的濾波參數(shù)。不同參數(shù)濾波后的干涉相位圖如圖5所示;不同參數(shù)濾波后的相位殘差點(diǎn)統(tǒng)計(jì)見(jiàn)表2。

圖5 不同參數(shù)濾波后的干涉相位圖Fig.5 Interferogram using different filtering parameters

表2 不同參數(shù)濾波后的相位殘差點(diǎn)統(tǒng)計(jì)Tab.2 Statistics of phase residual using different filtering parameters

通過(guò)圖5的干涉相位圖的直觀分析和表2相位殘差點(diǎn)的定量統(tǒng)計(jì)可以看出，α=0.5時(shí)Goldstein頻域?yàn)V波效果最好，此時(shí)近環(huán)狀相位條紋的內(nèi)層紋理最為清晰;α＜0.5時(shí)相位殘差點(diǎn)個(gè)數(shù)稍有增多，沒(méi)有達(dá)到最佳濾波效果;α＞0.5時(shí)圖像分辨率有所下降，邊緣信息略顯模糊，因此對(duì)應(yīng)的相位殘差點(diǎn)個(gè)數(shù)大幅度減少，導(dǎo)致圖像細(xì)節(jié)有所損失。所以，選擇α=0.5進(jìn)行Goldstein濾波最為適宜。

改進(jìn)的梯度自適應(yīng)濾波中k值取10倍平均梯度時(shí)平滑尺度較高，一些有小幅度梯度變化的像素點(diǎn)被當(dāng)作是邊緣像素點(diǎn)而不參與平滑，導(dǎo)致濾波窗口調(diào)為上限后模板內(nèi)參與平滑的像素點(diǎn)數(shù)仍然較少的情況增多，所以平滑后得到的相位值會(huì)出現(xiàn)較大誤差，多次迭代后甚至可能擴(kuò)大誤差范圍，形成塊狀噪斑。從圖5可以看到，k取10倍梯度時(shí)，圖像中塊狀分布的相位區(qū)域較多。如果k值過(guò)高，達(dá)到50倍平均梯度，部分低幅度邊緣會(huì)參與平滑，給周圍像素點(diǎn)的濾波都帶來(lái)影響，破壞邊緣信息的情況更為嚴(yán)重，表2顯示相位殘差點(diǎn)數(shù)也最多。由于噪聲與邊緣的混疊，導(dǎo)致圖5中顯示10倍和50倍的k值濾波后邊緣都有一定程度的破壞，k取30時(shí)干涉相位圖中的近環(huán)狀邊緣形態(tài)保持較好。因此，經(jīng)過(guò)以上分析，本文α取0.5，k取30作為濾波參數(shù)能得到最好的干涉圖濾波效果。

上文幾種濾波方法效果的對(duì)比分析表明，將Goldstein頻域?yàn)V波和改進(jìn)的梯度自適應(yīng)濾波相結(jié)合的綜合濾波方式應(yīng)用于礦區(qū)沉降監(jiān)測(cè)干涉圖濾波具有較明顯的優(yōu)勢(shì)。采用這種綜合的濾波方法對(duì)獲得的礦區(qū)形變干涉數(shù)據(jù)進(jìn)行濾波，可為礦區(qū)的沉降形變量和形變范圍分析提供更加清晰的相位圖，同時(shí)提高干涉圖像結(jié)合應(yīng)用光學(xué)圖像進(jìn)行礦區(qū)定位的精度。

4 結(jié)論

1)針對(duì)礦區(qū)沉降形成的SAR干涉相位模式所體現(xiàn)出的獨(dú)特特征，本文對(duì)傳統(tǒng)的基于梯度的濾波算法做出了改進(jìn)，并提出了一種Goldstein頻域?yàn)V波和改進(jìn)的梯度自適應(yīng)濾波相結(jié)合的綜合濾波方法。選取PALSAR干涉相位數(shù)據(jù)進(jìn)行了實(shí)驗(yàn)，對(duì)不同濾波方法得到的結(jié)果進(jìn)行了對(duì)比和定量分析，得到了這種綜合濾波算法的最佳濾波參數(shù)。

2)實(shí)驗(yàn)結(jié)果表明，將Goldstein濾波和改進(jìn)的梯度自適應(yīng)的濾波方法相結(jié)合，可以使因礦區(qū)沉降造成的形變相位環(huán)更加清晰，在較好地保持相位分辨率的同時(shí)能夠最大限度地降低噪聲干擾，發(fā)揮了兩種濾波方法的優(yōu)勢(shì)。因此，相比其他的單一濾波手段，本文提出的綜合濾波方法更適合應(yīng)用于包含礦區(qū)沉降信息的形變干涉圖的濾波。

志謝:感謝JAXA提供PALSAR數(shù)據(jù)。

［1］ Graham L C，Synthetic Interferometer Radar for Topographic Mapping［J］.Proc IEEE，1974，62(6):763 －768.

［2］ Rodriguez E，Martin J M.Theory and Design of Interferometric Synthetic Aperture Radars［J］.IEE Proc:F，1992，139(2):147 －159.

［3］ Oravecz K I.Improved Prediction of Surface Subsidence Using the Infuence Function Approach［C］∥ Sangorm Symposium:The Effect of Underground Mining on Surface，1986:73 －80.

［4］ Alejano L R，Ranirez － Oyanguren P，Taboada J，et al.Numerical Prediction of Subsidence Phenomena Due to Flat Coal Seam Mining［J］.InternationalJournalofRock Mechanicsand Mining Sciences，1998，35(4):440 －441.

［5］ Strozzi T，Wegmuller U，Werner C L，et al.JERS SAR Interferometry for Land Subsidence Monitoring［J］.IEEE Transactions on Geoscience and Remote Sensing，2003，41(7):1702 －1703.

［6］ Chang Z Q，Gong H L，Zhang J F，et al.A Feasible Approach for Improving Accuracy of Ground Deformation Measured by DInSAR［J］.Journal of China University of Mining and Technology，2007，17(2):262 －266.

［7］ Chang H C，Ge L L，Rizos C.DInSAR for Mine Subsidence Monitoring Using Multi－ source Satellite SAR Images［C］∥Geoscience and Remote Sensing Symposium.Seoul，Korea，2005.

［8］ Perski Z.Application of SAR Imagery and SAR Interferometry in Digital Geological Cartography［M］∥Ostaficzuk S.The Current role of GeologicalMapping in Geoscience.Netherlands Springer，Nato Science Series，2005，56:225 －244.

［9］ Wright P，Stow R.Detecting Mining Subsidence from Space［J］.International Journal of Remote Sensing，1999，20(6):1183 －1188.

［10］ Lee J S，Papathanassiou K P，Ainsworth T L，et al.A New Technique for Noise Filtering of SAR Interferometric Phase Images［J］.IEEE Transactions on Geoscience and Remote Sensing，1998，36(5):1456－1465.

［11］ Joo T H，Held D N.SAR Speckle Noise Reduction Using Wiener filter［J］.Proc of the NASA Symp on Math Pattern Recognition and Image Analysis，1983，295 －310.

［12］ Tates M，Nasrabadi N，Kwon H，et al.Wiener Filter－ based Change Detection for SAR Imagery［J］.Proceedings of SPIE，2006，6237:62370N1 －62370N8.

［13］王曉丹，吳崇明.基于MATLAB的系統(tǒng)分析與設(shè)計(jì)［M］.西安:西安電子科技大學(xué)出版社，2000:9－10.

［14］ Gabbouj M，Coyle E J，Gallagher N C.An Overview of Median and Stack Filtering［J］.Circuit Systems Signal Processing，1992，11(1):7－45.

［15］ Rees W G，Satchell M J F.The Effect of Median Filtering on Synthetic Aperture Radar Images［J］.International Journal of Remote Sensing，1997，11(3):2887 －2893.

［16］ Eichel P H，Ghiglia D C，et al.Spotlight SAR Inter－ ferometric Change Detection［R］.Sandia National Labs Tech Report，1993:2539－2546.

［17］林 卉，趙長(zhǎng)勝，杜培軍，等.InSAR干涉圖濾波方法研究［J］.測(cè)繪學(xué)報(bào)，2005，34(2):113 －117。

［18］ Goldstein R M，Werner C L.Radar Interferogram Filtering for Geophysical Applications［J］.Geophysical Research Letters，1998，25(21):4035－4038.

［19］ Perski Z.Application of SAR Imagery and SAR Interferometry in Digital Geological Cartography［J］.The Current Role of Geological Mapping in Geoscience，Netherlands:Springer，NatoScience Series，2005.

［20］ Liu G，Guo H D，Ramon H，et al.The Application of InSAR Technology to Mining Area Subsidence Monitoring［J］.Remote Sensing for Land and Resources，2008(2):51 －55.

［21］ Glodstein R M，Zebker H A，Werner C L.Satellite Radar Interferometry:Two － dimensional Phase Unwrapping［J］.Radio Science，1988，23(4):713 －720.

［22］桑會(huì)勇，郭華東，韓春明，等.一種基于梯度信息的小波SAR圖像濾波方法［J］.測(cè)繪通報(bào)，2005(2):17－21.

Research on an Improved Integrated Filtering Algorithm of SAR Interferogram

SONG Rui1，LIU Guang2，PERSKI Zbigniew3，GUO Hua －dong2
(1.School of Instrumentation Science and Opto－ electronics Engineering，Beijing University of Aeronautics＆ Astronautics，Beijing 100191，China;2.Center for Earth Observation and Digital Earth，Chinese Academy of Sciences，Beijing 100094，China;3.Polish Geological Institute－National Research Institute，Carpathian Branch 31560，Poland)

The interferogram with SAR phase information is the key factor that directly affects the further analysis of the deformation in the experimental area.Effective filter algorithm can suppress phase noise in the interferogram and improve the precision of the interferometry.In the obtained phase images，the characteristics of nearly annular stripe show the surface subsidence caused by coal mining.The authors therefore improve the original gradientbased filter algorithm firstly，and then combine Goldstein filter with the improved gradient－ based adaptive filter，making it perform more effectively in SAR interference phase mode from mining subsidence.In the experiment the authors selected the PALSAR interferogram data of the Fengfeng coal mine in Hebei Province to evaluate the detailed performance of such filtering methods.Experimental results show that the edge of deformation stripe in phase image is clearer，suggesting the validity of the integrated filtering algorithm which is characterized by good denoising effect and nice preservation of phase resolution.

SAR interferogram;Goldstein filter;improved gradient－based adaptive filter;mining area subsidence

TP 79

A

1001－070X(2012)01－0028－08

10.6046/gtzyyg.2012.01.06

2011－06－02;

2011－07－15

科技部中國(guó)與波蘭政府間科技合作項(xiàng)目(編號(hào):3417)、國(guó)家自然科學(xué)基金資助項(xiàng)目(編號(hào):41001264)以及中國(guó)科學(xué)院對(duì)地觀測(cè)與數(shù)字地球科學(xué)中心數(shù)字地球科學(xué)平臺(tái)重大項(xiàng)目(編號(hào):DESP01－04－10)共同資助。

宋 瑞(1990－)，男，主要研究方向?yàn)槔走_(dá)遙感。E－mail:songruibarry@qq.com。

(責(zé)任編輯:李 瑜)