基于經(jīng)驗模式分解的大氣邊界層高度提取

張娜　姜茂仁　李進杰

【摘 要】研究大氣邊界層的實時高度是研究天氣、氣候和大氣環(huán)境污染的迫切需要。激光雷達成為探測邊界層時空演變特征的最有效手段。目前國際上用于反演激光雷達數(shù)據(jù)獲取大氣邊界層高度的主要方法有梯度法、標準偏差法、曲線擬合法等。本文采用偏振拉曼—米散射激光雷達的實測數(shù)據(jù)，介紹梯度法中的二階求導反演激光雷達數(shù)據(jù)提取大氣邊界層高度，并采用經(jīng)驗模式分解（EMD）方法處理數(shù)據(jù)。通過對結(jié)果進行分析，發(fā)現(xiàn)提取的大氣邊界層高度值與利用診斷式求得的高度值的吻合性比較理想。

【關(guān)鍵詞】大氣邊界層；激光雷達；EMD

The Applidation of Emperical Mode Decomposition Method in Calculation Atmospheric Boundary Layer Height

ZHANG Na JIANG Mao-ren LI Jin-jie

（Naval aeronautical engineering academy qingdao branch， Qingdao Shandong 266041， China）

【Abstract】Research for the real-time atmospheric boundary layer height is important to the weather， climate and atmospheric pollution. The determination of atmospheric boundary layer height based on lidar data is the most effective tool. Currently in the world， the main methods to determine atmospheric boundary layer height from lidar data are three kinds of commonly used methods， such as the gradient method， standard deviation method， fitting mehord. This paper expounded the gradient method calculating atmospheric boundary layer height from lidar experiment data， and attempted the Empirical Mode Decomposition method which was first made use of processing data， with the continuously measured data of polarization-Raman-Mie scattering lidar. The results indicate that the height from EMD accord with that from diagnostic equation.

【Key words】Atmospheric boundary layer； Lidar； Empirical Mode Decomposition

0 引言

大氣邊界層（Atmospheric Boundary Layer， ABL）是指靠近地球表面，受地面摩擦阻力等表面強迫力影響并且響應時間尺度約為1小時或更小的大氣層[1]。它的高度隨天氣過程、地形和地表特征而變化，一般在2公里范圍內(nèi)。大氣邊界層是與人類關(guān)系密切的大氣層，是大氣污染物的主要集中層，邊界層高度的大小，決定了排入大氣中污染物的濃度值。所以針對氣候及大氣環(huán)境的研究需要，開展對大氣邊界層高度的觀測勢在必行。同時，對邊界層高度的觀測要求也越來越高。

激光雷達作為一種新型大氣主動遙感探測工具[2-4]，已經(jīng)廣泛應用于大氣氣溶膠和污染物、云等領(lǐng)域的測量和研究工作。激光雷達利用氣溶膠作為示蹤粒子，所測量到的光功率信號正比于大氣的氣溶膠含量。與此相對應，大氣邊界層中含有豐富的氣溶膠（較自由對流層的氣溶膠含量更大），這是因為氣溶膠粒子主要來源于地球表面，從而導致對激光的更大散射，因此可以利用激光雷達較容易的探測到這兩層之間的邊界。

本文利用激光雷達探測系統(tǒng)餓實測數(shù)據(jù)，采用了梯度法中的二階求導法來提取大氣邊界層高度，并在分析處理中應用了經(jīng)驗模式分解（EMD）方法，給出了典型天氣情況下的提取結(jié)果，并與診斷式求得的高度值進行比對，吻合度較高。

1 梯度法反演大氣邊界層高度

因為邊界層頂與自由對流層氣溶膠濃度差距較大，所以回波信號的梯度最大值點對應的高度邊界層厚度[5-7]。常用的梯度計算方法為一階導數(shù)法或二階導數(shù)法。計算激光雷達回波距離校正信號關(guān)于高度的二階導數(shù)，并計算其絕對最小值，最小值對應的高度定義為過渡層，即混合層和自由對流層交界的中間位置。一階導數(shù)的最小值對應高度是混合層的頂部。

激光雷達接收到的米散射回波信號可以式（1）所示的激光雷達方程表示：

P（z）=P■CY（z）z■[β■（z）+β■（z）]T■■（z）T■■（z）（1）

式中：P（z）是激光雷達接收到距離z處氣溶膠粒子和空氣分子的后向散射回波信號；P■是激光發(fā)射功率；C是激光雷達系統(tǒng)常數(shù)；Y（z）是幾何因子；β■（z）和β■（z）分別是大氣氣溶膠粒子和大氣分子后向散射系數(shù)，T■■（z）=exp-■？鄣■（z′）dz′和T■■（z）=exp-■？鄣■（z′）dz′分別是大氣氣溶膠粒子與大氣分子透過率；？鄣■（z）和？鄣■（z）分別為大氣氣溶膠粒子和大氣分子消光系數(shù)。

