[论文]河流最大径流量问题探究.doc
课程设计（论文）课程名称： 应用随机过程 设计题目： 河流最大径流量问题探究 院 系： 电子与信息技术研究院 班 级： 通信工程一班 设 计 者： 学 号： 指导教师： 设计时间： 2009-12-20 河流最大径流量问题研究数学模型预测方法主要有自回归模型(AR)、滑动平均模型(MA)和自回归滑动平均模型(ARMA)等，这些线性预测模型考虑因素较简单。自回归滑动平均模型(ARMA)计算简单，易于实时更新数据。河流的最大径流量是一种典型的时间序列，实际河流每年最大径流量的大小是一个依时间变化的过程，在这里我们取1年作为一个时间段来测量数据。下面是某条河流上的一个水文站从1915年到1973年记录的年最大径流量见表1的栏，共59个数据。将原始样本数据经过处理后变成时间序列，具体1156006931319900123128960291327310-135931040017313390403714106001931347310-13595108002131358850181698801211367840-829798501181371070020318109002231386190-2479988101413996109411099601291407580-1089111220035314199901321127510-1159426150-2519138640-29438250-419146380-2289446030-2639156810-1859458980311168820151466180-248917144005731479630961187440-1229489490821197240-1429492340-6329206430-22395011100243121110002331515090-3579227340-132952109002231239260591536490-2179245290-337954126003931259130461556640-2029267480-1189567430-1239276980-1689576760-190928965098158100001331297260-140959930063130875081表1 原始数据的计算过程如下所示：(1)求取误差时间序列,先计算令 则的计算数据如表1所示。(2)计算自协方差基函数的值根据，由于，则，表示样本数，表示计算的步骤。根据公式：上式中表示样本的数量，表示第个样本的自协方差函数值，计算如下：(3)计算样本的自相关函数值，根据公式 ，k=0，1，2上式中表示第个样本的自相关函数值，表示第个样本的自协方差基函数值，计算如下：(4)计算样本的偏相关函数值根据公式的计算如下表2所示，自相关函数值和偏相关函数值如表3所示k012345675020385-11569491470173-8171021421137-453251097564388677k891011121314-21049261836124109430039-370638157190-232551表2 基函数值kk1-0.23-0.2380.00-0.0720.290.2590.05-0.023-0.16-0.06100.02-0.0240.280.19110.09-0.025-0.010.1312-0.07-0.1160.220.15130.03-0.1070.080.1914-0.05-0.04表3 计算自相关函数值和偏相关函数值(5)判断和的截尾性和拖尾性。根据计算获得样本数据自相关函数值和偏相关函数值，绘制和的曲线趋势来判断。图1 自相关函数值的变化曲线趋势图图2偏相关函数值的变化曲线趋势图对AP()模型，已知当很大时，有，于是当>时，平均20个中至多有一个使，那么可认为截在=处。对MA()模型，同样有类似的性质，当>时，平均20个中至多有一个使，那么可认为截在=处。从图1中可看出拖尾，而从图2中看出截尾，又，>2时，所以可以认为尾巴截断在=2处。这里，取截尾在=2处较为合理。(6)确定模型的一般形式根据（5）得到的模型的阶数，则可以确定该模型为AR (2)。可得到AR(2)模型的一般形式如下：AR预测模型的阶数确定后，模型中的参数还需要进一步的求解其中各参数的值可通过以下3个公式确定：经过计算可得=-0.172，=0.250，=44553433。所以AR的预测模型如下：至此，我们已经得到了预测河流最大径流量的随机线性模型。参考文献：1. 田波平：应用随机过程讲义，20092. 刘波：Matlab信号处理，电子工业出版社，20063. 汪荣鑫：随机过程，西安交通大学出版社，2001附录：Matlab源程序close all;clear all;Zm=8669;n=59;K=14;Wt=load('data.txt');gama=zeros(1,K+1);for k=0:K for j=1:n-k gama(k+1)=Wt(j)*Wt(j+k)+gama(k+1); end gama(k+1)=gama(k+1)/n;endgama=round(gama);ro=zeros(1,K); for k=1:K ro(k)=gama(k+1)/gama(1);endro=round(ro*100)/100;fi=zeros(K,K);fi(1,1)=ro(1);fai=zeros(1,K);fai(1)=1;fai(2)=fi(1,1);for k=1:K-1 sro=0;srof=0; for j=1:k sro=ro(k+1-j)*fi(k,j)+sro; srof=ro(j)*fi(k,j)+srof; end fi(k+1,k+1)=(ro(k+1)-sro)/(1-srof); fai(k+2)=fi(k+1,k+1); for j=1:k fi(k+1,j)=fi(k,j)-fi(k+1,k+1)*fi(k,k-j+1); endendfai=round(fai*100)/100;ro=1,ro;x=zeros(1,K+1);plot(0:K,ro);xlabel('k');ylabel('自相关函数'); set(gca,'xtick',0:K);set(gca,'yminortick','on');hold on;plot(0:K,x,'k')set(gca,'ticklength',0.05 0.025);set(gca,'tickdir','out');grid onfigure;plot(0:K,fai);xlabel('k');ylabel('偏相关函数');set(gca,'xtick',0:K);set(gca,'yminortick','on');hold on;plot(0:K,x,'k')set(gca,'ticklength',0.05 0.025);set(gca,'tickdir','out');grid on
