趙超英，王寶行

1. 長安大學地質工程與測繪學院，陜西 西安 710054； 2. 地理國情監(jiān)測國家測繪地理信息局工程技術研究中心，陜西 西安 710054

合成孔徑雷達干涉測量(synthetic aperture radar interferometry，InSAR)已廣泛應用于火山、冰川、地面沉降、地震和滑坡等多類地表形變的監(jiān)測中，成為現(xiàn)代地球學科研究的重要技術之一。InSAR結果的監(jiān)測精度直接依賴于干涉相位的精度，然而由于熱噪聲、配準、體散射等失相干因素導致InSAR干涉圖往往存在嚴重噪聲，因此需要對干涉圖降噪(或稱濾波)以獲得高精度的測量成果。在干涉圖濾波方面，國內外研究大致可分為空間域和頻率域兩類方法?？臻g域濾波包括：中值濾波[1]及非局部濾波方法[2]等。頻率濾波的方法是針對信號高低頻特性進行平滑處理，最常用的有基于傅里葉變換的Goldstein濾波[3]和基于正余弦變換的小波濾波[4]，前者應用較為廣，且發(fā)展出諸多改進算法[5-8]。

這些改進主要集中在Goldstein濾波參數(shù)的估計方面，即自適應的濾波參數(shù)的確定。文獻[5]采用相干值代替了經典的單一濾波參數(shù)；文獻[6]兼顧了視數(shù)和相位標準差對相干值的影響，自適應地估計Goldstein濾波參數(shù)；文獻[7]以一種不顧及強度信息的偽相干圖來確定濾波參數(shù)；文獻[8]利用同質點計算的無偏相干性并采用穩(wěn)健估計的非線性濾波參數(shù)來進行濾波的。但是在噪聲很強的低相干地區(qū)，即使最新發(fā)展的穩(wěn)健估計的濾波參數(shù)也難以得到令人滿意的結果。本文研究從相位時間域信息來進行噪聲的濾除。

為克服噪聲及失相干的影響，文獻[9]提出SqueeSAR算法，分析農田、植被區(qū)等分布式目標(distributed scatterers，DS)的統(tǒng)計分布特性，對相位進行優(yōu)化，通過提高信噪比獲取更多的點目標，獲取更為完整的形變信息。該方法首先采用KS(Kolmogorov-Smirno)提取的服從同一統(tǒng)計分布的同質像元點，然后采用最大似然估計的方法優(yōu)化相位值。但是由于SqueeSAR的計算過程不僅要對相干矩陣求逆前進行插入阻尼因子防止負的或過小的特征值，而且需要一個最小的非線性目標函數(shù)進行最優(yōu)化相位的求解，計算比較復雜。文獻[10]利用同質點對干涉圖進行自適應多視降噪，在非城市區(qū)域取得不錯的結果。文獻[11]首先采用像元分類技術排除水體等非分布式目標，然后采用AD(Anderson-Darling)檢驗的方法精化同質點，最后采用同質濾波(自適應多視)的方法提高信噪比獲取了更多的監(jiān)測點目標，在煤礦區(qū)域獲取了更多的形變信息。文獻[12]提出了改進的非局部濾波方法，即時空同質濾波。首先采用KS(Kolmogorov-Smirno)檢驗的方法確定同質點，考慮相位在同質地區(qū)具有相關性，然后采用了同質點進行空間的非局部濾波，加權值采用了常規(guī)的基于歐氏距離的指數(shù)衰減模型來確定。文獻[13]提出了CAESAR方法，采用協(xié)方差矩陣特征分解的方法對干涉圖濾波，其計算效率優(yōu)于SqueeSAR。但是該方法采用的協(xié)方差矩陣分解方法容易受到協(xié)方差矩陣的不穩(wěn)定性影響，從而影響后續(xù)的相位解算精度。

本文對干涉圖濾波時，首先參考PS與DS方法解算中對每一像元進行時間域協(xié)方差矩陣進行穩(wěn)健估計，然后通過協(xié)方差矩陣特征值分解，選取最大特征值對應的特征向量(相位)來代表該像元的相位以進行SAR干涉圖的降噪，最后采用真實SAR數(shù)據(jù)驗證了該方法的優(yōu)越性。

1 協(xié)方差矩陣的穩(wěn)健估計

一般認為SAR信號的實部和虛部服從均值為零的復高斯分布，地物目標往往具有空間相關性，服從復高斯分布[9,14]，則任一像元p在N景SAR影像中的觀測向量為d(p)=[d1(p)d2(p)…dN(p)]T，其概率密度函數(shù)可表示為

C(p)-1d(p))

(1)

(2)

