基于集合卡尔曼滤波在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

...全文
14758 13 打赏 收藏 举报
写回复
13 条回复
切换为时间正序
当前发帖距今超过3年,不再开放新的回复
发表回复

您好,请问您的问题解决了吗

  • 打赏
  • 举报
回复 1
yyyy6yy1 03-29
@qq_49547341 请问您的问题解决了吗
  • 举报
回复
qq_37946279 2021-04-03
可以共享下数据吗 让孩子 看看格式
  • 打赏
  • 举报
回复
显然不显然 2020-10-26
引用 9 楼 cfttony 的回复:
[quote=引用 8 楼 显然不显然 的回复:]您好,我运行了您的代码认为问题应该出在“对背景场集合及观测场集合添加扰动”这一节。观测场不是通过观测值得到的,应该要用Y = H*Xb+RandomGC计算,最后结果就是很接近背景场了。


8楼的兄弟,你好。如果按照你说的Y = H*Xb+RandomGC,那观测值就不起作用了吧。那还要观测值干啥呢?[/quote]是我理解错了,非常抱歉
  • 打赏
  • 举报
回复
cfttony 2020-09-12
引用 8 楼 显然不显然 的回复:
您好,我运行了您的代码认为问题应该出在“对背景场集合及观测场集合添加扰动”这一节。观测场不是通过观测值得到的,应该要用Y = H*Xb+RandomGC计算,最后结果就是很接近背景场了。
8楼的兄弟,你好。如果按照你说的Y = H*Xb+RandomGC,那观测值就不起作用了吧。那还要观测值干啥呢?
  • 打赏
  • 举报
回复
显然不显然 2020-08-07
您好,请问您已经解决结果接近观测值的问题了吗?我想参考一下您编写的正确代码,希望您能将代码中有误的地方指出。十分感谢!
  • 打赏
  • 举报
回复
闫文慧 2020-02-22
引用 3 楼 Qisphs 的回复:
[quote=引用 2 楼 Qisphs 的回复:]
我也测试了一下代码,应该就是X1 和 X2 的那个问题


兄弟,你好像理解错了数据同化的意思[/quote]您好,冒昧打扰,请教一个问题,我完全按照上面的程序和文件,但运行几秒就停止了,也不显示结果,我不知道是什么原因。希望收到您的回复。十分感谢。
  • 打赏
  • 举报
回复
闫文慧 2020-02-22
引用 2 楼 Qisphs 的回复:
我也测试了一下代码,应该就是X1 和 X2 的那个问题
您好,冒昧打扰,请教一个问题,我完全按照上面的程序和文件,但运行几秒就停止了,也不显示结果,我不知道是什么原因。希望收到您的回复。十分感谢。
  • 打赏
  • 举报
回复
闫文慧 2020-02-22
引用 1 楼 weixin_41443838 的回复:
x1 = data(:,2); %背景场降水数据 x2 = data(:,3); %观测场降水数据 是不是这个问题呀
您好,冒昧打扰,请教一个问题,我完全按照上面的程序和文件,但运行几秒就停止了,也不显示结果,我不知道是什么原因。希望收到您的回复。十分感谢。
  • 打赏
  • 举报
回复
Qisphs 2019-12-28
引用 2 楼 Qisphs 的回复:
我也测试了一下代码,应该就是X1 和 X2 的那个问题


兄弟,你好像理解错了数据同化的意思
  • 打赏
  • 举报
回复
Qisphs 2019-12-26
我也测试了一下代码,应该就是X1 和 X2 的那个问题
  • 打赏
  • 举报
回复
weixin_41443838 2019-12-21
x1 = data(:,2); %背景场降水数据 x2 = data(:,3); %观测场降水数据 是不是这个问题呀
  • 打赏
  • 举报
回复
@weixin_41443838 同问
  • 举报
回复
相关推荐
发帖
其他

238

社区成员

其他产品/厂家
社区管理员
  • 其他
加入社区
帖子事件
创建了帖子
2019-04-18 04:48
社区公告
暂无公告