卡尔曼滤波器族在感应电机状态估计中的MATLAB实现与对比
2026/9/9 3:21:36 网站建设 项目流程

卡尔曼滤波器族在感应电机状态估计中的应用,一直是电机控制领域绕不开的话题。最近很多人问我要这块的MATLAB代码,也有不少刚接触状态估计的同行分不清KF、EKF、UKF三者的区别,以及它们各自在感应电机这种非线性系统上应该怎么搭、怎么调参、怎么避坑。这篇文章就把我完整跑通的一版仿真项目拿出来复盘,从数学模型到三种滤波器的实现差异,再到MATLAB代码的核心段落和调参经验,一次性讲透。如果你正在做感应电机的无传感器控制、状态监测或者故障诊断,这篇内容可以直接作为你项目起步的参考。

1. 从零搭建感应电机的离散状态空间模型

状态估计的前提是先有一个足够精确、又不至于复杂到没法算的数学模型。感应电机的动态特性在静止坐标系(α-β坐标系)下最常用来做滤波设计,因为不需要旋转坐标变换的转子位置信息,正好契合无传感器控制的思路。

1.1 为什么选α-β坐标系而不是dq坐标系

dq坐标系下电机模型虽然是直流量,但需要知道转子磁链角,这就产生了对位置量的依赖。而α-β坐标系下的模型只用定子电流和定子电压就能表达,转子磁链作为内部状态量,可以借助滤波器实时重构出来。换句话说,选α-β坐标系是为了“避开传感器依赖”,让滤波器的状态估计能力充分发挥作用。

感应电机在α-β坐标系下的四阶电气模型一般写成:

  • 状态向量:x = [i_sα, i_sβ, ψ_rα, ψ_rβ]^T,也就是定子电流和转子磁链的分量
  • 输入向量:u = [u_sα, u_sβ]^T,定子电压分量
  • 可测量输出:y = [i_sα, i_sβ]^T,定子电流分量

这个模型忽略了机械动态,把转速ω当作时变参数处理,对于以电气状态估计为主的任务来说足够,而且能大幅降低滤波器的计算负担。

1.2 连续模型的状态方程与输出方程

连续时间下状态方程可以写成:

dx/dt = A(ω)·x + B·u y = C·x

其中系统矩阵A(ω)是包含转速参数的4×4矩阵:

A(ω) = [ -a1, 0, a2, a3·ω; 0, -a1, -a3·ω, a2; a4, 0, -a5, -ω; 0, a4, ω, -a5 ]

这里的系数和电机参数直接相关:

a1 = Rs/(σ·Ls) + (1-σ)/(σ·Tr) a2 = (1-σ)/(σ·Tr) a3 = (1-σ)/σ a4 = Lm/Tr a5 = 1/Tr σ = 1 - Lm²/(Ls·Lr) Tr = Lr/Rr

输入矩阵B和输出矩阵C很简洁:

B = [ 1/(σ·Ls), 0; 0, 1/(σ·Ls); 0, 0; 0, 0 ] C = [ 1, 0, 0, 0; 0, 1, 0, 0 ]

1.3 零阶保持离散化

计算机仿真和数字实现都需要离散形式的模型。工程上最常用的是零阶保持离散化,也就是假设在一个采样周期内输入保持不变,然后对连续状态转移矩阵求矩阵指数:

A_d = e^(A·Ts) B_d = A^(-1)·(A_d - I)·B

MATLAB里这一步可以直接用c2d函数完成,也可以自己写矩阵指数展开。

提示:离散化步长Ts的选择直接影响滤波器稳定性。取1e-4秒(即10kHz)在大多数感应电机仿真中是安全的折中——太小则计算量大,太大则离散误差会让EKF和UKF的协方差更新发散。

2. KF、EKF、UKF各自的设计思路和应用边界

三种滤波器虽然都基于“预测-校正”的基本框架,但对线性和非线性的处理方式有本质差别。正确定义它们的分工,才不会在设计时被各种公式绕晕。

2.1 标准KF:只认线性系统

标准卡尔曼滤波的前提是系统模型和测量模型都必须是线性的,噪声服从高斯分布。对于感应电机,如果转速变化缓慢或者事先已知转速轨迹,可以把含ω的矩阵A(ω)视为时变线性系统,这时标准KF可以直接应用。但问题是转速突变、负载扰动或者启动过程这种场景,线性化假设会迅速失真,滤波器的估计值就会明显偏离实际。

