一種基于相位補償?shù)母缮鎴D自適應濾波算法

譚宏偉，周芳，賀豐收

(1.合肥工業(yè)大學 計算機與信息學院，合肥 230009；2.工業(yè)安全與應急技術(shù)安徽省重點實驗室，合肥 230601；3.中國航空工業(yè)集團公司雷華電子技術(shù)研究所，江蘇 無錫 214063)

0 引言

合成孔徑雷達干涉測量技術(shù)[1-2](interferometric synthetic aperture radar，InSAR)利用覆蓋同一地面區(qū)域的兩幅合成孔徑雷達圖像生成干涉圖，進而反演獲得地面的高程或形變。InSAR具有觀測范圍大、測量精度高、可以全天時全天候?qū)Φ孛孢M行觀測的優(yōu)點[3-4]，近年來受到人們的廣泛關(guān)注。然而，受各種去相干因素的影響，InSAR生成的干涉圖相位中往往含有大量噪聲，嚴重影響相位展開、高程重建等后續(xù)處理的性能[5-6]，因此，在干涉測量過程中需要對干涉圖進行濾波。

目前，研究者已經(jīng)提出多種干涉圖濾波算法，這些算法大致可以分為兩類：空間域濾波和變換域濾波?？臻g域濾波中多視濾波[7]是應用最早的一種算法，它通過將均值濾波模板在干涉圖上逐像素移動卷積，從而獲得每一像素的濾波結(jié)果。多視濾波的噪聲抑制效果是以犧牲圖像分辨率為代價的，不利于相位細節(jié)的保持。在多視濾波基礎(chǔ)上，Trouve等[8]提出坡度相位補償濾波算法，通過補償濾波窗口內(nèi)的線性坡度相位，對殘余相位進行窗口大小自適應的多視濾波，從而能更好地去除噪聲并降低相位信息損失。然而受線性相位模型的限制，當局部區(qū)域地形復雜或相位噪聲嚴重時，坡度相位補償濾波的濾波效果不佳。因此，Wang等[9]提出主相位成分補償濾波算法，它通過對局部干涉圖的頻譜進行加權(quán)截斷處理來獲得反映地形輪廓的主要相位成分，然后對補償了主要相位成分的殘余相位進行自適應方向窗口濾波。主相位成分補償濾波不受線性相位模型的約束，在一定程度上提高了濾波的效果，但該算法提取的主要相位分量存在過量噪聲殘留的情況，并且殘余相位濾波時濾波像元數(shù)目固定，因此存在過濾波或欠濾波的情況。比較常用的空間域濾波算法還包括Lee濾波[10-12]、非局部均值濾波[13-14]等，這些算法在濾波處理時均存在一些問題，如Lee濾波中方向模板窗口不能適應所有的相位條紋模式，在干涉圖相位變化復雜的區(qū)域會破壞條紋結(jié)構(gòu)；非局部均值濾波計算復雜度高，且在相位噪聲嚴重的區(qū)域難以有效估計像素間的相似性。

Goldstein濾波[15]是變換域濾波算法中最為經(jīng)典的一種算法，它通過對局部干涉圖的幅度譜進行平滑加權(quán)處理從而得到濾波器的響應。Goldstein濾波能有效保持相位的條紋結(jié)構(gòu)，然而其噪聲去除效果受濾波參數(shù)的影響較大，當參數(shù)設(shè)置不當時，干涉圖會出現(xiàn)過度平滑或者大量噪聲殘留的情況。此后，學者們針對Goldstein濾波的不足提出一系列改進算法，如Baran等[16]基于干涉圖局部相干性均值設(shè)置頻譜加權(quán)參數(shù)的大小，Suo等[17]結(jié)合坡度相位補償以及自適應頻譜平滑窗口大小和頻譜加權(quán)參數(shù)大小進行濾波。這些算法的濾波效果雖然較Goldstein濾波有所提升，然而均難以說明參數(shù)設(shè)置的最優(yōu)性。干涉圖小波濾波[18-21]能對干涉圖像進行多分辨率分析，因此能有效避免分辨率損失。然而，小波變換不能精確地表示圖像中的邊緣、條紋信息[22]，因此不利于條紋結(jié)構(gòu)的保持。

