Matleb多变量EOF（MV_EOF）程序讨论

2daidai

我理解的多变量EOF分解就是将多个变量合成为一个变量然后进行EOF分解
1.理论学习部分
MV_EOF相关介绍在http://bbs.06climate.com/forum.php?mod=viewthread&tid=20871说的特别详细
下面是我参考的MV_EOF的matlab的例子
话不多少直接上图
额一贴图他就是显示字数超了看下PDF吧。。。抱歉QQ1,QQ2,QQ3是重点内容的截图
2实践操作
2.1SST和U850处理
使用的是OISST和NCEP纬向风月平均资料

1. <font size="4">clc
   
2. clear
   
3. close all;
   
4. ncid = netcdf.open('E:\study\monsoon\uwnd.mon.mean.nc','NOWRITE');
   
5. ncid = netcdf.open('E:\study\monsoon\sst.mnmean.nc','NOWRITE');
   
6. ncdisp E:\study\monsoon\uwnd.mon.mean.nc;
   
7. ncdisp E:\study\monsoon\sst.mnmean.nc;
   
8. u= ncread('E:\study\monsoon\uwnd.mon.mean.nc','uwnd');
   
9. %U 144*73*17*823  2.5*2.5    1948.1.1-2016.7.1   1000, 925, 850, 700, 600, 500, 400, 300, 250, 200, 150, 100, 70, 50, 30, 20, 10
   
10. sst=ncread('E:\study\monsoon\sst.mnmean.nc','sst');
    
11. % SST 360*180*423 1*1 1981.12.1-2017.2.1
    
12. 
    
13. %%
    
14. %拼成一个新矩阵
    
15. %取U850 1981.1.1-2016.7.1
    
16. U=squeeze(u(17:49,29:45,3,(1982-1948)*12+1:end));
    
17. U=double(U);
    
18. sst1=sst(41:121,70:111,2:416);
    
19. %将lon,lat化为一维
    
20. for i=1:1:415;
    
21.     sst2=squeeze(sst1(:,:,i));
    
22.     sst3(:,i)=reshape(sst2,1,3402);
    
23. end
    
24. %标准化
    
25. sst4=zeros(3402,415);
    
26. for i=1:12;
    
27.       sstn(:,i)=nanmean(sst3(:,i:12:end),2);
    
28. end
    
29. for j=1:3402;
    
30.     for i=1:12
    
31.         sst4(j,i:12:end)=(sst3(j,i:12:end)-sstn(j,i))/std(sst3(j,:));
    
32.     end
    
33. end
    
34. %插值
    
35. xi=40:2.5:120;
    
36. yi=-20:2.5:20;
    
37. x=40.1:1:120.1;
    
38. y=-20.5:1:20.5;
    
39. x=double(x);
    
40. y=double(y);
    
41. xi=double(xi);
    
42. yi=double(yi);
    
43. [XX,YY]=meshgrid(yi,xi);%想要变成的格式
    
44. [X,Y]=meshgrid(y,x);
    
45. for i=1:415;
    
46.     U1(:,:,i)=griddata(XX,YY,U(:,:,i),X,Y,'v4');
    
47. end;
    
48. %将lon,lat化为一维
    
49. for i=1:1:415;
    
50.     U2=squeeze(U1(:,:,i));
    
51.     U3(:,i)=reshape(U2,1,3402);
    
52. end
    
53. %标准化
    
54. U4=zeros(3402,415);
    
55. for i=1:12;
    
56.       Un(:,i)=nanmean(U3(:,i:12:end),2);
    
57. end
    
58. for j=1:3402;
    
59.     for i=1:12
    
60.         U4(j,i:12:end)=(U3(j,i:12:end)-Un(j,i))/std(j,:);
    
61.     end
    
62. end
    
63. SST_U=zeros(2,3402,415);
    
64. SST_U(1,:,:)=sst4;
    
65. SST_U(2,:,:)=U4;
    
66. save  SST_U SST_U</font>

