随机水文学课件

随机水文学

第一章绪论本章主要内容第一节我国水资源情况第二节问题的提出第三节随机水文模型第四节随机水文的发展第五节随机水文模型的应用第一节我国水资源情况据2004年3月22日第十二届“世界水日”统计我国淡水的总量为28000亿立方米，居世界第四位（前三位为巴西、俄罗斯、加拿大），但人口平均占有量仅为2632m3，还不到世界平均的四分之一。居世界第108位左右。不仅人均水量低，而且水量在地区和时间上分布极不均匀，夏秋水量有余（汛期），冬春水量不足，北方水源亏缺，南方水量稍好。由于时空分布不均在我国经常出现旱涝灾害。

2003年淮河大水完工仅一个月的入海水道在了排洪中取到了巨大的作用1991年6、7月江苏、安徽、浙江、河南大水，为百年未遇，全国损失700亿元，死亡2200人，单江苏就损失100亿。

1998年长江、东北嫩江大水持续2个多月，死4000余人，损失2500亿元。

91年700亿元，94年1800亿元，95年1700亿元，96年2800亿元，1986年7、8月吉林连降十二场暴雨，松花江、东辽河、洮儿河、饮马河、辉发河、浑江三歧水40%农田受灾，损失45亿。

1985年6、8、9号台风袭击辽宁，10小时雨量达620mm，沈阳、鞍山、丹东、大连、盘锦、营口等50余县受灾，为了保住辽河大堤，炸沙山子村套堤保住了辽河大堤。

1975年8月河南林庄发生特大暴雨，板桥石漫滩水库溃坝，板桥，，京广铁路中断。

77.8内蒙古大暴雨，63.8海河大水危及到天津的安全。

81年长江上游洪水使四川、陕南遭受了严重损失，汉中机场变成一片汪洋,重现期只相当。同年黄河上游发生了50年一遇大水,上游龙羊峡正在施工。

83年汉江发生大水。汉江干流和支流洪峰同时在安康会合，3日暴雨仅为50年一遇而洪峰高达300年一遇，晚8点域门进水全城被淹。解放前1915年珠江发生200年一遇大洪水，珠江三角洲堤防崩溃，广州被淹7天7夜。31年长江淮河并流，共死亡22万人。32年松花江发生大洪水，哈尔滨市内水深达5米。38年花园口事件，国民党扒开黄河大堤，死亡89万人。干旱问题同样很严重，我国半数以上国土属干旱半干旱地区，年降雨量小于400mm，华北地区10年7旱。陕西关中地区曾出现连续十一年枯水年。天津1952-1968年来水90亿、海河入海51亿而70年代入海仅3.4亿，工业用水增加70倍、生活用水增加7倍。大量地下水开采使地面下沉。上海1921-1965年在150km2范围内下沉1.67米，最高达2.63米，66年后才基本得到控制。

解决水害的有效途径就是建水利工程，变水害为水利。长江三峡工程跨世纪工程、南水北调、各水库如小浪底、三门峡、天桥、青龙峡、盐锅峡、八盘峡、刘家峡、李家峡、渡、陈村、梅山、龙羊峡、大峡、黑三峡、葛洲坝、新安江、陆浑、丹江口、安康、石泉、乌江、佛子岭等，前年刚上马的云南小湾水电站，目前我们正在帮助设计的云南虎跳峡水电站，湖北江坪河水电站，以解决水资源的时空分布不均的问题。

中国万元工业产值用水量是发达国家的7倍，水重复利用率为35-40%，而发达国家达70-80%，节水是一个重要课题。

1975年8月，两千米长的板桥水库大坝被洪水掘根冲垮了三百七十多米。

小浪底水利工程新安江刘家峡三门峡葛洲坝富春江都江堰三峡可以说是万里长城又一巨大工程，坝高175米，装机1800万千瓦，回水长400—700公里，10个城市被完全淹没，8个城市（包括重庆）将受到影响，淹没万亩地，迁移100万人，各种迁移费35亿元，建设费130亿美元。防洪发电将取得巨大效益。

第二节问题的提出

就水文水利计算而言，目前采用数理统计方法来解决设计洪水、设计暴雨、设计年月径流以及设计库容、装机容量等问题。此类方法视水文变量为相互独立的随机变量，如Qm、WT、Xn、Q年等，此类方法或模型也称之为随机水文模型，但只是纯随机水文模型，模型中没有任何确定的成分。

20~30年代西方开始用这种纯随机模型来解决水文问题，最早是在1896年就有人将Qm作为独立随机变量，采用频率计算方法推求Qp。但这种方法求得的设计值往往误差很大，在实际运用中出现一系列问题。美国内布斯加州共和河1935年大洪水相当于万年一遇远远超过了过去用39年实测系列计算的设计标准。乌拉圭波德尔水电站用27年系列设计，而设计后59年特大洪水相当于50万年一遇。77.8内蒙相当万年2倍多。75.8林庄暴雨相当于千年一遇，83年安康洪水相当于300年一遇，龙羊峡相当于50年一遇。当然设计洪水算大了实际没出现这样大，问题还不算大，算小了出现大了，就会使建成的水工建筑物遭到破坏。

出现以上一系列的原因在于所用的实测样本系列太短，不足以求出精度较高的设计值。在水文统计中我们知道样本系列n愈短犯两类错误的概率愈大。

