李世金 張書畢 張秋昭 高延?xùn)|

（中國礦業(yè)大學(xué)環(huán)境與測繪學(xué)院，江蘇徐州221116）

合成孔徑雷達(dá)干涉測量技術(shù)（Interferometric synthethic aperture radar，InSAR）是一種用于提取地表三維變化信息的空間對地觀測技術(shù)，是傳統(tǒng)合成孔徑雷達(dá)遙感技術(shù)與無線電磁波干涉測量技術(shù)的有效融合。該技術(shù)不僅具有成本低、近連續(xù)性的特點(diǎn)和遠(yuǎn)程遙感觀測的能力［1］，還能夠全天候、全天時、高效率地獲取大面積的地面高精度三維信息。因此，InSAR技術(shù)在近些年得到了快速發(fā)展，在礦區(qū)沉陷監(jiān)測［2］、城市沉降［3］、地震形變［4］、火山運(yùn)動［5］、山體滑坡［6-7］及冰川運(yùn)動［8］等方面也展現(xiàn)出了強(qiáng)大的技術(shù)優(yōu)勢。在InSAR技術(shù)的實(shí)際應(yīng)用中，相位解纏的結(jié)果將直接決定著地面高程信息的獲取精度，而相位解纏的結(jié)果又與干涉圖的質(zhì)量密切相關(guān)，因此，研究干涉圖的濾波算法對于進(jìn)一步提高干涉圖的質(zhì)量乃至InSAR技術(shù)的應(yīng)用效果大有裨益。

InSAR干涉圖濾波算法主要有空間域?yàn)V波及頻率域?yàn)V波2類。其中，空間域?yàn)V波最常用的是圓周期均值濾波與圓周期中值濾波，主要根據(jù)干涉圖的周期性特點(diǎn)，利用局部統(tǒng)計(jì)特性來實(shí)現(xiàn)濾波。在此基礎(chǔ)上，圓周期加權(quán)中值濾波算法被提出，該算法融合了上述2種算法的優(yōu)點(diǎn)，濾波效果有了一定程度提升。上述算法的濾波效果與濾波窗口的尺寸存在著密切的關(guān)系，隨著濾波窗口的增大，噪聲去除效果越來越好，但相應(yīng)的圖像分辨率損失嚴(yán)重。為此，郭交等［9］將基于局部區(qū)域增長算法的窗口自適應(yīng)獲取方法應(yīng)用到濾波算法中，通過選取均值像素的濾波樣本，以提高濾波的自適應(yīng)效果；易輝偉等［10］提出了改進(jìn)的梯度自適應(yīng)濾波算法，根據(jù)不同方向的相位梯度信息來獲取相應(yīng)像素點(diǎn)的加權(quán)系數(shù)，進(jìn)而實(shí)現(xiàn)自適應(yīng)濾波。頻率域?yàn)V波算法則是依據(jù)傅里葉變換原理將干涉圖從空間域轉(zhuǎn)換至頻率域，在頻域內(nèi)依據(jù)相位噪聲和信號的頻譜特性進(jìn)行濾波處理，極大提升了相位信息的保留能力。最經(jīng)典的濾波算法為Goldstein濾波［11］，該算法對具有一定重疊效果的相位塊在頻率域采用平滑濾波器處理后，并用固定的濾波參數(shù)對其功率譜進(jìn)行處理。鑒于濾波參數(shù)的選取具有一定的主觀性，Baran等［12］利用相干值均值代替濾波參數(shù)，使其能夠根據(jù)干涉圖的相干性實(shí)現(xiàn)自適應(yīng)濾波；Suo等［13］為了進(jìn)一步提升該算法相位噪聲的抑制能力，利用相干值的平方值代替相干值均值，在相同的濾波窗口下進(jìn)一步增強(qiáng)了濾波強(qiáng)度；為進(jìn)一步提升自適應(yīng)效果，Zhao等［14］利用偽相干值均值代替相干值均值作為濾波參數(shù)，并進(jìn)行迭代處理，在一定程度上避免了偏置樣本一致性對濾波參數(shù)的影響；Liu等［15］提出了分窗口模型濾波算法，將改進(jìn)后的經(jīng)驗(yàn)?zāi)Ｐ头纸猓‥mpirical mode decomposition，EMD）算法［16］與Goldstein濾波算法［17］作為不同窗口的濾波模型。