由（1）式可得：

P（z）z■=P■CY（z）[β■（z）+β■（z）]T■■（z）T■■（z）（2）

激光雷達距離平方校正回波信號P（z）z■在一定程度上反映了大氣氣溶膠濃度隨探測高度變化的情況。由于覆蓋逆溫的作用，大量的大氣氣溶膠粒子富集在大氣邊界層以內(nèi)，這樣大氣邊界層大自由大氣層之間的大氣氣溶膠濃度就會發(fā)生變化。P（z）z■的梯度D（z）的變化即代表著大氣氣溶膠濃度梯度的變化[8]。

D（z）定義為D（z）=d[P（z）z■]/dz該梯度變化的最大位置就是大氣邊界層的高度。

2 EMD方法在信號處理中的應用

經(jīng)驗模式分解（EMD）方法的本質(zhì)是對信號進行平穩(wěn)化處理，根據(jù)特征時間尺度把一個復雜的數(shù)據(jù)序列分解成有限項內(nèi)在模函數(shù)（Instrinsic Mode Functions， IMF）[9-10]。最低頻率的IMF分量通常情況下代表原始數(shù)據(jù)的整體趨勢或均值，確定IMF的標準為：（1）整個數(shù)據(jù)集中，極值點的數(shù)目和跨零點的數(shù)目必須相等或至多只差一個；（2）任一點上，由局部極值定義的上下包絡線的均值必須等于零。

作為一種應用，EMD 分解方法可以有效地提取一個數(shù)據(jù)序列的趨勢或去掉該數(shù)據(jù)序列的均值。IMFS是通過篩選過程從數(shù)據(jù)中提取出的：找出原序列x（t）的各個局部極大值、極小值，用三階樣條函數(shù)分別進行插值，相應得到原序列x（t）的上包絡序列值xmax（t）及下包絡序列值xmin（t）；上下包絡線的均值被指定為m10；原序列x（t）和均值m10的差值h10是第一個組成。

x（t）-m10=h10

理想情況下，h10應該是一個IMF，但由于非線性數(shù)據(jù)包絡線的平均值可能不同于實際的局部均值，導致非對稱波的存在，因此重復這一篩選過程就非常必要了。把h10當作原序列，重復以上處理過程k次，直至滿足內(nèi)在模函數(shù)的定義，求出第一個內(nèi)在模函數(shù)為止。

h1k=h1（k-1）-m1k；

c1=h1k；

第一個IMF分量c1應包含原始信號中最短的周期分量。第二個IMF函數(shù)c2通過對余數(shù)項r1=x（t）-c1進行篩選得到。通常第n個IMF分量通過篩選rn=rn-1-cn得到。篩選過程直到下列任一標準滿足：（1）當分量cn或剩余分量rn比預定值小；（2）當剩余分量rn變成單調(diào)函數(shù)時，從中不能再篩選出基本模式分量。對于有趨勢的數(shù)據(jù)，其剩余分量就應該是該趨勢。如果把分解后的各分量合并起來，就得到原序列x（t）=■ci+rn。

EMD分解方法的特點是在分解的過程中，內(nèi)在模函數(shù)保留了原始數(shù)據(jù)本身的特性，各項內(nèi)在模函數(shù)的和等于原始數(shù)據(jù)序列，這就允許可以只用IMFS來分析和處理數(shù)據(jù)，保留了數(shù)據(jù)本身的特性。

3 EMD方法在信號處理中的應用

本次實驗中，激光雷達以脈沖能量50mJ出射激光，重復頻率1Hz，所接收到的大氣回波信號經(jīng)6分鐘時間的平均后存儲一次，距離分辨率為15m。

圖1為2014年11月03日晚20：00的大氣回波信號（未校準），考慮到大氣邊界層的高度通常為2km以下，通常我們只需處理從地面到2km這一段高度的數(shù)據(jù)。但本文為了使比較結(jié)果更形象明顯，故對采樣數(shù)據(jù)全部進行處理。

圖1 激光雷達回波信號（2014/11/03 20：00）

