CN102221629A - 一种基于小波变换的车辆横摆角速度滤波测量方法 - Google Patents

Abstract

本发明公开了一种基于小波变换的车辆横摆角速度滤波测量方法，适用于前轮为转向轮、后轮为非转向轮的汽车横摆角速度的测量。首先测出两个非转向轮的轮速，利用中值滤波法去除轮速粗大误差，由平面复合运动的运动关系，将通过中值滤波后的轮速进一步推算出未滤波的车辆横摆角速度，然后根据小波变换理论，选用‘sym8’小波基，对未滤波的车辆横摆角速度进行七层小波分解，并用软阈值法对高频系数进行阈值量化处理，最后用第七层低频系数和软阈值法量化处理后的高频系数进行逆小波变换，即可重构得到精度较高的滤波后的横摆角速度。

Description

一种基于小波变换的车辆横摆角速度滤波测量方法

技术领域

本发明涉及一种基于‘sym8’小波的车辆横摆角速度的滤波测量方法，其目的在于通过低成本的轮速传感器实现对横摆角速度的较高精度测量，以便评价车辆的操纵稳定性，属于汽车道路试验与测试领域。

背景技术

操纵稳定性是车辆主动安全性的重要评价指标，是现代车辆的重要使用性能之一。车辆横摆角速度，是衡量车辆操纵稳定性的重要标志。在汽车操稳试验和ESP电子稳定程序等领域都是一项重要参数。在GB/T6326-1994《汽车操纵稳定性试验方法》和QC/T480-1999《汽车操纵稳定性指标限值和评价方法》中都作了明确的规定，多项汽车操纵稳定性基本评价内容都将横摆角速度作为重要指标：横摆角速度是转向盘角阶跃输入下的响应、横摆角速度频率响应特性和回正性等基本评价内容的主要评价参数之一；典型行驶工况性能中的蛇形性能、移线性能、双移线性能，也都将横摆角速度纳为评价参量。此外，ESP电子稳定程序(Electronic Stability Program)是当今车辆主动安全控制系统的典型代表，这个系统的重要输入量和控制量之一就是横摆角速度。

根据国内外现有的研究情况，横摆角速度通常采用陀螺仪直接测量得到。但是，陀螺仪的价格昂贵，微机械陀螺仪一般价格在几千元人民币，而性能更好的光纤陀螺仪价格更是高达数万甚至几十万元人民币；并且陀螺仪安装比较麻烦，要尽量安装在车辆质心位置，限制了它的广泛使用。

实际上，横摆角速度还可以通过适当的运动关系由轮速等信号间接测量得到。这种方法一方面将有效的降低成本，另一方面在引入适当的滤波方法后还可以提高测量精度；即使在保留现有陀螺仪的情况下，也可以通过横摆角速度的估计实现对传感器的冗余故障诊断，因此得到了日益广泛的应用。

轮速信号的传感回路是弱电系统，容易受外界强电磁的干扰；而且车辆运行过程中由于路面的影响，检测装置也容易发生震动，这些都可能导致干扰信号产生，并与正常信号叠加在一起。由于这种横摆角速度测量方法主要通过轮速信号推算得到，因此，对测得的横摆角速度信号进行滤波，抑制和防止干扰，提高测量精度，具有重要的意义。信号的滤波方法主要有用傅立叶变换对信号进行时频变换方法、卡尔曼滤波和神经网络滤波等。由于干扰信号随机性较强，频率成分比较复杂，当信号和噪声的频带相互重叠时，傅立叶变换对信号进行时-频变换的方法无法在频域对其进行有效的区分，不能很好的滤除干扰信号的影响。卡尔曼滤波需要建立准确的系统状态和观测模型，并且需要知道精确的系统噪声和观测噪声，在实际操作中这些很难准确得到，导致卡尔曼滤波的精度不高。神经网络滤波需要较多的先验数据，并且需要较长的训练时间，这些制约了神经网络滤波的应用。

为克服上述方法的局限性，本发明将针对汽车操稳性试验中横摆角速度的测量要求，提供一种低成本、较高精度的横摆角速度测量方法。

发明内容

本发明提出了一种低成本、较高精度的基于‘sym8’小波的车辆横摆角速度信号的滤波测量方法，适用于前轮为转向轮、后轮为非转向轮的汽车操纵稳定性试验中横摆角速度的测量。

技术方案

本发明为实现上述目的，采用如下的技术方案：

本发明提出了一种基于小波变换的车辆横摆角速度滤波测量方法，适用于前轮为转向轮、后轮为非转向轮的汽车横摆角速度的测量，其特征在于包含如下步骤：

步骤1)首先利用汽车轮速传感器测出右后轮的轮速V_rr(i)和左后轮的轮速V_lr(i)(i＝1，2，…，n，i表示采样序数，n为数据的长度)，右后轮和左后轮均为非转向轮；

步骤2)对测得的右后轮轮速数据V_rr(i)先进行中值滤波，i＝1，2，…，n，除去粗大误差，保持V_rr(1)、V_rr(2)、V_rr(n-1)和V_rr(n)不变，设数据的采样间隔为T秒，取阈值为0.5gT，g为重力加速度，如果V_rr(j)≥V_rr(j-1)+0.5gT或者V_rr(j)≤V_rr(j-1)-0.5gT，则V_rr(j)取median{V_rr(j-2)，V_rr(j-1)，V_rr(j)，V_rr(j+1)，V_rr(j+2)}，否则V_rr(j)不变，3≤j≤n-2，median表示取中位数；

再对测得的左后轮轮速数据V_lr(i)进行中值滤波，i＝1，2，…，n，保持V_lr(1)、V_lr(2)、V_lr(n-1)和V_lr(n)不变，取阈值为0.5gT，如果V_lr(j)≥V_lr(j-1)+0.5gT或者V_lr(j)≤V_lr(j-1)-0.5gT，则V_lr(j)取median{V_lr(j-2)，V_lr(j-1)，V_lr(j)，V_lr(j+1)，V_lr(j+2)}，否则V_lr(j)不变，3≤j≤n-2；

步骤3)根据步骤2)中得到的滤波后的轮速数据，由公式

得到未滤波的横摆角速度ω(i)，i＝1，2，…，n，根据Mallat算法，并选用‘sym8’小波，对横摆角速度数据ω(i)进行7层小波分解，得到第七层低频系数

