气象统计方法实习报告

1 气象统计方法课程实践内容2013

目录

实习一 求500hPa高度场气候场、距平场和均方差场 ------------------------------------------------- 3

1、资料介绍 ------------------------------------------------------------------------------------------------- 3 2．要求 -------------------------------------------------------------------------------------------------------- 3 3、实习结果 ------------------------------------------------------------------------------------------------- 3 1）、FORTRAN源程序 ---------------------------------------------------------------------------------- 3

(2)、grads文件 --------------------------------------------------------------- 5 （3）、实习结果------------------------------------------------------------------------------------- 6

实习二 计算给定数据资料的简单相关系数和自相关系数 ------------------------------------------- 11

1、资料介绍 ------------------------------------------------------------------------------------------------ 11 2、要求 ------------------------------------------------------------------------------------------------------- 11 3、实习结果 ------------------------------------------------------------------------------------------------ 12

（1）、Fortran源程序 ----------------------------------------------------------------------------- 12 （2）、程序运行结果： -------------------------------------------------------------------------- 14

实习三 分析中国夏季降水线性趋势的分布特征 -------------------------------------------------------- 15

1.资料介绍及要求： -------------------------------------------------------------------------------------- 15 2.实习结果 ------------------------------------------------------------------------------------------------- 15

（1）.matlab程序 ----------------------------------------------------------------------------------- 15 （2）.程序运行结果 -------------------------------------------------------------------------------- 16

实习四 求给定数据的一元线性回归方程 ----------------------------------------------------------------- 17

1、资料介绍及要求 --------------------------------------------------------------------------------------- 17 2、实习结果 ------------------------------------------------------------------------------------------------ 18

第 1 页 共 36 页

2 气象统计方法课程实践内容2013

（1）、MATLAB程序 ---------------------------------------------------------------------------- 18 （2）、程序运行结果 ------------------------------------------------------------------------------ 18 （3）、结果分析 ------------------------------------------------------------------------------------ 19

实习五 对给定的海温数据进行EOF分析 ---------------------------------------------------------------- 21

1、资料介绍 ------------------------------------------------------------------------------------------------ 21 2、要求 ------------------------------------------------------------------------------------------------------- 21 3、实习结果： --------------------------------------------------------------------------------------------- 21

（1）、FORTRAN源程序 ------------------------------------------------------------------------ 21 （2）空间场和时间序列的ctl文件 ------------------------------------------------------------ 23 （3）运行结果 --------------------------------------------------------------------------------------- 24 （4）分析 --------------------------------------------------------------------------------------------- 26

实习三（附加） 计算给定数据的11年滑动平均和累积距平 -------------------------------------- 28

1、资料介绍 ------------------------------------------------------------------------------------------------ 28 2、要求 ------------------------------------------------------------------------------------------------------- 28 3、实习结果 ------------------------------------------------------------------------------------------------ 28 实习四（附加） 求给定数据的多元线性回归方程 ----------------------------------------------------- 31

1、说明 ------------------------------------------------------------------------------------------------------- 31 2、要求 ------------------------------------------------------------------------------------------------------- 31 3、实习结果： --------------------------------------------------------------------------------------------- 31

（1）Matlab源程序 -------------------------------------------------------------------------------- 31 （2）运行结果 --------------------------------------------------------------------------------------- 35

第 2 页 共 36 页

3 气象统计方法课程实践内容2013

实习一 求500hPa高度场气候场、距平场和均方差场

1、资料介绍

有一500hPa高度场资料，文件名h500.dat，范围：60~150E，0~40N. 时段：1982.1~1985.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_grads

implicit none

第 3 页 共 36 页

4 气象统计方法课程实践内容2013

integer,parameter::nx=37,ny=17,nz=4,nt=12 integer i,j,iz,it

real 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:\\gradsdata\\h500.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) enddo enddo

1000 format(2i7) 2000 format(37f6.2) 3000 format(37f8.1) 4000 format(37f7.2) close(10) !Output

open(16,file='g:\\gradsdata\\h500.grd',form='binary') do iz=1,nz do it=1,nt

write(16) ((var(i,j,iz,it),i=1,nx),j=1,ny) enddo enddo

!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/4 enddo enddo

enddo

! Output climate-file

open(12,file='g:\\gradsdata\\climate.grd',form='binary') do it=1,nt