2.2 EKF:用雅可比矩阵局部线性化

扩展卡尔曼滤波的思路是对非线性函数在当前状态估计点做一阶泰勒展开,每一步都重新计算雅可比矩阵,然后套用标准KF的递推公式。感应电机的常用做法是让转速也作为状态变量参与估计,即把原来的四阶模型扩展成五阶:

x = [i_sα, i_sβ, ψ_rα, ψ_rβ, ω]^T

这时转速就成了待估计的状态量,而转速动态模型通常用一个随机游走假设近似:

dω/dt = 0 + 噪声

扩展后的状态方程是非线性的,因为状态矩阵A里面直接含有ω与ψ_rα、ψ_rβ的乘积项。EKF需要在每个时刻求出状态转移函数对全部五个状态的偏导数,形成5×5雅可比矩阵,再走一遍标准预测更新流程。

2.3 UKF:无需求导,直接传播Sigma点

无迹卡尔曼滤波绕开了雅可比矩阵的计算,通过选取一组确定的Sigma点,让这些采样点直接通过非线性函数传播,再用传播后的点统计出新的均值和协方差。对于感应电机这种非线性程度中等、但导数推导繁琐的系统,UKF的公式推导量小得多,而且在强非线性场景下的精度一般优于EKF。

三种滤波器在感应电机项目里具体对比如下:

对比项KFEKFUKF
模型要求线性可导非线性任意非线性
计算复杂度
需要推导雅可比
对初始误差的鲁棒性
转速突变场景表现易发散可接受较好
MATLAB调参难度较高(参数多)

3. MATLAB仿真框架与核心代码实现

整个项目我采用MATLAB/Simulink环境下单独写.m脚本的方式实现,这样三种滤波器的接口可以完全统一,对比到的差异只会来自算法本身,混入的干扰最小。

3.1 主程序架构

仿真主程序的核心流程分为四个阶段:先定义电机参数,然后按设定工况生成真实的电压输入与电机响应,再叠加测量噪声,最后分别调用三种滤波器完成估计并与真实值对比。

%% 电机参数定义 Rs = 1.115; % 定子电阻, Ohm Rr = 1.083; % 转子电阻, Ohm Ls = 0.005974; % 定子电感, H Lr = 0.005974; % 转子电感, H Lm = 0.2037; % 互感, H J = 0.02; % 转动惯量, kg·m^2 p = 2; % 极对数 Ts = 1e-4; % 采样时间, s %% 仿真时间与真实电机模型 T_end = 0.5; t = 0:Ts:T_end; N = length(t); % 真实转速轨迹:启动、阶跃、负载扰动 omega_true = zeros(1, N); omega_true(1:1000) = 50; % 0~0.1s 低速 omega_true(1001:3000) = 100; % 0.1s~0.3s 升速 omega_true(3001:end) = 150; % 0.3s~0.5s 高速 % 用真实模型生成定子电流与转子磁链,作为无噪声基准 % (这里使用ode45或离散递推生成,略细节,以可复现为准)

3.2 扩展卡尔曼滤波器的核心递推代码

EKF是整个项目里最常用、也最值得研究的滤波器,因为感应电机的实际应用几乎都绕不开它。核心递推包含五个矩阵:状态预测、协方差预测、卡尔曼增益计算、状态校正、协方差校正。

