成信工统计天气预报实验指导

统计天气预报实验指导实验一回归分析一、目的和要求：回归分析是统计天气预报中最基础和最重要的方法，是必须掌握的方法之一。该方法既能作为天气分析模型，又能作为天气预报模型，在实际工作中有极强的实用意义。通过该实验，掌握运用物理统计方法建立区域降水预测的基本步骤，运用提供的资料和方法子程序，编写或补充程序当中的部分片断，了解区域降水的预报方法及其建立过程，输出实验要求的相应结果，并就方法对区域降水的拟合及实验预报效果进行分析，从而深刻理解回归统计模型的意义，理解和掌握统计分析中的核心概念如平均值、方差、协方差、相关系数等的分析计算，掌握多元回归模型和逐步回归模型的计算步骤、原理和计算差异、显著性检验的意义，为统计结果的分析、统计模型的实际应用打下一定的基础。二、实验的主要内容：我国属于季风气候，降水量的季节差异很大，年降水的大部分集中在夏半年。影响我国夏季降水的物理因子很多，根据多年的研究，在业务预测中常考虑的有副热带高压和季风活动、中高纬阻塞与冷空气的活动、海温特别是与ENSO相联系的赤道太平洋海温等。这里选用河北省10个地区500hPa高度距平和（1月、2～3月、4月），石家庄地面WSW-WNW风速合计（12月下旬一次年1月下旬）作为预报因子，采用多元回归分析和逐步回归分析，建立区域降水的预报方程。具体如下：2.1利用河北省10个地区500hPa高度距平和（1月、2～3月、4月），石家庄地面WSW-WNW风速合计（12月下旬一次年1月下旬）作为预报因子；2.2运用多元回归方法，对1951-1968年河北省6～8月降水量建立预报方程；2.3对1969-1971年进行多元回归预测实验；2.4运用逐步回归方法，对1951-1968年河北省6～8月降水量建立预报方程；2.5对1969-1971年进行逐步回归预测实验；2.6完成实验报告。三、步骤：3.1熟悉资料方法3.1.1资料资料预报因子和降水量，依次按列排放，时间从1951年-1971年，其中因子序列包括1月、2～3月、4月500hPa高度距平和，地面WSW-WNW风速合计共4项；降水量为河北省逐年6～8月。3.1.2多元回归方法回归分析是用来寻找若干变量之间统计关系的一种方法，利用所找到的统计关系对某一变量作出未来时刻的估计，称为回归预报值。回归作为一种统计模型，包括线性和非线性回归。由于线性回归比较简单，理论上也比较严谨，在气象中不少气象要素之间可以近似地存在这种线性的关系，因此，这种线性回归应用得比较多。预报量常常与多个因子之间存在线性相关关系，因而大部分气象统计预报中的回归分析都是多元回归。设预报量与m个因子的关系是线性的，则多元线性回归方程可表示为：其中Y是预报量，是预报因子，e为残差误差，，，，…为估计的回归系数。3.1.3逐步回归方法在气象预报中，对预报量的预报常常需要从可能影响预报y的许多因素中挑选一批关系较好的作为预报因子，应用多元线性回归的方法建立回归方程来作预报。逐步回归分析是在已选定的一批因子中得到“最优”的回归方程的一种常用方法。以下是几种实际计算方案来选择“最优”回归方程。逐步剔除方案：从包含全部变量的回归方程中逐步剔除不显著的因子。但此方案计算量大。逐步引进方案：在一批待选的因子中，考察它们对预报量y的方差贡献，挑选所有因子中方差贡献最大者，经统计检验是显著的，进入回归方程。但此方案不一定保证最后的方程是“最优”的。双重检验的逐步回归方案：将因子一个个引入，引入因子的条件是该因子的方差贡献是显著的。同时，每引入一个新因子后，要对老因子逐个检验，将方差贡献变为不显著的因子剔除。这种方案是利用求解线性方程中求逆同时并行的方法。此方案综合前两种方案的优点克服其缺点，所以本实验选择此方案。3.2效果分析3.2.1回归拟合效果的参数分析(1)残差平方和(SSR)(2)标准差(3)复相关系数其中当R近似等于1，则相对误差将近似0，说明回归效果好。(4)偏相关系数其中当偏相关系数越大，回归中因子的重要性越大。(5)回归方差(6)多元回归的F检验(7)方程中已引入l个因子时,其中第k个因子的方差贡献(8)对第k个因子作F检验3.2.2预测实验表1采用多元回归分析河北省6～8月降水量预测实验结果年份196919701971观测值预测值表2采用逐步回归分析河北省6～8月降水量预测实验结果年份196919701971观测值预测值3.3编写程序根据本实验的具体情况，参考主程序读入资料和预报量文件，建立回归方程并输出相关的参数与回归系数；然后根据拟合的系数，编写程序计算1951-1968年拟合估计值和对1969-1971年进行预测实验的预测值。3.3.1多元回归程序ccccccccccccccccccccccccccccccc多元回归子程序cccccccccccccccccccccccccccccccccsubroutineDYHG(x,y,m,mm,n,a,q,s,r,v,u,b,dyy) dimensionx(m,n),y(n),a(mm),b(mm,mm),v(m) b(1,1)=n do20j=2,mm b(1,j)=0.0do10i=1,n10 b(1,j)=b(1,j)+x(j-1,i)b(j,1)=b(1,j)20continuedo50i=2,mm do40j=i,mm b(i,j)=0.0 do30k=1,n30b(i,j)=b(i,j)+x(i-1,k)*x(j-1,k)b(j,i)=b(i,j)40continue50continuea(1)=0.0do60i=1,n60a(1)=a(1)+y(i)do80i=2,mm a(i)=0.0 do70j=1,n70a(I)=a(i)+x(i-1,j)*y(j)80continuecallCHOLESKY(b,mm,1,a,l) yy=0.0 do90i=1,n90yy=yy+y(i)/nq=0.0 dyy=0.0 u=0.0ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccdo110i=1,n p=a(1) do100j=1,m100p=p+a(j+1)*x(j,i)q=q+(y(i)-p)*(y(i)-p) dyy=dyy+(y(i)-yy)*(y(i)-yy) u=u+(yy-p)*(yy-p)110continuecccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc s=sqrt(q/n) r=sqrt(1.0-q/dyy) do150j=1,m p=0.0 do140i=1,n pp=a(1) do130k=1,m if(k.ne.j)pp=pp+a(k+1)*x(k,i)130continuep=p+(y(i)-pp)*(y(i)-pp)140continuev(j)=sqrt(1.0-q/p)150continuereturn endC!PerformtheCHOLESKYDecompositionsubroutineCHOLESKY(a,n,m,d,l) dimensiona(n,n),d(n,m) l=1 if(a(1,1)+1.0.eq.1.0)then l=0 write(*,30) return endif a(1,1)=sqrt(a(1,1)) do10j=2,n10 a(1,j)=a(1,j)/a(1,1)do100i=2,n do20j=2,i20a(i,i)=a(i,i)-a(j-1,i)*a(j-1,i)if(a(i,i)+1.0.eq.1.0)then l=0 write(*,30) return endif30format(1x,'fail')a(i,i)=sqrt(a(i,i)) if(i.ne.n)then do50j=i+1,n do40k=2,i40a(i,j)=a(i,j)-a(k-1,i)*a(k-1,j)50a(i,j)=a(i,j)/a(i,i)endif100continuedo130j=1,m d(1,j)=d(1,j)/a(1,1) do120i=2,n do110k=2,i110d(i,j)=d(i,j)-a(k-1,i)*d(k-1,j)d(i,j)=d(i,j)/a(i,i)120continue130continuedo160j=1,m d(n,j)=d(n,j)/a(n,n)do150k=n,2,-1 do140i=k,n140d(k-1,j)=d(k-1,j)-a(k-1,i)*d(i,j)d(k-1,j)=d(k-1,j)/a(k-1,k-1)150continue160continuereturn end 参考主程序(范例)C样本长度：N；因子个数：K；PROGRAMMAININTEGER,PARAMETER::N=18 INTEGER,PARAMETER::K=4REAL,DIMENSION(K,N)::X REAL,DIMENSION(N)::Y REAL,DIMENSION(K+1)::A REAL,DIMENSION(K+1,K+1)::B REAL,DIMENSION(K)::V REALQ,S,R,UCOPENTHEINPUTDATAFILEOPEN(10,FILE='x-1.TXT')OPEN(20,FILE='y.TXT')CREADTHEDATADOI=1,N READ(10,*)X(1,I),X(2,I),X(3,I),X(4,I)READ(20,*)Y(I) ENDDO CLOSE(10) MM=K+1C调用回归子程序callDYHG(X,Y,K,MM,N,A,Q,S,R,V,U,B,DYY) write(*,88)A(1)88format(/1x,'b0='f9.5)do89j=2,MM89write(*,100)j-1,A(j)100format(1x,'b',i2,'=',f9.5)ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccwrite(*,20)Q,S,R20format(1x,'Q=',f13.6,3x,'S=',f13.6,3x,'R=',f13.6)write(*,22)U,DYY22format(1x,'U=',f13.6,3x,'DYY=',f13.6)write(*,30)(i,V(i),i=1,K)30format(1x,'V(',i2,')=',f13.6)write(*,40)U40format(1x,'U=',f13.6)C输出回归系数及相关参数open(6,file='table.txt') write(6,180)180format(/2x,'regressioncoefficients:')write(6,88)A(1) do189j=2,MM189write(6,100)j-1,A(j)write(6,200)200format(1x,'GenericAnalysisofVarianceTablefortheMultiple&LinearRegression')write(6,202)202format(1x,'-----------------------------------------------------------&-------------------') write(6,204)204format(3x,'SourcedfSSMS')write(6,202) write(6,206)N-1,DYY206format(1x,'Totaln-1=',i2,'SST=',f13.4)u2=U/real(K) write(6,280)K,U,U2280format(1x,'RegressionK=',i2,'SSE=',f13.4,'MSR=SS&E/k=',f13.4) q2=q/real(n-k-1) write(6,209)n-k-1,q,q2209format(1x,'Residualn-k-1=',i2,'SSR=',f13.4,'MSE=SS&R/(n-k-1)=',f13.4) f=(U/real(K))/(Q/real(N-K-1))write(6,220)f220format(1x,'F=MSR/MSE=',f13.4)write(6,202) close(6) stop end 3.3.2逐步回归程序（略）3.4结果分析(范例)采用多元回归分析方法输出数据结果及图表，并用结果说明方法拟合预测的效果。以下为4个因子，18年样本的多元回归分析示例，仅供参考。图1输出的回归系数和相关参数图2预报量的观测值和预测值（应用EXCEL或Grads绘图，柱形图为预报量的预测值，折线图为预报量的观测值）逐步回归分析结果（略）3.5实验报告内容3.5.1实验目的3.5.2实验内容（含实验资料数据介绍）3.5.3计算方法和使用软件3.5.4计算程序清单3.5.5输出反映回归效果的参数及回归系数，并就相关参数分析回归效果。3.5.6输出预报实验的预测值3.5.7预报量与回归方程计算的预测值和观测值的历年曲线变化图(图2)，并附简单说明。

