气象统计实习报告

气象统计方法实习报告

注：本人上传实习报告，并不想为大家抄袭提供方便，只是想给大家一个参考和验证，我认为每一个用心之人都会做得比我更好，本文档只起抛砖之效。

专业：大气科学

班级：2010级x班

学号：xxxxxxxxxx

姓名：xx

目录

实习一气候场、距平场和均方差场的计算 (2)

实习二相关系数及自相关系数的计算 (6)

实习三一元线性回归 (9)

实习四滑动平均 (10)

实习五多元线性回归 (12)

实习六ＥＯＦ分析 (15)

实习心得 (26)

实习一气候场、距平场和均方差场的计算

要求：有一500hPa高度场资料，文件名h500.dat，范围：60~150E，0~40N.时段：1982.1~1985.12共48个月。水平分辨率：2.5*2.5，格点数：37*17。编fortran程序，求500hPa 高度场的

（1）气候场；

（2）距平场；

（3）均方差场。

并能用Grads做出图形,实习报告中气候场、距平场、均方差场任意给出两张图，图注要清楚，即要注明是哪个时间的图形，并做简单说明。

注：h500.For给出了如何用fortran读取ASCII码资料h500.dat.

编写的fortran程序如下：

program one

implicit none

integer nx,ny,mo,yr

parameter(nx=37,ny=17,mo=12,yr=4)

real datas(nx,ny,mo,yr)

real a(nx,ny,mo),aa(nx,ny,mo,yr),bb(nx,ny,mo)

integer i,j,t,it,ii,mi,irec

open(5,file='e:\h500.dat')

do ii=1,4

do mi=1,12

read(5,1000)

read(5,3000)((datas(i,j,mi,ii),i=1,nx),j=1,ny)

end do

end do

close(5)

!气候场

do t=1,12

do j=1,ny

do i=1,nx

a(i,j,t)=0

do it=1,4

a(i,j,t)=a(i,j,t)+datas(i,j,t,it)

end do

a(i,j,t)=a(i,j,t)/4

end do

end do

end do

!求距平场

do t=1,12

do j=1,ny

do i=1,nx

aa(i,j,t,1)=0

do it=1,4

aa(i,j,t,it)=datas(i,j,t,it)-a(i,j,t)

end do

end do

end do

end do

!求均方差场

do t=1,12

do j=1,ny

do i=1,nx

bb(i,j,t)=0

do it=1,4

bb(i,j,t)=bb(i,j,t)+(datas(i,j,t,it)-a(i,j,t))**2

end do

bb(i,j,t)=sqrt(bb(i,j,t)/4)

end do

end do

end do

open(11,file='d:\jp.grd',form='unformatted',access='direct',recl=nx*ny)

irec=0

do it=1,4

do t=1,12

irec=irec+1

write(11,rec=irec)((aa(i,j,t,it),i=1,nx),j=1,ny)

end do

end do

close(11)

open(10,file='d:\qhc.grd',form='unformatted',access='direct',recl=nx*ny) irec=0

do t=1,12

irec=irec+1

write(10,rec=irec)((a(i,j,t),i=1,nx),j=1,ny)

end do

close(10)

open(12,file='d:\jfc.grd',form='unformatted',access='direct',recl=nx*ny) irec=0

do t=1,12

irec=irec+1

write(12,rec=irec)((bb(i,j,t),i=1,nx),j=1,ny)

end do

close(12)

1000format(2i7)

2000FORMAT(37F6.2)

3000format(37f8.1)

4000format(37f7.2)

end program one

生成图片的gs文件

这里只给出生成气候场的gs文件

'reinit'

'open:e\cl i.ctl'

'enable print e:\cli.gmf'

'set t1'

*'set t7'

'd h'

'print'

'disable print'

;

图11月气候场

图27月气候场

图31982年1月距平场

图41982年7月距平场

图51月均方差场

图67月均方差场

实习二相关系数及自相关系数的计算

根据所给年平均气温和冬季平均气温的等级数据（数据略）进行下列计算：年份年平均气温

1970~1979 3.40 3.30 3.20 2.90 3.40 2.80 3.60 3.00 2.80 3.00 1980~1989 3.10 3.00 2.90 2.70 3.50 3.20 3.10 2.80 2.90 2.90

年份冬季平均气温