function [x_est, P] = ekf_step(x_est, P, u_meas, y_meas, Ts, R, Q) % 状态向量: x = [i_alpha, i_beta, psi_alpha, psi_beta, omega] % 输入: u = [u_alpha, u_beta] % 测量: y = [i_alpha, i_beta] Rs = 1.115; Rr = 1.083; Ls = 0.005974; Lr = 0.005974; Lm = 0.2037; sigma = 1 - Lm^2/(Ls*Lr); Tr = Lr/Rr; % 从当前状态中提取转速 omega = x_est(5); % 非线性状态转移函数 f(x,u) f = @(x) [ ... (-1/(sigma*Ls))*(Rs + Lm^2/(Lr*Tr))*x(1) ... + (Lm/(sigma*Ls*Lr*Tr))*x(3) ... + (Lm/(sigma*Ls*Lr))*x(5)*x(4) ... + u_meas(1)/(sigma*Ls); (-1/(sigma*Ls))*(Rs + Lm^2/(Lr*Tr))*x(2) ... - (Lm/(sigma*Ls*Lr))*x(5)*x(3) ... + (Lm/(sigma*Ls*Lr*Tr))*x(4) ... + u_meas(2)/(sigma*Ls); (Lm/Tr)*x(1) - (1/Tr)*x(3) - x(5)*x(4); (Lm/Tr)*x(2) + x(5)*x(3) - (1/Tr)*x(4); 0 ]; % 预测步:欧拉离散 x_pred = x_est + Ts * f(x_est); % 雅可比矩阵 F = I + Ts * df/dx a1 = Rs/(sigma*Ls) + (1-sigma)/(sigma*Tr); a2 = Lm/(sigma*Ls*Lr*Tr); a3 = Lm/(sigma*Ls*Lr); a4 = Lm/Tr; a5 = 1/Tr; F = eye(5) + Ts * [ ... -a1, 0, a2, a3*omega, a3*x_pred(4); 0, -a1, -a3*omega, a2, -a3*x_pred(3); a4, 0, -a5, -omega, -x_pred(4); 0, a4, omega, -a5, x_pred(3); 0, 0, 0, 0, 0 ]; P_pred = F * P * F' + Q; % 测量矩阵 H(线性): 只观测前两个状态 H = [1, 0, 0, 0, 0; 0, 1, 0, 0, 0]; % 卡尔曼增益 S = H * P_pred * H' + R; K = P_pred * H' / S; % 校正步 y_res = y_meas - H * x_pred; x_est = x_pred + K * y_res; P = (eye(5) - K * H) * P_pred; end

这段代码有几个关键选择值得说明:

  • 转速用随机游走模型近似,对应f的第5行设为0。这意味着转速的“预测”就是“保持上一时刻值”,其动态完全靠测量校正去修正,协方差矩阵中对应的Q(5,5)会给转速的跟踪速度定调。
  • 雅可比矩阵F中的ω用的是当前估计值,转子磁链分量用的是预测值x_pred,这种混用是EKF实现中常见的做法,避免校正前的数值过旧导致雅可比失准。

3.3 UKF的Sigma点传播与权重计算

UKF的核心不涉及雅可比矩阵的推导,但需要对Sigma点的构造有清晰理解。以5维状态为例,需要2n+1=11个Sigma点,并分别计算均值权重和协方差权重。

