气象统计方法实习报告

1、%象统计方法课程实践内容20132第 页共66页第 页共66页气象统计方法实习报告目录实习一求500hPa高度场气候场.距平场和均方TOC o 1-5 h z差场51、资料介绍52.要求53、实习结果61)FORTRAN源程序6.grads文件10.实习结果11实习二计算给定数据资料的简单相关系数和自相关系数19资料介绍19要求20实习结果20.Fortran源程序20.程序运行结果：25实习三分析中国夏季降水线性趋势的分布特征26资料介绍及要求:262实习结果26(1).matlab程序26(2)程序运行结果27实习四求给定数据的一元线性回归方程29TOC o 1-5 h z资料介绍及要求2

2、9实习结果30(1).MATLAB程序30、程序运行结果31.结果分析32实习五对给定的海温数据进行EOF分析353535要求TOC o 1-5 h z实习结果:35（1）.FORTRAN源程序35（2）空间场和时间序列的ctl文件-39（3）运行结果4145（4）分析实习三（附加）计算给定数据的11年滑动平均和累积距平48资料介绍49要求4950实习结果实习四（附加）求给走数据的多元线性回归方54TOC o 1-5 h z说明54要求54实习结果:54(1)Matlab源程序54(2)运行结果66实习一求500hPa高度场气候场.距平场和均方差场有一500hPa高度场资料,文件名h500.d

3、at,范围:60150E,040N时1982.1-1985.12共48个月。水平分辨率：25*25,格点数：37*17o2.要求编fortran程序，求SOOhPa高度场的(1)气候场;(2)距平场;(3)均方差场。并能用Grads做出图形,实习报告中气候场.距平场.均方差场任意给出两张图，图注要清楚，即要注明是哪个时间的图形，并做简单分析。注：hSOO.For给出了如何用fortran读取ASCII码资料h500.dat.3.卖习结果1人FORTRAN源程序programex_gradsimplicitnoneinteger,parameter:nx=37,ny=17,nz=4,nt二12i

4、ntegeri,j,iz,itrealvar(nx,ny,nz,nt),cl(nx,ny,nt),sum,jp(nx,ny,nz,nt),jfc(nx,ny,nt)!Openingfileopen(10,file二g:gradsdatah500.dat,)doiz=l,nzdoit=l,ntread(10,1000)read(10,3000)(var(i,j,iz,it),i=l,nx),j=l,ny)enddoenddo1000format(2i7)%象统计方法课程实践内容20132%象统计方法课程实践内容20132第 页共66页第 页共66页2000format(37f6.2)3000fo

5、rmat(37f&1)4000format(37f7.2)close(10)!Outputopen(16,file=,g:gradsdatah500grd,form=,binary，)doiz=l,nzdoit=l,ntwrite(16)(var(i,j,iz,it),i=l,nx),j=l,ny)enddoenddo!CalculatingtheClimatologicalFielddoit=l,ntdoi=l,nxdoj=l,nysum=0doiz=l,nzsum=sum+var(i,j,iz,it)enddocl(i,j,it)二sum/4enddoenddoenddo!Outputcl

6、imate-fileopen(12,file=,g:gradsdataclimategrd,form=,binary)doit=l,ntwrite(12)(cl(i,j,it),i=l,nx),j=l,ny)enddoICalculatingtheAnomalydoiz=l,nzdoit=l,ntdoi=l,nxdoj=l,nyjp(i,j,iz,it)=var(i,j,iz,it)-cl(i,j,it)enddoenddoenddoenddoopen(13,file=,g:gradsdataanomaly.grd,form=,binary,)!Outputanomaly-filedoiz=l

7、,nzdoit=l,ntwrite(13)(jp(i,j,iz,it),i=l,nx),j=l,ny)enddoenddo!CalculatingtheMean-squareDeviationdoit=l,ntdoi=l,nxdoj=l,nysum=0doiz=l,nzsum=sum+(jp(i,j,iz,it)*2enddojfc(i,j,it)=sqrt(sum/4)%象统计方法课程实践内容2013统计方法课程实践内容2013第 页共66页第 页共66页enddoenddoenddo!Outputmean-squaredeviation-fileopen(14,file=,g:gradsd

8、atadeviationgrdJ,form=,binary,)doit=l,ntwrite(14)(jfc(i,j,it),i=l,nx),j=l,ny)enddoend(2)grads文件openg:gradsdata*.ctT(*为所求场对应的ctl文件名)setlat040setIon60150,setlev500enableprintg:gradsdata*gmf,(*为所求场名称）1=1while（imax_y)THENmax_y=rxy_ty(ty)k=tyENDIFENDDOPRINT“C全年平均气温绝对值最大自相关系数rxy_ty=,f7.4,/,1滞后时间长度k=I2),rx

