1、原理介绍

研究的对象是房间的温度吗？根据经验，房间的温度是25℃也就是说，期望值是25℃。

受空气循环和阳光的影响，房间温度为25℃左右波动。

我们用温度计取样室内温度一分钟。

假设过程噪声为 W ( k ) W(k) W(k)，方差大小为 Q = 0.01 Q=0.01 Q=0.01。一步转移矩阵 A = 1 A=1 A=1，噪声驱动矩阵 G = 1 G=1 G=1。

因此，系统状态方程和测量方程为： X ( k ) = A X ( k ? 1 ) G W ( k ? 1 ) X(k)=AX(k-1) GW(k-1) X(k)=AX(k?1) GW(k?1)

Z ( k ) = H X ( k ) + V ( k ) Z(k)=HX(k)+V(k) Z(k)=HX(k)+V(k)

假设初始房间温度为25.1℃，即 X ( 1 ) = 25.1 X(1)=25.1 X(1)=25.1。

初始值方差为 P ( 1 ) = 0.01 P(1)=0.01 P(1)=0.01，

初始测量值为 Z ( 1 ) = 24.9 Z(1)=24.9 Z(1)=24.9。

Z ( 1 ) Z(1) Z(1)可以作为初始估计值， X k f ( 1 ) = Z ( 1 ) = 24.9 X_kf(1)=Z(1)=24.9 Xk​f(1)=Z(1)=24.9。

2、MATLAB仿真代码

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%功能说明：Kalman滤波用在一维温度数据测量系统中
function Kalman_Temperature
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 120;  %采样点的个数为120，时间单位是分钟，可以理解为进行了120分钟测量
CON = 25; %根据经验判断，室内温度理论值应该是25度左右，在这个25的基础上会有波动

%对状态和测量初始化
X_except = CON*ones(1,N); %期望温度是25度，但真实温度不可能是这样
X = zeros(1,N);           %房间各个时刻真实温度值
X_kf = zeros(1,N);        %Kalman滤波处理器的状态，也叫估计值
Z = zeros(1,N);           %温度计测量值
P = zeros(1,N);           %

%赋初始值
X(1) = 25.1;    %假设房间温度初始值是25.1度
P(1) = 0.01;    %初始值的协方差
Z(1) = 24.9;    
X_kf(1) = Z(1); %初始测量值是24.9度，可以作为滤波器的初始估计状态

%噪声
Q = 0.01;   
R = 0.25;   
W = sqrt(Q)*randn(1,N); %方差决定噪声大小
V = sqrt(R)*randn(1,N); %方差决定噪声大小

%系统矩阵
F = 1;  %一步转移阵
G = 1;  %过程噪声驱动矩阵
H = 1;  %量测矩阵
I = eye(1); %本系统状态为一维
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%模拟房间温度和测量过程，并滤波
for k = 2:N
    %step1：随着时间的推移，房间温度真实波动变化
    %k时刻房间内的真实温度，对于温度计来说，这个真实值是不知道的，但它的存在是客观事实
    X(k) = F*X(k-1) + G*W(k-1); %动态方程
    %step2：随着时间的推移，获取实时数据
    %温度计对k时刻房间温度的测量，Kalman滤波是站在温度计角度进行的，
    %它不知道此时的真实状态X(k)，
    %只能利用本次测量值Z(k)和上一次估计值X_kf(k)来做处理，
    %最大限度的降低量测噪声R的影响，尽可能的逼近X(k)，这也是Kalman滤波目的所在
    Z(k) = H*X(k) + V(k);
    %step3：Kalman滤波
    %有了k时刻的观测Z(k)和k-1时刻的状态，那么就可以进行滤波了
    X_pre = F*X_kf(k-1);                 %状态预测
    P_pre = F*P(k-1)*F' + Q;             %协方差预测
    Kg = P_pre*inv(H*P_pre*H' + R)       %计算Kalman增益
    e = Z(k) - H*X_pre;                  %新息
    X_kf(k) = X_pre + Kg*e;              %状态更新
    P(k) = (I-Kg*H)*P_pre;               %协方差更新
end

%计算误差
Error_Messure = zeros(1,N);         %测量值与真实值之间的误差
Error_Kalman = zeros(1,N);          %Kalman估计与真实值之间的偏差
for k = 1:N
    Error_Messure(k) = abs(Z(k)-X(k));       
    Error_Kalman(k) = abs(X_kf(k)-X(k));    
end

t = 1:N;
%figure('Name','Kalman Filter Simulation','NumberTitle','off');
figure %画图显示
plot(t,X_except,'-b',t,X,'-r.',t,Z,'-ko',t,X_kf,'-g*');
legend('期望值','真实值','观测值','Kalman滤波值');
xlabel('采样时间/s');
ylabel('温度值/℃');

figure %误差分析图
plot(t,Error_Messure,'-b.',t,Error_Kalman,'-k*');
legend('测量偏差','Kalman滤波偏差');
xlabel('采样时间/s');
ylabel('温度偏差值/℃');
 

3、仿真结果

figure1

figure2