复制代码

2.2进行MV-EOF分解

1. <font size="4"><p><font size="5">clc
   
2. clear
   
3. close all;
   
4. load SST_U.mat
   
5. 
   
6. for i=1:1:415
   
7.     SST_U1=squeeze(SST_U(:,:,i));
   
8.     SST_U2(:,i)=reshape(SST_U1,1,6804);
   
9. end
   
10. SST_U2(SST_U2<-1000000)=0;
    
11. SST_U_nan=isnan(SST_U2);
    
12. i=0;
    
13. for j=1:6804;
    
14.     if sum(SST_U_nan(j,:))==0;
    
15.         i=i+1;
    
16.         SST_U3(i,:)=SST_U2(j,:);
    
17.     end
    
18. end
    
19. X=zeros(6804,415);                             
    
20. for j=1:1:6804
    
21.     for i=1:1:12
    
22.     y(j,i)=mean(SST_U3(j,i:12:end));
    
23.     X(j,i:12:end)=SST_U3(j,i:12:end) - y(j,i);
    
24.    end
    
25. end
    
26. %EOF分解
    
27. S=X'*X;
    
28. [v,d]=eig(S);
    
29. v=fliplr(v);                                           % 矩阵作左右翻转
    
30. d=rot90(d,2);                                          % 矩阵上下翻转后再左右翻转(查看生成的对角阵是由小到大排列的~此指令可使其由大到小排列~fliplr(flipud(d))=rot90(d,2))
    
31. diagonal=diag(d);                                      %d为特征跟%diag(v,k)以向量v的元素作为矩阵X的第k条对角线元素，当k=0时，v为X的主对角线；当k>0时，v为上方第k条对角线；当k<0时，v为下方第k条对角线。
    
32. LR=X*v;                                                %LR为时间系数矩阵
    
33. for i=1:1:415;
    
34.     LR(:,i)=LR(:,i)/sqrt(diagonal(i));
    
35. end
    
36. Ym=LR'*X;%时间系数
    
37. 
    
38. sum_d=sum(diagonal);
    
39. for i=1:415;
    
40.     G(i)=diagonal(i)/sum_d;                           % G2(i)是方差贡献率
    
41. end
    
42. SST_U_L=zeros(6804,415);
    
43. SST_U_L(:,:)=NaN;
    
44. i=1;
    
45. for j=1:6804;
    
46.    if sum(SST_U_nan(j,:))==0;
    
47.       SST_U_L(j,:)=LR(i,:);%空间系数
    
48.       i=i+1;
    
49.    end
    
50. end
    
51. %将空间模态分开两个变量,只考虑第一模态
    
52. K=reshape(SST_U_L(:,1),2,3402);</font></p><p><font size="5"> KK=reshape(K(1,:),81,42);
    
53. KK2=reshape(K(2,:),81,42);
    
54. %参考PPT的分解是下面的
    
55. eSST=SST_U_L(1:3402,:);
    
56. eU=SST_U_L(3403:end);
    
57. %将经纬度分开
    
58. eeSST=reshape(eSST(:,1),81,42);
    
59. save SST_U_L SST_U_L</font>

复制代码

3，结果讨论
只画了第一模态的图
图1是按PPT方法画出来的
图2，3是我自己想的
结果和周磊老师分析的差别非常大，不知道是什么原因
希望有人可以给看看
周磊老师的图也在附件中

2daidai

忘了说我编程超级烂，没正统学过什么编程原理，所以程序看起来可能有点乱，我尽量把重要步骤做了备注。。。希望有人能给看看

高空急流

谢谢，学习了，赞一个

zwcshitoutou

学习了，哈哈

zwcshitoutou

学习了，哈哈

lgz

谢谢，学习了，赞一个

VAS111

这个不错，学习下哈！

夏夜

楼主的图贴的有点乱呢，前三个图的title都一样，怎么pattern差这么多?

zysunhan

谢谢分享，楼主好人

zh10eng

楼主讲得挺不错的