tibetrain.zip_ATT6_ZNGR_anythingqkz_降水_青藏高原极端降水的线性趋势

%青藏高原站点降水的线性拟合 a=textread('D:\西藏整编降水资料\1\51886.txt','%f'); a1=reshape(a,[15,19480])'; a2=a1(:,[1 5:7 9 13]);%读取指定列的数据，分别为区站号 年 月 日 平均气温 日降水 n=size(a2(:,1)); a4=[]; j=1; for i=1:n(1) if a2(i,2)>1960 && a2(i,6)>=1 && a2(i,6)<2000 a3(1,:)=a2(i,2:6); a4=[a4;a3]; j=j+1; end end %将相同年份的资料存储在同一个数组中 k=1961; c=cell(1,50);%建立多元数组 for k=1961:2011 b1=find(a4(:,1)==k);%查找每年的数据的位置 i=b1(1);%开始的位置 j=b1(end);%结束的位置 c{k-1960}=a4(i:j,:);%将每年的数据赋值给不同的数组 end %将每年的数据升序排列，提取每年降水的最大值,每年降水前十的累计值 c1=cell(51,1); c2=[]; for i=1:51 c1{i,1}=sortrows(c{1,i},-5); c2(i,1)=c1{i,1}(1,1); c2(i,2)=c1{i,1}(1,5);%top_1极端降水 c2(i,3)=sum(c1{i,1}(1:10,5));%top_10极端降水 end %x线性拟合 timeSeries=c2(:,1); dataSeries=c2(:,2); dataCount=length(dataSeries); %% 正序列计算 SK=zeros(size(dataSeries)); % 定义SK UFK=zeros(size(dataSeries)); % 定义UFK s=0; for i=2:dataCount for j=1:i if dataSeries(i)>dataSeries(j) s=s+1; else s=s+0; end end SK(i)=s; E=i*(i-1)/4; % SK均值 Var=i*(i-1)*(2*i+5)/72; % SK方差 UFK(i)=(SK(i)-E)/sqrt(Var); end %% 逆序列计算 oppSK=zeros(size(dataSeries)); UBK=zeros(size(dataSeries)); s=0; opphumiditySeries=flipud(dataSeries); for i=2:dataCount for j=1:i if opphumiditySeries(i)>opphumiditySeries(j) s=s+1; else s=s+0; end end oppSK(i)=s; E=i*(i-1)/4; % SK均值 Var=i*(i-1)*(2*i+5)/72; % SK方差 UBK(i)=0-(oppSK(i)-E)/sqrt(Var); end oppUBK=flipud(UBK); xlswrite(outputdata,UFK,'Sheet2','A1');% 输出结果 xlswrite(outputdata,oppUBK,'Sheet2','B1'); %% 制图 figure(3); plot(timeSeries,UFK,'r-','linewidth',1.5); hold on plot(timeSeries,oppUBK,'b-.','linewidth',1.5); plot(timeSeries,1.96*ones(dataCount,1),':','linewidth',1); axis([min(timeSeries),max(timeSeries),-5,5]); legend('UF统计量','UB统计量','0.05显著水平'); xlabel('t (year)','FontName','TimesNewRoman','FontSize',12); ylabel('统计量','FontName','TimesNewRoman','Fontsize',12); hold on plot(timeSeries,0*ones(dataCount,1),'-.','linewidth',1); plot(timeSeries,1.96*ones(dataCount,1),':','linewidth',1); plot(timeSeries,-1.96*ones(dataCount,1),':','linewidth',1);