汉江安康站河流最大和最小径流量的数学模型.doc

1、汉江安康站最大和最小径流量的数学模型饶梦彬 杨婕妤 欧林海摘要： 河流的径流量是一个非平稳的时间序列。由于径流量的变化受主要受气候的影响，而影响气候的因素又是复杂多样的，对于这样一个复杂的系统，常规的思路如回归分析往往误差较大。基于时间序列的角度出发，我们分别建立了Holt双指数平滑模型和自回归平滑模型ARMA(p,q)。在ARMA(p,q)模型中，我们采用动态参量法在Matlab中得到拟合效果最佳的p和q的值分别为4和3。并通过估计参数得到了预测方程，从1991到1997年的汉江安康站最大、最小径流量实际值和预测值来看，预测的误差较小，同时用该方程得到1998至2002年的最大、最小径流量的

2、预测值。人工神经网络由于能够实现各种复杂的功能而备受各领域学者的关注，在模型三中，我们建立了一个以S型函数为传递函数的三层BP神经网络，并把原始数据归一化至0.2，0.8，使用Levenberg-Marquardt算法对1943到1990年的48个样本进行训练，达到了较高的精度，并且得到的1991至1997年的预测值与实际值比较接近。利用这个BP神经网络，我们得到了1998到2002年汉江安康站最大、最小径流量进行预测，并给出了模型的误差分析。关键词：时间序列、Holt双指数平滑、ARMA(p,q)、BP神经网络23一、问题重述：气候是重要的环境因素，研究我国干旱和半干旱地区的气候变化规律，对

3、确定陕西的经济发展战略，制定发展规划具有重要意义。1请根据陕南汉江安康站统计的最大、最小径流量数据（附录1），分析这些数据之间的关系；2建立最大、最小径流量适当的数学模型，并检验模型的合理性；3利用您所建立的模型，对1998，1999，2000，2001，2002年汉江安康站的最大、最小径流量进行预报，并与实际情况进行比较。二、问题分析：由题目中给出的汉江安康站1943年至1997年55年最大、最小径流量的数据，我们判断这是属于利用统计预测方法建立模型的问题。再者，由于数据的无规则性，以及与随机因素影响很大，可以其归为时间序列预测的一类。1对分析数据之间的关系的理解由于该问题中只有最大、最小径

4、流量的数据，我们判断可以通过分析数据的平均趋势、自相关系数、偏相关系数来观察最大、最小径流量的总体趋势以及数据的横、纵关系。2 对问题2的理解为了建立汉江安康站最大、最小径流量的数学模型，我们从时间序列的理论方法出发。根据相关理论建立较好的预测方程，在得到预测方程的基础上，分析残差，检验出模型的误差值和可靠度。从而确定模型的可行性。3 关于19982002年径流量的预测根据建立的数学模型，我们可以预测1998至2002汉江安康站的最大、最小径流量。同时为了预测效果，我们可以从已经给出的55组数据中，留下一小部分来验证模型，即预测值与实际值进行比较。由上述分析，本文在关于建立最大、最小径流量模型

5、的过程中采用三个模型，逐步精确。三、符号说明：、为平滑指数；为实际观测值；为外推预测时期数 4次移动平均数P（双指数模型）时序项数最大径流量时间序列最小径流量时间序列对预处理后所得时间序列对预处理后所得时间序列对零均值后的时间序列对零均值后的时间序列对四次差分处理后的时间序列对四次差分处理后的时间序列 时间序列的最小值和最大值y为映射后的数据。M为训练样本总数;n为输入层单元数;H为隐层单元数四、模型假设：（1）题目中给出的关于汉江安康最大、最小径流量的数据是准确无误的，能够反映实际情况；（2）不考虑1997年以后人类活动以及其他原系统以外因素对汉江安康及其上游径流量的影响。五、初始数据分析1

