基于经验模式分解的大气边界层高度提取

张娜　姜茂仁　李进杰

【摘 要】研究大气边界层的实时高度是研究天气、气候和大气环境污染的迫切需要。激光雷达成为探测边界层时空演变特征的最有效手段。目前国际上用于反演激光雷达数据获取大气边界层高度的主要方法有梯度法、标准偏差法、曲线拟合法等。本文采用偏振拉曼—米散射激光雷达的实测数据，介绍梯度法中的二阶求导反演激光雷达数据提取大气边界层高度，并采用经验模式分解（EMD）方法处理数据。通过对结果进行分析，发现提取的大气边界层高度值与利用诊断式求得的高度值的吻合性比较理想。

【关键词】大气边界层；激光雷达；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公里范围内。大气边界层是与人类关系密切的大气层，是大气污染物的主要集中层，边界层高度的大小，决定了排入大气中污染物的浓度值。所以针对气候及大气环境的研究需要，开展对大气边界层高度的观测势在必行。同时，对边界层高度的观测要求也越来越高。

激光雷达作为一种新型大气主动遥感探测工具[2-4]，已经广泛应用于大气气溶胶和污染物、云等领域的测量和研究工作。激光雷达利用气溶胶作为示踪粒子，所测量到的光功率信号正比于大气的气溶胶含量。与此相对应，大气边界层中含有丰富的气溶胶（较自由对流层的气溶胶含量更大），这是因为气溶胶粒子主要来源于地球表面，从而导致对激光的更大散射，因此可以利用激光雷达较容易的探测到这两层之间的边界。

本文利用激光雷达探测系统饿实测数据，采用了梯度法中的二阶求导法来提取大气边界层高度，并在分析处理中应用了经验模式分解（EMD）方法，给出了典型天气情况下的提取结果，并与诊断式求得的高度值进行比对，吻合度较高。

1 梯度法反演大气边界层高度

因为边界层顶与自由对流层气溶胶浓度差距较大，所以回波信号的梯度最大值点对应的高度边界层厚度[5-7]。常用的梯度计算方法为一阶导数法或二阶导数法。计算激光雷达回波距离校正信号关于高度的二阶导数，并计算其绝对最小值，最小值对应的高度定义为过渡层，即混合层和自由对流层交界的中间位置。一阶导数的最小值对应高度是混合层的顶部。

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

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

式中：P（z）是激光雷达接收到距离z处气溶胶粒子和空气分子的后向散射回波信号；P■是激光发射功率；C是激光雷达系统常数；Y（z）是几何因子；β■（z）和β■（z）分别是大气气溶胶粒子和大气分子后向散射系数，T■■（z）=exp-■？鄣■（z′）dz′和T■■（z）=exp-■？鄣■（z′）dz′分别是大气气溶胶粒子与大气分子透过率；？鄣■（z）和？鄣■（z）分别为大气气溶胶粒子和大气分子消光系数。

由（1）式可得：

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

激光雷达距离平方校正回波信号P（z）z■在一定程度上反映了大气气溶胶浓度随探测高度变化的情况。由于覆盖逆温的作用，大量的大气气溶胶粒子富集在大气边界层以内，这样大气边界层大自由大气层之间的大气气溶胶浓度就会发生变化。P（z）z■的梯度D（z）的变化即代表着大气气溶胶浓度梯度的变化[8]。

D（z）定义为D（z）=d[P（z）z■]/dz该梯度变化的最大位置就是大气边界层的高度。

2 EMD方法在信号处理中的应用

经验模式分解（EMD）方法的本质是对信号进行平稳化处理，根据特征时间尺度把一个复杂的数据序列分解成有限项内在模函数（Instrinsic Mode Functions， IMF）[9-10]。最低频率的IMF分量通常情况下代表原始数据的整体趋势或均值，确定IMF的标准为：（1）整个数据集中，极值点的数目和跨零点的数目必须相等或至多只差一个；（2）任一点上，由局部极值定义的上下包络线的均值必须等于零。

作为一种应用，EMD 分解方法可以有效地提取一个数据序列的趋势或去掉该数据序列的均值。IMFS是通过筛选过程从数据中提取出的：找出原序列x（t）的各个局部极大值、极小值，用三阶样条函数分别进行插值，相应得到原序列x（t）的上包络序列值xmax（t）及下包络序列值xmin（t）；上下包络线的均值被指定为m10；原序列x（t）和均值m10的差值h10是第一个组成。

x（t）-m10=h10

