一種適用于低頻InSAR的改進空變維納預濾波算法

夏路福，馮 東，李建鵬，安道祥，周智敏

（國防科技大學電子科學學院，湖南長沙 410073）

0 引 言

干涉合成孔徑雷達是一項綜合了合成孔徑成像技術和干涉測量的遙感技術［1］，它通過對不同的下視角或者時間觀測同一塊地區獲取的兩幅或多幅SAR 復圖像進行干涉處理，得到數字高程模型（Digital Elevation Model, DEM）［2］。InSAR 系統根據工作頻段不同分為高頻InSAR 系統和低頻InSAR 系統，高頻InSAR 系統發展較早，相關技術較為成熟［3］；而低頻InSAR 系統發展較晚，相關技術的研究有所欠缺［4］，由于其獨有的性質，近年來受到研究人員的關注［5］。

低頻InSAR 系統工作在VHF∕UHF∕L 等低頻波段［6］，其信號能夠穿透林葉，該系統可以完成對葉簇覆蓋區域的地形精確測繪，在民用和軍用領域用途廣泛。然而由于信號波長較長，可以穿透植被或者地表淺層等，導致所接收到的回波中雜波較多，噪聲較大，干涉處理過程中主輔圖像間相干性較差，難以獲得高質量的干涉相位圖。而高相干性的SAR復圖像對是InSAR 干涉測量的前提［7］，在InSAR 處理流程中，可以通過預濾波消除一些去相干因素的影響［8］，提高圖像間的相干性，改善干涉圖的質量［9］。

預濾波算法研究中，幾乎都是基于高頻InSAR數據開展研究，最常用的是公共帶預濾波算法［10］，該算法認為去相干引起了頻譜的偏移，圖像頻譜中公共部分具有相干性，非公共部分被認為是噪聲。公共帶預濾波算法在圖像頻域上進行處理，可以提高圖像對之間的相干性，但提升效果較小，且由于一部分頻譜被削減，導致圖像分辨率降低，進而影響干涉相位圖的質量，后續對該算法的研究與改進也沒有改變頻譜削減的事實［11-13］。Fornaro 等提出了空變維納預濾波算法［14］，它與公共帶預濾波算法不同，在圖像域上進行處理，不會削減頻譜。該算法對相干性的提升較大，但也存在明顯的缺點［15］，在大尺寸圖像上處理時計算量較大，耗時較長，且算法需要未知的干涉相位值作為已知條件，限制了算法的應用。

關于低頻InSAR 預濾波算法方面的研究很少，為了得到高質量的干涉相位圖［16］，本文針對低頻InSAR 的回波特點，結合公共帶預濾波理論，提出一種適合低頻InSAR 的公共帶預濾波算法；同時針對空變維納預濾波算法的缺點，簡化算法結構，優化參數估計，并結合公共帶預濾波的原理提出一種改進的空變維納預濾波算法。通過仿真和實測數據驗證了改進算法優于傳統算法，可以顯著提升圖像的相干性，改善干涉相位圖的質量。

1 低頻InSAR公共帶預濾波算法

1.1 公共帶預濾波算法

公共帶預濾波算法處理流程可以分為兩步，分別為距離向預濾波處理和方位向預濾波處理，其算法原理如圖1所示。

圖1 公共帶預濾波算法原理示意圖

圖1（a）為距離向公共帶預濾波算法原理示意圖，其中fc為系統載頻，B為系統帶寬，θ1和θ2為兩次飛行過程中的入射角。由于基線的存在，造成入射角的差異，導致主輔圖像地距頻譜產生偏移，通過保留公共地距頻譜部分，濾除非公共部分的方式完成距離向預濾波。

圖1（b）表示方位向公共帶預濾波算法原理示意圖，其中fη1、fη2分別為主輔圖像的多普勒中心頻率，Ba1和Ba2分別為主輔圖像的多普勒帶寬，圖像方位向頻譜與雷達斜視角和速度有關，公共帶預濾波算法濾除主輔圖像方位向頻譜中非公共區域，留下公共區域。

1.2 適用于低頻InSAR的公共帶預濾波算法

以條帶模式，正側視成像為例，SAR 系統接收到的回波二維頻譜范圍可以表示為［17］

式中，fr為距離向頻率，B為系統帶寬，fa為方位向頻率，va為雷達飛行速度，fc為系統載頻，c為光速，Δθ為方位向積累角。