1970~1979 3.24 3.14 3.26 2.38 3.32 2.71 2.84 3.94 2.75 1.83 1980~1989 2.80 2.81 2.63 3.20 3.60 3.40 3.07 1.87 2.63 2.47

1）计算两个气温之间的相关系数。

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

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

求相关系数r的fortran程序：

program two

implicit none

integer i

real x(20),y(20)

real sum1,sum2,av1,av2,f1,f2,f3,rxy

data(x(i),i=1,20)/3.40,3.30,3.20,2.90,3.40,2.80,3.60,3.00,2.80,3.00,3.10,3.00,2.90,2.70,3.50,3.20,3. 10,2.80,2.90,2.90/

data(y(i),i=1,20)/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/

!计算均值

sum1=0

sum2=0

do i=1,20

sum1=sum1+x(i)

sum2=sum2+y(i)

end do

av1=sum1/20

av2=sum2/20

!计算方差的相关系数r

f1=0

do i=1,20

f1=f1+(x(i)-av1)*(y(i)-av2)

end do

f2=0

f3=0

do i=1,20

f2=f2+(x(i)-av1)**2

f3=f3+(y(i)-av2)**2

end do

rxy=f1/(sqrt(f2*f3))

print*,'R=',rxy

!!compute the R above

end program two

运行程序，得到结果如下图所示：

图7相关系数R的计算结果

计算自相关系数的程序其中j取1~7：

program two_2

implicit none

integer i,j

real x(20),y(20),s1(7),s2(7),r1(7),r2(7)

real sum1,sum2,av1,av2,f1,f2,f3,rxy

data(x(i),i=1,20)/3.40,3.30,3.20,2.90,3.40,2.80,3.60,3.00,2.80,3.00,3.10,3.00,2.90,2.70,3.50,3.20,3.10,2.80,2.90,2. 90/

data(y(i),i=1,20)/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/

!open(5,file='d:\sedata.txt')

!read(5,100)x

!read(5,100)y

!do i=1,20

!print*,x(i)

!print*,y(i)

!end do

sum1=0

sum2=0

do i=1,20

sum1=sum1+x(i)

sum2=sum2+y(i)

end do

av1=sum1/20

av2=sum2/20

!print*,av1,av2

f2=0

f3=0

do i=1,20

f2=f2+(x(i)-av1)**2

f3=f3+(y(i)-av2)**2

end do

f2=f2/20

f3=f3/20

!print*,f2,f3

do j=1,7

s1(j)=0

s2(j)=0

do i=1,20-j

s1(j)=s1(j)+(x(i)-av1)*(x(i+j)-av1)

s2(j)=s2(j)+(y(i)-av2)*(y(i+j)-av2)

!print*,s1(j),s2(j)

end do

end do

do j=1,7

s1(j)=s1(j)/(20-j)

s2(j)=s2(j)/(20-j)

!print*,s1(i),s2(i)

end do

do j=1,7

r1(j)=s1(j)/f2

r2(j)=s2(j)/f3

end do

open(10,file='d:\result.txt')

!do i=1,7

write(10,'(7f10.5)')(r1(i),i=1,7)

write(10,'(7f10.5)')(r2(i),i=1,7)

!end do

!write(10,*)r2(i)

close(10)

write(*,'(7f10.5)')(r1(i),i=1,7)

write(*,'(7f10.5)')(r2(i),i=1,7)

end

运行结果如下：

图8自相关系数的计算结果

由结果可知，年平均气温在滞后长度j=7、冬季序列在j=4最大。

实习三一元线性回归

根据下列数据，编写一元线性回归程序，当x=14.5时，y的值应该为多少？

X012345678910111213141516 Y3029.128.428.12827.727.527.22726.826.526.326.125.725.324.824要求：实习报告中附出一元线性回归程序，并给出相应原数据及回归数据图形（一张图）

线性回归

程序如下：

program linear

implicit none

integer i

real x(17),y(17)

real f1,f2,b,avx,avy,b0,y1

data(x(i),i=1,17)/0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16/

data(y(i),i=1,17)/30,29.1,28.4,28.1,28,27.7,27.5,27.2,27,26.8,26.5,26.3,26.1,25.7,25.3,24.8,24/ f1=0

do i=1,17

f1=f1+x(i)*y(i)

end do

avx=sum(x)/17!求均值

avy=sum(y)/17

f2=0

do i=1,17

f2=f2+x(i)**2

end do