write(12) ((cl(i,j,it),i=1,nx),j=1,ny) enddo

第 4 页 共 36 页

5 气象统计方法课程实践内容2013

!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) enddo enddo enddo enddo

open(13,file='g:\\gradsdata\\anomaly.grd',form='binary')

!Output anomaly-file do iz=1,nz do it=1,nt

write(13) ((jp(i,j,iz,it),i=1,nx),j=1,ny) enddo enddo

!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 enddo

enddo

!Output mean-square deviation-file

open(14,file='g:\\gradsdata\\deviation.grd',form='binary') do it=1,nt

write(14) ((jfc(i,j,it),i=1,nx),j=1,ny) enddo end

(2)、grads文件

'open g:\\gradsdata\\***.ctl' (***为所求场对应的ctl文件名) 'set lat 0 40'

第 5 页 共 36 页

6 气象统计方法课程实践内容2013

'set lon 60 150'

'set lev 500'

'enable print g:\\gradsdata\\***.gmf' （***为所求场名称） i=1

while(i<=48(或 12)) 'set t 'i 'd h' 'print' 'c' i=i+1 endwhile

'disable print' ;

（3）、实习结果

①、原始场 1982年1月

1982年7月

第 6 页 共 36 页

7 气象统计方法课程实践内容2013

结果分析：

冬季（此处以1月为代表）等高线分布整体平缓，表明高度场分布相对均匀，且北部接近极地位势高度低，赤道地区位势高度高，这与太阳直射点在1月在南半球，使北半球整体的辐射吸收随纬度增加而减小有关。北半球的气块受热随纬度递增而递减，因而膨胀率递减，故南方接近赤道地区的气体膨胀大，位势高，而北方近极地地区气体膨胀小，位势低。

夏季（此处以7月为代表），大洋上空出现副热带高压（588线位置），在东亚地区存在较为明显的位势高度槽，即东亚大槽。夏季在30°N以北的地区位于西风带中，从图中可看出明显的西风带长波特征。除东亚大槽外，在中亚地区也存在一长波槽，这些槽线发生长波调整时，会在部分地区发生较剧烈的天气变化。此外，从图中可以看出，1982年7月副热带高压脊线的平均位置位于25°N ,125°E附近，我国华东地区位于副高北侧西南气流控制下，西南气流为水汽输送的主要通道，该地区发生降水较为频繁。

从图中还可以看出冬夏季的位势高度分布存在明显差异，这与太阳直射点的年纪变化密切相关。

②、气候场 2月

第 7 页 共 36 页

8 气象统计方法课程实践内容2013

8月

结果分析：

气候场是多年数据中同时段的平均值序列，表征了区域内多年平均的位势高度变化。 从气候场图形可以看出多年平均的500hPa高度场中，冬季（此处以多年平均的2月气候场为代表）等高线较平直，大陆上等高线较稀疏，而海洋上等高线较密集，这表明大陆上空的位势高度变化率较海洋小。此外，冬季等位势高度线分布平直，还说明冬季的天气变化显著。

夏季（此处以多年平均的8月气候场为代表）在太平洋上有副热带高压，副高是深厚的系统，所以可以一直延伸到500hPa高度。东亚地区存在明显的西风带长波槽，即东亚大槽，

第 8 页 共 36 页

9 气象统计方法课程实践内容2013

东亚大槽的移动和变化配合副热带高压和夏季风的影响，会使我国大部分地区天气造成强烈

变化，如形成大范围降水或强对流天气。同时，孟加拉湾处存在一低涡。由于高原的阻挡作用，这一系统对我国的影响并不显著。图中副高脊线8月的多年平均位置位于30°N以北，且东部长波槽位于110°E附近，故由气候场分析，华北地区位于长波槽前，又槽线受到副高的阻挡作用，因而华北地区容易形成降水。 ③、距平场1982年1月

1982年7月

结果分析：

距平场指示了位势高度的震荡趋势，因距平的平均值为0，则大于0的值表明位势高度偏高，小于0的值表示位势高度偏低。

从图中可看出冬季（此处以1982年1月距平场做代表）在大陆位势高度为正距平，而在大洋则存在明显的负距平。则由距平场的性质得，冬季在大洋上位势高度偏高，在大陆上位势高度偏低。其原因是，海水的比热大于陆地，则冬季海洋温度比陆地高，所以海洋上气块膨胀更多，位势高度更高。