6、原始数据的时间序列图根据汉江安康站19431997年的最大、最小径流量数据，在EXCEL分别得到它们的时间序列图：图（1）图（2）从图（1）、图（2）中可以看出，汉江安康站记录的最大径流量和最小径流量在55年里的变化是不平稳的，波动幅度较大。在某些年份，它的数值会突然变得很大（如1983年的最大径流量和1984年的最小径流量）。主要的原因是河流流量受气候的影响，而对气候产生影响的因素是复杂多样的，从而导致径流量的随机波动较大。2 时间序列的长期趋势分析为了降低时间序列中随机扰动成分的影响，我们在Minitab软件中对最大径流进行3年移动平均处理，得到下图：图（3）在3年移动平均的分析报告中，我

7、们得到它的平均百分误差（MAPE）为57，进一步表明，径流量是一个非平稳的序列。因此，对于这种动态数据而言，用一般的回归模型进行按拟合往往误差较大。3 最大径流量与最小径流量的相关系数为了考察最大径流量与最小径流量的相关性，我们可以计算它们之间的相关系数。两个变量的相关系数计算公式为：经计算，我们得到最大径流量与最小径流量的相关系数为0.237。说明两者之间存在较弱的正相关。4 变量的自相关系数自相关分析法是进行时间序列分析的有效方法，它简单易行、较为直观，根据绘制的自相关分析图和偏自相关分析图，我们可以初步地识别平稳序列的模型类型和模型阶数。利用自相关分析法可以测定时间序列的随机性和平稳性，

8、以及时间序列的季节性。自相关函数的定义：滞后期为k的自协方差函数为：，则的自相关函数为：，其中。当序列平稳时，自相关函数可写为：。样本自相关函数为：，其中，它可以说明不同时期的数据之间的相关程度，其取值范围在-1到1之间，值越接近于1，说明时间序列的自相关程度越高。用Minitab求得的最大径流量与最小径流量的自相关系数如表1。表1 自相关系数表滞后最大径流量ACF最小径流量的ACF10.2090930.55508620.0899430.47210130.1795730.39355840.1437430.29321550.1854640.2662486-0.0691060.197657-0.0

9、22150.1295918-0.1010030.07712990.0411130.075182100.0460840.03446611-0.2644120.00107112-0.1384840.03282313-0.1580050.00347914-0.2354550.1925965 变量的偏自相关分析样本的偏自相关函数： 其中，。同样的，由Minitab软件可求偏相关系数如表2。表2 偏相关系数表滞后最大径流量PACF最小径流量PACF10.2090930.55508620.0483360.23700830.1588690.09502240.08029-0.01582650.1378190.

10、0471316-0.175141-0.0198957-0.025105-0.0501868-0.165399-0.044890.1126880.03635100.033517-0.02259911-0.220116-0.03580612-0.0695160.06021513-0.123085-0.01085314-0.2038590.2917596 互相关系数（CCF） 最大、最小径流量序列的互相关系数有下式给出： 在Minitab中，我们得到它的互相关函数图（图4）。图（4）六、模型的建立与求解模型一：双指数（霍尔特指数）模型对于为汉江安康站最大、最小径流量建立适合的数学模型，本文首先建立双

11、指数模型。双指数模型又称为霍尔特（HOLT）指数。霍尔特（HOLT）指数平滑法是一种线性指数平滑方法。这种方法最突出的优点是对具有趋势变动的时间数列，不用二次指数平滑，而是对趋势数据直接进行平滑并对原时间数列进行预测。这种方法因具有很大的灵活性而被广泛地使用着。1 移动平均处理为了建立合适的双指数模型，我们先对数据做移动平均处理。移动平均法，是指根据时间数列资料，逐项递推移动，依次计算包含一定项数的扩大时距。设序时项数为P，则移动平均数见（1-1）和（1-2）：当P为奇数是有： 当P为偶数时有： （1-1） （1-2）预测模型：趋势值预测：若P为奇数项：将每一移动平均值置于序时项数末期，然后再

