張繼超，李繼虎，趙鵬飛，林泓全

(1.遼寧工程技術(shù)大學(xué) 測(cè)繪與地理科學(xué)學(xué)院，遼寧 阜新 123000；2.遼寧工程技術(shù)大學(xué) 軟件學(xué)院，遼寧 葫蘆島 125105)

0 引言

合成孔徑雷達(dá)干涉測(cè)量技術(shù)廣泛應(yīng)用于大區(qū)域的地表高程反演和長(zhǎng)期形變監(jiān)測(cè)，其原理是通過(guò)對(duì)同一地表區(qū)域多次觀(guān)測(cè)的SAR影像進(jìn)行干涉，從而獲取較高精度的干涉相位信息，干涉相位圖的質(zhì)量直接決定干涉相位測(cè)量的精度。在InSAR干涉測(cè)量中，由于地表環(huán)境復(fù)雜、大氣環(huán)境因素導(dǎo)致的傳感器接收地表信號(hào)不佳和受到熱噪聲去相干、時(shí)間去相干、基線(xiàn)去相干以及信號(hào)去相干等多種因素影響而產(chǎn)生噪聲復(fù)雜區(qū)域，在濾波效果較差的情況下會(huì)導(dǎo)致解纏過(guò)程中由于產(chǎn)生較多奇異點(diǎn)使得相位解纏精度下降，甚至無(wú)法解纏，而且還會(huì)導(dǎo)致后續(xù)衍生產(chǎn)品質(zhì)量不佳[1-2]。因此，為了減小噪聲對(duì)干涉測(cè)量結(jié)果產(chǎn)生的負(fù)面影響，并保證干涉測(cè)量結(jié)果的準(zhǔn)確，在相位解纏之前對(duì)干涉相位進(jìn)行適合的濾波是必不可少的。

現(xiàn)階段，用于干涉相位濾波的方法大致分為頻率域?yàn)V波和空間域?yàn)V波。空間域?yàn)V波主要有均值濾波、中值濾波和圓周期均值等方法，其原理簡(jiǎn)單、運(yùn)算速度快，但沒(méi)有考慮干涉相位圖的實(shí)際含義，濾波效果不盡人意，而且在濾波窗口的選擇上有一定難度，在干涉條紋密集區(qū)域容易破壞細(xì)節(jié)，出現(xiàn)失真的情況。隨著空間域?yàn)V波的發(fā)展，又有學(xué)者將自適應(yīng)濾波和自適應(yīng)窗口大小的濾波方法等相關(guān)技術(shù)應(yīng)用于干涉相位濾波中，例如精細(xì)Lee濾波，通過(guò)自適應(yīng)的方法可以根據(jù)圖像的紋理特征相似性來(lái)更好地消除噪聲并保留圖像中的細(xì)節(jié)，但其計(jì)算復(fù)雜，而且仍舊無(wú)法很好保留干涉跳變信息[3-4]。頻率域?yàn)V波是將圖像變換到頻率域，通過(guò)去除高頻噪聲部分從而實(shí)現(xiàn)濾波[5]。1998年，Goldstein和Werner提出了一種在頻率域中的自適應(yīng)干涉相位濾波方法，這一具有里程碑意義的濾波算法是采用傅里葉變換將干涉相位變換到頻率域平滑處理，再將平滑后的頻率域干涉相位做傅里葉逆變換回到空間域。經(jīng)典Goldstein濾波是典型的干涉相位濾波方法，它在去噪與保持影像分辨率的效果較好，但在選取參數(shù)主觀(guān)因素較強(qiáng)。在2003年，Baran等[6]將Goldstein濾波的濾波參數(shù)根據(jù)兩張SAR影像的相干性系數(shù)相聯(lián)系，使之可以根據(jù)相關(guān)性的好壞自動(dòng)地改變?yōu)V波參數(shù)，從而改變?yōu)V波強(qiáng)度，使得傳統(tǒng)濾波效果進(jìn)一步提升。2019年，Hensley[7]用Goldstein與Werner濾波結(jié)合的方法來(lái)解決因長(zhǎng)時(shí)間、長(zhǎng)基線(xiàn)引起的對(duì)噪聲的影響，但該方法在噪聲復(fù)雜區(qū)域去噪效果不佳，使得在高噪聲區(qū)域?yàn)V波質(zhì)量仍然不高。

