陳冬,王江安,康圣
(海軍工程大學(xué)電子工程學(xué)院光電研究所,湖北武漢 430033)
在激光遙感領(lǐng)域,研究者的注意力主要集中在信號反演方法上而忽視了降噪處理算法的深入研究。通常采用的僅僅是中值濾波、窄帶濾波等傳統(tǒng)方法,也有研究者采用多次測量取均值法。近年才有研究者應(yīng)用小波閾值法進(jìn)行濾波并取得了較好效果[1-2]。
對脈沖激光雷達(dá)信號頻譜進(jìn)行分析可知,其近場信號頻譜成分分布廣泛而遠(yuǎn)場高頻分量能量較低而低頻分量起主導(dǎo)作用。因而較之傳統(tǒng)的濾波方法,時頻分析工具等現(xiàn)代信號處理方法更適合于該課題。經(jīng)驗?zāi)B(tài)分解理論是美籍華人黃鄂教授提出的一種基于時頻分析的新算法,它與小波分解算法的區(qū)別在于對瞬時頻率定義的差別。本文采用經(jīng)驗?zāi)B(tài)分解(EMD)濾波對激光雷達(dá)信號進(jìn)行降噪處理,并與擴(kuò)展卡爾曼算法(EKF)和小波濾波進(jìn)行比較。將3種方法在不同大氣條件下的處理結(jié)果進(jìn)行對比,從而得到最優(yōu)的降噪方法。
仿真驗證時如果采用脈沖激光雷達(dá)回波實測數(shù)據(jù),無法對原始數(shù)據(jù)和濾波后數(shù)據(jù)的信號進(jìn)行評價,因而必須構(gòu)建仿真理想信號。嚴(yán)格地說,每種大氣條件下的信號模型都是不同的。它受到大氣粒子尺寸分布,粒子形狀分布、濃度,粒子性質(zhì)及激光雷達(dá)發(fā)射和接收視場角等多種因素影響。研究者通常把大氣分為3種情況進(jìn)行建模,把非常稀薄和非常稠密的粒子群分布單獨考慮,分別采用單粒子散射近似和漫反射近似,而最普遍的情況是采用一級多次散射近似進(jìn)行建模。大氣脈沖激光雷達(dá)方程(式(1))正是在一級多次散射近似條件下推導(dǎo)出來的,適用于絕大多數(shù)大氣回波情況,也與實測信號最為接近。
式中:Pr(R)為距離R處的回波功率;βπ(R)為后向散射系數(shù);σ(R)為消光系數(shù);Cq為儀器常數(shù),它又可以表示為:
式中:Ar為接收鏡頭面積;Tt和Tr分別為發(fā)射/接收光學(xué)系統(tǒng)透過率;E0為發(fā)射脈沖能量;τp為脈沖寬度;c為光速;η(R)為幾何重疊因子。
后向散射系數(shù)是激光雷達(dá)方程中的重要參數(shù),根據(jù)Klett的研究[3],在米散射起主導(dǎo)作用時,消光系數(shù)與后向散射系數(shù)的關(guān)系為β(z)=B0σk(z),B0為后向散射比,k的取值和波長與粒子尺度的相對值有關(guān),通常取值區(qū)間為[0.66,1]。
實際的激光雷達(dá)接收信號要通過光電探測器、高頻放大器、模數(shù)轉(zhuǎn)換器等器件。根據(jù)實際器件參數(shù),將信號模型進(jìn)一步修正。設(shè)放大器的偏置電平為lv,AD轉(zhuǎn)換器的采樣頻率為f,同時假定第n-1個采樣點與第n個采樣點范圍內(nèi)采樣體積內(nèi)的平均消光系數(shù)為,回波模型可以修正為:
圖1 仿真、實測回波信號比較Fig.1Compare between real and simulated signal
為了檢驗?zāi)P偷恼_性,將真實信號與同等條件仿真波形進(jìn)行比較。圖1中虛線為某型激光雷達(dá)在某機(jī)場晴朗均勻大氣中的實測回波信號,實線是采用式(3)擬合的理想回波信號??梢钥闯?,仿真信號對遠(yuǎn)場回波信號擬合較好,近場信號存在差別,這是幾何重疊因子的差異造成的。儀器常數(shù)的諸多因子中,幾何重疊因子是難以確定的。有研究人員假設(shè)光束能量高斯分布通過激光光束發(fā)射與接收視場截面的大小來確定幾何重疊因子的大?。?]。但激光光束能量的分布特性的隨機(jī)性大,近場能量強(qiáng)、多次散射強(qiáng)烈使得這種方法意義不大。因此,研究中只關(guān)注視場完全重疊后的區(qū)域,近場擬合的精確性對仿真效果并無影響。
激光雷達(dá)的噪聲按形成原因可分為散粒噪聲、熱噪聲、背景光噪聲及1/f噪聲等[2]。其中除了1/f噪聲外,其他噪聲都可認(rèn)為是加性噪聲,并且符合高斯分布。因此研究中采用高斯噪聲模型。
評價信號質(zhì)量通常采用的參數(shù)是信噪比,但信號信噪比的定義有很多,本文選用信號幅值均方根與噪聲幅度均方根之比來定義。
式中:Pi為信號幅度;Ni為噪聲幅度;n為信號長度。
在激光雷達(dá)領(lǐng)域,研究者將擴(kuò)展卡爾曼算法當(dāng)作1種大氣狀態(tài)參數(shù)的反演方法[5-6]。這樣難以對算法降噪性能加以評估,也無法驗證反演模型的準(zhǔn)確性。因此筆者傾向于把EKF法修正為1種降噪算法。
將消光系數(shù)作為信號主要參量,信號模型為
其中:v(n)為消光系數(shù)波動參數(shù),它服從正態(tài)分布,方差為R,均值為0。
偏置電平很容易濾除,去掉放大器偏置電平并考慮加性噪聲,式(3)可化簡為式(6),其中w(n)均值為0,方差為Q。
初始條件為P(n0)=0,n0為發(fā)射與接收視場的完全重疊位置,此時η(n0)=1。
小波濾波采用多尺度正交小波基對信號進(jìn)行分解,將信號和噪聲在不同時頻域進(jìn)行分離,然后在各時頻范圍內(nèi)對依據(jù)閾值進(jìn)行濾波,最后重構(gòu)信號從而達(dá)到降噪的目的。離散小波基的定義為:
任何函數(shù)都可以分解為小波基的加權(quán)和,
其中,加權(quán)值即為小波系數(shù):
小波理論產(chǎn)生于20世紀(jì)80年代,學(xué)者們對小波降噪方法的研究廣泛且深入,因此本文不再討論。對于激光雷達(dá)信號降噪問題,小波基和閾值的選取是最重要的2個方面。不同小波基的特性不同,選取不同的小波基對濾波結(jié)果影響很大。Db5小波具有雙正交性、緊支撐性和近似對稱性,十分適合對指數(shù)衰減型曲線進(jìn)行濾波。閾值的選取對濾波效果影響也很大,常用的降噪閾值選取方法有長度對數(shù)閾值(Sqtwolog)、最小極大方差閾值(Minimaxi)、無偏風(fēng)險閾值(Rigrsure)、啟發(fā)式閾值(Heursure)等[8-10]。這些方法的主要思想是根據(jù)小波基長度與噪聲均方根構(gòu)造最佳閾值,它們都可以改善信號信噪比,但降噪能力有限。對于本課題來說,仿真表明固定閾值降噪效果更好,當(dāng)閾值選擇為二倍噪聲均方根時可使濾波效果達(dá)到最佳。
EMD分解濾波也是一種基于時頻分析的算法,它與小波分解算法的區(qū)別在于對瞬時頻率定義的差別。EMD算法具有自適應(yīng)性,它會為每條信號曲線產(chǎn)生不同的正交基。
經(jīng)驗?zāi)B(tài)分解(EMD),是通過固有模態(tài)函數(shù)(IMF)提取來實現(xiàn)的[11]。1個函數(shù)可以分解為n個固有模態(tài)函數(shù)ci(t)與1個剩余分量r(t)的和,即
固有模態(tài)函數(shù)IMF必須滿足以下2個條件:一是信號持續(xù)時間內(nèi),其過0點的個數(shù)和極值點的個數(shù)相等或最多相差1個;二是信號極大值或極小值構(gòu)成的包絡(luò)線的均值為0。
固有模態(tài)函數(shù)采用以下步驟進(jìn)行提取:尋找信號f(t)的局部極值點;在極大值之間進(jìn)行插值擬合處理形成最大值包絡(luò)曲線maxi(t),在極小值之間進(jìn)行插值擬合處理形成最小值包絡(luò)曲線mini(t);求包絡(luò)均值曲線m(t)=[max(t)+min(t)]/2;信號減去包絡(luò)均值曲線;多次迭代后,剩余信號滿足停止準(zhǔn)則時篩選結(jié)束。篩選門限見式(15),篩選門限一般在0.2~0.3之間。
經(jīng)過EMD提取之后,信號被分解為多個固有模態(tài)函數(shù),對固有函數(shù)進(jìn)行閾值濾波,將濾波后的數(shù)據(jù)重構(gòu)就可以達(dá)到降噪的目的。
仿真表明,篩選門限取值2.5時濾波效果較好,軟閾值法的效果比硬閾值法濾波效果好,當(dāng)閾值為噪聲均方根時可使降噪效果達(dá)到最佳。
仿真主要針對2種大氣條件下的回波信號:均勻大氣和局部突變大氣。均勻大氣主要針對無風(fēng)或微風(fēng)條件下的自然大氣狀態(tài),這時消光系數(shù)比較穩(wěn)定隨距離起伏不大。局部突變大氣針對存在風(fēng)切變、湍流及云霧或煙塵等天氣條件,它們會改變大氣的局部分布,造成局部消光系數(shù)突變。
常規(guī)大氣條件下,消光系數(shù)均勻,在均值附近上下起伏。因而選擇不同的消光系數(shù)均值和方差構(gòu)造激光脈沖回波,并疊加高斯噪聲。如圖2和圖3所示。
仿真依據(jù)消光系數(shù)從低至高進(jìn)行,每個消光系數(shù)取值又采用不同的噪聲幅度進(jìn)行仿真,仿真結(jié)果如表2所示。
表2中第1行為仿真消光系數(shù)值,第1列為噪聲均方根(δ)。表中數(shù)據(jù)依次為原始信號EKF濾波小波濾波EMD濾波對應(yīng)的信噪比。
此種條件時,絕大部分消光系數(shù)分布與均勻大氣下相同,只是局部空氣中消光系數(shù)陡增驟減。當(dāng)消光系數(shù)如圖4所示分布時,回波信號波形如圖5所示,仿真結(jié)果如表3所示。
表3中第1行為仿真消光系數(shù)值,第1列為噪聲均方根(δ)。表中數(shù)據(jù)依次為原始信號EKF濾波小波濾波EMD濾波對應(yīng)的信噪比。
擴(kuò)展卡爾曼濾波的降噪能力有限,但可以將信號直接反演為消光系數(shù),初始值的選取對這種方法的濾波效果影響很大。消光系數(shù)較高時,EKF濾波后的失真較大;在消光系數(shù)較低時的濾波效果較好。在均勻大氣條件下,EKF法反演精度較高;而當(dāng)大氣中有局部突變時,若消光系數(shù)陡增幅度過大會破壞這種方法的成立條件,導(dǎo)致誤差較大。小波濾波的效果在多數(shù)情況下是幾種方法中的最佳方法,信噪比較高時,小波濾波的優(yōu)勢并不明顯,并且有可能造成信號失真;中低信噪比條件下,這種方法對信號質(zhì)量的改善最為明顯;在突變大氣條件下,盡管這種方法的濾波效果降低了,但仍然是最佳方法。在信噪比較高時EMD濾波效果明顯好于其他2種方法。目前采用EMD法對激光雷達(dá)回波信號進(jìn)行濾波的研究仍處于起步階段,在降噪閾值和篩選門限這2個方面進(jìn)一步研究會使該方法的優(yōu)勢得到進(jìn)一步體現(xiàn)。
綜上所述,本文根據(jù)激光雷達(dá)性能參數(shù)和典型的大氣條件建立了1種更為實用的回波信號仿真模型。本文工作證明EMD濾波法在某些大氣條件下可以達(dá)到最佳濾波效果,是現(xiàn)有方法的必要補充。綜合采用3種方法才能使濾波達(dá)到最佳效果。均勻大氣條件下,如果噪聲幅度較高,應(yīng)采用小波濾波;噪聲幅度較低時,如果消光系數(shù)較小,應(yīng)采用EKF濾波;如果消光系數(shù)較大,就應(yīng)采用EMD算法。突變大氣條件下,如果噪聲幅度較高,宜采用小波濾波;如果噪聲幅度較低則應(yīng)采用EMD濾波。將該方法應(yīng)用到激光雷達(dá)消光系數(shù)測量工程中,得到了較好的效果。
[1]FANG Hai-tao,HUANG De-shuang.Noise reduction in lidar signal based on discrete wavelet transform[J].Optic communication,2004,233:67-76.
[2]Aime Lay-Ekuakille,TROTTAA.Comparisonbetween adjustablewindowtechniqueandwaveletmethodin processing backscattering lidar sinal[C].SPIE,2004 (5240):72-82.
[3]KLETT J D.Stableanalyticalinversionsolutionfor processing lidar returns[J].Applied Optics.1981,20(2): 211-210.
[4]王治華,賀應(yīng)紅,左浩毅,等.基于高斯光束特性的Mie散射大氣激光雷達(dá)回波近場信號校正研究[J].物理學(xué)報,2006,55(6):3188-3192.
WAN Zhi-hua,HE Ying-hong,ZUO Hao-yi,et al.The correction of short-range Mie scattering laser lidar returns based on the Gaussian character of laser beam[J].Acta Physica Sinica,2006,55(6):3188-3192.
[5]LAINIOTIS D G,PAPAPARASKEVA P,KOTHAPALLI G,et al.Adaptive filter applications to LIDAR:return power and log power estimation[J].IEEE Transactons on Geoscience and Remote Sensing,1996,34(4):886-891.
[6]José M.Bioucas Dias,et al.Reconstruction of backscatter and extinction coeffictents in lidar:A stochastic filtering approach[J].IEEE Transactions on Geoscience and Remote Sensing,2004,42(2):443-456.
[7]HUANG N E,ZHENG Shen,STEVEN R L,et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J].The Royal Society London A,1998,454:903-995.
[8]張斌,王彤,谷傳綱,等.改進(jìn)的小波閾值消噪法在湍流信號處理中的應(yīng)用[J].工程熱物理學(xué)報,2009,30(3): 407-410.
ZHANG Bin,WANG Tong,GU Chuan-gang,et al.Applcation of improved wavelet threshold de-noising in turbulent singnal processing[J].Journal of Engineering Thermophysics,2009,30(3):407-410.
[9]DONOHO D L.Denoising by soft-thresholding[J].IEEE Trans on Information Theory,1995,41(3):613-627.
[10]DONOHO D L,JOHNSTONE I M.Ideal spatial adaptation by wavelet shrinkage[J].Journal of the Royal Statistical Society,1997,59(2):319-351.