式中，M是與任一點p具有相同分布的同質像元個數(shù)；m表示p點的同質像元。因此，對于任一像元，識別其同質點是進行協(xié)方差矩陣估計的前提。經典的識別同質點方法有基于復觀測值本身的散射數(shù)據(jù)(如強度或者幅度值)進行非參數(shù)檢驗法，如:KS(Kolmogorov-Smirno)檢驗[9]、AD(Anderson-Darling)檢驗[15]和BWS(Baumgartner-Weiβ-Schindler)檢驗[15]等。但由于這類方法選擇同質點時計算速度慢，尤其對于大數(shù)據(jù)集時其效率較差。本文采用基于統(tǒng)計區(qū)間估計的方法[16]，即

(3)

假設待估參數(shù)向量X，進行了n次觀測，即觀測值L=[l1l2…lN]T，f是L的概率密度函數(shù)，由極大似然估計有

(4)

胡貝爾將其廣義化為

(5)

同樣，聯(lián)合式(1)和式(2)，式(5)可具體化為

(6)

(7)

(8)

式中，w(·)表示等價權函數(shù)。一般地，等價權函數(shù)w(·)可由式(4)得

(9)

在實用中，由于SAR影像中與p點異質的點服從長尾分布，可用復多維t分布表示[17]

(10)

式中，v是t分布的自由度；Γ(·)是gamma函數(shù)。將式(10)代入式(9)，并考慮到復數(shù)t分布向實數(shù)t分布的變換，可得[18-19]

(11)

(12)

(13)

最終可得任一點的穩(wěn)健的協(xié)方差矩陣由加權的同質點強度來估計，即

(14)

該方法又稱為符號協(xié)方差矩陣，不需迭代也能達到抗差的效果[20]。

2 協(xié)方差矩陣的特征值分解

干涉圖每一像元的信號矢量中包含多類型的散射信號，圖1所示為SAR像元的不同散射機制模型。圖1(a)表示任意散射機制模型，是SAR影像中常見的散射模型；圖1(b)為永久散射體，即存在單一主導散射機制模型，在強相干區(qū)域很常見；圖1(c)為分布式目標的單一主導散射機制模型，主要存在于高分辨率像元內；圖1(d)為分布式目標的多主導散射機制模型，易出現(xiàn)在低分辨率像元內。

由于SAR時間序列數(shù)據(jù)的像元協(xié)方差矩陣的特征值對應不同強度的散射信號，可以通過對每一像元的協(xié)方差矩陣進行特征值分解來分離不同散射特征的信號，即

(15)

(16)

圖1 SAR像元的散射機制模型[21]Fig.1 Phase scattering mechanism model of one SAR pixel[21]

因此，對穩(wěn)健估計的協(xié)方差矩陣進行特征值分解后，找出最大的特征值和特征向量，即主導散射體，從而達到對干涉圖像元相位降噪的目的。由于本文采用的是高分辨率的TerraSAR-X的Strip模式數(shù)據(jù)，距離向和方位向的分辨率分別為1.7和0.9 m。分辨率單元內包含了較少的地物目標，所以假設分辨率單元地物類型單一，符合圖1(b)、(c)散射機制模型。因此本文基于穩(wěn)健估計的協(xié)方差矩陣分解的SAR干涉圖降噪只考慮單一主導的散射相位。

圖2為穩(wěn)健協(xié)方差矩陣分解SAR干涉圖降噪流程。該方法首先對多期SAR影像進行配準，并對強度圖進行輻射定標；然后基于定標后的SAR強度圖根據(jù)式(3)逐項元選取同質點；之后采用同質點構建待估像元的穩(wěn)健協(xié)方差矩陣；最后對協(xié)方差矩陣進行特征值分解，選取最大特征值對應的特征向量(相位)作為該像元的最終相位值，實現(xiàn)干涉圖像元降噪的目的。

圖2 穩(wěn)健協(xié)方差矩陣分解SAR干涉圖降噪流程Fig.2 Flowchart of SAR interferogram denoising based on robust covariance matrix decomposition

3 試驗結果與分析

3.1 試驗數(shù)據(jù)

本文采用8景覆蓋山西省清徐縣地面沉降區(qū)域的TerraSAR-X數(shù)據(jù)進行干涉圖去噪試驗，其參數(shù)列表見表1，影像中心入射角為23.87°。待去噪干涉圖由主影像20151013與從影像20150717組成，時間間隔為88 d，垂直基線為-45.63 m，干涉圖大小在SAR坐標系下為4000×4000像素。該區(qū)域由于受到地下水開采的影響產生嚴重的地面沉降和地裂縫等地質災害等現(xiàn)象[22-24]，而地面沉降主要發(fā)生在農田等低相干區(qū)域。采用常規(guī)InSAR濾波方法難以得到有效的形變信息。

