張 斌, 胡慶榮, 韋立登, 李 爽

(1. 北京無線電測量研究所, 北京 100854; 2. 中國航天科工集團(tuán)第二研究院, 北京 100854)

0 引 言

干涉合成孔徑雷達(dá)(interferometric synthetic aperture radar, InSAR)通過利用兩幅或多幅SAR圖像的相位信息可以獲取大范圍、高精度的地表三維信息[1-2]。其利用相位信息反演高程信息,但其干涉相位圖仍然存在較為嚴(yán)重的噪聲(雷達(dá)系統(tǒng)的熱噪聲、SAR處理過程中的相干斑、干涉基線引起的相位噪聲、雷達(dá)陰影、配準(zhǔn)誤差等),而噪聲(尤其斑點噪聲)的存在使得干涉相位圖出現(xiàn)很多孤立的殘差點,同時相位的連續(xù)性和一致性受到影響,給后續(xù)相位解纏造成困難,也影響InSAR的測高精度[3-5]。因此在InSAR處理過程中,有必要對干涉相位圖進(jìn)行濾波操作,最大限度消除殘差點對相位的影響,提高信噪比和相位解纏的精度及效率。

目前,干涉圖濾波方法主要包括空間域濾波、頻率域濾波、時頻域濾波和幾何學(xué)濾波等4大類。圓周期均值或中值濾波[6]作為典型的空間域濾波方法,其算法較為簡單,運算速度快,但其在條紋密集區(qū)域和條紋稀疏區(qū)域時效果差別較大,條紋密集區(qū)域容易出現(xiàn)細(xì)節(jié)丟失現(xiàn)象。繼而出現(xiàn)的魯棒加權(quán)圓周濾波[7]結(jié)合干涉圖中條紋的特點,使濾波效果得到提高,但其實現(xiàn)較為復(fù)雜,運算量大;而梯度自適應(yīng)濾波[8]利用像素的噪聲統(tǒng)計特性,對于干涉條紋邊緣能起到較好的保護(hù)作用,但不可避免的是部分噪聲區(qū)域因為梯度值大于k值(保護(hù)邊緣的作用)而被保留,使得該處噪聲濾波效果不佳。在頻域方法中,典型的算法為Goldstein算法[9],該算法通過數(shù)據(jù)分塊和頻域加權(quán)進(jìn)行濾波,效果優(yōu)良,但其受分塊大小和頻域加權(quán)參數(shù)的影響較大,并且有過度濾波的傾向。時頻域方法中,文獻(xiàn)[10]提出的離散小波變換相位濾波方法能較好地保持圖像的干涉條紋信息,但去噪效果不好。

而數(shù)學(xué)形態(tài)學(xué)[11]作為幾何學(xué)方法具有描述圖像幾何形態(tài),保持圖像紋理信息的性能,故基于形態(tài)學(xué)的圖像濾波算法適合于干涉相位圖的濾波。形態(tài)學(xué)的濾波方法中,文獻(xiàn)[12]針對干涉相位圖的特點,提出一種修正的灰度形態(tài)學(xué)開閉操作濾波器,濾波性能優(yōu)于中值濾波、模態(tài)濾波器等,但其濾波效果仍有待提高;文獻(xiàn)[13]將干涉圖進(jìn)行量化處理,再運用二值形態(tài)學(xué)濾波處理,最后將量化后的濾波圖進(jìn)行融合得到最終的干涉濾波圖,雖然該算法的處理效果一定程度上優(yōu)于Candeias算法,但其量化融合的范圍受到限制;文獻(xiàn)[14]提出了基于形態(tài)學(xué)開閉重建運算的融合算法(modified version of alternate sequential filters with reconstruction, MASF),從理論和實測數(shù)據(jù)兩方面進(jìn)行了驗證,并與局部頻率估計濾波方法[15]進(jìn)行對比顯示出其方法的優(yōu)勢,但其在不同噪聲環(huán)境下的穩(wěn)定性和局部細(xì)節(jié)信息表現(xiàn)仍有待提高?；谏鲜龇治?本文提出了一種改進(jìn)的自適應(yīng)梯度形態(tài)學(xué)重建算法,對梯度估計進(jìn)行修正求取更優(yōu)的線結(jié)構(gòu)元素,并借鑒了Candeias中改進(jìn)的形態(tài)學(xué)開閉運算和MASF方法的基于膨脹和腐蝕重建操作的思想進(jìn)行濾波,能在去除噪聲的同時有效保護(hù)相位的連續(xù)性。文章通過相位圖像矢量化(正余弦化)和改進(jìn)自適應(yīng)結(jié)構(gòu)元素使得濾波更加合理,算法的性能與兩種形態(tài)學(xué)方法(Candeias方法、MASF方法)以及中值圓周濾波方法相比相位的連續(xù)性得到提升,殘差點數(shù)目也更少。

