自適應偽相干值改進的InSAR干涉圖Goldstein濾波算法

李世金 張書畢 張秋昭 高延東

（中國礦業大學環境與測繪學院，江蘇徐州221116）

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

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

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

1 干涉圖濾波算法

1.1 Goldstein濾波

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

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

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

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

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

式中?i為第i個像素的復數相位值；n為像素點數量。

1.2 改進后的濾波算法

為進一步提高Goldstein濾波算法的自適應性，本研究對獲取的偽相干值進行進一步優化，首先對大于閾值T的偽相干值進一步增大，對小于該閾值的偽相干值進一步減小；然后將當前濾波窗口中改進后的偽相干值均值作為濾波參數，并根據設定的閾值自適應確定平滑窗口尺寸。為此，構建的改進后的偽相干系數的取值模型為

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

當Pc≥T時（0＜β＜1），模型為

式（5）、式（6）中T值選取依據為：當 -Pc較大時，即干涉圖質量整體較好，濾波重點應為提升相位細節信息，故T值應相對較小；當-Pc較小時，濾波重點應為提升去噪效果，故T值應相對較大。本研究經過仿真試驗，構建的λ取值模型為

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

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

2 算法試驗

2.1 TerraSAR數據試驗

本研究采用TerraSAR數據進行濾波分析，結果見圖1、表1。

注：原始含噪聲干涉圖中殘差點有26 137個。

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

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

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

2.2 DEM反演分析

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

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

對不同濾波算法處理后的干涉圖運用GAMMA軟件處理后獲取的地距結構干涉高程圖如圖4所示。依據圖4提取的高程信息與SRTM高程數據的RMSE值如表3所示。

注：原始含噪聲干涉圖的殘差點數量有33 096個。

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

3 結語

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