ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

卡尔曼滤波在Simulink中的建模实现与参数整定

卡尔曼滤波在Simulink中的建模实现与参数整定 卡尔曼滤波这四个字在控制、导航、信号处理圈子里几乎是个绕不开的坎。我最早正儿八经去碰它是在做车辆纵向车速估计的时候手里只有轮速传感器和加速度计想要一个靠谱的车速信号一阶低通滤波压不住噪声直接积分又怕漂移最后把希望全押在卡尔曼滤波上。模型搭在Simulink里从建立系统方程、定状态转移矩阵到调Q和R参数、对比波形整个过程踩了不少坑也把很多教科书上讲得模模糊糊的细节彻底理顺了。这篇文章就把整套东西重新梳理一遍卡尔曼滤波到底在算什么、五个核心公式怎么理解以及在Simulink里怎么把系统模型真正搭起来跑通。无论你是刚接触卡尔曼滤波的在校学生还是要在工程里做状态估计的工程师这套思路和模型都能直接拿过去改改用。1. 卡尔曼滤波到底在解决什么问题1.1 从一阶低通滤波说起为什么它不够用很多人第一次听到滤波两个字第一反应就是低通滤波尤其是一阶RC低通公式简单、Simulink里一个Transfer Fcn模块就解决了。但实际做状态估计的时候一阶低通有两个很致命的毛病。第一个是滞后。只要是低通就必然有相位滞后输入信号一变输出总要慢半拍。我之前做车速估计的时候用一阶低通滤轮速信号正常巡航没什么问题但一到急加速或者重刹滤波后的车速明显跟不上真实车速的变化趋势这一慢后面ESP或者扭矩控制逻辑就会误判。第二个问题是它完全不利用系统模型的信息。一阶低通只根据当前测量值做平滑它不知道系统本身的运动规律也不知道控制输入比如电机扭矩、加速度对状态的影响。换句话说低通滤波是一个盲滤波器。卡尔曼滤波的思路完全不一样。它把系统的状态方程拿进来先根据上一时刻的状态和控制输入去预测当前状态再用当前测量值去修正这个预测。每一步都同时利用模型规律和观测数据所以它既能压制噪声又不会像低通那样产生明显的滞后这就是它被称为最优估计的底气所在。1.2 核心数学思想预测与校正的闭环卡尔曼滤波的数学思想可以浓缩成一句话把不确定性问题拆成两步先按物理规律往前走一步再用观测数据拉回来一步然后反复迭代。这个过程像极了你在一个陌生的城市里导航手机GPS告诉你位置但GPS有漂移你手里有地图和你的运动模型知道往前走一百米大概会到哪。你把模型预测的位置和GPS测到的位置按各自的置信度加权合并得到一个比两者都更准的位置。对应到离散时间系统里系统的状态方程是x(k) A * x(k-1) B * u(k-1) w(k-1)观测方程是z(k) H * x(k) v(k)其中w是过程噪声v是测量噪声它们都被假设为零均值的高斯白噪声。A是状态转移矩阵B是控制输入矩阵H是观测矩阵。卡尔曼滤波不直接处理原始测量信号而是处理状态估计值也就是从带噪声的观测量里把真实状态一步步递推出来。很多人一上来就被矩阵和高斯分布吓住了其实抛开这些形式核心就是预测步用模型给你一个先验估计更新步用观测给你一个后验估计而两者之间的权重由卡尔曼增益K来动态分配。噪声大的系统K就倾向于模型预测模型不可靠K就倾向于测量值。1.3 五个核心公式一次讲透卡尔曼滤波的计算流程可以压缩为五个公式这也是我建议所有人在搭Simulink模型之前手推一遍的内容。预测步x_pred A * x_est B * uP_pred A * P_est * A^T Q更新步K P_pred * H^T * (H * P_pred * H^T R)^(-1)x_est x_pred K * (z - H * x_pred)P_est (I - K * H) * P_pred每个量的含义要理清楚P是状态估计的误差协方差矩阵它衡量的是我对当前状态估计到底有多信任Q是过程噪声协方差矩阵表示模型本身的不确定性R是测量噪声协方差矩阵表示传感器测量值的不确定性K是卡尔曼增益它决定了在预测值和测量值之间各信多少。P越小说明估计越自信K就越小这时更新步给的修正量也小R越大说明测量越不可信K同样越小。换句话说卡尔曼滤波的智能全在这个K的自动调节上它不是一个固定参数而是随着每一步的P和R实时算出来的。这五个公式的输入输出关系可以用一个很直观的比喻来理解P_pred是我对预测结果的把握x_pred是我觉得现在应该在哪K是在多大程度上相信外部观测。(z - H*x_pred)叫新息也就是实际观测和我预测的差异有多大差异越大如果K也大修正量就越大。整个滤波器就这么一圈一圈滚下去直到P收敛到一个稳态值。2. Simulink建模前的状态方程与参数准备2.1 系统模型怎么定以纵向车速估计为例要在Simulink里搭卡尔曼滤波模型第一步不是画模块而是把状态方程和观测方程写清楚。拿我熟悉的车辆纵向车速估计举例假设整车可以简化为一维运动模型状态量取车速v控制输入取加速度a那么连续时间状态方程就是dv/dt a w其中w代表未建模的加速度扰动。离散化之后取采样周期Ts可以得到v(k) v(k-1) Ts * a(k-1) w(k-1)也就是说状态转移矩阵A1控制输入矩阵BTs。观测来源用轮速传感器设轮速直接测车速那么H1观测方程就是z(k) v(k) v_k这里v_k是传感器噪声。如果你要处理多维系统比如同时估计车速和加速度那么状态向量就是[v; a]状态转移矩阵就会变成2x2矩阵A [1 Ts; 0 1]控制输入矩阵B可以继续保留也可以把加速度并进状态里做随机游走。总之建模这一步的关键是把真实的物理过程离散化并且每个矩阵的维度要逐一核对Simulink里的矩阵乘法模块对这些是最敏感的维度一错直接报错。2.2 Q和R的工程整定思路Q和R的设定是卡尔曼滤波建模里最玄学、也最影响效果的环节。很多教程喜欢说根据经验试凑但实际工程里是有章法的。R代表测量噪声方差这个最容易估计。你拿传感器静置一段时间采集一组数据算一下方差基本就是R的近似值。比如轮速传感器静态采集数据噪声方差如果是0.1 (m/s)^2那R就定在0.1附近。Q代表过程噪声它没有直接测量手段通常要靠递推试验来标定。一个比较有效的初值方法是Q的数值级从R的百分之一到十分之一开始试。如果滤波响应太迟钝曲线像被焊死了一样老半天不动说明Q给得偏小过程噪声权重太低如果滤波输出跟着测量噪声剧烈跳动说明Q给得偏大模型预测权重太低。另外要留意Q并不是越大越好R也不是越小越好。二者之间需要一个平衡最终看P是否收敛、滤波后的曲线是否平顺且跟踪及时。我习惯把Q、R做成Simulink模型里的Constant参数或者MATLAB脚本里的变量每次仿真改数值不需要动模型结构只需要在workspace里重新赋值然后运行调试效率会提高很多。2.3 滤波模型架构选型纯模块、MATLAB Function还是C Function在Simulink里实现卡尔曼滤波通常有三条路。第一条路是用最基本的模块搭增益、加法器、矩阵乘法、单位延迟。这种方案的好处是每个环节都可视化特别适合教学和理解算法流程但缺点是模型会很乱尤其是多维系统连线多到眼晕而且矩阵运算模块配置起来也挺繁琐。第二条路是写成MATLAB Function或者MATLAB FunctionLevel-2把五步迭代直接写在脚本里代码直观、易于调试是目前大多数人的首选。第三条路是C Function或S-Function适合需要代码生成、嵌入到实际控制器里的场景尤其是在MBD开发流程中算法最终要生成C代码那就得提前用C语言实现。我个人建议只是学习验证用选第二条路要做产品部署直接从第三条路开始写。很多工程师习惯先在MATLAB Function里调通再翻译成C这中间反而容易因为语言差异引入新bug。倒是直接写C Function再在Simulink里做单元测试一步到位后续生成代码几乎没有迁移成本。3. Simulink系统模型的具体搭建与核心实现3.1 一维卡尔曼滤波的模块级实现流程为了把卡尔曼滤波五个公式和Simulink模块一一对应起来可以先从一维系统入手。我以楼层加速度估计为例搭建一个纯模块版本。你需要准备这些模块Constant输入真实值并叠加Random Number模拟带噪声观测Gain模块若干Add模块若干Unit Delay模块一个以及Scope用来观察波形。第一步搭建预测支路。从Unit Delay中获取上一时刻的状态估计值x_est通过Gain模块乘上A一维时就是1再加上B*u这个案例里如果没有控制输入就省略得到预测状态x_pred。同时从上一个误差协方差P_est出发通过Gain乘上A^2再加上Q得到预测协方差P_pred。第二步搭建更新支路。用预测协方差P_pred除以P_predR得到卡尔曼增益K。然后用测量值z减去预测状态x_pred得到新息乘以K再加上x_pred得到当前时刻的状态估计值x_est。这一步的输出接回Unit Delay的输入完成递推回路。同时计算P_est(1-K)*P_pred接入另一个Unit Delay更新协方差。这里有一个特别容易踩的坑如果直接把状态估计信号引到前面的计算环节而不经过Unit DelaySimulink会报代数环Algebraic Loop错误。解决办法就是在反馈回路上放一个单位延迟模块让它把上一拍的值保存下来下一拍再用。这也是卡尔曼滤波天然适合数字递推的原因所在每个采样周期只依赖上一拍的状态。3.2 多维状态与数组读取Selector、Demux怎么选当状态量从一维变成多维问题就来了状态向量是一个数组或者向量信号你要在回路里取某个分量用于观测更新或者要拆开送入不同的计算支路这时就需要用到信号选择模块。Simulink里处理这类问题最常见的模块是Demux和Selector。Demux会把向量信号拆成多个标量输出简单直接但它对输入端口的维度是静态匹配的状态向量维度一变整个模型都要改灵活性很差。Selector则要灵活得多它可以从一个向量或矩阵信号里按索引抽取需要的行、列或者子块。以三轴姿态估计为例状态向量是[横滚角; 俯仰角; 横滚角速率]如果你只需要第一个分量和第三分量去做观测更新用Selector的Index Vector模式设置Index为[1 3]输出就是一个二维向量正好接给后续的矩阵运算。使用Selector的时候要注意它的Port设置如果你要从多维数组里抽数据需要把Index Mode设为Index Vector或者Starting and Ending Indices并正确设置输入维度。还有一个习惯性技巧在总线信号里传卡尔曼滤波的状态估计配合Bus Selector取总线里的各个分量模型会清晰很多避免了一大堆连线绕在一起。另外从Workspace读取测试数据时我建议直接用From Workspace模块加上数组信号数据格式用Structure或者Timeseries都行。如果你在信号线上看到维度标识是一个[3x1]的向量但接收端需要[1x3]别慌加一个Squeeze模块或者用Reshape模块调整维度即可。这类维度问题在Simulink里极其常见尤其是卡尔曼滤波这种大量矩阵运算的模型排查维度错误的时候最好的办法是在关键线路上右键调出Signal Dimensions显示一眼就能看出信号是几维的比猜快得多。3.3 用C Function实现卡尔曼滤波核心算法在实际工程中我强烈推荐在Simulink模型里嵌入C Function来实现卡尔曼滤波。这样做的直接好处是单元测试很干净代码生成很顺滑而且最终的C代码可以直接给嵌入式工程师用不需要二次移植。下面是一维卡尔曼滤波C Function的参考实现输入是测量值z、上一个状态x_pre、上一个协方差p_pre输出是更新后的x_est和p_estvoid kalman_1d_step(double z, double q, double r, double *x_pre, double *p_pre, double *x_est, double *p_est) { double x_pred *x_pre; double p_pred *p_pre q; double k p_pred / (p_pred r); *x_est x_pred k * (z - x_pred); *p_est (1.0 - k) * p_pred; *x_pre *x_est; *p_pre *p_est; }在Simulink的C Function块里你需要定义输入端口z、q、r和输出端口x_est、p_est同时把x_pre和p_pre作为持久化的内部状态。一个比较常见的做法是使用DWork向量来保存这两个值在每个采样周期调用函数时读入旧值、更新后写回。实现过程中有两个细节容易被忽略一是要确认采样周期的设置C Function块必须在离散采样模式下运行否则系统会认为它连续执行导致递推关系错乱二是代码中要避免在函数内部printf或写文件这类语句在普通仿真没问题但一旦做代码生成或静态代码检查很可能报错或直接被优化掉。如果你的模型要导入外部C头文件比如已经有现成的卡尔曼滤波库或者惯性导航算法库可以在C Function块的Custom Code选项卡里配置头文件路径和源文件路径。配置之后模型里直接调用库函数即可。这一步在MBD基于模型设计流程里非常常见我建议所有做算法开发的人都尽早熟悉C Function的配置方式因为它能让你基于Simulink做完整的单元测试而不是每次都要等整套系统联调完再去查问题。3.4 仿真对比卡尔曼滤波与一阶低通滤波模型搭完之后一定要做对比验证。我建议在Simulink里建一个对比模型信号源输出一个斜坡加正弦的复合真实信号叠加白噪声后分别进入一阶低通滤波模块和卡尔曼滤波子系统然后用Scope同时观察三条曲线。这里的关键是给两个滤波器设置合理的参数一阶低通的时间常数和卡尔曼滤波的Q/R要做到大致同一水平否则对比不公平。从波形上你通常会发现两个现象。第一卡尔曼滤波输出的滞后明显小于一阶低通。这是因为它有模型预测这个前馈机制不是一味地平滑。第二卡尔曼滤波的稳态噪声抑制能力也更好。一阶低通在滞后和噪声抑制之间是此消彼长的关系时间常数调小了噪声大调大了滞后明显卡尔曼滤波通过动态调整增益能同时兼顾两者。如果对比波形不理想不要急着怀疑卡尔曼滤波本身先检查Q和R的比例。我曾经遇到过目标跟踪场景里卡尔曼滤波输出几乎跟随噪声、没有滤波效果的情况后来发现是把R设得太小了测量值被滤波器当成了高度可信自然不肯滤波。把R调大一个数量级波形立刻恢复正常。这是所有调参者都要经历的一课。4. 工程落地中的常见问题与调试技巧4.1 模型发散不收敛先查这三件事卡尔曼滤波模型最常见的故障就是输出直接飞了无穷大或者NaN。遇到这种问题我建议按顺序排查三个地方。第一检查状态转移矩阵A是不是写错了。尤其是多维矩阵很容易把转置搞混。在Simulink里可以用Display模块直接看矩阵乘法的输出如果某个元素出现NaN基本就是矩阵维度不匹配或者相乘出错。第二检查P矩阵的递推是否保持在正定范围。由于浮点数舍入误差P可能在长时间运行后退化为非正定矩阵导致K计算异常。解决办法是定期把P强制修正为对称矩阵比如P(PP^T)/2或者使用平方根卡尔曼滤波的变体来提升数值稳定性。第三检查采样时间是否有混叠。卡尔曼滤波是离散递推算法如果A矩阵里包含Ts而模型步长设置和Ts不一致那预测步的物理意义就全错了。我见过有人把Ts写成了模型步长结果模型跑得越快滤波效果越差调了半天才发现是时间单位不一致。检查代数环也是一个重点。纯模块搭建时如果反馈路径上没有延迟单元Simulink诊断窗口会提示存在代数环这种模型即使能仿真结果也可能不符合预期。看到Algebraic Loop的告警不要犹豫直接把对应的反馈路径上加上Unit Delay或者考虑改用函数实现。4.2 从仿真到硬件外部模式、C代码生成与FMU导出卡尔曼滤波算法在Simulink里调通之后下一步就是往工程落地走。这里有几个关键的Simulink功能很多人不知道或者没用好。外部模式External Mode是一个非常实用的功能。它允许你在一台主机上运行Simulink模型同时通过串口或者以太网连接目标硬件实时调整卡尔曼滤波的R、Q等参数不用重新编译就能观测波形变化。我在调试实际传感器数据时经常先用External Mode跑一遍在线把噪声方差测出来再填进模型参数里效率非常高。不过要注意外部模式对硬件的实时性要求比较高如果模型本身计算量太大或者通信链路有延迟波形会出现毛刺这时候可以尝试降低通信采样频率。代码生成方面如果最终要部署到MCU或者嵌入式Linux环境建议用Embedded Coder生成C代码。在模型配置里选好目标硬件设置好求解器为离散定步长然后生成代码。生成的代码里要特别关注卡尔曼滤波函数是否被内联优化掉了有时编译器会把看似无用的调试变量删掉但如果你把这些变量标记为Output可避免这个问题。另外建议做一次静态代码检查比如用Polyspace或者MISRA C检查重点排查矩阵运算时数组越界的问题卡尔曼滤波算法因为用了大量一维数组当矩阵用索引稍微写错编译期不会报错运行期却会踩内存这种问题越早发现越省事。再一个值得提的功能是FMU导出。Simulink模型可以导出为FMUFunctional Mock-up Unit这样就能和其他工具链做联合仿真比如跟Carsim联合、跟Amesim联合或者导入到别的仿真环境。卡尔曼滤波模块如果封装成FMU最大的好处是接口固定、模型保密可以在不泄露算法细节的前提下让别人集成。导出FMU的时候需要在模型里设置好输入输出端口并选择支持FMI 2.0的导出选项整个过程并不复杂。4.3 扩展卡尔曼滤波与惯性导航的衔接卡尔曼滤波只适用于线性系统但工程里大多数系统是非线性的比如惯性导航中的姿态解算、四旋翼的飞行控制、车辆的运动学模型。这时候就要用扩展卡尔曼滤波EKF它的核心思想是在每一步把非线性系统在当前状态附近做一阶泰勒展开得到雅可比矩阵然后用这个线性化后的模型继续套用卡尔曼滤波的五步迭代。在惯性导航场景里状态量通常包括位置、速度、姿态角、陀螺仪零偏和加速度计零偏维度经常达到15维以上。这时候Q矩阵和R矩阵的维度也会膨胀到15x15手动设定已经不太现实工程上通常用Allan方差分析陀螺和加速度计的噪声特性再确定Q的对角元素。Simulink里搭建EKF没有新的魔法依然是用MATLAB Function写雅可比矩阵计算然后复用卡尔曼滤波的更新链路。一个很常见的误区是认为EKF的雅可比矩阵只要初始给得好就行实际上雅可比需要在每一步重新计算因为它是状态依赖的。如果你在Simulink里偷懒把雅可比设为常量那模型可能在初始状态附近表现正常一旦状态跑远滤波输出就会迅速发散。这也是很多人在做四旋翼姿态解算时在地面调试一切正常一上天就出问题的原因之一。4.4 联合仿真场景中卡尔曼滤波的位置很多朋友会在Carsim和Simulink联合仿真里用到卡尔曼滤波。典型的场景是Carsim输出车辆动力学状态比如纵向速度、横摆角速度但由于传感器噪声和模型误差直接拿到控制算法里可能不够可靠。这时候卡尔曼滤波就夹在Carsim输出的真值噪声和下游控制算法之间充当一个状态估计器。需要注意的一点是Carsim和Simulink联合仿真的采样步长往往不一样Carsim通常是连续动力学仿真Simulink控制算法是离散的。卡尔曼滤波的步长必须和控制算法的离散步长一致而不能直接采用Carsim的输出步长否则滤波器递推的频率与状态方程离散化频率对不上效果会大打折扣。另外CARSim输出的量单位通常已经换算成国际单位制但有些版本会用km/h这在搭建系统模型时特别容易忽略。卡尔曼滤波对量纲极其敏感如果你把km/h当m/s代入状态方程状态转移矩阵里的系数就会差2.78倍整个模型看起来滤波有效实际输出却是错的。我的经验是每次在联合仿真里接传感器信号先做一次单位换算并验证量级再进入卡尔曼滤波器。4.5 常见问题速查表问题现象可能原因排查方向滤波输出发散为NaNP矩阵失去正定性或A矩阵错误检查矩阵维度、强制P对称化滤波曲线滞后严重Q设置太小模型预测权重过低增大Q值或检查离散化系数滤波曲线跟随噪声毛刺多R设置太小测量值权重过高增大R值或重新计算测量噪声方差模型报代数环错误反馈回路缺少单位延迟在反馈路径上插入Unit Delay仿真和代码生成结果不一致采样时间或代码生成优化选项设置不一致统一求解器为离散定步长禁止优化内联关键函数外部模式下波形卡顿通信速率低或模型计算量大降低通信采样频率简化Scope采样频率多维数组维度不匹配报错信号维度设置错误显示Signal Dimensions检查Selector索引最后说一个我一直沿用的习惯每次搭完一个卡尔曼滤波模型我都会刻意把测量信号停掉只用预测步跑一段看看状态估计在没有观测的情况下怎么变化。这个开环测试能帮你第一时间发现状态转移矩阵和离散化系数是否写对也能直观感受Q对预测步的影响。这招帮我省了无数次在闭环里反复排查的麻烦。卡尔曼滤波的坑确实不少但只要你把模型、参数、调试手段这三样东西理顺了它在Simulink里绝对是你做状态估计最趁手的工具。
返回列表