夏季（此处以1982年7月距平场做代表）相反，大洋上位势高度偏低，而陆地偏高。在70°E~90°E附近位势偏高的原因在于夏季青藏高原接受太阳辐射使之在对流层中层形成热源，位势高度因比大陆其他位置高。由此可见，位势高度的变化不仅与海陆差异有关，同时

第 9 页 共 36 页

10 气象统计方法课程实践内容2013

也与地形有关。在海洋上副热带高压所在的位置存在证据平值。

④、均方差场 6月

12月

结果分析：

均方差场反应同一时段内的位势高度变化幅度的大小。

由图可以看出，整体位势高度在大陆上的变化幅度比海洋小， 且海洋上冬季的变化幅度比夏季大，而陆地上相反，冬季的变化幅度比夏季小。因为陆地的比热小，所以陆地在夏季白天与夜间的温差大于冬季，对应的高度场震荡就比冬季剧烈。海洋上的位势高度变化幅度的影响方面温差为次要因素，其主要受到副热带高压，西风带长波槽脊影响，位势高度根据天气系统的移动而变化，所以震荡幅度较大。除受到天气系统影响外，海洋上的位势高度场还受到洋流等因素的影响。

第 10 页 共 36 页

11 气象统计方法课程实践内容2013

实习二 计算给定数据资料的简单相关系数和自相关系数

1、资料介绍

根据下表中年平均气温和冬季平均气温的等级数据进行下列计算： 1）计算两个气温之间的简单相关系数。

2）分别找出两个气温数据自相关系数绝对值最大的滞后时间长度。（滞后长度τ最大取10）

2、要求

实习报告中附出简单相关系数或自相关系数程序。

答案：r=0.47

年平均气温在滞后长度j=3、冬季序列在j=4最大。

第 11 页 共 36 页

12 气象统计方法课程实践内容2013

3、实习结果

（1）、Fortran源程序

PROGRAM EXAM

IMPLICIT NONE

INTEGER,PARAMETER::N=20 INTEGER i,j,k,ty,tw,tyw

REAL::avr_y=0,avr_w=0,sy=0,sw=0,rxy=0,max_y=0,max_w=0,max_yw=0 REAL y(N),w(N)

DATA y/3.4,3.3,3.2,2.9,3.4,2.8,3.6,3.0,2.8,3.0,3.1,3.0,2.9,2.7,3.5,3.2,3.1,2.8,2.9,2.9/ DATA

w/3.24,3.14,3.26,2.38,3.32,2.71,2.84,3.94,2.75,1.83,2.80,2.81,2.63,3.20,3.60,3.40,3.07,1.87,2.63,2.47/

REAL syy(N),sww(N),r(N),rty(N),rtw(N),rtyw(N),rxy_ty(N),rxy_tw(N),rxy_tyw(N)

!求两数组平均值 DO i=1,N

avr_y=avr_y+y(i) avr_w=avr_w+w(i) END DO

avr_y=avr_y/N avr_w=avr_w/N !简单相关系数 DO j=1,N

syy(j)=(y(j)-avr_y)**2 sy=sy+syy(j)

sww(j)=(w(j)-avr_w)**2 sw=sw+sww(j) END DO

sy=sqrt(sy/N) sw=sqrt(sw/N) DO j=1,N

r(j)=((y(j)-avr_y)/sy)*((w(j)-avr_w)/sw) rxy=rxy+r(j) END DO rxy=rxy/N

PRINT \年全年平均气温与冬季平均气温的简单相关系数rxy=',f5.2)\k=0

!自相关系数 DO ty=1,N/2 DO i=1,N-ty

第 12 页 共 36 页

13 气象统计方法课程实践内容2013

rty(i)=((y(i)-avr_y)/sy)*((y(i+ty)-avr_y)/sy)

rxy_ty(ty)=rxy_ty(ty)+rty(i) END DO

rxy_ty(ty)=rxy_ty(ty)/(N-ty) rxy_ty(ty)=ABS(rxy_ty(ty)) IF(rxy_ty(ty)>max_y) THEN max_y=rxy_ty(ty) k=ty END IF END DO

PRINT \全年平均气温绝对值最大自相关系数rxy_ty=',f7.4,/,'滞后时间长度 k=',I2)\k=0

