锂离子电池的SOC荷电状态估算是BMS里最有技术含量也最让人头疼的一块。安时积分会累积误差开路电压法又没法在线用工程上主流的精度方案基本都落在卡尔曼滤波这一族。我之前在纯C语言环境里实现过一套锂电池SOC仿真模型把扩展卡尔曼滤波EKF和容积卡尔曼滤波CKF完整各跑了一遍并且分别在固定参数和带遗忘因子的在线参数识别FFRLS两种配置下做了交叉验证从模型搭建、公式推导、C代码实现到调参踩坑整个过程都亲自走了一遍。这篇文章就把这些经验完整整理出来适合正在做BMS算法验证、准备把卡尔曼滤波从仿真环境真正搬到嵌入式C代码里的工程师参考。1. 为什么把这套模型做成“双滤波双参数模式”1.1 选型逻辑EKF是行业基准线CKF是更好的非线性工具先聊一个最基础的问题为什么是EKF和CKF而不是更主流的UKF或者号称吊打一切的粒子滤波EKF的思路是把非线性函数做一阶泰勒展开用雅可比矩阵代替真实导数然后套用标准卡尔曼滤波的框架。它的优势很直接——计算量小、公式简单、几十年工程验证几乎所有BMS相关的文献都拿它当基准。我项目里的EKF部分目的就是给整套模型立一个“最低可用”的参照系任何新算法如果连EKF都跑不过那就不值得继续投入。CKF是2009年前后出现的球形-径向容积积分方法核心思想是用2n个等权重的容积点去逼近非线性函数的概率分布n是状态维数。和UKF最大的区别在于UKF的sigma点需要自己选alpha、beta、kappa三个系数选不好效果差别很大而CKF的容积点和权重是固定推导出来的——每个点的权重都是1/2n没有任何需要拍脑袋的自由参数。在状态维数低比如我们这里n3时就是6个点计算量完全可以接受。至于粒子滤波需要几百上千个粒子才能保证统计精度在MCU级别的平台上实时跑SOC估算根本不现实。这个项目是仿真模型不假但目标从一开始就是“能在嵌入式硬件上复现”所以一切选型都要兼顾算力预算。1.2 “定参”和“遗忘因子”分别解决什么问题标题里有个“定参与FF”完整点说就是固定参数定参与遗忘因子递推最小二乘FFRLS两种工作模式。固定参数模式最简单粗暴先用HPPC混合脉冲功率特性测试离线辨识出电芯的欧姆内阻R0、极化电阻R1/R2、极化电容C1/C2然后把这个模型参数写死滤波器在运行过程中不再修改。好处是实时性好、代码逻辑简单坏处是电池会老化、温度会变化R0和极化参数会漂移参数一旦失配滤波器的精度就会明显下滑。FFRLS模式解决的就是这个漂移问题。它用带遗忘因子的递推最小二乘实时在线辨识模型参数遗忘因子λ决定“过去的数据以多快的速度被遗忘”λ越接近1估计越平稳但跟踪越慢λ越小跟踪越快但越容易受噪声干扰。后面的实跑结果里这一项对误差的影响比换EKF还是CKF本身还要大。这也是我为什么坚持把两个维度交叉起来做测试——只看单变量对比很容易得出错误的结论。2. 电池等效模型与状态方程离散化2.1 采用二阶RC等效电路模型的理由电池SOC估算的模型有很多从最简单的理想电压源模型到电化学模型都有。理想模型太粗糙用一个电阻加一个电压源根本描述不了极化效应电化学模型精度虽然高但对SOC估算来说参数过多、实现过重工程上性价比太差。我做这个项目时用的是一阶RC的改进版——二阶RC等效电路模型。二阶RC模型由一个欧姆内阻R0和两个RC网络串联组成表达式如下U_t OCV(SOC) - U_1 - U_2 - I * R_0其中U_1对应电化学极化时间常数一般在几秒到几十秒U_2对应浓度极化时间常数更长。两个RC网络能更好地拟合电池在动态工况下的电压响应。实测下来一阶RC在脉冲放电结束瞬间的电压预测误差通常有几十毫伏二阶RC能把误差压到十毫伏级别而三阶RC的提升就非常有限了反而明显增加计算量——这是我的经验值具体选几阶要看你手上的电芯数据和算力上限。状态向量取x [SOC, U_1, U_2]^T输入u为电流I放电为正观测量为端电压U_t。2.2 SOC-OCV曲线的标定与拟合细节SOC-OCV关系是整个观测方程的“心脏”。我这边用0.1C倍率做完整的充放电静置实验取充电和放电OCV的平均值来消除迟滞效应再把数据点拟合成高次多项式OCV(SOC) c0 c1SOC c2SOC^2 ... c6*SOC^6拟合阶数不是越高越好我试过9阶和11阶在数据区间内拟合误差小了但在曲线两端出现了明显的龙格振荡反而让dOCV/dSOC在SOC接近0和1时失真滤波器直接受到误导。最后停在6阶在SOC5%到95%的区间内最大拟合误差能控制在5mV以内够用了。这里有个必须提醒的坑如果你的电芯是磷酸铁锂它的OCV曲线在中段SOC 30%-70%非常平缓dOCV/dSOC可能只有几微伏每百分之一SOC而噪声却有几毫伏。这意味着观测方程在这段区间几乎“看不见”SOC滤波修正能力极弱这时候无论EKF还是CKF都很吃力因为问题出在系统可观测性上而不是滤波算法本身。对这个仿真模型来说最直接的对策是保证OCV拟合在平台区的导数足够平滑不要因为拟合振荡产生虚假的修正方向。2.3 连续方程离散化与雅可比矩阵推导二阶RC模型的连续状态方程是dSOC/dt -η * I / (3600 * Qn) dU_1/dt -U_1/(R_1C_1) I/C_1 dU_2/dt -U_2/(R_2C_2) I/C_2其中η是库仑效率充电取1放电约等于1Qn是额定容量Ah。用前向欧拉离散化采样周期Ts0.1s得到SOC[k1] SOC[k] - ηTs/(3600Qn) * I[k] U_1[k1] a1 * U_1[k] R_1*(1-a1) * I[k] U_2[k1] a2 * U_2[k] R_2*(1-a2) * I[k]其中a1 exp(-Ts/(R_1C_1))a2 exp(-Ts/(R_2C_2))。写成矩阵形式A矩阵是diag(1, a1, a2)B矩阵是[-ηTs/(3600Qn), R_1*(1-a1), R_2*(1-a2)]^T。观测方程的雅可比矩阵H也在这里一并算出来H [dOCV/dSOC, -1, -1]这个H矩阵就是EKF的核心后面调用OCV多项式求导数即可。注意H的第一个元素是随SOC实时变化的不能写成常数这是很多初版实现最容易写错的地方。3. EKF在C语言下的完整实现3.1 单观测量的一个关键简化卡尔曼增益退化成标量EKF的实现绕不开矩阵运算但这里有个很多人没注意的简化我们的观测量只有一个——端电压所以创新协方差S HPH^T R是一个1x1的矩阵本质上就是个标量。卡尔曼增益K P*H^T / S从“矩阵求逆”退化成了一次标量除法。这个简化对C语言实现影响巨大意味着我完全不需要写复杂的矩阵求逆函数一个除法就搞定。我在PC端验证时发现有些同学一上来就引入通用的矩阵求逆库不仅代码臃肿还引入了数值不稳定的风险。正确做法是先算S判断它接近0的情况意味着测量信息完全不可信再决定是否更新。3.2 五个核心公式对应的C代码骨架EKF的标准五个公式我按固定流程拆成预测和更新两步。以3x3矩阵和3维向量为最小结构定义如下typedef struct { float d[3][3]; } Mat3x3; typedef struct { float d[3]; } Vec3; void mat_mul(const Mat3x3 *a, const Mat3x3 *b, Mat3x3 *out); void mat_transpose(const Mat3x3 *a, Mat3x3 *out); void mat_axpy(const Mat3x3 *a, float alpha, const Mat3x3 *b, Mat3x3 *out); void mat_sym(Mat3x3 *a); /* 将矩阵对称化: a (a a) / 2 */预测部分的代码非常直观/* 状态预测: x_pred A*x B*I */ x[0] x[0] - ETA * Ts / (3600.0f * Qn) * I; x[1] a1 * x[1] R1 * (1.0f - a1) * I; x[2] a2 * x[2] R2 * (1.0f - a2) * I; /* 协方差预测: P_pred A*P*A Q */ A.d[0][0] 1.0f; A.d[1][1] a1; A.d[2][2] a2; /* 其余为0 */ mat_axpy(A, P, A); /* temp A*P*A */ mat_axpy(temp, 1.0f, Q, P_pred); /* P_pred temp Q */更新部分先算H和S再算增益和状态修正H[0] ocv_deriv(x[0]); /* dOCV/dSOC 多项式求导 */ H[1] -1.0f; H[2] -1.0f; /* S H*P_pred*H R, 标量 */ float S R; for (int i 0; i 3; i) for (int j 0; j 3; j) S H[i] * P_pred.d[i][j] * H[j]; /* K P_pred*H / S, 3x1向量 */ for (int i 0; i 3; i) { K[i] 0; for (int j 0; j 3; j) K[i] P_pred.d[i][j] * H[j]; K[i] / S; } /* 新息 z U_meas - (OCV(x0) - x1 - x2 - R0*I) */ float z U_meas - (ocv_lookup(x[0]) - x[1] - x[2] - R0 * I); /* 状态修正 */ for (int i 0; i 3; i) x[i] K[i] * z; /* 协方差更新: P (I - K*H)*P_pred */ /* 3x3手写展开避免调用矩阵乘的开销 */这版代码在采样周期0.1s下单步耗时在常见的ARM-M系列内核上只有几十微秒量级实时性完全无压力。3.3 P矩阵对称化一个避免发散的廉价处理EKF调得久了你会发现一个常见问题由于浮点舍入误差P矩阵在多次迭代后逐渐变得不对称甚至对角线出现负值然后滤波器就莫名其妙发散了。解决方式非常简单暴力——每次更新完P之后强制做一次对称化P.d[0][1] P.d[1][0] 0.5f * (P.d[0][1] P.d[1][0]); P.d[0][2] P.d[2][0] 0.5f * (P.d[1][2] P.d[2][1]); P.d[1][2] P.d[2][1] 0.5f * (P.d[1][2] P.d[2][1]);同时检查对角线是否为负如果出现负值就加上一个极小的正量比如1e-12强制恢复正定。这个操作几乎不增加时间开销但对长期稳定的影响非常大。我实际跑过连续1小时的动态工况不做对称化的版本在30分钟左右就开始出现协方差负值做了对称化之后全程稳定。4. CKF从原理到C语言实现4.1 容积点规则为什么它没有待调参数CKF的数学基础是球形-径向容积规则对n维高斯分布取2n个容积点每个点权重都是1/(2n)。第i个容积点的位置是sqrt(n)乘以单位矩阵的第i列正负各一个。对我们的3维状态就是6个点ξ_i sqrt(3) * [±1, 0, 0], [0, ±1, 0], [0, 0, ±1]每个点都带着相等的权重1/6经非线性函数传播后再加权求和得到非线性变换后的均值和协方差。这比起UKF的好处是干干净净不需要调参。UKF的alpha、beta、kappa选不好协方差可能为负而CKF在理论推导上就保证了权值非负数值稳定性天然更好。4.2 Cholesky分解CKF绕不开的那道坎CKF在一个关键环节和EKF不同协方差预测前需要把P分解成S*S^T的形式即Cholesky分解然后用S去生成容积点。所以必须实现一个3x3的Cholesky分解函数int chol3(const Mat3x3 *P, Mat3x3 *S) { int i, j, k; float sum; for (i 0; i 3; i) { for (j 0; j i; j) { sum P-d[i][j]; for (k 0; k j; k) sum - S-d[i][k] * S-d[j][k]; if (i j) { if (sum 0.0f) return -1; /* 非正定失败 */ S-d[i][j] sqrtf(sum); } else { S-d[i][j] sum / S-d[j][j]; } } } return 0; }这里有一个必须处理的异常分支如果P在数值上不再是正定矩阵Cholesky会失败函数返回-1。我在实现里做了降级策略——返回失败时本次循环退化为EKF的更新步骤同时在下一次迭代前强制对P做上述对称化和正定恢复。实际调试中这个降级分支用得比想象中频繁尤其是在滤波器收敛到稳态之后P的某些元素会变得极小浮点误差很容易让Cholesky失败。4.3 CKF完整流程的C代码框架CKF的预测阶段由于状态方程是线性的容积点传播的结果其实和EKF的线性预测完全一致所以代码复用EKF的预测部分即可。真正的区别在测量更新——先分解协方差生成容积点再把容积点通过非线性观测函数OCV拟合多项式传播/* 1. P_pred Cholesky分解 */ if (chol3(P_pred, S) ! 0) { /* 降级: 走EKF更新 */ return ekf_update(x, P_pred, U_meas, I); } /* 2. 生成6个容积点 x_i x_pred sqrt(3) * S的第i列 */ Vec3 xi[6]; for (int i 0; i 3; i) { xi[i] x_pred; xi[i].d[i] 1.732f * S.d[i][i]; xi[3i] x_pred; xi[3i].d[i] - 1.732f * S.d[i][i]; } /* 3. 容积点通过观测方程: z_i OCV(SOC_i) - U1_i - U2_i - R0*I */ float z[6], z_mean 0; for (int i 0; i 6; i) { z[i] ocv_lookup(xi[i].d[0]) - xi[i].d[1] - xi[i].d[2] - R0 * I; z_mean z[i] / 6.0f; } /* 4. 计算互协方差和量测协方差然后照常算增益和更新 */这段代码的注释已经说明了核心流程。可以看到CKF对非线性观测函数的处理是“用6个点去打探真实曲线”而不是像EKF那样“在一条直线上硬求斜率”——这正是它在OCV曲线曲率大的区域更有优势的根本原因。4.4 实测计算量对比EKF vs CKF我针对这个模型做了粗略的运算量统计以浮点乘加次数计不涉及具体硬件项目EKFCKF预测阶段约30次乘加约35次乘加含Cholesky更新阶段约40次乘加约90次乘加含6点传播总耗时比例1.0x1.6x-1.9xCKF大约是EKF的1.6到1.9倍但绝对量仍然很小。在1kHz的MCU主频下单步计算远小于1ms无论放在100ms还是10ms的控制周期里都绰绰有余。所以结论很明确在状态维数只有3时CKF的额外成本完全可以接受可以放心选用。5. 仿真工况设计与结果对比5.1 用什么工况验证才有说服力仿真验证最忌讳只用恒流放电来测精度那是自欺欺人。我做了一套完整流程先用HPPC工况离线辨识二阶RC模型的参数——在SOC的几个关键点分别做脉冲充放电静置足够时间后取电压响应曲线拟合出R0、R1、C1、R2、C2。然后跑两种动态工况验证一种是BBDST北京动态应力测试工况电流变化剧烈且带有明显脉冲特征对滤波器的动态响应要求很高另一种是匀速缓变的充放电循环考察长时间尺度上的误差累积。最后专门做了初值偏差测试把SOC初始估计强制设置为真实值偏移20%比如真值70%估计90%观测滤波器能否快速收敛回来。所有仿真里的SOC真值用电流积分法获得——因为仿真环境里电流是精确已知的安时积分没有误差可以作为参考基准。这比实验台架上用昂贵设备测SOC要方便得多也是仿真模型存在的最大价值。5.2 Q和R矩阵怎么调才不踩坑卡尔曼滤波里最“玄学”的一步就是调Q和R。我的经验是先物理后试错R是测量噪声方差可以从传感器数据手册和实际采集的静态电压噪声推算。比如端电压采集噪声标准差约3mV那R就取1e-5(0.003)^2如果ADC量化步进是5mV那R至少要取到2.5e-5量级。Q反映的是模型误差没有直接物理对应只能凭经验估。我惯用的起点是这组数值针对一个2.8Ah的三元锂电芯参数初值建议调节方向Q(1,1)SOC过程噪声1e-6偏大→估计跳动偏小→收不动初值偏差Q(2,2)、Q(3,3)RC电压过程噪声1e-5动态工况下适当调大R端电压测量噪声1e-5传感器噪声大就调大P0初始协方差diag(1e-2, 1e-2, 1e-2)不确信初值就调大帮助快速收敛一个实操技巧在模型里故意引入10%的初始SOC偏差如果滤波器在2分钟内还没收敛到2%误差以内大概率是Q(1,1)偏小而不是R的问题。先固定R逐步放大Q看响应找到“收敛快但不过度振荡”的临界点。实测下来初值偏差20%的情况下Q(1,1)1e-6时约需2分钟收敛Q(1,1)1e-4时十几秒就收回来了但稳态误差抖动明显增大。5.3 定参模式与FFRLS模式的结果对比下面这组数据来自我的仿真模型以2.8Ah三元锂电芯、BBDST工况运行60分钟为例。定参模式的参数在25℃下离线辨识仿真温度设为35℃人为制造模型失配配置组合MAE%RMSE%最大误差%定参 EKF0.921.152.8定参 CKF0.740.932.1FFRLS EKF0.580.721.6FFRLS CKF0.460.581.2两个结论很清晰第一在固定参数下CKF确实优于EKF但优势幅度有限约15%-20%的RMSE下降第二切换到FFRLS之后无论配哪种滤波器误差都大幅下降——FFRLS对失配参数的修正贡献远远大于滤波算法本身的选择。所以如果你的模型参数来自标称值而非实测辨识请优先考虑上在线参数识别那才是收益最大的一步。5.4 遗忘因子λ的选取经验FFRLS的遗忘因子λ是另一个关键旋钮。我扫了λ0.90到0.998之间的一组值λ0.95时参数跟踪最快R0估计能在几十秒内跟上工况变化但辨识出的参数波动很大导致滤波器的预测电压噪声明显。λ0.995时参数几乎不更新退化成定参模式。最后折中在λ0.98——既能在一分钟内捕捉到参数慢漂移又不会让辨识结果大幅抖动。具体数值跟你电流动态的剧烈程度有关动态越强λ要越小才能及时跟踪但要提醒的是λ太小时电芯参数的物理连续性会被破坏出现参数跳变反而干扰SOC滤波。另外提一句如果做FFRLS输入电流数据必须先做低通滤波。我最初直接把原始电流送进递推最小二乘结果R0辨识结果被电流的毛刺噪声污染得很厉害λ怎么调都不稳。把电流做了一阶低通之后辨识结果立刻稳定下来。6. 高频问题速查与排错实录6.1 我实际遇到过的问题清单仿真调试过程中我把遇到过的高频问题整理成了一张表方便定位现象可能原因解决办法CKF在运行中途出现NaNCholesky分解失败P失去正定性P矩阵对称化、加小量正定恢复失败时降级为EKFSOC估计跑到0%以下或100%以上R调得太小导致增益过大、修正过冲增加SOC限幅适当增大R电流突变瞬间误差尖峰明显模型中R0不准电压预测出现系统性偏差改用FFRLS在线修正R0或在电流突变时暂时增大Q长时间运行后误差缓慢变大电池发热导致参数漂移定参模型失效增加温度查表参数或者开启FFRLS滤波器收敛极慢初值偏差拉不回来Q(1,1)太小或P0初始方差太小增大Q(1,1)和P0优先保证收敛速度高压平台期估算值抖动OCV曲线平台段dOCV/dSOC过小观测修正不可信提高OCV拟合精度考虑联合安时积分控制修正权重这些坑我在复现过程中基本都踩过一遍前两个是最容易让人怀疑人生的。特别是Cholesky失败它不像普通bug会直接报错而是表现为滤波结果在某几秒突然跳到离谱值然后恢复。排查方法也简单在Cholesky函数返回值非0时打一个日志计数器看是不是周期性出现如果是说明P在每个循环里都在接近奇异需要从源头上改善正定性。6.2 从仿真代码移植到嵌入式代码的几个要点这套模型虽然是仿真我一开始就是按嵌入式风格写的所以移植比较顺利。这里分享几个移植时必须注意的点第一数据结构全部用固定大小的数组和结构体绝不使用动态内存分配。MCU上的内存碎片问题是隐性的但调试起来万分痛苦。3x3矩阵、3维向量都用静态声明。第二浮点运算注意单精度与双精度的取舍。ARM-M系列默认开启FPU后float运算速度远快于double内存也减半。但如果你的MCU没有FPU纯软件浮点的话double会让你吃大亏建议直接用float并在关键位置如Cholesky和协方差更新做定点数或Q格式的可选适配。第三OCV拟合多项式在嵌入式上建议改为查表加上线性插值。多项式求值在仿真环境里很省事但移植后每个周期都要算6次乘加功耗不划算。做一个129点的OCV表每0.5%SOC一个点二分查找加线性插值精度几乎无损失速度更快。第四务必在主循环里加一个执行时间监测。滤波函数的前后各读取一次硬件定时器算出差值连续N次超时再触发报警。嵌入式系统最怕的不是算法本身变慢而是被中断或其他任务抢占后悄悄超时这个监测能帮你早期发现调度问题。最后聊几句个人体会整套模型跑完之后我最大的感受是EKF和CKF在SOC估算上的差距并没有论坛上吹得那么玄乎。在模型参数准确的情况下两者都能把误差稳定控制在2%以内真正把精度拉开的是模型失配的场景——电池老化、温度变化、工况突变。这时候有没有遗忘因子在线修正参数远比纠结滤波算法本身更值得投入精力。如果你也正在做类似的仿真我的建议是先把OCV曲线标定做扎实再用HPPC把等效模型参数辨识干净最后才轮到滤波算法的迭代。参数识别和滤波算法的配合才是整个SOC估算链路里性价比最高的优化点。