在高頻InSAR 系統中，由于相對帶寬小，波束角小，圖像的二維頻譜近似為矩形，且方位向頻率范圍和距離向頻率近似不相關，公共帶預濾波時可以分開對距離向和方位向進行處理。而低頻In-SAR 系統相對帶寬較大，方位向波束角寬度大，其圖像的二維頻譜為梯形，方位向頻譜范圍與載頻以及距離向頻率有關。不同InSAR 系統的回波二維頻譜如圖2所示。

圖2 不同InSAR系統回波的二維頻譜示意圖

由于低頻InSAR 的回波方位向頻率與距離向頻率相關，對其進行預濾波處理時，距離向與方位向不能分開處理，本文結合低頻InSAR 系統回波的特點提出一種適用于低頻InSAR 的公共帶預濾波算法，該算法原理依舊是保留頻譜的公共部分，濾除非公共部分，與高頻InSAR 公共帶預濾波算法的區別在于其直接對圖像的二維頻譜處理，不再分開對距離向頻譜和方位向頻譜進行處理。

2 適用于低頻InSAR的改進空變維納預濾波算法

2.1 空變維納預濾波算法

以一維距離向為例，經過成像、配準后的信號對可以表示為

式中s1、s2分別為主、輔圖像信號，r'為距離坐標軸，z(r)為地面散射系數，φ1(r)、φ2(r)為相位，f1、f2為脈沖響應函數，n1、n2為加性熱噪聲。

假設一條距離線的長度為2N+ 1，將式（2）轉化為數學矩陣形式：

以信號S1為例，則其中的參數可以表示為

假設地面散射系數、系統熱噪聲服從零均值的高斯分布，且三者互不相關。則它們的自相關矩陣可以表示為

式中I為單位矩陣，I2N+1表示(2N+ 1) ×(2N+ 1)的單位矩陣。根據正交化原理，可以推導出濾波器的矩陣形式，如下所示：

在式（11）和式（12）中，

式中?inf為干涉相位，可通過主輔圖像干涉處理獲得。

上述式（11）和式（12）中濾波器的構建，需要知道4個參數即可，分別為距離向雷達脈沖響應函數（f1和f2）、地面散射系數的自相關矩陣σ0I2N+1、噪聲功率、干涉相位?inf。其中干涉相位?inf可通過主輔圖像干涉處理獲得。距離向雷達脈沖響應函數一般為sinc型函數。地面散射系數的自相關矩陣可以通過信號的互相關求解，如下所示：

同理，噪聲功率可以根據信號的自相關求解，如下所示：

通過式（15）求解出參數σ0I2N+1后，代入式（16）中可以求解出求解步驟與求解步驟類似。

空變維納預濾波算法適用性廣，對相干系數提升效果較大，但在實際應用中會產生兩個問題：第一，使用該算法進行預濾波時是對一整條距離線處理，構建濾波器以及參數估計時會直接出現大型矩陣的相乘和求逆的情況，計算量較大，計算效率不高，限制了算法的使用；第二，濾波器中干涉相位直接來自未進行處理的干涉相位，其所含噪聲較多，干涉相位不準確，最終會影響預濾波效果。

2.2 改進的空變維納預濾波算法

針對空變維納預濾波算法出現的問題，本文對該算法的結構進行變化，采用單個像素點進行濾波，可以使濾波器結構從(2N+ 1) ×(2N+ 1)矩陣轉變為1 × 1 的矩陣，參數估計也隨之發生變化，這樣降低了運算的復雜度。

信號對的表達式依舊如式（2）所示，將其轉變為矩陣形式：

在式（17）中，以信號S11中各參數為例，可知

上述式子中i表示一條距離線上的某一點，即某個像素點。最終濾波器的結構如下所示：

式中I1為1 × 1 的單位矩陣。式（21）和式（22）的前半段、后半段都為1 × 1 的矩陣，計算量有所減小。同時參數的估計也作出調整，其中距離向雷達脈沖響應函數和干涉相位的求解沒有變化，地面散射系數可以直接通過SAR復圖像強度獲得：

式中A為SAR 復圖像強度。則地面散射系數的自相關矩陣σ0I1可直接通過式（9）求解。

對于噪聲功率，需要估計每個像素點的噪聲功率，上述估計方法不再適用，查閱文獻可知［18］

針對第二個問題，空變維納預濾波算法需要干涉相位作為輸入，一般直接采用未進行預濾波處理的主輔圖像共軛相乘得到的干涉相位，這樣處理會引入含有較多噪聲的干涉相位，進而影響濾波器的準確性。而經過預濾波處理后的干涉相位所含噪聲減小，因此本文引入迭代操作，以上一次像素點空變維納預濾波處理后的干涉相位、主輔圖像作為下一次像素點空變維納預濾波的輸入，經過一次預濾波處理后得到的干涉相位所包含的噪聲較少，干涉相位更加準確，避免了使用含噪聲較高的數據，提高了預濾波效果。