k₇＝1，2，...，m₇，

以及第一层高频系数k₁＝1，2，...，m₁，

第二层高频系数

k₂＝1，2，...，m₂，

第三层高频系数k₃＝1，2，...，m₃，第四层高频系数

k₄＝1，2，...，m₄，

第五层高频系数

k₅＝1，2，...，m₅，

第六层高频系数

k₆＝1，2，...，m₆，

和第七层高频系数

k₇＝1，2，...，m₇，

floor表示下取整运算，T_w为后轴左右轮的轮距，分解的具体过程如下：首先用初始序列

表示未滤波的横摆角速度ω(i)序列，i＝1，2，…，n，k₀＝1，2，…，m₀，m₀＝n，对序列

进行长度为16的对称延拓，得到

k′₀＝1，2，...，n+32，将延拓后的

与低通分解滤波器LO_D作卷积，然后将卷积后的结果首、尾部各去掉16个数，再对剩余的元素进行二元下抽样即可得到第一层低频系数

k₁＝1，2，...，m₁，

将延拓后的

与高通分解滤波器HI_D做卷积，然后再将卷积后的结果首、尾部各去掉16个数，再对剩余的元素进行二元下抽样即可得到第一层高频系数

k₁＝1，2，...，m₁，

同理，接着对

k₁＝1，2，...，m₁，

进行长度为16的对称延拓得到

k′₁＝1，2，...，m₁+32，将延拓后的

与低通分解滤波器LO_D做卷积，然后再将卷积后的结果首、尾部各去掉16个数，再对剩余的元素进行二元下抽样即可得到第二层低频系数

k₂＝1，2，...，m₂，将延拓后的与高通分解滤波器HI_D作卷积，然后再将卷积后的结果首、尾部各去掉16个数，再对剩余的元素进行二元下抽样即可得到第二层高频系数

k₂＝1，2，...，m₂，

按上述过程继续进行分解，依次得到第三层高频系数

k₃＝1，2，...，m₃，

第四层高频系数k₄＝1，2，...，m₄，第五层高频系数

k₅＝1，2，...，m₅，

第六层高频系数

k₆＝1，2，...，m₆，

第七层高频系数

k₇＝1，2，...，m₇，

和第七层低频系数

k₇＝1，2，...，m₇，其中低通分解滤波器LO D和高通分解滤波器HI_D的取值如下：

LO_D＝{-0.0034 -0.0005 0.0317 0.0076 -0.1433 -0.0613 0.48140.7772 0.3644 -0.0519 -0.0272 0.0491 0.0038 -0.0150 -0.00030.0019}，

HI_D＝{-0.0019 -0.0003 0.0150 0.0038 -0.0491 -0.0272 0.05190.3644 -0.7772 0.4814 0.0613 -0.1433 0.0076 0.0317 0.0005-0.0034}，

步骤4)小波分解高频系数的阈值量化：保留第七层低频系数

不变，k₇＝1，2，…，m₇，对分解得到的高频系数序列d_k，l用软阈值法进行量化处理，k＝1，2，...，7，l＝1，2...，m_k，即当d_k，l≥λ时，

当d_k，l≤-λ时，

当|d_k，l|＜λ时，其中，k＝1，2，...，7，l＝1，2...，m_k，

k₁＝1，2，...，m₁；经量化处理得到量化处理后的第一层高频系数

k₁＝1，2，...，m₁，量化处理后的第二层高频系数

k₂＝1，2，...，m₂，量化处理后的第三层高频系数

k₃＝1，2，...，m₃，量化处理后的第四层高频系数

k₄＝1，2，...，m₄，量化处理后的第五层高频系数k₅＝1，2，...，m₅，量化处理后的第六层高频系数

k₆＝1，2，...，m₆，和量化处理后的第七层高频系数

k₇＝1，2，...，m₇；

步骤5)进行逆小波变换，重构滤波后的信号：由第七层低频系数

k₇＝1，2，...，m₇，量化处理后的第一层高频系数k₁＝1，2，...，m₁，量化处理后的第二层高频系数k₂＝1，2，...，m₂，量化处理后的第三层高频系数

k₃＝1，2，...，m₃，量化处理后的第四层高频系数

k₄＝1，2，...，m₄，量化处理后的第五层高频系数k₅＝1，2，...，m₅，量化处理后的第六层高频系数

k₆＝1，2，...，m₆，以及量化处理后的第七层高频系数k₇＝1，2，...，m₇，做逆小波变换即可得到滤波后的横摆角速度信号，具体过程如下：首先让序列

k₇＝1，2，...，m₇，再将

和量化处理后的第七层高频系数

分别进行长度为16的对称延拓得到

和

其中，k₇＝1，2，...，m₇，k′₇＝1，2，…，m₇+32，再对延拓后的和延拓后的

分别做二元上抽样，k′₇＝1，2，…，m₇+32，

表示

二元上抽样的结果，

表示

二元上抽样的结果，然后将

与低通重构滤波器LO_R作卷积得到序列

p₁＝1，2，…，(2m₇+80)，

与高通重构滤波器HI_R作卷积得到序列

并将序列与序列

中相对应元素分别进行相加，得到新序列

p₁＝1，2，…，(2m₇+80)，找到

的第

项

如果步骤4)中的m₆为奇数，则在

中从

的前一项开始向前取

项，从的后一项开始向后也取

项，保留

这m₆个数据按原来的先后顺序排列就得到重构后的第六层低频系数

k₆＝1，2，...，m₆；如果步骤4)中的m₆为偶数，则在

中从

的前一项开始向前取

项，从

的后一项开始向后取项，保留

这m₆个数据按原来的顺序排列就得到重构后的第六层低频系数

k₆＝1，2，...，m₆，同理，将

和量化处理后的第六层高频系数

分别进行长度为16的对称延拓得到

和

k₆＝1，2，...，m₆，k′₆＝1，2，…，m₆+32，再分别对延拓后的

和延拓后的做二元上抽样，k′₆＝1，2，…，m₆+32，

表示

二元上抽样的结果，

表示

二元上抽样的结果，然后将

与低通重构滤波器LO_R作卷积得到序列

p₂＝1，2，…，(2m₆+80)，