b=(17*f1-sum(x)*sum(y))/(17*f2-sum(x)**2)!求b

b0=avy-b*avx

y1=b*14.5+b0

print*,'b=',b

print*,'b0=',b0

print*,'y=',b0,b,'x'

!!compute the b above

print*,'当x=14.5时，y=',y1

!100format(20F5.2)

end program

程序运行结果如下：

图9线性回归系数的计算结果

图10线性回归拟合结果

实习四滑动平均

利用数据ma.dat，编写11点滑动平均的程序，ma.for给出了阅读资料的fortran程序。数据在文件夹中单独给出。

要求：实习报告中附出程序，并给出原数据和滑动后数据的图形（一张图）

滑动平均程序如下：

program ma

!dimension x(1000),x1(1000),nny1(1000)

!write(*,10)

!10format(5x,'n=?,ih=?,nyear=?')

!Read(*,*)n,ih,nyear

!**********************************************

!*n:sample size of the time series*

!*ih:moving length*

!*nyear:first year of the series*

!*x(n):oroginal time series*

!*x1(n-ih+1):moved series*

!**********************************************

integer i,n,ih,nyear

parameter(n=85,ih=11,nyear=1922)

real x(n),x1(n-ih+1)

open(2,file='d:\ma.dat')

read(2,*)(x(i),i=1,n)

close(2)

do i=1,n-ih+1

x1(i)=sum(x(i:i+ih-1))/ih

end do

print*,x1

open(10,file='d:\moveaverage.txt')

write(10,'(f10.6)')(x1(i),i=1,n-ih+1)

close(10)

end

运行结果如图：

图11平滑结果

画出的图像为

图12滑动平均后与原始数据的比较

实习五多元线性回归

求给定数据的多元线性回归方程

要求：实习报告中附出多元线性回归程序，并给出回归数据图形（一张图）

数据说明：参考教材中数据进行程序验证.

以下是就是魏凤英书中多元线性回归中例子的数据，本实习也以此为验证对象（数据见下表）。

X1203250220212212172225220217238199 195222224229216201199227195192215 203213221216228208202239204240209 234247234222202211212241220225219

X2109-6-61-272-18033-647983-85-78 97-69547-102-9186-95-1331-66 44215881-8-3333-1081348-18 138166120-16-13-48918164129

X2105959010014911510810899105105 109140133149149105109140134105149 1131491051009910010510099110105 1459010512011510590909090105

Y 4461000576496477349367398373570378 416477364342388778526395370527497 548407605269458775371550662512370 481558478606433552466688522601705

程序如下：

program duoyuangxianxinghuigui

implicit none

integer n,i,j,t

parameter(n=44)

real x1(n),x2(n),x3(n),x(3,n)

real s(3,3),sy(3),c,d,avx(3),avy,b1,b2,b3,y(n),b0,yy(n),a1,a2,a3

!数据来自己魏凤英书籍中的相关章节

data

x1/203,250,220,212,212,172,225,220,217,238,199,195,222,224,229,216,201,199,227,195,192,215,203,213,221,21 6,228,208,202,239,204,240,209,234,247,234,222,202,211,212,241,220,225,219/

data

x2/109,-6,-61,-272,-180,33,-64,79,83,-85,-78,97,-69,5,47,-102,-91,86,-95,-133,1,-66,44,2,158,81,-8,-33,33,-108,13 4,8,-18,138,166,120,-16,-13,-48,91,8,16,41,29/

data

x3/105,95,90,100,149,115,108,108,99,105,105,109,140,133,149,149,105,109,140,134,105,149,113,149,105,100,9 9,100,105,100,99,110,105,145,90,105,120,115,105,90,90,90,90,105/

data

y/446,1000,576,496,477,349,367,398,373,570,378,416,477,364,342,388,778,526,395,370,527,497,548,407,605,26 9,458,775,371,550,662,512,370,481,558,478,606,433,552,466,688,522,601,705/

avx(1)=sum(x1)/n

avx(2)=sum(x2)/n

avx(3)=sum(x3)/n

avy=sum(y)/n

x(1,1:n)=x1

x(2,1:n)=x2

x(3,1:n)=x3

do i=1,3

do j=1,3

c=0

d=0

do t=1,n

c=c+x(i,t)*x(j,t)

d=d+x(i,t)

end do

s(i,j)=c-sum(x(i,1:n))*sum(x(j,1:n))/n sy(i)=d-sum(x(i,1:n))*sum(y)/n