DO 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 IF END DO

PRINT \冬季平均气温绝对值最大自相关系数rxy_tw=',f7.4,/,'滞后时间长度 k=',I2)\k=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 IF END DO

PRINT \全年平均温度与冬季平均气温之间的落后交叉相关系数rxy_tyw=',f7.4,/,'滞后时间长度 k=',I2)\

第 13 页 共 36 页

14 气象统计方法课程实践内容2013

END

（2）、程序运行结果：

第 14 页 共 36 页

15 气象统计方法课程实践内容2013

实习三 分析中国夏季降水线性趋势的分布特征

1.资料介绍及要求：

利用数据160zhan-rainfall-summer.txt，编写求1982-2006年中国160站各站夏季降水线性倾向率，给出分布图，并进行简单分析。read-rain.for给出了阅读资料的fortran程序。数据在文件夹中单独给出。 2.实习结果

（1）.matlab程序

%编写求1982-2006年中国160站各站夏季降水线性倾向率

clear all clc

fid=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); end

b=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)

第 15 页 共 36 页

16 气象统计方法课程实践内容2013

xlabel('精度');ylabel('纬度');title('1982-2006年中国160站各站夏季降水

线性倾向率分布图')

（2）.程序运行结果

1982-2006年中国160站各站夏季降水线性倾向率分布图55 504540度纬35302520 8090100110120130经度)

第 16 页 共 36 页

151050-5-10

17 气象统计方法课程实践内容2013

实习四 求给定数据的一元线性回归方程

1、资料介绍及要求

利用下表数据，以环流指标为预报因子，气温为预报量，计算气温和环流指标之间的一元线性回归方程，并对回归方程进行检验。

年份 1951 1952 1953 1954 1955 1956 1957 1958 1959 1960 1961 1962 1963 1964 1965 1966 1967 1968 1969 1970

气温T 0.9 1.2 2.2 2.4 -0.5 2.5 -1.1 0 6.2 2.7 3.2 -1.1 2.5 1.2 1.8 0.6 2.4 2.5 1.2 -0.8 环流指标 32 25 20 26 27 24 28 24 15 16 24 30 22 30 24 33 26 20 32 35 答案：

? ?7.5-0.23x yF=20.18>Fα=4.41，回归方程显著

第 17 页 共 36 页

18 气象统计方法课程实践内容2013

2、实习结果

（1）、MATLAB程序

%实习四 求给定数据的一元线性回归方程

ClimateData=xlsread('F:气象统计方法\\实验数据\\气象统计 实验四 数据.xls'); %从Excel文件读取数据

x=ClimateData(:,3); %提取ClimateData的第三列，即环流指标 y=ClimateData(:,2); %提取ClimateData的第三列，即气温T

xdata=[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.1253] ans =

'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]

第 18 页 共 36 页

19 气象统计方法课程实践内容2013

[ 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.2606] ans =

'判定系数' 'F统计量的观测值' '检验的P值' '误差方差的估计值' [ 0.5313] [ 20.4045] [2.6673e-004] [ 1.5134]

（3）、结果分析

从输出的结果看，常数项和回归系数的估计值分别为7.5095和-0.2343，从而可以写

出线性回归方程为

??7.5095?0.2343x y回归系数估计值的置信区间为 [-0.3433,-0.1253]。

对回归直线进行显著性检验，原假设和对立假设分别为

H0:?1?0,H1:?1?0

-4

检验P的值为2.6673×10<0.01，可知显著性水品α=0.01下应拒绝原假设H0,可认为y(环流指数)与x（气温T）的线性关系是显著的。

原始数据散点图与回归直线图

>>plot(x,y,'k.','Markersize',15) %原始数据散点图 >>hold on

>>plot(x,yhat,'linewidth',3) %回归直线图 >>xlabel('环流指标(x)')%标注x轴 >>ylabel('气温(y)')%标注y轴

>>legend('原始散点','回归直线')%加标注框

第 19 页 共 36 页

20 气象统计方法课程实践内容2013

7654原始散点回归直线 气温(y)3210-1-2 152025环流指标(x)3035

第 20 页 共 36 页

21 气象统计方法课程实践内容2013

实习五 对给定的海温数据进行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:\\sstpx\\sstpx.grd',form='unformatted',access='direct',recl=nx*ny)