与高通重构滤波器HI_R作卷积得到序列

并将序列

与序列中相对应元素分别进行相加，得到序列

p₂＝1，2，…，(2m₆+80)，找到的第

项如果步骤4)中的m₅为奇数，则从

的前一项开始向前取

项，从

的后一项开始向后也取

项，保留

这m₅个数据按原来的顺序排列就得到重构后的第五层低频系数

k₅＝1，2，...，m₅；如果步骤4)中的m₅为偶数，则从

的前一项开始向前取

项，从

的后一项向后取

项，保留

这m₅个数据按原来的顺序排列就得到重构后的第五层低频系数

k₅＝1，2，...，m₅，按上述过程继续进行重构，最终可得到

k₀＝1，2，…，m₀，m₀＝n，即为滤波后的横摆角速度数据，其中低通重构滤波器LO_R和高通重构滤波器HI_R的取值如下：

LO_R＝{0.0019 -0.0003 -0.0150 0.0038 0.0491 -0.0272 -0.05190.3644 0.7772 0.4814 -0.0613 -0.1433 0.0076 0.0317 -0.0005-0.0034}，

HI_R＝{-0.0034 0.0005 0.0317 -0.0076 -0.1433 0.0613 0.4814-0.7772 0.3644 0.0519 -0.0272 -0.0491 0.0038 0.0150 -0.0003-0.0019}。

有益效果

1本发明提出一种基于小波变换的横摆角速度滤波测量方法，适用于前轮为转向轮、后轮为非转向轮的车辆的横摆角速度测量。

2本发明提供的横摆角速度滤波测量方法，通过小波变换能够有效地滤除横摆角速度测量值中的干扰，测量精度高且数据平滑。

3本发明通过轮速数据来间接推算横摆角速度，并利用提出的基于小波变换的滤波方法来提高横摆角速度的测量精度。具体而言，采用小波收缩阈值法，不但能够有效地抑制噪声，而且可以很好的保留横摆角速度信号的特征峰值点，滤波过程中采用‘sym8’小波基进行7层小波分解，能将低频信号和高频的噪声信号有效地分离开，且计算量小。

4本发明提供的横摆角速度滤波测量方法成本低、实现简便、效率高。

附图说明

图1是本滤波方法的流程图；

图2是车辆三自由度模型；

图3是小波分解图；

图4是Mallat算法小波分解图；

图5是小波重构图；

图6是Mallat算法小波重构图；

图7是未滤波的横摆角速度、小波滤波后的横摆角速度和陀螺仪测得的横摆角速度对比图；

图8是未滤波的横摆角速度和陀螺仪测得的横摆角速度对比图；

图9是小波滤波后的横摆角速度和陀螺仪测得的横摆角速度对比图；

图10～图11是图7的局部放大图。

具体实施方式

实施实例1

操纵稳定性是汽车主动安全性的重要评价指标，是现代汽车的重要性能之一。在GB/T6326-1994《汽车操纵稳定性试验方法》和QC/T480-1999《汽车操纵稳定性指标限值和评价方法》中对汽车操纵稳定性试验方法、指标限值和评价方法都作了明确的规定。据此，可以看出汽车操纵稳定性主要评价方法有以下几个特点：

(1)转向盘角阶跃输入下的稳态响应利用横摆角速度与前轮转角之比，即稳态横摆角速度增益来反映转向灵敏度。其瞬态响应也以横摆角速度波动的无阻尼圆频率作为主要评价参数。

(2)横摆角速度频率响应特性是指转向盘转角正弦输入下，频率由0→∞变化时，汽车横摆角速度与转向盘转角的振幅比及相位差的变化规律，是另一个表征汽车操纵稳定性的重要特性。

(3)转向回正性能试验中，按松开方向盘3s时的残留横摆角速度绝对值和横摆角速度总方差作为评价指标。

(4)典型行驶工况性能中的蛇形性能、移线性能、双移线性能，都将横摆角速度纳为评价参量。

根据上述特点，可以看出横摆角速度是多项汽车操纵稳定性的主要评价指标之一。

根据国内外现有的研究情况，横摆角速度通常采用陀螺仪直接测量得到。但是，陀螺仪价格昂贵，微机械陀螺仪一般价格在几千元人民币，而性能更好的光纤陀螺仪价格更是高达数万甚至几十万元人民币；并且陀螺仪安装比较麻烦，要尽量安装在车辆质心位置，限制了它的广泛使用。实际上，横摆角速度还可以通过适当的运动关系由轮速等信号间接测量得到。具体的思路如下：

操纵稳定性试验通常在比较平整、坡度较小的路面上进行的，针对前轮是转向轮、后轮是非转向轮的车辆，此时整车的运动学模型可以简化为一个三自由度的平面复合运动模型，如图2所示，即包括沿OX方向的前向运动，沿OY方向的横向运动和横摆运动(横摆角速度为图中ω)。图2中，V是对地绝对平面速度，V_lf、V_rf、V_lr、V_rr分别代表左前轮、右前轮、左后轮、右后轮的线速度。β为侧滑角，ω为横摆角速度，T_w为后轴左右轮的轮距。根据平面复合运动的运动关系，可得

$V_{\mathrm{lr}}=V\cos \mathrm{\β}-\frac{T_w}{2}\mathrm{\ω}---\left(1\right)$

$V_{\mathrm{rr}}=V\cos \mathrm{\β}+\frac{T_w}{2}\mathrm{\ω}---\left(2\right)$

由式(1)和(2)，可以得到

$\mathrm{\ω}=\frac{V_{\mathrm{rr}}-V_{\mathrm{lr}}}{T_w}---\left(3\right)$

