基于局部條紋特征的頻率估計方法

合成孔徑干擾測量技術是利用合成孔內數據提取的相位信息為信息源獲得三維信息和變化信息的技術。干涉合成孔徑雷達/聲納通過兩幅接收陣獲取地表同一場景的復圖像對,經共軛復乘后得到干涉圖。干涉圖中包含了斜距向上的點到兩接收陣之間精確的距離差信息。但是,由于受系統噪聲、去相干、配準誤差和地表變化等因素的影響,實際獲得的干涉圖會存在不同程度的噪聲,掩蓋了真實值的同時,也影響了后續相位解纏的速度和精度。為了降低噪聲的影響,需要對干涉圖進行濾波。目前干涉圖濾波方法有很多,濾波可以在空域或頻域進行。空域濾波主要有圓周期均值濾波,圓周期中值濾波和Lee濾波等。空域濾波器實現簡單,運算速度快,但容易破壞陡峭地形產生的高頻干涉條紋。頻域濾波的代表為Goldstein濾波,該算法利用干涉圖自身的頻譜冪指數進行頻域加權實現濾波,對相位噪聲較小的干涉圖具有良好的降噪性能,但在相位噪聲大的區域濾波效果下降。無論是空域還是頻域濾波,都很少具有根據局部條紋特征調整濾波器的能力。近些年,一些文獻提出利用局部頻率估計的空頻域自適應濾波算法,相比傳統方法,這類算法對局部區域的濾波更有針對性,能夠在保持干涉條紋的前提下更好地去除相位噪聲。但存在的問題是:濾波性能依賴局部信噪比的大小。當信噪比較低時,可能產生錯誤的局部頻率估計結果,進而致使相位濾波出現錯誤。為了解決局部頻率估計類算法存在的問題,本文提出了一種改進的干涉圖濾波方法。其主要思想是:在提取局部條紋方向和頻率特征的基礎上,利用方位和頻率平穩的特征,對二維頻率估計結果進行分析和判定,對認定為頻率估計出現錯誤的點,通過適當調整估計窗口改善信噪比,實現對錯誤頻率的修正,然后利用修正后的頻率作為濾波器的參數完成濾波。1反映地形輪廓相位t-c假設單視主輔復圖像s1和s2為s1=C+n1s2=Ce-j?Τ+n2(1)s1=C+n1s2=Ce?j?T+n2(1)式中:C為主輔圖像的公共部分,n1和n2為熱噪聲,?T為真實的干涉相位。s1和s2經共軛相乘后得到的干涉信號可以表示為I=s1s*2=|C|2ej?T+Cn*2+C*n1ej?T+n1n*2(2)式(2)可以進一步寫為I=ej?c{|C|2ej(?T-?c)+e-j?cEJ?T+N1N*2}(3)式中:?c為反映地形輪廓的相位,?T-?c為反映地形細節的相位,方括號內的項可視為噪聲。根據式(3)濾波過程可分兩步進行,首先提取反映地形輪廓的相位,然后對剩余地形細節相位進行提取。?c在局部區域內(窗口)可以近似認為線性不變,利用平面近似,即ej?c≈ej(4)式中:m和n表示相對窗口中心點(x,y)的位置偏移,?0為初始相位,可根據文獻中的方法進行估計;fx和fy為點(x,y)對應的二維空間頻率。在提取了地形輪廓對應的干涉相位后,再對細節相位?T-?c進行提取,可以通過計算相位均值的方法獲得。反映地形輪廓相位?c的提取對于該濾波算法非常重要,因此二維頻率fx和fy的參數估計就成為了濾波的關鍵。最常用的二維頻率估計方法是最大似然法,它通過二維信號周期圖的峰值位置來估計二維頻率fx和fy,可以采用二維離散傅立葉變換(DiscreteFouriertransform,DFT)實現[?fx,?fy]=maxfx,fy{|Μ∑m=1Ν∑n=1Ι(x,y)e-j(mfx+nfy)|}(5)[f?x,f?y]=maxfx,fy{∣∣∣∣∑m=1M∑n=1NI(x,y)e?j(mfx+nfy)∣∣∣∣}(5)式中:M×N為選取的樣本窗口尺寸,?ff?x和?ff?y分別為fx和fy的估計值。為了提高估計精度,需要對二維DFT進行插值。插值可以采用補零實現,但由于圖像尺寸通常很大,直接補零插值計算量巨大。文獻提出先利用DFT對局部條紋頻率粗估計,然后再采用線性調頻Z變換(Chirp-Ztransform,CZT)在粗估計值附近細化頻譜的方法,大大減少了二維頻率估計所需的計算量。2干涉圖信噪比較低時頻率估計錯誤的發生局部頻率估計干涉圖濾波利用局部相位特征進行濾波,濾波時很有針對性,濾波能力強。但不足的是,由于該濾波方法過于依賴頻率估計結果,一旦局部頻率估計出現較大誤差或錯誤,就會破壞原有的干涉圖,或在干涉圖中引入額外的奇異點。局部頻率估計的錯誤主要是由噪聲引起的,當局部區域干涉圖信噪比較低時,可能出現的錯誤情況有兩種:一是頻譜中相位噪聲的峰值超過真實條紋所對應的頻率峰值,如圖1(a)所示。在采用頻譜峰值位置估計條紋頻率時,錯誤地選擇了噪聲峰值。另一種頻率估計錯誤是由DFT頻率粗估計產生的,如圖1(b)所示。雖然干涉條紋對應的峰值比相位噪聲的峰值高,但超過幅度并不大,當DFT的點數較少時,就可能出現因頻率采樣點少,導致頻率采樣后的噪聲頻率峰值相比條紋頻率峰值高,產生頻率粗估計錯誤,從而導致后續的CZT頻譜細化的是噪聲頻率。雖然干涉圖信噪比較低時容易發生頻率估計錯誤,但并非完全無法克服。條紋信號的頻率峰值仍然存在,只是被噪聲峰值掩蓋了。只要能夠對局部頻率估計錯誤值進行判定及識別,對頻譜估計值進行修正,仍然可以得到精確的濾波結果。3條紋方向和條紋頻率干涉條紋可以看成一種方向流場,并且其局部可認為是平穩的,這一特征與激光干涉條紋和指紋條紋的局部特征非常相似,可采用類似方法對局部干涉條紋進行建模。局部區域內的干涉條紋可以用個兩個參數來描述:條紋方向和條紋頻率。其中方向性是局部干涉條紋最重要的特征。相對于干涉條紋的頻率,濾波對干涉條紋的方向更加敏感。若條紋方向錯誤,即使條紋頻率估計正確,得到的濾波結果也將完全偏離真實值。條紋頻率是局部干涉條紋重要的信息,決定濾波后相位誤差大小。下文給出條紋方向和條紋頻率的定義。假設圖像中任意點(x,y)的相位為?(x,y),則干涉條紋方向定義為相位場方向導數為零的方向,即???l=???xcosθ′+???ysinθ′=Gxcosθ′+Gysinθ′=0(6)???l=???xcosθ′+???ysinθ′=Gxcosθ′+Gysinθ′=0(6)式中:θ′為l到x軸的轉角,稱為條紋的方向角,求解后可寫為θ′=arctan(Gy/Gx)-π/2(7)需要注意的是θ′取值范圍為[0,π),而實際條紋方向取值范圍應為),經修正后的方向角θ可表示為【MATH244Z】由于二維頻率表示為【MATH245Z】因此方向圖估計也可以用二維頻率估計結果表示,即【MATH246Z】條紋的頻率定義為沿條紋方向的相位變化率,可以直接計算二維頻率矢量的模值得到頻率估計【MATH247Z】4維頻率估計的判別進行頻率估計修正,首先需要對頻率估計錯誤的像素點進行判定。由于局部區域干涉條紋是平穩的,條紋頻率的變化比較平緩,局部區域干涉條紋方向和頻率不會存在大的跳變。因此,可以通過考察相鄰像素點的方向或頻率的連續性,對當前像素點二維頻率估計的可信度進行判斷。可信度可以分為方向估計可信度和頻率估計可信度,分別從方向圖和頻率圖中獲得。圖2給出了一組仿真數據的干涉相位圖、方向圖和頻率圖。其中,圖2(a)為真實相位為圓錐體的干涉相位圖,信噪比為2dB,圖2(b)和(c)分別為利用式(8)和(11)計算得到的方向圖和頻率圖,部分像素點方向和頻率估計存在明顯的不連續,需要進行修正。本文以條紋方向和頻率不連續程度作為測度提出一種可信度評價標準。假設像素點(x,y)的條紋方向和頻率分別為θ(x,y)和f(x,y),則其方向和頻率估計可信度Cd(x,y)和Cf(x,y)可表示為Cd(x,y)=L∑u=1Κ∑v=1∑u=1L∑v=1K|arg{exp[j?θ(x,y)-j?θ(x-u,y-v)]}|(12)Cf(x,y)=L∑u=1Κ∑v=1|?f(x,y)-?f(x-u,y-v)|(13)式中:u和v表示相對中心點(x,y)的位置偏移,L×K為相鄰樣本點數。當干涉圖信噪比較高時,局部區域的干涉信號近似平穩,相鄰點的條紋方向和頻率基本保持不變,Cd和Cf值小,方向和頻率估計可信度高。當干涉圖信噪比較低時,局部區域的非平穩性增強,Cd和Cf值增大,方向和頻率估計可信度降低。由于有可能出現頻率相同方向不同或方向相同頻率不同的情況,為了對頻率估計錯誤點判定更加準確,可信度CA和Cf之間需要進行互補,這時干涉圖中需要修正的區域R可以表示為R=Rd∪Rf(14)Rd={(x,y)|Cd(x,y)>DTd}(15)Rf={(x,y)|Cf(x,y)>DTf}(16)需修正區域R和非需修正區域ˉR可以視為二模態的,因此,最佳閾值DTd和DTf可利用最大類別方差法確定。下文從更直觀的角度,以當前像素點及其四鄰域點為例說明二維頻率估計的判定過程。如圖3(a)所示,中心灰色A0為待修正點,其鄰點分別為A1、A2、A3和A4。當信噪比較高時,各點的二維頻率估計是準確的,方向和頻率變化都很小,A0頻率估計的可信度小于閾值,無需修正。當信噪比較低時,由于局部噪聲是非平穩的,空間位置不同點的噪聲頻率峰值也是不確定的,因此必然會導致A0與4個相鄰點的頻率估計值相關性變差,可信度大于閾值,需要進行頻率修正。頻率修正采用增大空間估計窗口長度的辦法,通過將估計窗口由W擴大為W′,提高窗口內干涉條紋信噪比。需要注意的是,W′尺寸選擇不能過大,以免破壞干涉圖局部區域平穩的假設,同時也要確保窗口內信噪比能明顯提高。另外,為減小由于窗口增大可能引入的頻率估計誤差的影響,可采取以下步驟:將W′窗口的頻率估計結果(?fx′,?fy′)作為初步頻率估計。在W窗口內的信號二維頻譜中,搜索(?fx′,?fy′)附近的頻譜峰值,如圖3(b)所示,中間虛線框內為峰值搜索窗口,將最接近(?fx′,?fy′)的峰值所對應的頻率作為最終估計頻率(?fx,?fy)。若W窗口內所有的頻譜峰值離(?fx′,?fy′)都比較遠,則采用(?fx′,?fy′)作為A0區域的頻率估計值。5研究結果與分析為了驗證方法的濾波性能,本文采用幾種經典的濾波算法進行比較。分別為圓周期均值濾波、Goldstein濾波和基于局部頻率估計的地形自適應濾波。采用仿真數據和實際數據分別進行驗證。仿真的理想干涉圖如圖4(a)所示。為了比較各算法在不同信噪比下的濾波性能,將干涉圖分為上下兩部分加入相位噪聲,上半部分信噪比為9dB,下半部分信噪比為0.5dB,加入噪聲后的干涉圖如圖4(b)所示。各濾波算法的樣本窗口大小均為11×11,Goldstein濾波算法中的冪指數α取值0.5。圖4(c)～(f)依次給出了采用圓周期均值濾波、Goldstein濾波、基于局部頻率估計的地形自適應濾波以及本文方法的處理結果。由圖4可見,在兩種不同信噪比下,圓周期均值濾波對四周比較稀疏的條紋濾波可以取得比較滿意的結果,但是中間條紋密集的區域效果不好,即使噪聲較小的區域也影響了干涉條紋的致密性。Goldstein算法只在高信噪比區域有較好的濾波效果。基于局部頻率估計的地形自適應濾波方法在低信噪比條件下,濾波性能較前兩種方法有較大改善,但仍然存在明顯的小塊相位噪聲,這是由于頻率估計錯誤所致,采用本文方法可獲得更優的濾波效果。為進一步檢驗本文方法的可靠性和有效性,采用航天飛機合成孔徑雷達X-SAR對意大利Enta火山的實測數據進行驗證,處理結果如圖5所示。各濾波算法的樣本窗口大小均為15×15,Goldstein濾波算法中的冪指數α取值0.9。從圖5可以看出,圓周期均值濾波雖然可以顯著去除相位噪聲,但是在干涉條紋密集區域其致密性受到嚴重影響,Goldstein算法濾波后仍存在大量相位噪聲,局部頻率估計方法相比前兩種算法,處理性能有明顯改善,但在中間區域還是出現部分明顯錯誤。而采用本文方法處理得到的干