卡尔曼滤波算法在MEMS 热式风速风向传感器中的应用∗

赵 天田 润易真翔秦 明

(东南大学电子科学与工程学院，江苏 南京210097)

风速风向是气象测量中重要参数之一，能够直观地表征空气流动及气压分布的情况，因此，风速风向的测量在气象预测、航天、军事等领域均有重要意义[1]。 传统的风速风向传感器主要包括风杯式和超声式，但是此类传感器具有维护成本高、小型化困难的缺点[2]。 MEMS 技术的发展为解决上述问题提供了新的技术方案。 MEMS 热式风速风向传感器主要由加热元件和温度敏感元件组成，基于强迫对流原理，通过监测芯片表面气体流动引起的热场梯度变化实现风速和风向测量[3-5]。 1974 年，荷兰代尔夫特理工大学的Van Putten 和Middelhoek 首次采用标准的硅加工技术成功研制了热流量传感器[6]。自此之后，世界各大高校纷纷开始了MEMS 风速风向传感器的研究[7-11]，通过近半个世纪的理论和技术积累，该传感器已逐渐可应用于实际生产生活中[12-16]。 相较于传统的风速风向仪，MEMS 热式风速风向传感器具有体积小、功耗低、精度高、便于集成和可批量化制备的优点[17]。

但是，在实际的长期应用中，本征噪声的存在容易造成MEMS 风速风向传感器输出不稳定等性能退化的现象。 由于本征噪声源自环境或传感器本身存在的分子的随机运动，属于白噪声。 这部分噪声是传感器本征固有的，无法通过电路设计解决这一问题[18]。 在此之前，实验室研究的MEMS 风速风向传感器系统均采用加权滑动均值滤波算法以减小本征白噪声[18，19]，但该算法存在滤波效果差、内存消耗多等问题。 因此，为了解决上述问题，本文首次将Kalman 滤波算法引入到MEMS 风速风向传感器中，并通过仿真及实验证明通过该数字滤波算法可以有效提高响应速度，降低测量误差。

1 Kalman 滤波原理及算法

1.1 算法简介

Kalman R E 在20 世纪60 年代初提出了Kalman滤波算法[20]。 Kalman 滤波算法是一种时域算法，其将信号处理的过程视为一个线性系统，将本征噪声视为系统的一个激励，用状态方程来描述输入和输出的关系，以此形成滤波算法。 Kalman 滤波只涉及时域变量，可以对平稳的一维随机过程如风速传感器的输出进行估计。 此外，由于算法的递推性，可以在计算机上实现即时处理。 自卡尔曼滤波算法被发明以来，其已在各类工程实践如惯性导航，制导系统，全球定位系统，通信与信号处理，航空航天等领域被广泛应用[21]。

Kalman 滤波算法的核心原理是综合考虑上一时刻的状态最优估计Xk-1以及当前时刻的测量结果Zk，并根据系统噪声和测量噪声的统计特点对当前时刻的状态的最优估计Xk进行预测，具体的推导过程利用到了线性系统原理以及射影定理等，较为复杂，在此不再赘述[21]。 基于以上理论，可以给出以下状态空间模型，如式(1)所示:

在该模型中，Xk和Xk-1均是n×1 维的状态向量；φk/k-1是n×n维的状态转移矩阵；Γk/k-1是n×l维的系统噪声分配矩阵；Zk是m×1 维的量测向量；Hk是m×n维的量测矩阵；此外Wk-1和Vk分别是l×1维的系统噪声向量和m×1 维的测量噪声向量，且这两个向量均为期望为零的白噪声序列，彼此互不相关，即满足式(2):

Kalman 滤波算法主要分为以下五步，如式(3)～式(7)所示:

①计算状态一步预测:

②计算状态一步预测方差:

③计算滤波增益:

④计算状态估计:

⑤计算新状态估计均方误差:

通过以上五步循环迭代即可实现滤波。 具体的计算过程可以用图1 来表示。

图1 Kalman 滤波迭代框图[22]

可以看到对于这个系统来说，量测信息Zk是主要的输入，并且对于非时变系统，系统参数包括φk/k-1，Γk/k-1，Hk以及噪声参数Qk-1，Rk均为定值，这样在实际计算中只需要存储Xk，Pk两个变量，输入也仅有Zk，对于单片机来说内存消耗小，计算较简单。 以下将介绍本系统中的具体参数设置。

1.2 参数设置

对于本文实验中所涉及的风速风向传感器而言，其系统状态可以设置为

其中Xk表示系统状态，Tk表示当前时刻测温元件处的温度，Ak表示测温元件处温度的变化速率，因此对于下一时刻的状态预测，可以认为状态保持不变，即:

因此，可以得到状态转移矩阵:

其中Δt为采样间隔或采样周期。