例：有一车苹果

原假设H0：苹果质量好

第一类错误：拒绝H0/H0为真拿一箱打开全是烂的，不买了，实际上其他全是好的。

第二类错误：接受H0/H0不真拿一箱打开全是好的，买回去其他全是烂的。要使P（拒绝H0/H0为真）P（接受H0/H0不真）同时减小，只有增加n。

增加样本容量的一般途径有历史调查洪水、插补延长。另外一条途径就是建立随机水文模型来模拟生成水文系列，来达到用传统计算方法不能达到的目的。

Q1,…,Qn实测ModelQ1’,Q2’,…Qn’模拟系列

1927年Suldel为确定库容分布生成1000年系列。

1984年陈文辉在华水学报上发表了“随机生成系列不能减少样本抽样误差”在全国引起了很大反响。但随机模拟方法对现行水文计算方法的研究无疑是非常有用的。如参数估计方法研究、插补延长方法的研究、洪水地区组成方法的研究、站网规划方面的研究、水文中长期预报、实时洪水预报、洪水保险风险分析等。

“随机生成系列不能减少样本抽样误差”x1,x2,…,x3000xp0(x1,x2,…,x30)1xp1(x1,x2,…,x30)2xp2(x1’,x2’,…,x1000’)1xp1’(x1’,x2’,…,x1000’)2xp2’………(x1,x2,…,x30)100xp100(x1’,x2’,…,x1000’)100xp100’总体

CvCs争议表明更接近，更小本课程专门研究模型，模型如何建立、参数如何估计、系列如何模拟等方面内容。

如自回归模型AR(P)：

1.模型如何确定，即为什么采用自回归模型AR(P)；2.模型参数如何确定,阶数P、φP,0、φP,1…φP,P即随机项εt的参数σε和Csε如何确定；3.水文系列的模拟。

水文现象是随时间而变化的过程，水文过程包含两种成份：一是确定成分、二是随机成分。

模拟水文过程的模型纯随机模型（不存在确定性因素和相依性）随机模型(确定性因素和随机因素共存)确定性模型（包含周期和非周期因素）平稳的（只有纯随机和相依成分参数不随时间变化）非平稳的（确定与平稳共存参数随时间变化）第三节随机水文模型模拟水文过程，可用三种数字模型：

纯随机模型：无确定性因素和相依性（如自相关系数ρ(1)=0…），水文计算中频率计算方法就是采用纯随机模型来研究水文变化的问题。只研究某事件的概率而不研究过程。如去年发生的洪水与今年无关等。确定性模型：不存在任何随机性因素，但包括周期和非周期成分，是一个确定的函数关系，知道去年过程就可确定今年水文过程。随机模型：以上两种极端情况不能真正反映水文的复杂过程，引入随机模型。两种因素在模型中均得到反映。第四节随机水文的发展

本学科是在近20年才成为一门完整学科。海森（Hazen）和塞德芬（Sudlen）早期用数理统计方法来分析径流序列，五十年代初赫斯特（Hurst）研究了径流和地球物理现象的长期序列。巴尔纳（Barnes）54年应用正态随机数表模拟径流。60年初塞姆思（Thomas）、费营（Fiering）和叶菲耶维奇（Yevjevich）考虑到水文现象的相依性而引用了马尔柯夫模型（自回归模型），并用以模拟年月径流序列，现在广泛用于水文预报、水库调节计算中。

从70年代起引用贝叶斯方法和卡尔曼滤波方法来估计模型参数，模型的参数亦随之变化，所以特别称这样的模型为变参数模型。

近来在广泛应用的线性正态模型基础上提出了非线性偏态模型，尚能反映水文变量P-Ⅲ型分布的特征。

60年代末产生了自回归滑动平均求和模型、解集模型、散粒噪声模型、分数高斯噪声模型、快速分数高斯噪声模型以及折线模型等。目前广泛用的是前两种模型，这种模型是数学模型。

总之，随机水文学的发展史从马尔柯夫模型到非马尔柯夫模型，从纯数学模型到有一定物理基础的模型，从定参数模型到变参数模型，从单变量模型到多变量模型，从线性正态模型到非线性偏态模型。随机水文学将有更大更快的发展。

第五节随机水文模型的应用（一）水文水利计算方面1.在系统分析计算中的应用如：年月径流序列的随机模拟随机模型年月径流系列系统分析与计算模拟序列系统响应

通过随机模型，借助统计试验法可获得大量的统计分析与计算所必需的模拟序列。黄河上游这样复杂的系统在进行水库联合优化调度时，作为系统输入多站模拟序列是必不可少的。其区间也通过模拟，不必用上下同频区间相应的方法求区间。龙刘盐八青2.

在插补序列中的应用

传统插补延长序列的方法，只是用回归方程而没有考虑随机误差，是插补变量的方差偏小。用来估计

插补后y的Cv值偏小，求的设计值偏小，m值偏大，偏小愈甚，ρ愈小插补效果愈差。

xαE(y/xα)xf(y/x)y假定、ρ0、Cvx0、Csx0、、Cvx0、Csx0、yp0x31…x40x1x2…x30y1y2…y30矩法、Cvy’、Csy’、yp’求b0、b建立y1…y30y31…y40矩法、Cvy’’、Csy’’、yp’’………100组自回归模型插补延长的统计试验研究

