我用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
uj5u.com熱心網友回復：
x1 = data(:,2); %背景場降水資料x2 = data(:,3); %觀測場降水資料
是不是這個問題呀
uj5u.com熱心網友回復：
我也測驗了一下代碼，應該就是X1 和 X2 的那個問題uj5u.com熱心網友回復：
，兄弟，你好像理解錯了資料同化的意思
uj5u.com熱心網友回復：
您好，冒昧打擾，請教一個問題，我完全按照上面的程式和檔案，但運行幾秒就停止了，也不顯示結果，我不知道是什么原因。希望收到您的回復。十分感謝。uj5u.com熱心網友回復：
您好，冒昧打擾，請教一個問題，我完全按照上面的程式和檔案，但運行幾秒就停止了，也不顯示結果，我不知道是什么原因。希望收到您的回復。十分感謝。uj5u.com熱心網友回復：
您好，冒昧打擾，請教一個問題，我完全按照上面的程式和檔案，但運行幾秒就停止了，也不顯示結果，我不知道是什么原因。希望收到您的回復。十分感謝。uj5u.com熱心網友回復：
您好，請問您已經解決結果接近觀測值的問題了嗎？我想參考一下您撰寫的正確代碼，希望您能將代碼中有誤的地方指出。十分感謝！uj5u.com熱心網友回復：
您好，我運行了您的代碼認為問題應該出在“對背景場集合及觀測場集合添加擾動”這一節。觀測場不是通過觀測值得到的，應該要用Y = H*Xb+RandomGC計算，最后結果就是很接近背景場了。轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/14622.html
標籤：其他
上一篇：日本IT派遣狀況（東京）