基于集合卡尔曼滤波在matlab中实现数据同化的问题

weixin_42855030 2019-04-18 04:48:32

我用matlab写的集合卡尔曼滤波（ENKF）的程序用于对遥感降水进行数据同化/融合，但是为什么我最后得到的融合降水却几乎变成了观测数据（也是一种遥感降水数据，我把它当做观测数据用于数据同化/融合），我的代码写错了吗？



clear all;

clc;



%%  卫星遥感降水数据的导入及参数的设定

Dec = 0.25; %降水数据空间分辨率

Site_number = 27;   %实测站点个数

A = 1;  %状态转移矩阵

H = 1;  %观测算子

data = xlsread('.\ENKF\1.xlsx');    %导入背景场降水数据和融合数据

Coord = xlsread('.\ENKF\流域内格点经纬度坐标.xlsx');  %导入尼洋河流域内网格点经纬度坐标

x1 = data(:,1); %背景场降水数据1

x2 = data(:,2); %观测场降水数据

R = 0.7;

Random_count = 1000;    %设定生成随机数个数

N=50;   %设定从随机数中抽取的个数

n=80;   %设定观测数



%%  生成背景场集合及观测场集合的随机扰动

for i=1:Site_number

    random = 0+randn(1,Random_count);    %随机生成符合均值为0，方差为1的共1000个随机数（背景场）

    randomGC = 0+sqrt(R)*randn(1,Random_count);  %随机生成符合均值为0，方差为R的共1000个随机数（观测场）

    randomcopy(i,:) = random;   %记录此行的所有随机数（背景场）

    randomcopyGC(i,:) = randomGC;   %同理（观测场）

    random_choose(i,:) = randsample(Random_count,N,false);  %从之前生成的随机数中随机选取其中N个数作为列的索引（背景场）

    random_chooseGC(i,:) = randsample(Random_count,n,false);    %同理（观测场）

    Random(i,:) = random(random_choose(i,:));    %按照上面的索引获取N个随机数的值（背景场）

    RandomGC(i,:) = randomGC(random_chooseGC(i,:));    %同理（观测场）

end



%%  对背景场集合及观测场集合添加扰动

x1copy = kron(x1,ones(1,N));    %将初始背景场x1复制为每列都相同的数组

x2copy = kron(x2,ones(1,n));    %将初始观测场x2复制为每列都相同的数组

Xb = x1copy+Random; %生成带扰动的背景场集合

XbGC = x2copy+RandomGC; %生成带扰动的观测场集合



%%  计算背景场及观测场的误差协方差矩阵

Xbmean = mean(Xb,2);    %求集合平均值

Xbmeancopy = kron(Xbmean,ones(1,N));    %将集合平均值复制为每列都相同的n列数组

XbRD = Xb-Xbmeancopy;   %计算背景场集合扰动项

Q = var(XbRD,0,2);  %计算背景场集合扰动项的方差

PbInitial = (XbRD*XbRD')/(N-1);    %计算初始背景场误差协方差矩阵

Ob = (RandomGC*RandomGC')/(n-1);   %计算观测场误差协方差矩阵



%%  计算局地裁剪函数并更新背景场误差协方差矩阵

for j=1:size(data,1)

    for m=1:size(data,1)

        Z(j,m) = Distance(Coord(j,3),Coord(j,4),Coord(m,3),Coord(m,4));   %计算距离Z的矩阵

        LCF(j,m) = Local_clipping(Z(j,m),0.75); %计算局地裁剪函数的矩阵

    end

    

end



Pb = LCF.*PbInitial;    %计算局地化之后的背景场误差协方差矩阵



%% ENKF数据融合

for ii=1:27

    X(1,:,ii) = Xb(ii,:);   %初始最优估计值

    P(ii,1) = Pb(ii,ii);    %初始误差协方差

    for jj=2:n+1     %迭代次数jj

%%  1、预测方程

        X(jj,:,ii) = A*X(jj-1,:,ii);    %计算模式预测值

        P1(ii,jj-1) = A*P(ii,jj-1)*A'+Q(ii);   %计算误差协方差预测值,代表第ii个站点的背景场误差协方差

%%  2、更新方程

        Kg(ii,jj-1) = P1(ii,jj-1)*H'/(H*P1(ii,jj-1)*H'+Ob(ii,ii));  %计算此状态卡尔曼系数，代表第ii个站点迭代到第jj次的Kg

        XGC = kron(XbGC(ii,:)',ones(1,N));

        X(jj,:,ii) = X(jj,:,ii)+Kg(ii,jj-1)*(XGC(jj-1,:)-H*X(jj,:,ii));    %计算此状态最优估计值,代表第ii个站点迭代到第jj次的值

        P(ii,jj) = (1-Kg(ii,jj-1)*H)*P1(ii,jj-1);    %计算此状态下最优协方差,代表第ii个站点迭代到第jj次的最优协方差

    end 

    Xmean(ii,:) = mean(X(:,:,ii),2)';   %其最后一列的数据为最终融合后的降水数据

    

end

function [ Z ] = Distance( lat1,lon1,lat2,lon2 )

    % Distance函数用于计算两个网格点的经纬度坐标之间的距离

Z = sqrt((lat1-lat2)^2 + (lon1-lon2)^2);



end



function [ P ] = Local_clipping( z,c )

    %Local_clipping函数是局地裁减函数，用于减少长距离伪相关的影响

if z>=0 && z<c

    P=1-(5/3)*(z/c)^2+(5/8)*(z/c)^3+(1/2)*(z/c)^4-(1/4)*(z/c)^5;

elseif z>=c && z<=2*c

    P=-(2/3)*(z/c)^-1+4-5*(z/c)+(5/3)*(z/c)^2+(5/8)*(z/c)^3-(1/2)*(z/c)^4+(1/12)*(z/c)^5;

else

    P=0;

end

...全文