針對現(xiàn)有干涉圖濾波算法存在噪聲去除不充分或相位細節(jié)保持困難的問題，本文提出一種更加有效的自適應濾波算法。該算法通過對干涉圖預濾波并對預濾波干涉圖進行頻譜截斷的Goldstein平滑處理，從而提取更加精確的地形輪廓相位。在對補償了地形輪廓相位的殘余相位濾波時，基于像元相干系數(shù)設(shè)計了一種自適應像元數(shù)目的多視濾波方法，在有效去除噪聲的同時充分保持相位細節(jié)。

1 基礎(chǔ)理論分析

1.1 相位噪聲統(tǒng)計特性

干涉圖相位噪聲的嚴重程度可以由生成它的兩幅SAR圖像的相干系數(shù)來衡量，相干系數(shù)越大，相位噪聲越小。Franceschetti等[23]給出了相位噪聲標準差與相干系數(shù)以及多視視數(shù)的關(guān)系，如式(1)所示。

(1)

式中：σφ為相位噪聲標準差；γ為像元相干系數(shù)；L為多視視數(shù)；φ為干涉相位；φ0為相位的數(shù)學期望；ρ(φ，γ，L，φ0)代表相位的概率密度函數(shù)，其解析表達如式(2)所示。

(2)

式中：Γ(g)表示伽瑪函數(shù)；β=γcos(φ-φ0)；F(g)表示高斯超幾何函數(shù)。

結(jié)合式(1)和式(2)，可以計算不同相干系數(shù)以及多視視數(shù)情況下的相位噪聲標準差。當視數(shù)為1、2、5、10、15、20時，噪聲標準差隨相干系數(shù)的變化曲線如圖1所示。從圖1可以看出，相位噪聲標準差隨著相干系數(shù)的增大而逐漸減小，并且增大多視視數(shù)可以有效降低噪聲。

圖1 相位噪聲標準差與相干系數(shù)以及多視視數(shù)的關(guān)系

研究表明，當多視視數(shù)大于4且相干系數(shù)大于0.2時，相位噪聲標準差的近似表示如式(3)所示[24]。

(3)

1.2 相位補償濾波分析

相位補償濾波中的干涉圖相位模型如式(4)所示。

φz=φm+φr+v

(4)

式中：φz為測量相位；φm為地形變化引起的相位或反應地形輪廓的相位；φr為殘余信號相位；v為噪聲相位，v均值為零且與φm和φr相互獨立。相位補償濾波算法通過在局部窗口內(nèi)估計并補償?shù)匦巫兓辔换虻匦屋喞辔?，然后對殘余相位進行更加有效的濾波處理，從而能有效增強濾波的去噪效果并減少相位信息損失。

坡度相位補償濾波假設(shè)局部窗口內(nèi)的相位近似滿足線性模型，即式(4)中地形變化相位φm可以表示為式(5)。

φm(x，y)=2π(xfx+yfy)

(5)

式中：(x，y)為像素在窗口內(nèi)的位置坐標；fx和fy分別為方位向和距離向的相位頻率。通過估計頻率fx和fy可以得到線性坡度相位，進而對殘余相位進行窗口大小自適應的多視濾波。

然而在坡度相位補償濾波中當局部區(qū)域地形復雜時，干涉相位在濾波窗口內(nèi)難以用線性模型近似，濾波會存在細節(jié)損失嚴重的問題。此外，當相位噪聲較大時，相位頻率fx和fy的估計精度會降低，濾波可能引入虛假相位條紋。殘余相位濾波時，由于不同大小的窗口內(nèi)像元數(shù)目不連續(xù)變化，因此坡度相位補償濾波不能實現(xiàn)濾波強度自適應。在此基礎(chǔ)上，文獻[9]提出了主相位成分補償濾波算法，算法在更大的局部窗口內(nèi)通過對干涉圖的頻譜進行加權(quán)截斷處理，獲得反映地形輪廓的主要相位成分，避免了線性模型近似，并對殘余相位進行自適應方向窗口濾波。主相位成分補償濾波中主要相位分量的計算方法如式(6)所示。

φm(x，y)=arg{FFT2-1{S(u，v)Appc(u，v)}}

(6)

式中：arg{g}為求復數(shù)相位；FFT2-1{g}為二維逆傅里葉變換；(u，v)表示頻率坐標；S(u，v)代表局部窗口內(nèi)的干涉圖頻譜；Appc(u，v)代表頻譜截斷后的干涉圖幅度譜。在文獻[9]中，頻譜截斷的閾值設(shè)為半功率點。