12、外推一期（P为奇数）或两期（P为偶数）作为预测趋势值。本文利用Minitab求解的结果如下表1表3 原始数据四年移动平均后的结果最大径流量四年移动平均值最小径流量四年移动平均值年份数据年份数据年份数据年份数据1945796019718025194555.4875197174.5519469546.25197210728.75194660.4125197266.85194710183.75197313700194769.125197357.037519488685197414766.25194889.3125197456.119499715197514410194998.825197559.11

13、95011055197613062.5195094.2375197657.0125195111876.25197712225195192.075197755.1875195212821.25197813133.75195289.575197851.0375195312982.5197914740195386.325197948.65195413195198015400195483.8625198051.4875195513436.25198117150195582.65198155.0875195614763.75198219587.5195676.2125198273.12519571386

14、1.25198319637.5195777.575198391.8375195812322.5198418095195881.7625198493.1625195912333.75198515377.5195980.3875198588.7875196011195198611888.75196084.725198667.325196112058.7519878803.75196183.675198740.75196213047.519887680196286.8875198826.26125196313948.7519896252.5196393.2625198916.417519641454

15、3.7519905011.25196491.17519908.6762519651300019916223.75196589.5519915.465196613112.519926723.75196685.77519926.62875196712142.519936445196780.419937.06375196810932.519946265196881.27519947.3087519691088019955308.75196986.62519956.8219708805197083.33752 双指数模型的建立霍尔特指数平滑方法有二个基本平滑公式和一个预测公式。二个平滑公式分别对时间数

16、列的二种因素进行。它们是： （1-3） （1-4）预测公式是： （1-5）、为平滑常量；为实际观测值；为外推预测时期数。式（1-3）是对时间数趋势因素的平滑式；式（1-4）是对趋势增量的平滑式。3 、的估计为了估计、的值，本文利用步长加速法原理，并用计算机求解，原理如下：设目标函数：:即寻找使实际值与向前一步预测值（）离差平方和最小的平滑参数。为此，我们任选一基点，算出此点的目标函数值。然后对某个平滑参数按某一步长进行探索，即比较的目标函数值，以目标函数值最小的点为临时矢点；再由此点出发沿对另一参数进行同样的探索，如能得到比以前更好的点，就以该点代替前面的点作为新的临时矢点。如此对各个参数轮流

17、探索一遍，并选这一轮探索最好的点（最后的临时矢点）为第二个基点。由第一个基点到第二个基点构成了第一个模矢。对第一个基点来说，这是使目标函数得以改善的最有利的移动方向，沿这一方向前进，目标函数下降最快（就 而言）。显然，这一方向近似于目标函数的负梯度方向。下一步就把第一个模矢加速一倍，进行第二轮探索。余此类推。经过若干轮探索，即可以得到越来越好的目标函数下降点。如果探索进行到某一步时，得不出新的下降点，则应缩小步长以进行更精细的探索。当步长已缩小到某一精度要求，但仍得不到新的下降点时，即可将该点作为所求的近似最优点，就此停止迭代。用计算机求得的值如表4。表4 双指数平滑模型系数最大径流量双指数平

18、滑模型最小径流量双指数平滑模型水平系数1.98154水平系数0.8趋势系数0.0044趋势系数0.8以及最大径流量的拟合效果图（图5）。图（5）4 拟合误差分析利用Minitab软件，我们得到它们的拟合误差情况（表5）。表5 双指数平滑模型误差分析误差项最大径流量最小径流量平均百分误差（MAPE)819.689平均绝对误差8727.426平均偏差平方和1151520133.634从上表来看，Holt双指数平滑模型能够比较好地拟合汉江安康站所记录的多年的径流量数据。但由于径流量时间序列中的不规划成分占的比例较大和双指数平滑模型本身的特性，导致了它的预测效果不太好，因此，我们无法用该模型来对未来几