均值愈接近愈好，均方差愈小愈好。计算不同种参数用不同的m。ρ>0.8m<nn为样本长m插补延长长度水文参数估计方法的统计试验研究

假定P-Ⅲ分布EX0、Cv0、Cs0、XP0矩法xp’极大似然法xp’’适线法xp’’’各种经验频率公式不偏的经验频率公式(x1,…,x30)1(x1,…,x30)2………(x1,…,x30)100均值愈接近yP0愈好，均方差愈小愈好

（二）水文预报方面1.预报误差的处理

2.随机模型预报3.卡尔曼滤波实时预报

（三）在水文测验和站网布设方面

划为三块每块随机抽一个，三个雨量站时段平均雨量与所有雨量站时段平均雨量作比较，后者为总体面雨量，如差别大多分几块再随机抽重复，最后确定雨量站个数和位置，以误差控制。

沂沭泗流域雨量站网优化规划采用0-1整数规划优选雨量站个数和位置，是优选出雨量站时段（40）雨量与总体相差以内。雨量站为xi(i=1,2,…,n)xi=0不选；xi=1选中

目标函数∶min=

约束条件：时段面雨量约束i=1,2,…,40——时段aij——i号站j时段雨量xi——i号站（取值0或1）

上海黄浦江建闸方案研究（1998-2000年）潮位H0由天文潮HP和风暴潮ΔH组成，两者独立。

H0=HP+ΔH

HP样本HP1,HP2,…,HPn→求参数、Cv、Cs→生成HP1

，HP2

，…，HP1000000

ΔH=H0-HP样本ΔH1,ΔH2,…,ΔHn→求参数、Cv、Cs→生成ΔH1，ΔH2，…，ΔH1000000相加可得H01

，H02

，…，H01000000，可得百万个H0，由大到小排在第100个为万分之一设计值。

课题中的应用黄浦江吴淞高潮位、上游太湖来水、上海市区间来水碰头的概率三维联合分布，随机模拟系列求课题中的应用具体做法：

求三维联合分布，随机模拟一百万组系列，求三组碰头的概率(1)对H吴年最大值相对应的太湖30天雨量X太，相应上海一日雨量X区三系列分别取对数（lnH0吴、lnX区、lnX太），使其近似为正态分布；

(2)求均值和均方差矩阵课题中的应用(3)由样本推求以下系数矩阵其中(4)则符合相应参数的正态随机变量可由以下式子生成：y1、y2、y3为服从标准正态分布的随机数N(0,1)课题中的应用(5)分别对H0吴、X太、X区进行频率计算,求出设计值H0吴P、X太P、X区P参数(6)模拟一百万组H0吴(x1)、X太(x2)、X区(x3)(7)统计同时发生的组数(8)求碰头的概率P=发生组数/一百万长江口、杭州湾上游耒水、潮位和波浪相遭遇的概率

求三维联合分布，随机模拟一百万组系列，求三项均满足上述的组数，除一百万即得均满足上述的概率暴雨过程的随机模拟方法

（云南虎跳峡水电站）

（1）根据1962-2000年39年年最大15日逐日雨量序列，采用适线法估计出各日的参数和相关系数∶EXi、Cvxi、Csxi、ρij（2）由各日的参数，根据（7）—（11）式估计各自的ai，bi，σi,rij七个参数

（3）生成15个[0,1]均匀分布的随机数u1、u2、u3、u4………u15（4）利用下式生成独立标准正态分布的随机数n1、n2、n3………n15

（5）利用Y=AZ其中A=（aij），（11）得到y1=a11z1y2=a21z1+a22z2…………（12）

yn=an1z1+an2z2+……+annzn

其中∶n=15（6）由xi=exp（σiyi-bi）+aij；i=1～n（13）得到（x1,x2,……，xn）就是具有给定参数ai,bi,σi及（或EXi,Cvxi，Csxi,ρij）的对数正态变量（X1,X2,……,Xn）的抽样值。重复上述步骤1000次，即可得到1000组，每组15日的暴雨过程。沪宁高速加宽水文计算平原河网地区我们采用的推求设计水位的方法

H桥P=H水文站P+（H桥99-H水文站99）用二维随机模拟方法论证我们采用的方法的合理性第一节随机水文过程的概念t1流量x(t)X1(t)X2(t)Xn(t)X1(t1)X2(t1)Xn(t1)时间t例：某水文站各年的日流量过程可用图来表示，有几年就有几条日流量过程，我们称一年的日流量过程，为一次试验的结果（一个现实或一个样本函数）。一次试验的结果就是日流量随时间t而变化的函数。其函数各次不同，且事先无法确定（即每年日流量过程各次不同，事先无法确定），对于多次试验的结果是一族时间t的函数，我们称这族时间t的函数为随机函数。而每次试验结果，即族中每一个函数称为随机函数的一个现实或样本函数。可见，随机函数就是所有现实或样本函数的集合。第一节随机水文过程的概念第一节随机水文过程的概念

当随机函数随时间t连续的取有限区间内的值时，称此随机函数为随机过程。而随机函数随时间t取离散值称为随机序列（现在统称随机过程）。如月径流过程是以月为参数的12维随机变量，日径流过程是以日为参数的365维随机变量。一般随机过程用x(t)表示，各个现实用xi(t)表示。