根据式(3)可以看出，通过测量右后轮的轮速和左后轮的轮速可以推算得到横摆角速度。据此，首先利用汽车轮速传感器测出右后轮的轮速V_rr(i)和左后轮的轮速V_lr(i)(i＝1，2…，n，i表示采样序数，n为数据的长度)，右后轮和左后轮均是非转向轮。n的确定方法为：试验过程的持续时间除以轮速数据的采样间隔T，然后再对除法的结果取整即为n。由于受道路情况、轮速检测装置震动等因素的影响，测得的原始数据含有较多的粗大误差，为除掉粗大误差，需要进行中值滤波。一般车辆在进行操纵稳定性试验时，前进方向的加速度通常小于0.5g(g为重力加速度)，根据这个特点，当采样间隔为T秒时，选取中值滤波的阈值为0.5gT。中值滤波过程如下：对测得的右后轮轮速数据V_rr(i)(i＝1，2，…，n)进行中值滤波，即保持V_rr(1)、V_rr(2)、V_rr(n-1)和V_rr(n)不变，如果V_rr(j)≥V_rr(j-1)+0.5gT或者V_rr(j)≤V_rr(j-1)-05gT，则V_rr(j)取median{V_rr(j-2)，V_rr(j-1)，V_rr(j)，V_rr(j+1)，V_rr(j+2)}，否则V_rr(j)(3≤j≤n-2)不变，其中median表示取中位数；

再对测得的左后轮轮速数据V_lr(i)进行中值滤波，i＝1，2，…，n，保持V_lr(1)、V_lr(2)、V_lr(n-1)和V_lr(n)不变，取阈值为0.5gT，如果V_lr(j)≥V_lr(j-1)+0.5gT或者V_lr(j)≤V_lr(j-1)-0.5gT，则V_lr(j)取median{V_lr(j-2)，V_lr(j-1)，V_lr(j)，V_lr(j+1)，V_lr(j+2)}，否则V_lr(j)不变，3≤j≤n-2。

由上述中值滤波后的轮速，根据公式(3)可以计算得到未滤波的横摆角速度ω(i)(i＝1，2，…，n)。中值滤波只能滤除粗大误差，对于随机噪声等干扰还无法滤除。目前，信号的滤波方法主要有傅立叶变换对信号进行时频变换方法，卡尔曼滤波和神经网络滤波等。由于干扰信号随机性较强，频率成分比较复杂，当信号和噪声的频带相互重叠时，傅立叶变换对信号进行时频变换的方法无法在频域对其进行有效的区分，不能很好的滤除干扰信号的影响。卡尔曼滤波需要建立准确的系统状态和观测模型，并且需要知道精确的系统噪声和观测噪声，在实际操作中这些很难准确得到，导致卡尔曼滤波的精度不高。神经网络滤波需要较多的先验数据，并且需要较长的训练时间，这些制约了神经网络滤波的应用。

为了克服这些滤波方法的局限性，针对操纵稳定性试验的需求，本发明提出一种基于小波变换的滤波方法。小波变换属于时频分析的一种，它具有多分辨率的特点，而且在时域和频域都具有表征信号局部特征的能力，是一种具有灵活多变的时间-频率窗的时频局部化分析方法。由于其具有低熵性、多分辨率、去相关性和小波基选择的多样性等特点，使它在信号滤波领域得到了广泛的研究并取得了非常好的应用效果。用小波变换滤波的基本方法有：模极大值法，基于各尺度下小波系数的相关性法，平移不变量小波法和小波收缩阈值法等几种方法。这几种方法的特点如下：

(1)模极大值法。主要适用于信号中混有白噪声且含有较多奇异点的情况。该方法在滤波的同时，能有效地保留信号的奇异点信息，并且滤波后的信号没有多余的震荡。但是该方法的计算速度较慢，并且小波分解的尺度的选择非常重要，不好把握。

(2)基于各尺度下小波系数的相关性法。这种方法的的去噪效果比较稳定，比较适合高信噪比信号。不足之处是计算量大，并且需要估计噪声方差，主要应用在边缘检测和图像增强等方面。

(3)平移不变量小波法。该方法主要适用于信号中含有白噪声且含有若干个不连续点的情况，与阈值法相比，它可以得到更小的均方根误差，而且信噪比也能得到一定程度的提高，但其缺点是计算速度很慢。

(4)小波收缩阈值法。这种方法的优点是噪声几乎完全可以得到抑制，且反应原始信号的特征尖峰点得到很好的保留。该方法具有广泛的适应性，因而是众多小波滤波方法中最为广泛应用的一种。

针对横摆角速度滤波测量的实际要求，结合上述各种方法的特点，本发明选用小波收缩阈值法来对轮速推算得到的横摆角速度进行滤波。小波收缩阈值法的主要理论依据是，小波变换特别是正交小波变换具有很强的去数据相关性，它能使信号的能量在小波域集中在一些大的小波系数中；而噪声的能量却分布于整个小波域内，因此，经小波变换后，信号的小波系数幅值要大于噪声的系数幅值，可以认为，幅值比较大的小波系数一般以信号为主，而幅值比较小的系数在很大程度上是噪声。因此，采用阈值的办法可以把有用信号系数保留，而使大部分噪声系数减小至零。

小波收缩阈值法的滤波效果主要依赖于所选的小波基函数。选用怎样的小波基函数可由待处理的数据的特性来决定，理想的小波基函数应该具有下列性质：线性相位性，可以减少或消除重构信号在边缘处的失真；紧支集特性，支集越短，小波变换计算的复杂度越低，便于快速实现；消失矩特性，决定小波变换后能量集于低频分量的程度，消失矩越高越好。

常用的小波基函数有‘Haar’小波，‘biorNr·Nd’(Nr＝1，2，…，6，Nd＝1，2，…，8)小波系，‘coifN’(N＝1，2，3，4，5)小波系，‘dbN’(N＝1，2，…，10)小波系和‘symN’(N＝2，3，…，8)小波系。‘Haar’小波是最早用到的一个具有紧支撑集的正交小波函数，但是它是非连续的，不适合用来逼近光滑函数；‘biorNr·Nd’小波系的最主要特点是具有线性相位性，它主要应用在信号和图像的重构中；‘coifN’小波系中，相同阶数的小波函数比其他系列小波函数有较长的的支撑长度，导致其小波变换的复杂度高；‘dbN’(N＝1，2，…，10)小波系是由Daubechies提出的紧支撑标准正交小波，其大多数不具有对称性，但对于给定的支撑长度和消失矩比相同阶数的其他小波好；‘symN’小波系是对‘dbN’小波系的一种改进，与‘dbN’小波系相比，‘symN’小波系有更好的对称性，可以减小计算量。针对各种小波基的特点，结合横摆角速度的滤波测量要求，本发明选用‘sym8’小波基。为能将低频信号和高频的噪声信号有效地分离开，并为了保证横摆角速度信号的重要奇异点不会丢失，小波分解的层数选为7层。