9、y_ty(k),kk=0DOtw=l,N/2DOi=l,N-twrtw(i)=(w(i)-avr_w)/sw)：i：(w(i+tw)-avr_w)/sw)rxy_tw(tw)=rxy_tw(tw)+rtw(i)ENDDOrxy_tw(tw)=rxy_tw(tw)/(N-tw)rxy_tw(tw)=ABS(rxy_tw(tw)IF(rxy_tw(tw)max_w)THENmax_w=rxy_tw(tw)k=twENDIFENDDOPRINTC冬季平均气温绝对值最大自相关系数rxy_tw=f7.4,/；滞后时间长度k=I2)H,rxy_tw(k),kk=0!落后交叉相关系数DOtyw=l,N/2D

10、Oi=l9N-tywrtyw(i)=(y(i)-avr_y)/sy)：i：(w(i+tyw)-avr_w)/sw)rxy_tyw(tyw)=rxy_tyw(tyw)+rtyw(i)ENDDOrxy_tyw(tyw)=rxy_tyw(tyw)/(N-tyw)rxy_tyw(tyw)=ABS(rxy_tyw(tyw)IF(rxy_tyw(tyw)max_yw)THENmax_yw=rxy_tyw(tyw)k=tywENDIFENDDOPRINTC全年平均温度与冬季平均气温之间的落后交叉相关系数rxy_tyw=,f7.4,/,滞后时间长度k=,I2)rxy_tyw(k),kEND(2).程序运行结果

11、:i-G:MSDEV98MyProjects22Debugshixi2.exe*xy_tyw=0.4066相温巳自气rlu大均ti最平or)?值4季3CEt-V=对=4=to無绝k绝kurk9H温度温度度度ke舟气长温长y98坷间坷间坷间ar)T平时平时平时lnesn?19全夂希全win实习三分析中国夏季降水线性趙势的分布特征1姿蚪介绍及要求：用数据160zhan-rainfall-summer.txt,求1982-2006年中罔160北各秋夏孚阵水线性倾向率，给岀分布图并进行简单分析oread-rain.for给出了闻快竇料飽fortran程序。救据蛊丈件夾中单独给出。2卖习姑果(1).ma

12、tlab程序%编写求1982-2006年中国160站各站夏季降水线性倾向率clearallclcfid=fopen(1E:/160zhan-rainfall-summertxtTA1rtf);tline=fgets(fid);datal=fscanf(fidAT%f%f%f%f%f%28,160);data2=datal1;fclose(fid);fori=l:160;j(iA1:25)=data2(iA4:28);nl=1982:l:2006;pp(i,:)=polyfit(nl,j(i,l:25),1);endb=pP(:J)；jd=data2(:,3);wd=data2(:f2);jdc