第 21 页 共 36 页

22 气象统计方法课程实践内容2013

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 do

end do ! output file

open(2,file='g:\\sstpx\\standard.grd',form='unformatted',access='direct',recl=nx*ny) irec=1

第 22 页 共 36 页

23 气象统计方法课程实践内容2013

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.ctl

dset g:\\sstpx\\evf.grd title Coads SSTA E undef -999.0

xdef 18 linear 120 10 ydef 12 linear -27.5 5 zdef 1 linear 1000 1

tdef 2 linear 1jan1948 1month vars 1

S 0 99 Coads SST anomaly interperated using endvars tcf.ctl

dset g:\\sstpx\\tcf.grd title Coads SSTA T undef -999.0

xdef 1 linear 120 10 ydef 1 linear -27.5 5 zdef 1 linear 1000 1

tdef 516 linear 1jan1948 1month vars 2

a 0 99 time coefficient 1 b 0 99 time coefficient 2 endvars

第 23 页 共 36 页

24 气象统计方法课程实践内容2013

（3）运行结果

第一模态

空间场

时间系数

第 24 页 共 36 页

25 气象统计方法课程实践内容2013

第二模态

空间场

时间系数

第一特征向量

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

第 25 页 共 36 页

26 气象统计方法课程实践内容2013

-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）分析

第一模态

空间场

第 26 页 共 36 页

27 气象统计方法课程实践内容2013

时间系数

此次试验EOF分析中的前两个特征向量最大限度地表征了海平面温度场的主要结构。

第一特征向量所描绘的第一经验正交函数的特征场（即第一模态）具有海表面气温516个样本的最相似的特征。若其可以解释为516个月的标准化特征，它指示出海表温度变化的扰动。其对应的时间系数可以表示为第一模态空间场的时间权重。

从第一模态的空间特征场可以看出，受到大尺度环流影响，整场的空间变化基本全为负值。而其值乘以时间权重后均变为负值。也就是大的时间系数乘以空间特征值对应海表温度的低值，而小的时间系数乘以空间特征值则对应高值。海表温度的低值对应了气象上的拉尼娜年，而海表温度的高值对应了厄尔尼诺年。

厄尔尼诺现象泛指赤道附近的东部太平洋表层海水温度上升引起的气候异常现象，它是热带海洋洋流与大气互作用的产物。其基本特征是太平洋沿岸的海面水温异常升高，海水水位上涨，并形成一股暖流向南流动。它使原属冷水域的太平洋东部水域变成暖水域，结果引起海啸和暴风骤雨，造成一些地区干旱，另一些地区又降雨过多的异常气候现象。所以，在空间特征场乘以时间系数后的高值表示厄尔尼诺年。

拉尼娜现象是指海洋中的赤道的中部和东部太平洋，东西上万公里，南北跨度上千公里的范围内，海洋温度比正常温度东部和中部海面温度偏低0.2摄氏度，并持续半年（与厄尔尼诺现象正好相反），东南信风将表面被太阳晒热的海水吹向太平洋西部，致使西部比东部海平面增高将近60厘米，西部海水温度增高，气压下降，潮湿空气积累形成台风和热带风暴，东部底层海水上翻，致使东太平洋海水变冷的现象。所以，在空间特征场乘以时间系

数后的低值表示拉尼娜年。

第 27 页 共 36 页

28 气象统计方法课程实践内容2013

实习三（附加） 计算给定数据的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） 3.072727 2.936364 2.745455 2.736364 2.763637 3.045455 2.900000 2.754546 2.772728 2.727273 2.990910 2.854546 2.727273 2.772728 2.681818 2.954546 2.781819 2.727273 2.763637 2.663637 2.918182 2.736364 2.754546 2.772728 2.663637

第 28 页 共 36 页

29 气象统计方法课程实践内容2013

2.618182

2.572727 2.536364 2.527273 2.536364 2.572727 3.009091 2.972727 2.990909 3.054546 3.072727 3.045455

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.172728 3.218183 3.227273 3.181819 3.163637 3.136364 3.100001 3.063637 3.018182 3.072727 3.054546 3.081818 3.009091 3.000000

第 29 页 共 36 页

30 气象统计方法课程实践内容2013

② 、原始数据及实验滑动平均拟合曲线图表(图表由matlab画出)