19、年的径流量作合理预测。模型二：自回归平滑模型（ARMA）1 数据预处理为了使数据满足ARMA（p,q）模型的应用条件，我们对数据做如下处理：（1）为了消除偶然因素的影响做三年平滑处理 = （2）为了压缩数据起伏度，对数据开方 （3）求零均值 （4）四阶差分处理 数据预处理的到的数据，见附表（1）。2 模型的建立与求解 （1）模型的建立 通过对数据的预处理，得时间序列，且E（W）=0，D（W）=和E（Z）=0，D（Z）=。因此它们符合达到平稳化要求，可建立ARMA（p,q）模型，得到ARMA(p,q)模型： （2） p,q的确定根据动态参数法，通过matlab和minitab软件得p=4,q=3

20、.（3）参数估计 模型的估计就是对参数的有效估计。本文采用Yule-Walker法（简称Y-W法），得到了精确度较高的参数，具体步骤如下：1）Y-W法估计参数序列的自相关系数与，有如下关系：则由可确定参数，。结果为（）2）又因为 （k=0,1kq,kq）由此式可得，。具体如表6。表6 最大径流量ARMA参数估计类型系数系数标准误TPAR11.56270.17538.910AR2-0.91950.3062-30.005AR30.34520.30821.120.269AR4-0.08020.1722-0.470.644移动平均10.57140.12144.710移动平均2-0.72730.081-

21、8.980移动平均31.08120.16966.370常量-0.10790.1928-0.560.579平均值-1.1742.098表7 最小径流量ARMA模型估计参数类型系数系数标准误差TPAR12.05250.22958.940AR2-1.57940.4029-3.920AR30.44330.38081.160.251AR40.02790.18710.150.882移动平均10.73860.21473.440.001移动平均2-0.73230.1288-5.690移动平均30.8870.18914.690常量-0.030050.01288-2.330.025平均值-0.53940.2311

22、3 模型的检验本文利用白噪声的随机性对ARMA(P,A)模型进行检验。对所建模型计算到其误差项所构成的时间序列的白自相关系数，对进行随机性检验，又因为修正后的Box-pierce卡方统计量如下：表8 最大径流量修正Box-Pierce(Ljung-Box)卡方统计量滞后12243648卡方2030.239.246自由度4162840P值00.0170.0780.238表9 最小径流量修正Box-Pierce(Ljung-Box)卡方统计量滞后12243648卡方6.512.824.531.8自由度4162840P值0.1630.6870.6560.819若数据误差项的白自相关系数落在95%的置

23、信区间之内，就说明误差项是随机的，由Q=n，计算最大、最小径流量的，均小于，说明模型是合理的，可以用来预测。4 模型的预测利用所建模型对1998-2002年的最大最小径流量进行预测：ARMA（p,q）序列的递推预测式由上分别得出ARMA模型的最大、最小径流量预测值以及还原后的最大、最小径流量的预测，如下：表10 ARMA模型的最大径流量预测值周期预测95%置信区间下限上限55-18.4777-39.34342.388156-3.5683-32.948425.8118578.3107-32.4949.11135811.0286-30.947453.0046599.7344-32.556752.0

24、255表11 还原后的最大径流量预测值周期年份预测值5519972741.45619984732.757199944395820004937.65920013958.16020023192.8表12 ARMA模型的最小径流量预测值周期预测95%置信区间下限上限55-1.4456-2.79005-0.1011456-2.19993-4.419920.0200657-1.64298-4.976691.6907258-0.5759-4.326273.17447590.3673-3.406394.14098600.84375-2.994664.68216表13 还原后的最小径流量预测值周期年份预测值5

25、519973.15619985.35719995.75820004.95920013.86020023.4模型三：BP-NN网络模型对于水资源系统而言，由于年径流量的随机复杂性与年径流量的各种因素的关联性，很难用一个准确的数学解析式来描述，我们知道以神经网络为代表的智能计算方法中，单个神经元有较简单却又反映非线性本质特征的非线性核，通过这些基本的非线性核自组织复合，使其能够重建任意的非线性连续函数。在时间序列预测中，前馈网络是最常使用的网络。在这种情形下，从数学角度看，网络成为输入输出的非线性函数。记一个时间序列为，对其进行预测可用下式描述:时间序列预测方法即是用神经网络来拟合函数，然后预测未