实验二经验正交函数分解一、目的和要求：经验正交函数分解（EOF）是统计天气分析中气象要素场最基础的研究模型，是必须理解和掌握的方法之一，是后续课程中许多气象要素场的计算结果的理解的基础理论，也是毕业设计和论文中的基本分析方法。该方法用个数较少的几个空间分布模态来描述环流形势，而且基本涵盖环流场的信息，既能作为天气分析模型，其方法的延拓又能作为天气预报模型，在实际工作中也有极强的实用意义。通过该实验，深刻理解气象要素场的统计模型的意义，掌握气象要素场分析的基本方法，为实际预报业务和科研工作打下一定的基础。二、实验的主要内容：对（-，-）850hPa高度场进行经验正交展开（EOF.FOR），输出分析主要参数指标；绘制环流型图和相应的时间系数序列图，并加以分析。三、步骤：3.1熟悉资料方法3.1.1资料提供的资料为NCEP/NCAR60年（1948年-2007年）逐年1～12月的850hPa高度场资料，资料范围为（-，-），网格距为2.5*2.5，纬向格点数为144，经向格点数为73。资料为NC格式，资料从南到北、自西向东排列，每月为一个记录，按年逐月排放，注意读取方式以及记录长度。本次实验应用NCEP/NCAR（-，-）58年（1948年-2005年）逐年7月的850hPa高度场资料，纬向格点数为73，经向格点数为37。3.1.2方法（经验正交函数分解EOF）EOF（经验正交函数分解）是针对气象要素场进行的，其基本原理是把包含p个空间点（变量）的场随时间变化进行分解。设抽取样本容量为n的资料.则场中任一空间点i和任一时间点j的距平观测值可看成由p个空间函数和时间函数(k=1,2,…,p)的线性组合,表示成EOF功能是从一个气象场多次观测资料中识别出主要空间型及其时间演变规律。EOF展开就是将气象变量场分解为空间函数（V）和时间函数（T）两部分的乘积之和：X=VT。应用步骤：资料预处理(距平或标准化处理)计算协方差矩阵用Jacobi方法或迭代法计算协方差矩阵的特征值与特征向量将特征值从大到小排列计算特征向量的时间系数计算每个特征向量的方差贡献结果输出3.2编写程序要求编写主程序，其中包括资料读入，范围截取，子程序调用。注意：EOF的资料输入，时间场一维，空间场一维。*********************（附程序，对关键部分标志出）**********************EOF程序C**********************************************************************C*CPROGRAMNOTES*C*CTHISPROGRAMUSESEOFTOANALYSISTIMESERIES*COFMETEOROLOGICALFIELD*C*C**********************************************************************C*C********ParameterTable**********C*CMt===>LENTHOFTIMESERIES*CN===>NUMBEROFGRID-POINTS(orSTATIONS)*CKS=-1,SELF;KS=0,DEPATURE;KS=1,STANDERDLIZEDDEPATURE*CKV=NUMBEROFEIGENVALUESWILLBEOUTPUT*CKVT=NUMBEROFEIGENVECTORSANDTIMESERIESWILLBEOUTPUT*CMNH=Minimum(Mt,N)*CEGVT===>EIGENVECTORS,ECOF===>TIMECOEFFICIENTSFOREGVT*CER(KV,1)====>LAMDA;LAMDA===>EIGENVALUE*CER(KV,2)====>ACCUMULATELAMDA*CER(KV,3)====>THESUMOFCOMPONENTSVECTORSPROJECTEDONTO*CEIGENVACTOR.*CER(KV,4)====>ACCUMULATEER(KV,3)*C*C**********************************************************************PARAMETER(N=73*37,MT=58,MNH=58)PARAMETER(KS=1,KV=10,KVT=10)REALF(N,MT),AVF(N),DF(N),ER(MNH,4)REALA(MNH,MNH),S(MNH,MNH),V(MNH)c**************************************************************************cINFN－输入数据文件名；OUTERA－输出特征值及方差贡献、累积方差贡献的文件名（文本）;cOUTTC1－输出时间系数文件（文本）；OUTTC2－输出时间系数文件（二进制）；cOUTTEVT－输出特征向量文件（二进制）；c************************************************************************** CHARACTER*50INFN,OUTERA,OUTTC1,OUTTC2,OUTEVT DATAINFN/'hgt8501948-2005july.grd'/ DATAOUTERA/'hgt_XT03ER3.DAT'/ DATAOUTTC1/'hgt_XT03TC13.DAT'/ DATAOUTTC2/'hgt_XT03TC23.DAT'/ DATAOUTEVT/'hgt_XT03VT3.DAT'/ C----------------ReadORIGINALDATA----------------------------write(*,*)'Nowisreadingprimativefield!'OPEN(8,FILE=INFN,FORM='UNFORMATTED',ACCESS='DIRECT',RECL=N)DOIT=1,MTREAD(8,REC=IT)(F(IS,IT),IS=1,N)ENDDOpauseC**************STARTTORUNEOFPROGRAM******************************WRITE(*,*)write(*,*)'FIRSTSTEP'write(*,*)'formingtheinitialmatrix(F)byusingTRANSF!'CALLTRANSF(N,Mt,F,AVF,DF,KS)WRITE(*,*)write(*,*)'STEP2'write(*,*)'achivingthecovariancematrixbyusingtheFORMA!'CALLFORMA(N,Mt,MNH,F,A)WRITE(*,*)write(*,*)'STEP3'write(*,*)'caculatingtheeigenvalueandeigenvectors'WRITE(*,*)'byusingJacobmethod!'CALLJCB(MNH,A,S,0.001)WRITE(*,*)write(*,*)'STEP4'write(*,*)'arrangetheeigenvalueandeigenvectors' WRITE(*,*)'byusingARRANG!'CALLARRANG(MNH,A,ER,S)WRITE(*,*)write(*,*)'STEP5'write(*,*)'thecaculationofstandardeigenvectors'WRITE(*,*)'byusingTCOEFF!' CALLTCOEFF(KVT,N,Mt,MNH,S,F,V,ER)write(*,*)write(*,*)'STEP6'write(*,*)'outputingeigenvalueandaccumulationusingOUTER!'CALLOUTER(MNH,ER,OUTERA)WRITE(*,*)WRITE(*,*)'STEP7'write(*,*)'outputingthetimecoefficientoftheeigenvecters!'CALLOUTVT1(KVT,N,Mt,MNH,S,F,OUTTC1,OUTTC2)WRITE(*,*)WRITE(*,*)'STEP8'write(*,*)'outputingtheeigenvecters!' CALLOUTVT3(KVT,N,Mt,MNH,S,F,OUTEVT)ENDC***************FINISHTHEMAINPROGRAM*****************************C *CSUBROUTINEFUNCTION*C*CTHISSUBROUTINEPRINTSARRAYER*CER(KV,1)FORSEQUENCEOFEIGENVALUEFROMBIGTOSMALL*CER(KV,2)FOREIGENVALUEFROMBIGTOSMALL*CER(KV,3)FORSMALLLO=(LAMDA/TOTALVARIANCE)*CER(KV,4)FORBIGLO=SUMOFSMALLLO)*C*********************************************************************C--------SAVINGTHEEIGENVALUEANDERROR---------------------------*SUBROUTINEOUTER(MNH,ER,OUTERA)DIMENSIONER(MNH,4) CHARACTER*50OUTERAopen(30,file=OUTERA)WRITE(30,510)WRITE(30,520)WRITE(30,530)(IS,(ER(IS,J),J=1,4),IS=1,MNH)CLOSE(30)510FORMAT(25X,'EIGENVALUEANDANALYSISERROR')520FORMAT(5X,'N',8X,'LAMDA',10X,'SLAMDA',11X,'PH',12X,'SPH')530FORMAT(I6,2E15.6,2F15.5) RETURNENDC**********************************************************************CSUBROUTINEFUNCTION*C*CTHISSUBROUTINEPRINTSSTANDARDEIGENVECTORS*CANDITSTIME-COEFFICENTSERIES*C**********************************************************************C-------------savetime-coeffivcentseriedofS.E.------------SUBROUTINEOUTVT1(KVT,N,M,MNH,S,F,OUTTC1,OUTTC2)DIMENSIONF(N,M),S(MNH,MNH) CHARACTER*50OUTTC1,OUTTC2OPEN(31,file=OUTTC1)OPEN(32,file=OUTTC2,FORM='UNFORMATTED',ACCESS='DIRECT',RECL=KVT)WRITE(31,400)WRITE(31,200)(IS,IS=1,KVT)DOJ=1,MIF(M.GE.N)THENWRITE(31,300)J,(F(IS,J),IS=1,KVT)WRITE(32,REC=J)(F(IS,J),IS=1,KVT)ELSEWRITE(31,300)J,(S(J,IS),IS=1,KVT)WRITE(32,REC=J)(s(J,IS),IS=1,KVT)ENDIF ENDDOCLOSE(31)200FORMAT(3X,10I15)300FORMAT(I5,10E15.7)400FORMAT(30X,'TIME-COEFFICENTSERIESOFS.E.')RETURNENDC---------savestandardeignvectors------------------SUBROUTINEOUTVT3(KVT,N,M,MNH,S,F,OUTEVT)DIMENSIONF(N,M),S(MNH,MNH) CHARACTER*50OUTEVTOPEN(33,file=OUTEVT,FORM='UNFORMATTED',ACCESS='DIRECT',RECL=N)DOJS=1,KVTIF(M.GE.N)THENWRITE(33,REC=JS)(S(I,JS),I=1,N)ELSEWRITE(33,REC=JS)(F(I,JS),I=1,N)ENDIF ENDDOCLOSE(33) RETURN ENDC***********************************************************************CSUBROUTINEFUNCTION*C*CTHISSUBROUTINEPROVIDESINITIALFBYKS(optionalparameter)*Cks=-1,0,or1accordingtoprimativefield*C***********************************************************************SUBROUTINETRANSF(N,M,F,AVF,DF,KS)REALF(N,M),AVF(N),DF(N)IF(KS)30,10,1010DOI=1,NAVF(I)=0.0DF(I)=0.0ENDDODOI=1,NDOJ=1,MAVF(I)=AVF(I)+F(I,J)ENDDOAVF(I)=AVF(I)/MDOJ=1,MF(I,J)=F(I,J)-AVF(I)ENDDOENDDOIF(KS.EQ.1)THENDOI=1,NDOJ=1,MDF(I)=DF(I)+F(I,J)*F(I,J)ENDDODF(I)=SQRT(DF(I)/M)DOJ=1,MF(I,J)=F(I,J)/DF(I) ENDDOENDDOENDIF30CONTINUE RETURNENDC-----------------FORMA-----------------------------SUBROUTINEFORMA(N,M,MNH,F,A)REALF(N,M),A(MNH,MNH)IF(M-N)40,50,5040DOI=1,MNHDOJ=1,IA(I,J)=0.0DOIS=1,NA(I,J)=A(I,J)+F(IS,I)*F(IS,J)ENDDOA(J,I)=A(I,J) ENDDO ENDDORETURN50DOI=1,MNHDOJ=1,IA(I,J)=0.0DOJS=1,MA(I,J)=A(I,J)+F(I,JS)*F(J,JS)ENDDOA(J,I)=A(I,J) ENDDO ENDDO RETURNENDc***********************************************************************CSUBROUTINEFUNCTION*C*CTHISSUBROUTINECOMPUTSEIGENVALUESANDEIGENVECTORSOFA*c***********************************************************************SUBROUTINEJCB(N,A,S,EPS)DIMENSIONA(N,N),S(N,N)DOI=1,NDOJ=1,NIF(I.EQ.J)THENS(I,J)=1.0 ELSE S(I,J)=0.0 ENDIFENDDOENDDOG=0.0DOI=2,NI1=I-1DOJ=1,I1G=G+2.*A(I,J)*A(I,J) ENDDO ENDDOS1=SQRT(G)S2=EPS/FLOAT(N)*S1S3=S1L=050S3=S3/FLOAT(N)60DO130IQ=2,NIQ1=IQ-1DO130IP=1,IQ1IF(ABS(A(IP,IQ)).LT.S3)GOTO130L=1V1=A(IP,IP)V2=A(IP,IQ)V3=A(IQ,IQ)U=0.5*(V1-V3)IF(U.EQ.0.0)G=1.IF(ABS(U).GE.1E-10)G=-SIGN(1.,U)*V2/SQRT(V2*V2+U*U)ST=G/SQRT(2.*(1.+SQRT(1.-G*G)))CT=SQRT(1.-ST*ST)DOI=1,NG=A(I,IP)*CT-A(I,IQ)*STA(I,IQ)=A(I,IP)*ST+A(I,IQ)*CTA(I,IP)=GG=S(I,IP)*CT-S(I,IQ)*STS(I,IQ)=S(I,IP)*ST+S(I,IQ)*CTS(I,IP)=GENDDODOI=1,NA(IP,I)=A(I,IP)A(IQ,I)=A(I,IQ) ENDDOG=2.*V2*ST*CTA(IP,IP)=V1*CT*CT+V3*ST*ST-GA(IQ,IQ)=V1*ST*ST+V3*CT*CT+GA(IP,IQ)=(V1-V3)*ST*CT+V2*(CT*CT-ST*ST)A(IQ,IP)=A(IP,IQ)130CONTINUEIF(L-1)150,140,150140L=0GOTO60150IF(S3.GT.S2)GOTO50 RETURNENDc**********************************************************************CSUBROUTINEFUNCTION*C*CTHISSUBROUTINEPROVIDESASERIESOFEIGENVALUES*CFROMMAXTOMIN*C**********************************************************************SUBROUTINEARRANG(MNH,A,ER,S)DIMENSIONA(MNH,MNH),ER(MNH,4),S(MNH,MNH)TR=0.0DOI=1,MNHTR=TR+A(I,I) ER(I,1)=A(I,I) ENDDOMNH1=MNH-1DOK1=MNH1,1,-1DOK2=K1,MNH1IF(ER(K2,1).LT.ER(K2+1,1))THENC=ER(K2+1,1)ER(K2+1,1)=ER(K2,1)ER(K2,1)=CDOI=1,MNHC=S(I,K2+1)S(I,K2+1)=S(I,K2)S(I,K2)=C ENDDOENDIF ENDDO ENDDOER(1,2)=ER(1,1)DOI=2,MNHER(I,2)=ER(I-1,2)+ER(I,1) ENDDODOI=1,MNHER(I,3)=ER(I,1)/TRER(I,4)=ER(I,2)/TR ENDDO RETURNENDC**********************************************************************CTHISSUBROUTINEPROVIDESSTANDARDEIGENVECTORS*C(M.GE.N,SAVEDINS;M.LT.N,SAVEDINF)ANDITSTIMECOEFFICENTS*CSERIES(M.GE.N,SAVEDINF;M.LT.N,SAVEDINS)*C**********************************************************************SUBROUTINETCOEFF(KVT,N,M,MNH,S,F,V,ER)DIMENSIONS(MNH,MNH),F(N,M),V(MNH),ER(MNH,4)DOJ=1,MNHC=0.0DOI=1,MNH C=C+S(I,J)*S(I,J) ENDDOC=SQRT(C)DOI=1,MNH S(I,J)=S(I,J)/C ENDDO ENDDOIF(M.GE.N)THENDOJ=1,MDOI=1,NV(I)=F(I,J)F(I,J)=0.0 ENDDODOIS=1,KVTDOI=1,N F(IS,J)=F(IS,J)+V(I)*S(I,IS) ENDDO ENDDO ENDDOELSEDOI=1,NDOJ=1,MV(J)=F(I,J)F(I,J)=0.0 ENDDODOJS=1,KVTDOJ=1,MF(I,JS)=F(I,JS)+V(J)*S(J,JS) ENDDO ENDDO ENDDODOJS=1,KVTDOJ=1,MS(J,JS)=S(J,JS)*SQRT(ER(JS,1)) ENDDODOI=1,NF(I,JS)=F(I,JS)/SQRT(ER(JS,1)) ENDDO ENDDOENDIF RETURNEND3.3结果输出(范例)范例所用资料为NCEP/NCAR（-，-）58年（1948年-2005年）逐年7月的850hPa高度场资料，方差贡献及典型场分布（图1～图6）仅供参考。EIGENVALUEANDANALYSISERRORNLAMDASLAMDAPHSPH10.490105E+050.490105E+050.312850.3128520.223947E+050.714052E+050.142950.4558030.133197E+050.847249E+050.085020.5408340.109573E+050.956822E+050.069940.6107750.931490E+040.104997E+060.059460.6702360.665244E+040.111650E+060.042460.7127070.587009E+040.117520E+060.037470.7501780.463021E+040.122150E+060.029560.77972图1EOF第一特征向量场图2第一特征向量场对应的时间系数序列图3EOF第二特征向量场图4第二特征向量场对应的时间系数序列图5EOF第三特征向量场图6第三特征向量场对应的时间系数序列第一模态分析第一模态表现出大体以东经130o为界,西低东高的形势,即在东亚大陆气压较低而北太平洋地区气压较高,另从时间系数序列中可以看出这个形势总体趋势在减弱。3.4实验报告内容3.4.1实验目的3.4.2实验内容（含实验资料数据介绍）3.4.3计算方法和使用软件3.4.4计算程序清单3.4.5绘制前三个特征向量场3.4.6用所学知识分析图表资料

