MATLAB EOF处理 缺测值问题

最近要用到EOF处理数据，根据魏凤英老师在《现代气候统计诊断与预测技术》中提出的那样：

变量场输入的三种形式决定了分析重点，关于变量的距平场和标准化场，ncl中给出了直接方便的处理方式，但是ncl中没有给出关于原始变量场这一项的选择，导致我在想研究其平均状况时受阻（这种情况不多）。所以打算用matlab来解决，翻看一些博主的文章，发现很多都已经给出了EOF函数，但是关于气象中缺测值的问题没有提（很可能是我看的太少），所以做点补充。

%model_bias等都是自己定义的数组，进行替换即可
%若EOF区域跨越的纬度较大，例如从30-70N，需要将model_bias乘以一个系数再进行EOF，以弥补极地面积减小的变化。
%经度，纬度，模式个数（时间维），季节时次;
%经度由小到大，纬度由小到大;
model_bias=zeros(192,144,46,5);
wgt = zeros(144);
pi   = atan(1.0)*4;
for i=1:144
wgt(i)  = sqrt(cosd( lat(i)*pi/180.0 ));
model_bias(:,i,:,:) = bias(:,i,:,:)*wgt(i);%192x144x46x5
end


%因为气象数据时常有缺测，需要先剔除掉缺测数据，否则进行EOF计算会报错
tic;
num=0;
for l=1:192%经度时次
disp(l)
for k=1:144%纬度时次
disp(k)
if (ismember(1,isnan(model_bias(l,k,:,:)))==0)
num=num+1;%将二维数组转化为一维数组，便于计数
disp(num);
new_model_bias(num,:,:) = model_bias(l,k,:,:);
new_lon(num) = lon(l,1);
new_lat(num) = lat(k,1);
end
end
end
toc;%计时


%原始数据
%原始数据处理维持原状
%计算距平数据，运行下一行
%new_model_bias=new_model_bias-mean(new_model_bias,2)

%EOF计算
[v_ann,pc_ann,lambda_ann,rvc_ann]=fEOF(new_model_bias(:,:,1));
[v_mam,pc_mam,lambda_mam,rvc_mam]=fEOF(new_model_bias(:,:,2));
[v_jja,pc_jja,lambda_jja,rvc_jja]=fEOF(new_model_bias(:,:,3));
[v_son,pc_son,lambda_son,rvc_son]=fEOF(new_model_bias(:,:,4));
[v_djf,pc_djf,lambda_djf,rvc_djf]=fEOF(new_model_bias(:,:,5));


%标准化时间序列;
pc_ann = zscore(pc_ann(1,:),0)%计算样本标准差
pc_mam = zscore(pc_mam(1,:),0)
pc_jja = zscore(pc_jja(1,:),0);
pc_son = zscore(pc_son(1,:),0);
pc_djf = zscore(pc_djf(1,:),0);


%复原空间场
tic;
num=0;
for l=1:192%经度时次
disp(l)
for k=1:144%纬度时次
disp(k)
if (ismember(1,isnan(model_bias(l,k,:,:)))==0)
%将一维数组复原成空间场
num=num+1;
disp(num);
mode_ann(l,k) =v_ann(num,1);
mode_mam(l,k)= v_mam(num,1);
mode_jja(l,k)=v_jja(num,1);
mode_son(l,k)=v_son(num,1);
mode_djf(l,k)=v_djf(num,1);
else
mode_ann(l,k) =nan;
mode_mam(l,k)=nan;
mode_jja(l,k)=nan;
mode_son(l,k)=nan;
mode_djf(l,k)=nan;
end
end
end
toc;%计时


%年和四季第一模态的方差贡献率
first_var(1,1)=rvc_ann(1,1);
first_var(1,2)=rvc_mam(1,1);
first_var(1,3)=rvc_jja(1,1);
first_var(1,4)=rvc_son(1,1);
first_var(1,5)=rvc_djf(1,1);

%年和四季的第一模态的时间系数
first_pc(1,:)=pc_ann;
first_pc(2,:)=pc_mam;
first_pc(3,:)=pc_jja;
first_pc(4,:)=pc_son;
first_pc(5,:)=pc_djf;

%年和四季的第一模态
first_mode(1,:,:)=mode_ann;
first_mode(2,:,:)=mode_mam;
first_mode(3,:,:)=mode_jja;
first_mode(4,:,:)=mode_son;
first_mode(5,:,:)=mode_djf;


%这个函数是在网站上下载的，由 2019/5/3 Ke Li 编写，感谢！
function [v,pc,lambda,rvc]=fEOF(x)
% EOF analysis
% supposing the input data x has been standardized
% of course you can standardize the data by anomaly or std-anomaly
% it all depends on your study demand
% for time-space matrix, 
% be sure the row denotes the spatial dimension and column for time
%
% return varibles:
% v: spatial mode or pattern for each column (eigenvector)
% pc: time factors for each row (principle component)
% lambda: eigenvalue
% rvc: the ratio of variance contribution
%
% 2019/5/3 editd by Ke Li

[~,n]=size(x);

% covariance matrix
c=x*x'/n;

% eigenvalue and eigenvector
[v,e]=eig(c,'nobalance');   % c*v - v*e is much closer to 0, 
                            % so the 'nobalance' option produces more accurate results in this case.
lambda=diag(e);

% inverse the result
v=fliplr(v);
lambda=flipud(lambda);

% time factors (principle components)
pc=v'*x;

% the ratio of variance contribution
rvc=lambda./sum(lambda);

这样子处理的话，我们只要在进行fEOF之前对数组进行不同形式的处理就可以得到我们想要的结果。令人疑惑的是，我用MATLAB进行的距平场EOF和ncl的距平场EOF，模态和时间系数都一致，但是方差贡献有些许差异，不知道是什么原因。

若有错误，感谢大家指正！