26、来值。数学上己经证明多层前馈网络具有很强的函数映射功能，1989年Robert-Nielson证明一个隐层的BP网可逼近闭区间内任意一个连续函数。理论上已经证明:具有偏差和至少一个S型隐层加上一个线性输出层的网络，能够逼近任何有理函数。大量BP网络模型在水文水资源中的应用表明，具有一个隐层就能满足要求，也就是说三层BP网络模型可刻化水文水资源研究对象，本文采用三层BP网络。图（6）为了优化BP网络，首先，我们必须选择合适的网络结构以减少训练次数而加快训练速度。其次，由于训练速度还与数据的规范化方法、训练速率、冲量项的引入以及联接权的初始值等因素有关，因而本文首先通过调整这些因素加快训练速度。本

27、文采用的BP网络以S形函数作为转换函数，该函数的值域为0，1区间，通常采用标准的归一化方法来实现。S函数的表达式如下（3-1），图形见图3-1。 （3-1） 1 数据归一化处理由于规范后每一个输出的教师值序列至少有一个值为0，一个值为1，恰是S形函数的极小值和极大值，要求联接权足够大才能使网络的输出值与其匹配，从而需要相当多的训练次数来不断修正权值，导致训练速度缓慢。为了避免这种现象，本节中将数据归一到02，0.8区间，归一化公式如式（3-2）： （3-2）2 隐含层单元数的确定用于预测的三层BP网络存在一个最佳隐层单元数问题。通常采用“试错法”(Trial-and- Ermr)确定隐层单元数

28、。首先给定较小初始隐单元数，构成一个结构较小的BP网络，进行训练。如果训练次数很多或则在规定的训练次数内没有满足收敛条件，停止训练，逐渐增加隐层单元数形成新的网络重新训练。另外，Weigend等人根据图象识别问题建议用下式（3-3）确定隐层单元数，即： （3-3）现在为提取较好的隐层单元数，本文综合（3-3）以及经验来确定。在神经网络学习时，必须选定迭代计算过程中的两个重要参数:学习速率和动量参数a。一般情况下，和a的取值目前无一致的意见，在变步长改进算法的同时，通常也是通过“试错法”反复训练，提取其中较为理想的值。随着科学技术的发展，一定量的训练次数己经不成问题，本节在初步试错的基础上，主要

29、采用了不断增加训练次数，以使拟合和预测结果能达到实际要求。权重和闭值通过网络学习最后确定。为对汉江安康站最大、最小径流量进行预测，本文引入BP网络模型进行预测。建模所采用的资料见问题重述部分。模型的数据本文采用从1943年到1994共48年的径流量资料。所建的模型为三层BP网络模型，4个输入，1个输出，隐含层用5个神经元来建立。3 模型的计算机实现首先，由公式（3-1）对最小年径流量做归一化处理，结果见表14。表14 汉江安康站最小径流量归一化处理后的数据年份归一化值年份归一化值年份归一化值年份归一化值19430.430719550.364719670.441819790.355319440.

30、47119560.4819680.500219800.36219450.339219570.473319690.443219810.36219460.226819580.394519700.508619820.403219470.516919590.553819710.511619830.401219480.516919600.403219720.342219840.819490.52719610.513219730.365319850.42619500.580719620.47119740.379819860.374719510.418319630.449219750.377119870.3

31、12419520.492519640.594119760.451519880.31319530.493519650.493519770.336519890.219540.547119660.434819780.352619900.212其次，进行样本训练，将以上样本导入网络，设置好参数，样本输入分为4组，构建一个451（输入层隐含层输出层）的三层BP-NN;对最小径流量训练次数达891次，精确度达0.0037.训练结果如表15。表15 最小径流量的训练结果年份预测值实际值绝对误差年份预测值实际值绝对误差19430.46390.4307-0.033219670.47390.4418-0.0321