然而，由于噪聲分布在干涉圖頻譜的所有頻帶上[25]，因此主相位成分補償濾波算法提取的主要相位分量存在過量噪聲殘留的情況。并且在實際數(shù)據(jù)處理時，頻譜截斷閾值設(shè)為半功率點時提取的主要相位分量往往無法反映干涉相位的輪廓結(jié)構(gòu)，而且殘余相位濾波時沒有進行自適應處理，因此濾波存在噪聲去除不足或相位細節(jié)損失嚴重的問題。

2 改進的相位補償濾波算法

針對干涉圖相位補償濾波中存在的上述問題，本文提出一種改進算法，以進一步提高濾波的性能。提出算法將濾波過程分為以下四步。

1)對干涉圖預濾波，降低相位噪聲對干涉圖頻譜的影響。

2)將預濾波干涉圖劃分為一系列重疊的圖像塊，并采用頻譜截斷的Goldstein濾波方法提取每一圖像塊的地形輪廓相位。

3)從相應位置的噪聲干涉圖中去除圖像塊的地形輪廓相位，然后基于像元相干系數(shù)對殘余相位進行自適應濾波。

4)將步驟3)中的濾波后相位與步驟2)中提取的地形輪廓相位相加，得到每一圖像塊的濾波結(jié)果，并對圖像塊的邊緣重疊部分進行加權(quán)平滑處理。

算法流程如圖2所示。

圖2 濾波算法流程圖

2.1 干涉圖預濾波

噪聲會影響地形輪廓相位的提取精度，因此有必要對干涉圖預濾波。本文采用坡度相位補償濾波與多視濾波相結(jié)合的方法進行預濾波處理。

坡度相位補償濾波可以表示為式(7)。

(7)

(8)

式中：max{g}表示求函數(shù)最大值；P、Q為估計窗口的方位向和距離向像素半徑；j為虛部；|g|表示求復數(shù)的模。

由1.2節(jié)分析可知，坡度相位補償濾波在干涉圖噪聲嚴重的區(qū)域濾波效果較差，嚴重時可能因為相位頻率的錯誤估計而引入虛假相位條紋。因此本文在低相干區(qū)域采用多視濾波方法進行預濾波處理，多視濾波計算方法如式(9)所示。

(9)

預濾波在初步去除噪聲的同時應該避免過度模糊相位，因此本文設(shè)置較小的濾波窗口，如坡度相位補償濾波的濾波窗口大小設(shè)為5像素×5像素，多視濾波的濾波窗口大小設(shè)為3像素×3像素，以充分保留相位信息。

2.2 局部區(qū)域地形輪廓相位提取

局部干涉圖頻譜中地形相位分量主要分布在較窄的頻帶范圍內(nèi)，而噪聲分量則分布在所有頻帶上[25]，本文利用干涉圖頻譜的這一特點從預濾波干涉圖中提取地形輪廓相位。

首先，將預濾波干涉圖劃分為一系列一定大小的圖像塊，并且為了避免濾波后相鄰圖像塊間相位條紋不連續(xù)，本文采用與Goldstein濾波方法相同的策略，設(shè)置圖像塊間有一定的重疊，在濾波完成后，對重疊部分的濾波干涉圖加權(quán)平滑處理。

對每一塊預濾波干涉圖進行二維傅里葉變換，得到其頻譜S(u，v)。由于預濾波初步降低了相位噪聲，因此信號頻譜中混疊的噪聲分量降低，更有利于地形輪廓相位的提取。為了進一步避免提取的相位中包含噪聲，本文對預濾波干涉圖進行頻譜截斷的Goldstein濾波處理，即只保留頻譜中幅度超過閾值的窄帶信號分量并對其加權(quán)平滑以進一步降低噪聲的影響。提取的輪廓相位如式(10)所示。

φ′m(k，l)=arg{FFT2-1{SM(|S(u，v)|)T(S(u，v))}}

(10)

式中：SM(g)表示幅度譜平滑，通常采用3像素×3像素大小的均值窗口進行平滑處理；T(g)表示頻譜閾值截斷，如式(11)所示。

(11)

式中：b為頻譜截斷閾值。經(jīng)過預濾波處理，地形相位分量的窄帶特性進一步增強。該頻譜會在較少頻率上形成明顯的峰值，因此本文將b設(shè)為幅度譜均值的3倍，即b=mean(|S(u，v)|)×3，經(jīng)過測試提取的地形輪廓相位較為平滑且能精確地反映干涉相位的結(jié)構(gòu)。

2.3 殘余相位自適應濾波

將地形輪廓相位從噪聲干涉圖中去除，得到殘余相位，如式(12)所示。

φR(k，l)=φr(k，l)+v(k，l)=arg{φ(k，l)e-jφ′m(k，l)}