表1 SAR干涉圖參數(shù)

注：*表示主影像。

3.2 同質點的選取試驗

利用式(3)采用8景SAR影像的強度圖進行同質點的選取，選取同質點的試驗區(qū)域如圖3(a)中P點箭頭所指示為某一池塘的邊緣。圖3(b)中中心白點代表待估計的參考點，本文同質點識別窗口選取15×15個像元，如圖3(b)灰色點陣所示。該窗口對應的實際距離約為30×30 m，充分考慮到形變區(qū)域的形變特征(特別是形變梯度)與協(xié)方差矩陣的穩(wěn)定性。圖3(c)為最終確定與該點同質的點。

圖3 同質點識別示意圖Fig.3 Sketch map of homogeneous point identification

圖4(a)為主影像采用同質點進行強度圖的降噪圖，該圖將用于計算協(xié)方差矩陣中的強度值，圖4(b)為8景SAR影像的平均強度圖，圖4(c)為主影像降噪前的強度圖，可見同質點濾波大大抑制了SAR強度圖的噪聲。

圖4 強度圖噪聲濾波前后對比圖Fig.4 The intensity maps comparison before and after noise filtering

為比較本文的濾波方法，采用文獻[8]提出的非線性濾波參數(shù)法(本文稱為改進的Goldstein濾波)來進行，即

(17)

式中，γ為無偏估計相干值，位于0～1；L為干涉圖視數(shù)。

3.3 干涉圖降噪對比分析

圖5—圖7為原始干涉圖、改進的Goldstein濾波和本文方法濾波后的干涉圖、相干圖及相干直方圖。圖5(a)為采用30 m分辨率的外部DEM差分后的原始差分干涉圖，圖中左上角的條紋為地面沉降信息，但由于失相干的影響，干涉圖噪聲很嚴重。由圖5(b)、(c)可以看出干涉圖相干性很低，特別在具有形變的農田區(qū)域，因此需要進行濾波才能獲取有用的地表形變信息。

對比圖5，圖6中改進的Goldstein濾波后的干涉圖干涉條紋明顯變清晰了，特別在城市等高相干地區(qū)。但由于低相干地區(qū)相干性改善不明顯，干涉條紋仍出現(xiàn)不連續(xù)現(xiàn)象。圖6(c)比圖5(c)的整體相干性大大提高，但兩個直方圖的分布特征類似。如圖7所示，本文方法降噪后的干涉圖的條紋變更加清晰，條紋的連續(xù)性得到增強，特別是在低相干區(qū)域，噪聲得到有效的抑制。圖7(c)的相干統(tǒng)計直方圖的分布有明顯改善，整體相干值顯著提升。

圖5 原始干涉圖及其質量圖Fig.5 Original interferogram and its quality map

圖6 改進的Goldstein濾波后的干涉圖及其質量圖Fig.6 Filtered interferogram with the improved Goldstein filer and its quality map

圖7 本文方法濾波后的干涉圖及其質量圖Fig.7 Denoised interferogram with the proposed method and its quality map

為了對比協(xié)方差矩陣穩(wěn)健估計的效果，對傳統(tǒng)協(xié)方差矩陣(式(2))分解后的相干值和經過穩(wěn)健估計的協(xié)方差矩陣(式(14))分解后的相干值進行統(tǒng)計，如圖8所示(橫軸代表相干性，縱軸代表像元的個數(shù))。由圖8可以看出，經過穩(wěn)健估計得到的協(xié)方差矩陣分解后的干涉圖相干性整體有所提升，從8.7e+06提升到9.0e+06。

以下選取3個典型的散射點區(qū)域，如圖3(a)p1、p2和p3分別對應低相干點(農田)、中等程度的相干點(池塘邊)和高相干點(房頂)，通過采用本文濾波方法前后相干性矩陣圖的變化來展示本文方法降噪的效果，如圖9—圖11所示。圖中(a)對應原始干涉圖的相干性矩陣，圖中(b)對應經過本文降噪后的相干性矩陣，由于本文總共采用8景SAR數(shù)據(jù)，所以每個點的相干矩陣為8×8。由3個相干性矩陣在濾波前后的對比分析可見，在3類散射區(qū)域，本文的降噪方法均能很好地抑制噪聲水平，特別對于中低相干區(qū)域(見圖9與圖10)。本文算法濾波后，相干性均得到顯著提升，這對于增加有效的地表監(jiān)測點信息至關重要。

為定量分析本文方法的降噪能力，采用相位梯度[25]和相位殘差點[9]作為定量評定的指標。相位梯度主要判斷相位的平滑程度，利用鄰域像元之間的相位差值，計算每個像元八鄰域的梯度APD，如式(18)所示，最后計算整個干涉圖的梯度和(SPD)，如式(19)

