CN101017201A - 基于经验模态分解的激光雷达信号处理方法 - Google Patents

Abstract

本发明公开了一种基于经验模态分解的激光雷达信号处理方法，属于激光雷达信号处理领域。针对激光雷达信号和经验模态分解方法的特点，结合信号的小波变换理论和统计检测理论，利用中心极限定理和莱以特准则对信号强度起伏的随机性进行检测；检测出的信号逆变部分视为信号的发展趋势保留不参与经验模态分解方法分解，其余部分则使用经验模态分解方法处理。

Description

基于经验模态分解的激光雷达信号处理方法

技术领域

本发明涉及激光雷达信号处理方法。尤其是，适用于激光雷达信号存在逆变的情况。

技术背景

光电探测器接收到的激光雷达信号，激光雷达方程决定了信号的主要因素，激光雷达方程可以表示为如下形式：

$P\left(r\right)=P_t{\mathrm{kr}}^{-2}[{\mathrm{\β}}_m\left(r\right)+{\mathrm{\β}}_a\left(r\right)]\exp \{-2{\∫}_0^r[{\mathrm{\α}}_m\left(r^{\′}\right)+{\mathrm{\α}}_a\left(r^{\′}\right)\left]{\mathrm{dr}}^{\′}\right\}---\left(1\right)$

式中P_t为激光发射功率(W)，k是激光雷达系统常数(W·km³·Sr)，β_a(r)和β_m(r)分别是距离r处大气气溶胶粒子和空气分子后向散射系数(km^-1Sr^-1)，α_a(r)和α_m(r)分别为距离r处大气气溶胶粒子和空气分子消光系数(km^-1)。在总体上激光雷达信号呈现与距离的平方反比衰减趋势；而在局部，则由于大气不稳定性以及噪声(光电探测器的各种噪声以及天空背景辐射等)的影响引起一定的随机不可预测的起伏。对直接探测激光雷达来说，可以通过多发脉冲累计的办法削弱信号的局部起伏，提高信噪比；但在高层，由于回波信号较弱，信噪比仍然偏低。同时，由于可能的大气粒子分层分布明显和云层的存在，导致信号的局部逆变现象。当系统探测距离分辨率不够高时，这种现象尤为明显。定义：

p_corr(r)＝P(r)r²    (2)

p_corr(r)称为距离校正信号(下文简称信号)。由于空间大气的非均匀、非稳定性，激光雷达信号一般属于非平稳信号。

目前，对该类信号的滤波方法主要有中值滤波和小波变换滤波。中值滤波方法是一种非线性滤波方法(也是较普遍采用的方法)，提供了信号的简单、粗略滤波。它的基本原理是把离散信号序列中一点的值，用该点邻域中各点值的中值来替代。存在的问题是，若选择的邻域过小，滤波效果差；邻域过大，势必将信号的逆变部分平滑。小波变换是较新的信号分析理论，将信号进行“数学显微镜”式的层层分解。根据具体信号选择合适的小波基函数和分解层数，并选择合适的阈值规则对分解信号重构，实现信号滤波的目的。当信号存在较强的逆变时，会产生pseudo-Gibbs现象(即重构后的信号在突变位置附近的领域内出现局部振荡)，小波分析方法存在最优基函数的选择问题，自适应性差。由于小波基函数尺度的有限性，会导致能量泄漏问题。经验模态分解(Empirical ModeDecomposition，EMD)方法具有小波多尺度(多分辨率)分析的特点，信号分解时不需要构造基函数，其基函数由数据自身构造，是一种特殊的自适应信号分解方法。EMD方法分解的基本原理是将一列信号X(t)分解为一系列本征模态函数(IMF)。这些函数满足：1)极大值和极小值个数之和与过零点的个数之差不超过1。2)分别由极大值和极小值构成的上、下包络的平均值应处处等于0或者接近0。第一个条件类似于平稳窄带高斯过程要求；第二个条件，修正了全局性要求，以保证瞬时频率(Instantaneous Frequency)不包含不对称波形造成的不必要的波动。每一个本征模可以认为是一个新的信号。EMD方法分解算法如下：

1.初始化：R(t)＝X(t)，i＝1；

2.提取第i个IMF(IMF_i(t))：

(a).初始化：HF(t)＝R(t)；

(b)TEMP(t)＝HF(t)：

