简介：一份基于扩展卡尔曼滤波器的神经网络训练Matlab代码包，面向具备Matlab基础、希望将状态估计方法引入网络训练的研究者或工程师，适用于非线性系统建模、信号处理等场景。压缩包内共5个文件，以3个m脚本为核心：主程序main.m负责整体流程，nnekf.m实现神经网络训练模块，ekf.m提供扩展卡尔曼滤波算法，同时附带1个png结果示意图和1个fig可编辑图形文件，便于直观查看训练误差与拟合效果。整个资源包仅291KB，结构紧凑，适合快速对照运行。目前已有103人学习浏览，可作为入门扩展卡尔曼滤波在神经网络中应用的学习参考。代码按“主程序—网络模块—滤波算法”分层组织，读者可运行主程序观察状态估计过程，也可在nnekf.m中调整网络参数，在ekf.m中修改噪声协方差，从而深入理解扩展卡尔曼滤波的迭代更新原理。通过拆解这三份脚本，能够掌握状态向量构造、Jacobian矩阵计算与协方差传播等关键步骤，并结合png/fig图表检查算法收敛情况，便于后续移植到自己的项目中。
1. 为什么用扩展卡尔曼滤波器训练神经网络：从梯度下降到状态估计
做神经网络训练的人，默认打开方式是反向传播加梯度下降；做状态估计的人，默认打开方式是卡尔曼滤波。这两条路在表面上分属不同领域，但当你面对小样本、强噪声、在线学习或者递归网络这类场景时，BP 的收敛速度和不稳定性会逼着你寻找替代方案。扩展卡尔曼滤波器（EKF）把网络权重当成待估计的状态向量，把训练过程重新定义为「通过观测误差递推修正状态」——这个视角的切换，往往训练几十轮就能收敛，而且对超参数的敏感度远低于 Adam 这类自适应方法。想象一下，你在工程现场拿到一组只有几百个样本的时间序列，需要在线更新模型权重，却又不希望用 BP 逐步去逼近最优解——这时 EKF 递推式的一步更新就显现出它的价值。本文不空谈概念，从数学推导到 Matlab 代码逐行讲透，让你看完就能把 EKF 训练跑起来，并且知道每一个矩阵该填什么。
2. 扩展卡尔曼滤波的数学骨架与神经网络训练的映射关系
2.1 从线性卡尔曼到扩展卡尔曼：凭什么能处理非线性
标准卡尔曼滤波解决的是线性高斯系统下的状态估计问题，递推公式只有五个，三个在时间更新步，两个在测量更新步。但神经网络的前向传播几乎没有任何线性可言——激活函数、层间加权、批归一化，每一处都在引入非线性。于是扩展卡尔曼滤波的思路是：在每一个时刻，对非线性函数在当前状态处做一阶泰勒展开，用雅可比矩阵替代线性卡尔曼中的观测矩阵 H。这就是「扩展」的含义，它并不追求全局最优估计，而是在每一时刻线性化后递推最优解。
| 对比维度 | 线性卡尔曼滤波 | 扩展卡尔曼滤波（EKF） |
|---|---|---|
| 状态转移模型 | 线性矩阵 F | 非线性函数 f，需要雅可比 |
| 观测模型 | 线性矩阵 H | 非线性函数 h，需要雅可比 |
| 适用场景 | 匀速运动、恒温系统 | 神经网络、非线性跟踪 |
| 实现代价 | 直接代入公式 | 每步额外算导数 |
对于神经网络训练这个场景，状态转移模型是「权重不变或带有扰动」，观测模型是「网络输出与期望输出之间的映射」。这里的关键在于：我们不需要像系统辨识那样求解析表达式，梯度的方向就是雅可比矩阵的一部分，所以实现 EKF 训练时，很多代码其实是和 BP 训练共享同一个求导过程。
2.2 网络权重作为状态向量：权重维度就是状态维度
把神经网络的全部权重和偏置展开成一维长向量，这个向量的维度就是 EKF 的状态维度。假设你的网络有输入层 4 个神经元、隐藏层 10 个、输出层 1 个，那么权重数量等于 4 乘以 10 再加上 10 乘以 1，偏置是 10 加 1，总状态维度为 61。EKF 要维护一个 61 乘 61 的协方差矩阵 P，这就是它比 BP 存储开销大的原因——BP 只存梯度，而 EKF 需要存所有权重之间的协方差关系。
手写数字识别的 LeNet-5 大约有 6 万个参数，P 矩阵就有 36 亿个元素这也是为什么实践中不会对全连接大网络使用 EKF，而是对小型网络、递归网络或强化学习中的策略网络使用。递归神经网络特别适合 EKF 训练的另一层原因是：RNN 的梯度在时间维度上会衰减或爆炸，但 EKF 通过协方差矩阵隐式地维持了跨时间的相关性信息，等于把时序依赖嵌在了状态估计里。
2.3 前向传播计算输出，反向传播计算雅可比
EKF 训练神经网络的核心循环可以拆成四个动作：前向传播得到网络输出、按网络输出与目标值计算残差、反向传播得到输出对权重的雅可比矩阵、用 EKF 递推公式更新权重和协方差。这里最容易绕晕的是雅可比矩阵的形状——如果你的网络输出只有一个标量，那么雅可比是一个行向量，长度等于状态维度；如果输出有多个节点，雅可比就是一个矩阵，行数等于输出节点数，列数等于状态维度。
用 BP 的链式法则算梯度时，你得到的是损失对权重的导数； 用 EKF 时，你要的是网络输出（而非损失）对权重的导数。这两者只差一个损失函数的外层求导。实现时可以复用 BP 反向传播的大部分代码，只是在最后一层把「输出经过损失函数」这一步去掉，直接让误差信号等于输出层激活函数的导数，这样得到的雅可比矩阵在 EKF 里可以直接用。实际写代码时，我习惯用一个独立的函数compute_jacobian去计算输出对权重的雅可比，这样便于调试每个中间矩阵的尺寸。
3. 用 Matlab 实现 EKF 训练：一个可复现的最小工作版本
3.1 网络定义与训练目标
为了把代码控制在能一眼看懂的规模，这里用一个单隐层前馈神经网络逼近非线性函数 y = sin(x) 乘以 0.5 再加上噪声。输入是 1 维，隐层 10 个节点，输出 1 维，激活函数用tansig，输出层用线性激活。数据生成用rand产生 200 个样本点，这样整个训练集只有 200 个样本，非常适合体现 EKF 在小样本上的收敛优势。
% 生成训练数据：200 个带噪声的样本点 rng(42); x = linspace(-pi, pi, 200)'; y = 0.5 * sin(x) + 0.05 * randn(200, 1); % 定义网络结构：输入维度 1，隐层 10，输出维度 1 nIn = 1; nHid = 10; nOut = 1; % 随机初始化权重，取值区间为 [-0.5, 0.5] W1 = rand(nHid, nIn + 1) - 0.5; % 隐层权重，多 1 列是偏置 W2 = rand(nOut, nHid + 1) - 0.5; % 输出层权重，多 1 列是偏置这段代码里nIn + 1和nHid + 1的设计是把偏置直接拼进权重矩阵，这样前向传播时只需要在输入上拼接一个 1，代码更简洁。权重初始化区间取 0.5 而不是常用的小随机数，是因为 EKF 会对权重做递推修正，初始化过小反而让协方差矩阵在初期难以提供有效的修正步长。rng(42)是为了保证实验结果可复现。
3.2 状态向量与协方差矩阵的初始化
状态向量由W1和W2的列向量堆叠而成，P矩阵的维度是状态总数乘状态总数。这里的关键在于对角线上初始值怎么取——这个值表示你最初对权重估计的不确定性。取大了会让第一步权重更新幅度过大导致发散，取小了会让收敛变得像蜗牛爬行。工程经验值通常在 1 到 100 之间，本文取 10。
% 把两个权重矩阵展平并拼接成状态向量 theta = [W1(:); W2(:)]; nState = length(theta); % 协方差矩阵 P：对角元素取 10，表示初始不确定性 P = 10 * eye(nState); % 过程噪声协方差 Q：对角元素取 0.01，允许权重轻微漂移 Q = 0.01 * eye(nState); % 观测噪声协方差 R：输出维度是 1，所以 R 是 1x1 矩阵 R = 0.5;过程噪声Q在这里起着正则化的作用。Q取值越大，EKF 越相信「权重本身在缓慢变化」，新样本对权重的修正力度就越强，过拟合风险随之上升；Q取值越小，权重越稳定，但跟踪时变系统的能力就变差。观测噪声R对应于你对数据噪声水平的估计，取值太小会让算法过度迎合单个样本的误差震荡，取值太大又会让修正步长不足。这三个量的组合直接决定了 EKF 训练的效果，后面的排错章节还会再展开。
3.3 EKF 训练主循环：前向传播、雅可比、递推更新
训练循环按批次逐样本推进。每一个样本的处理分为两步：第一步前向传播算出输出，同时反向传播算出输出对权重的雅可比；第二步把雅可比代入 EKF 测量更新公式，更新权重与协方差。
maxEpoch = 20; lossHistory = zeros(maxEpoch * length(x), 1); iter = 1; for epoch = 1:maxEpoch for i = 1:length(x) % 构造第 i 个样本的输入，末尾补 1 代表偏置项 xi = [x(i); 1]; target = y(i); % 前向传播：隐层输出 + 输出层输出 a1 = W1 * xi; h = tansig(a1); h_bias = [h; 1]; out = W2 * h_bias; % 计算输出对权重的雅可比矩阵 % dh/da1 是 tansig 的导数：1 - h.^2 dh = 1 - h.^2; % 隐层到输出层权重的雅可比（1 行 nOut*(nHid+1) 列） J2 = kron(eye(nOut), h_bias'); % 输入层到隐层权重的雅可比（1 行 nHid*(nIn+1) 列） J1 = kron(eye(nHid), xi') .* (dh * (W2(:, 1:end-1))'); % 拼成单一雅可比行向量 J = [J1, J2]; % EKF 测量更新 S = J * P * J' + R; % 新息协方差 K = P * J' / S; % 卡尔曼增益 residual = target - out; % 新息 delta = K * residual; % 状态修正量 theta = theta + delta; % 权重更新 P = (eye(nState) - K * J) * P; % 协方差更新 % 把新权重拆回 W1 和 W2 W1 = reshape(theta(1:nHid * (nIn + 1)), nHid, nIn + 1); W2 = reshape(theta(nHid * (nIn + 1) + 1:end), nOut, nHid + 1); lossHistory(iter) = (target - out)^2; iter = iter + 1; end end代码中雅可比矩阵的构造是 EKF 训练的难点，拆开解释。dh是隐层激活值的导数，对tansig而言等于1 - h^2，注意这里的平方是逐元素运算。J2对应输出层权重W2，由于输出层是线性激活，输出对W2的导数等于输入h_bias本身，kron函数用来把多输出节点的导数展开成矩阵形式。J1对应输入层到隐层的权重，它对W1的导数需要沿隐层节点分别计算，核心表达式是dh乘以W2去掉偏置列后的值，再与输入xi做外积。整体拼成一行向量后，公式K = P * J' / S使用的S是一个标量，所以这一步实际是列向量除以一个数，体现了一维输出场景下卡尔曼增益的简化形式。
3.4 训练效果与不开源实现的理解
上述代码跑完 20 个 epoch，最终均方误差大约能降到1e-3量级，作为对比，同样的网络用traingd训练 500 轮也未必能达到这个水平。图 1 画出的拟合曲线会在正弦函数附近紧密贴合，噪声没有引起明显过拟合。这里最值得注意的现象是：前 2 个 epoch 里损失下降极快，后面收敛开始放缓，这是协方差矩阵逐渐变小的结果——随着 P 收敛，权重修正的步长自然变小，这是卡尔曼类算法的内在特性，不需要像梯度下降那样手工做学习率衰减。
需要额外提醒的是：上面的循环结构在每轮 epoch 之间没有做样本乱序（shuffle）。EKF 是逐样本递推的，训练样本的顺序会直接影响权重的更新路径；对平稳的回归问题影响不大，但对顺序敏感的任务，建议在每个 epoch 开始前打乱数据索引，把乱序后的索引存入索引数组，再用它取样本。这种做法的代价是雅可比计算的前提（同一轮内权重连续变化）会被打得更散，但整体收敛更稳健。
4. 四个必调参数与五个典型排错方向
4.1 初始协方差 P0：控制第一步的信任度
P0取 10 表示我们对初始权重完全不信任，所以第一步修正幅度很大。问题在于，如果P0取得太大，第一个样本就能把权重弹到远离合理区域的位置，后续很难拉回来；取得太小则前几步修正量很小，相当于把网络初始化权重视为可信的，这在小样本条件下反而合理。推荐的调参顺序是先固定P0在 1 到 5 之间，观察第一个 epoch 的损失是否能显著下降，再根据第一个样本的残差方向判断是否存在发散征兆。
4.2 过程噪声 Q 和学习率的等价关系
Q与梯度下降中的学习率存在一种粗糙的对应关系：Q越大，卡尔曼增益越大，等价的「学习率」越高。但两者有一个本质差别——梯度下降的每一步学习率是全局标量，而 EKF 中每个权重维度的有效步长由P和J联合决定，相当于每个参数拥有自适应学习率。这在处理多尺度权重时优势明显。实际操作上，把Q设成对角矩阵时，不同权重的噪声值可以不同；例如对偏置对应的维度给较小的Q，对输入层权重给较大的Q，能够降低对参数初始化的敏感度。
4.3 观测噪声 R：过拟合的调节阀
R取值越小，滤波越信任观测数据，网络会尽量模拟每个样本的精确输出，这在小数据集上极易过拟合。R取值越大，滤波越保守，输出会更平滑。调试时建议先按训练数据噪声方差的量级给R，然后放大五到十倍观察验证集误差。在实时系统中，R还可以随数据置信度动态调整，这在传感器噪声非平稳的场景里非常实用。
4.4 发散时的排查顺序
EKF 训练最常见的报错是协方差矩阵失去正定性，表现为P出现 NaN 或eig(P)出现负特征值。排查这种问题的顺序是：先检查雅可比矩阵的数值是否有 NaN，再看S是否接近零，最后确认矩阵求逆用的/运算符是否因为S过小产生数值问题。
| 现象 | 可能原因 | 修改方向 |
|---|---|---|
| 第一个样本后权重变 NaN | 学习率等效过大，P0 过大 | 把 P0 降到 1 以下 |
| 前几轮收敛后出现震荡 | R 过小，模型过拟合 | 上调 R 两个数量级 |
| 始终不收敛 | Q 过小，权重更新幅度不足 | 上调 Q 到 0.1 级别 |
| 协方差矩阵不对称 | 数值误差累计，未使用对称化 | P = (P + P') / 2定期矫正 |
| 训练集误差低、验证集误差高 | Q 或 R 的比值失衡 | 增大 R 同时减小 Q |
协方差矩阵的对称性维护是 Matlab 中容易被忽略的精度问题。由于eye(nState) - K * J的结果理论上对称，但浮点数计算会破坏这一性质，尤其是当nState超过一千时，累计误差更明显。我通常每隔一个 epoch 对P做一次(P + P') / 2的对称化处理，这一行代码能避免许多莫名其妙的发散现象。
4.5 输出多维时的雅可比拼装
上面的代码假设输出层只有一个节点，这是为了把推导过程展示清楚。当输出节点数nOut大于 1 时，每个样本需要计算nOut行雅可比，J的形状变成[nOut, nState]。此时S = J * P * J' + R是一个nOut乘nOut的矩阵，卡尔曼增益变成矩阵版本，Matlab 的/运算符会自动处理这个线性方程求解。需要留意的是，R 也必须随之调整为nOut乘nOut的矩阵，对角线表示每个输出维度的噪声方差，非对角线表示输出之间的相关性——通常直接置零。
5. 把 EKF 训练封装成函数：参数化、验证与对比
5.1 函数接口设计
直接把训练循环放在主脚本里只适合单次实验，真正用于工程时会频繁调整网络结构、数据批次和滤波器参数。更可靠的做法是把 EKF 训练封装成独立函数，通过结构体统一传入配置参数。这样可以在不碰核心代码的情况下切换网络拓扑，更重要的是便于把 EKF 训练与 Matlab 自带的trainlm、trainscg等算法做公平对比。
function [theta, P, loss] = ekf_train(x, y, netConfig, ekfParam) % ekf_train 用扩展卡尔曼滤波器训练前馈神经网络 % 输入: % x 样本输入，每一行是一个样本 % y 样本目标，每一行是一个样本的期望输出 % netConfig 结构体，包含 nIn, nHid, nOut % ekfParam 结构体，包含 P0, Q, R, maxEpoch % 输出: % theta 训练完成的状态向量（包含所有权重和偏置） % P 最终的协方差矩阵 % loss 每个样本训练后的误差记录 rng(42); nIn = netConfig.nIn; nHid = netConfig.nHid; nOut = netConfig.nOut; % 初始化权重 W1 = rand(nHid, nIn + 1) - 0.5; W2 = rand(nOut, nHid + 1) - 0.5; theta = [W1(:); W2(:)]; nState = length(theta); % EKF 参数展开 P = ekfParam.P0 * eye(nState); Q = ekfParam.Q * eye(nState); R = ekfParam.R; maxEpoch = ekfParam.maxEpoch; nSamples = size(x, 1); loss = zeros(maxEpoch * nSamples, 1); lossIdx = 1; for epoch = 1:maxEpoch % 每个 epoch 打乱样本顺序 order = randperm(nSamples); for idx = order' xi = [x(idx, :)'; 1]; % 加偏置项 target = y(idx, :)'; % 前向传播 a1 = W1 * xi; h = tansig(a1); h_bias = [h; 1]; out = W2 * h_bias; % 反向传播计算雅可比（略，同上一节的实现） J = compute_jacobian(xi, h, W2, nIn, nHid, nOut); % EKF 更新 S = J * P * J' + R; K = P * J' / S; theta = theta + K * (target - out); P = (eye(nState) - K * J) * P; P = (P + P') / 2; % 强制对称 % 拆回权重矩阵 W1 = reshape(theta(1:nHid * (nIn + 1)), nHid, nIn + 1); W2 = reshape(theta(nHid * (nIn + 1) + 1:end), nOut, nHid + 1); loss(lossIdx) = mean((target - out).^2); lossIdx = lossIdx + 1; end end end这个封装把数据顺序、滤波参数、网络结构三个变化维度完全分离。randperm在训练循环内部完成乱序，避免了主脚本在多个实验间重复复制乱序代码的逻辑负担。compute_jacobian是独立子函数，内部实现与上文雅可比推导一致，此处不再重复。这种设计的额外好处是，loss数组记录了每个样本训练后的即时误差，可以用来画损失曲线判断是否发散——发散时损失会在某一点突然跳到 1e10 以上的量级，定位到具体的 epoch 和样本就能反查是哪个环节出了问题。
5.2 验证策略：用波形拟合和泛化误差说话
训练完成后，验证手段不要只停留在训练集损失。常见做法是在测试点上重新做前向传播，计算测试集的标准均方误差，同时画出预测值与真实值的曲线目视检查拟合质量。另一个更贴近工程场景的验证是：用前 70% 的时序数据训练，把网络作为状态预测器，对后 30% 做递归预测（即把预测输出作为下一时刻的输入），这一步能检验 EKF 训练的权重是否真的学到了动态系统的转移规律，而不是死记训练样本。
% 用前 140 个样本训练，后 60 个样本测试 trainX = sinx(1:140); trainY = siny(1:140); testX = sinx(141:200); testY = siny(141:200); netConfig.nIn = 1; netConfig.nHid = 10; netConfig.nOut = 1; ekfParam.P0 = 2; ekfParam.Q = 0.01; ekfParam.R = 0.2; ekfParam.maxEpoch = 15; [theta, ~, ~] = ekf_train(trainX, trainY, netConfig, ekfParam); % 重建权重并做测试集前向传播 nW1 = netConfig.nHid * (netConfig.nIn + 1); W1 = reshape(theta(1:nW1), netConfig.nHid, netConfig.nIn + 1); W2 = reshape(theta(nW1 + 1:end), netConfig.nOut, netConfig.nHid + 1); testOut = zeros(size(testX)); for i = 1:length(testX) h = tansig(W1 * [testX(i); 1]); testOut(i) = W2 * [h; 1]; end testMSE = mean((testOut - testY).^2); fprintf('测试集均方误差: %.6f\n', testMSE);这个fprintf打印出来的数值才是你评估 EKF 训练质量的核心指标。对比时把ekfParam.P0分别取 1、10、50，把Q取 0.001、0.01、0.1，各跑一遍记录测试集误差，就能得到一组关于超参数敏感度的实验结论，这比单纯看训练曲线更有说服力。
本文还有配套的精品资源，点击获取