最后結合適用于低頻InSAR 的公共帶預濾波算法，保留公共頻譜，濾除非公共頻譜可以進一步去除噪聲，提升相干性。在不影響圖像分辨率的情況下，選擇合適的頻譜范圍作為公共部分，其中頻譜公共區域范圍可以表示為

式中F2d1為主圖像的二維頻譜范圍，F2d2為輔圖像的二維頻譜范圍。經過公共帶濾波處理后，公共區域的頻率范圍包含了回波的信號頻率，保證圖像的分辨率不會受損，同時去除非公共帶中不相干的噪聲，可以進一步提升預濾波效果，提高相干性。本文算法的算法流程圖如圖3所示。

圖3 本文算法的算法流程圖

其具體算法步驟為：

步驟1：分別估計主輔圖像的噪聲功率和IRF矩陣。

步驟2：計算干涉相位，根據公式估計地面散射系數自相關矩陣。

步驟3：構建濾波器矩陣形式并進行像素點空變維納預濾波。

步驟4：將得到的主輔圖像作為輸入，重復步驟1～3。

步驟5：根據式（26）以及雷達和場景參數求解合適的公共頻譜范圍。

步驟6：對主輔圖像的二維頻譜分別進行公共帶濾波，得到最終的結果。

從圖3可以看出，所提算法為迭代像素點空變維納預濾波算法與公共帶預濾波算法的結合，其中迭代像素點空變維納預濾波算法與輸入信號的頻率無關，而公共帶預濾波算法的處理與信號頻率相關，當采用高頻信號時，將公共帶預濾波算法變化為適用高頻數據的即可，即分別進行距離向和方位向預濾波處理。

3 仿真數據驗證

為驗證本文所提算法的有效性，本小節使用仿真數據驗證，采用Matlab 軟件進行信號仿真，仿真參數如表1 所示，仿真場景如圖4 所示，為一個位于平面中心的理想圓錐體。

表1 仿真參數表

圖4 仿真場景

經過回波信號的仿真后，采用后向投影（Back-Projection,BP）成像算法處理，使用地距成像網格，網格參考高度為0 m。隨后對主、輔圖像添加高斯白噪聲，經過配準得到處理后的主輔圖像，如圖5所示，它們的二維頻譜如圖6所示。

圖5 經過配準后的主、輔圖像對

圖6 仿真數據二維頻譜圖

對主輔圖像對分別采用不同的預濾波算法處理，其中公共帶預濾波算法中公共帶區域如圖7所示，本文所提算法中噪聲功率估計數值如圖8所示。

圖7 公共帶區域

圖8 噪聲功率估計數值

經過預濾波處理后得到實驗結果，其中圖9所示為經過不同預濾波算法處理后得到的相干系數分布圖，圖10 所示為經過預濾波算法處理后得到的干涉相位圖。為了觀察方便，干涉相位圖經過了去平地相位處理。

圖9 經過預濾波算法處理后得到的相干系數分布圖

圖10 經過預濾波算法處理后得到的干涉相位圖

分析對比圖9 中的4 幅圖可以發現，經過預濾波處理后，主輔圖像間的相干性得到提升，相對于平地區域，圓錐體區域的相干性提升效果較弱，主要由于陡峭的地形與坡度一定程度上會影響預濾波的效果，圖9（d）為經過本文所提算法處理后得到的相干系數分布圖，可以看出相干系數提升最大，不管是平地區域還是圓錐體區域，相干性的提升十分明顯。

觀察圖10 可以發現，不同預濾波算法處理后干涉相位圖的條紋都有所改善，對比圓錐體區域的右側條紋可以發現經過公共帶預濾波和空變維納預濾波處理后，條紋上面的噪點依舊較多，條紋的清晰度提升不明顯；而經過本文所提算法處理后，圖像整體變得清晰，相位圖上的噪聲水平明顯比其他算法的低，圓錐體區域的條紋很明顯。整體上看，本文所提算法的處理結果目視效果最優，條紋最清晰。

表2 為經過不同預濾波算法處理后的數值結果，觀察平均相干系數這一欄可以發現經過公共帶預濾波、空變維納預濾波、本文所提算法處理后相干系數分別提高了0.087 3，0.116 1，0.194 3，客觀數據上可以看出本文所提算法對主輔圖像間相干系數的提升最大。對比算法的耗時，空變維納預濾波算法耗時為9.71 s，而本文所提算法的耗時只有5.47 s，可以證明本文所提算法縮短了耗時，提升了算法效率，具有一定的可行性。