(c)分别提取TEMP(t)的极大值和极小值。采用三次样条插值求出极大值和极小值构成的上、下包络。计算上、下包络的平均值M(t)

(d).HF(t)＝TEMP(t)-M(t)；

(e).如果满足本征模的条件，置IMF_i(t)＝HF(t)，否则回到步骤(b)；

3.R(t)＝R(t)-IMF_i(t)；

4..如果R(t)最多只有一个极值或者满足本征模分解个数的要求，则退出分解，R(t)是残余项；否则返回步骤2，并置i＝i+1

一般可用下面的公式作为IMF筛分的判据：

$\mathrm{SD}=\underset{t=0}{\overset{T}{\mathrm{\Σ}}}\frac{{[H_{j-1}\left(t\right)-H_j\left(t\right)]}^2}{{H_{j-1}}^2\left(t\right)}---\left(3\right)$

门限SD一般取0.2～0.3。若要将所用的IMF分解出来，则要求R(t)中的局部极值个数不超过2个，作为提取IMF个数的终止标准。重构信号可由下式表示：

$X\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{IMF}}_i\left(t\right)+R\left(t\right)---\left(4\right)$

对给定的信号ψ(t)，对其进行Hilbert变换：

$\mathrm{\γ}\left(t^{\′}\right)=\frac{1}{\mathrm{\π}}P\underset{R}{\∫}\frac{\mathrm{\ψ}\left(t^{\′}\right)}{t-t^{\′}}{\mathrm{dt}}^{\′---\left(5\right)}$

P为柯西主值，构造信号ψ(t)的解析信号x(t)：

x(t)＝ψ(t)+iγ(t)＝a(t)e^iθ(t)    (6)

式中，

$a\left(t\right)=\sqrt{{\mathrm{\ψ}}^2\left(t\right)+{\mathrm{\γ}}^2\left(t\right)}---\left(7\right)$

$\mathrm{\θ}\left(t\right)=\arctan \left(\frac{\mathrm{\γ}\left(t\right)}{\mathrm{\ψ}\left(t\right)}\right)---\left(8\right)$

a(t)和θ(t)分别是x(t)的瞬时振幅和位相。因此，x(t)表征了信号ψ(t)瞬时特征，定义瞬时频率ω(t)：

$\mathrm{\ω}\left(t\right)=\frac{\mathrm{d\θ}\left(t\right)}{\mathrm{dt}}---\left(9\right)$

ω(t)描述了位相θ(t)的变化情况，反映了信号ψ(t)的频谱随时间的变化情况，Hilbert时频谱反映了信号的重要特征。高频IMF提取了原始信号的最精细成分，噪声一般主要包含在高频IMF中，因此，可以通过减去分解得到的高频IMF，从而实现滤波的目的。

但是当信号存在较强的逆变时，或者是信号的特征尺度发生跃变时，直接使用EMD方法分解信号，会产生模态混叠现象。模态混叠现象是指分解得到一个IMF中包含以后筛分得到的IMF的信息，由于模态混叠现象的存在，使得本征模不能够清晰的表现出信号的频率过程和信号的内在性质。因此，当产生混态模现象时，高频IMF中必定包含了逆变高频信号的成分，直接使用EMD方法分解信号会造成重构信号的失真。可能存在的逆变信号(例如卷云)会影响EMD方法在激光雷达信号滤波的效果。

发明内容

本发明解决了EMD方法应用于激光雷达信号滤波所带来的上述缺点，尤其是当信号存在较强的逆变情况下，提供了一种即能够取得较好的滤波效果又能够较好的保留逆变信息的信号处理方法。

本发明的技术方案如下：

一种基于经验模态分解的激光雷达信号处理方法，其特征在于包括以下步骤：

(f)将经A/D转换采集到激光雷达信号乘以距离的平方得到距离校正信号f(t)；

(g)选择合适的小波基函数和伸缩尺度因子，利用

${\mathrm{WT}}_f(a,b)=\frac{1}{\sqrt{\left|a\right|}}\underset{R}{\∫}f\left(t\right){\mathrm{\ψ}}^{*}\left(\frac{t-b}{a}\right)\mathrm{dt}$

其中，ψ(t)为小波函数，a为伸缩尺度因子，b为平移量，

为ψ(t)经伸缩平移后的共轭。

选择合适的尺度因子a₀，如果有：