1 數(shù)學(xué)形態(tài)學(xué)的相位濾波

在干涉圖濾波方法中,數(shù)學(xué)形態(tài)學(xué)能夠保持相位的跳變特性和相關(guān)紋理,其利用結(jié)構(gòu)元素在相位圖像上不斷移動來確定相位之間的關(guān)系,通過其幾何特征濾除疊加在相位圖像上的噪聲。本文算法中在數(shù)學(xué)形態(tài)學(xué)中的應(yīng)用算法主要為改進(jìn)的形態(tài)學(xué)開閉操作和基于膨脹和腐蝕的重建操作。

1.1 改進(jìn)的形態(tài)學(xué)開閉操作

(1)

εB(φ)(x)=min{φ(y)|y∈Bx∩E}

(2)

在相位濾波中,由于相位值在[-π,π]的限制,如果直接使用膨脹和腐蝕操作進(jìn)行濾波,會造成局部濾波結(jié)果不合理。而Candeias[12]對膨脹和腐蝕操作進(jìn)行了修改,提出了擴(kuò)大和縮小操作,使得其在進(jìn)行運算時更符合真實值。擴(kuò)大和縮小操作定義[12]為

(3)

min{φ(y)+L(v-φ(y))|y∈Bx∩E}

(4)

式中,+L表示模L加,此處L=2π;v表示取值范圍的中心值,此處v=0。

圖1(a)表示了一維相位值中的相鄰三點,其相位模糊如圖1(b)所示,圖1(c)和圖1(d)分別顯示了腐蝕操作和縮小操作的結(jié)果。

圖1 腐蝕和縮小操作對比結(jié)果圖Fig.1 Comparison of erosion and contraction

圖1(c)和圖1 (d)對比發(fā)現(xiàn),3個像素中的中間像素經(jīng)過縮小操作的結(jié)果比腐蝕操作更為符合實際值。

圖2(a)表示了一維相位值中的相鄰三點,其相位模糊如圖2(b)所示,圖2(c)和圖2(d)分別顯示了膨脹操作和擴(kuò)大操作的結(jié)果。

圖2 膨脹和擴(kuò)大操作對比結(jié)果圖Fig.2 Comparison of dilation and expansion

圖2(c)和圖2(d)對比發(fā)現(xiàn),3個像素中的中間像素經(jīng)過擴(kuò)大操作的結(jié)果比膨脹操作更為接近真實值。因此,改進(jìn)的形態(tài)學(xué)開操作和閉操作可表示為

(5)

(6)

由于改進(jìn)的形態(tài)學(xué)開操作和閉操作基于擴(kuò)大和縮小操作,因此在相位濾波中進(jìn)行基于改進(jìn)的開閉操作濾波時,其結(jié)果比傳統(tǒng)的形態(tài)學(xué)開閉操作更為接近真實相位值,對相位的保護(hù)更好。

1.2 基于膨脹和腐蝕的重建操作

基于膨脹和腐蝕的形態(tài)學(xué)重建運算分別可表示[16]為

(7)

(8)

(9)

(10)

式中,Rd(φ′,φ)表示基于膨脹的形態(tài)學(xué)重建操作;δB(φ′)(x)表示膨脹操作;∧表示點方式的最小算子(與相同位置處φ值的最小值);Re(φ′,φ)表示基于腐蝕的形態(tài)學(xué)重建結(jié)果;εB(φ′)(x)表示腐蝕算子;∨表示點方式的最大算子(與相同位置處φ值的最大值);φ表示模板圖像,式(7)和式(9)中的k值分別表示關(guān)于φ的測地膨脹和測地腐蝕反復(fù)迭代達(dá)到穩(wěn)定值時的迭代次數(shù)。

