基于經(jīng)驗(yàn)?zāi)Ｊ椒纸獾拇髿膺吔鐚痈叨忍崛?/h1>
2016-06-13 18:34:34張娜姜茂仁李進(jìn)杰
科技視界 2016年13期
張娜 姜茂仁 李進(jìn)杰
【摘 要】研究大氣邊界層的實(shí)時(shí)高度是研究天氣、氣候和大氣環(huán)境污染的迫切需要。激光雷達(dá)成為探測(cè)邊界層時(shí)空演變特征的最有效手段。目前國(guó)際上用于反演激光雷達(dá)數(shù)據(jù)獲取大氣邊界層高度的主要方法有梯度法、標(biāo)準(zhǔn)偏差法、曲線擬合法等。本文采用偏振拉曼—米散射激光雷達(dá)的實(shí)測(cè)數(shù)據(jù)，介紹梯度法中的二階求導(dǎo)反演激光雷達(dá)數(shù)據(jù)提取大氣邊界層高度，并采用經(jīng)驗(yàn)?zāi)Ｊ椒纸猓‥MD）方法處理數(shù)據(jù)。通過(guò)對(duì)結(jié)果進(jìn)行分析，發(fā)現(xiàn)提取的大氣邊界層高度值與利用診斷式求得的高度值的吻合性比較理想。
【關(guān)鍵詞】大氣邊界層；激光雷達(dá)；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）是指靠近地球表面，受地面摩擦阻力等表面強(qiáng)迫力影響并且響應(yīng)時(shí)間尺度約為1小時(shí)或更小的大氣層[1]。它的高度隨天氣過(guò)程、地形和地表特征而變化，一般在2公里范圍內(nèi)。大氣邊界層是與人類關(guān)系密切的大氣層，是大氣污染物的主要集中層，邊界層高度的大小，決定了排入大氣中污染物的濃度值。所以針對(duì)氣候及大氣環(huán)境的研究需要，開(kāi)展對(duì)大氣邊界層高度的觀測(cè)勢(shì)在必行。同時(shí)，對(duì)邊界層高度的觀測(cè)要求也越來(lái)越高。
激光雷達(dá)作為一種新型大氣主動(dòng)遙感探測(cè)工具[2-4]，已經(jīng)廣泛應(yīng)用于大氣氣溶膠和污染物、云等領(lǐng)域的測(cè)量和研究工作。激光雷達(dá)利用氣溶膠作為示蹤粒子，所測(cè)量到的光功率信號(hào)正比于大氣的氣溶膠含量。與此相對(duì)應(yīng)，大氣邊界層中含有豐富的氣溶膠（較自由對(duì)流層的氣溶膠含量更大），這是因?yàn)闅馊苣z粒子主要來(lái)源于地球表面，從而導(dǎo)致對(duì)激光的更大散射，因此可以利用激光雷達(dá)較容易的探測(cè)到這兩層之間的邊界。
本文利用激光雷達(dá)探測(cè)系統(tǒng)餓實(shí)測(cè)數(shù)據(jù)，采用了梯度法中的二階求導(dǎo)法來(lái)提取大氣邊界層高度，并在分析處理中應(yīng)用了經(jīng)驗(yàn)?zāi)Ｊ椒纸猓‥MD）方法，給出了典型天氣情況下的提取結(jié)果，并與診斷式求得的高度值進(jìn)行比對(duì)，吻合度較高。
1 梯度法反演大氣邊界層高度
因?yàn)檫吔鐚禹斉c自由對(duì)流層氣溶膠濃度差距較大，所以回波信號(hào)的梯度最大值點(diǎn)對(duì)應(yīng)的高度邊界層厚度[5-7]。常用的梯度計(jì)算方法為一階導(dǎo)數(shù)法或二階導(dǎo)數(shù)法。計(jì)算激光雷達(dá)回波距離校正信號(hào)關(guān)于高度的二階導(dǎo)數(shù)，并計(jì)算其絕對(duì)最小值，最小值對(duì)應(yīng)的高度定義為過(guò)渡層，即混合層和自由對(duì)流層交界的中間位置。一階導(dǎo)數(shù)的最小值對(duì)應(yīng)高度是混合層的頂部。
激光雷達(dá)接收到的米散射回波信號(hào)可以式（1）所示的激光雷達(dá)方程表示：
P（z）=P■CY（z）z■[β■（z）+β■（z）]T■■（z）T■■（z）（1）
式中：P（z）是激光雷達(dá)接收到距離z處氣溶膠粒子和空氣分子的后向散射回波信號(hào)；P■是激光發(fā)射功率；C是激光雷達(dá)系統(tǒng)常數(shù)；Y（z）是幾何因子；β■（z）和β■（z）分別是大氣氣溶膠粒子和大氣分子后向散射系數(shù)，T■■（z）=exp-■？鄣■（z′）dz′和T■■（z）=exp-■？鄣■（z′）dz′分別是大氣氣溶膠粒子與大氣分子透過(guò)率；？鄣■（z）和？鄣■（z）分別為大氣氣溶膠粒子和大氣分子消光系數(shù)。
由（1）式可得：
P（z）z■=P■CY（z）[β■（z）+β■（z）]T■■（z）T■■（z）（2）
激光雷達(dá)距離平方校正回波信號(hào)P（z）z■在一定程度上反映了大氣氣溶膠濃度隨探測(cè)高度變化的情況。由于覆蓋逆溫的作用，大量的大氣氣溶膠粒子富集在大氣邊界層以內(nèi)，這樣大氣邊界層大自由大氣層之間的大氣氣溶膠濃度就會(huì)發(fā)生變化。P（z）z■的梯度D（z）的變化即代表著大氣氣溶膠濃度梯度的變化[8]。
D（z）定義為D（z）=d[P（z）z■]/dz該梯度變化的最大位置就是大氣邊界層的高度。
2 EMD方法在信號(hào)處理中的應(yīng)用
經(jīng)驗(yàn)?zāi)Ｊ椒纸猓‥MD）方法的本質(zhì)是對(duì)信號(hào)進(jìn)行平穩(wěn)化處理，根據(jù)特征時(shí)間尺度把一個(gè)復(fù)雜的數(shù)據(jù)序列分解成有限項(xiàng)內(nèi)在模函數(shù)（Instrinsic Mode Functions， IMF）[9-10]。最低頻率的IMF分量通常情況下代表原始數(shù)據(jù)的整體趨勢(shì)或均值，確定IMF的標(biāo)準(zhǔn)為：（1）整個(gè)數(shù)據(jù)集中，極值點(diǎn)的數(shù)目和跨零點(diǎn)的數(shù)目必須相等或至多只差一個(gè)；（2）任一點(diǎn)上，由局部極值定義的上下包絡(luò)線的均值必須等于零。
作為一種應(yīng)用，EMD 分解方法可以有效地提取一個(gè)數(shù)據(jù)序列的趨勢(shì)或去掉該數(shù)據(jù)序列的均值。IMFS是通過(guò)篩選過(guò)程從數(shù)據(jù)中提取出的：找出原序列x（t）的各個(gè)局部極大值、極小值，用三階樣條函數(shù)分別進(jìn)行插值，相應(yīng)得到原序列x（t）的上包絡(luò)序列值xmax（t）及下包絡(luò)序列值xmin（t）；上下包絡(luò)線的均值被指定為m10；原序列x（t）和均值m10的差值h10是第一個(gè)組成。
x（t）-m10=h10
理想情況下，h10應(yīng)該是一個(gè)IMF，但由于非線性數(shù)據(jù)包絡(luò)線的平均值可能不同于實(shí)際的局部均值，導(dǎo)致非對(duì)稱波的存在，因此重復(fù)這一篩選過(guò)程就非常必要了。把h10當(dāng)作原序列，重復(fù)以上處理過(guò)程k次，直至滿足內(nèi)在模函數(shù)的定義，求出第一個(gè)內(nèi)在模函數(shù)為止。
h1k=h1（k-1）-m1k；
c1=h1k；
第一個(gè)IMF分量c1應(yīng)包含原始信號(hào)中最短的周期分量。第二個(gè)IMF函數(shù)c2通過(guò)對(duì)余數(shù)項(xiàng)r1=x（t）-c1進(jìn)行篩選得到。通常第n個(gè)IMF分量通過(guò)篩選rn=rn-1-cn得到。篩選過(guò)程直到下列任一標(biāo)準(zhǔn)滿足：（1）當(dāng)分量cn或剩余分量rn比預(yù)定值小；（2）當(dāng)剩余分量rn變成單調(diào)函數(shù)時(shí)，從中不能再篩選出基本模式分量。對(duì)于有趨勢(shì)的數(shù)據(jù)，其剩余分量就應(yīng)該是該趨勢(shì)。如果把分解后的各分量合并起來(lái)，就得到原序列x（t）=■ci+rn。
EMD分解方法的特點(diǎn)是在分解的過(guò)程中，內(nèi)在模函數(shù)保留了原始數(shù)據(jù)本身的特性，各項(xiàng)內(nèi)在模函數(shù)的和等于原始數(shù)據(jù)序列，這就允許可以只用IMFS來(lái)分析和處理數(shù)據(jù)，保留了數(shù)據(jù)本身的特性。
3 EMD方法在信號(hào)處理中的應(yīng)用
本次實(shí)驗(yàn)中，激光雷達(dá)以脈沖能量50mJ出射激光，重復(fù)頻率1Hz，所接收到的大氣回波信號(hào)經(jīng)6分鐘時(shí)間的平均后存儲(chǔ)一次，距離分辨率為15m。
圖1為2014年11月03日晚20：00的大氣回波信號(hào)（未校準(zhǔn)），考慮到大氣邊界層的高度通常為2km以下，通常我們只需處理從地面到2km這一段高度的數(shù)據(jù)。但本文為了使比較結(jié)果更形象明顯，故對(duì)采樣數(shù)據(jù)全部進(jìn)行處理。
圖1 激光雷達(dá)回波信號(hào)（2014/11/03 20：00）
對(duì)該段數(shù)據(jù)先進(jìn)行預(yù)處理，接著用上述介紹過(guò)的EMD方法分解，按照頻率由高到低得到16個(gè)IMFS，見(jiàn)圖2。從各個(gè)內(nèi)在模函數(shù)可以看出，其分解出的前6項(xiàng)IMF代表原始數(shù)據(jù)中較高頻的組成，對(duì)信號(hào)的能量貢獻(xiàn)最少，我們可以考慮將其消除，只處理低頻的部分。圖3為消除高頻項(xiàng)后的信號(hào)與原始信號(hào)的比較。通過(guò)比較可發(fā)現(xiàn)，二者基本吻合，消除高頻項(xiàng)后的信號(hào)較原始信號(hào)平滑的多，且仍能顯著地描述出原始信號(hào)的變化趨勢(shì)。將經(jīng)過(guò)以上處理的信號(hào)采用Menut求取二階導(dǎo)數(shù)的方法，認(rèn)為二階導(dǎo)數(shù)最小值所對(duì)應(yīng)的高度為大氣邊界層高度。從圖4可以看出，由EMD預(yù)先處理再求二階導(dǎo)數(shù)，信號(hào)強(qiáng)度最小值對(duì)應(yīng)的高度為1290m。
圖2 距離校準(zhǔn)信號(hào)中分解出的IMFS
在此本文把通過(guò)以上方法求得的邊界層高度值與根據(jù)常規(guī)氣象資料計(jì)算出的值比較，邊界層的計(jì)算公式采用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。利用青島氣象臺(tái)提供的氣象資料及以上關(guān)系式，我們得到11月03日18：00時(shí)的大氣邊界層厚度值為1053m。
下面，我們把2014年10月和11月測(cè)量的部分結(jié)果在表中作比較，同時(shí)表中亦列出僅根據(jù)Menut的距離校準(zhǔn)信號(hào)的二階微分最小確定的大氣邊界層厚度。
表4 兩種方法提取邊界層厚度的結(jié)果比較
從以上實(shí)驗(yàn)結(jié)果可計(jì)算出：激光雷達(dá)回波信號(hào)用EMD處理結(jié)合二階導(dǎo)數(shù)最小值對(duì)應(yīng)的大氣邊界層高度與氣象資料所確定的高度值的相關(guān)性高達(dá)76.32%，依據(jù)二階微分最小確定的大氣邊界層高度與氣象資料所確定的高度值的相關(guān)性為52.11%。顯然，前者要絕對(duì)優(yōu)于后者。也就是說(shuō)，用EMD結(jié)合二階微分最小方法確定大氣邊界層高度相比較與僅利用二階微分最小的方法來(lái)確定大氣邊界層高度，結(jié)果得到明顯改善，進(jìn)一步證明了經(jīng)驗(yàn)?zāi)Ｊ椒纸夥椒ㄔ诖髿膺吔鐚雍穸鹊募す饫走_(dá)提取應(yīng)用中的可行性。同時(shí)，由于邊界層高度的診斷關(guān)系式受探測(cè)時(shí)大氣邊界層穩(wěn)定度的影響，所以與EMD結(jié)合二階微分最小確定邊界層高度的方法相比，要相對(duì)復(fù)雜得多。
EMD結(jié)合二階微分最小確定邊界層高度的方法，僅根據(jù)激光雷達(dá)所接收到的回波信號(hào)（與大氣中氣溶膠含量成正比），利用經(jīng)驗(yàn)?zāi)Ｊ椒纸馀cMenut的二階微分求導(dǎo)二者的結(jié)合，即可得到大氣邊界層高度值，方法簡(jiǎn)單而有效。
【參考文獻(xiàn)】
[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]王珍珠，李炬，鐘志慶，等.激光雷達(dá)探測(cè)北京城區(qū)夏季大氣邊界層[J].應(yīng)用光學(xué)，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].
[責(zé)任編輯：楊玉潔]
張娜 姜茂仁 李進(jìn)杰
【摘 要】研究大氣邊界層的實(shí)時(shí)高度是研究天氣、氣候和大氣環(huán)境污染的迫切需要。激光雷達(dá)成為探測(cè)邊界層時(shí)空演變特征的最有效手段。目前國(guó)際上用于反演激光雷達(dá)數(shù)據(jù)獲取大氣邊界層高度的主要方法有梯度法、標(biāo)準(zhǔn)偏差法、曲線擬合法等。本文采用偏振拉曼—米散射激光雷達(dá)的實(shí)測(cè)數(shù)據(jù)，介紹梯度法中的二階求導(dǎo)反演激光雷達(dá)數(shù)據(jù)提取大氣邊界層高度，并采用經(jīng)驗(yàn)?zāi)Ｊ椒纸猓‥MD）方法處理數(shù)據(jù)。通過(guò)對(duì)結(jié)果進(jìn)行分析，發(fā)現(xiàn)提取的大氣邊界層高度值與利用診斷式求得的高度值的吻合性比較理想。
【關(guān)鍵詞】大氣邊界層；激光雷達(dá)；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）是指靠近地球表面，受地面摩擦阻力等表面強(qiáng)迫力影響并且響應(yīng)時(shí)間尺度約為1小時(shí)或更小的大氣層[1]。它的高度隨天氣過(guò)程、地形和地表特征而變化，一般在2公里范圍內(nèi)。大氣邊界層是與人類關(guān)系密切的大氣層，是大氣污染物的主要集中層，邊界層高度的大小，決定了排入大氣中污染物的濃度值。所以針對(duì)氣候及大氣環(huán)境的研究需要，開(kāi)展對(duì)大氣邊界層高度的觀測(cè)勢(shì)在必行。同時(shí)，對(duì)邊界層高度的觀測(cè)要求也越來(lái)越高。
激光雷達(dá)作為一種新型大氣主動(dòng)遙感探測(cè)工具[2-4]，已經(jīng)廣泛應(yīng)用于大氣氣溶膠和污染物、云等領(lǐng)域的測(cè)量和研究工作。激光雷達(dá)利用氣溶膠作為示蹤粒子，所測(cè)量到的光功率信號(hào)正比于大氣的氣溶膠含量。與此相對(duì)應(yīng)，大氣邊界層中含有豐富的氣溶膠（較自由對(duì)流層的氣溶膠含量更大），這是因?yàn)闅馊苣z粒子主要來(lái)源于地球表面，從而導(dǎo)致對(duì)激光的更大散射，因此可以利用激光雷達(dá)較容易的探測(cè)到這兩層之間的邊界。
本文利用激光雷達(dá)探測(cè)系統(tǒng)餓實(shí)測(cè)數(shù)據(jù)，采用了梯度法中的二階求導(dǎo)法來(lái)提取大氣邊界層高度，并在分析處理中應(yīng)用了經(jīng)驗(yàn)?zāi)Ｊ椒纸猓‥MD）方法，給出了典型天氣情況下的提取結(jié)果，并與診斷式求得的高度值進(jìn)行比對(duì)，吻合度較高。
1 梯度法反演大氣邊界層高度
因?yàn)檫吔鐚禹斉c自由對(duì)流層氣溶膠濃度差距較大，所以回波信號(hào)的梯度最大值點(diǎn)對(duì)應(yīng)的高度邊界層厚度[5-7]。常用的梯度計(jì)算方法為一階導(dǎo)數(shù)法或二階導(dǎo)數(shù)法。計(jì)算激光雷達(dá)回波距離校正信號(hào)關(guān)于高度的二階導(dǎo)數(shù)，并計(jì)算其絕對(duì)最小值，最小值對(duì)應(yīng)的高度定義為過(guò)渡層，即混合層和自由對(duì)流層交界的中間位置。一階導(dǎo)數(shù)的最小值對(duì)應(yīng)高度是混合層的頂部。
激光雷達(dá)接收到的米散射回波信號(hào)可以式（1）所示的激光雷達(dá)方程表示：
P（z）=P■CY（z）z■[β■（z）+β■（z）]T■■（z）T■■（z）（1）
式中：P（z）是激光雷達(dá)接收到距離z處氣溶膠粒子和空氣分子的后向散射回波信號(hào)；P■是激光發(fā)射功率；C是激光雷達(dá)系統(tǒng)常數(shù)；Y（z）是幾何因子；β■（z）和β■（z）分別是大氣氣溶膠粒子和大氣分子后向散射系數(shù)，T■■（z）=exp-■？鄣■（z′）dz′和T■■（z）=exp-■？鄣■（z′）dz′分別是大氣氣溶膠粒子與大氣分子透過(guò)率；？鄣■（z）和？鄣■（z）分別為大氣氣溶膠粒子和大氣分子消光系數(shù)。
由（1）式可得：
P（z）z■=P■CY（z）[β■（z）+β■（z）]T■■（z）T■■（z）（2）
激光雷達(dá)距離平方校正回波信號(hào)P（z）z■在一定程度上反映了大氣氣溶膠濃度隨探測(cè)高度變化的情況。由于覆蓋逆溫的作用，大量的大氣氣溶膠粒子富集在大氣邊界層以內(nèi)，這樣大氣邊界層大自由大氣層之間的大氣氣溶膠濃度就會(huì)發(fā)生變化。P（z）z■的梯度D（z）的變化即代表著大氣氣溶膠濃度梯度的變化[8]。
D（z）定義為D（z）=d[P（z）z■]/dz該梯度變化的最大位置就是大氣邊界層的高度。
2 EMD方法在信號(hào)處理中的應(yīng)用
經(jīng)驗(yàn)?zāi)Ｊ椒纸猓‥MD）方法的本質(zhì)是對(duì)信號(hào)進(jìn)行平穩(wěn)化處理，根據(jù)特征時(shí)間尺度把一個(gè)復(fù)雜的數(shù)據(jù)序列分解成有限項(xiàng)內(nèi)在模函數(shù)（Instrinsic Mode Functions， IMF）[9-10]。最低頻率的IMF分量通常情況下代表原始數(shù)據(jù)的整體趨勢(shì)或均值，確定IMF的標(biāo)準(zhǔn)為：（1）整個(gè)數(shù)據(jù)集中，極值點(diǎn)的數(shù)目和跨零點(diǎn)的數(shù)目必須相等或至多只差一個(gè)；（2）任一點(diǎn)上，由局部極值定義的上下包絡(luò)線的均值必須等于零。
作為一種應(yīng)用，EMD 分解方法可以有效地提取一個(gè)數(shù)據(jù)序列的趨勢(shì)或去掉該數(shù)據(jù)序列的均值。IMFS是通過(guò)篩選過(guò)程從數(shù)據(jù)中提取出的：找出原序列x（t）的各個(gè)局部極大值、極小值，用三階樣條函數(shù)分別進(jìn)行插值，相應(yīng)得到原序列x（t）的上包絡(luò)序列值xmax（t）及下包絡(luò)序列值xmin（t）；上下包絡(luò)線的均值被指定為m10；原序列x（t）和均值m10的差值h10是第一個(gè)組成。
x（t）-m10=h10
理想情況下，h10應(yīng)該是一個(gè)IMF，但由于非線性數(shù)據(jù)包絡(luò)線的平均值可能不同于實(shí)際的局部均值，導(dǎo)致非對(duì)稱波的存在，因此重復(fù)這一篩選過(guò)程就非常必要了。把h10當(dāng)作原序列，重復(fù)以上處理過(guò)程k次，直至滿足內(nèi)在模函數(shù)的定義，求出第一個(gè)內(nèi)在模函數(shù)為止。
h1k=h1（k-1）-m1k；
c1=h1k；
第一個(gè)IMF分量c1應(yīng)包含原始信號(hào)中最短的周期分量。第二個(gè)IMF函數(shù)c2通過(guò)對(duì)余數(shù)項(xiàng)r1=x（t）-c1進(jìn)行篩選得到。通常第n個(gè)IMF分量通過(guò)篩選rn=rn-1-cn得到。篩選過(guò)程直到下列任一標(biāo)準(zhǔn)滿足：（1）當(dāng)分量cn或剩余分量rn比預(yù)定值小；（2）當(dāng)剩余分量rn變成單調(diào)函數(shù)時(shí)，從中不能再篩選出基本模式分量。對(duì)于有趨勢(shì)的數(shù)據(jù)，其剩余分量就應(yīng)該是該趨勢(shì)。如果把分解后的各分量合并起來(lái)，就得到原序列x（t）=■ci+rn。
EMD分解方法的特點(diǎn)是在分解的過(guò)程中，內(nèi)在模函數(shù)保留了原始數(shù)據(jù)本身的特性，各項(xiàng)內(nèi)在模函數(shù)的和等于原始數(shù)據(jù)序列，這就允許可以只用IMFS來(lái)分析和處理數(shù)據(jù)，保留了數(shù)據(jù)本身的特性。
3 EMD方法在信號(hào)處理中的應(yīng)用
本次實(shí)驗(yàn)中，激光雷達(dá)以脈沖能量50mJ出射激光，重復(fù)頻率1Hz，所接收到的大氣回波信號(hào)經(jīng)6分鐘時(shí)間的平均后存儲(chǔ)一次，距離分辨率為15m。
圖1為2014年11月03日晚20：00的大氣回波信號(hào)（未校準(zhǔn)），考慮到大氣邊界層的高度通常為2km以下，通常我們只需處理從地面到2km這一段高度的數(shù)據(jù)。但本文為了使比較結(jié)果更形象明顯，故對(duì)采樣數(shù)據(jù)全部進(jìn)行處理。
圖1 激光雷達(dá)回波信號(hào)（2014/11/03 20：00）
對(duì)該段數(shù)據(jù)先進(jìn)行預(yù)處理，接著用上述介紹過(guò)的EMD方法分解，按照頻率由高到低得到16個(gè)IMFS，見(jiàn)圖2。從各個(gè)內(nèi)在模函數(shù)可以看出，其分解出的前6項(xiàng)IMF代表原始數(shù)據(jù)中較高頻的組成，對(duì)信號(hào)的能量貢獻(xiàn)最少，我們可以考慮將其消除，只處理低頻的部分。圖3為消除高頻項(xiàng)后的信號(hào)與原始信號(hào)的比較。通過(guò)比較可發(fā)現(xiàn)，二者基本吻合，消除高頻項(xiàng)后的信號(hào)較原始信號(hào)平滑的多，且仍能顯著地描述出原始信號(hào)的變化趨勢(shì)。將經(jīng)過(guò)以上處理的信號(hào)采用Menut求取二階導(dǎo)數(shù)的方法，認(rèn)為二階導(dǎo)數(shù)最小值所對(duì)應(yīng)的高度為大氣邊界層高度。從圖4可以看出，由EMD預(yù)先處理再求二階導(dǎo)數(shù)，信號(hào)強(qiáng)度最小值對(duì)應(yīng)的高度為1290m。
圖2 距離校準(zhǔn)信號(hào)中分解出的IMFS
在此本文把通過(guò)以上方法求得的邊界層高度值與根據(jù)常規(guī)氣象資料計(jì)算出的值比較，邊界層的計(jì)算公式采用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。利用青島氣象臺(tái)提供的氣象資料及以上關(guān)系式，我們得到11月03日18：00時(shí)的大氣邊界層厚度值為1053m。
下面，我們把2014年10月和11月測(cè)量的部分結(jié)果在表中作比較，同時(shí)表中亦列出僅根據(jù)Menut的距離校準(zhǔn)信號(hào)的二階微分最小確定的大氣邊界層厚度。
表4 兩種方法提取邊界層厚度的結(jié)果比較
從以上實(shí)驗(yàn)結(jié)果可計(jì)算出：激光雷達(dá)回波信號(hào)用EMD處理結(jié)合二階導(dǎo)數(shù)最小值對(duì)應(yīng)的大氣邊界層高度與氣象資料所確定的高度值的相關(guān)性高達(dá)76.32%，依據(jù)二階微分最小確定的大氣邊界層高度與氣象資料所確定的高度值的相關(guān)性為52.11%。顯然，前者要絕對(duì)優(yōu)于后者。也就是說(shuō)，用EMD結(jié)合二階微分最小方法確定大氣邊界層高度相比較與僅利用二階微分最小的方法來(lái)確定大氣邊界層高度，結(jié)果得到明顯改善，進(jìn)一步證明了經(jīng)驗(yàn)?zāi)Ｊ椒纸夥椒ㄔ诖髿膺吔鐚雍穸鹊募す饫走_(dá)提取應(yīng)用中的可行性。同時(shí)，由于邊界層高度的診斷關(guān)系式受探測(cè)時(shí)大氣邊界層穩(wěn)定度的影響，所以與EMD結(jié)合二階微分最小確定邊界層高度的方法相比，要相對(duì)復(fù)雜得多。
EMD結(jié)合二階微分最小確定邊界層高度的方法，僅根據(jù)激光雷達(dá)所接收到的回波信號(hào)（與大氣中氣溶膠含量成正比），利用經(jīng)驗(yàn)?zāi)Ｊ椒纸馀cMenut的二階微分求導(dǎo)二者的結(jié)合，即可得到大氣邊界層高度值，方法簡(jiǎn)單而有效。
【參考文獻(xiàn)】
[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]王珍珠，李炬，鐘志慶，等.激光雷達(dá)探測(cè)北京城區(qū)夏季大氣邊界層[J].應(yīng)用光學(xué)，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].
[責(zé)任編輯：楊玉潔]