③、分析

图中滑动平均值很好的描述了变量x随时间的变化趋势。滑动平均值滤掉了较大的震荡，使趋势更加明显。从图中看出，x的整体趋势被体现而震荡的极大值却在趋势线中没有显著表示。所以，滑动平均在分析数据时可以更好的体现变化趋势但无法体现较大的异常值。在分析异常时，不宜使用滑动平均。

第 30 页 共 36 页

31 气象统计方法课程实践内容2013

实习四（附加） 求给定数据的多元线性回归方程

1、说明

x1-x4为四个预报因子，y为预报量；样本个数n=13

2、要求

选取预报因子1、2、4，求预报量的标准化回归方程。 i x1 x2 x3 x4 1 7 2 1 3 11 4 5 6 11 7 3 8 1 9 2 10 11 12 21 1 11 13 10 68 8 12 11 7 26 29 56 31 52 55 6 15 8 8 6 9 71 31 54 47 40 66 17 22 18 6 4 23 9 60 52 20 47 33 22 44 22 26 34 12 y 78.5 74.3 104.3 87.6 95.9 109.2 102.7 72.5 93.1 115.9 83.8 113.3 109.4

? =0.5679x1 +0.4323x2?0.2613x4 答案：标准化变量回归方程：y附加：利用上表资料，尝试编写逐步回归程序。

3、实习结果：

（1）Matlab源程序