上述改進(jìn)濾波算法雖然在一定程度上提升了算法的濾波性能，但總體上無法有效兼顧相位噪聲抑制及相位信息保留兩方面。因此，本研究采用偽相干值對Goldstein算法的濾波參數(shù)進(jìn)行進(jìn)一步優(yōu)化，使其能夠根據(jù)不同閾值選取不同的自適應(yīng)模型，進(jìn)而自適應(yīng)確定最優(yōu)濾波參數(shù)，并根據(jù)閾值自適應(yīng)確定平滑窗口尺寸，對干涉圖的功率譜進(jìn)行平滑處理。

1 干涉圖濾波算法

1.1 Goldstein濾波

Goldstein濾波［11］首先將干涉圖劃分為若干個相互重疊的相位塊，然后依據(jù)傅里葉變換將其轉(zhuǎn)換至頻率域，對其功率譜進(jìn)行相應(yīng)的平滑處理：

式中，S{*}為平滑算子；F（u，v）和H（u，v）分別為濾波前后的頻率域干涉圖；（u，v）為頻率域中干涉圖中某像素點(diǎn)的坐標(biāo)值；α為濾波參數(shù)，α∈[0，1]。

當(dāng)α=0時，S{|F(u,v)|}=1，表明無濾波效果；當(dāng)α=1時，S{|F(u,v)|}表示由矩形平滑窗口構(gòu)成的低通濾波器，表明濾波效果非常強(qiáng)。此外，較大的相位塊尺寸和高濾波參數(shù)α值，有利于增強(qiáng)干涉圖的相干性；相位塊相互重疊（重疊度不宜小于75%）有助于減弱干涉圖邊界處的不連續(xù)性。

考慮到干涉圖中相位噪聲分布的不均勻性，如果對整個干涉圖采用固定的α值進(jìn)行濾波，易降低濾波的自適應(yīng)性，可能會在高相干地區(qū)存在過濾波，而在低相干地區(qū)存在欠濾波的現(xiàn)象，因此難以選擇合理的全局濾波參數(shù)α值。對此，Baran等［12］采用相干值的均值代替α值，使得干涉圖非相干性區(qū)域的濾波程度強(qiáng)于相干性區(qū)域：

Zhao等［14］利用偽相干值均值代替α值，相對于Baran濾波算法在一定程度上避免了偏置樣本一致性對α參數(shù)估計(jì)的影響，其偽相干值（Pc）可通過下式求取：

式中?i為第i個像素的復(fù)數(shù)相位值；n為像素點(diǎn)數(shù)量。

1.2 改進(jìn)后的濾波算法

為進(jìn)一步提高Goldstein濾波算法的自適應(yīng)性，本研究對獲取的偽相干值進(jìn)行進(jìn)一步優(yōu)化，首先對大于閾值T的偽相干值進(jìn)一步增大，對小于該閾值的偽相干值進(jìn)一步減??；然后將當(dāng)前濾波窗口中改進(jìn)后的偽相干值均值作為濾波參數(shù)，并根據(jù)設(shè)定的閾值自適應(yīng)確定平滑窗口尺寸。為此，構(gòu)建的改進(jìn)后的偽相干系數(shù)的取值模型為

為有效兼顧濾波算法在噪聲抑制及相位信息保留兩方面的效果，需要有效確定α值，本研究構(gòu)建了隨偽相干值（Pc）自適應(yīng)變化的α取值模型。

當(dāng)Pc≥T時（0＜β＜1），模型為

式（5）、式（6）中T值選取依據(jù)為：當(dāng) -Pc較大時，即干涉圖質(zhì)量整體較好，濾波重點(diǎn)應(yīng)為提升相位細(xì)節(jié)信息，故T值應(yīng)相對較??；當(dāng)-Pc較小時，濾波重點(diǎn)應(yīng)為提升去噪效果，故T值應(yīng)相對較大。本研究經(jīng)過仿真試驗(yàn)，構(gòu)建的λ取值模型為

為進(jìn)一步提升算法的自適應(yīng)濾波效果，依據(jù)上述確定的閾值進(jìn)一步自適應(yīng)確定平滑窗口尺寸。當(dāng)默認(rèn)平滑窗內(nèi)的偽相干值均值于小閾值時，對其平滑窗口尺寸進(jìn)行增大處理，直至窗口內(nèi)的偽相干值均值大于閾值為止，以提高算法在低相干區(qū)域的降噪能力。隨后采用二維傅里葉變換后的平滑窗口對獲取的相位塊的功率譜進(jìn)行平滑處理［13］，并運(yùn)用改進(jìn)后偽相干值進(jìn)一步確定濾波參數(shù)，最終獲取濾波后的干涉相位。改進(jìn)后的Goldstein濾波算法可表述如下

