用卡尔曼滤波器预测PM2.5浓度问题
这是我用卡尔曼滤波器去预测PM2.5浓度的一段代码，但是我的结果出现了如图片中预测值后面变成一条平直线的问题，请问我的代码哪里出现了问题，应该怎么改正% 清空环境
clc;
clear;
close all;
%% 1. 数据预处理
% 读取数据
= xlsread('dalian_pm2.5.xlsx', 'Sheet1', 'B2:B201');
observed_pm25 = data';
% 数据标准化（提升模型稳定性）
mu = mean(observed_pm25);
sigma = std(observed_pm25);
normalized_data = (observed_pm25 - mu) / sigma;
% 分割数据集
train_data = normalized_data(1:150);
test_data = normalized_data(151:end);
num_predictions = length(test_data);
%% 2. ARIMA模型优化
% --- 平稳性检验与自动定阶 ---
% ADF检验
= adftest(train_data);
if h == 0
fprintf('数据非平稳，自动差分...\n');
d = ndiffs(train_data); % 自动选择差分阶数
train_data_diff = diff(train_data, d);
else
d = 0;
train_data_diff = train_data;
end
% 通过AIC准则自动选择(p,q)
max_order = 3; % 最大尝试阶数
best_aic = Inf;
for p = 0:max_order
for q = 0:max_order
model = arima(p, d, q);
try
[~, ~, logL] = estimate(model, train_data_diff', 'Display', 'off');
current_aic = aicbic(logL, p + q);
if current_aic < best_aic
best_aic = current_aic;
best_p = p;
best_q = q;
end
catch
continue
end
end
end
fprintf('最优模型: ARIMA(%d,%d,%d)\n', best_p, d, best_q);
% 拟合最终模型
model = arima(best_p, d, best_q);
fit = estimate(model, train_data_diff');
%% 3. 改进的卡尔曼滤波实现
% --- 基于AR模型的扩展状态空间 ---
% 从ARIMA模型中提取AR系数
ar_coeff = cell2mat(fit.AR); % 获取AR项系数
% 构建状态空间
p = length(ar_coeff); % 状态维度=AR阶数
A = [ar_coeff;
eye(p-1), zeros(p-1, 1)]; % 状态转移矩阵
H = ; % 观测矩阵
% 初始化参数
Q = 0.001 * eye(p);% 过程噪声（需调优）
R = 0.5; % 观测噪声（需调优）
x_hat = '; % 初始状态
P = eye(p);
% --- 预测与滤波联合过程 ---
= forecast(fit, num_predictions, 'Y0', train_data_diff');
% 卡尔曼滤波递归修正
filtered_diff = zeros(1, num_predictions);
for k = 1:num_predictions
% 预测步骤
x_hat_minus = A * x_hat;
P_minus = A * P * A' + Q;
% 更新步骤（使用ARIMA预测作为观测输入）
K = P_minus * H' / (H * P_minus * H' + R);
z = forecast_diff(k); % ARIMA预测值
x_hat = x_hat_minus + K * (z - H * x_hat_minus);
P = (eye(p) - K * H) * P_minus;
filtered_diff(k) = H * x_hat;
end
% --- 逆标准化与差分还原 ---
% 还原差分
if d > 0
last_obs = train_data(end-d+1:end);
forecast_values = cumsum();
filtered_values = cumsum();
else
forecast_values = forecast_diff;
filtered_values = filtered_diff;
end
% 截取有效部分
forecast_values = forecast_values(end-num_predictions+1:end);
filtered_values = filtered_values(end-num_predictions+1:end);
% 逆标准化
forecast_values = forecast_values * sigma + mu;
filtered_values = filtered_values * sigma + mu;
test_data_orig = test_data * sigma + mu; % 原始尺度的测试数据
%% 4. 结果分析与可视化
% 检查维度
forecast_values_size = size(forecast_values);
test_data_orig_size = size(test_data_orig);
% 确保维度一致
if forecast_values_size ~= test_data_orig_size
% 根据需要调整维度
if forecast_values_size(2) == 1
forecast_values = forecast_values';
end
if test_data_orig_size(2) == 1
test_data_orig = test_data_orig';
end
end
% 计算RMSE
rmse_arima = sqrt(mean((forecast_values - test_data_orig).^2));
% 计算误差指标
rmse_arima = sqrt(mean((forecast_values - test_data_orig).^2));
rmse_kf = sqrt(mean((filtered_values - test_data_orig).^2));
fprintf('ARIMA RMSE: %.2f\n卡尔曼滤波RMSE: %.2f\n', rmse_arima, rmse_kf);
% 动态预测效果可视化
figure;
t = 1:num_predictions;
plot(t, test_data_orig, 'k-o', 'LineWidth', 1.5, 'MarkerSize', 4);
hold on;
plot(t, forecast_values, 'b--s', 'LineWidth', 1.5, 'MarkerSize', 4);
plot(t, filtered_values, 'r-.d', 'LineWidth', 1.5, 'MarkerSize', 4);
% 图形美化
title(sprintf('PM_{2.5}浓度预测对比 (ARIMA RMSE=%.2f, KF RMSE=%.2f)', rmse_arima, rmse_kf));
xlabel('时间（天）');
ylabel('PM_{2.5}浓度 (μg/m?)');
legend('实际值', 'ARIMA预测', '卡尔曼滤波修正', 'Location', 'best');
grid on;
set(gca, 'FontSize', 12, 'XTick', 1:5:num_predictions);
% 残差分析
figure;
subplot(2,1,1);
plot(t, test_data_orig - forecast_values, 'b', 'LineWidth', 1.5);
hold on;
plot(t, test_data_orig - filtered_values, 'r', 'LineWidth', 1.5);
title('预测残差序列');
legend('ARIMA残差', 'KF残差', 'Location', 'best');
grid on;
subplot(2,1,2);
histogram(test_data_orig - forecast_values, 'BinWidth', 5, 'FaceColor', 'b');
hold on;
histogram(test_data_orig - filtered_values, 'BinWidth', 5, 'FaceColor', 'r');
title('残差分布直方图');
legend('ARIMA', 'KF');
你的卡尔曼滤波预测结果变成一条平直线，可能是由于以下几个原因导致的：
### 可能的问题：
1. **状态转移矩阵 A 未正确定义**
你的代码中 `A = ;` 这一行是空的，状态转移矩阵 `A` 需要正确设置，否则状态不会正确更新，导致预测值趋于平稳。
2. **观测矩阵 H 未正确定义**
`H = ;` 也是空的，`H` 负责将状态变量映射到观测值，如果 `H` 设错，滤波器可能无法正确更新状态。
3. **初始状态 x_hat 为空**
`x_hat = ';` 这一行也是错误的，`x_hat` 需要初始化为合理的值，比如 `x_hat = zeros(p, 1);`。
4. **过程噪声 Q 和观测噪声 R 可能不合适**
你的 `Q = 0.001 * eye(p);` 和 `R = 0.5;` 可能需要调整，过小的 `Q` 可能导致滤波器对新观测值的调整过慢，导致预测值趋于平稳。
5. **ARIMA 预测值未正确传递给卡尔曼滤波器**
`z = forecast_diff(k);` 这一行可能有问题，确保 `forecast_diff` 变量正确存储了 ARIMA 预测值。
### 解决方案：
1. **正确定义状态转移矩阵 A**
如果你的 ARIMA 模型是 AR(p)，那么 `A` 应该是：
A = ;
这样可以保证状态向量正确更新。
2. **正确定义观测矩阵 H**
H = ;
这样 `H * x_hat` 就能正确映射到观测值。
3. **初始化状态变量**
x_hat = zeros(p, 1);
4. **调整噪声参数**
你可以尝试增大 `Q`，比如：
Q = 0.01 * eye(p);
这样可以让滤波器对新观测值更加敏感。
5. **确保 `forecast_diff` 变量正确**
在 `forecast` 之后，检查 `forecast_diff` 是否正确：
disp(forecast_diff);
确保它不是全零或异常值。
### 结论：
你的问题很可能是 `A`、`H`、`x_hat` 这些关键变量未正确定义，导致卡尔曼滤波器无法正确更新状态，最终预测值变成一条直线。尝试按照上述方法修正代码，并调整 `Q` 和 `R` 以优化滤波效果。
以上回复来自 FishC AI 小助理，如未能正确解答您的问题，请继续追问。 FishC 发表于 2025-3-27 16:10
你的卡尔曼滤波预测结果变成一条平直线，可能是由于以下几个原因导致的：
### 可能的问题：
更改后依然平稳不变化 蒋文德 发表于 2025-03-27 16:23
FishC 发表于 2025-3-27 16:10
你的卡尔曼滤波预测结 ...
如果更改后预测值仍然趋于平稳，可能的原因包括：
### **可能的核心问题**
1. **卡尔曼增益 K 过小**
- 如果 `Q`（过程噪声协方差）过小，或者 `R`（观测噪声协方差）过大，卡尔曼增益 `K` 可能会变得很小，导致滤波器对新观测值的调整幅度很小，最终预测值趋于平稳。
- **解决方案**：尝试增大 `Q` 或减小 `R`，例如：
Q = 0.1 * eye(p);% 增大过程噪声
R = 0.1; % 减小观测噪声
2. **状态转移矩阵 A 可能不稳定**
- `A` 需要正确反映 AR 过程的动态特性。如果 `A` 过于接近单位矩阵（即 `A ≈ I`），状态更新可能会变得缓慢，导致预测值趋于平稳。
- **解决方案**：确保 `A` 由 ARIMA 模型的 AR 系数正确构造，例如：
A = ;
3. **观测矩阵 H 可能不正确**
- `H` 需要正确映射状态变量到观测值。如果 `H` 设为单位矩阵 `I`，但状态变量的维度大于 1，可能会导致观测值更新不合理。
- **解决方案**：确保 `H` 只提取状态向量的第一维：
H = ;
4. **初始状态 x_hat 可能不合理**
- 如果 `x_hat` 初始化为全零，且 `A` 过于平稳，状态可能会一直保持接近零。
- **解决方案**：尝试用训练数据的均值初始化：
x_hat = mean(train_data) * ones(p, 1);
5. **ARIMA 预测值未正确传递**
- `forecast_diff(k)` 可能是一个常数或错误的值，导致滤波器的更新无效。
- **解决方案**：在 `forecast` 之后检查 `forecast_diff`：
disp(forecast_diff);
### **调试步骤**
1. **检查 `K` 是否过小**
在 `K = P_minus * H' / (H * P_minus * H' + R);` 之后，添加：
disp(K);
如果 `K` 过小（接近零），说明滤波器对新观测值的调整幅度太小。
2. **检查 `x_hat` 是否变化**
在 `x_hat = x_hat_minus + K * (z - H * x_hat_minus);` 之后，添加：
disp(x_hat);
如果 `x_hat` 变化很小，说明状态更新不充分。
3. **检查 `forecast_diff` 是否合理**
在 `z = forecast_diff(k);` 之后，添加：
disp(z);
确保 `z` 不是一个常数或异常值。
### **总结**
- **增大 `Q`，减小 `R`**，让滤波器更敏感。
- **确保 `A` 和 `H` 正确构造**，避免状态更新过慢。
- **检查 `K`、`x_hat` 和 `forecast_diff`**，确保滤波器在每一步都在更新状态。
尝试这些方法后，看看预测值是否仍然趋于平稳。
以上回复来自 FishC AI 小助理，如未能正确解答您的问题，请继续追问。 蒋文德 发表于 2025-3-27 16:23
更改后依然平稳不变化
% 原代码：
for k = 1:length(X_test)
% --- 预测步骤 ---
x_hat_minus = A * x_hat;
P_minus = A * P * A' + Q;
% --- 更新步骤 ---
K = P_minus * H' / (H * P_minus * H' + R);
x_hat = x_hat_minus + K * (X_test(k) - H * x_hat_minus);
P = (eye(3) - K * H) * P_minus;
predicted_pm25(k) = x_hat(1);
end
% 修改为：
for k = 1:length(X_test)
% --- 预测步骤 ---
x_hat_minus = A * x_hat;
P_minus = A * P * A' + Q;
% 保存预测值（使用先验估计）
predicted_pm25(k) = x_hat_minus(1);
% --- 不进行更新步骤，直接使用预测状态 ---
x_hat = x_hat_minus;% 状态更新为先验估计
P = P_minus; % 协方差更新为先验协方差
end
我在将这一部分修改后出现未定义函数或变量 'X_test'。这个问题应该怎么修改
蒋文德 发表于 2025-03-27 16:41
蒋文德 发表于 2025-3-27 16:23
更改后依然平稳不变 ...
`X_test` 未定义的错误通常是因为变量 `X_test` 在代码中没有正确赋值或命名不一致。你可以按照以下步骤检查并修正：
### **可能的原因及解决方案**
1. **检查 `X_test` 是否在代码中定义**
在 `for` 循环之前，确保 `X_test` 已经正确赋值。例如，如果 `X_test` 是测试数据，你可能需要这样定义：
X_test = test_data;% 确保 test_data 已经正确赋值
如果 `test_data` 也未定义，请检查数据预处理部分。
2. **检查变量命名是否一致**
代码中可能使用了不同的变量名，比如 `test_data` 或 `forecast_diff`，但 `for` 循环中使用的是 `X_test`。请确保变量名一致，例如：
for k = 1:length(test_data)
x_hat_minus = A * x_hat;
P_minus = A * P * A' + Q;
predicted_pm25(k) = x_hat_minus(1);
x_hat = x_hat_minus;
P = P_minus;
end
如果 `X_test` 其实是 `forecast_diff`，请改为：
for k = 1:length(forecast_diff)
x_hat_minus = A * x_hat;
P_minus = A * P * A' + Q;
predicted_pm25(k) = x_hat_minus(1);
x_hat = x_hat_minus;
P = P_minus;
end
3. **检查 `X_test` 是否在代码执行前被清除**
如果 `X_test` 在 `clear` 之后没有重新赋值，MATLAB 会报错。请确保 `X_test` 在 `clear` 之后被正确赋值。
### **调试方法**
在 `for` 循环之前，添加：
disp('X_test 内容:');
disp(X_test);
如果 `X_test` 为空或未定义，说明数据未正确传递，需要回溯数据处理部分。
### **总结**
- 确保 `X_test` 在 `for` 循环前正确赋值。
- 检查变量命名是否一致，例如 `X_test` 是否应该是 `test_data` 或 `forecast_diff`。
- 确保 `X_test` 没有在 `clear` 之后被误删除。
- 使用 `disp(X_test);` 检查变量是否正确赋值。
以上回复来自 FishC AI 小助理，如未能正确解答您的问题，请继续追问。 鱼油把 forecast_diff(k)（即 ARIMA 多步预测）当作了卡尔曼滤波的观测量 z，而并没有用实际的测试数据（也就是真实观测到的 PM2.5）来更新卡尔曼滤波器。
这样一来，卡尔曼滤波器就只能围绕着 ARIMA 的预测值进行修正，而 ARIMA 的多步预测在稳态情况下往往会向平均值收敛，所以滤波输出也就“跟着”收敛为一条平线。
也就是说，在卡尔曼滤波的“更新”步骤里，应该使用真实测试集（即 test_data，如果差分过就用相应的 test_data_diff）来作为观测值 z，而不是 forecast_diff(k)。
参考：
for k = 1:num_predictions
% (1) 预测步骤
x_hat_minus = A * x_hat;
P_minus = A * P * A' + Q;
% (2) 更新步骤 —— 用真实观测来更新
K = P_minus * H' / (H * P_minus * H' + R);
z = test_data_diff(k);% 用实际的测试数据(差分后)作为观测
x_hat = x_hat_minus + K*(z - H*x_hat_minus);
P = (eye(p) - K*H)*P_minus;
filtered_diff(k) = H*x_hat;
end 不是，鱼C论坛都能讨论matlab了吗
页:
[1]