用卡尔曼滤波器预测PM2.5浓度问题

蒋文德

这是我用卡尔曼滤波器去预测PM2.5浓度的一段代码，但是我的结果出现了如图片中预测值后面变成一条平直线的问题，请问我的代码哪里出现了问题，应该怎么改正


% 清空环境
clc;
clear;
close all;

%% 1. 数据预处理
% 读取数据
[data, ~, ~] = 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检验
[h, ~] = 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 = [1, zeros(1, p-1)];        % 观测矩阵

% 初始化参数
Q = 0.001 * eye(p);  % 过程噪声（需调优）
R = 0.5;             % 观测噪声（需调优）
x_hat = [train_data_diff(end:-1:end-p+1)]'; % 初始状态
P = eye(p);

% --- 预测与滤波联合过程 ---
[forecast_diff, ~, ~] = 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([last_obs, forecast_diff]);
    filtered_values = cumsum([last_obs, filtered_diff]);
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');

FishC

你的卡尔曼滤波预测结果变成一条平直线，可能是由于以下几个原因导致的：

### 可能的问题：
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 = [ar_coeff; eye(p-1), zeros(p-1,1)];

   这样可以保证状态向量正确更新。

2. **正确定义观测矩阵 H**  

H = [1, zeros(1, p-1)];

   这样 `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
你的卡尔曼滤波预测结果变成一条平直线，可能是由于以下几个原因导致的：

### 可能的问题：

更改后依然平稳不变化

FishC

蒋文德 发表于 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 = [ar_coeff; eye(p-1), zeros(p-1,1)];

3. **观测矩阵 H 可能不正确**  
   - `H` 需要正确映射状态变量到观测值。如果 `H` 设为单位矩阵 `I`，但状态变量的维度大于 1，可能会导致观测值更新不合理。  
   - **解决方案**：确保 `H` 只提取状态向量的第一维：
   

H = [1, zeros(1, p-1)];

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'。这个问题应该怎么修改

FishC

蒋文德 发表于 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了吗