根据小波收缩阈值法，本发明中的横摆角速度滤波方法由以下三步构成：首先，选用‘sym8’小波基，对未滤波的横摆角速度ω(i)(i＝1，2，…，n)进行7层小波分解，得到第7层低频系数和第1层到第7层高频系数；其次，对高频系数进行非线性阈值量化处理，即保留第7层低频系数不变，选择一个阈值，对从第1层到第7层的每一层高频系数进行阈值量化处理；最后，根据小波分解的第7层的低频系数和经过量化处理后的第1层到第7层的高频系数，进行小波重构，得到小波滤波后的横摆角速度。实现步骤如下：

第一步，选用‘sym8’小波基，对横摆角速度ω(i)(i＝1，2，…，n)进行7层小波分解。对信号进行小波分解与重构，可以通过Mallat算法实现。Mallat算法是信号的塔式多分辨分解与重构的著名算法，类似于快速傅里叶变换(FFT)在经典傅里叶变换中的地位。根据Mallat算法，本发明对横摆角速度信号进行小波分解的基本步骤为：用初始序列(k₀＝1，2，…，m₀，m₀＝n，n为数据的长度)序列表示未滤波的横摆角速度ω(i)(i＝1，2，…，n)序列，可将

分解为第一层低频系数

(k₁＝1，2，…，m₁，

)和第一层高频系数

(k₁＝1，2，…，m₁，

)；同理，对

(k₁＝1，2，…，m₁)可进一步分解为第二层低频系数(k₂＝1，2，…，m₂，

)和第二层高频系数

(k₂＝1，2，…，m₂，

)，依次进行7次，即可以实现对

的7层小波分解。floor表示下取整运算，即若x表示任一实数，则floor(x)表示不超过x的整数中最大的一个。迭代分解过程如图3所示。

上述小波分解的具体实现过程如下(如图4所示)：首先用初始序列

(k₀＝1，2，…，m₀，m₀＝n，n为数据的长度)序列表示未滤波的横摆角速度ω(i)(i＝1，2，…，n)序列，对

(k₀＝1，2，…，m₀，m₀＝n，n为数据的长度)进行长度为16的对称延拓，得到

k′₀＝1，2，…，n+32)，将延拓后的

与低通分解滤波器LO_D作卷积，然后将卷积后的结果首、尾部各去掉16个数，再对剩余的元素进行二元下抽样即可得到第一层低频系数(k₁＝1，2，…，m₁，)；将延拓后的与高通分解滤波器HI_D做卷积，然后再将卷积后的结果首、尾部各去掉16个数，再对剩余的元素进行二元下抽样即可得到第一层高频系数

(k₁＝1，2，…，m₁，

)。同理，对

(k₁＝1，2，…，m₁，

)进行长度为16的对称延拓得到

(k′₁＝1，2，…，m₁+32)，将延拓后的

与低通分解滤波器LO_D做卷积，然后再将卷积后的结果首、尾部各去掉16个数，再对剩余的元素进行二元下抽样即可得到第二层低频系数

(k₂＝1，2，…，m₂，)；将延拓后的与高通分解滤波器HI_D作卷积，然后再将卷积后的结果首、尾部各去掉16个数，再对剩余的元素进行二元下抽样即可得到第二层高频系数

(k₂＝1，2，…，m₂，

)。按上述过程继续进行分解，依次得到第三层高频系数

(k₃＝1，2，…，m₃，

)，第四层高频系数

(k₄＝1，2，…，m₄，)，第五层高频系数(k₅＝1，2，…，m₅，)，第六层高频系数

(k₆＝1，2，…，m₆，

)，第七层高频系数

(k₇＝1，2，…，m₇，

)和第七层低频系数

(k₇＝1，2，…，m₇，

)。对称延拓是指信号延拓部分的数据与信号延拓端数据对称的一种延拓方式，长度为对称延拓的元素个数，如对{x₀，x₁，x₂，x₃，x₄，x₅}进行长度为3的对称延拓后得到{x₂，x₁，x₀，x₀，x₁，x₂，x₃，x₄，x₅，x₅，x₄，x₃}。上述小波分解的具体过程如图4所示，图4中M＝0，1…，6，图中的‘延拓’是指长度为16的对称延拓，‘截取’是指将一组数(卷积后的结果)首、尾部各去掉16个数后取剩余的数据。其中‘↓2’代表二元下抽样，它用来从一个向量中每隔一个元素抽取一个元素组成一个新的向量，如{x₀，x₁，x₂，x₃，x₄，x₅}经二元下抽样后的向量为{x₁，x₃，x₅}。

第二步，对高频系数进行非线性阈值量化处理。获得小波分解后的第七层低频系数

(k₇＝1，2，…，m₇)，第一层高频系数

(k₁＝1，2，…，m₁)，第二层高频系数

(k₂＝1，2，…，m₂)，第三层高频系数

(k₃＝1，2，…，m₃)，第四层高频系数

(k₄＝1，2，…，m₄)，第五层高频系数

(k₅＝1，2，…，m₅)，第六层高频系数(k₆＝1，2，…，m₆)以及第七层高频系数(k₇＝1，2，…，m₇)之后，就要进行对高频系数的非线性阈值量化处理。阈值的选择对滤波的效果有着很重要的影响，通常有两种方法，即硬阈值法和软阈值法。

为保持信号的整体形状不变，保留第七层低频系数

(k₇＝1，2，…，m₇)不变。取阈值

(其中n为数据的长度，函数ln(n)表示对n取自然对数)，对每个高频系数d_k，l(k＝1，2，...，7，l＝1，2，...，m_k)，采用如下硬阈值或者软阈值方法进行处理：

硬阈值法：

软阈值法：

其中

(k＝1，2，...，7，l＝1，2，...，m_k)表示阈值量化处理过的高频系数，σ是噪声的标准方差，需要进行估计。可以看出，硬阈值法虽然在均方误差意义上优于软阈值，但是所得到的估计信号会产生附加震荡，出现震铃和伪吉布斯现象；软阈值法估计得到的小波高频系数整体连续性好，不会使估计信号产生附加震荡，本发明选择软阈值法。取