(18)

(19)

圖8 穩(wěn)健估計(C-M)與傳統(tǒng)協(xié)方差矩陣(C-MLE)分解后相干性統(tǒng)計圖Fig.8 Coherent statistical analysis afterdecomposition of robust covariance matrix estimation and traditional covariance matrix estimation

圖9 低相干點的相干矩陣對比圖(農田，如圖4(a)中p1所示)Fig.9 Coherence matrix of low coherence points (in the farmland, p1 in Fig.4(a))

對本文提出的干涉圖降噪方法和改進的Goldstein濾波兩種方法分別進行了干涉圖濾波對比試驗，其濾波相位質量統(tǒng)計結果見表2。首先，相位梯度和(SPD)指標表明，傳統(tǒng)的協(xié)方差矩陣分解法和改進的Goldstein濾波后的相位梯度和均值分別為0.40 rad和0.53 rad，減少比例比改進的Goldstein濾波方法從15.47%提高到35.27%，而本文具有穩(wěn)健估計的方法比傳統(tǒng)的協(xié)方差矩陣分解法從35.27%提高到43.92%，而相位梯度和均值降低至0.35 rad，相比改進的Goldstein濾波，降噪能力提升了2.8倍。其次從相位殘差點的統(tǒng)計結果可以看出，本文方法使得干涉圖相位殘差點由2 635 542個減少到667 200個，減少比例為74.68%，約為改進的Goldstein濾波干涉圖的殘差點的1/3，可見在低相干區(qū)域相位不連續(xù)性得到很大的改善。進一步發(fā)現(xiàn)穩(wěn)健估計方法比傳統(tǒng)的協(xié)方差矩陣分解法的殘差點少近1/3，說明穩(wěn)健的協(xié)方差矩陣估計對于干涉圖噪聲的減弱效果明顯。總之，基于協(xié)方差矩陣分解的干涉圖降噪法比改進的Goldstein濾波在相干性提高與有效目標點的增加方面均有顯著效果，特別是在低相干區(qū)域的農田等區(qū)域。進一步試驗發(fā)現(xiàn)穩(wěn)健的協(xié)方差矩陣估計法比傳統(tǒng)的協(xié)方差矩陣分解法的整體相干性又有一定的提高。

圖10 中等相干點的相干矩陣對比圖(池塘邊,如圖4(a)中p2所示)Fig.10 Coherence matrix of moderate coherence points(near the pool, p2 in Fig.4(a))

圖11 高相干點的相干矩陣對比圖(房頂，如圖4(a)中p3所示)Fig.11 Coherence matrix of high coherence points(on the roof, p3 in Fig.4(a))

表2 干涉圖去噪質量統(tǒng)計

同樣本文采用相位剖線對比原始相位、改進的Goldstein濾波和本文方法，即3種方法的降噪能力及相位精度，如圖7(a)中的P剖線。該條剖線的相位值基本穩(wěn)定在-1 radian左右，如圖12所示，其中淺色圓點代表原始相位，黑色圓點代表改進的Goldstein濾波，星號點代表本文方法。圖12表明本文方法使得纏繞相位更加平滑，相位噪聲得到有效控制，而原始相位基本處于一種隨機飄動的狀態(tài)，改進的Goldstein濾波相位同樣噪聲表現(xiàn)嚴重。

圖12 3個干涉圖的相位剖線比較圖(剖線P位置如圖7(a)所示)Fig.12 The comparison of phases of three interferograms(the profile P shown in Fig.7(a))

4 結 論

本文運用基于穩(wěn)健估計的協(xié)方差矩陣分解方法對干涉圖進行降噪。該方法基于同質點估計協(xié)方差矩陣，對少量引入的異質點采用穩(wěn)健估計的方法進行處理，理論上保證了協(xié)方差矩陣的無偏性；在此基礎上對無偏的協(xié)方差矩陣進行特征值分解，選取最大特征值對應的特征向量作為單一主導散射體的像元或只含一種散射機制的分布式散射體像元的最終相位，從而達到對干涉圖降噪的目的。通過覆蓋山西清徐地表形變區(qū)域的8景Strip模式的TerraSAR-X數(shù)據(jù)試驗，表明本文降噪方法比改進的Goldstein濾波得到的干涉圖條紋更加清晰；定量的相位梯度和與相位殘差點對比發(fā)現(xiàn)本文提出的降噪方法濾波效果更優(yōu)，特別在低相干區(qū)域的相位相干性得到了明顯提高，增加了低相干區(qū)域的監(jiān)測點密度，這對于低相干區(qū)域InSAR技術的應用具有重要作用。