32、19440.46390.4710.007119680.47840.50020.021819450.46390.3392-0.124719690.48470.4432-0.041619460.46390.2268-0.237119700.49360.50860.01519470.46390.51690.05319710.50550.51160.006119480.46390.51690.05319720.34220.3422019490.46390.5270.063119730.35580.36530.009619500.46390.58070.116819740.37930.37980.000

33、519510.46390.4183-0.045619750.40560.3771-0.028519520.46390.49250.028519760.43310.45150.018419530.46390.49350.029519770.33650.3365019540.4640.54710.083219780.34110.35260.011519550.4640.3647-0.099319790.36130.3553-0.00619560.4640.480.01619800.37660.362-0.014619570.46410.47330.009219810.38440.362-0.022

34、419580.46420.3945-0.069719820.37180.40320.031419590.46440.55380.089519830.40120.4012019600.46460.4032-0.061319840.790.80.0119610.46490.51320.048319850.46650.426-0.040519620.46540.4710.005619860.31370.37470.061119630.46610.4492-0.016919870.27370.31240.038719640.46710.59410.12719880.26560.3130.0475196

35、50.46860.49350.024819890.26430.2-0.064319660.47080.4348-0.03619900.26440.212-0.0524根据上表15用MATLAB绘制图形如下：图（7）最后，用上述模型做预测，把样本模拟的网络导入，得到表16。表16 BP神经网络对19911997年的最小径流量预测值与实际值对比年份1991199219931994199519961997实际值5.764.518.59.45.936.32.8预测值5.70 5.70 5.71 5.79 4.88 4.64 3.59 绝对误差-0.06 1.19 -2.79 -3.61 -1.05 -

36、1.66 0.79 对训练后的拟和、预测结果作相应的反归一化处理。处理后得出汉江安康站最小径流量时间序列的模拟预测结果如表16，从表16中可以看出，应用BP-ANN模型对汉江安康站最小径流量进行预测具有较高的精度，这是因为神经网络具有很强的非线性拟和能力。由BP网络模型我们得到汉江安康最小径流量从19982002的预测结果见表17。表17 BP神经网络对汉江安康19982002年最小径流量的预测值年份19981999200020012002预测值5.70 5.70 5.79 4.44 3.59 用同样的方法，我们用MATLAB得到最大径流量的训练结果见表18。表18 BP神经网络对最大径流量的

37、预测值与实际值对比（归一化数据）年份预测值实际值绝对误差年份预测值实际值绝对误差19430.41220.4307-0.112719670.40970.44180.005919440.41220.471-0.154919680.40530.50020.168519450.41220.3392-0.046419690.3940.4432-0.081719460.41220.22680.023519700.36790.5086-0.035819470.41220.5169-0.119119710.32220.51160.003319480.41220.51690.105619720.30560.34

38、220.025919490.41220.527-0.204719730.45880.3653-0.029119500.41220.5807-0.058219740.61460.37980.033219510.41220.41830.127619750.51220.3771-0.026619520.41220.49250.073519760.32830.45150.013919530.41220.4935-0.041219770.36860.3365-0.006719540.41220.5471-0.070419780.48890.35260.010819550.41220.36470.1656

39、19790.50160.3553-0.001919560.41220.480.069519800.49440.362-0.020719570.41220.47330.001519810.47290.3620.014819580.41220.39450.099619820.47970.4032019590.41220.5538-0.148919830.79430.40120.005719600.41220.40320.137619840.57250.8-0.008819610.41220.5132-0.064819850.40050.4260.005119620.41210.471-0.016519860.34470.3747-0.0319630.41210.44920.105619870.33150.31240.198219640.41210.59410.041619880.32860.313-0.053319650.41190.49350.175919890.3280.2-0.12819660.41130.4348-0.160819900.32780.2120.0125以及1998年2002年最大径流量的预测数据表19。表19 BP神经网络对汉江安康站1