这是由Dohono在1994年提出的，是基于最大最小估计得出的最优阈值。由于噪声主要分布在最高分辨级，所以可以利用高频系数

(k₁＝1，2，…，m₁)来估计噪声标准差σ，采用Donoho提出的鲁棒性中值估计：

$\mathrm{\σ}=\mathrm{median}\left(\right|d_{1,k_1}\left|\right)/0.6745(k_1=\mathrm{1,2},...,m_1)$

经上述量化处理可以得到量化处理后的第一层高频系数

(k₁＝1，2，...，m₁)，量化处理后的第二层高频系数(k₂＝1，2，...，m₂)，量化处理后的第三层高频系数

(k₃＝1，2，...，m₃)，量化处理后的第四层高频系数(k₄＝1，2，...，m₄)，量化处理后的第五层高频系数

(k₅＝1，2，...，m₅)，量化处理后的第六层高频系数(k₆＝1，2，...，m₆)，和量化处理后的第七层高频系数

(k₇＝1，2，...，m₇)。

第三步，进行小波重构，即逆小波变换，得到滤波后的横摆角速度。根据Mallat算法，进行小波重构的基本步骤为：由第七层低频系数

(k₇＝1，2，…，m₇)和量化处理后的第七层高频系数

(k₇＝1，2，…，m₇)可以得到重构后的第六层低频系数

(k₆＝1，2，…，m₆)，由(k₆＝1，2，…，m₆)和量化处理后的第六层高频系数

(k₆＝1，2，…，m₆)可以重构得到重构后的第五层低频系数

(k₅＝1，2，…，m₅)，依次进行上述过程即可重构得到滤波后的数据(k₀＝1，2，…，m₀，m₀＝n，n为数据的长度)。迭代重构过程如图5所示。

上述重构的具体实现过程为：首先让序列(k₇＝1，2，…，m₇)，再将

(k₇＝1，2，…，m₇)和量化处理后的第七层高频系数

(k₇＝1，2，…，m₇)分别进行长度为16的对称延拓得到

(k′₇＝1，2，…，m₇+32)和(k′₇＝1，2，…，m₇+32)，再对延拓后的(k′₇＝1，2，…，m₇+32)和延拓后的

(k′₇＝1，2，…，m₇+32)分别做二元上抽样，

表示

二元上抽样的结果，

表示

二元上抽样的结果，然后将

与低通重构滤波器LO_R作卷积得到序列

(p₁＝1，2，…，(2m₇+80))，

与高通重构滤波器HI_R作卷积得到序列(p₁＝1，2，…，(2m₇+80))，并将序列

与序列中的相对应元素分别进行相加，得到新序列

(p₁＝1，2，…，(2m₇+80))。找到的第

项

如果第二步中的m₆为奇数，则在中从的前一项开始向前取

项，从

的后一项开始向后也取

项，保留

这m₆个数据按原来的先后顺序排列就得到重构后的第六层低频系数(k₆＝1，2，…，m₆)；如果第二步中的m₆为偶数，则在

中从

的前一项开始向前取

项，从

的后一项开始向后取

项，保留

这m₆个数据按原来的顺序排列就得到重构后的第六层低频系数

(k₆＝1，2，…，m₆)。同理，将

(k₆＝1，2，…，m₆)和量化处理后的第六层高频系数(k₆＝1，2，…，m₆)分别进行长度为16的对称延拓得到

(k′₆＝1，2，…，m₆+32)和(k′₆＝1，2，…，m₆+32)，再分别对延拓后的

(k′₆＝1，2，…，m₆+32)和延拓后的

(k′₆＝1，2，…，m₆+32)做二元上抽样(

表示二元上抽样的结果，

表示

二元上抽样的结果)，然后将与低通重构滤波器LO_R作卷积序列

(p₂＝1，2，…，(2m₆+80))，

与高通重构滤波器HI_R作卷积序列

(p₂＝1，2，…，(2m₆+80))，并将序列

与序列

中的相对应元素分别进行相加，得到新序列(p₂＝1，2，…，(2m₆+80))。找到

的第

项

如果第二步中的m₅为奇数，则从的前一项开始向前取

项，从

的后一项开始向后也取

项，保留

这m₅个数据按原来的顺序排列就得到重构后的第五层低频系数

(k₅＝1，2，…，m₅)；如果第二步中的m₅为偶数，则从

的前一项开始向前取项，从

的后一项向后取

项，保留

这m₅个数据按原来的顺序排列就得到重构后的第五层低频系数

(k₅＝1，2，…，m₅)。按上述过程依次继续进行重构，最终可得到

(k₀＝1，2，…，m₀，m₀＝n)，即为滤波后的横摆角速度数据。上述过程如图6所示，图6中M＝6，5，…，0，图6中的‘延拓’即为长度为16的对称延拓，‘中心截取’是指对

和得到的序列

(p_7-M＝1，2，…，(2m_M+1+80))，找到第项

如果第二步中的m_M为奇数，则从

的前一项开始向前取

项，从

的后一项开始向后也取

项，保留

这m_M个数据按原来的顺序排列就得到重构后的第M层低频系数

(k_M＝1，2，…，m_M)；如果第二步中的m_M为偶数，则从

的前一项开始向前取

项，从

的后一项开始向后取项，保留这m_M个数据按原来的顺序排列就得到重构后的第M层低频系数

(k_M＝1，2，…，m_M)。‘↑2’代表二元上抽样，它用来从一个向量中每隔一个元素填充一个0元素组成一个新的向量，如{x₀，x₂，x₄}经二元上抽样得到的向量为{0，x₀，0，x₂，0，x₄，0}。m₀＝n，即为原始数据的长度。

由于本发明选用‘sym8’小波基来实现小波分解与重构，因此上述低通分解滤波器LO_D、高通分解滤波器HI_D、低通重构滤波器LO_R和高通重构滤波器HI_R的取值分别为：

LO_D＝{-0.0034 -0.0005 0.0317 0.0076 -0.1433 -0.0613 0.48140.7772 0.3644 -0.0519 -0.0272 0.0491 0.0038 -0.0150 -0.00030.0019}；

HI_D＝{-0.0019 -0.0003 0.0150 0.0038 -0.0491 -0.0272 0.05190.3644 -0.7772 0.4814 0.0613 -0.1433 -0.0076 0.0317 0.0005-0.0034}；