第一节随机水文过程的概念对于每一个固定时刻，如t=t1（5月15日）x(t1)是一个随机变量，而x1(t1),x2(t1),…xn(t1)是随机变量X(t1)在每次试验中的取值，我们称X(t1)是随机过程X(t)在时刻t=t1的截口或称为随机过程X(t)在时刻t=t1时的状态。

第二节随机过程的分布函数

X(t)是一个随机过程，对于每一个固定的t，X(t)是一个随机变量，它的分布一般与t有关，记为随机过程的分布。

称为随机过程X(t)的一维分布函数。则称为随机过程X(t)的一维概率密度。如t=2，则X(t)是表示2月径流的随机变量，则表示2月份径流小于等于50m3/s的概率。表示2月份月径流的概率密度（为一维的）。

第二节随机过程的分布函数同理二维或多维表示为：

第三节随机水文过程的数字特征一、均值函数

X(t)tnXm(t)X2(t)X1(t)μ(t)随机过程在固定t时X(t)为一随机变量，它的均值或数学期望一般与t有关。称μ(t)为随机过程X(t)的均值函数（简称均值），它是普通t的函数。

第三节随机水文过程的数字特征

二、方差函数

称为随机过程X(t)的方差函数（方差），它也是个t的普通函数。σ(t)称为随机过程X(t)的均方差或标准差。

第三节随机水文过程的数字特征三、自协方差函数

为随机过程X(t)的自协方差函数（或称协方差函数）。如t1=2,t2=4表示了二月份与四月份月径流的协方差。

第三节随机水文过程的数字特征四、自相关系数

称为随机过程X(t)的标准化协方差函数或自相关系数。表示了随机过程任意两个不同时刻截口之间的线性相依关系程度。

第四节

自相关和互相关系数的估计

一、自相关分析

由实测样本xt（其容量为n）估计的自相关系数rk

式中：

——样本协方差

——样本均方差

——样本均方差

第四节

自相关和互相关系数的估计式中均值

第四节

自相关和互相关系数的估计代入整理后：

第四节

自相关和互相关系数的估计当n较大时，和都用

代替；同时当n增大及k较小时，因此有滞时R=0,1,2,…,m，，当n>50时，可取，常取m在n/10左右；当n<50时，取m在n/4左右某数值；有的取m<n-10，即参加计算的至少在10项以上。

第四节

自相关和互相关系数的估计Kr(k)乌江武隆站年径流自相关图

第四节

自相关和互相关系数的估计二、互相关分析（两个随机过程）

由实测序列xt和yt（样本容量为n）估计的互相关函数（互相关系数）为rk(x,y)

式中：——用样本估计的互协方差

——样本估计的xt和yt的均方差

第四节

自相关和互相关系数的估计式中均值

第四节

自相关和互相关系数的估计代入最后可得出前式相仿的公式

第四节

自相关和互相关系数的估计n较大时，当k=0时互相关系数r0(x,y)是普通相关系数。注意：k与-k求出的rk(x,y)与r-k(x,y)是不相同的。

第四节

自相关和互相关系数的估计月rk(x,y)程家河站月降水量（流域平均值）和月径流量互相关图

k=10表示月降水量（如2月）与滞时为10个月（12月）月径流量互相关系数；k=-10表示月降水量（如2月）与滞时为-10个月（去年4月）月径流量互相关系数

第五节随机独立性的假设检验一、持续性检验二、趋势性检验

三、随机特性检验

1.序列自相关系数检验

2.序列秩号相关系数检验

1.秩号和检验

2.秩号与序号相关系数检验

1.游程检验

2.转折点检验

一、持续性检验(1)序列自相关系数检验

原假设H0：一阶自相关系数

统计量：

式中：N～样本长度；

r1～样本一阶相关系数

统计量t符合自由度为的t分布如则接受，认为样本自相关系数为0反之拒绝，认为不等于0

检验方法一、持续性检验检验方法(2)序列秩号相关系数检验

所谓秩号是对原数列按大小顺序排列后的序号

原假设H0：一阶秩号相关系数

统计量：

式中：s1～秩号一阶相关系数

统计量t符合自由度的t分布，检验方法同前二、趋势性检验检验方法(1)秩号和检验

按时间顺序将项序列分成相等二段，其项数分别为和（n1+n2=N）再统计每段内的秩号总和。如样本具有趋势性，则秩号相近的事件将集中在一起出现，从而和差别较大，反之较小。

原假设H0：R1和R2差别不显著

统计量：

式中：(1)秩号和检验检验方法可以证明U符号正态分布，其参数

同样在置信限α给定的情况下可做出接受和拒绝H0的决定

(2)秩号与序号相关系数检验检验方法序列从大到小排列后所得到的秩号和反映时间序列的序号之间的相关系数反映了序列的趋势性。

原假设H0：秩号与序号相关系数等于0

统计量：

式中：

统计量t符合自由度df=N-2的t分布。与上法相同，对各站样本可作出接受或拒绝H0的判断

三、随机特性检验检验方法(1)游程检验

将时间序列从中间分成项数相等的A、B两项(n1和n2)，再全部按大小排列，并标出所属的组号，连续出现同一组号的元素为一个游程，检验全部序列总的游程数K，若样本中大值或小值成组出现，则K小，表明样本存在持续性；反之则K大，表明样本存在周期性；只有当游程数K接近其期望值EK时随机性才较好。