式中，|*|m×m為傅里葉變換后的平滑窗算子；m為窗口尺寸。

2 算法試驗(yàn)

2.1 TerraSAR數(shù)據(jù)試驗(yàn)

本研究采用TerraSAR數(shù)據(jù)進(jìn)行濾波分析，結(jié)果見圖1、表1。

注：原始含噪聲干涉圖中殘差點(diǎn)有26 137個。

分析圖1可知：在高相干性區(qū)域4種濾波算法的濾波效果基本相似，而在中部及右下角相干性較低的區(qū)域，本研究算法的去噪效果明顯優(yōu)于其余3種算法。

分析表1可知：4種濾波算法的噪聲抑制能力均較理想，其中Goldstein濾波與Zhao濾波處理后的干涉圖中剩余的殘差點(diǎn)數(shù)量較接近，相應(yīng)的噪聲去除率分別為80.74%、80.40%；Baran濾波效果略優(yōu)于前2種算法，噪聲抑制率約82.24%；本研究算法的噪聲抑制率達(dá)到88.64%，并且該算法的EPI明顯高于其余3種算法，表明該算法不僅濾波性能較好，而且具有良好的相位保持能力。

在上述分析的基礎(chǔ)上，采用GAMMA軟件中的枝切樹相位解纏模塊對不同濾波算法處理后的干涉圖分別開展了相位解纏工作，結(jié)果如圖2所示。分析圖2可知：由于殘余殘差點(diǎn)的影響，易導(dǎo)致解纏結(jié)果中出現(xiàn)大量的未解纏區(qū)域，經(jīng)過本研究算法濾波后的干涉圖的整體解纏效果明顯優(yōu)于Goldstein濾波、Baran濾波以及Zhao濾波，尤其是圖中標(biāo)記區(qū)域，未解纏區(qū)域明顯少于其余3種算法；此外，Zhao濾波算法處理后的干涉圖的解纏結(jié)果中出現(xiàn)了大量誤差傳遞現(xiàn)象（圖2（c）），本研究算法則無此現(xiàn)象，進(jìn)一步證明了該算法的有效性。

2.2 DEM反演分析

為進(jìn)一步驗(yàn)證本研究算法對于實(shí)測數(shù)據(jù)的處理效果，采用1組5 000×5 000的TeraaSAR/TanDEM數(shù)據(jù)進(jìn)行了DEM反演分析。首先對比分析不同濾波算法對于該數(shù)據(jù)的濾波效果，然后采用GAMMA軟件對其進(jìn)一步處理，對不同濾波算法處理后獲取的高程精度進(jìn)行比較分析。

本研究采用去地勢后的干涉圖進(jìn)行試驗(yàn)。分析圖3、表2可知：4種濾波算法對圖3（a）中標(biāo)定區(qū)域的濾波效果均較好，噪聲去除率依次為86.31%、86.70%、86.62%、90.51%；本研究算法濾波后的干涉圖中的殘余相位噪聲點(diǎn)明顯少于其余3種算法，并且該算法處理后的干涉圖的EPI指數(shù)明顯大于其余3類算法。

對不同濾波算法處理后的干涉圖運(yùn)用GAMMA軟件處理后獲取的地距結(jié)構(gòu)干涉高程圖如圖4所示。依據(jù)圖4提取的高程信息與SRTM高程數(shù)據(jù)的RMSE值如表3所示。

注：原始含噪聲干涉圖的殘差點(diǎn)數(shù)量有33 096個。

分析圖4及表3可知：經(jīng)過Goldstein濾波、Baran濾波及Zhao濾波算法處理后獲取的干涉高程圖均出現(xiàn)了大量空洞現(xiàn)象，主要是由于殘余的噪聲點(diǎn)影響了后續(xù)的枝切樹相位解纏效果，導(dǎo)致相位解纏后出現(xiàn)了“孤島”現(xiàn)象，最終影響了高程信息的提取精度；本研究算法處理后獲取的干涉高程圖未出現(xiàn)該現(xiàn)象，并且相應(yīng)的RMSE值明顯小于其余3種算法。

3 結(jié)語

采用偽相干值對Goldstein濾波算法的濾波參數(shù)及濾波窗口進(jìn)行了優(yōu)化，TerraSAR-X/TanDEM-X數(shù)據(jù)濾波分析以及DEM反演分析表明，改進(jìn)后的濾波算法在噪聲抑制及條紋相位信息保持方面優(yōu)于Zhao算法、Baran算法以及Goldstein算法，并且在高程信息提取精度方面也有明顯優(yōu)勢。