SAR干涉圖降噪的穩(wěn)健協(xié)方差矩陣分解法

趙超英，王寶行

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

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

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

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

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

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

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

C(p)-1d(p))

(1)

(2)

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

(3)

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

(4)

胡貝爾將其廣義化為

(5)

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

(6)

(7)

(8)

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

(9)

在實(shí)用中，由于SAR影像中與p點(diǎn)異質(zhì)的點(diǎn)服從長尾分布，可用復(fù)多維t分布表示[17]

(10)

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

(11)

(12)

(13)

最終可得任一點(diǎn)的穩(wěn)健的協(xié)方差矩陣由加權(quán)的同質(zhì)點(diǎn)強(qiáng)度來估計(jì)，即

(14)

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

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

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

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

(15)

(16)

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

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

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

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

3 試驗(yàn)結(jié)果與分析

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

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

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

注：*表示主影像。

3.2 同質(zhì)點(diǎn)的選取試驗(yàn)

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

圖3 同質(zhì)點(diǎn)識別示意圖Fig.3 Sketch map of homogeneous point identification

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

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

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

(17)

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

3.3 干涉圖降噪對比分析

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

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

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

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

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

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

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

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

(18)

(19)

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

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

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

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

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

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

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

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

4 結(jié) 論

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