原假设

统计量

Ki——游程数

(1)游程检验检验方法K符合正态分布，参数如下：

将统计量值K标准化后，在给定的置信限α下，可作出拒绝或接受原假设的决定。

(2)转折点检验

检验方法时间序列波动的极大值（Xi-1<Xi>Xi+1）或极小值（Xi-1>Xi<Xi+1）总数可作为统计量进行检验。当统计量P较大时表示序列呈锯齿出现，则样本随机性差。当P较小时，序列连续上升或下降具有持续性。因此，P接近其期望值EP，则以为样本随机性较好。

原假设

统计量值

式中：Pi—第i项转折点数

第一节概述

水文序列包括确定成分和随机成分，后者包括相依成分和纯随机成分。欲随机模拟水文序列，需先随机模拟序列中的纯随机成分，再依时序将其叠加在其它成分上。

第一节概述例：随机模拟某站符合P-Ⅲ分布的年最大流量Qm序列（即纯随机序列）PiQmiQmPP-ⅢEXCsCv1.由样本先估计总体分布Qm-P频率曲线

2.随机模拟Pi(i=1,2,…)，由Qm-P关系查出，由此可见纯随机模拟方法之一，就是要解决如何模拟Pi，以及由Pi如何转换为指定分布序列两个问题。

第一节概述统计试验法（也叫作蒙特卡洛法（M-K））

首先模拟[0，1]区间上均匀分布的纯随机序列（简称随机数），即相当Pi；再转换为指定分布的纯随机序列。

第二节随机数的模拟

目前均采用乘同余法，乘同余法递推公式为：n=0,1,…

式中：x0——为初值λ——为乘子

M——为模

这些都是非负整数，而且λ<M，上式表示为xn+1是λxn被M整除后的余数，而。

就是[0，1]区间上均匀分布的随机数。

第二节随机数的模拟例：设x0=1，λ=7，M=103则有

按这公式模拟的序列到一定长度后，最终要出现重复循环，循环的长度称为周期L。关键选择参数λ和M使周期足够长，通称模拟随机数为“伪随机数”。

第三节指定分布纯随机序列的模拟

一、正态分布纯随机序列的模拟

ξ1、ξ2为相互独立的标准正态分布N（0，1）的随机数

则正态分布的伪随机数

第三节指定分布纯随机序列的模拟二、对数正态分布纯随机序列的模拟

假定随机变量x为对数正态分布，y为正态分布，它们之间有如下关系：

第三节指定分布纯随机序列的模拟式中：，Cv，Cs为的参数可由样本用矩法估计

三、皮尔逊Ⅲ型分布纯随机序列的模拟模拟均匀分布随机数u1,u2,…,uα’α’=[α][α]小于等于α的最大整数计算模拟伪随机数uα’+1uα’+2uα’+3计算计算计算xt是否舍选法

舍选法t=1,2,…即生成t个P-Ⅲ型分布纯随机数

在模拟时必须否则舍去，重新生成一对伪随机数重新计算再判断，直到满足要求为止即可生成一个P-Ⅲ型分布随机数。

式中：

为第三节指定分布纯随机序列的模拟P-Ⅲ型分布纯随机序列的模拟，还可以利用Φ值表。模型中的参数一般均采用矩法估计，此算法计算简便，适用于计算机大量计算。

先抽ut（即P值）由Cs查出Φt

为P-Ⅲ型分布随机数

用舍选法模拟一个xt需用多个伪随机数，特别Cs较小时，可用下列公式模拟P-Ⅲ型分布纯随机序列

第三节指定分布纯随机序列的模拟Φt为标准化P-Ⅲ型分布随机数（即离均系数）

其中ξt为标准正态分布

xt——为P-Ⅲ型分布

第三节指定分布纯随机序列的模拟…

…

其中y1…yn为标准正态分布随机数，按前面的公式进行模拟

四、多维正态分布随机变量的模拟四、多维正态分布随机变量的模拟四、多维正态分布随机变量的模拟参数可用矩法估计

五、多维对数正态随机数模拟1.先由样本估计参数

Exi、Cvxi、Csxi(i为维数)

31312424úúûùêêëé++--úúûùêêëé++=iiiiiCsxCsxCsxCsxh五、多维对数正态随机数模拟2.生成[0,1]均匀分布随机数ui，i=1~n（二维就模拟两个）

3.由ui生成N（0，1）标准正态分布随机数Zi，i=1~n

4.由公式生成yi系列(i=1~n)…

…

其中

5.对数正态分布随机数

第四节纯随机序列模拟的应用例：假定X、Y均服从对数正态分布，本课题研究由X（参证站）来插补Y（设计站），插补效果的研究。即ρxy要多大（是否>0.8）？插补长度m不宜超过多长（是否m<n）？

xynm假定XY总体参数EX0Cvx0Csx0EY0Cvy0Csy0ρxy0求设计值真值yp0p=1%p=1‰随机模拟x1,x2,…,x30y1,y2,…,y30估计样本参数EXCvxCsx

EYCvyCsxrxy0求设计值（未插补）yp=EY(CvYΦP+1)随机模拟x31’…x40’m=10建立回归方程y=b0+bx插补yy31’…y40’由样本y1…y30,y31’…y40’用矩法求EY’Cvy’Csy’求设计值yp’=EY’(CvYΦP+1)(插补后)重复100次k=k+1k<100分别求设计值均值和方差（插补前后）