${\frac{\∂{\mathrm{WT}}_f(a_0,b)}{\∂b}|}_{b=b_0}=0$

则称小波变换WT_f(a₀，b)在(a₀，b₀)有局部极值。若在b₀的某一领域δb₀，b∈δb₀，均有：

|WT_f(a₀，b)|≤|WT_f(a₀，b₀)|

并且在左领域或右领域严格满足：

|WT_f(a₀，b)|＜|WT_f(a₀，b₀)|

称|WT_f(a₀，b₀)|是小波变换|WT_f(a₀，b)|在尺度a₀下的模极大值点。

将(a)中的离散数据代入上述公式中计算，得到反映激光雷达信号细节信息的WT_f(a₀，b)数据序列；

(h)经中心极限定理，即由(b)中WT_f(a₀，b)数据序列组成序列{x_n}，定义部分和

$S_j=\underset{i=1}{\overset{j}{\mathrm{\Σ}}}{x}_i(j=\mathrm{1,2},...n...)$

定义归一化随机变量：

$z_n=\frac{S_n-E\left(S_n\right)}{\mathrm{\σ}\left(S_n\right)}$

式中，E(·)为平均运算符；σ(S_n)为序列{S_n}的标准偏差。

计算，将(b)中的WT_f(a₀，b)数据序列变换为一个趋于正态分布的{z_n}序列；

(i)采用莱以特准则对由步骤(c)得到的{z_n}进行随机性检测，如果有异常值，记录其在对应的步骤(b)中激光雷达信号细节信息序列中的位置；同时剔除该位置处所对应的小波系数，组成的新序列替换由步骤(b)得到的细节信息序列，回到步骤(c)；若无异常值，则转到(e)；

(j)若经(d)检测，没有异常值，则直接对f(t)采用EMD方法处理。若有异常值，则找到其对应的模极值所对应的局部领域，模极值所在的局部领域内的信号不参与EMD方法处理，信号的其余部分采用EMD方法处理。

本发明相关原理和实施步骤如下：

1.信号逆变信息检测原理介绍

(a)连续小波变换与模极大值

小波变换是将信号与一个时域和频域都具有局部化特征的平移伸缩小波基函数进行卷积，将信号分解成位于不同时频带的各个成分，这种变换有利于提取信号的本质特征。小波基函数的定义为，设ψ(t)为平方可积函数，即ψ(t)∈L²(R)。如果满足：

$C_{\mathrm{\&psi;}}=\underset{R}{\&Integral;}\frac{{\left|\widehat{\mathrm{\&psi;}}\left(\mathrm{\&omega;}\right)\right|}^2}{\left|\mathrm{\&omega;}\right|}\mathrm{d\&omega;}<+\&infin;---\left(10\right)$

则称ψ(t)为基本小波函数或母小波函数。将任意L²(R)空间中的函数(信号)f(t)与小波函数的卷积，称为f(t)的连续小波变换：

${\mathrm{WT}}_f(a,b)=\frac{1}{\sqrt{\left|a\right|}}\underset{R}{\&Integral;}f\left(t\right){\mathrm{\&psi;}}^{*}\left(\frac{t-b}{a}\right)\mathrm{dt}---\left(11\right)$

其中，a为伸缩尺度因子，b为平移量。

由小波函数的空间局部化特征可知，WT_f(a，b)的值主要决定于信号f(t)在b处领域内的值，也即信号在某点处的小波变换在尺度a下完全由该点附近的局部信息所确定；并且a越小，领域区间也越小。因此，在合适的尺度(a)上，WT_f(a，b)提供了信号的局部信息。

给定尺度a₀，如果有：

${\frac{{\&PartialD;\mathrm{WT}}_f(a_0,b)}{\&PartialD;b}|}_{b=b_0}=0---\left(12\right)$

则称小波变换WT_f(a₀，b)在(a₀，b₀)有局部极值。若在b₀的某一领域δb₀，∈δb₀，均有：

|WT_f(a₀，b)|≤|WT_f(a₀，b₀)|    (13)

并且在左领域或右领域严格满足：

|WT_f(a₀，b)|＜|WT_f(a₀，b₀)|    (14)

称|WT_f(a₀，b₀)|是小波变换|WT_f(a₀，b)|在尺度a₀下的极大点。

