clear all;close all;clc
%% 读取数据
lat=ncread('NCEP.sst.mean.nc','lat');
lon=ncread('NCEP.sst.mean.nc','lon');
sst=ncread('NCEP.sst.mean.nc','sst');
time=ncread('NCEP.sst.mean.nc','time');
sst(sst == -9.969209968386869e+36) = NaN;
%% 插值
[Lat Lon]=meshgrid(lat,lon);
Latq=25;Lonq=130;
for i=1:226
sst1=squeeze(sst(:,:,i));
sst2(i)=interp2(Lat,Lon,sst1,Latq,Lonq);
end
%% 去除线性趋势
dsst=detrend(sst2);
trendline=sst2-dsst;
sst2=dsst;
%% 气候态月平均sst
time1=datestr(time,30);
mm=str2num(time1(:,5:6));
for j=1:12
mm_index=find(mm==j);
month_average(j)=mean(sst2(mm_index));%求每年同一月份的平均
circle_sst(mm_index)=month_average(j);%把月平均数赋予给每年的月份
anomaly_sst(mm_index)=sst2(mm_index)-circle_sst(mm_index);%
end
%% 一维自回归预报
n=226;
for i=6:n
X(i-5,:)=anomaly_sst(i-1:-1:i-5);
Y(i-5)=anomaly_sst(i);
end
Y=Y';
B=regress(Y,X);
for j=6:n
anomaly_sst1(j)=B(1)*anomaly_sst(j-1)+B(2)*anomaly_sst(j-2)+B(3)*anomaly_sst(j-3)+B(4)*anomaly_sst(j-4)+B(5)*anomaly_sst(j-5);
end
%% 加上趋势项
sst2=trendline+sst2;
forecast=anomaly_sst1+circle_sst+trendline;
circle_sst=circle_sst+trendline;
%% 画图
subplot(211)
set(gcf,'color','w')
plot(time,sst2)
hold on
plot(time,circle_sst,'k-','linewidth',2)
plot(time,forecast,'linewidth',2)
xlabel('year')
ylabel('sst(\circC)')
datetick('x','yyyymm')
legend('时间序列','年周期','预报值')
subplot(212)
plot(time,anomaly_sst,'k-');
hold on
plot(time,anomaly_sst1,'r--');
xlabel('year')
ylabel('sst(\circC)')
datetick('x','yyyymm')
legend('去气候态平均','预报值')
%%
L=5;
eil=anomaly_sst1(L+1:end)-anomaly_sst(L+1:end);
lamdal=sqrt(sum(eil.^2)/length(eil));
lamdat=std(anomaly_sst);
g=lamdal/lamdat
matlab 准平稳时间序列分析预报方法
最新推荐文章于 2023-12-17 11:02:53 发布