K≥100上海黄浦江防洪问题苏通大桥

课题中的应用(一)应用H0HHPΔHt天文潮和增水的随机组合

HP——天文潮可由调和分析求出

ΔH——由实测H0-HP得到

应用求H0P

(1)可由年最大H0经频率计算求出(2)通过如下随机模拟方法求出

HP1HP2…HPn

CvHPCsHPHP1…HP100000ΔH1ΔH2…ΔHnΔH1…ΔH100000

CvΔHCsΔH+H01…H0100000H0P(P=0.01%0.1%1%)应用课题中的应用(二)推求苏通大桥涨落潮流量

应用杨林徐六泾天生港苏通大桥长江青龙港桥位及各测站分布示意图

应用由徐六泾1984-2001年实测和通过1997-2001年整点潮位资料采用非恒定流计算的桥址潮流量资料，得如下表：

徐六泾站潮位要素与潮流量关系

（1）（2）推求苏通大桥涨落潮流量应用Z低ΔZ涨ΔZ落由ΔZ徐涨序列求参数随机模拟代入（1）求Q涨P

由ΔZ徐落和Z天低序列用二维对数正态分布求参数随机模拟代入（2）求Q落P

从而通过水位要素求设计潮流量

推求苏通大桥涨落潮流量数值积分法

通过数值积分求Q落P

假定Q落，则由（2）式可得出：

应用数值积分法应用ΔZ徐落a2a1Z天低假定ΔZ徐涨(y)和Z天低(x)为两维对数正态分布，则a1为Z天低最小取值，a2为ΔZ徐落最小取值，由积分区域可求出相应的Q落的P

应用数值积分法则可得到不同的Q落和不同的P求P=1/300、1/100、1/50相应的Q落P

数值积分法应用其中：二维对数正态联合分布

安全修正值问题

Cs=4.03.53.0BP其中：a在0.7左右

由上式可算出B，再假定不同参数不同的P，再算可点出上图。

流程图见下页：假定P-ⅢCv0Cs0Xp0

(x1,x2,…,x30)1(x1,x2,…,x30)100(CvCsXp)1(CvCsXp)100流程图……关于偏态分布的曲线插补问题

xαEy/xα=b0+bxαxf(y/x)y正态分布：

偏态分布：为非线性导出对数正态曲线的插补公式，并用随机模拟进行，证明曲线插补更优

第一节平稳随机过程

定义：过程的统计特性随时间t的变化而改变称为非平稳随机过程，反之称为平稳随机过程。

平稳随机过程的特点：过程的统计特性不随时间的变化而变化。（如年径流过程）

第二节线性平稳模型的一般概念

一、自回归模型模型AR(P)一般形式：

式中：μ——x的均值

——自回归系数

εt——一般假定为独立随机序列（白噪声序列）且假定服从均值为零、方差为的正态分布N(0,σε)p——自回归阶数第二节线性平稳模型的一般概念二、自回归滑动平均模型ARMA(p，q)

一般形式：

式中：Q1…Qq——滑动平均系数

q——滑动平均阶数

第二节线性平稳模型的一般概念当q=0时，ARMA(p,q)模型为AR(p)模型

当p=0时，为MA(q)模型，称之为滑动平均模型

以上两个模型称之为线性平稳模型

第一节一般自回归模型

一般自回归模型表示为：

又可表示为：

式中：

第一节一般自回归模型自回归模型用AR(p)来表示，有时称为滞时为p的自回归模型和滞时为p的马尔柯夫模型。

ΦP,0

ΦP,1…ΦP,P——回归系数

μ——序列均值xt——自回归序列简称AR(p)序列P——阶数εt——独立随机变量在水文中通常假定为P-Ⅲ分布，均值为零，其参数仅有两个即标准差σε与偏态系数，没有Cvε。

第一节一般自回归模型将上述中心化得到尤尔—沃尔克（Yule-Walker）方程：

式中：

建立上述模型主要解决：1）模型最大阶数P；2）参数估计

第二节自回归模型的参数估计

自回归模型参数包括：自回归系数ΦP,0,ΦP,1,…,ΦP,P,μ,σε和Csε

自回归系数的均值与方差可由矩法来估计

第二节自回归模型的参数估计一、回归系数的求解方法

递推解法

当P=1时

当P=2时

递推解法当P=3时

递推解法当P=K时

其中：j=1,2,…,k-1二、随机项εt的统计参数

式中：r1,r2,…,rp为系列xi的1阶,2阶,…,P阶相关系数，可由实测序列采用矩法估计

二、随机项εt的统计参数三、模型的识别—最大自相关阶数P的选定

实用选定最大阶数P的方法是根据一般经验公式左右

对于AR(P)序列而言，偏相关系数Φk,k是P步截尾。N为径流系列的样本系列长度，因此，从理论上说如当k>p后，其中，则可选定自回归的最大阶数P。

尤尔—沃尔克（Yule-Walker）方程组中各阶系数中Φ1,1,Φ2,2,Φ3,3,…,Φp,p叫做偏相关系数。可以证明，当径流系列为P阶自回归模型时，则当k<p时Φkk≠0

呈截尾状。