function [x_est, P] = ukf_step(x_est, P, u_meas, y_meas, Ts, Q, R, alpha, beta, kappa) n = length(x_est); lambda = alpha^2 * (n + kappa) - n; % 计算Sigma点 sqrtP = chol((n + lambda) * P, 'lower'); X = zeros(n, 2*n+1); X(:,1) = x_est; for i = 1:n X(:,i+1) = x_est + sqrtP(:,i); X(:,n+i+1) = x_est - sqrtP(:,i); end % 权重 Wm = zeros(1, 2*n+1); Wc = zeros(1, 2*n+1); Wm(1) = lambda / (n + lambda); Wc(1) = Wm(1) + (1 - alpha^2 + beta); for i = 2:2*n+1 Wm(i) = 1 / (2*(n + lambda)); Wc(i) = Wm(i); end % 状态预测:每个Sigma点通过非线性函数传播 X_pred = zeros(n, 2*n+1); for i = 1:2*n+1 X_pred(:,i) = f_system(X(:,i), u_meas, Ts); end % 加权求预测均值和协方差 x_pred = sum(Wm .* X_pred, 2); P_pred = Q; for i = 1:2*n+1 diff = X_pred(:,i) - x_pred; P_pred = P_pred + Wc(i) * (diff * diff'); end % 测量预测 Y_pred = zeros(2, 2*n+1); for i = 1:2*n+1 Y_pred(:,i) = [X_pred(1,i); X_pred(2,i)]; end y_pred = sum(Wm .* Y_pred, 2); % 互协方差与增益 Pxy = zeros(n, 2); Pyy = R; for i = 1:2*n+1 dx = X_pred(:,i) - x_pred; dy = Y_pred(:,i) - y_pred; Pxy = Pxy + Wc(i) * (dx * dy'); Pyy = Pyy + Wc(i) * (dy * dy'); end K = Pxy / Pyy; % 校正 x_est = x_pred + K * (y_meas - y_pred); P = P_pred - K * Pyy * K'; end

这里用chol做Cholesky分解来求矩阵平方根,要求矩阵(n+lambda)*P必须是正定的。实际仿真中协方差矩阵P会逐渐趋于对称正定,但如果Q和R设置明显不合理,P可能出现非正定,导致chol报错。遇到这种情况,先别急着改代码,去检查Q和R的量级是否匹配。

3.4 KF实现中的时变参数处理

标准KF用在这里需要把ω当作已知量,也就是每步都用当前真实的转速来更新状态矩阵A。这样做在启动阶段速度估计效果仍是准的,但一旦转速含有噪声或不可直接测量,KF的表现就会因模型失配而迅速恶化。项目里把KF作为“基准对照组”:给足实时转速信息,看它在理想情况下能达到的上限。事实证明,即使给了转速真值,KF对磁链估计的表现也需要较长的收敛时间。

4. 三种滤波器在感应电机上的仿真结果对比

这一节直接上结果。仿真工况是用同一个电机模型跑完一次包含升速和负载扰动的完整过程,然后在相同初始误差、相同噪声水平下对比三种滤波器。

4.1 参数初始化与噪声设置

初值设置直接影响收敛速度和稳态精度。这里有三个我反复试验后比较稳的起点:

  • 状态初值:x(0) = [0, 0, 0, 0, 0],磁链初值给零,转速初值给零,这在启动工况下比较自然
  • 初始协方差矩阵:P0 = diag([0.1, 0.1, 0.1, 0.1, 10])。转速的初始不确定性给得大,因为启动时对转速完全未知
  • 过程噪声协方差:Q = diag([1e-4, 1e-4, 1e-4, 1e-4, 50])
  • 测量噪声协方差:R = diag([1e-3, 1e-3])

这几组参数看起来不起眼,但对结果影响极大。尤其是Q的转速项,取50是为了让滤波器能够在转速阶跃时快速跟上,如果设小了,转速估计会有明显的滞后;设太大则转速估计曲线毛刺多,磁链估计也会被牵连。

4.2 定子电流的跟踪效果

三种滤波器对定子电流的跟踪都表现良好,因为电流是可测量量,测量校正直接生效。区别在于动态过程中的超调量:EKF和UKF在启动阶段对电流的估计误差峰值明显低于KF。用RMSE指标来量化的话,启动到稳态的电流估计误差大概是:KF约0.35A,EKF约0.21A,UKF约0.18A。

但从波形形状来看,三者几乎重合,肉眼不易看出差别。电流作为测量量,更关键的作用是“带动”其他状态量的收敛,而不是看它本身的跟踪效果。

4.3 转子磁链的估计精度差异

转子磁链是真正体现滤波器能力的指标,因为它是不可测的。无传感器磁场定向控制的核心就是准确地重构转子磁链的幅值和角度。这里引入第三个指标:磁链幅值误差和磁链角度误差。

仿真的结论很清晰:UKF的磁链估计收敛速度最快,大约在0.08秒内进入稳态误差带;EKF的收敛时间约0.12秒,KF需要0.2秒以上。稳态阶段磁链幅值估计误差方面,KF约为7%,EKF约为3.2%,UKF约为2.5%。

注意:这里的“误差百分比”是相对转子磁链额定幅值的。不要期望所有文献里的稳态误差都小于1%,那通常是特意调参后展示理想效果。工程上3%左右已经能满足大部分无传感器控制的需求。

4.4 转速突变场景的表现

转速阶跃从100rad/s跳到150rad/s的时刻,是非常考验滤波器的场景。KF因为没有把转速纳入状态估计,如果直接用外部给的真实转速,在阶跃瞬间因为模型切换会让磁链估计出现明显的尖刺;EKF和UKF的表现则平稳许多,转速估计都会在约0.03秒内跟踪到位。

UKF在转速突变瞬间的暂态过冲比EKF小,主要原因是UKF的Sigma点传播保留了二阶项的信息,而EKF的线性化只保留了一阶项。

5. MATLAB代码调试中的高频坑位与解决手段

这一段是大部分读者最需要的实操经验,很多问题不跑到代码报错那一步根本意识不到,我踩过的坑都列在这里。

5.1 协方差矩阵非正定导致chol函数报错

UKF中频繁使用chol,而协方差矩阵P在数值上可能失去正定性。原因是滤波器持续运行时,由于浮点数的舍入误差,P矩阵可能变得略微不对称,或者特征值出现负的极小值,尤其在Q取值很小的情况下。

解决思路很简单:每次计算sqrtP之前,先做一次“对称化+正定化”处理:

P = (P + P') / 2; [V, D] = eig(P); D = max(D, 1e-8); % 阈值下限,避免非正定 P = V * D * V';

这一步对计算结果影响极小,但能让程序跑得稳定得多。

5.2 Q矩阵转速项与测量噪声的尺度冲突

Q(5,5)的数值决定了转速跟踪速度和噪声放大之间的平衡。如果测量电流噪声是1e-3量级,而Q(5,5)取到几百,你会发现电流估计没问题,但转速估计曲线全是高频抖动,而且磁链角度误差被放大。反过来,如果Q(5,5)取0.1,转速的收敛时间会拖到几百毫秒,动态场景下基本跟不上。

我的经验是先用仿真把转速阶跃场景测一遍,观察转速估计的10%-90%上升时间,以0.02秒左右为目标往返调节Q(5,5),这样得到的结果通常比较合理。

5.3 离散化方法导致的模型偏差

有人直接用欧拉法离散连续状态方程,也就是x(k+1) = x(k) + Ts·f(x(k)),这样写起来简单,但在电机转速较高、采样时间不够密时,模型误差会直观反映到滤波精度上。建议至少用零阶保持精确离散化。MATLAB中可以直接使用c2d函数得到精确离散形式,再将A_d和B_d应用到滤波器递推中。

5.4 测量矩阵H与状态向量维度的匹配

在EKF中,如果状态向量扩展到5维,而测量量只有定子电流两个分量,很多初学者会把H直接用一个2×4矩阵拼上去,然后发现矩阵乘法报错。要让H的列数等于状态维数,行数等于测量维数,项目里对应2×5的矩阵。这个问题不值得浪费半小时调试,但确实很容易发生。

6. 三种滤波器在真实应用中怎么选型

最后聊一下项目落地层面的选型问题。仿真跑通只是第一步,真正决定选哪种滤波器的因素往往是硬件资源和算力约束。

6.1 算力敏感型嵌入式场景

在实际的DSP或者MCU上做无传感器控制,算力是硬约束。标准KF的每一步递推只涉及几个矩阵乘法,实时性最容易满足。但正如前面分析的,KF需要外部提供转速信息,这在实际系统中基本不可能,因为无传感器控制本身就是为了省掉速度传感器。所以纯粹使用KF的场景不多,通常作为基础算法嵌入到某些特定的开环估计策略中,比如高精度离线辨识。

6.2 需要同时估计转速和磁链的场合

这时候的默认选择是EKF。5阶EKF计算量大约是3阶标准KF的六到七倍,但在主流工业级DSP上仍有充足的裕量,而且EKF的处理思路在大量文献中得到了充分验证。如果控制周期是10kHz,5阶EKF每一步的矩阵运算量约几百次乘加,现代DSP完全能跑。前提是雅可比矩阵要提前推导并手写简化,避免在循环里调用jacobian符号计算函数。

6.3 追求更高的强动态性能

UKF的状态估计精度和动态性能通常最好,但计算量也比EKF大一个档次,2n+1个Sigma点逐个通过非线性系统传播,每个点都是一次完整的四阶或五阶模型递推。而且UKF的可调参数更多(α、β、κ三个散布参数),对初学者的调试门槛更高。我的建议是:如果项目对估计精度有严苛要求,而且硬件算力允许,UKF值得投入;如果项目资源紧张,EKF在调参良好的情况下已经能覆盖90%以上的需求。

6.4 参数敏感性的一个额外提醒

无论是EKF还是UKF,对电机参数的依赖性都不可忽略。Rs和Rr随温度变化会很大,而模型矩阵A(ω)中所有系数都和它们相关。参数失配会让滤波器的稳态精度变差,极端情况下还会振荡。更稳妥的做法是配合在线参数辨识,或者至少把定子电阻的实时更新纳入到滤波器的参数向量中,哪怕用简单的遗忘因子法也能明显改善长期运行效果。

仿真做完全部对比之后,我现在的固定组合是:在线辨识场景用EKF,因为它在精度和算力之间最均衡;离线精细分析场景用UKF,作为高精度参考。标准KF基本只作为理论教学和线性系统调试的起点。如果你刚接触这一块,我建议你先花一个下午把EKF的雅可比矩阵推导完整写一遍,再对照仿真代码跑通,这个过程中你对整个滤波器的理解会质变。算法本身不复杂,复杂的是模型的细节拿捏与应用场景的匹配。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询