對該段數(shù)據(jù)先進行預處理，接著用上述介紹過的EMD方法分解，按照頻率由高到低得到16個IMFS，見圖2。從各個內(nèi)在模函數(shù)可以看出，其分解出的前6項IMF代表原始數(shù)據(jù)中較高頻的組成，對信號的能量貢獻最少，我們可以考慮將其消除，只處理低頻的部分。圖3為消除高頻項后的信號與原始信號的比較。通過比較可發(fā)現(xiàn)，二者基本吻合，消除高頻項后的信號較原始信號平滑的多，且仍能顯著地描述出原始信號的變化趨勢。將經(jīng)過以上處理的信號采用Menut求取二階導數(shù)的方法，認為二階導數(shù)最小值所對應的高度為大氣邊界層高度。從圖4可以看出，由EMD預先處理再求二階導數(shù)，信號強度最小值對應的高度為1290m。

圖2 距離校準信號中分解出的IMFS

在此本文把通過以上方法求得的邊界層高度值與根據(jù)常規(guī)氣象資料計算出的值比較，邊界層的計算公式采用Benkley and Schulman提供的中性和穩(wěn)定情況下的關(guān)系式：h=a×■，式中u？鄢為摩擦速度（m/s）；f為地轉(zhuǎn)參數(shù)（1/s）；a為常數(shù)，由氣象資料決定，不同研究者給出的a值不同。本文中，我們根據(jù)u？鄢值的范圍選擇a，文中有兩種情況，a=0.31或a=0.56。利用青島氣象臺提供的氣象資料及以上關(guān)系式，我們得到11月03日18：00時的大氣邊界層厚度值為1053m。

下面，我們把2014年10月和11月測量的部分結(jié)果在表中作比較，同時表中亦列出僅根據(jù)Menut的距離校準信號的二階微分最小確定的大氣邊界層厚度。

表4 兩種方法提取邊界層厚度的結(jié)果比較

從以上實驗結(jié)果可計算出：激光雷達回波信號用EMD處理結(jié)合二階導數(shù)最小值對應的大氣邊界層高度與氣象資料所確定的高度值的相關(guān)性高達76.32%，依據(jù)二階微分最小確定的大氣邊界層高度與氣象資料所確定的高度值的相關(guān)性為52.11%。顯然，前者要絕對優(yōu)于后者。也就是說，用EMD結(jié)合二階微分最小方法確定大氣邊界層高度相比較與僅利用二階微分最小的方法來確定大氣邊界層高度，結(jié)果得到明顯改善，進一步證明了經(jīng)驗模式分解方法在大氣邊界層厚度的激光雷達提取應用中的可行性。同時，由于邊界層高度的診斷關(guān)系式受探測時大氣邊界層穩(wěn)定度的影響，所以與EMD結(jié)合二階微分最小確定邊界層高度的方法相比，要相對復雜得多。

EMD結(jié)合二階微分最小確定邊界層高度的方法，僅根據(jù)激光雷達所接收到的回波信號（與大氣中氣溶膠含量成正比），利用經(jīng)驗模式分解與Menut的二階微分求導二者的結(jié)合，即可得到大氣邊界層高度值，方法簡單而有效。

【參考文獻】

[1]Stull R B. Introduction to Boundary Layer Meteorology[M]. Norwell： Kluwer Academic Publisher， 1988

[2]Menut， L.， Flamant， C.， Pelon， J.， Flamant， P.H.， Urban boundary-layer height determination from lidar measurements over the Paris area， Applied Optics， Vol.38，No.6， 945-954，1999[Z].

[3]Steyn D G， Bottenheim J W， Thomson R B. Over-view of tropospheric ozone in the lower fraser valley and the Pacific 93 field study[J]. J. Atmos. Environ.， 1997， 31：2025-2035.

[4]Davis K J， Len Schow D H，Oncley S P， et al. Role of entrainment in surface-atmosphere interaction over the boreal forest[J]. J.Geophy Res， 1997， 102（D24）： 29219-29230.

[5]Arya S P S. Parameterizing the height of the stable atmospheric boundary layer. JAM.，20，1981： 1192-1202[Z].

[6]Estournel C，Guedalia D. Improving the diagnostic relation for the nocturnal boundary layer height. BLM， 53， 1990： 191-198[Z].

[7]Surridge A D，Swanepoel D J.On the evolution of the height and temperature difference across the nocturnal stable Bondary Layer. BLM， 40， 1987： 87[Z].

[8]王珍珠，李炬，鐘志慶，等.激光雷達探測北京城區(qū)夏季大氣邊界層[J].應用光學，2008，29（1）：96-100.

[9]Norden E.Huang The empirical mode decomposition and Hilbert spectrum for nonlinear and non-stationary time series analysis. Proc. R. Soc. Lond A 1998， 454（1971）： 903-995[Z].

[10]Norden E. Huang， et all. A confidence limit for the empirical mode decomposition and Hilbert spectral analysis， Proc. Roy. Soc. London A， Vol. 459， （2003）： 2317-2345[Z].

[責任編輯：楊玉潔]