為了更有效抑制噪聲復(fù)雜區(qū)域，并保持條紋密集區(qū)域的信息不被破壞，以達(dá)到濾波結(jié)果有利于解纏和提升干涉測(cè)量精度的目的，本文提出了一種融合頻譜細(xì)化和自適應(yīng)Goldstein相位濾波的方法，將圖像轉(zhuǎn)換至頻率域后，通過(guò)線(xiàn)性調(diào)頻Z變換(CZT)的方法進(jìn)行局部的頻率估計(jì)，提高計(jì)算的頻率分辨率，使得在高噪聲區(qū)域的頻譜得以細(xì)化以達(dá)到更好區(qū)分噪聲與實(shí)際有用相位信息的目的；通過(guò)自適應(yīng)的方式，可根據(jù)相干影像的相關(guān)性改變?yōu)V波強(qiáng)度，以達(dá)到去噪與保留條紋密集區(qū)信息的兼顧。將本文方法與Boxcar濾波、非局部均值濾波(non-local means，NL-Means)和自適應(yīng)Goldstein濾波進(jìn)行對(duì)比，通過(guò)仿真實(shí)驗(yàn)與實(shí)際數(shù)據(jù)驗(yàn)證本文方法的有效性。

1 融合頻譜細(xì)化的自適應(yīng)濾波原理

1.1 線(xiàn)性調(diào)頻Z變換頻譜細(xì)化算法原理

如果說(shuō)N點(diǎn)的快速傅里葉變換是在Z平面的單位圓上的等間距取N個(gè)點(diǎn)，那么CZT則是Z平面螺旋線(xiàn)的周線(xiàn)上Z變換取等間距N個(gè)點(diǎn)[8]。CZT頻譜細(xì)化可以計(jì)算在N點(diǎn)中輸入數(shù)列為x(n)的給定一條路徑的Z變換，采用增加輸出點(diǎn)數(shù)，使得輸出點(diǎn)多余輸入點(diǎn)數(shù)就可以達(dá)到對(duì)指定頻譜細(xì)化的目的[9-10]。

假設(shè)采樣序列x(n)，n=0，1，…，N-1中，采樣頻率為fs，其中f1～f2為需要細(xì)化的頻譜條帶，對(duì)信號(hào)進(jìn)行N點(diǎn)快速傅里葉變換，故有原有信號(hào)的頻率分辨率為Δf=fs/N，所選需要細(xì)化條帶中心頻率為f0=(f1+f2)/2，帶寬為fw=f2-f1。用X(z)來(lái)表示x(n)序列的Z變換，通過(guò)CZT算法就可以計(jì)算給定點(diǎn)Zk上的Z變換X(Zk)。

Zk=AW-k，k=0，1，…，M-1

(1)

式中：A=A0e-jθ0；W=W0e-jφ0；M為所要分析復(fù)數(shù)頻譜的采樣點(diǎn)數(shù);j為虛數(shù)。其中，A0和θ0分別表示第一個(gè)采樣點(diǎn)Z0的矢量半徑的大小與相角；W0為Z平面螺旋線(xiàn)的伸展率；φ0表示相鄰采樣點(diǎn)間的角度差。根據(jù)式(1)所給的Zk值，推導(dǎo)可得Z變換X(Zk)，如式(2)所示。

(2)

令

(3)

(4)

可得

(5)

令

(6)

可得

(7)

在實(shí)際過(guò)程中，頻帶點(diǎn)的采樣點(diǎn)數(shù)M需要根據(jù)整個(gè)頻譜的頻譜分辨率Δf′來(lái)確定，它們之間的關(guān)系滿(mǎn)足公式：Δf′=fw/(M-1)。由于此次頻譜分析在單位圓上進(jìn)行CZT算法的實(shí)現(xiàn)，因此A0與W0取1，θ0=2πf1/fs，φ0=2πΔf′/fs。

1.2 自適應(yīng)Goldstein算法模型

Goldstein濾波算法是基于一種傅里葉變換的自適應(yīng)濾波方法，將干涉相位轉(zhuǎn)換至復(fù)數(shù)形式，有效避免干涉相位圖在相位跳變處的影像，并在空間域進(jìn)行處理濾波。其濾波算法過(guò)程[11-13]如下所示。

將干涉相位圖(x，y)轉(zhuǎn)換為復(fù)數(shù)形式的矢量空間E(x，y)，如式(8)所示。

E(x，y)=exp(φ(x，y))=cos(φ(x，y))+j·sin(φ(x，y))

(8)

