融合CZT頻譜細(xì)化算法的InSAR相位濾波方法

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

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

0 引言

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

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

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

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

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

如果說N點(diǎn)的快速傅里葉變換是在Z平面的單位圓上的等間距取N個(gè)點(diǎn)，那么CZT則是Z平面螺旋線的周線上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á)到對指定頻譜細(xì)化的目的[9-10]。

假設(shè)采樣序列x(n)，n=0，1，…，N-1中，采樣頻率為fs，其中f1～f2為需要細(xì)化的頻譜條帶，對信號進(jìn)行N點(diǎn)快速傅里葉變換，故有原有信號的頻率分辨率為Δf=fs/N，所選需要細(xì)化條帶中心頻率為f0=(f1+f2)/2，帶寬為fw=f2-f1。用X(z)來表示x(n)序列的Z變換，通過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平面螺旋線的伸展率；φ0表示相鄰采樣點(diǎn)間的角度差。根據(jù)式(1)所給的Zk值，推導(dǎo)可得Z變換X(Zk)，如式(2)所示。

(2)

令

(3)

(4)

可得

(5)

令

(6)

可得

(7)

在實(shí)際過程中，頻帶點(diǎn)的采樣點(diǎn)數(shù)M需要根據(jù)整個(gè)頻譜的頻譜分辨率Δf′來確定，它們之間的關(guā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)行處理濾波。其濾波算法過程[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)

選定合適濾波窗口大小，對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)為過濾器參數(shù)。當(dāng)過濾器參數(shù)α=0時(shí)，不對其進(jìn)行平滑濾波，當(dāng)α的值過大時(shí)，濾波結(jié)果將會在相位跳變密集區(qū)過度濾波，使得精度下降。為了提高濾波效果，α的取值將根據(jù)兩張SAR影像的干涉相關(guān)系數(shù)，即α=1-r(0≤r≤1)。最后將經(jīng)過濾波后的頻率域H(u，v)做二維快速傅里葉逆變換得到濾波后的結(jié)果。

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

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

由于傳統(tǒng)的頻域?yàn)V波方法在相對較低的頻域分辨率中分辨有用信息與噪聲的能力有限，為了能夠?qū)⑦@些噪聲復(fù)雜區(qū)域進(jìn)行良好的濾波，同時(shí)保持干涉條紋緊湊區(qū)域的信息不被破壞，本節(jié)結(jié)合了上述CZT頻譜細(xì)化和自適應(yīng)Goldstein濾波方法的特點(diǎn)，提出了CZT頻譜細(xì)化與自適應(yīng)濾波相融合的方法。通過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í)盡量減小對條紋密集區(qū)細(xì)節(jié)的丟失。其具體算法流程如下。

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

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

(10)

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

(11)

(12)

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

ih=0，1，…，Nr-1

(13)

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

(14)

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

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

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

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

通過仿真實(shí)驗(yàn)與實(shí)測數(shù)據(jù)來驗(yàn)證本文提出的方法。其中，仿真數(shù)據(jù)來源于MATLAB軟件中peak函數(shù)生成的理想干涉圖(圖4(a))與添加了真實(shí)SAR影像噪聲的仿真干涉圖(圖4(b))，用來模擬真實(shí)噪聲的影響。實(shí)測數(shù)據(jù)來源于將Sentinel-1A的單視復(fù)數(shù)影像通過ENVI軟件中的SARscape模塊通過裁剪、配準(zhǔn)等預(yù)處理所產(chǎn)生的干涉圖(圖7(a))。為了驗(yàn)證本文濾波方法的有效性，將本文方法與Boxcar、NL-mean和自適應(yīng)Goldstein濾波方法作對比，從目視結(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 評價(jià)指標(biāo)

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

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

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

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

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

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

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

此外，通過對比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ì)圖可看出本文研究方法相對于其他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ì)對比

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

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

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

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

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

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

4 結(jié)束語

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