四、自回归模型AR(P)的计算过程

(一)自回归模型AR（P）的一般形式

ΦP,0ΦP,1…ΦP,P——自回归系数

εt——白噪声服从P-Ⅲ分布独立随机变量

四、自回归模型AR(P)的计算过程(二)计算实测系列的统计参数

四、自回归模型AR(P)的计算过程(三)按公式选定最大阶数P

(四)用递推公式计算自回归系数同前递推解法；当P=K时

其中：j=1,2,…,k-1

四、自回归模型AR(P)的计算过程(五)计算ΦP，0

(六)计算εt的参数σε和Csε

例

江西鹤洲站资料见课本P31表5-1，建立自回归模型AR（2），并生成10年年径流系列。

例解：一、一般形式

二、计算样本参数

三、确定阶数P=2

例四、计算回归系数

例当P=1时

当P=2时

例五、计算Φ2，0

六、计算εt的参数

例七、模型的确定

八、由模型模拟10年径流系列

采用

例1.先模拟10个εt

2.模拟年径流序列

徐六泾年最大落潮流量频率计算因苏通长江公路大桥河床冲刷模型计算需要，分析连续五年出现大于20年一遇徐六泾落潮流量的概率，本次计算建立了徐六泾年最大落潮流量的自回归模型—AR(P)模型，模拟了500万年年最大落潮流量，统计连续五年出现大于20年一遇落潮流量的次数，并计算分析出现的频率。计算说明

1模型的确定

自回归模型—AR（P）模型是水文常用来模拟水文系列的模型，模型形式如下式：2模型阶数的确定

计算说明

阶数的确定可以采用AIC准则或通过自相关系数的分析做出。本次采用对自相关系数的分析确定阶数P。计算相关系数的方程如下：徐六泾年最大落潮流量自相关系数表

3模型参数的计算

计算说明

徐六泾年最大落潮流量序列和参数见表：3.1回归系数的确定

3模型参数的计算计算说明

3.2随机变量ε的参数估计

3.3模型的建立

徐六泾年最大落潮流量AR（1）模型如下∶

4随机模拟

计算说明

为统计连续五年出现大于20年一遇落潮流量的频率，本次采用上述模型模拟500万年徐六泾年最大落潮流量，统计连续五年出现大于20年一遇落潮流量的次数，计算出现的频率。

4.1随机项的生成计算说明

先生成0-1均匀分布的伪随机数ui，由ui和并用Φ值表采用三点插值查出对应ui(pi)，Csε(pi)的Φ值，由下式模拟出皮尔逊Ⅲ型分布随机项。

4.2徐六泾年最大落潮流量序列的模拟

计算说明

由建立的徐六泾站年最大落潮流量序列自回归模型AR(1)：

由2000年徐六泾站最大落潮流量代入上页公式中，再生成一个随机数代入上页公式中，则可算出模拟2001年的年最大落潮流量。如此反复计算500万次，得到模拟500万年徐六泾年最大落潮流量序列。

主要内容：第一节滑动平均模型第二节自回归滑动平均模型

第三节模型的识别和检验

第一节滑动平均模型一、一般滑动平均模型MA(q)滑动平均模型的一般形式MA(q)其中心化形式为：

式中：θq,1,θq,2,…,θq,q——各阶权重又称权重系数是待求的参数

εt——白噪声系列，在水文上仍假定为P-Ⅲ分布，均值为零，其参数仅有两个即标准差σε与偏态系数，没有Cvε

二、参数估计式中：样本序列的和r1,r2,…,rq（各阶自相关系数）均可由样本计算出。

由一般形式经推导可得到如下方程组：

三、的计算一般形式经变形可得：

其中：x1’,x2’,…,xn’为已知中心化径流系列θq,1,θq,2,…,θq,q，已由上面解出

…

……

则

三、的计算（2）

利用这样计算得到εt系列，系列ε1，ε2，…，εn按照矩法估计：

四、一阶滑动平均模型以中心化变量表示的MA(1)模型形式如下：可由上述方法计算得到

五、二阶滑动平均模型以中心化变量表示的MA(2)模型形式如下：

其中

当r2>0时取+，当r2<0时取-六、例：

同前，江西塞塘站有15年（1964-1978年）年径流资料，试建立MA（1）模型，并生成10年符合P—Ⅲ型分布年径流系列。

例计算步骤1．MA（1）模型的一般形式2．计算样本统计参数同前例计算步骤3．用公式计算参数例计算步骤4．由实测样本经中心化后分别计算εt值……

例计算步骤5．计算Csε6．模型的确定例计算步骤7．由模型模拟10项新系列(1)先生成0~1均匀分布伪随机数10个作为P值

(2)由上述计算Csε查P—Ⅲ离均系数Φ值表，得到相应的Φ值

(3)代入上面公式可生成xt+1

同理生成xt+2,…,xt+10

计算结束第二节

自回归滑动平均模型一、一般自回归滑动平均模型一般自回归滑动平均模型ARMA(p,q)，若以中心化变量xt’来表示，则：

当时，即q=0时为自回归模型AR(p)

当时，即p=0时为自回归模型MA(q)

二、ARMA(1,1)模型的参数估计方法

经推算可得到如下方程：

方程中

二、ARMA(1,1)模型的参数估计方法则前式可由上述三个用实测样本估计的值算出Φ1,1，θ1,1，σε2