選定合適濾波窗口大小，對(duì)E(x，y)進(jìn)行二維快速傅里葉變換到頻率域Z(u，v)(其中u，v表示空間頻率)，并進(jìn)行Goldstein濾波處理得到H(u，v)，如式(9)所示。

H(u，v)=S{|Z(u，v)|}α·Z(u，v)

(9)

式中：S{}為濾波平滑算子；α(0≤α≤1)為過(guò)濾器參數(shù)。當(dāng)過(guò)濾器參數(shù)α=0時(shí)，不對(duì)其進(jìn)行平滑濾波，當(dāng)α的值過(guò)大時(shí)，濾波結(jié)果將會(huì)在相位跳變密集區(qū)過(guò)度濾波，使得精度下降。為了提高濾波效果，α的取值將根據(jù)兩張SAR影像的干涉相關(guān)系數(shù)，即α=1-r(0≤r≤1)。最后將經(jīng)過(guò)濾波后的頻率域H(u，v)做二維快速傅里葉逆變換得到濾波后的結(jié)果。

這種使用相干性系數(shù)作為過(guò)濾器的平滑參數(shù)的特點(diǎn)是在相干性好的區(qū)域可以防止過(guò)度濾波，從而使干涉圖的分辨率降低，而在相干性差、噪聲多的區(qū)域增強(qiáng)濾波效果，達(dá)到對(duì)噪聲更好的過(guò)濾。

1.3 基于融合CZT頻譜細(xì)化的自適應(yīng)濾波方法

由于傳統(tǒng)的頻域?yàn)V波方法在相對(duì)較低的頻域分辨率中分辨有用信息與噪聲的能力有限，為了能夠?qū)⑦@些噪聲復(fù)雜區(qū)域進(jìn)行良好的濾波，同時(shí)保持干涉條紋緊湊區(qū)域的信息不被破壞，本節(jié)結(jié)合了上述CZT頻譜細(xì)化和自適應(yīng)Goldstein濾波方法的特點(diǎn)，提出了CZT頻譜細(xì)化與自適應(yīng)濾波相融合的方法。通過(guò)CZT頻譜細(xì)化的方法細(xì)化相位干涉頻率域的頻率分辨率，區(qū)別噪聲頻率與有用信息相位的頻率界限，并根據(jù)兩幅干涉的主副影像在噪聲復(fù)雜區(qū)域相關(guān)性較差，條紋密集區(qū)有著較高相關(guān)性的特點(diǎn)改變?yōu)V波參數(shù)的大小實(shí)現(xiàn)在噪聲復(fù)雜區(qū)域?yàn)V波良好的同時(shí)盡量減小對(duì)條紋密集區(qū)細(xì)節(jié)的丟失。其具體算法流程如下。

步驟1：首先將干涉相位轉(zhuǎn)換至矢量空間，使得跳變的不連續(xù)的相位變成連續(xù)的復(fù)相位數(shù)據(jù)，通過(guò)式(8)得到復(fù)干涉相位E(x，y)。在一個(gè)復(fù)干涉相位的小窗口內(nèi)(窗口大小為Na×Nr)，相位變化的坡度可以認(rèn)為是近似相等的常數(shù)，即通過(guò)快速傅里葉變換將復(fù)干涉相位轉(zhuǎn)換到頻率域中的僅有的一個(gè)主頻，這個(gè)主頻對(duì)應(yīng)的相位如式(10)所示。

φ(m+p,n+q)=2πp·fa+2πq·fr+φ(m,n)

(10)

步驟2：針對(duì)復(fù)干涉相位轉(zhuǎn)換到頻率域中的頻譜峰值作為主瓣始點(diǎn)進(jìn)行頻譜細(xì)化，通過(guò)式(1)至式(7)，對(duì)計(jì)算窗口的方位向與距離向取3×32個(gè)采樣點(diǎn)做CZT變換，對(duì)主頻譜瓣進(jìn)行32倍采樣估計(jì)，可得式(11)。

(11)

(12)

則在這個(gè)小窗口內(nèi)的線(xiàn)性相位細(xì)化頻譜如式(13)所示。

ih=0，1，…，Nr-1

(13)

步驟4：根據(jù)細(xì)化后的復(fù)干涉相位頻帶做自適應(yīng)濾波(式(14))。

(14)

根據(jù)上述步驟原理，其流程圖如圖1所示。

圖1 本文干涉相位濾波方法流程