(12)

由于殘余信號相位φr主要是干涉相位的細節(jié)分量，因此殘余相位濾波應能在有效去除噪聲的同時充分避免細節(jié)損失。本文基于像元相干系數(shù)設(shè)計了一種自適應像元數(shù)目的多視濾波方法，避免濾波在低相干區(qū)域噪聲去除不足或者在高相干區(qū)域相位細節(jié)損失嚴重。

由1.1節(jié)的分析可知，當相干系數(shù)大于0.2且多視視數(shù)大于4時，相位噪聲標準差與相干系數(shù)以及多視視數(shù)間的關(guān)系可以表示為式(3)。因此，可以綜合考慮噪聲去除以及細節(jié)保持兩方面因素，設(shè)定濾波后的噪聲標準差值，根據(jù)此設(shè)定值就可以得到不同相干系數(shù)的像元應該采用的視數(shù)。另一方面，當像元相干系數(shù)很低時，需要較大視數(shù)才能將噪聲標準差降至期望值，然而視數(shù)太大可能使相位損失嚴重，因此本文定義最大視數(shù)值，當視數(shù)達到最大視數(shù)時不再增加。當像元相干系數(shù)較高時，得到的視數(shù)可能較小，為了充分降低噪聲，本文定義最小視數(shù)值，當視數(shù)降至最小視數(shù)時也不再減小。由此，可以得到多視視數(shù)，如式(13)所示。

(13)

式中：Nmax、Nmin、σ分別為設(shè)定的最大視數(shù)、最小視數(shù)以及濾波后的相位噪聲標準差；round(g)代表四舍五入。

由于已經(jīng)將輪廓相位從干涉圖中去除，殘余相位基本不含有相位趨勢，因此可以從濾波像元鄰域內(nèi)選擇相應視數(shù)的濾波樣本進行多視處理。殘余相位濾波的結(jié)果如式(14)所示。

(14)

式中：x為選擇的濾波樣本位置。

2.4 濾波結(jié)果

將提取的地形輪廓相位與濾波后的殘余信號相位相加，即可得到濾波結(jié)果，如式(15)所示。

φ(k，l)=wrap{φ′m(k，l)+φ′r(k，l)}

(15)

式中：wrap{g}表示相位纏繞操作。

3 實驗結(jié)果與分析

為驗證提出算法的有效性，分別采用仿真干涉圖和實測干涉圖進行濾波實驗，并將提出算法的濾波結(jié)果與Lee濾波、坡度相位補償濾波、Goldstein濾波、文獻[9]中的主相位成分補償濾波以及文獻[14]中的非局部均值濾波的濾波結(jié)果進行對比分析。

3.1 仿真干涉圖濾波

在MATLAB中仿真兩種具有不同相位變化特點的無噪干涉圖，圖像大小均為300像素×300像素，干涉相位如圖3(a)和圖3(c)所示。其中圖3(a)的相位頻率呈線性變化，頻率從上到下由1/28逐漸增大到1/8，而圖3(c)的條紋結(jié)構(gòu)較為復雜，相位非線性變化嚴重。為驗證算法在不同相干質(zhì)量條件下的濾波性能，在圖像上加入不同強度的相位噪聲，使干涉圖相干性從左到右由0.1增大到0.9，噪聲干涉相位分別如圖3(b)和圖3(d)所示。

圖3 仿真干涉相位圖

分別采用本文算法以及幾種對比算法對仿真的噪聲干涉圖進行處理，結(jié)果如圖4和圖5所示。在濾波參數(shù)設(shè)置時，本文將對比算法的參數(shù)大小設(shè)置為該類算法所采用的一般性大小，如Goldstein濾波、主相位成分補償濾波的圖像塊劃分尺寸設(shè)為32像素×32像素，重疊率設(shè)為50%，Goldstein濾波采用3像素×3像素大小的頻譜平滑窗口，頻譜加權(quán)參數(shù)α設(shè)為0.6。將本文算法與對比算法中具有相同含義的參數(shù)設(shè)置為一致，以觀察相同參數(shù)下的濾波效果。對于本文算法獨有的參數(shù)，本文則是參考相關(guān)文獻并結(jié)合濾波實驗進行的設(shè)置，如濾波后噪聲標準差是參考文獻[26]中的相關(guān)說明將其設(shè)為0.2 rad，最大視數(shù)、最小視數(shù)則是參考坡度相位補償濾波中的最大窗口、最小窗口設(shè)置，將其分別設(shè)為81和9。