LO_R＝{0.0019 -0.0003 -0.0150 0.0038 0.0491 -0.0272 -0.05190.3644 0.7772 0.4814 -0.0613 -0.1433 0.0076 0.0317 -0.0005-0.0034}；

HI_R＝{-0.0034 0.0005 0.0317 -0.0076 -0.1433 0.0613 0.4814-0.7772 0.3644 0.0519 -0.0272 -0.0491 0.0038 0.0150 -0.0003-0.0019}。

实施实例2

为检验本发明提出的横摆角速度滤波测量方法的实际效果，进行了实车实验。实验基本情况说明如下：

实验目的：检验本发明提出的车辆横摆角速度的滤波测量方法的效果。

实验系统组成：实验系统由软件程序和硬件设备共同组成。横摆角速度滤波程序是按照本发明提出的横摆角速度滤波方法编制的；主要硬件设备包括：计算机(AMD TK-53CPU、1G内存)，Buick实验用车，中星环宇ZX-VG MEMS陀螺仪，汽车轮速传感器，车载电源逆变器，PC-104工控机等。

实验设置：两个汽车轮速传感器分别装在汽车后轮上，数据采集模块通过无线方式将采集信号传到PC-104中，采样间隔T为10ms。陀螺仪固定于车顶上方的中心位置。

实验线路：在总装备部定远汽车试验场、南京江宁开发区等路面上进行了多次跑车实验。

实验结果：实验表明(见下面的滤波效果图)，本发明提出的横摆角速度滤波方法取得了良好的滤波效果，有很好的实际应用价值，具体体现在：①以陀螺仪测得的横摆角速度为参照，设未滤波的横摆角速度数据与陀螺仪测得的横摆角速度数据之差为y₁，滤波后的数据与陀螺仪数据之差为y₂，可得y₁的均值为0.0112rad/s，y₁的均方差为0.0754rad/s，y₂的均值为0.0112rad/s，y₂的均方差为0.0157rad/s，可以看出均方差大大减小了，大幅度提高了横摆角速度数据的精度；②能有效地滤除各种噪声，滤波后的数据很平滑，很接近陀螺仪测出的数据，具有较高的实际应用价值。

图7～11给出了某次实验结果图。图7是未滤波的横摆角速度、小波滤波后的横摆角速度以及陀螺仪直测的横摆角速度三种情形的对比图。为显示清楚，图8单独给出了未滤波的横摆角速度和陀螺仪直测的横摆角速度对比图，而图9给出了小波滤波后的横摆角速度和陀螺仪直测的横摆角速度对比图。为更进一步的说明滤波的效果，图10～11是图7的局部放大图。从上述实验效果图(图7～11)，可以看出未滤波的横摆角速度数据波动较大，含有较多的噪声，而滤波后的横摆角速度数据非常接近陀螺仪直测的横摆角速度数据，说明本发明提出的基于小波变换的车辆横摆角速度滤波测量方法能够有效地滤除噪声，取得了良好的滤波效果。

Claims (1)

1.一种基于小波变换的车辆横摆角速度滤波测量方法，适用于前轮为转同轮、后轮为非转向轮的汽车横摆角速度的测量，其特征在于包含如下步骤：

步骤1)首先利用汽车轮速传感器测出右后轮的轮速V_rr(i)和左后轮的轮速V_lr(i)(i＝1，2，...，n，i表示采样序数，n为数据的长度)，右后轮和左后轮均为非转向轮；

步骤2)对测得的右后轮轮速数据V_rr(i)先进行中值滤波，i＝1，2，...，n，除去粗大误差，保持V_rr(1)、V_rr(2)、V_rr(n-1)和V_rr(n)不变，设数据的采样间隔为T秒，取阈值为0.5gT，g为重力加速度，如果V_rr(j)≥V_rr(j-1)+0.5gT或者V_rr(j)≤V_rr(j-1)-0.5gT，则V_rr(j)取median{V_rr(j-2)，V_rr(j-1)，V_rr(j)，V_rr(j+1)，V_rr(j+2)}，否则V_rr(j)不变，3≤j≤n-2，median表示取中位数；

再对测得的左后轮轮速数据V_lr(i)进行中值滤波，i＝1，2，...，n，保持V_lr(1)、V_lr(2)、V_lr(n-1)和V_lr(n)不变，取阈值为0.5gT，如果V_lr(j)≥V_lr(j-1)+0.5gT或者V_lr(j)≤V_lr(j-1)-0.5gT，则V_lr(j)取median{V_lr(j-2)，V_lr(j-1)，V_lr(j)，V_lr(j+1)，V_lr(j+2)}，否则V_lr(j)不变，3≤j≤n-2；

步骤3)根据步骤2)中得到的滤波后的轮速数据，由公式得到未滤波的横摆角速度ω(i)，i＝1，2，...，n，根据Mallat算法，并选用‘sym8’小波，对横摆角速度数据ω(i)进行7层小波分解，得到第七层低频系数

k₇＝1，2，...，m₇，

以及第一层高频系数

k₁＝1，2，...，m₁，

第二层高频系数

k₂＝1，2，...，m₂，

第三层高频系数k₃＝1，2，...，m₃，

第四层高频系数

k₄＝1，2，...，m₄，

第五层高频系数

k₅＝1，2，...，m₅，

第六层高频系数k₆＝1，2，...，m₆，和第七层高频系数

k₇＝1，2，...，m₇，

floor表示下取整运算，T_w为后轴左右轮的轮距，分解的具体过程如下：首先用初始序列

表示未滤波的横摆角速度ω(i)序列，i＝1，2，...，n，k₀＝1，2，...，m₀，m₀＝n，对序列

进行长度为16的对称延拓，得到

k′₀＝1，2，...，n+32，将延拓后的

与低通分解滤波器LO_D作卷积，然后将卷积后的结果首、尾部各去掉16个数，再对剩余的元素进行二元下抽样即可得到第一层低频系数k₁＝1，2，...，m₁，

将延拓后的

与高通分解滤波器HI_D做卷积，然后再将卷积后的结果首、尾部各去掉16个数，再对剩余的元素进行二元下抽样即可得到第一层高频系数

k₁＝1，2，...，m₁，

同理，接着对

