气象统计方法实习报告.doc

36气象统计方法课程实践内容2013 目录实习一 求500hPa高度场气候场、距平场和均方差场31、资料介绍32要求33、实习结果31）、FORTRAN源程序3(2)、grads文件5（3）、实习结果6实习二 计算给定数据资料的简单相关系数和自相关系数111、资料介绍112、要求113、实习结果12（1）、Fortran源程序12（2）、程序运行结果：14实习三 分析中国夏季降水线性趋势的分布特征151.资料介绍及要求：15 2.实习结果15（1）.matlab程序15（2）.程序运行结果16实习四 求给定数据的一元线性回归方程171、资料介绍及要求172、实习结果18（1）、MATLAB程序18（2）、程序运行结果18（3）、结果分析19实习五 对给定的海温数据进行EOF分析211、资料介绍212、要求213、实习结果：21（1）、FORTRAN源程序21（2）空间场和时间序列的ctl文件23（3）运行结果24（4）分析26实习三（附加） 计算给定数据的11年滑动平均和累积距平281、资料介绍282、要求283、实习结果28实习四（附加） 求给定数据的多元线性回归方程311、说明312、要求313、实习结果：31（1）Matlab源程序31（2）运行结果35实习一 求500hPa高度场气候场、距平场和均方差场 1、资料介绍有一500hPa高度场资料，文件名h500.dat，范围：60150E，040N.时段：1982.11985.12共48个月。水平分辨率：2.5*2.5，格点数：37*17。2要求编fortran程序，求500hPa高度场的（1）气候场；（2）距平场；（3）均方差场。并能用Grads做出图形,实习报告中气候场、距平场、均方差场任意给出两张图，图注要清楚，即要注明是哪个时间的图形，并做简单分析。注：h500.For给出了如何用fortran读取ASCII码资料h500.dat.3、实习结果 1）、FORTRAN源程序program ex_gradsimplicit noneinteger,parameter:nx=37,ny=17,nz=4,nt=12integer i,j,iz,itreal var(nx,ny,nz,nt),cl(nx,ny,nt),sum,jp(nx,ny,nz,nt),jfc(nx,ny,nt)! Opening file open(10,file=g:gradsdatah500.dat) do iz=1,nz do it=1,nt read(10,1000) read(10,3000) (var(i,j,iz,it),i=1,nx),j=1,ny) enddoenddo1000 format(2i7)2000 format(37f6.2)3000 format(37f8.1)4000 format(37f7.2)close(10)!Outputopen(16,file=g:gradsdatah500.grd,form=binary) do iz=1,nz do it=1,nt write(16)(var(i,j,iz,it),i=1,nx),j=1,ny) enddoenddo!Calculating the Climatological Field do it=1,nt do i=1,nx do j=1,ny sum=0 do iz=1,nz sum=sum+var(i,j,iz,it) enddo cl(i,j,it)=sum/4enddo enddoenddo! Output climate-file open(12,file=g:gradsdataclimate.grd,form=binary) do it=1,ntwrite(12)(cl(i,j,it),i=1,nx),j=1,ny) enddo !Calculating the Anomaly do iz=1,nz do it=1,nt do i=1,nx do j=1,ny jp(i,j,iz,it)=var(i,j,iz,it)-cl(i,j,it) enddoenddo enddoenddoopen(13,file=g:gradsdataanomaly.grd,form=binary)!Output anomaly-filedo iz=1,nzdo it=1,ntwrite(13)(jp(i,j,iz,it),i=1,nx),j=1,ny) enddoenddo!Calculating the Mean-square Deviation do it=1,nt do i=1,nx do j=1,ny sum=0 do iz=1,nz sum=sum+(jp(i,j,iz,it)*2 enddo jfc(i,j,it)=sqrt(sum/4)enddo enddoenddo!Output mean-square deviation-file open(14,file=g:gradsdatadeviation.grd,form=binary) do it=1,ntwrite(14)(jfc(i,j,it),i=1,nx),j=1,ny) enddo end(2)、grads文件open g:gradsdata*.ctl (*为所求场对应的ctl文件名)set lat 0 40set lon 60 150set lev 500enable print g:gradsdata*.gmf （*为所求场名称）i=1while(imax_y) THEN max_y=rxy_ty(ty) k=ty END IFEND DOPRINT (全年平均气温绝对值最大自相关系数rxy_ty=,f7.4,/,滞后时间长度 k=,I2),rxy_ty(k),kk=0DO tw=1,N/2 DO i=1,N-tw rtw(i)=(w(i)-avr_w)/sw)*(w(i+tw)-avr_w)/sw) rxy_tw(tw)=rxy_tw(tw)+rtw(i) END DO rxy_tw(tw)=rxy_tw(tw)/(N-tw) rxy_tw(tw)=ABS(rxy_tw(tw) IF(rxy_tw(tw)max_w) THEN max_w=rxy_tw(tw) k=tw END IFEND DOPRINT (冬季平均气温绝对值最大自相关系数rxy_tw=,f7.4,/,滞后时间长度 k=,I2),rxy_tw(k),kk=0!落后交叉相关系数DO tyw=1,N/2 DO i=1,N-tyw rtyw(i)=(y(i)-avr_y)/sy)*(w(i+tyw)-avr_w)/sw) rxy_tyw(tyw)=rxy_tyw(tyw)+rtyw(i) END DO rxy_tyw(tyw)=rxy_tyw(tyw)/(N-tyw) rxy_tyw(tyw)=ABS(rxy_tyw(tyw) IF(rxy_tyw(tyw)max_yw) THEN max_yw=rxy_tyw(tyw) k=tyw END IFEND DOPRINT (全年平均温度与冬季平均气温之间的落后交叉相关系数rxy_tyw=,f7.4,/,滞后时间长度 k=,I2),rxy_tyw(k),kEND（2）、程序运行结果：实习三 分析中国夏季降水线性趋势的分布特征 1.资料介绍及要求：利用数据160zhan-rainfall-summer.txt，编写求1982-2006年中国160站各站夏季降水线性倾向率，给出分布图，并进行简单分析。read-rain.for给出了阅读资料的fortran程序。数据在文件夹中单独给出。2.实习结果（1）.matlab程序%编写求1982-2006年中国160站各站夏季降水线性倾向率clear allclcfid=fopen(E:/160zhan-rainfall-summer.txt,rt);tline=fgets(fid);data1=fscanf(fid,%f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f,28,160);data2=data1;fclose(fid);for i=1:160; j(i,1:25)=data2(i,4:28); n1=1982:1:2006; pp(i,:)=polyfit(n1,j(i,1:25),1);endb=pp(:,1);jd=data2(:,3);wd=data2(:,2);jdc=75:0.5:135;wdc=18:.5:55;bz=griddata(jd,wd,b,jdc,wdc,cubic);c=contour(jdc,wdc,bz)xlabel(精度);ylabel(纬度);title(1982-2006年中国160站各站夏季降水线性倾向率分布图)（2）.程序运行结果实习四 求给定数据的一元线性回归方程1、资料介绍及要求利用下表数据，以环流指标为预报因子，气温为预报量，计算气温和环流指标之间的一元线性回归方程，并对回归方程进行检验。年份气温T环流指标19510.93219521.22519532.22019542.4261955-0.52719562.5241957-1.128195802419596.21519602.71619613.2241962-1.13019632.52219641.23019651.82419660.63319672.42619682.52019691.2321970-0.835答案：F=20.18F=4.41，回归方程显著2、实习结果（1）、MATLAB程序%实习四 求给定数据的一元线性回归方程ClimateData=xlsread(F:气象统计方法实验数据气象统计实验四数据.xls); %从Excel文件读取数据x=ClimateData(:,3); %提取ClimateData的第三列，即环流指标y=ClimateData(:,2); %提取ClimateData的第三列，即气温Txdata=ones(size(x,1),1),x; %在原始数据x的左边加一列1，即模型包含常数项b,bint,r,rint,s=regress(y,xdata); %调用regress函数作一元线性回归yhat=xdata*b; %计算y的估计值%定义元胞数组，以元胞数组形式显示系数的估计值和估计值得95%置信区间head1=系数的估计值,估计值的95%置信下限,估计值的95%置信上限;head1;num2cell(b,bint)%定义元胞数组，以元胞数组形式显示y的真实值、y的估计值、残差和残差的95%置信区间head2=y的真实值,y的估计值,残差,残差的95%置信下限,残差的95%置信上限;%同时显示y的真实值，y的估计值、残差和残差的95%置信区间head2;num2cell(y,yhat,r,rint)%定义元胞数组，以元胞数组形式显示判定系数、F统计量的观测值、检验的P值和误差方差的估计值head3=判定系数,F统计量的观测值,检验的P值,误差方差的估计值;head3;num2cell(s)（2）、程序运行结果ans = 系数的估计值 估计值的95%置信下限 估计值的95%置信上限 7.5095 4.6554 10.3637 -0.2343 -0.3433 -0.1253ans = y的真实值 y的估计值 残差 残差的95%置信下限 残差的95%置信上限 0.9000 0.0123 0.8877 -1.5633 3.3388 1.2000 1.6523 -0.4523 -3.0331 2.1285 2.2000 2.8237 -0.6237 -3.1171 1.8696 2.4000 1.4180 0.9820 -1.5611 3.5251 -0.5000 1.1837 -1.6837 -4.1250 0.7576 2.5000 1.8866 0.6134 -1.9531 3.1800 -1.1000 0.9494 -2.0494 -4.4072 0.3084 0 1.8866 -1.8866 -4.2867 0.5136 6.2000 3.9952 2.2048 0.1971 4.2125 2.7000 3.7609 -1.0609 -3.3535 1.2317 3.2000 1.8866 1.3134 -1.1840 3.8108 -1.1000 0.4808 -1.5808 -3.9959 0.8342 2.5000 2.3552 0.1448 -2.4137 2.7034 1.2000 0.4808 0.7192 -1.8001 3.2385 1.8000 1.8866 -0.0866 -2.6717 2.4986 0.6000 -0.2220 0.8220 -1.5996 3.2437 2.4000 1.4180 0.9820 -1.5611 3.5251 2.5000 2.8237 -0.3237 -2.8318 2.1844 1.2000 0.0123 1.1877 -1.2301 3.6056 -0.8000 -0.6906 -0.1094 -2.4794 2.2606ans = 判定系数 F统计量的观测值 检验的P值 误差方差的估计值 0.5313 20.4045 2.6673e-004 1.5134（3）、结果分析 从输出的结果看，常数项和回归系数的估计值分别为7.5095和-0.2343，从而可以写出线性回归方程为 回归系数估计值的置信区间为 -0.3433,-0.1253。 对回归直线进行显著性检验，原假设和对立假设分别为 检验P的值为2.667310-4 plot(x,y,k.,Markersize,15) %原始数据散点图hold onplot(x,yhat,linewidth,3) %回归直线图xlabel(环流指标(x)%标注x轴ylabel(气温(y)%标注y轴legend(原始散点,回归直线)%加标注框实习五 对给定的海温数据进行EOF分析 1、资料介绍给出海表温度距平数据资料sstpx.grd，以及相应的数据描述文件sstpx.ctl，对其进行EOF分析，资料的时空范围可以根据sstpx.ctl获知。 数据在文件夹中单独给出，距平或者标准化距平处理后再进行EOF。Zhunsst.for给出了如何读取资料，Ssteof.for为对距平或者标准化距平处理后的资料进行EOF分析。2、要求实习报告中给出第一特征向量及其时间系数，并分析其时空特征。 3、实习结果：（1）、FORTRAN源程序! prepare data for eof analysis! the program is to normalize sea surface temperature(SST)! mt: the length of time series;! mo: the month numbers; my: the year numbers;! sst: sea surface temperature data;!sst3: the work array; avf: the average of SST;! df: the variance of SST; program main parameter(mo=12,my=43,nx=18,ny=12,mt=516) dimension avf(mo,nx,ny),df(mo,nx,ny) dimension sst(nx,ny,mt),sst3(nx,ny,mt) open(1,file=g:sstpxsstpx.grd,form=unformatted,access=direct,recl=nx*ny) irec=1 do it=1,mt read(1,rec=irec)(sst(i,j,it),i=1,nx),j=1,ny) irec=irec+1 end do! average do j=1,ny do i=1,nx do k=1,mo do it=k,mt,12 avf(k,i,j)=avf(k,i,j)+sst(i,j,it) end do avf(k,i,j)=avf(k,i,j)/my end do end do end do! variance do j=1,ny do i=1,nx do k=1,mo do it=k,mt,12 df(k,i,j)=df(k,i,j)+(sst(i,j,it)-avf(k,i,j)*2 end do df(k,i,j)=sqrt(df(k,i,j)/my) end do end do end do ! standardizing do j=1,ny do i=1,nx do k=1,mo do it=k,mt,12 if(sst(i,j,it)=-999.0)then sst3(i,j,it)=-999.0 else sst3(i,j,it)=(sst(i,j,it)-avf(k,i,j)/df(k,i,j) end if end do end do end doend do! output file open(2,file=g:sstpxstandard.grd,form=unformatted,access=direct,recl=nx*ny) irec=1 do it=1,mt write(2,rec=irec)(sst3(i,j,it),i=1,nx),j=1,ny) irec=irec+1 end do close(2) close(1)end分解后的时间系数写入tcf.grd文件中，空间场写入evf.grd文件中，特征值和分析误差写入sstpx文件夹下的eigenvalue.dat文件，特征向量写入eigenvactors.dat文件。由eigenvalues.dat 中的标准特征向量可得出一般特征值的前两个模态有效。用grads打开evf.ctl和tcf.ctl，分别画出海平面气温EOF分解后的空间场和时间序列。（2）空间场和时间序列的ctl文件evf.ctldset g:sstpxevf.grdtitle Coads SSTA Eundef -999.0xdef 18 linear 120 10ydef 12 linear -27.5 5zdef 1 linear 1000 1tdef 2 linear 1jan1948 1monthvars 1S 0 99 Coads SST anomaly interperated usingendvars tcf.ctldset g:sstpxtcf.grdtitle Coads SSTA Tundef -999.0xdef 1 linear 120 10ydef 1 linear -27.5 5zdef 1 linear 1000 1tdef 516 linear 1jan1948 1monthvars 2a 0 99 time coefficient 1b 0 99 time coefficient 2endvars （3）运行结果第一模态空间场时间系数第二模态空间场时间系数第一特征向量Eigenvalues.dat文件给出了EOF分析的第一特征向量值的216个值 0 0 0 -0.0220 0.0180 0.0430 0.0340 0.0640 0.0540 0.0600 0.0690 0.0550 0.0370 0.0100 -0.0190 -0.0210 -0.0460 0 0 0 0 -0.0230 0.0210 0.0500 0.0610 0.0540 0.0500 0.0260 0.0130 0.0070 -0.0600 -0.0580 -0.0690 -0.0530 -0.0550 0 -0.0130 0 -0.0010 0.0230 0.0280 0.0350 0.0510 0.0320 0.0080 -0.0550 -0.0730 -0.1070 -0.1180 -0.0990 -0.0760 -0.0680 -0.0780 0 -0.0140 -0.0140 0.0010 0.0360 0.0270 0.0180 -0.0010 -0.0560 -0.0610 -0.1050 -0.1230 -0.1230 -0.1220 -0.1140 -0.0970 -0.0870 -0.1110 0 -0.0130 -0.0050 0.0170 0.0490 -0.0140 -0.0510 -0.0540 -0.1050 -0.1120 -0.0970 -0.1010 -0.1400 -0.1410 -0.1270 -0.1270 -0.1070 -0.1180 0 -0.0080 0.0170 0.0490 0.0120 -0.0640 -0.0950 -0.1100 -0.1330 -0.1250 -0.1220 -0.1130 -0.1220 -0.1270 -0.1360 -0.1190 -0.1180 0 0 -0.0200 0.0120 0.0150 -0.0010 -0.0480 -0.1040 -0.1030 -0.1270 -0.1300 -0.1160 -0.1180 -0.1130 -0.1010 -0.1000 -0.1060 -0.1240 -0.1080 0 -0.0280 0.0070 0.0140 0.0120 -0.0020 -0.0340 -0.0710 -0.0810 -0.0800 -0.1020 -0.1150 -0.0980 -0.0950 -0.0890 -0.1080 -0.1380 -0.1120 0 -0.0320 0.0120 0.0250 0.0120 0.0010 -0.0190 -0.0080 -0.0440 -0.0620 -0.0770 -0.0810 -0.0630 -0.0530 -0.0810 -0.0770 -0.1310 -0.0780 -0.0640 -0.0210 -0.0090 0.0210 0.0310 0.0240 0.0040 -0.0090 -0.0370 -0.0610 -0.0560 -0.0640 -0.0650 -0.0720 -0.0840 -0.0850 -0.0520 -0.0660 -0.0780 -0.0110 -0.0260 -0.0100 -0.0110 0.0280 0.0180 0.0240 0 -0.0230 -0.0420 -0.0660 -0.0630 -0.0650 -0.0960 0.0160 -0.0340 -0.0480 -0.0630 -0.0080 -0.0110 -0.0070 -0.0050 0.0130 0.0350 0.0450 0.0600 0.0550 0.0280 -0.0230 -0.0590 -0.0950 -0.0630 0.0070 0.0040 0.0100 -0.0140（4）分析 第一模态空间场时间系数此次试验EOF分析中的前两个特征向量最大限度地表征了海平面温度场的主要结构。第一特征向量所描绘的第一经验正交函数的特征场（即第一模态）具有海表面气温516个样本的最相似的特征。若其可以解释为516个月的标准化特征，它指示出海表温度变化的扰动。其对应的时间系数可以表示为第一模态空间场的时间权重。从第一模态的空间特征场可以看出，受到大尺度环流影响，整场的空间变化基本全为负值。而其值乘以时间权重后均变为负值。也就是大的时间系数乘以空间特征值对应海表温度的低值，而小的时间系数乘以空间特征值则对应高值。海表温度的低值对应了气象上的拉尼娜年，而海表温度的高值对应了厄尔尼诺年。 厄尔尼诺现象泛指赤道附近的东部太平洋表层海水温度上升引起的气候异常现象，它是热带海洋洋流与大气互作用的产物。其基本特征是太平洋沿岸的海面水温异常升高，海水水位上涨，并形成一股暖流向南流动。它使原属冷水域的太平洋东部水域变成暖水域，结果引起海啸和暴风骤雨，造成一些地区干旱，另一些地区又降雨过多的异常气候现象。所以，在空间特征场乘以时间系数后的高值表示厄尔尼诺年。 拉尼娜现象是指海洋中的赤道的中部和东部太平洋，东西上万公里，南北跨度上千公里的范围内，海洋温度比正常温度东部和中部海面温度偏低0.2摄氏度，并持续半年（与厄尔尼诺现象正好相反），东南信风将表面被太阳晒热的海水吹向太平洋西部，致使西部比东部海平面增高将近60厘米，西部海水温度增高，气压下降，潮湿空气积累形成台风和热带风暴，东部底层海水上翻，致使东太平洋海水变冷的现象。所以，在空间特征场乘以时间系数后的低值表示拉尼娜年。实习三（附加） 计算给定数据的11年滑动平均和累积距平 1、资料介绍利用数据ma.dat，编写11点滑动平均的程序，ma.for给出了阅读资料的fortran程序。数据在文件夹中单独给出。2、要求实习报告中附出程序，并给出原数据和滑动后数据的图形（1张图）Matlab程序load g:MA.DAT x=MA; year=1922:1:2006; year2=year(1+(ih-1)/2:length(x)-(ih-1)/2); ih=11; for i=1:length(x)-10 avr(i)=sum(x(i:i+10)/ih; end plot(year,x,b:) hold on plot(year2,avr,r) save (Exam_4_output_data.txt,avr,-ascii) 3、实习结果、FORTRAN程序滑动平均计算值（已导入文件Exam_4 output file_DATA.dat） 第 36 页 共 36 页3.072727 3.045455 2.990910 2.954546 2.918182 2.936364 2.900000 2.854546 2.781819 2.736364 2.745455 2.754546 2.727273 2.727273 2.754546 2.736364 2.772728 2.772728 2.763637 2.772728 2.763637 2.727273 2.681818 2.663637 2.663637 2.618182 2.572727 2.536364 2.527273 2.536364 2.572727 2.609091 2.654546 2.681818 2.754545 2.790909 2.881818 3.009091 3.036364 3.072727 3.118182 3.090909 3.109091 3.109091 3.145455 3.127273 3.109091 3.081819 3.054546 3.081819 3.118182 3.136364 3.181819 3.154546 3.200001 3.1727