可采用图解法或迭代法求解θ1,1，σε2

三、Csε的计算

将一般形式变形为：

仿照MA(q)的计算εt系列的方法，由实测的系列x1’,x2’,…,xn’（中心化后的）及上面已计算出的Φ1,1,θ1,1，利用上式计算出ε1，ε2，，εn系列，再利用Csε的矩法公式：

计算Csε

第三节模型的识别和检验

对给定的水文序列，应选什么样的模型，是AR(p)，MA(q)还是ARMA(p,q)？若选出一种模型，例如AR(p)模型，又如何确定阶数P。这样的问题在时间序列分析中叫做模型类型的选择及其形式的确定，可以统称为模型的识别。通过识别确定模型，在由实测系列估计出模型中的参数，我们就初步得到了表征其随机变化的模型。对这种模型，尚需作进一步检验，以验证是否符合模型的一些基本假定。

一、模型初步认别的原理和方法

在模型初步认别中，主要的依据是样本序列的自相关函数和偏相关函数。主要特点归纳如下：

(1)AR(p)序列

自相关函数ρk，随滞时k的增大逐步变小，自相关图呈拖尾状。它的偏相关函数Φk,k呈截尾状，单调或波动衰减趋向于零。在k=p处出现一个截止点，即在时，当k>p时Φk,k=0。

(1)AR(p)序列（如下图所示）(2)MA(q)序列

自相关函数ρk呈截尾状，k=q时出现一个截尾点，即在时；当k>q时ρk=0。相反，它的偏相关函数随阶数的增加而逐渐变小，呈拖尾状。单调或波动衰减趋向于零。

(3)ARMA(p,q)序列自相关函数和偏相关函数都没有截尾点，均以拖尾状而逐渐变小，趋向于零。

缺陷：但由样本序列去估计自相关函数ρk和偏相关函数Φk,k抽样误差较大，难以直观判断，必须进行统计推断。例如，即使是出自AR(p)模型的序列，由于抽样的缘故在k>p时由样本序列估计出的也不全会为零，而在零上下起伏，而用统计检验方法，推断出它们和零的差异并不显著，则可认为截尾点是k=p，从而推断序列为AR(p)序列。下面分别叙述MA(q)，AR(p)和ARMA(p,q)模型识别的具体方法：

1．MA(q)模型由于k>q时ρk=0，rk在k>q时渐进服从的正态分布（其中为rk的方差D(rk)），根据正态分布的性质：1．MA(q)模型先假定q0如k<q0，rk应明显不为零（k=1,2,…,q0）而当q=q0时，k分别取q0+1，q0+2，…，q0+M（M可取左右）代入上式不等式，如果满足上式不等式，则统计满足个数占总数的比例，如果占总数M的68.3%或95%左右则可判断rk在q0处截尾，并初步认为序列是MA(q0)序列。2．AR(p)模型对AR(p)序列，当k>p时Φk,k=0；当k>p时将渐近服从正态分布，为的方差D()。

同样找到p0，在p0前（k取1,2,…,p0）明显不为零，而p=p0时，k分别取p0+1,p0+2,…,p0+M（M可取左右）代入上式不等式统计满足上式不等式的个数占总数M的比例，如在68.3%或95.0%左右可判断在p0处截尾，并初步认为序列是AR(p0)序列。3．ARMA(p,q)模型对于这种混合模型，自相关函数和偏相关函数均无截尾的性质，认别较困难，一般识别p,q的方法可以从低阶到高阶逐个尝试，(p,q)分别取(1,1)，(1,2)，(2,1)，(2,2)然后分别进行参数估计，定出模型在检验这个模型是否被接受。若被接受，则模型确定；否则重新调整(p,q)直到接受为止。

二、模型识别的AIC准则日本的赤池（Akaike）对于ARMA(p,q)模型中阶数p和q的确定提出了AIC准则：

式中：n为序列的长度；ln为自然对数，为残差的方差

q=0为AR(p)模型p=0为MA(q)模型使AIC达到最小值的模型被认为是可以接受的好模型。

对AR(p)模型，按AIC准则识别步骤：

（1）计算样本序列的自相关函数rk（2）按递推计算自回归系数Φp,1,Φp,2,…,Φp,p

按AIC准则识别步骤：（3）计算残差方差

（4）计算AIC

（5）根据不同P计算出相应的AIC(P)使AIC达到最小的P为我们所求

例

有序列n=45采用AR(p)模型，分别计算得到AIC(k)值，如下：试判断AR(p)的阶数p

p=14

三、模型的检验原假设H0：εt序列相互独立

统计量

M——最大滞时，一般m取n/4左右

Q服从自由度为（m-p-q）的Χ2分布

由样本ε1,ε2,…,εm计算出Q值，取α值由Χ2表查出Χα2，若则接受原假设，认为εt为独立序列的假设成立，反之不成立。

例

某站有44年径流序列，已识别出AR（1）模型。由样本序列，通过该模型算得εt序列，进而算出它的自相关关系：r1(ε)=-0.13，r2(ε)=0.14，r3(ε)=0.14，r4(ε)=-0.09，r5(ε)=0.04，r6(ε)=0.18，r7(ε)=-0.18，r8(ε)=0.05，r9(ε)=0.02，r10(ε)=-0.09。