基于膨脹的重建操作是利用結(jié)構(gòu)元素對φ′進(jìn)行膨脹操作并將結(jié)果圖像作為標(biāo)記圖像,然后取標(biāo)記圖像與原圖像中每個點對應(yīng)的最小值,以此進(jìn)行迭代直到圖像趨于穩(wěn)定,該操作能夠?qū)唿c噪聲有效地去除[14]。而基于腐蝕的重建操作用結(jié)構(gòu)元素對φ′進(jìn)行腐蝕操作并將結(jié)果圖像作為標(biāo)記圖像,然后取標(biāo)記圖像與原圖像每個點對應(yīng)的最大值,以此進(jìn)行迭代直到圖像趨于穩(wěn)定,該操作能夠彌合一些間斷的區(qū)域,保持條紋的連續(xù)性[14]。

2 改進(jìn)的自適應(yīng)梯度形態(tài)學(xué)重建濾波

本節(jié)對所提算法進(jìn)行詳盡的說明,首先是對原復(fù)數(shù)據(jù)進(jìn)行矢量化,接著估計相位梯度以確定形態(tài)學(xué)操作必需的線結(jié)構(gòu)元素,之后設(shè)計線結(jié)構(gòu)元素進(jìn)行基于改進(jìn)的灰度形態(tài)學(xué)開閉運算和膨脹與腐蝕的重建濾波操作,然后進(jìn)行公式反演得到濾波后的相位圖。

算法的流程圖如圖3所示。

圖3 算法流程圖Fig.3 Flow chart of the algorithm

2.1 復(fù)數(shù)據(jù)矢量化(正余弦變換)

為了避免復(fù)數(shù)干涉圖的強(qiáng)度信息(主輔圖像強(qiáng)度信息)對濾波精度和效果的影響,可以把干涉圖中的干涉相位進(jìn)行等權(quán)處理,采用相位矢量化(正余弦化)的處理方式。而MASF方法中則對原始相位進(jìn)行shifting操作(相位加π后對相位圖纏繞)作為第二幅圖像,兩幅相位圖像都存在相位的跳變,并沒有對相位起到有效的保護(hù)。本文相位正余弦化的實質(zhì)就是運用正余弦函數(shù)的周期性避免相位跳躍的影響,把相位圖像轉(zhuǎn)換為正弦圖像和余弦圖像,兩幅相位圖像不存在相位跳躍,有效保護(hù)相位條紋的連續(xù)性。后續(xù)的步驟在兩幅圖像(正、余弦圖像)中分別進(jìn)行操作。

2.2 中值濾波操作

在對相位進(jìn)行正余弦化后,分離開的正余弦數(shù)據(jù)仍有很多的噪聲,而本文的濾波方法中基于膨脹和腐蝕的形態(tài)學(xué)重建操作需要先確定線結(jié)構(gòu)元素,而噪聲的存在影響了梯度估計和線結(jié)構(gòu)元素的選取,因此有必要對正余弦圖像進(jìn)行預(yù)濾波操作,以增強(qiáng)梯度估計的準(zhǔn)確性。本文采用3×3大小的元素進(jìn)行中值濾波操作,既能夠?qū)D像進(jìn)行預(yù)濾波,也能有效保護(hù)相位條紋保護(hù)相位的連續(xù)性。

2.3 改進(jìn)的自適應(yīng)梯度估計

中值預(yù)濾波處理后,需要確定形態(tài)學(xué)處理所需的線結(jié)構(gòu)元素,而線結(jié)構(gòu)元素的選取由梯度估計得到。

對相位數(shù)據(jù)進(jìn)行分塊,在小塊里面進(jìn)行梯度估計,計算中心元素與窗口內(nèi)各點差值的絕對值,將差值絕對值最大的點與中心點作為梯度方向,為后續(xù)線結(jié)構(gòu)元素的選取提供依據(jù)。