ClimateData=xlsread('F:\\-Oの\\-Oの伍\\气象统计方法\\实验数据\\气象统计 实验

四 数据（附加题）.xls'); %从Excel文件读取数据

X=ClimateData(2:5,:); %提取ClimateData的第二到5行数据，即自变量观测矩阵X y=ClimateData(6,:); %提取ClimateData的第六行，即预报量 reglm(y,X) %调用reglm函数做4重线性回归，显示回归分析的方差分析表和参数估计表

第 31 页 共 36 页

32 气象统计方法课程实践内容2013

function stats = reglm(y,X,model,varnames) % 多重线性回归分析或广义线性回归分析

%

% reglm(y,X)，产生线性回归分析的方差分析表和参数估计结果，并以表格形式显示在屏幕上. 参

% 数X是自变量观测值矩阵，它是n行p列的矩阵. y是因变量观测值向量，它是n行1列的列向量. %

% stats = reglm(y,X)，还返回一个包括了回归分析的所有诊断统计量的结构体变量stats. %

% stats = reglm(y,X,model)，用可选的model参数来控制回归模型的类型. model是一个字符串， % 其可用的字符串如下

% 'linear' 带有常数项的线性模型（默认情况） % 'interaction' 带有常数项、线性项和交叉项的模型

% 'quadratic' 带有常数项、线性项、交叉项和平方项的模型 % 'purequadratic' 带有常数项、线性项和平方项的模型 %

% stats = reglm(y,X,model,varnames)，用可选的varnames参数指定变量标签. varnames

% 可以是字符矩阵或字符串元胞数组，它的每行的字符或每个元胞的字符串是一个变量的标签，它的行

% 数或元胞数应与X的列数相同. 默认情况下，用X1,X2,…作为变量标签. %

% 例:

% x = [215 250 180 250 180 215 180 215 250 215 215

% 136.5 136.5 136.5 138.5 139.5 138.5 140.5 140.5 140.5 138.5 138.5]';

% y = [6.2 7.5 4.8 5.1 4.6 4.6 2.8 3.1 4.3 4.9 4.1]'; % reglm(y,x,'quadratic') %

% ------------------------------------------------------------------------

方

差

分

析

表

% 方差来源 自由度 平方和 均方 F值 p值

% 回归 5.0000 15.0277 3.0055 7.6122 0.0219

% 残差 5.0000 1.9742 0.3948 % 总计 10.0000 17.0018 %

% 均方根误差(Root MSE) 0.6284 判定系数(R-Square) 0.8839

第 32 页 共 36 页

33 气象统计方法课程实践内容2013

% 因变量均值(Dependent Mean) 4.7273 调整的判定系数(Adj

R-Sq) 0.7678 %

% ----------------------------------------------------------------------

参

数

估

计

% 变量 估计值 标准误 t值 p值

% 常数项 30.9428 2011.1117 0.0154 0.9883

% X1 0.7040 0.6405 1.0992 0.3218

% X2 -0.8487 29.1537 -0.0291 0.9779

% X1*X2 -0.0058 0.0044 -1.3132 0.2461

% X1*X1 0.0003 0.0003 0.8384 0.4400

% X2*X2 0.0052 0.1055 0.0492 0.9626 %

% Copyright 2009 - 2010 xiezhh.

% $Revision: 1.0.0.0 $ $Date: 2009/12/22 21:41:00 $

if nargin < 2

error('至少需要两个输入参数'); end

p = size(X,2); % X的列数，即变量个数 if nargin < 3 || isempty(model)

model = 'linear'; % model参数的默认值 end

% 生成变量标签varnames

if nargin < 4 || isempty(varnames)

varname1 = strcat({'X'},num2str([1:p]'));

varnames = makevarnames(varname1,model); % 默认的变量标签 else

if ischar(varnames)

varname1 = cellstr(varnames); elseif iscell(varnames)

varname1 = varnames(:); else

error('varnames 必须是字符矩阵或字符串元胞数组'); end

if size(varname1,1) ~= p

第 33 页 共 36 页

34 气象统计方法课程实践内容2013

error('变量标签数与X的列数不一致');

else

varnames = makevarnames(varname1,model); % 指定的变量标签 end end

ST = regstats(y,X,model); % 调用regstats函数进行线性回归分析，返回结构体变量ST

f = ST.fstat; % F检验相关结果 t = ST.tstat; % t检验相关结果

% 显示方差分析表 fprintf('\\n');

fprintf('------------------------------------------------------------------------'); fprintf('\\n');

方

差

分

析

表

fprintf('%s%7sssss','方差来源','自由度','平方和','均方','F值','p值'); fprintf('\\n');

fmt = '%s.4f.4f.4f.4f.4f';

fprintf(fmt,'回归',f.dfr,f.ssr,f.ssr/f.dfr,f.f,f.pval); fprintf('\\n');

fmt = '%s.4f.4f.4f';

fprintf(fmt,'残差',f.dfe,f.sse,f.sse/f.dfe); fprintf('\\n');

fmt = '%s.4f.4f';

fprintf(fmt,'总计',f.dfe+f.dfr,f.sse+f.ssr); fprintf('\\n'); fprintf('\\n');

% 显示判定系数等统计量

fmt = '"s.4f%s.4f';

fprintf(fmt,'均方根误差(Root MSE)',sqrt(ST.mse),'判定系数(R-Square)',ST.rsquare); fprintf('\\n');

fprintf(fmt,'因变量均值(Dependent Mean)',mean(y),'调整的判定系数(Adj R-Sq)',ST.adjrsquare); fprintf('\\n'); fprintf('\\n');

% 显示参数估计及t检验相关结果

fprintf('----------------------------------------------------------------------'); fprintf('\\n');

参数估计

第 34 页 共 36 页

35 气象统计方法课程实践内容2013

fprintf('%8sssss','变量','估计值','标准误','t值','p值');

fprintf('\\n');

for i = 1:size(t.beta,1) if i == 1

fmt = '%8s .4f.4f.4f.4f\\n';

fprintf(fmt,'常数项',t.beta(i),t.se(i),t.t(i),t.pval(i)); else

fmt = 's .4f.4f.4f.4f\\n';

fprintf(fmt,varnames{i-1},t.beta(i),t.se(i),t.t(i),t.pval(i)); end end

if nargout == 1

stats = ST; % 返回一个包括了回归分析的所有诊断统计量的结构体变量 end

% -----------------子函数----------------------- function varnames = makevarnames(varname1,model) % 生成指定模型的变量标签 p = size(varname1,1); varname2 = []; for i = 1:p-1 varname2

[varname2;strcat(varname1(i),'*',varname1(i+1:end))]; end

varname3 = strcat(varname1,'*',varname1); switch model

case 'linear'

varnames = varname1; case 'interaction'

varnames = [varname1;varname2]; case 'quadratic'

varnames = [varname1;varname2;varname3]; case 'purequadratic'

varnames = [varname1;varname3]; end

（2）运行结果

有错误，未来得及仔细改

第 35 页 共 36 页

=

36 气象统计方法课程实践内容2013

第 36 页 共 36 页

本文来源：https://www.bwwdw.com/article/rfdg.html