本文中涉及的热式风速传感器是通过热温差的方式实现风速测量的，具体原理为，当有风流动时会使传感器表面的热场发生改变，导致传感器上下游出现温差ΔT，通过一对热电堆可以对该温差进行测量，测量是通过Seebeck 效应实现的，Seebeck 效应又称第一热电效应，是指半导体两端的温度差异而引起电压差的热电现象，对于同一种材料有简单的表达式如式(11)所示[23]:

根据以上原理，给出式(12)所示测量方法:

即有量测矩阵:

以上两式中α均为Seebeck 系数，此外在实际应用中根据测量原理不同，这里的温度T可以指当前温度，也可以指当前温度差或是其他适用于风速测量的物理参数。

考虑到利用Seebeck 效应的计算其实是一个线性过程，因此完全可以直接将热电堆的输出电压Vk和热电堆输出电压的变化速度Dk视为系统状态的组成，这样Seebeck 系数α就不会出现在表达式中，可以显著简化运算，如式(14)所示:

至于噪声参数，测量噪声V可以通过对传感器噪声测量得到，对于确定的传感器，其具有与传感器物理结构相关的本征噪声。 为了测量该部分噪声的大小，可将传感器置于风洞中，不经过数字滤波直接测量给定风速下传感器输出电压的变化范围，将其峰峰值作为测量噪声V的取值。

对于系统噪声W，在卡尔曼滤波算法中，该参数与状态一步预测结果的可信度有关，如果该值较大，则滤波结果将更接近于测量值，如果该值较小，则滤波结果将更趋近于预测值。 在实际使用中可以通过调整该参数大小以取得更好的滤波效果。

对于采样间隔Δt这个参数值的选择有几个需要考虑的地方。 首先采样过程相当于一个低通滤波器，采样间隔Δt越大，滤波器的截止频率越低。 如果增大采样间隔，很多噪声会在采样过程中被消除，但是如果采样间隔过大，就会导致信息的丢失。 其次，选取采样间隔时还应考虑单片机本身的性能，采样间隔至少要大于一次Kalman 循环所需要的计算时间，所以Δt必须根据实际需求进行权衡。

2 Kalman 滤波算法软件仿真

本节通过软件仿真验证Kalman 滤波算法的有效性，所选仿真软件为MATLAB 2016a。

2.1 滤波效果仿真

首先，利用示波器对8.0 m/s 风速下传感器的输出信号进行采集，去除有效信号提取出噪声信号，分别采用加权滑动均值滤波算法和Kalman 滤波算法，对比两种算法的滤波效果。 本次实验中共采集了50 s 的数据，采样间隔为0.000 8 s，因此共有62 500 个数据。 所得结果如图2 所示。

图2 Kalman 滤波及加权滑动均值滤波两算法滤波效果仿真结果对比

如图2 所示，原始噪声数据的峰峰值高达9.97 mV，使用加权滑动均值滤波后的峰峰值减小至1.39 mV，而使用Kalman 滤波后的峰峰值为4.72 mV。可以看到两种滤波算法都有不错的滤波效果，其中加权滑动均值滤波效果要更好一些，但是实际应用中，不能只看滤波效果，还应考虑算法的响应速度。

2.2 响应速度仿真

为了比较两种算法的响应速度，在25 s 时刻添加了一个大小为1 V 的阶跃信号，之后分别进行加权滑动均值滤波和Kalman 滤波，观察其响应速度，其结果如图3 所示。

图3 Kalman 滤波算法及加权滑动均值滤波算法响应速度仿真结果对比

此处选取首次到达目标电压的时间作为响应时间。 可以看到，加权滑动均值滤波的响应时间为7.5 s，而Kalman 滤波的响应时间仅为1.64 s，Kalman 滤波的响应速度明显占优。 实际上，根据国家标准，有效的风向风速测量系统需要3 s 内回传一组数据，即响应时间应不小于3 s。 因此，按照上述仿真的参数设置，加权滑动均值滤波的响应速度并未达到国家标准的要求，故无法实际使用。 与之形成对比，Kalman 滤波虽然滤波效果稍差，但响应速度上仍有很多余量，考虑到在数字滤波算法中，响应速度和滤波效果是一对相互矛盾的变量，因此在实际使用中可以通过调整参数，牺牲部分响应速度以加强滤波效果。

通过仿真结果可以发现Kalman 滤波虽然滤波效果逊于加权滑动均值滤波，但相比之下响应速度快很多，在实际使用中，需要针对实际指标要求判断两种滤波算法的优劣。

3 Kalman 滤波算法实验结果

3.1 实验条件

为了验证实际使用中Kalman 滤波的效果，将Kalman 滤波算法及加权滑动均值滤波算法移植到8位单片机C8051F410 上，在风速8 m/s，湿度25%，环境温度18 ℃的条件下进行测试。 实验装置及测控电路原理如图4 所示。