由于在小塊里面求取差值最大點,此時差值絕對值最大的點并不一定在邊緣處,也可能在中間元素中(見圖4)。而MASF方法中對梯度值估計時,其在小塊邊緣查找差值絕對值的最大值,使得所選取的差值最大點在小塊的邊緣上,由于相位及條紋方向的不確定性,最大值點并不能保證完全在邊緣處,此時造成了梯度估計的不準(zhǔn)確性,線結(jié)構(gòu)元素為1的值也均在邊緣,使得在進(jìn)行形態(tài)學(xué)處理時容易造成小塊內(nèi)細(xì)節(jié)信息的丟失。而本方法求取的差值最大值點擺脫了該限制,差值最大點既可能在邊緣上也可能在小塊內(nèi),從而使得梯度估計更為準(zhǔn)確,后續(xù)步驟線結(jié)構(gòu)元素選取時不僅在邊緣處有值,在小塊里面也有為1的值,使得后續(xù)形態(tài)學(xué)濾波得到的細(xì)節(jié)信息更為豐富。

2.4 改進(jìn)的形態(tài)學(xué)處理

2.4.1 線結(jié)構(gòu)元素選取

當(dāng)形態(tài)學(xué)結(jié)構(gòu)元素為線結(jié)構(gòu)元素時,其條紋保持能力要好于矩形平坦結(jié)構(gòu)元素(矩形結(jié)構(gòu)元素更平滑一些)。因此，本文在進(jìn)行形態(tài)學(xué)處理時選取線結(jié)構(gòu)元素作為結(jié)構(gòu)元。文中選取垂直于梯度方向的元素點位置(不止有邊緣點)作為線結(jié)構(gòu)元素(值為1)。由于梯度估計時最大值點不一定在邊緣處,使得小塊內(nèi)部也存在線結(jié)構(gòu)元素的點(見圖4中線結(jié)構(gòu)元素L1在小塊內(nèi)部有兩個點作為線結(jié)構(gòu)元素的一部分),從而形態(tài)學(xué)濾波更為合理,有效保護(hù)相位條紋的細(xì)節(jié)信息。此外為了系統(tǒng)的魯棒性考慮(盡量避免噪聲的影響),將與線結(jié)構(gòu)元素距離最近的兩條線結(jié)構(gòu)元素也作為后續(xù)的處理。如圖4所示,本小塊內(nèi)的線結(jié)構(gòu)元素為L1、L2、L3。

圖4 梯度估計Fig.4 Estimation of gradients

2.4.2 改進(jìn)的形態(tài)學(xué)灰度開閉操作和重建操作的濾波

線結(jié)構(gòu)元素選取之后,進(jìn)入改進(jìn)的基于灰度形態(tài)學(xué)開閉運算和膨脹及腐蝕的重建操作。改進(jìn)的形態(tài)學(xué)開操作能夠平滑相位圖像的輪廓,并能消除細(xì)小的孤立區(qū)域和孤立噪聲像素;改進(jìn)的閉操作能夠彌合細(xì)小的空洞,填補(bǔ)條紋線的斷裂[16];基于膨脹的重建操作能夠?qū)唿c噪聲有效去除;而基于腐蝕的重建操作能夠彌合一些間斷的區(qū)域,保持細(xì)節(jié)的連續(xù)性[14]。

首先利用3個線結(jié)構(gòu)元素(L1、L2、L3)依次對圖像中選取的小塊進(jìn)行改進(jìn)的形態(tài)學(xué)開運算,選取3個小塊中對應(yīng)像素的最大值,將其作為標(biāo)記圖像,并將該步驟前的相位圖作為模板圖像,進(jìn)行基于膨脹的形態(tài)學(xué)重建操作。該步驟可用公式表示為

(11)

式中,φ1為正余弦化和中值濾波后的相位圖;χi(·)表示用線結(jié)構(gòu)元素L1、L2、L3分別進(jìn)行改進(jìn)的形態(tài)學(xué)開操作運算，i=1,2,3;Rd表示基于膨脹的重建操作。

其次分別利用3個線結(jié)構(gòu)元素(L1、L2、L3)對形態(tài)學(xué)膨脹重建操作后的圖像進(jìn)行改進(jìn)的形態(tài)學(xué)閉運算操作,選取3個小塊中對應(yīng)像素的最小值,將其作為標(biāo)記圖像,并將該步驟前的相位圖作為模板圖像,進(jìn)行基于腐蝕的形態(tài)學(xué)重建操作。該步驟可用公式表示為

(12)

式中,γi表示用線結(jié)構(gòu)元素L1、L2、L3分別進(jìn)行改進(jìn)的形態(tài)學(xué)閉操作運算(i=1,2,3);Re表示基于腐蝕的重建操作。