end do

end do

do i=1,3

sy(i)=0

do t=1,n

sy(i)=sy(i)+(x(i,t)-avx(i))*(y(t)-avy) end do

end do

a1=s(2,1)/s(1,1)

a2=s(3,1)/s(1,1)

do i=1,3

s(2,i)=s(2,i)-s(1,i)*a1

s(3,i)=s(3,i)-s(1,i)*a2

end do

sy(2)=sy(2)-sy(1)*a1

sy(3)=sy(3)-sy(1)*a2

a3=s(3,2)/s(2,2)

do i=1,3

s(3,i)=s(3,i)-s(2,i)*a3

end do

sy(3)=sy(3)-sy(2)*a3

b3=sy(3)/s(3,3)

b2=(sy(2)-s(2,3)*b3)/s(2,2)

b1=(sy(1)-s(1,3)*b3-s(1,2)*b2)/s(1,1)

b0=avy-avx(1)*b1-avx(2)*b2-avx(3)*b3 print*,'b1=',b1,'b2=',b2,'b3=',b3,'b0=',b0 do i=1,n

yy(i)=b0+b1*x1(i)+b2*x2(i)+b3*x3(i) end do

open(3,file='e:\ad.txt')

!write(3,'(3f15.5)')((s(i,j),j=1,3),i=1,3) write(3,'(f15.5)')yy

write(3,'(f15.5)')y

close(3)

end

图13多元线性回归各系数的计算结果

所以得到回归方程为：

图14x1与y的效果图

图15多元线性回归拟合效果

实习六ＥＯＦ分析

给出海表温度距平数据资料sstpx.grd，以及相应的数据描述文件sstpx.ctl，对其进行EOF 分析，资料的时空范围可以根据sstpx.ctl获知。

数据在文件夹中单独给出，距平或者标准化距平处理后再进行EOF。

要求：实习报告中给出第一特征向量及其时间系数，并分析其时空特征。程序如下：

PROGRAM EOFPW

C THIS PROGRAM USES EOF FOR ANALYSING TIME SERIES

C OF METEOROLOGICAL FIELD

C M:LENTH OF TIME SERIES

C N:NUMBER OF GRID-POINTS

C KS=-1:SELF;KS=0:DEPATURE;KS=1:STANDERDLIZE

D DEPATURE

C KV:NUMBER OF EIGENVALUES WILL BE OUTPUT

C KVT:NUMBER OF EIGENVECTORS AN

D TIM

E SERIES WILL BE OUTPUT

C MNH=MIN(M,N)

C Evf=EIGENVACTORS(特征向量),tcF=TIME COEFFICIENTS FOR EGVT.

C ER(KV,1)=LAMDA,LAMDA EIGENVALUE

C ER(KV,2)=ACCUMULATE LAMDA

C ER(KV,3)=THE SUM OF COMPONENTS VECTORS PROJECTE

D ONTO

c EIGENVACTOR.

C ER(KV,4)=ACCUMULATE ER(KV,3)C

PARAMETER(M=516,N=216,MNH=216,KS=-1,KV=10,KVT=2)

PARAMETER(ff=-999.0,nx=18,ny=12)C

DIMENSION F(N,M),A(MNH,MNH),S(MNH,MNH),ER(mnh,4),

*DF(N),V(MNH),AVF(N),evf(N,KVT),tCF(M,KVT),

*data(Nx,ny),nf(N)

CCCCCCCCCCCCCCCCINPUT DATA

c OPEN(2,FILE='d:\eof\grd\eofz.DAT')

c read(2,*)F

c write(*,*)'ok1'

open(11,file='e:\eof\sstpx.grd',form='unformatted',

!access='direct',recl=nx*ny)

do132it=1,m

read(11,rec=it)((data(i,j),i=1,nx),j=1,ny)

do132jj=1,ny

do132ii=1,nx

kkkk=nx*(jj-1)+ii

f(kkkk,it)=data(ii,jj)

132continue

close(11)

CCCCCCCCCCCCCCCCINPUT DATA CCCCCCCCCCCCCCCCCCC ccccccccccccccccccccccccccccccccccccc

CALL Test1(n,m,ff,f,nf)

write(*,*)'ok2'

CALL TRANSF(N,M,F,nf,AVF,DF,KS)

write(*,*)'ok3'

CALL FORMA(N,M,MNH,F,A)

