4. 扩展卡尔曼滤波（EKF）：非线性系统的线性化处理，雅可比矩阵计算

各位同学，欢迎来到第四讲。

前面我们聊了标准卡尔曼滤波，它很完美，但有个前提——系统必须是线性的。可现实世界哪有那么多线性？电池的OCV-SOC曲线是弯的，RC网络的响应也是非线性的。这时候，标准KF就束手无策了。

怎么办？

我个人的习惯是：既然解决不了非线性，那就把它变成线性。这就是扩展卡尔曼滤波（EKF）的核心思想。

4.1 为什么需要EKF？

说白了，标准KF假设状态转移和观测都是线性关系。但电池模型里，你想想看：

• 状态方程：SOC的更新本身是线性的（安时积分），但RC网络的电压响应是指数形式的。
• 观测方程：端电压 = OCV(SOC) + RC电压 + 内阻压降。OCV和SOC的关系，就是那条著名的非线性曲线。

如果你强行用标准KF，结果就是发散。我在项目初期就踩过这个坑，直接用标准KF去估计SOC，结果没跑几分钟，估计值就飞到天上去了。嗯，从那以后，我再也不敢对非线性系统用标准KF了。

4.2 泰勒展开与线性化

回忆一下高数里的泰勒展开。对于一个非线性函数 \( f(x) \)，在 \( x = \hat{x} \) 处展开：

f(x) ≈ f(ˆx) + f'(ˆx) * (x - ˆx)

我们只取一阶项，忽略高阶项。这就是线性化的本质。

应用到电池模型里：

• 状态方程线性化：对 \( f(x_k, u_k) \) 求关于状态 \( x \) 的偏导，得到矩阵 \( F_k \)。
• 观测方程线性化：对 \( h(x_k) \) 求关于状态 \( x \) 的偏导，得到矩阵 \( H_k \)。

这两个偏导矩阵，就是传说中的雅可比矩阵。

4.3 雅可比矩阵计算实战

我们以二阶RC模型为例。状态向量为：

x = [SOC, V1, V2]^T

其中V1和V2是两个RC网络的电压。

状态方程：

SOC(k+1) = SOC(k) - (η * Δt / Q) * I(k)
V1(k+1) = V1(k) * exp(-Δt / τ1) + R1 * (1 - exp(-Δt / τ1)) * I(k)
V2(k+1) = V2(k) * exp(-Δt / τ2) + R2 * (1 - exp(-Δt / τ2)) * I(k)

求状态转移矩阵 \( F_k \)（即 \( \partial f / \partial x \)）：

F_k = [1,      0,      0;
       0,  exp(-Δt/τ1), 0;
       0,      0,  exp(-Δt/τ2)]

你看，SOC对V1、V2的偏导都是0，因为SOC更新只和电流有关。V1对SOC的偏导也是0。所以F_k是一个对角矩阵。

观测方程：

Vt(k) = OCV(SOC(k)) - V1(k) - V2(k) - R0 * I(k)

求观测矩阵 \( H_k \)（即 \( \partial h / \partial x \)）：

H_k = [dOCV/dSOC,  -1,  -1]

这里最关键的就是 dOCV/dSOC。它需要你对OCV-SOC曲线求导。

4.4 EKF的完整流程

有了雅可比矩阵，EKF的流程和标准KF几乎一样：

1. 预测步：
   • 状态预测：\( \hat{x}_{k|k-1} = f(\hat{x}_{k-1}, u_k) \)（用原始非线性函数）
   • 协方差预测：\( P_{k|k-1} = F_k P_{k-1} F_k^T + Q \)（用线性化的F_k）
2. 更新步：
   • 卡尔曼增益：\( K_k = P_{k|k-1} H_k^T (H_k P_{k|k-1} H_k^T + R)^{-1} \)
   • 状态更新：\( \hat{x}_k = \hat{x}_{k|k-1} + K_k (z_k - h(\hat{x}_{k|k-1})) \)（用原始非线性h）
   • 协方差更新：\( P_k = (I - K_k H_k) P_{k|k-1} \)

注意看，预测和更新里，状态和观测的传播用的是原始非线性函数，只有协方差传播用的是线性化的雅可比矩阵。这一点很重要，千万别搞混了。

4.5 知识体系与核心逻辑

下面这张图，是我自己总结的EKF在电池状态估计中的核心逻辑，你一看就明白：

4.6 代码实现要点

在Python里实现EKF，我建议你这样组织代码：

class BatteryEKF:
    def __init__(self, dt, R0, R1, C1, R2, C2, Q_capacity):
        # 初始化参数
        self.dt = dt
        self.R0 = R0
        # ... 其他参数
        self.x = np.array([0.5, 0, 0])  # [SOC, V1, V2]
        self.P = np.eye(3) * 0.1
        
    def state_transition(self, x, I):
        """原始非线性状态转移函数"""
        SOC, V1, V2 = x
        SOC_next = SOC - (self.eta * self.dt / self.Q) * I
        V1_next = V1 * np.exp(-self.dt/self.tau1) + self.R1 * (1 - np.exp(-self.dt/self.tau1)) * I
        V2_next = V2 * np.exp(-self.dt/self.tau2) + self.R2 * (1 - np.exp(-self.dt/self.tau2)) * I
        return np.array([SOC_next, V1_next, V2_next])
    
    def compute_F(self, x, I):
        """计算状态雅可比矩阵 F_k"""
        # 对于二阶RC模型，F_k是常数矩阵（对角）
        F = np.array([[1, 0, 0],
                      [0, np.exp(-self.dt/self.tau1), 0],
                      [0, 0, np.exp(-self.dt/self.tau2)]])
        return F
    
    def compute_H(self, x):
        """计算观测雅可比矩阵 H_k"""
        SOC, V1, V2 = x
        dOCV_dSOC = self.ocv_poly_deriv(SOC)  # 多项式求导
        H = np.array([[dOCV_dSOC, -1, -1]])
        return H

4.7 避坑指南

我曾经在EKF的实现上栽过跟头，这里分享几个经验：

• 雅可比矩阵必须在当前估计点计算。我见过有人把雅可比矩阵算成常数，结果EKF直接退化成标准KF，精度惨不忍睹。
• 注意数值稳定性。当协方差矩阵P变得非正定时，EKF会发散。我建议加一个np.linalg.eig检查，或者用平方根滤波。
• OCV-SOC曲线的导数要平滑。如果拟合的多项式阶数太高，导数会剧烈震荡。我一般用6阶多项式，效果不错。
• 初始化要保守。初始协方差P设大一点，让滤波器先“相信”测量值，再慢慢收敛。

好了，这一讲的内容就到这里。EKF是处理非线性系统的经典方法，虽然现在有UKF、粒子滤波等更高级的方法，但EKF的线性化思想是基础，必须掌握扎实。

---