圖4 仿真干涉圖1濾波結(jié)果

圖5 仿真干涉圖2濾波結(jié)果

由圖4和圖5可以看出，各種濾波算法在仿真干涉圖的高相干區(qū)域濾波結(jié)果相近，然而隨著相干性的降低，濾波算法的性能出現(xiàn)較大差異。Lee濾波在低相干且相位變化較快的區(qū)域造成嚴重的條紋結(jié)構(gòu)損失，這主要是由于Lee濾波中方向窗口不能適應各種條紋模式，并且噪聲會影響方向窗口檢測的準確性。坡度相位補償濾波在中高相干性的圖像區(qū)域濾波性能較好，然而當相干性較低時去噪效果較差。Goldstein濾波以及非局部均值濾波均能有效保持相位結(jié)構(gòu)，然而其去噪性能較差，在圖像左側(cè)殘留有大量的相位噪聲。主相位成分補償濾波在仿真干涉圖1上濾波效果較好，然而在仿真干涉圖2上存在低相干區(qū)域去噪不充分，干涉條紋復雜區(qū)域的相位結(jié)構(gòu)損失嚴重的情況。綜合來看，本文算法在不同相干性條件下均能有效去除相位噪聲并保持相位的結(jié)構(gòu)信息不被破壞，取得最好的濾波效果。

采用殘差點減少率以及濾波相位相對于真實相位的均方根誤差來客觀評價濾波算法的性能，殘差點減少率越大說明算法的噪聲抑制效果越顯著，均方根誤差越小說明算法的相位保持效果越好。各種濾波算法的濾波性能分別如表1和表2所示，可以看出，本文算法在兩幅干涉圖上均具有最優(yōu)的濾波性能。

表1 不同算法對仿真干涉圖1濾波的效果評價

表2 不同算法對仿真干涉圖2濾波的效果評價

3.2 實測干涉圖濾波

實測干涉圖由覆蓋重慶地區(qū)的兩幅ERS-1/ERS-2衛(wèi)星數(shù)據(jù)生成，主輔圖像成像日期分別為1996年4月4日和1996年4月5日，對應成像區(qū)域地形起伏劇烈，干涉相位具有條紋密集且噪聲分布不平穩(wěn)的特點。濾波干涉圖大小為256像素×256像素，干涉相位如圖6(a)所示，其對應的相干系數(shù)如圖6(b)所示。分別采用本文算法以及五種對比算法對干涉圖進行濾波處理，參數(shù)設(shè)置情況與前文一致，結(jié)果如圖6(c)～圖6(h)所示。

圖6 實測干涉圖濾波結(jié)果

從圖6可以看出，幾種對比算法在實測干涉圖上的濾波表現(xiàn)與模擬干涉圖濾波較為相似，Lee濾波造成嚴重的條紋結(jié)構(gòu)損失；坡度相位補償濾波、Goldstein濾波以及非局部均值濾波均有效保持了相位的條紋結(jié)構(gòu)，然而它們在圖像低相干區(qū)域濾波不足；主相位成分補償濾波在一些區(qū)域改變了相位的結(jié)構(gòu)，并且在低相干區(qū)域噪聲去除不夠充分。從目視效果來看，本文算法的濾波結(jié)果平滑度最好，恢復出的條紋結(jié)構(gòu)也最清晰。

由于沒有無噪干涉圖作為參考，無法計算濾波相位的均方根誤差。各濾波算法的殘差點減少率分析如表3所示。從表3可知，本文算法的殘差點減少率最大，噪聲去除效果最好。

表3 不同算法對實測干涉圖濾波的效果評價

4 結(jié)束語

本文提出一種基于相位補償?shù)母缮鎴D自適應濾波算法。算法通過對干涉圖預濾波以及頻譜截斷的Goldstein平滑濾波，解決了相位補償濾波中地形輪廓相位提取不精確的問題，并且在殘余相位濾波時基于像元相干系數(shù)進行濾波強度的自適應設(shè)置，從而使算法在噪聲去除以及細節(jié)保持兩方面取得平衡。實驗結(jié)果表明，相比于五種對比算法，本文算法的殘差點減少率最大，相位信息損失最小，具有最優(yōu)的濾波性能。然而本文算法存在最優(yōu)參數(shù)選擇困難的問題，在實驗部分參數(shù)的設(shè)置是根據(jù)經(jīng)驗進行的，如何從理論上得到自適應的最優(yōu)參數(shù)仍然有待研究。