write(*,*)'ok4'

CALL JCB(MNH,A,S,0.00001)

write(*,*)'ok5'

CALL ARRANG(KV,MNH,A,ER,S)

write(*,*)'ok6'

CALL TCOEFF(KVT,KV,N,M,MNH,S,F,V,evf,tcf,ER) write(*,*)'ok7'

call test3(N,ff,nf,evf,kvt)

write(*,*)'ok8'

open(21,file='e:\eof\evf.grd',form='unformatted',

!access='direct',recl=nx*ny)

irec=0

do668kk=1,kvt

irec=irec+1

668write(21,rec=irec)((evf(nx*(j-1)+i,kk),i=1,nx),j=1,ny) close(21)

open(21,file='e:\eof\tcf.grd',form='unformatted',

!access='direct',recl=kvt)

irec=0

do345it=1,m

irec=irec+1

345write(21,rec=irec)(tcf(it,iik),iik=1,kvt)

close(21)

106format(10f8.4)

c OPEN(16,FILE='D:\pteof\grd\dats.DAT')

c CALL OUTER(KV,ER,mnh)

open(21,file='e:\eof\dats.dat')

write(21,106)(er(iiii,3),iiii=1,kv)

close(21)

c CALL OUTVT(KVT,N,M,MNH,S,F,Evf,tcF)

STOP

END

SUBROUTINE Test1(n1,m,ff,f,nf)

DIMENSION F(N1,M)

DIMENSION nF(N1)

do i=1,n1

nf(i)=0.0

enddo

do i=1,n1

do k=1,m

if(f(i,k).eq.ff)then

f(i,k)=0.0

nf(i)=nf(i)+1

endif

enddo

return

end

SUBROUTINE TRANSF(n1,m,f,nf,avf,df,ks)

C THIS SUBROUTINE PROVIDES INITIAL F BY KS and kv.

DIMENSION F(N1,M),AVF(N1)

DIMENSION DF(N1)

DIMENSION nF(N1)

if(ks.eq.-1)then

goto30

endif

do i=1,n1

avf(i)=0.0

enddo

if(ks.eq.0)then

goto5

endif

do i=1,n1

df(i)=0.0

enddo

cccccccccccccccccc

5continue

DO141I=1,N1

if(nf(i).ne.0)goto141

do12j=1,m

12AVF(I)=A VF(I)+F(I,J)

AVF(I)=A VF(I)/M

DO14J=1,M

14F(I,J)=F(I,J)-A VF(I)

141CONTINUE

IF(KS.EQ.0)THEN

RETURN

ELSE

DO241I=1,N1

if(nf(i).ne.0)goto241

DO22J=1,M

22DF(I)=DF(I)+F(I,J)*F(I,J)

DF(I)=SQRT(DF(I)/M)

DO24J=1,M

24F(I,J)=F(I,J)/DF(I)

241CONTINUE

ENDIF

30CONTINUE

RETURN

SUBROUTINE FORMA(N,M,MNH,F,A)

C THIS SUBROUTINE FORMS A BY F

DIMENSION F(N,M),A(MNH,MNH)

IF(M-N)40,50,50

40DO44I=1,MNH

DO44J=I,MNH

A(I,J)=0.0

DO42IS=1,N

42A(I,J)=A(I,J)+F(IS,I)*F(IS,J)

A(J,I)=A(I,J)

44CONTINUE

RETURN

50DO54I=1,MNH

DO54J=I,MNH

A(I,J)=0.0

DO52JS=1,M

52A(I,J)=A(I,J)+F(I,JS)*F(J,JS)

A(J,I)=A(I,J)

54CONTINUE

RETURN

END

SUBROUTINE JCB(N,A,S,EPS)

C THIS SUBROUTINE COMPUTS EIGENVALUES AN

D EIGENVECTORS OF A

DIMENSION A(N,N),S(N,N)

DO30I=1,N

DO30J=1,I

IF(I-J)20,10,20

10S(I,J)=1.

GO TO30

20S(I,J)=0.

S(J,I)=0.

30CONTINUE

G=0.

DO40I=2,N

I1=I-1

DO40J=1,I1

40G=G+2.*A(I,J)*A(I,J)

S1=SQRT(G)

S2=EPS/FLOAT(N)*S1

S3=S1

L=0

50S3=S3/FLOAT(N)

60DO130IQ=2,N