13、=75:05:135;wdc=18:5:55;bz=griddata(jdAwdAbjdcAwdc1r1cubic=contour(jderwdcrbz)xlabel(1精度1);ylabel(1纬度1);title(1982-2006年中国160站各站夏季降水线性倾向率分布图J(2)程序运行结果%象统计方法课程实践内容2013%象统计方法课程实践内容20132第 页共66页第 页共66页1982-200拜中国160站各站真季降水线性倾向率分布图经度）实习四求给走数据的一元线性回归方程1资料介绍及要求利用下表数据，以环流指标为预报因子，气温为预报量，计算气温和Imn=lmi环流指标之间的一元线

14、性回归方程，并对回归方程进行检验。年份气温T环流指标19510.93219521.22519532.22019542.4261955-0.52719562.52419571128195802419596.21519602.71619613.2241962113019632.52219641.23019651.82419660.63319672.42619682.52019691.23219700.835答案：y=7.50.23xF=20.18F=4.41,回归方程显著、MATLAB程序%实习四求给定数据的一元线性回归方程ClimateData=xlsread(!F:气象统计方法从Excel文件

15、读取数据x=ClimateData(:,3);%提取ClimateData的第三列9即环流指标y=ClimateData(:f2);%提取ClimmteDmtm的第三列，即气温T实验数据气象统计实验四数据xls1);%xdata=ones(size(xA1)rl)Ax;%在原始数据x的左边加一列即模型包含常数项brbintrrrrintrs=regress(y,xdata)%调用regress函数作一元线性回归yhat=xdata*b;%计算y的估计值%定义元胞数组，以元胞数组形式显示系数的估计值和估计值得95%置信区间headl=1系数的估计值fJ估计值的95%置信下限，J估计值的95%置信

16、上限headl;num2cell(bAbint)%定义元胞数组，以元胞数组形式显示y的真实值、y的估计值、残差和残差的95$置信区间head2=y的真实值*Jy的估计值，残差1J残差的95$置信下限，J残差的95%置信上限&同时显示y的真实值，y的估计值、残差和残差的95%置信区间head2;num2cell(y,yhatArArint)$定义元胞数组，以元胞数组形式显示判定系数、F统计量的观测值、检验的P值和误差方差的估计值head3=1判定系数，5统计量的观测值，检验的P值，J误差方差的估计值，；head3;num2cell(s)(2)程序运行结果ans=系数的估计值7.5095-0.23

17、43估计值的95%3信下限估计值的95%置信上限4.655410.3637-0.34330.1253ans=y的真实值.y的估计值.残差0.90000.01230.88771.20001.6523-0.45232.20002.8237-0.62372.40001.41800.9820-0.50001.1837-1.6837残差的95%置信下限.黴差的95%置信上限-1.56333.3388-3.03312.1285-3.11711.8696-1.56113.5251-4.12500.7576ans=判定系数乍统计址的观测值检验的P值误差方差的估计值0.531320.40452.6673e-00

18、41.51342.5000-1.100006.20002.70003.2000-1.10002.50001.20001.80000.60002.40002.50001.2000-0.80001.88660.94941.88663.99523.76091.88660.48082.35520.48081.8866-0.22201.41802.82370.0123-0.69060.6134-2.0494-1.88662.2048-1.06091.3134-1.58080.14480.7192-0.08660.82200.9820-0.32371.1877-0.1094-1.9531-4.4072-4

19、.28670.1971-3.3535-1.1840-3.9959-2.4137-1.8001-2.6717-1.5996-1.5611-2.8318-1.2301-2.47943.18000.30840.51364.21251.23173.81080.83422.70343.23852.49863.24373.52512.18443.60562.2606(3).结果分析从输出的结果看，常数项和回归系数的估计值分别为7.5095和-02343,从而可以写出线性回归方程为$=7.50950.2343兀归系数估计值的置信区间为-0.3433,-0.1253。对回归直线进行显著性检验，原假设和对立假设

20、分别为矶：“=0,耳：“H0检验P的值为2.6673X104FORTRAN源程序!preparedataforeofanalysis!theprogramistonormalizeseasurfacetemperature(SST)!mt:thelengthoftimeseries;!mo:themonthnumbers;my:theyearnumbers;!sst:seasurfacetemperaturedata;!sst3:theworkarray;avf:theaverageofSST;!df:thevarianceofSST;programmainparameter(mo=12,my

21、=43,nx=18,ny=12,mt=516)dimensionavf(mo,nx,ny),df(mo,nx,ny)dimensionsst(nx,ny,mt),sst3(nx,ny,mt)open(l,file=,g:sstpxsstpxgrd,form=,unformatted,accessdirect/,recl=nx*ny)irec=ldoit=1,mtread(1,rec=irec)(sst(i,j,it),i=l,nx),j=l,ny)irec=irec+lenddo!averagedoj=l,nydoi=l,nxdok=l,modoit=k,mt,12avf(k,i,j)=avf

22、(k,i,j)+sst(i,j,it)enddoavf(k,i,j)=avf(k,i,j)/myenddoenddoenddo!variancedoj=l,nydoi=l,nxdok=l,modoit二k,mt,12df(k,i,j)=df(k,i,j)+(sst(i,j,it)-avf(k,i,j)*2enddodf(k,i,j)=sqrt(df(k,i,j)/my)enddo%象统计方法课程实践内容2013%象统计方法课程实践内容2013第 页共66页第 页共66页enddoenddo!standardizingdoj=l,nydoi=l,nxdok=l,modoit=k,mt,12辻(s

23、st(i,j,it)=-999.0)thensst3(i,j,it)=-999.0elsesst3(i,j,it)=(sst(i,j,it)-avf(k,i,j)/df(k,i,j)endifenddoenddoenddoenddo!outputfileopen(2,file=,g:sstpxstandardgrd,form=，unformatted,accessdirect,recl=nx*ny)irec=ldoit=l,mtwrite(2,rec=irec)(sst3(i,j,it),i=l,nx),j=1,ny)irec=irec+lenddoclose(2)close(1)end分解后

24、的时间系数写Atctgrd文件中，空间场写入evf.grd文件中，特征值和分析误差写入sstpx文件夹下的eigenvalue.dat文件，特征向量写入eigenvactors.dat文件。由eigenvalues.dat中的标准特征向量可得出一般特征值的前两个模态有效。用grads打开evf.ctl和tcf.ctl,分别画出海平面气温EOF分解后的空间场和时间序列。(2)空间场和时间序列的ctl文件evfctl%象统计方法课程实践内容20134%象统计方法课程实践内容20134第 页共66页第 页共66页dsetg:sstpxevfgrdtitleCoadsSSTAEundef-999.0

25、xdef18linear12010ydef12linear-27.55zdef1linear10001tdef2linear叮anl948lmonthvars1anomalyS099CoadsSSTinterperatedusingendvarstcf.ctldsetg:sstpxtcfgrdtitleCoadsSSTATundef-999.0 xdef1linear12010ydef1linear-27.55zdef1linear10001tdef516linear1janl948lmonthvars2a099timecoefficient1b099timecoefficient2endva

26、rs(3)运行结果第_模态空间场theSSTSpacialFieldforEOFAnalysis/NOJ时间系数theSSTTimeCoefficientforEOFAnalysis/NoJ19501955I960196519701975198019851990o-olouoxvGIAodnudn第二模态空间场theSSTSpacialFieldforEOFAnalysis/N0.225N-a15、-0.090.00-0.12M，-aas0.03-0.Q5-20S140E160WMOW80WI20WI00W二溜-012-0.09时间系数theSSTTimeCoefficientforEOFAn

27、alysis/No.21W01921501960197019751W01921930第一特征向量Eigenvalues.dat文件给出了EOF分析的第一特征向量值的216个值000-0.02200.01800.04300.03400.06400.05400.06000.06900.05500.03700.01000.0190-0.0210-0.046000000.02300.02100.05000.06100.05400.05000.02600.01300.0070-0.0600-0.0580-0.06900.0530-0.05500-0.01300-0.00100.02300.02800.0

28、3500.05100.03200.00800.0550-0.0730-0.1070-0.1180-0.0990-0.0760-0.0680-0.07800-0.0140-0.01400.00100.03600.02700.0180-0.0010-0.0560-0.0610-0.1050-0.12300.1230-0.12200.1140-0.0970-0.0870-0.111000.0130-0.00500.01700.0490-0.0140-0.0510-0.0540-0.10500.1120-0.0970-0.1010-0.1400-0.1410-0.1270-0.1270-0.1070-

29、0.11800-0.00800.01700.04900.0120-0.0640-0.0950-0.1100-0.1330-0.1250-0.1220-0.1130-0.1220-0.1270-0.1360-0.1190-0.1180000.02000.01200.0150-0.0010-0.0480-0.1040-0.1030-0.1270-0.1300-0.11600.1180-0.1130-0.1010-0.1000-0.10600.1240-0.10800-0.02800.00700.01400.0120-0.0020-0.0340-0.0710-0.0810-0.0800-0.1020

30、-0.1150-0.0980-0.0950-0.0890-0.1080-0.1380-0.11200-0.03200.01200.02500.01200.0010-0.0190-0.0080-0.0440-0.0620-0.0770-0.0810-0.0630-0.05300.0810-0.07700.1310-0.0780-0.06400.0210-0.00900.02100.03100.02400.0040-0.00900.0370-0.06100.0560-0.06400.0650-0.0720-0.0840-0.0850-0.0520-0.0660-0.0780-0.0110-0.02

31、60-0.0100-0.01100.02800.01800.02400-0.02300.0420-0.0660-0.0630-0.0650-0.09600.0160-0.0340-0.0480-0.0630-0.00800.0110-0.0070-0.00500.01300.03500.04500.06000.05500.0280-0.02300.0590-0.09500.06300.00700.00400.0100-0.0140(4)分析第_模态空间场theSSTSpociolFieldforEOFAnalysis/N0.1时间系数theSSTTimeCoefficientforEOFAna

32、lysis/No.1-201950191960196519701975IW1990此次试验EOF分析中的前两个特征向量最大限度地表征了海平面温度场的主要结构。第_特征向量所描绘的第一经验正交函数的特征场（即第一模态）具有海表面气温516个样本的最相似的特征。若其可以解释为516个月的标准化特征，它指示出海表温度变化的扰动。其对应的时间系数可以表示为第一模态空间场的时间权重。从第一模态的空间特征场可以看出，受到大尺度环流影响，整场的空间变化基本全为负值。而其值乘以时间权重后均变为负值。也就是大的时间系数乘以空间特征值对应海表温度的低值，而小的时间系数乘以空间特征值则对应高值。海表温度的低值对应了

33、气象上的拉尼娜年，而海表温度的高值对应了厄尔尼诺年。厄尔尼诺现象泛指赤道附近的东部太平洋表层海水温度上升引起的气候异常现象，它是热带海洋洋流与大气互作用的产物。其基本特征是太平洋沿岸的海面水温异常升高，海水水位上涨，并形成一股暖流向南流动。它使原属冷水域的太平洋东部水域变成暖水域，结果引起海啸和暴风骤雨，造成一些地区干旱，另一些地区又降雨过多的异常气候现象。所以，在空间特征场乘以时间系数后的高值表示厄尔尼诺年。拉尼娜现象是指海洋中的赤道的中部和东部太平洋，东西上万公里，南北跨度上千公里的范围内，海洋温度比正常温度东部和中部海面温度偏低0.2摄氏度，并持续半年（与厄尔尼诺现象正好相反），东南信风

34、将表面被太阳晒热的海水吹向太平洋西部，致使西部比东部海平面增高将近60厘米，西部海水温度增高，气压下降,潮湿空气积累形成台风和热带风暴，东部底层海水上翻，致使东太平洋海水变冷的现象。所以，在空间特征场乘以时间系数后的低值表示拉尼实习三（附加）计算给定数据的11年滑动平均和累积距平1脊料介绍利用ma.dat,编写11点滑动平均的程序，ma.for给出了阅读资料的fortran程序。数据在文件夹中单独给出。2.要求实习报告中附出程序，并给出原数据和滑动后数据的图形（1张图）Matlab程序load1g:MADAWx=MA;year=1922:l:20061;year2=year(1+(ih-1)/

35、2:length(x)-(ih-l)/2)；ih=ll;fori=l:length(x)-10avr(i)=sum(x(i:i+10)/ih;endplot(yearAxA)holdonplot(year2ravrA1r1)save(1Exam_4_output_datatxt1r1avr1rascii1)、FORTRAN程序滑动平均计算值（已导入文件Exam_4outputfiIe_DATA.dat)3.0727242.7727233.0363672.7454582.5363643.0454552.7636343.0727252.7545472.5727272.9909162.7272773

36、.1181802.7272732.6090922.9545432.6818113.0909062.7272782.6545492.9181832.6636363.1090922.7545472.6818112.9363662.6636383.1090942.7363672.7545412.9000042.6181853.1454502.7727222.7909052.8545482.5727293.1272762.7727272.8818132.7818182.5363683.1090992.7636343.0090912.7363672.5272713.081819913.054543.16

37、3632.972726773.081813.136362.990909493.118183.100003.054542163.136363.063633.072724773.181813.018183.045459253.154543.07272673.200003.05454163.172723.08181883.218183.00909313.227273.00000303.181813.00909、原始数据及实验滑动平均拟合曲线图表（图表由matlab画出）2010McnngAxvragainrti合|Datallpor初ovAvg、分析图中滑动平均值很好的描述了变量x随时间的变化趋势。

38、滑动平均值滤掉了较大的震荡，使趋势更加明显。从图中看出，x的整体趋势被体现而震荡的极大值却在趋势线中没有显著表示。所以，滑动平均在分析数据时可以更好的体现变化趋势但无法体现较大的异常值。在分析异常时，不宜使用滑动平均。%象统计方法课程实践内容20134%象统计方法课程实践内容20134%象统计方法课程实践内容2013第 页共66页第 页共66页第 页共66页实习四（附加）求给定数据的多元线性回归方程!1!为四个预报因子，丿为预报量；样本个数“13选取预报因子1、2、4,求预报量的标准化回归方程。112345678910111213X17111117113122111110X2262956315

39、25571315447406668X3615886917221842398兀4605220473322644222634121278.574.3104.387.695.9109.2102.772.593.1115.983.8113.3109.4答案：标准化变量回归方程：$=0.5679再+0.4323兀2-0.2613“附加：利用上表资料，尝试编写逐步回归程序。3.卖习结果：(1)Matlab源程序ClimateData=xlsread(F:-O(D-OCD伍气象统计方法实验数据气象统计实验四数据（附加题）xl）；UExuel文件读取数据X=ClimateData(2:5,:);%提取Clim

40、ateData的第二到5行数据，即自变量观测矩阵Xy=ClimateData(6,:);%提取ClimmteDmtm的第六行,即预报量reglm(yAX)%调用reglm函数做4重线性回归，显示回归分析的方差分析表和参数估计表functionreglm(ymodel丫varnames)%多重线性回归分析或广义线性回归分析%reglm(yAX),产生线性回归分析的方差分析表和参数估计结果，并以表格形式显示在屏幕上.参%数X是自变量观测值矩阵，它是n行P列的矩阵.y是因变量观测值向量，它是n行1列的列向量.%stats=reglm(yAX)9还返回一个包括了回归分析的所有诊断统计量的结构体变量st

41、ats=reglm(yAXAmodelAvarnames),用可选的varnames参数指定变量标签.vmmmmes%可以是字符矩阵或字符串元胞数组，它的每行的字符或每个元胞的字符串是一个变量的标签，它的行%数或元胞数应与X的列数相同.默认情况stats%stats=reglm(yAXAmodel),用可选的model参数来控制回归模型的类型.model是一个字符串，%其可用的字符串如下带有常数项%linear的线性模型(默认情况)带有常数%interaction*项、线性项和交叉项的模型带有常数%quadratic项、线性项、交叉项和平方项的模型带有常数%1purequadratic*项、线

42、性项和平方项的模型%下，用X1.X2,-作为变量标签.%例：%x=215250180250180215180215250215215%136513651365138513951385140.5140.5140.5138513851;%y=6.27.54.85.14.64.62.83.14.34.94.11;%reglm(yAxA*quadratic*)%分析表%方差来源自由度方和均方值p值平F5.000015.02773.0055761220.0219%残差5.00001974203948%总计10.000017.0018%均方根误差(RootMSE)0.6284判定系数(R-Square)0

43、8839%因变量均值(DependentMean)4.7273调整的判定系数(AdjR-Sq)07678%参数估计变量估计值130.94280.0154标准误%常数项2011.1109883TOC o 1-5 h z%XI0640503218%X229.153709779%X1*X20.004402461%X1*X10.000304400%X2*X20105509626%070401.0992-0.8487-0.0291-0.0058-1.31320.0003083840.00520.0492%Copyright2009-2010 xiezhh%$Revision:1000$Date:2009

44、/12/2221:41:00$ifnargin2，至少需要两个输入参数,endp=size(X,2);数%X的列数，即变量个ifnargin3|isempty(model)model=1linear1;%model参数的默认值end%生成变量标签varnamesifnargin4|isempty(varnames)varnamel=streat(1X1rnum2str(1:p1);varnames=makevarnames(varnamelAmodel);%默认的变量标签ifischar(varnames)varnamelcellstr(varnames);elseifiscell(varnames)varnamel=varnames(:);error(fvarnames必须是字符矩阵或字符串元胞数组T;endifsize(varnamelf1)=perr