2 實(shí)驗(yàn)方案與評(píng)價(jià)指標(biāo)

2.1 實(shí)驗(yàn)方案

通過(guò)仿真實(shí)驗(yàn)與實(shí)測(cè)數(shù)據(jù)來(lái)驗(yàn)證本文提出的方法。其中，仿真數(shù)據(jù)來(lái)源于MATLAB軟件中peak函數(shù)生成的理想干涉圖(圖4(a))與添加了真實(shí)SAR影像噪聲的仿真干涉圖(圖4(b))，用來(lái)模擬真實(shí)噪聲的影響。實(shí)測(cè)數(shù)據(jù)來(lái)源于將Sentinel-1A的單視復(fù)數(shù)影像通過(guò)ENVI軟件中的SARscape模塊通過(guò)裁剪、配準(zhǔn)等預(yù)處理所產(chǎn)生的干涉圖(圖7(a))。為了驗(yàn)證本文濾波方法的有效性，將本文方法與Boxcar、NL-mean和自適應(yīng)Goldstein濾波方法作對(duì)比，從目視結(jié)果、相位殘差點(diǎn)數(shù)量、相位導(dǎo)數(shù)標(biāo)準(zhǔn)差以及仿真實(shí)驗(yàn)誤差的角度進(jìn)行分析。實(shí)驗(yàn)方法流程如圖2所示。

圖2 實(shí)驗(yàn)方案流程

2.2 評(píng)價(jià)指標(biāo)

為了更好地對(duì)比相位濾波效果，本文除了在目視效果上作對(duì)比外，還加入了相位殘差點(diǎn)與相位導(dǎo)數(shù)標(biāo)準(zhǔn)差作為對(duì)比依據(jù)，并且在后續(xù)的仿真實(shí)驗(yàn)中對(duì)比了這幾種濾波方法的誤差效果。

相位殘差點(diǎn)的多少將直接影響相位解纏的質(zhì)量，濾波后相位殘差點(diǎn)的數(shù)量越少，對(duì)相位解纏的效果越有利，其解纏的精度越高[14-15]。

相位導(dǎo)數(shù)標(biāo)準(zhǔn)差是評(píng)價(jià)干涉圖質(zhì)量的有效方式，其值越大，說(shuō)明干涉效果越差。

3 實(shí)驗(yàn)結(jié)果分析

3.1 仿真實(shí)驗(yàn)

為了更準(zhǔn)確地對(duì)比實(shí)驗(yàn)結(jié)果，本節(jié)使用仿真模擬干涉相位來(lái)對(duì)比本文濾波方法和其他3種相位濾波的提升效果。本小節(jié)仿真模擬相位圖由MATLAB軟件的peaks函數(shù)生成，大小為512×512，如圖3(a)所示。模擬真實(shí)SAR影像噪聲的仿真干涉相位圖如圖3(b)所示。圖3(c)至圖3(f)是4種濾波結(jié)果圖。通過(guò)這4幅濾波結(jié)果圖，可以較為明顯地判斷出Boxcar濾波與NL-mean濾波去噪效果較差，自適應(yīng)Goldstein濾波與本文方法濾波效果較好。為了更好地顯示4種濾波方法在條紋密集區(qū)的濾波效果，圖3(g) 至圖3(j)給出了4種方法在圖上標(biāo)出的條紋密集區(qū)濾波結(jié)果放大圖?？梢钥闯?，本文方法在去噪效果與干涉密集條紋保持情況相較于其他3種濾波方法更好。

圖3 仿真數(shù)據(jù)及模擬干涉相位圖濾波結(jié)果與條紋密集區(qū)濾波放大圖

此外，通過(guò)對(duì)比4種方法濾波結(jié)果的相位殘差點(diǎn)數(shù)(表1)，不論在全局區(qū)域還是條紋密集區(qū)域，本文研究的濾波方法明顯少于另外3種濾波方法，這也意味著該方法更有益于后續(xù)的相位解纏階段。圖4表示4種濾波方法的相位導(dǎo)數(shù)標(biāo)準(zhǔn)差統(tǒng)計(jì)情況。從相位導(dǎo)數(shù)標(biāo)準(zhǔn)差統(tǒng)計(jì)圖可看出本文研究方法相對(duì)于其他3種方法的相位導(dǎo)數(shù)標(biāo)準(zhǔn)差值的數(shù)量更靠近0，表明本文方法濾波結(jié)果的干涉圖質(zhì)量相較于另外3種方法有不同程度的提升。從仿真實(shí)驗(yàn)的誤差統(tǒng)計(jì)圖(圖5)可以看出，這4種濾波方法的誤差結(jié)果都呈現(xiàn)近似正態(tài)分布，而本文濾波方法在全局區(qū)域和條紋密集區(qū)的誤差所呈現(xiàn)的正態(tài)分布尖峰更高，證明本文方法在全局區(qū)域與條紋密集區(qū)域的濾波精度相較于其他3種方法有所提升，證明本文方法在去噪效果與保持條紋密集區(qū)域的信息上有更好效果。