k₁＝1，2，...，m₁，进行长度为16的对称延拓得到

k′₁＝1，2，...，m₁+32，将延拓后的

与低通分解滤波器LO_D做卷积，然后再将卷积后的结果首、尾部各去掉16个数，再对剩余的元素进行二元下抽样即可得到第二层低频系数

k₂＝1，2，...，m₂，

将延拓后的

与高通分解滤波器HI_D作卷积，然后再将卷积后的结果首、尾部各去掉16个数，再对剩余的元素进行二元下抽样即可得到第二层高频系数

k₂＝1，2，...，m₂，

按上述过程继续进行分解，依次得到第三层高频系数

k₃＝1，2，...，m₃，

第四层高频系数

k₄＝1，2，...，m₄，

第五层高频系数

k₅＝1，2，...，m₅，

第六层高频系数

k₆＝1，2，...，m₆，

第七层高频系数

k₇＝1，2，...，m₇，和第七层低频系数

k₇＝1，2，...，m₇，其中低通分解滤波器LO_D和高通分解滤波器HI_D的取值如下：

LO_D＝{-0.0034 -0.0005 0.0317 0.0076 -0.1433 -0.0613 0.48140.7772 0.3644 -0.0519 -0.0272 0.0491 0.0038 -0.0150 -0.00030.0019}，

HI_D＝{-0.0019 -0.0003 0.0150 0.0038 -0.0491 -0.0272 0.0519 0.3644-0.7772 0.4814 0.0613 -0.1433 -0.0076 0.0317 0.0005 -0.0034}，

步骤4)小波分解高频系数的阈值量化：保留第七层低频系数

不变，k₇＝1，2，...，m₇，对分解得到的高频系数序列d_k，l用软阈值法进行量化处理，k＝1，2，...，7，l＝1，2...，m_k，即当d_k，l≥λ时，

当d_k，l≤-λ时，

当|d_k，l|＜λ时，

其中，k＝1，2，...，7，l＝1，2...，m_k，

k₁＝1，2，...，m₁；经量化处理得到量化处理后的第一层高频系数

k₁＝1，2，...，m₁，量化处理后的第二层高频系数

k₂＝1，2，...，m₂，量化处理后的第三层高频系数

k₃＝1，2，...，m₃，量化处理后的第四层高频系数

k₄＝1，2，...，m₄，量化处理后的第五层高频系数

k₅＝1，2，...，m₅，量化处理后的第六层高频系数

k₆＝1，2，...，m₆，和量化处理后的第七层高频系数k₇＝1，2，...，m₇；

步骤5)进行逆小波变换，重构滤波后的信号：由第七层低频系数

k₇＝1，2，...，m₇，量化处理后的第一层高频系数k₁＝1，2，...，m₁，量化处理后的第二层高频系数

k₂＝1，2，...，m₂，量化处理后的第三层高频系数

k₃＝1，2，...，m₃，量化处理后的第四层高频系数

k₄＝1，2，...，m₄，量化处理后的第五层高频系数

k₅＝1，2，...，m₅，量化处理后的第六层高频系数

k₆＝1，2，...，m₆，以及量化处理后的第七层高频系数

k₇＝1，2，...，m₇，做逆小波变换即可得到滤波后的横摆角速度信号，具体过程如下：首先让序列

k₇＝1，2，...，m₇，再将

和量化处理后的第七层高频系数

分别进行长度为16的对称延拓得到

和

其中，k₇＝1，2，...，m₇，k′₇＝1，2，...，m₇+32，再对延拓后的和延拓后的

分别做二元上抽样，k′₇＝1，2，...，m₇+32，

表示二元上抽样的结果，表示

二元上抽样的结果，然后将与低通重构滤波器LO_R作卷积得到序列

p₁＝1，2，...，(2m₇+80)，

与高通重构滤波器HI_R作卷积得到序列

并将序列

与序列

中相对应元素分别进行相加，得到新序列p₁＝1，2，...，(2m₇+80)，找到的第

项

如果步骤4)中的m₆为奇数，则在

中从

的前一项开始向前取

项，从的后一项开始向后也取

项，保留这m₆个数据按原来的先后顺序排列就得到重构后的第六层低频系数

k₆＝1，2，...，m₆；如果步骤4)中的m₆为偶数，则在

中从

的前一项开始向前取项，从的后一项开始向后取项，保留这m₆个数据按原来的顺序排列就得到重构后的第六层低频系数

k₆＝1，2，...，m₆，同理，将

和量化处理后的第六层高频系数分别进行长度为16的对称延拓得到和

k₆＝1，2，...，m₆，k′₆＝1，2，...，m₆+32，再分别对延拓后的

和延拓后的

做二元上抽样，k′₆＝1，2，...，m₆+32，

表示二元上抽样的结果，

表示二

元上抽样的结果，然后将

与低通重构滤波器LO_R作卷积得到序列

p₂＝1，2，...，(2m6+80)，

与高通重构滤波器HIR作卷积得到序列

并将序列

与序列

中相对应元素分别进行相加，得到序列

p₂＝1，2，...，(2m₆+80)，找到

的第

项如果步骤4)中的m5为奇数，则从

的前一项开始向前取

项，从

的后一项开始向后也取

项，保留

这m₅个数据按原来的顺序排列就得到重构后的第五层低频系数

k₅＝1，2，...，m₅；如果步骤4)中的m₅为偶数，则从

的前一项开始向前取

项，从

的后一项向后取

项，保留

这m₅个数据按原来的顺序排列就得到重构后的第五层低频系数

k₅＝1，2，...，m₅，按上述过程继续进行重构，最终可得到

k₀＝1，2，...，m₀，m₀＝n，即为滤波后的横摆角速度数据，其中低通重构滤波器LO_R和高通重构滤波器HI_R的取值如下：

LO_R＝{0.0019 -0.0003 -0.0150 0.0038 0.0491 -0.0272 -0.0519 0.36440.7772 0.4814 -0.0613 -0.1433 0.0076 0.0317 -0.0005 -0.0034}，

HI_R＝{-0.0034 0.0005 0.0317 -0.0076 -0.1433 0.0613 0.4814 -0.77720.3644 0.0519 -0.0272 -0.0491 0.0038 0.0150 -0.0003 -0.0019}。

Images