白噪声估值器及其在信号处理中的应用

在上一篇博客中，我讲到信息融合涉及的领域中包含了系统白噪声的估计问题，即根据系统测量值来估计系统的观测误差和系统误差的大小。下面针对上述白噪声估计算法进行讲解。

定理1：

针对如下状态空间模型：

                                                           

                                                                

存在最优输入白噪声估值器

                                                                        

                                                                          （N<0）

                                                       （N>0）

估计误差方差阵如下：

                                                                                 （N<=0）

                                             

证明：

根据射影定理，有：

                （1）

其中有：

                                                       （2）

所以有：

                                                                          （3）

将式2和式3带入式1，有：

                                                            （4）

下面来证明误差方差阵：

                                         （5）

                         （6）

综上，定理1得证。

定理2：

最优输入白噪声递推平滑器如下所示：

其中：

误差方差阵如下：

证明：

根据射影定理可得：

（7）

因为有：

                                               （8）

（9）

根据式8递推表达式可知：

（10）

将式10带入式8求期望可得：

            （11）

其中：

将式11带入式7即可得到定理中式子。

对于求误差方差阵，利用定理1中的证明方式即可得到，这里不再赘述。

定理3：最优观测白噪声估值器如下

初值为：

其中：

误差方差阵如下：

证明：

因为：

（12）

同定理2的证明过程，由式8和式9可得上述定理表达式。

仿真例子

考虑随机系统

其中有：

b(t)为取值为1或者0的白噪声，概率为：

且有Qg=7,Qv=0.1,则仿真效果图如下：

仿真代码如下：

clc;
clear;
N=50;
Lamda=0.3;
Qg=7;
Q=Lamda*Qg;
R=0.1;
RandnumBt=rand(1, N);
for i=1:N
    if RandnumBt(i)<0.5
       RandnumBt(i)=0; 
    else
       RandnumBt(i)=1; 
    end
end
RandNum=normrnd(0,sqrt(Qg),1,N);
for i=1:N
    W(i)=RandNum(i)*RandnumBt(i);
end
RandNum1=normrnd(0,sqrt(R),1,N);
%%卡尔曼滤波参数
A=[1,0;0.3,-0.5];
G=[-1;2];
H=[1,1];
x(:,1)=[0;0];
y(1)=H*x(:,1)+RandNum1(1);
for i=2:N
    x(:,i)=A*x(:,i-1)+G*W(i-1);
    y(i)=H*x(:,i)+RandNum1(i);
end
Pk(:,:,1)=[Q,0;0,R];
Pkk(:,:,1)=zeros(2,2);
K(:,1)=zeros(2,1);
X(:,1)=x(:,1);
Xk(:,1)=x(:,1);
I=eye(2);
for i=2:N
    Xk(:,i)=A*X(:,i-1);
    Pkk(:,:,i)=A*Pk(:,:,i-1)*A'+G*Q*G';
    K(:,i)=Pkk(:,:,i)*H'/(H*Pkk(:,:,i)*H'+R);
    e(i)=y(i)-H*Xk(:,i);
    X(:,i)=Xk(:,i)+K(:,i)*e(i);
    Pk(:,:,i)=(I-K(:,i)*H)*Pkk(:,:,i);
end
N1=3;
for i=1:N
    Persai(:,:,i)=A*(I-K(:,i)*H);
    Qe(i)=H*Pkk(:,:,i)*H'+R;
end
for i=1:N-2
    M1e(i)=Q*G'*H'/Qe(i);
    M2e(i)=Q*G'*Persai(:,:,i+1)'*H'/Qe(i+2);%
    M3e(i)=Q*G'*(Persai(:,:,i+1)*Persai(:,:,i+2))'*H'/Qe(i+2);
end
for i=1:N-5
    wjian(1,i)=M1e(i)*e(i+1);%%一阶平滑滤波
    wjian(2,i)=wjian(1,i)+M2e(i)*e(i+2);%%二阶平滑滤波
    wjian(3,i)=wjian(2,i)+M3e(i)*e(i+3);%%三阶平滑滤波
end
wjian(3,N)=W(N);
wjian(3,N-1)=W(N-1);
wjian(3,N-2)=W(N-2);
t=1:N;
figure(1);
subplot(2,2,1);
plot(t,x(1,:),'r',t,X(1,:));
subplot(2,2,2);
plot(t,x(2,:),'r',t,X(2,:));
figure(2);
subplot(2,2,1);
plot(t,W,'*r',t,wjian(1,:),'*k');
title('一次滤波');
subplot(2,2,2);
plot(t,W,'*r',t,wjian(2,:),'*k');
title('二次滤波');
subplot(2,2,3);
plot(t,W,'*r',t,wjian(3,:),'*k');
title('三次滤波');