改進的形態學干涉圖濾波方法

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

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

0 引 言

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

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

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

1 數學形態學的相位濾波

在干涉圖濾波方法中,數學形態學能夠保持相位的跳變特性和相關紋理,其利用結構元素在相位圖像上不斷移動來確定相位之間的關系,通過其幾何特征濾除疊加在相位圖像上的噪聲。本文算法中在數學形態學中的應用算法主要為改進的形態學開閉操作和基于膨脹和腐蝕的重建操作。

1.1 改進的形態學開閉操作

(1)

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

(2)

在相位濾波中,由于相位值在[-π,π]的限制,如果直接使用膨脹和腐蝕操作進行濾波,會造成局部濾波結果不合理。而Candeias[12]對膨脹和腐蝕操作進行了修改,提出了擴大和縮小操作,使得其在進行運算時更符合真實值。擴大和縮小操作定義[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)分別顯示了腐蝕操作和縮小操作的結果。

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

圖1(c)和圖1 (d)對比發現,3個像素中的中間像素經過縮小操作的結果比腐蝕操作更為符合實際值。

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

圖2 膨脹和擴大操作對比結果圖Fig.2 Comparison of dilation and expansion

圖2(c)和圖2(d)對比發現,3個像素中的中間像素經過擴大操作的結果比膨脹操作更為接近真實值。因此,改進的形態學開操作和閉操作可表示為

(5)

(6)

由于改進的形態學開操作和閉操作基于擴大和縮小操作,因此在相位濾波中進行基于改進的開閉操作濾波時,其結果比傳統的形態學開閉操作更為接近真實相位值,對相位的保護更好。

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

基于膨脹和腐蝕的形態學重建運算分別可表示[16]為

(7)

(8)

(9)

(10)

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

基于膨脹的重建操作是利用結構元素對φ′進行膨脹操作并將結果圖像作為標記圖像,然后取標記圖像與原圖像中每個點對應的最小值,以此進行迭代直到圖像趨于穩定,該操作能夠將斑點噪聲有效地去除[14]。而基于腐蝕的重建操作用結構元素對φ′進行腐蝕操作并將結果圖像作為標記圖像,然后取標記圖像與原圖像每個點對應的最大值,以此進行迭代直到圖像趨于穩定,該操作能夠彌合一些間斷的區域,保持條紋的連續性[14]。

2 改進的自適應梯度形態學重建濾波

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

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

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

2.1 復數據矢量化(正余弦變換)

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

2.2 中值濾波操作

在對相位進行正余弦化后,分離開的正余弦數據仍有很多的噪聲,而本文的濾波方法中基于膨脹和腐蝕的形態學重建操作需要先確定線結構元素,而噪聲的存在影響了梯度估計和線結構元素的選取,因此有必要對正余弦圖像進行預濾波操作,以增強梯度估計的準確性。本文采用3×3大小的元素進行中值濾波操作,既能夠對圖像進行預濾波,也能有效保護相位條紋保護相位的連續性。

2.3 改進的自適應梯度估計

中值預濾波處理后,需要確定形態學處理所需的線結構元素,而線結構元素的選取由梯度估計得到。

對相位數據進行分塊,在小塊里面進行梯度估計,計算中心元素與窗口內各點差值的絕對值,將差值絕對值最大的點與中心點作為梯度方向,為后續線結構元素的選取提供依據。

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

2.4 改進的形態學處理

2.4.1 線結構元素選取

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

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

2.4.2 改進的形態學灰度開閉操作和重建操作的濾波

線結構元素選取之后,進入改進的基于灰度形態學開閉運算和膨脹及腐蝕的重建操作。改進的形態學開操作能夠平滑相位圖像的輪廓,并能消除細小的孤立區域和孤立噪聲像素;改進的閉操作能夠彌合細小的空洞,填補條紋線的斷裂[16];基于膨脹的重建操作能夠將斑點噪聲有效去除;而基于腐蝕的重建操作能夠彌合一些間斷的區域,保持細節的連續性[14]。

首先利用3個線結構元素(L1、L2、L3)依次對圖像中選取的小塊進行改進的形態學開運算,選取3個小塊中對應像素的最大值,將其作為標記圖像,并將該步驟前的相位圖作為模板圖像,進行基于膨脹的形態學重建操作。該步驟可用公式表示為

(11)

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

其次分別利用3個線結構元素(L1、L2、L3)對形態學膨脹重建操作后的圖像進行改進的形態學閉運算操作,選取3個小塊中對應像素的最小值,將其作為標記圖像,并將該步驟前的相位圖作為模板圖像,進行基于腐蝕的形態學重建操作。該步驟可用公式表示為

(12)

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

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

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

2.5 公式轉換

(13)

3 試驗結果分析

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

3.1 仿真分析

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

表1 算法仿真實驗比較

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

為了定量分析上述方法的仿真濾波效果,現采用濾波前后的RMSE和殘差點數目[18]作為判斷依據,各方法的統計結果如表1所示。

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

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

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

3.2 實測數據處理及比較

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

圖6 實測實驗結果圖Fig.6 Results of filtering in real data

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

為了定量分析上述方法的濾波效果,現采用殘差點數目[18]、濾波后平均偽相干值[19]來評價濾波效果,相關統計結果如表2所示。

表2 相關評價指標對比結果

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

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

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

圖7 實測實驗結果圖Fig.7 Filtering results in real data

為了定量分析上述方法的濾波效果,采用殘差點數目[18]、濾波后平均偽相干值[19]來進行評價濾波效果,相關統計結果如表3所示。

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

表3 相關評價指標對比結果

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

4 結 論

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