实验三功率谱分析和滤波一、目的和要求：周期是大气运动的基本理念之一，周期分析也是统计天气分析中的基础理论。功率谱分析的模型和滤波思想是最基本的周期分析统计模型，是必须理解和掌握的基本理论方法，是后续课程中许多天气分析事实理解的基础，也是毕业设计和论文中的基本分析方法。该方法作为天气分析模型，在实际业务工作和科研工作中有极强的实用意义。通过该实验，能使学生深刻理解气象时间序列周期分析的意义和周期分析的方法，为实际的业务和科研工作打下一定的基础。二、实验的主要内容：对原序列滤波前（或后）利用功率谱检测序列的周期，并进行显著性检验。本实验利用功率谱方法检测原序列的周期，并进行显著性检验；利用带通滤波器提取所需要的周期分量，再利用功率谱检验所提取周期的显著性。三、步骤：3.1熟悉资料方法3.1.1资料本次实验使用国家气候中心提供58年（1951年-2008年）夏季（6-8月）逐月降水资料。取得是长江中下游及江淮地区降水场第一主要因子的时间序列为例。3.1.2功率谱分析方法谱分析是时间序列在频域上进行分析的方法，亦称频谱分析或波谱分析。主要侧重在时间序列波动分析上。功率谱的概念是针对功率有限信号的，所表现的是单位频带内信号功率随频率的变换情况。功率谱密度是指对于具有连续频谱和有限平均功率的信号或噪声，表示其频谱分量的单位带宽功率的频率函数。气象上常利用功率谱作周期分析，功率谱图显示不同频率振动的功率大小，同时也是方差贡献的大小，因而可以从谱曲线中的谱值最大来确认主要振动及其对应的周期。功率谱估计是指对一实测时间序列现实（t=1，2，…，n)的估计，可有两种估计功率谱方法：离散功率谱估计和连续功率谱估计，本实验选择连续功率谱估计。具体求连续功率谱估计方法步骤如下：(1)计算样本落后自相关系数，m为最大步长，或最大落后时间长度，本实验m取18，通常m宜选择在n/3～n/10之间。（为时间落后步长）(2)求粗谱估计波数与圆频率的关系在取样点数为m时有：波数k最大可取到m/2。习惯上令l=2k，波数l取值点可从0到m。(3)计算平滑功率谱也可（式中l=0，1，…，m）(4)作谱图，以波数l为横轴，平滑功率谱密度估计值为纵坐标作图。波数l与周期关系为连续功率谱检验包括红色噪音过程和白色噪音过程，本实验选择白色噪音过程。3.1.3滤波方法对气象要素作谱分析过程中，一些规则周期占有很大的分量，但这种周期是众所周知的，它的存在，压低了其它周期的表现，一旦把它去掉之后，则可突出地表现其它周期的成分，即把感兴趣的周期成分从原来的序列中识别和提取出来，这一过程就是滤波过程。实际上是原始序列经一定的变换转化为另一序列的过程。滤波的总体过程：输入输出x(t)滤波器y(t)（数字处理）其中x(t)为原始序列，y(t)为新的时间序列。要求：过滤系统具有时间不变性和稳定性。滤波器可分为三类：用于从原序列中滤出低于某个频率成分的滤波器，称为低通滤波器；用于滤出高于某个频率成分的滤波器称为高通滤波器；滤出两个频率之间一个频带成分的滤波器称为带通滤波器。根据本实验的具体情况，选择带通滤波器，滤波范围为2～8年。3.2编写程序要求编写主程序，其中包括资料读入，范围截取，子程序调用。Cccccccccccccccccccccccccccc（附程序，对关键部分标志出）ccccccccccccccccccccccccc连续功率谱程序 parameter(n=50,m=n/3-1,ta=1.645,xx=12.59,pi=3.141592653589793) realx(n),r(0:m),ss(0:m),sr(0:m),sss(0:m),ssss(0:m),vv,ssr,sum,rc&tta,ave open(1,file='soi-1fil.txt') open(3,file='soi-1filsl1.txt')open(5,file='soi-1filtest1.txt') doi=1,n read(1,*)x(i) enddo close(1) **********************求方差***************** ex=0.0 doi=1,n ex=ex+x(i) enddoex=ex/n s=0.0 doi=1,n s=s+(x(i)-ex)**2 enddo s=sqrt(s/n)**********************求自相关系数***************** doj=0,m r(j)=0.0 doi=1,n-j r(j)=r(j)+(x(i)-ex)*(x(i+j)-ex)/(s*s) enddo r(j)=r(j)/(n-j)enddo*********************求sl***************************doi=0,m ssr=0.0dok=1,m-1 ssr=ssr+r(k)*(1+cos((pi*k)/m)*cos(i*pi*k/m))enddosr(i)=ssrif((i.ne.0).or.(i.ne.m))thenss(i)=(1.0/m)*(r(0)+sr(i))elsess(i)=1.0/2*m*(r(0)+sr(i)) endif enddodok=0,m write(3,'(f8.4)')ss(k) enddo close(2) **********检验******************* sum=0.0 doi=0,m sum=sum+ss(i) enddo ave=sum/mtta=n-2rc=(-1+ta*sqrt(tta))/(n-1)vv=(2.0*n-m/2.0)/m wri