陳冬,王江安,康圣
(海軍工程大學電子工程學院光電研究所,湖北武漢 430033)
在激光遙感領域,研究者的注意力主要集中在信號反演方法上而忽視了降噪處理算法的深入研究。通常采用的僅僅是中值濾波、窄帶濾波等傳統方法,也有研究者采用多次測量取均值法。近年才有研究者應用小波閾值法進行濾波并取得了較好效果[1-2]。
對脈沖激光雷達信號頻譜進行分析可知,其近場信號頻譜成分分布廣泛而遠場高頻分量能量較低而低頻分量起主導作用。因而較之傳統的濾波方法,時頻分析工具等現代信號處理方法更適合于該課題。經驗模態分解理論是美籍華人黃鄂教授提出的一種基于時頻分析的新算法,它與小波分解算法的區別在于對瞬時頻率定義的差別。本文采用經驗模態分解(EMD)濾波對激光雷達信號進行降噪處理,并與擴展卡爾曼算法(EKF)和小波濾波進行比較。將3種方法在不同大氣條件下的處理結果進行對比,從而得到最優的降噪方法。
仿真驗證時如果采用脈沖激光雷達回波實測數據,無法對原始數據和濾波后數據的信號進行評價,因而必須構建仿真理想信號。嚴格地說,每種大氣條件下的信號模型都是不同的。它受到大氣粒子尺寸分布,粒子形狀分布、濃度,粒子性質及激光雷達發射和接收視場角等多種因素影響。研究者通常把大氣分為3種情況進行建模,把非常稀薄和非常稠密的粒子群分布單獨考慮,分別采用單粒子散射近似和漫反射近似,而最普遍的情況是采用一級多次散射近似進行建模。大氣脈沖激光雷達方程(式(1))正是在一級多次散射近似條件下推導出來的,適用于絕大多數大氣回波情況,也與實測信號最為接近。

式中:Pr(R)為距離R處的回波功率;βπ(R)為后向散射系數;σ(R)為消光系數;Cq為儀器常數,它又可以表示為:

式中:Ar為接收鏡頭面積;Tt和Tr分別為發射/接收光學系統透過率;E0為發射脈沖能量;τp為脈沖寬度;c為光速;η(R)為幾何重疊因子。
后向散射系數是激光雷達方程中的重要參數,根據Klett的研究[3],在米散射起主導作用時,消光系數與后向散射系數的關系為β(z)=B0σk(z),B0為后向散射比,k的取值和波長與粒子尺度的相對值有關,通常取值區間為[0.66,1]。
實際的激光雷達接收信號要通過光電探測器、高頻放大器、模數轉換器等器件。根據實際器件參數,將信號模型進一步修正。設放大器的偏置電平為lv,AD轉換器的采樣頻率為f,同時假定第n-1個采樣點與第n個采樣點范圍內采樣體積內的平均消光系數為,回波模型可以修正為:


圖1 仿真、實測回波信號比較Fig.1Compare between real and simulated signal
為了檢驗模型的正確性,將真實信號與同等條件仿真波形進行比較。圖1中虛線為某型激光雷達在某機場晴朗均勻大氣中的實測回波信號,實線是采用式(3)擬合的理想回波信號。可以看出,仿真信號對遠場回波信號擬合較好,近場信號存在差別,這是幾何重疊因子的差異造成的。儀器常數的諸多因子中,幾何重疊因子是難以確定的。有研究人員假設光束能量高斯分布通過激光光束發射與接收視場截面的大小來確定幾何重疊因子的大小[4]。但激光光束能量的分布特性的隨機性大,近場能量強、多次散射強烈使得這種方法意義不大。因此,研究中只關注視場完全重疊后的區域,近場擬合的精確性對仿真效果并無影響。
激光雷達的噪聲按形成原因可分為散粒噪聲、熱噪聲、背景光噪聲及1/f噪聲等[2]。其中除了1/f噪聲外,其他噪聲都可認為是加性噪聲,并且符合高斯分布。因此研究中采用高斯噪聲模型。
評價信號質量通常采用的參數是信噪比,但信號信噪比的定義有很多,本文選用信號幅值均方根與噪聲幅度均方根之比來定義。

式中:Pi為信號幅度;Ni為噪聲幅度;n為信號長度。
在激光雷達領域,研究者將擴展卡爾曼算法當作1種大氣狀態參數的反演方法[5-6]。這樣難以對算法降噪性能加以評估,也無法驗證反演模型的準確性。因此筆者傾向于把EKF法修正為1種降噪算法。
將消光系數作為信號主要參量,信號模型為

其中:v(n)為消光系數波動參數,它服從正態分布,方差為R,均值為0。
偏置電平很容易濾除,去掉放大器偏置電平并考慮加性噪聲,式(3)可化簡為式(6),其中w(n)均值為0,方差為Q。

初始條件為P(n0)=0,n0為發射與接收視場的完全重疊位置,此時η(n0)=1。
小波濾波采用多尺度正交小波基對信號進行分解,將信號和噪聲在不同時頻域進行分離,然后在各時頻范圍內對依據閾值進行濾波,最后重構信號從而達到降噪的目的。離散小波基的定義為:

任何函數都可以分解為小波基的加權和,

其中,加權值即為小波系數:

小波理論產生于20世紀80年代,學者們對小波降噪方法的研究廣泛且深入,因此本文不再討論。對于激光雷達信號降噪問題,小波基和閾值的選取是最重要的2個方面。不同小波基的特性不同,選取不同的小波基對濾波結果影響很大。Db5小波具有雙正交性、緊支撐性和近似對稱性,十分適合對指數衰減型曲線進行濾波。閾值的選取對濾波效果影響也很大,常用的降噪閾值選取方法有長度對數閾值(Sqtwolog)、最小極大方差閾值(Minimaxi)、無偏風險閾值(Rigrsure)、啟發式閾值(Heursure)等[8-10]。這些方法的主要思想是根據小波基長度與噪聲均方根構造最佳閾值,它們都可以改善信號信噪比,但降噪能力有限。對于本課題來說,仿真表明固定閾值降噪效果更好,當閾值選擇為二倍噪聲均方根時可使濾波效果達到最佳。
EMD分解濾波也是一種基于時頻分析的算法,它與小波分解算法的區別在于對瞬時頻率定義的差別。EMD算法具有自適應性,它會為每條信號曲線產生不同的正交基。
經驗模態分解(EMD),是通過固有模態函數(IMF)提取來實現的[11]。1個函數可以分解為n個固有模態函數ci(t)與1個剩余分量r(t)的和,即

固有模態函數IMF必須滿足以下2個條件:一是信號持續時間內,其過0點的個數和極值點的個數相等或最多相差1個;二是信號極大值或極小值構成的包絡線的均值為0。
固有模態函數采用以下步驟進行提取:尋找信號f(t)的局部極值點;在極大值之間進行插值擬合處理形成最大值包絡曲線maxi(t),在極小值之間進行插值擬合處理形成最小值包絡曲線mini(t);求包絡均值曲線m(t)=[max(t)+min(t)]/2;信號減去包絡均值曲線;多次迭代后,剩余信號滿足停止準則時篩選結束。篩選門限見式(15),篩選門限一般在0.2~0.3之間。

經過EMD提取之后,信號被分解為多個固有模態函數,對固有函數進行閾值濾波,將濾波后的數據重構就可以達到降噪的目的。
仿真表明,篩選門限取值2.5時濾波效果較好,軟閾值法的效果比硬閾值法濾波效果好,當閾值為噪聲均方根時可使降噪效果達到最佳。
仿真主要針對2種大氣條件下的回波信號:均勻大氣和局部突變大氣。均勻大氣主要針對無風或微風條件下的自然大氣狀態,這時消光系數比較穩定隨距離起伏不大。局部突變大氣針對存在風切變、湍流及云霧或煙塵等天氣條件,它們會改變大氣的局部分布,造成局部消光系數突變。
常規大氣條件下,消光系數均勻,在均值附近上下起伏。因而選擇不同的消光系數均值和方差構造激光脈沖回波,并疊加高斯噪聲。如圖2和圖3所示。

仿真依據消光系數從低至高進行,每個消光系數取值又采用不同的噪聲幅度進行仿真,仿真結果如表2所示。

表2中第1行為仿真消光系數值,第1列為噪聲均方根(δ)。表中數據依次為原始信號EKF濾波小波濾波EMD濾波對應的信噪比。
此種條件時,絕大部分消光系數分布與均勻大氣下相同,只是局部空氣中消光系數陡增驟減。當消光系數如圖4所示分布時,回波信號波形如圖5所示,仿真結果如表3所示。

表3中第1行為仿真消光系數值,第1列為噪聲均方根(δ)。表中數據依次為原始信號EKF濾波小波濾波EMD濾波對應的信噪比。
擴展卡爾曼濾波的降噪能力有限,但可以將信號直接反演為消光系數,初始值的選取對這種方法的濾波效果影響很大。消光系數較高時,EKF濾波后的失真較大;在消光系數較低時的濾波效果較好。在均勻大氣條件下,EKF法反演精度較高;而當大氣中有局部突變時,若消光系數陡增幅度過大會破壞這種方法的成立條件,導致誤差較大。小波濾波的效果在多數情況下是幾種方法中的最佳方法,信噪比較高時,小波濾波的優勢并不明顯,并且有可能造成信號失真;中低信噪比條件下,這種方法對信號質量的改善最為明顯;在突變大氣條件下,盡管這種方法的濾波效果降低了,但仍然是最佳方法。在信噪比較高時EMD濾波效果明顯好于其他2種方法。目前采用EMD法對激光雷達回波信號進行濾波的研究仍處于起步階段,在降噪閾值和篩選門限這2個方面進一步研究會使該方法的優勢得到進一步體現。
綜上所述,本文根據激光雷達性能參數和典型的大氣條件建立了1種更為實用的回波信號仿真模型。本文工作證明EMD濾波法在某些大氣條件下可以達到最佳濾波效果,是現有方法的必要補充。綜合采用3種方法才能使濾波達到最佳效果。均勻大氣條件下,如果噪聲幅度較高,應采用小波濾波;噪聲幅度較低時,如果消光系數較小,應采用EKF濾波;如果消光系數較大,就應采用EMD算法。突變大氣條件下,如果噪聲幅度較高,宜采用小波濾波;如果噪聲幅度較低則應采用EMD濾波。將該方法應用到激光雷達消光系數測量工程中,得到了較好的效果。
[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]王治華,賀應紅,左浩毅,等.基于高斯光束特性的Mie散射大氣激光雷達回波近場信號校正研究[J].物理學報,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].工程熱物理學報,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.