表2 不同預濾波算法處理結果

圖11 為相干系數分布的統計直方圖，不同預濾波算法處理后的相干系數分布統計曲線分別用不同顏色表示，對比未預濾波處理的曲線，其他曲線明顯后移，圖像的相干性整體上得到了提升；與其他預濾波算法對比，本文所提算法處理后得到相干系數大部分分布在0.5～0.8 之間，明顯高于其他算法的結果，顯著提升了整幅圖像的平均相干系數水平。

圖11 相干系數分布的統計直方圖

4 實測數據驗證

為了驗證所提算法的適用性，本小節采用國防科技大學在2014年通過外場飛行實驗獲取的機載重軌P 波段InSAR 實測數據，雷達平臺高度為3 800 m，基線長度約100 m，雷達離場景中心的最近地距為9 000 m。通過BP算法成像，成像網格為地距網格，選擇一個2 305×4 113 像素的典型區域，其中橫軸代表地距向，縱軸代表方位向。經過成像與配準后的圖像如圖12 所示，主輔圖像的二維頻譜圖如圖13所示。

圖12 實測數據圖像

圖13 實測數據二維頻譜圖

對主輔圖像進行預濾波處理，其中公共帶預濾波算法中公共帶區域如圖14 所示，本文所提算法中噪聲功率估計數值如圖15所示。

圖14 公共帶區域

圖15 噪聲功率估計數值

經過預濾波后的實驗結果如圖16、圖17所示，其中圖16 為相干系數圖，圖17 為干涉相位圖，為了方便觀察，干涉相位圖經過了去平地相位處理。

圖16 經過預濾波算法處理后得到的相干系數圖

圖17 經過預濾波算法處理后得到的干涉相位圖

從圖16 可以看出，經過公共帶預濾波算法處理后相干系數分布圖的顏色沒有明顯變化，相干系數提升不大，而經過本文所提算法處理后，相干系數分布圖的顏色加深，相干系數提升較大。觀察藍色的低相干區域可以看出，經過公共帶預濾波和空變維納預濾波處理后，提升效果不顯著，經過本文所提算法處理后，這些區域的相干性有了明顯提升，提升效果優于其他預濾波算法。

對比圖17 的4 幅圖可以發現，公共帶預濾波算法和空變維納預濾波算法處理后的干涉相位圖噪聲有所減少，但是條紋依舊不清晰，上面的噪聲變化不明顯，而經過本文所提算法處理后的干涉相位圖上的條紋更加明顯，圖像整體的噪聲水平下降，其中高噪聲水平區域的條紋仍然不清晰。

表3 所示經過不同預濾波算法處理后的數值結果，可以看出預濾波處理后圖像對之間的平均相干系數都得到提升，其中經過公共帶預濾波、空變維納預濾波處理后平均相干系數分別提高了0.023 8 和0.081 3，而本文所提算法對平均相干系數的提升達到了0.170 6，提升效果較大。對比空變維納預濾波算法的耗時，本文所提算法耗時大幅縮減，計算耗時在30 min 左右，計算效率較高，處理較大尺寸圖像有一定的實用性。

表3 不同預濾波算法處理結果

圖18表示經過不同預濾波算法處理后的相干系數分布的統計直方圖，可以看出經過本文所提算法處理后，相干系數分布曲線變化明顯，低相干系數的像素點數量大幅減少，高相干系數的像素點數量增多，相干系數統計曲線明顯后移，大部分像素點對應的相干系數在0.6～0.9 之間，相干系數整體有所提升。

圖18 相干系數分布的統計直方圖

5 結束語

本文針對低頻InSAR 成像特點提出了一種改進的空變維納預濾波算法，為了增強空變維納預濾波算法的計算效率以及預濾波效果，本文改進了濾波器的結構，改變了參數估計的方法，避免了大型矩陣的相乘和求逆，提升了算法效率；引入迭代處理，避免直接采用含噪聲較多的干涉相位構建濾波器的問題，從而提升了預濾波效果；最后聯合公共帶預濾波算法，在不損失圖像分辨率的情況下去除多余的噪聲，進一步提升了主輔圖像的相干性。后續通過仿真數據和實測數據的驗證，定性和定量分析表明，本文所提算法的性能優于其他預濾波算法，可以顯著提升主輔圖像間的相干性，改善干涉相位圖的質量，具有一定的實用性。