ARTICLE DETAIL

资讯详情

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

从矩阵分析到卡尔曼滤波:RM电控状态估计实战指南

从矩阵分析到卡尔曼滤波:RM电控状态估计实战指南 做RM电控的同学十个里有九个第一次翻开卡尔曼滤波资料时脑子里是空白的。明明“预测”和“更新”这两步用嘴讲都能懂可一到公式满屏的矩阵符号直接把你劝退。标题里的“卡尔曼滤波前瞻-矩阵分析基础”其实就是干这件事的在真正碰卡尔曼滤波之前先把线性代数里那些和它强相关的矩阵分析工具啃下来。这篇文章不讲抽象的数学定理而是从RM电控里真实会碰到的状态估计问题出发把状态向量、转移矩阵、协方差矩阵、雅可比矩阵这些东西挨个拆开。适合谁看适合已经会写一点C语言、准备上卡尔曼滤波但数学没底的同学也适合被扩展卡尔曼滤波折磨过的老队员。1. 为什么卡尔曼滤波会让电控人栽在数学上1.1 一个真实场景云台瞄准怎么调都“飘”我以前调步兵云台视觉识别给出一组带噪声的目标角度。直接把角度送给云台Yaw轴跟踪瞄准线像喝醉了酒一样左右摆。换低通滤波噪声是压住了但滞后特别严重对面装甲板一加速弹丸全打在尾巴上。后来决定上卡尔曼滤波。按网上的教程写出状态方程和观测方程x_k A x_{k-1} B u_k w_kz_k H x_k v_k然后就开始头疼了。A矩阵是什么H矩阵怎么构造P、Q、R四个矩阵的初始值给多少这几个矩阵里任何一个理解不到位滤波器要么发散要么比低通还钝。这个场景我相信很多RM队员都经历过尤其是电控刚入门、视觉那边又催着要落地的赛季中期。1.2 卡尔曼的五个公式里每一步都在做矩阵运算标准卡尔曼滤波更新流程已经很固定了预测x̂_k⁻ A x̂_{k-1} B u_kP_k⁻ A P_{k-1} Aᵀ Q更新K_k P_k⁻ Hᵀ (H P_k⁻ Hᵀ R)⁻¹x̂_k x̂_k⁻ K_k (z_k - H x̂_k⁻)P_k (I - K_k H) P_k⁻如果不懂矩阵乘法、转置、逆矩阵这里面每一个符号都是天书。其实逻辑本身不复杂A把上一时刻的状态推到当前P代表对状态信任程度的不确定性Q是模型自己带的过程噪声H把状态映射到传感器观测的空间K则是“该信预测还是信测量的权重”。但矩阵不是数乘法交换律不成立求逆也不是简单除一下这才是卡壳的地方。1.3 学矩阵分析前你要先转换心态很多同学问“我是不是要把线性代数整本书学完才能看卡尔曼”我的经验是不需要。你只需要把矩阵当成“一组数按规则组织起来的容器”然后掌握乘法、转置、逆、行列式、特征值、矩阵求导这些工具就够了。卡尔曼用到的矩阵分析是应用工具箱不是数学系的抽象空间。重点是学会什么时候用哪个运算而不是背证明。等你真的用滤波解决过几个问题再回头补理论效率会高非常多。2. 矩阵分析第一课把卡尔曼的“状态”装进矩阵2.1 状态向量与状态转移矩阵从坐标变换说起RM里最常见的状态估计是云台角度。假设我们估计两个量当前角度θ和角速度ω。状态向量可以写成x [θ, ω]ᵀ忽略控制输入时匀速旋转模型的物理关系是θ_k θ_{k-1} ω_{k-1} · dtω_k ω_{k-1}把它们打包成矩阵形式就是 x_k A x_{k-1}其中A [[1, dt], [0, 1]]如果认为电机有一个角加速度控制量 a那么还需要一个控制输入项 B u_k这里 B [[0], [dt]]u a。这个例子想说明A矩阵并不是什么神秘的东西它就是把上一时刻的多路物理量线性组合到当前时刻。对RM电控来说最常用的状态量是角度、角速度、加速度、位置、速度这些。建模时只要把“下一时刻等于什么”写清楚矩阵自然就出来了。很多同学一上来就抄别人的A矩阵抄完不知道每个元素代表什么滤波效果差了也没法调。2.2 向量与矩阵的维度为什么乘法顺序不能乱普通数乘法满足交换律a乘b和b乘a结果一样。矩阵不满足。A B和B A大多数时候不但数值不同连维度都可能不匹配。卡尔曼滤波里常用的符号维度是这样的符号维度说明xn×1状态向量An×n状态转移矩阵Pn×n协方差矩阵Qn×n过程噪声协方差zm×1观测向量Hm×n观测矩阵Rm×m测量噪声协方差Kn×m卡尔曼增益为什么预测协方差是 P_k⁻ A P_{k-1} Aᵀ而不是 Aᵀ P A这里不光要维度对物理意义也要对。P是状态空间里的不确定性A从上一时刻的状态空间映射到当前时刻所以左边乘A就不要了不对矩阵乘法不交换P两边都要和A发生关系。A在左、Aᵀ在右正好保证结果仍然是n×n而且保留了状态演化对不确定性的影响。另一个重点是H矩阵。假设用编码器测角度但角速度不可直接观测那么m1H [[1, 0]]。K矩阵就是n×m也就是2×1。这个维度匹配关系写代码的时候最好用注释标出来“此处矩阵乘法左矩阵列数必须等于右矩阵行数”。在ARM上跑滤波一不留神就会数组越界轻则滤波乱跳重则HardFault。2.3 单位矩阵、转置和逆矩阵卡尔曼公式中的“还原”操作单位矩阵I是矩阵里的“1”任何矩阵乘I不变。转置Aᵀ把行和列互换很多公式里出现Aᵀ是因为要把误差从观测空间映射回状态空间。逆矩阵A⁻¹相当于矩阵里的“除法”用来解线性方程组。卡尔曼增益里的这个组合 (H P⁻ Hᵀ R) 是m×m矩阵。H P⁻ Hᵀ 把状态协方差变换到观测空间R加上测量自身的噪声整个矩阵描述了“预测测量值的不确定性”。对它求逆就能把修正量折算回状态空间。多传感器融合时R从标量变成矩阵这个逆必须用矩阵求逆算法实现。大一大二学的线性代数里2×2矩阵可以直接套公式n3的时候要写LU分解或者用现成矩阵库。RM比赛通常建议自己维护一套轻量矩阵库把所有运算封装好后期调LQR控制器也能复用。这里有一个容易忽略的点单位矩阵I的维度是n×n但在P_k (I - K_k H) P_k⁻ 里K_k H 是n×nI也必须取n×n不是默认的2×2或者3×3。维度随状态向量数量变化建议用宏定义。3. 协方差矩阵RM赛场上误差传播的数学账本3.1 方差到协方差单个量到多个量的不确定性先回顾一个变量的方差σ²它表示这个量自身波动的剧烈程度。两个量一起波动时还需要知道它们是同向还是反向这就是协方差。把所有状态量两两之间的协方差排成一个方阵就是协方差矩阵P。对角线是每个状态量自己的方差非对角线是不同状态量之间的相关性。卡尔曼滤波每一步都在更新这个矩阵因为它代表着我们对当前状态估计的“信任程度”。Q和R不是随便填的。Q大说明你觉得模型本身不靠谱于是P会被冲大滤波更愿意相信测量R大说明你觉得传感器噪声大滤波会更平滑但反应迟钝。RM里常见错误是一律把Q和R设成0.01结果滤波效果和普通低通差不多白写一堆矩阵运算。正确做法是先估一下传感器噪声方差编码器读数在静止时的抖动范围再反推RQ则结合模型误差设定比如忽略摩擦阻尼时Q可以稍微给大一点。3.2 卡尔曼滤波里的P矩阵为什么要不断更新预测步骤里P_k⁻ A P_{k-1} Aᵀ Q表示老的不确定性在状态演化中传播再加上过程噪声带来的新不确定性。更新步骤里P_k (I - K H) P_k⁻表示观测给了一部分信息后不确定性缩小了。有人为了省事把P矩阵从头到尾设成一个固定对角矩阵不更新。短期看也能工作但一旦状态突然发生大幅度变化比如云台被弹丸砸了一下固定的P会让滤波反应不过来。P就像是滤波器的“自信心账本”它一直在线你才能真正自适应。举一维弹丸测速的例子摩擦轮速度传感器偶尔丢包模型又有打滑扰动。如果P更新不及时滤波对突变的响应会非常慢。一维还好办到二维状态时P的非对角线开始反映“角度估计误差”和“角速度估计误差”之间的相关程度。比如角度偏了很可能角速度也偏了同一方向这种相关性能在后续预测中自动传递给下个时刻。这是低通滤波永远做不到的。3.3 对角化与相关性问题IMU数据融合的一个例子IMU姿态估计里经常把状态设计成[姿态角角速度加速度计漂移]。加速度计和陀螺仪都有误差它们对姿态角的影响会耦合在一起P矩阵的非对角项就会非零。如果非对角项一直很大说明系统状态之间存在强耦合不能简单地拆开单独处理。矩阵分析里的“对角化”概念可以帮助理解通过特征向量变换可以把一组耦合的变量旋转到一个新坐标系里让它们的误差解耦。实际调试中你不必真的去对角化但理解这个过程能帮你读懂论文里“白噪声化”“解耦”的说法。有个新队员把协方差矩阵P的非对角项全部强制清零理由是“我状态量看起来没关系”。结果滤波精度下降不少。除非状态量在物理上完全独立否则不要手动删除非对角项。就算你在建模时觉得无关数据里的相关性也可能来自控制耦合或者机械安装误差P矩阵本来就应该把它们体现出来。4. 矩阵求导与雅可比矩阵非线性系统里的“线性化”工具4.1 雅可比矩阵到底是什么RM里很多模型不是线性的。云台转角到像素坐标的变换里面有旋转矩阵和相机投影弹道估计要考虑空气阻力这些方程没法写成固定的A矩阵。这个时候要用扩展卡尔曼滤波EKF核心思想是把非线性函数在当前估计点附近做一阶泰勒展开展开后的导数矩阵就是雅可比矩阵。雅可比矩阵本质就是“多输入多输出函数对每个输入求偏导然后排成一个表”。假设状态是[p, q]ᵀ测量方程是两组非线性函数z₁ f₁(p, q)z₂ f₂(p, q)那么观测雅可比矩阵就是H [[∂f₁/∂p, ∂f₁/∂q], [∂f₂/∂p, ∂f₂/∂q]]举个RM里的直觉例子镜头看到的装甲板像素位置与云台转角的关系。当目标在画面中心时像素偏移和云台转角近似线性当目标在画面边缘畸变变大非线性增强。EKF每一帧都在算一个局部的线性化近似等于沿着一根弯曲的路径不断画切线每走几步就画一条新的。4.2 用雅可比矩阵做卡尔曼滤波的近似EKF的思路EKF的预测和更新与线性卡尔曼基本相同只是把A和H替换为当前估计值处函数对应的雅可比矩阵 A_j 和 H_j。因为每次估计值都变雅可比矩阵也要跟着更新。这个近似只有在当前点附近才准。如果系统强非线性比如云台在极限角度附近出现机械限位打滑EKF可能发散。这时可以考虑无迹卡尔曼UKF或者粒子滤波但计算量会上去。对RM大部分应用来说EKF足够解决云台和底盘的问题。写EKF时矩阵求导结果要在嵌入式上实时计算。解析推导太复杂时可以用数值差分近似。比如求∂f/∂x就用(f(xε)-f(x-ε)) / (2ε)。ε不能太小否则浮点截断误差会占主导也不能太大否则线性化误差太明显。实际调试时建议先用Python把解析式算出来和数值差分结果对比确认解析实现没有手滑。4.3 实战摩擦轮测速的温度漂移补偿中用得到的导数模型RM电控里有一个容易被忽视的场景超级电容放电时电压会持续下降摩擦轮电机的响应会跟着变化。如果只把摩擦轮转速当成状态电压变化对转速的影响就会被当成随机噪声滤掉。但如果你把电压也写进状态方程那“电压变化对转速的影响系数”就是一个偏导数它会出现在雅可比矩阵里。例如状态里有转速v和电压U状态方程写成v_k v_{k-1} c·U_k·dt这里c是固定系数但如果考虑摩擦阻力随温度变化c会变成v和温度的函数所以需要重新求导。链式法则是关键复合函数求雅可比等于中间雅可比矩阵相乘。矩阵乘法不满足交换律所以推导的时候每一步顺序都要小心最好在纸上画变量依赖图再对照矩阵乘法顺序。实际项目里我不建议手推特别复杂的雅可比容易错。先用MATLAB/Python写符号推导再把结果整理成C代码然后用一个真实的log数据离线跑一遍看滤波值和原始传感器值是否贴合。这样能省一个赛季的调试时间。5. 特征值与特征向量判断系统收敛性的一扇窗5.1 矩阵作用于向量时方向不变意味着什么矩阵A表示一个线性变换。如果有非零向量v满足 A v λ v那么v是A的特征向量λ是特征值。意思是这个矩阵对向量v只做长度伸缩不改变方向。在卡尔曼滤波里A矩阵的特征值能告诉我们系统开环的稳定性。比如云台匀速旋转模型A[[1,dt],[0,1]]特征值是两个1说明这个模型本身临界稳定——既不会自己收敛也不会很快发散。角度和角速度会按照模型一直推算下去真正让误差变小的是后面的测量更新也就是“反馈”。5.2 离散系统稳定性与卡尔曼滤波收敛的关系对于离散系统 x_k A x_{k-1}如果所有特征值的模都小于1系统会收敛大于1会发散等于1刚好处在边界。卡尔曼滤波的预测步骤里如果A本身不稳定P矩阵会越来越大更新步骤又会用测量修正把P拉回来。滤波能不能收敛还和系统的可观测性有关也就是H矩阵能不能通过测量把所有状态量都“看透”。RM里的一个现象滤波输出疯狂增长很多人的第一反应是Q、R没调好。其实更快的方法是检查A矩阵的写法。比如把dt写成0.01但实际控制周期是0.002意味着每个周期里模型多推了5倍的角度和角速度特征值自然会偏离真实。A矩阵写错时特征值往往已经大于1离线分析很容易抓出来。5.3 用特征值分析避免滤波器发散写代码实现卡尔曼滤波后建议先做一件事把测量更新禁用掉只保留预测步骤看P矩阵会不会发散。如果预测步骤里P越变越大说明模型本身不可观测或者A取错了。然后再启用K更新看P是否收敛到一个比较小的稳态值。当P收敛到稳态后卡尔曼增益K会趋于常数。这个稳态增益对应的是代数黎卡提方程的解。实际工程里如果主控算力吃紧可以先在线算一段时间等P收敛到阈值后切换成常数K省掉每一步的矩阵求逆。RM主控资源有限这个技巧很划算。但要注意比赛过程中车速、射速、环境温度都会变常数K并不总是最优。折中办法是每隔一段时间重新打开在线计算或者用一组特征值变化很敏感的标志量触发重新初始化。矩阵特征值还有一个用途判断滤波器的响应速度。特征值模长越接近0滤波衰减越快接近1滤波越“钝”。这和低通滤波器的截止频率在概念上很像只不过卡尔曼的“截止频率”会随噪声环境自适应变化。6. 从矩阵到卡尔曼RM电控里最常见的五个数学坑6.1 矩阵维度不匹配运行时断言崩溃C语言里没有运行时维度检查矩阵乘法维度不对会访问越界表现往往是滤波结果随机跳变或者程序卡死。解决办法是建立统一的矩阵结构体比如typedef struct { uint8_t rows; uint8_t cols; float data[8][8]; } Matrix;在矩阵乘法函数里先检查左矩阵列数是否等于右矩阵行数不相等就返回错误码。调试阶段打开断言让程序在出错时立刻停在对应行然后再定位是哪个矩阵的维度设计错了。这样能省下大量查数组越界的时间。6.2 逆矩阵不存在奇异问题卡尔曼增益里要对 S H P Hᵀ R 求逆。如果R是全零矩阵而H P Hᵀ又奇异S就不可逆程序直接算出无穷大。实际使用中R对角线要保证为正也就是每个传感器都至少有一点噪声估计这样S才是正定矩阵逆一定存在。有些同学把P0初始化成全零矩阵这也不合适。P0全零代表“我完全相信初始状态”滤波器在后续很长时间里不太敢修正收敛很慢。建议P0对角线取比较大的值比如角度方差给1角速度方差给100表示“我其实不确定初值”这样滤波器会快速往真实值靠拢。6.3 数值误差累积P矩阵非对称浮点运算经过大量迭代之后P矩阵会因为舍入误差逐渐失去对称性甚至出现对角线为负的“负方差”。负方差显然没有物理意义会导致滤波发疯。解决办法很粗暴每次更新完P之后强制做一次对称化P (P Pᵀ) / 2这个操作几乎不增加计算量但对稳定性帮助很大。矩阵求逆和特征值计算之前尤其要保证P对称。很多标准库内部会默认输入是对称阵如果你的P已经漂移得不对称结果会不可预测。6.4 初始化不当导致滤波收敛慢卡尔曼滤波需要x0和P0。x0可以用第一帧传感器读数初始化比如角度直接读编码器角速度读陀螺仪。P0代表初值的置信度。如果P0设得特别小滤波会长时间“迷信”x0实际运动起来要几十毫秒才能跟上如果P0设大一点滤波器会快速修正。RM项目上电瞬间云台角度是确定的但角速度往往未知我习惯把角速度的方差给大一些让滤波器初期更信任陀螺仪的数据来修正。6.5 忘记时间戳状态转移在不同频率下的问题A矩阵里的dt必须和实际控制周期一致。有的同学在主循环里写死dt0.002但被视觉算法拖一下实际循环变成0.003甚至0.005模型和真实物理系统就错位了。更靠谱的做法是每一步读取系统时间计算真实的dt然后动态更新A矩阵。这个坑听起来不算矩阵分析但其实是矩阵应用里最要命的。因为A矩阵一旦和真实时间不匹配特征值偏离协方差的传播也不对后面矩阵运算做得再漂亮都没用。这也是为什么卡尔曼滤波在RM赛场上“落地难”的常见原因之一。7. 自测与建议的学习路径7.1 五个小练习带你过一遍矩阵分析与其看一堆理论不如动手自测。这五个练习都是我带新人时用过的难度递增但都和卡尔曼相关给定A[[1,0.01],[0,1]]P0[[1,0],[0,4]]手算 A P0 Aᵀ写出结果。默写2×2矩阵[[a,b],[c,d]]的逆矩阵公式并说明什么时候不存在。用一组IMU静止数据计算加速度计x、y两轴的方差和协方差填进2×2协方差矩阵。给定函数 f(θ,ω) θ ω·dt 0.5·sin(θ)·dt²对θ和ω求偏导构造雅可比矩阵。求矩阵A[[0.9,0.1],[0,0.95]]的两个特征值判断系统是否稳定。做完这几个练习再回头看卡尔曼五个公式你会觉得流畅很多。7.2 推荐资料与使用顺序网上经常搜到“矩阵分析史荣昌pdf”这类资源但我不建议一上来就啃教材正文。对RM电控来说先建立几何直觉比刷证明重要。可以先看3Blue1Brown的线性代数本质系列把“线性变换”“特征向量”“矩阵乘法”的几何意义搞清楚然后找一份带例子的卡尔曼滤波教程用Python跟着实现一遍最后如果想去证明收敛性、推导EKF的雅可比矩阵再去翻矩阵分析教材。按这个顺序走效率最高。7.3 我的个人经验边写代码边补数学效果最好我带新人最有效的办法不是先讲完线性代数再讲卡尔曼而是直接给一套能跑的卡尔曼代码把P和K的更新中间量打日志出来每个矩阵用表格打印。让他们观察角度误差如何随P变化增益如何根据噪声调整。矩阵运算在数值上的表现看得多了再回头补理论速度会快很多。数学重要但要带着问题学。比如“为什么K矩阵要乘一个Hᵀ”光看书可能记不牢。一旦你代码里预测值和观测值维度对不上调试到半夜之后这个问题会刻进DNA。RM电控的学习节奏本来就很紧与其按部就班啃教材不如直接从实际现象出发把矩阵分析当成工具用到哪学到哪。最后分享一个我自己的习惯每次在嵌入式上写一段矩阵运算我都会先用电脑脚本跑一遍同样的运算对比数值结果。卡尔曼滤波里的矩阵分析不是考试是你手里的扳手。把这一课补完后面再看LQR、状态观测器、IMU姿态解算你会发现到处都有它的影子。祝各位在赛场上滤波不发散调车不顺的时候少一点。
返回列表