理想情况下，h10应该是一个IMF，但由于非线性数据包络线的平均值可能不同于实际的局部均值，导致非对称波的存在，因此重复这一筛选过程就非常必要了。把h10当作原序列，重复以上处理过程k次，直至满足内在模函数的定义，求出第一个内在模函数为止。

h1k=h1（k-1）-m1k；

c1=h1k；

第一个IMF分量c1应包含原始信号中最短的周期分量。第二个IMF函数c2通过对余数项r1=x（t）-c1进行筛选得到。通常第n个IMF分量通过筛选rn=rn-1-cn得到。筛选过程直到下列任一标准满足：（1）当分量cn或剩余分量rn比预定值小；（2）当剩余分量rn变成单调函数时，从中不能再筛选出基本模式分量。对于有趋势的数据，其剩余分量就应该是该趋势。如果把分解后的各分量合并起来，就得到原序列x（t）=■ci+rn。

EMD分解方法的特点是在分解的过程中，内在模函数保留了原始数据本身的特性，各项内在模函数的和等于原始数据序列，这就允许可以只用IMFS来分析和处理数据，保留了数据本身的特性。

3 EMD方法在信号处理中的应用

本次实验中，激光雷达以脉冲能量50mJ出射激光，重复频率1Hz，所接收到的大气回波信号经6分钟时间的平均后存储一次，距离分辨率为15m。

图1为2014年11月03日晚20：00的大气回波信号（未校准），考虑到大气边界层的高度通常为2km以下，通常我们只需处理从地面到2km这一段高度的数据。但本文为了使比较结果更形象明显，故对采样数据全部进行处理。

图1 激光雷达回波信号（2014/11/03 20：00）

对该段数据先进行预处理，接着用上述介绍过的EMD方法分解，按照频率由高到低得到16个IMFS，见图2。从各个内在模函数可以看出，其分解出的前6项IMF代表原始数据中较高频的组成，对信号的能量贡献最少，我们可以考虑将其消除，只处理低频的部分。图3为消除高频项后的信号与原始信号的比较。通过比较可发现，二者基本吻合，消除高频项后的信号较原始信号平滑的多，且仍能显著地描述出原始信号的变化趋势。将经过以上处理的信号采用Menut求取二阶导数的方法，认为二阶导数最小值所对应的高度为大气边界层高度。从图4可以看出，由EMD预先处理再求二阶导数，信号强度最小值对应的高度为1290m。

图2 距离校准信号中分解出的IMFS

在此本文把通过以上方法求得的边界层高度值与根据常规气象资料计算出的值比较，边界层的计算公式采用Benkley and Schulman提供的中性和稳定情况下的关系式：h=a×■，式中u？鄢为摩擦速度（m/s）；f为地转参数（1/s）；a为常数，由气象资料决定，不同研究者给出的a值不同。本文中，我们根据u？鄢值的范围选择a，文中有两种情况，a=0.31或a=0.56。利用青岛气象台提供的气象资料及以上关系式，我们得到11月03日18：00时的大气边界层厚度值为1053m。

下面，我们把2014年10月和11月测量的部分结果在表中作比较，同时表中亦列出仅根据Menut的距离校准信号的二阶微分最小确定的大气边界层厚度。

表4 两种方法提取边界层厚度的结果比较

从以上实验结果可计算出：激光雷达回波信号用EMD处理结合二阶导数最小值对应的大气边界层高度与气象资料所确定的高度值的相关性高达76.32%，依据二阶微分最小确定的大气边界层高度与气象资料所确定的高度值的相关性为52.11%。显然，前者要绝对优于后者。也就是说，用EMD结合二阶微分最小方法确定大气边界层高度相比较与仅利用二阶微分最小的方法来确定大气边界层高度，结果得到明显改善，进一步证明了经验模式分解方法在大气边界层厚度的激光雷达提取应用中的可行性。同时，由于边界层高度的诊断关系式受探测时大气边界层稳定度的影响，所以与EMD结合二阶微分最小确定边界层高度的方法相比，要相对复杂得多。

EMD结合二阶微分最小确定边界层高度的方法，仅根据激光雷达所接收到的回波信号（与大气中气溶胶含量成正比），利用经验模式分解与Menut的二阶微分求导二者的结合，即可得到大气边界层高度值，方法简单而有效。

【参考文献】

[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]王珍珠，李炬，钟志庆，等.激光雷达探测北京城区夏季大气边界层[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].

[责任编辑：杨玉洁]