表1 不同濾波方法結(jié)果的相位殘差點(diǎn)數(shù)量

圖4 4種濾波方法的相位導(dǎo)數(shù)標(biāo)準(zhǔn)差統(tǒng)計(jì)對(duì)比

圖5 4種濾波方法結(jié)果及條紋密集區(qū)濾波結(jié)果的誤差統(tǒng)計(jì)對(duì)比

3.2 實(shí)測(cè)數(shù)據(jù)實(shí)驗(yàn)

本節(jié)使用的實(shí)測(cè)數(shù)據(jù)為哨兵1號(hào)SAR影像數(shù)據(jù)，選用2018年12月12日與2018年12月23日福建省泉州市某區(qū)域裁剪后的影像進(jìn)行干涉，截取了600像素×600像素大小的干涉區(qū)域作為實(shí)驗(yàn)對(duì)象。所選區(qū)域?yàn)槎嗌角伊謪^(qū)密集，干涉形成的噪聲較為復(fù)雜且干涉條紋比較密集有利于進(jìn)一步驗(yàn)證本文算法良好的去噪能力與保持條紋細(xì)節(jié)信息的能力。圖6為實(shí)測(cè)數(shù)據(jù)和4種方法的濾波結(jié)果?？梢钥闯觯谙嗤笮?3×13的濾波窗口條件下，Boxcar與NL-mean濾波方法在去噪能力上較弱，自適應(yīng)Goldstein濾波無(wú)法較好地保持條紋細(xì)節(jié)，而本文濾波方法在良好去除噪聲的同時(shí)保證了緊湊條紋信息能夠較好保存。此外，在表1所示不同方法濾波結(jié)果所含殘差點(diǎn)數(shù)目中，本文濾波結(jié)果后的殘差點(diǎn)最少，更利于相位解纏。圖7和圖8分別對(duì)比了4種方法的相位導(dǎo)數(shù)標(biāo)準(zhǔn)差及其統(tǒng)計(jì)圖，本文濾波方法由于相位導(dǎo)數(shù)標(biāo)準(zhǔn)差值更多地趨向于0，因此，本文濾波方法優(yōu)于其他3種濾波方法。

圖6 實(shí)測(cè)數(shù)據(jù)和4種濾波方法結(jié)果

圖7 實(shí)測(cè)數(shù)據(jù)和4種濾波方法干涉結(jié)果的相位導(dǎo)數(shù)標(biāo)準(zhǔn)差圖對(duì)比

圖8 4種濾波方法干涉結(jié)果的相位導(dǎo)數(shù)標(biāo)準(zhǔn)差統(tǒng)計(jì)圖對(duì)比

4 結(jié)束語(yǔ)

本文將頻譜細(xì)化原理與自適應(yīng)濾波相結(jié)合，通過(guò)提升圖像變換為頻率域的頻率分辨率的方式，更好地區(qū)分干涉相位中的噪聲信息來(lái)提升濾波效果與濾波質(zhì)量，在保證良好的去除噪聲能力的同時(shí)又保證了條紋細(xì)節(jié)不被破壞。通過(guò)仿真實(shí)驗(yàn)與實(shí)測(cè)數(shù)據(jù)，用本文方法與Boxcar濾波、非局部均值濾波以及自適應(yīng)Goldstein濾波方法作實(shí)驗(yàn)對(duì)比，在相位殘差點(diǎn)數(shù)量和相位導(dǎo)數(shù)標(biāo)準(zhǔn)差中，本文的濾波方法相對(duì)于其他3種方法效果更好。在仿真實(shí)驗(yàn)誤差分析中，本文濾波方法對(duì)數(shù)據(jù)濾波后的準(zhǔn)確性較高，有利于提高干涉測(cè)量的精度。通過(guò)以上實(shí)驗(yàn)內(nèi)容證明了本文濾波方法的有效性。