小波变换的模极大与信号的局部突变点(奇异性)有关。而信号局部突变主要由各种噪声和局部逆变信息引起的。因此，可以通过结合中心极限定理和莱以特准则来检测WT_f(a，b)以及模极大值，从而可以检测逆变信息。

(b)中心极限定理

设{x_n}是独立随机变量序列，取：

$S_j=\underset{i=1}{\overset{j}{\mathrm{\&Sigma;}}}{x}_i(j=\mathrm{1,2},...n...)---\left(15\right)$

定义归一化随机变量：

$z_n=\frac{S_n-E\left(S_n\right)}{\mathrm{\&sigma;}\left(S_n\right)}---\left(16\right)$

式中，E(·)为平均运算符；σ(S_n)为序列{S_n}的标准偏差。对于有限均值和方差的统计独立同分布的随机变量，中心极限定理可以表述为：在独立随机变量序列中，每个随机变量x_i对归一化随机变量z_n的影响足够小，当n充分大时，{z_n}是收敛于标准正态分布的随机变量序列。中心极限定理的必然结果是：如果一个物理过程(例如电路噪声)为许多独立作用之和，并且满足独立同分布，均值和方差有限的条件，那么这个过程就趋于正态(高斯)过程。

(3)莱以特准则(3σ准则)

对于一个正态分布的随机变量序列，可以采用莱以特准则对序列的正态性进行检测(当序列的长度大于500时，检测效果更好)。由正态分布的相关理论，变量的残余误差落在±3σ以外的概率约为0.3％。计算随机变量序列中变量的残余误差，如果残余误差大于3，则可认为该变量是异常值。

2.具体实施步骤

(k)选择合适的小波基函数和伸缩尺度因子，利用小波变换(11)式提取信号的细节信息；

(l)经中心极限定理，由(15)式和(16)式，将步骤(a)中由小波变换得到的信号细节信息变换为一个趋于正态分布的序列；

(m)采用莱以特准则对由步骤(b)得到的序列进行随机性检测，如果有异常值，记录其在对应的步骤(a)中细节信息序列中的位置；同时剔除该位置处所对应的小波系数，组成的新序列替换由步骤(a)得到的细节信息序列，回到步骤(b)；若无异常值，则转到(d)；

(n)若无异常值，采用EMD方法对信号进行处理；若有异常值，由(12)式、(13)式和(14)式，找到异常值对应的模极值对应的逆变信息局部领域。局部领域内的信号不参与EMD方法分解，信号按照局部领域分段采用EMD方法处理。

对激光雷达信号采用分段处理的目的有三个：(1)由于云层回波信号和边界层回波信号往往比低信噪比信号的强度高几个量级，若直接使用EMD方法对信号进行滤波处理，在信号的强逆变领域可能产生严重的混态模现象，使得重构信号失真，滤波失去意义。(2)EMD方法在信号分解时，需要运用三次样条插值来获得信号的上、下包络，在数据比较多的情况下需要运算高阶矩阵，占用较多的内存，耗时多。(3)信号的分段处理有利于保留信号逆变位置携带的信息。

本发明结合小波分析理论和统计检测理论，利用中心极限定理和莱以特准则对由连续小波变换分解信号得到的系数进行随机性检测，检测出的逆变信息局部领域(一般是分层部分以及云层回波信号)所对应的信号可视为信号的发展趋势保留不参与EMD方法分解，而其余部分则分段使用EMD方法分解。本发明即能够取得较好的滤波效果，又能够较好的保留信号的细节信息和逆变信息。另外，EMD方法在信号分解时，需要运用三次样条插值来获得信号的上、下包络，在数据比较多的情况下需要运算高阶矩阵，占用较多的内存，耗时多。本发明通过分段处理信号，减少矩阵的阶数，提高处理速度。

附图说明

图1为实测的激光雷达信号(有卷云)。

图2为实测的激光雷达信号(无卷云)。

图3为图2所示的信号通过EMD方法分解得到的IMF和残余项R(t)。

图4为图3所示的IMF经Hilbert变换得到的时频谱。

图5为图3所示的IMF经离散Fourier变换得到的频谱。

由图3和图5可知，高频IMF主要体现了噪声的特征。

图6为无卷云时，信号经EMD方法滤波后的效果图。

图7为当信号存在较强的逆变信号时，有关模态混叠现象的示意图。