式(11)和式(12)類似于基于開(閉)操作的形態(tài)學(xué)重建(先腐蝕(膨脹)輸入圖像,后進(jìn)行基于膨脹(腐蝕))的重建)[16],其有效地保留了相位圖像腐蝕(膨脹)后留下的相位分量的形狀,減弱噪聲的影響,同時其利用改進(jìn)的形態(tài)學(xué)開(閉)運算和3個線結(jié)構(gòu)元素的優(yōu)化處理提高了相位濾波的魯棒性。

通過式(11)和式(12)的基于膨脹和腐蝕的形態(tài)學(xué)重建濾波操作,得到正余弦化后的相位濾波圖φsin和余弦化后的相位濾波圖φcos。

2.5 公式轉(zhuǎn)換

(13)

3 試驗結(jié)果分析

本文同時進(jìn)行仿真和實測數(shù)據(jù)處理試驗,并與中值圓周濾波[6]、Candeias方法[12]、MASF方法[14]方法進(jìn)行對比驗證所提方法的濾波效果。

3.1 仿真分析

文章將Matlab生成的peaks面作為相位仿真曲面來驗證所提算法的濾波效果。相位纏繞后的圖像如圖5(a)所示,數(shù)據(jù)大小為512像素×512像素。圖5(b)為帶噪聲的纏繞相位圖,圖中所加噪聲為1視相干系數(shù)為0.9、標(biāo)準(zhǔn)差為0.422 5的均勻噪聲[17]。上述4種方法分別對該干涉圖進(jìn)行濾波處理,相應(yīng)的濾波結(jié)果圖如圖5(c)～圖5(f)所示,統(tǒng)計信息(殘差點數(shù)目[18]和濾波前后的均方根誤差(root mean square error, RMSE))如表1所示。

表1 算法仿真實驗比較

從圖5可以看出MASF方法濾除效果較差,對噪聲的濾除較差,且對條紋造成了一定程度的破壞。中值圓周濾波和Candeias方法能濾除一部分噪聲,但仍有噪聲未能去除。而本文方法能有效濾除噪聲,且條紋保持性能較好。

為了定量分析上述方法的仿真濾波效果,現(xiàn)采用濾波前后的RMSE和殘差點數(shù)目[18]作為判斷依據(jù),各方法的統(tǒng)計結(jié)果如表1所示。

從表1可以看出本文方法濾波后的殘差點數(shù)目較少,與原始相位圖的均方根誤差也較少。

圖5 仿真實驗結(jié)果圖Fig.5 Results of simulation with different methods

仿真分析可以看出,本文方法的噪聲濾除效果較好,且條紋保持能力強(qiáng)于上述比較的方法。

3.2 實測數(shù)據(jù)處理及比較

為了驗證其在實測數(shù)據(jù)中的濾波效果,為此選取毫米波InSAR的數(shù)據(jù)進(jìn)行處理,并與其他方法進(jìn)行比較。圖6(a)所示為機(jī)載毫米波InSAR系統(tǒng)獲取的中國西南某地區(qū)的干涉條紋圖,大小為4 096像素×4 096像素。其相應(yīng)的相干系數(shù)圖如圖6(b)所示,圖6(b)右側(cè)坐標(biāo)代表相干系數(shù)的大小。從相干系數(shù)圖(計算窗口33×33)可以看出該相位數(shù)據(jù)有較多的低相干區(qū)域,對干涉處理造成較大的影響。4種方法的濾波結(jié)果如圖6(c)～圖6(f)所示。

圖6 實測實驗結(jié)果圖Fig.6 Results of filtering in real data

從圖6可以看出中值圓周濾波和Candeias方法對高相干區(qū)域有過度濾波的傾向,使得條紋相對變暗,條紋保持性能相對較差。而MASF方法對低相干區(qū)域的濾除效果相對不足,與本文方法相比,有更多的低相干區(qū)域并未濾除。而本文方法的噪聲濾除效果較好,既能在高相干區(qū)域保持較好的條紋,也能最大限度地濾除低相干區(qū)域的噪聲。但由于部分噪聲是陰影造成,為無效區(qū)域,無法進(jìn)行濾除,因此各方法均在一定區(qū)域都無法對噪聲進(jìn)行濾除。