实验过程中，将风速传感器与支架固定后放入风洞内，当风洞工作时，风速传感器表面的热场会根据风速大小发生不同程度改变，风速信号会首先被转化为热信号，再通过Seebeck 效应将热信号转化为电压信号，之后通过PCB 上的数据采集电路(包括滤波、放大电路) 将数据传入单片机并进行Kalman 滤波、数据处理得到风速风向信息，最后通过UART 串口发送至PC 上显示，此外PCB 板上还有一块加热电路使芯片表面电路略高于环境温度，便于测量温场变化。

图5 为实际测试中所使用的MEMS 热式风速风向传感器，PCB 电路及测试用风洞的实物照片。

图4 实验装置及测控电路原理示意图

图5 MEMS 热式风速风向传感器，PCB 电路及测试用风洞的实物照片

在本次实验中，根据1.2 节中参数设置的规则，参数值设置如下:

系统噪声方差:

测量噪声方差:

该值为风速8 m/s 时传感器低频噪声(0.1 Hz～10.0 Hz)的测量值。

测量间隔:

3.2 实验结果

根据上述参数，进行实验，得到输出电压，结果如图6 所示。

图6 Kalman 滤波算法及加权滑动均值滤波算法实际滤波效果对比

首先考虑滤波效果，可以看到原始信号噪声的峰峰值高达56.2 mV，加权滑动均值滤波后的峰峰值为19.7 mV，而Kalman 滤波后的峰峰值仅为9.7 mV，相较于原始信号降低了82.7%，相较于加权滑动均值滤波的结果降低了50.8%。

为了测试两种算法的响应速度，先将传感器启动后放置于风洞中，设置风速为0，放置15 min 待其输出稳定后将风速设置为8 m/s，等待5 min 后，收集数据，结果如图7 所示。

图7 Kalman 滤波算法及加权滑动均值滤波算法实际响应速度对比

不难看出，加权滑动均值滤波的响应时间为5.76 s，而Kalman 滤波的响应时间仅为2.56 s，相较于加权滑动滤波快了55.6%，满足3 s 回传时间的国家标准。 值得注意的是，图7 中Kalman 滤波结果存在一段过冲，但考虑到过过冲电压较小，为0.04 V，满足本项目中指标要求(在该风速段中换算到电压值为0.048 V)，因此，仍取首次到达目标电压值的时间为响应时间。

之后为计算滤波算法对测量误差的影响，考虑到无法直接将电压值与风洞设置的准确值进行对比，需要先将输出电压转化为风速数据。 由于表征输出电压与风速关系的数学表达式较复杂[24]，因此在嵌入式代码中直接使用线性插值的方法将电压与风速一一对应起来，这样就可以对测量误差进行表征。

图8(a)和8(b)分别为不同风速下经过滤波、插值后所得风速结果较风洞输入值的绝对误差和相对误差，很明显，风速测量的绝对误差随风速增大而增大，相对误差则随风速增大而减小。 此外，应用Kalman 滤波算法后， 最大风速绝对误差从±0.34 m/s 缩小至±0.25 m/s，减小了27%。 绝对误差在风速为11.2 m/s 时减小最多，由±0.20 m/s 缩小为±0.08 m/s，减幅高达60%。 相对误差在风速为0.8 m/s 时减小最多，由2.7%减小为1.42%。

可以看出，在上述参数设置条件下，Kalman 滤波无论在滤波效果、响应速度和误差上都明显优于加权滑动均值滤波，而且响应时间尚有一定余量，如有必要仍可进一步强化滤波效果。 因此，Kalman 滤波算法可以有效应用于本文所涉及的MEMS 热式风速风向传感器。

图8 Kalman 滤波算法及加权滑动均值滤波算法误差对比

4 总结与展望

本论文就Kalman 滤波应用于MEMS 热式风速风向传感器的可行性进行了讨论。 首先，介绍了Kalman 滤波算法的原理，逐步分析了Kalman 滤波算法的流程，详细给出了实际应用时的参数设置方法和注意事项。 然后，利用MATLAB 对Kalman 滤波算法进行仿真，并与传统的加权滑动均值滤波算法进行对比，从仿真结果可以看出，Kalman 滤波算法在滤波效果上稍逊于加权滑动均值滤波，但响应速度有显著优势。 最后，将Kalman 滤波算法写入单片机，直接采集传感器输出信号进行滤波，观察滤波结果并与加权滑动均值滤波的结果进行对比，通过调整参数可以使Kalman 滤波的滤波效果、响应时间以及误差都优于加权滑动均值滤波，由此确定Kalman 滤波算法可以应用于MEMS 热式风速风向传感器，该算法能有效抑制测量中出现的低频白噪声，提高传感器的精确度和稳定性。