图8和图9分别为无卷云和有卷云时对逆变信号的检测结果。

图10为有卷云时，采用EMD方法和分段处理信号的示意图。

具体实施方式

参见附图。

如图1所示，信号强度分层分布明显，边界层回波信号强，起伏少，信噪比高；2～3km和4～5km之间分别有一个干净层，7km的地方出现了一层薄薄的卷云；由于卷云的作用，7km之后的区域基本上是噪声。

如图4所示，在Hilbert时频谱中，能够同时在时域和频域中考察信号的局部和整体特征。由图可见，激光雷达信号属于低频信号，图2所示的信号其能量主要分布在边界层至6km范围内。

图7中的(a)图显示的是信号经EMD方法滤波后的效果，由图可见，由于卷云的存在，使用EMD方法分解信号时，产生了混态模现象如(b)图所示，导致了重构信号的失真，滤波失去意义。

在信号逆变区域，小波变换得到的系数均具有良好的对应性。若认为信号强度起伏变化独立随机，信号中逆变信息的出现属于小概率事件，噪声则属于大概率事件。由于激光发射脉冲能量的有限性，将信号经小波变换的系数通过中心极限定理转换，再结合莱以特准则和模极大值的条件，对逆变信息进行检测。无卷云和有卷云的检测效果如图8和图9所示。

由于采用分段及多尺度EMD分解方法处理信号，信号的逆变位置可能携带的重要信息得到了很好的保留，同时又取得了满意的滤波效果；典型的处理效果如图10所示。这对反演空间大气随高度分布的直接、间接特征参数具有重要的意义。

Claims (1)

1.一种基于经验模态分解的激光雷达信号处理方法，其特征在于包括以下步骤：

(a)将经A/D转换采集到激光雷达信号乘以距离的平方得到距离校正信号f(t)；

(b)选择合适的小波基函数和伸缩尺度因子，利用

$W{T}_f(a,b)=\frac{1}{\sqrt{\left|a\right|}}\underset{R}{\&Integral;}f{\left(t\right)\mathrm{\&psi;}}^{*}\left(\frac{t-b}{a}\right)\mathrm{dt}$

其中，ψ(t)为小波函数，a为伸缩尺度因子，b为平移量， 为ψ(t)经伸缩平移后的共轭。

选择合适的尺度因子a₀，如果有：

$\frac{\&PartialD;W{T}_f(a_0,b)}{\&PartialD;b}{|}_{b=b_0}=0$

则称小波变换WT_f(a₀，b)在(a₀，b₀)有局部极值。若在b₀的某一领域δb₀，b∈δb₀，均有：

|WT_f(a₀，b)|≤|WT_f(a₀，b₀)|

并且在左领域或右领域严格满足：

|WT_f(a₀，b)|＜|WT_f(a₀，b₀)|

称|WT_f(a₀，b₀)|是小波变换|WT_f(a₀，b)|在尺度a₀下的模极大值点。

将(a)中的离散数据代入上述公式中计算，得到反映激光雷达信号细节信息的WT_f(a₀，b)数据序列；

(c)经中心极限定理，即由(b)中WT_f(a₀，b)数据序列组成序列{x_n}，定义部

分和

$S_j=\underset{i=1}{\overset{j}{\mathrm{\&Sigma;}}}{x}_i,(j=\mathrm{1,2},...n...)$

定义归一化随机变量：

$z_n=\frac{S_n-E\left(S_n\right)}{\mathrm{\&sigma;}\left(S_n\right)}$

式中，E(·)为平均运算符；σ(S_n)为序列{S_n}的标准偏差。计算，将(b)中的WT_f(a₀，b)数据序列变换为一个趋于正态分布的{z_n}序列；

(d)采用莱以特准则对由步骤(c)得到的{z_n}进行随机性检测，如果有异常值，记录其在对应的步骤(b)中激光雷达信号细节信息序列中的位置；同时剔除该位置处所对应的小波系数，组成的新序列替换由步骤(b)得到的细节信息序列，回到步骤(c)；若无异常值，则转到(e)；

(e)若经(d)检测，没有异常值，则直接对f(t)采用EMD方法处理。若有异常值，则找到其对应的模极值所对应的局部领域，模极值所在的局部领域内的信号不参与EMD方法处理，信号的其余部分采用EMD方法处理。

Images