為了定量分析上述方法的濾波效果,現(xiàn)采用殘差點數(shù)目[18]、濾波后平均偽相干值[19]來評價濾波效果,相關(guān)統(tǒng)計結(jié)果如表2所示。

表2 相關(guān)評價指標(biāo)對比結(jié)果

從表2中也可看出,MASF方法的平均偽相干系數(shù)最小,殘差點數(shù)目也最多,濾波并不充分;中值圓周濾波、Candeias濾波方法較為適中,其偽相干系數(shù)較高,殘差點數(shù)目也較少,而本文方法的濾除效果相對最好,其平均為相干系數(shù)最高,殘差點也最少,可以看出表中統(tǒng)計結(jié)果與圖6實測實驗結(jié)果圖相一致。

為了進(jìn)一步驗證其該方法在去噪及相位條紋保持方面的性能,選取第二組實測數(shù)據(jù)進(jìn)行實驗,并與其他方法進(jìn)行比較。圖7(a)所示為機(jī)載毫米波InSAR系統(tǒng)獲取的中國東南某地區(qū)的干涉條紋圖,大小為512像素×512像素,由于該系統(tǒng)基線較長,故條紋比較密集。其相應(yīng)的相干系數(shù)圖如圖7(b)所示,圖7(b)右側(cè)坐標(biāo)代表相干系數(shù)的大小。從相干系數(shù)圖(計算窗口9×9)可以看出該相位數(shù)據(jù)左側(cè)區(qū)域有較多的低相干區(qū)域,對干涉處理造成較大的影響。4種方法的濾波結(jié)果如圖7(c)～ 圖7(f)所示。

從圖7(e)可以看出MASF方法左側(cè)區(qū)域仍有部分噪聲尚未濾除,右側(cè)區(qū)域濾波效果較好。圖7(c)和圖7(d)中的兩種方法中值圓周濾波和Candeias方法在高相干區(qū)域濾波效果略好于MASF方法,而在低相干區(qū)域的濾波較差,使得不連續(xù)區(qū)域較多,條紋保持性能相對較差。而本文方法不僅在右側(cè)相干系數(shù)較好的區(qū)域噪聲濾除效果較好,條紋更為明亮清晰,在左側(cè)區(qū)域也能較好地實現(xiàn)濾波,同時細(xì)節(jié)信息也比較豐富,相位條紋保持較為良好。

圖7 實測實驗結(jié)果圖Fig.7 Filtering results in real data

為了定量分析上述方法的濾波效果,采用殘差點數(shù)目[18]、濾波后平均偽相干值[19]來進(jìn)行評價濾波效果,相關(guān)統(tǒng)計結(jié)果如表3所示。

從表3中也可看出,4種方法的平均偽相干值和殘差點數(shù)目預(yù)濾波前相比有明顯改善。其中,MASF方法的平均偽相干系數(shù)最小,殘差點數(shù)目也最多,濾波并不充分;中值圓周濾波、Candeias濾波方法較為適中,其偽相干系數(shù)較高,殘差點數(shù)目也較少,而本文方法的濾除效果相對最好,其平均偽相干系數(shù)最高,殘差點也最少,可以看出表中統(tǒng)計結(jié)果與圖7實測實驗結(jié)果圖相一致。

表3 相關(guān)評價指標(biāo)對比結(jié)果

兩組實測數(shù)據(jù)實驗分析可以看出,本文方法不僅能較好地濾除噪聲,且能有較好的條紋細(xì)節(jié)保持能力。

4 結(jié) 論

文中提出了一種改進(jìn)的自適應(yīng)梯度形態(tài)學(xué)重建濾波算法。算法基于改進(jìn)的灰度形態(tài)學(xué)開閉運算和基于膨脹與腐蝕重建濾波操作,利用復(fù)數(shù)據(jù)的矢量化和改進(jìn)的梯度估計有效地保護(hù)相位的跳躍特性及連續(xù)性,改進(jìn)的自適應(yīng)的線結(jié)構(gòu)元素使得細(xì)節(jié)信息不被破壞。仿真分析和實測數(shù)據(jù)對比處理表明該方法的濾波效果較好,細(xì)節(jié)表現(xiàn)力也較強(qiáng)。