径流随机模拟

1、径流随机模拟由于实际水文资料往往比较短，难于满足实际水文程作随机模拟。这种随机模拟的目的 之一在于充分利用是用来延长资料长度。当所建模型及参数准确时，这种年月径流随机模拟 方法，对于多站及更深入的随机模拟一、随机过程基本知识( 一) 随机过程和时间序列的定义在实际问题中，常涉及试验过程随某个参变量的变的流量、水位是随时间变化的随机变 量，气温是随时间这种随机变量为随机函数，并称以时间为参数的随机函数为随机过程，记 为 (t),t T，T是t 变化的范围。随机过程在一次试验或观测中所得结果，称为随机过程的一个实现。若时间参变量 T 是连续时刻的集合，则称这种随创机过程为连续参数随机过程，如水位

2、过程、流量过程等。若时间参变量 T 是程为离散时刻的集合，则称这种随机过离散参数随机 过程，也称为随机序列或时间序列。如年、月径流程，年最大流量过程都是时间序列，也称 水文时间序列。( 二 ) 随机过程的数字特征随机过程 (t) 在任一固定时刻的状态是随机变量，因此可按与前述随机变量同样的方法定义 随机过程的数学期望和方差。定义如下数学期望 (t) E (t) (811)12)方差2(t) E (t) (t)2(8为了规划随机变量两个不同时刻状态间关系的密切程度，可定义随机变量的自相关函数 为R(t1,t2)E (t1) (t1) (t2) (t2)(t1) (t2) 8-13 )( 三) 随

3、权过程基本分类l 。按统计性质的稳定性分类 按随机过程的统计性质是否随时间而变化，可分成平稳和非平稳过程。若随机过程统计 数字特征不随时间的平移而变化，则称为平稳过程，否则为非平稳过程。2 按不同时刻状态间的关系分类 可分成独立过程和马尔柯夫过程。若过程各状态相互独立，则称为独立随机过程。在非独立随机过程中，最重要的一类是马尔柯夫过程，其特点是t n时刻状态只与 tn 1时刻有关，而与 tn 1 以前各时刻无关。以上各个概念的严格定义，可参见有关文献。二、径流随机模拟一船步骤图 87 给出了单站年径流随机模拟的一般步骤：时间序列组成分析；模型的建立； 序列的生成；模型及生成系列的检验。三、水文

4、时间序列的组成分析水文序列 Qt 一殷可按下式表示：Qt Tt Ct Pt St （8-14 ） Tt，Ct，Pt，St分别为趋势项、跳跃项、周期项和随机项。当水文序列 Qt 中不含 Tt，Ct，Pt等确定性成分时， 如Qt St ，即仅包括随机成分的序列。 对年径流序列而言，这种情况是出较常见的。但月径流序列因存在明显的年周期，所以不是 仅包括随机成分的序列。趋势项指的是水文变量的统计参数 （主要是均值 ）的长期系统性升降， 如图 88 即存在明 显趋势 （ 增长 ） 。一般是由于气候因子或下垫面因子逐步改变而引起的缓慢变化。对实测水文序列，可用假设检验或滑动平均的方法查明是否存在趋势。若存

5、在趋势，呈 线性变化时，常用线性方程拟合，然后从序列中将趋势滤掉。跳跃项是指水文序列急剧变化的一种形式， 当水文序列从一种状态过渡到另一种状态时 表现出来。见图 8-9 。跳跃是人为或自然原因造成的。如建库后面积增大，蒸发量等损失增 加，有可能出现跳跃， 并反映在年径流序列的均值等参数。 突变可看作跳跃的一种特殊情况， 如由于地震塌方，拦截江河，形成水库以后溃坝，这样引起流量的突变，随着临时水坝的冲 毁，又恢复到原来状态。跳跃是否存在于序列中， 多用分割样本的方法检验。 若存在较显著的跳跃成分， 应从 序列中排除掉，使得剩余序列具有原始状态或一致条件。周期项 （含近似周期 ）是由于天体运动的周

6、期性影响造成的， 如地球公转、 自转引起以 年和日为周期的变化，以及太阳黑子活动引起的旱涝多年变化。通常可用谐波分析的方法析 出，再从序列中滤掉。随机成分是由于不规则及随机振荡引起的， 一般由相依成分和纯随机成分组成。 严格 地讲，几乎所有水文变量应是非平稳过程。不过除了人为影响及自然灾变外，水文环境的变 化在数十年或几百年期间都相当小，因此，从实用观点，常把水文序列中的随机成分看作平 稳过程。图 88 序列趋势变化示意图图 89 跳跃成分及包含跳跃成分的合成序列四、单站年月径流随机模型的建立对年月径流序列建立随机模型一般是对原始年月径流序列排除趋势跳跃等确定性成分后 的随机成分而建立的模型。

7、设用于建模的年径流序列为Q1， Q2，Qn。径流序列随机模型是实测时间序列一种概化的数学式子。目前常见的随机模型有：线性 平稳模型、非线性平稳模型，马尔柯夫模型及非马尔柯夫模型等。对年径流序列而言，常采 用线性平稳模型或马尔柯夫模型。建立随机模型的一般步骤为：选择模型；确定阶数；估计模型参数(一)单站年径流随机模型的建立通常采用线性自回归模型，即马尔柯夫模型1 线性自回归模型的一船形式Qt Qp,1(Qt 1 Q)p,2(Qt 2 Q)p,p(Qt p Q) t (8-15 )Qt 第 t 年的年径流量， t=1 ，2，常称式( 8-15 )的 Qt 为自回归系列。Q Qt 序列的平均值；p,

8、1， p, p 自回归系数或偏相关系数，反映Qt 在时间上相依性大小；t 模型残差项，纯随机成分， t 与 Qt 1， Qt ，无关且是独立随机变量，其均值为0，方差为 2t 。由于 2t与 Qt 的方差 Q2有确定关系， 因此，一般自回归模型中参数有： Q 、 Q和 p,1 ，p,p，共 p十 2个参数。该模型说明第 t 年年径流量仅依赖于第 t 一 1 年、第 t 一 2 年、第 t 一 p 年的年径流量和一个纯随机变量 t 。若令 yt Qt Q ，则式 （8 15） 变为：ytp,1yt 1p,2yt 2p,p yt p t (8 16) 式（8 16） 是中心化变量表示的自回归模型。

9、2模型参数的估计Q 1Qini1n(Qt Q) 2 i1?QsQr11 rp 2r2r1rp 3rp 1r1rp 2r21rp8-19 )其中 K阶样本自相关系数 rK 在 n较大， K较小时，计算公式为rK n KnK(Qt Q)(Qt K Q) n t 1n2(Qt Q )2t18-20 )据推导?2t?Q2(1 ?p,1r1 ?p,2r2?p,prp )8-21)在数学上序列，一般假定 t 为正态分布，故不需计算其偏态系数。但对于具有偏态的水文系列，一般把 t当作 P-型分布，因此还必须估计t的偏态系数 Cs t 。Cs t1(N P 3) ? tn( t?t )tP18-22 )t（t

10、=P+1，P+2，n）是根据估计出的以上 P+2个参数及观测序列 Qt ，利用公式（8-15 ） 反推得到的。平均值 ?t 是根据反推序列 t 利用矩法估计的。3常见 AR（1）及 AR（ 2）模型参数估计公式AR（ 1）模型形式为：Qt Q1,1（Qt 1 Q） t模型参数估计公式为：?1,1 r18-23)? t?Q 1 r12nCs t1 3 ( t ?)2s t (n 4) ?3t t 2 tAR( 2)模型形式为：?2,1 r1(1 r2 )/(1 r12)8-26 )?2,2 (r2 r12 )/(1 r12)?Q 1 ?2,1r1Cs1(n 5) ?3tn( t ?)3t38-3

11、1 )4模型阶数 P 的确定对于 AR（ P）序列，可以证明：它的自相关系数随滞时增大而减小，呈拖尾状，而偏相关系数 K,K 则呈截尾状， 在 K P时出现一个截止点， 即在 K P时， K,K 0，当 K P时， K,K 0。因此从理论上讲，可以通过计算不同的?K,K 进行模式识别。例如，当从样本序列估计 ?K,K在 K 3时具有明显的截尾现象， 那么可以推断该水文序列 P 3即适合于 AR（ 3）模型。但是由于实际水文样本容量较小， 故统计量 ?K,K 抽样误差较大， 即使是 AR（P） 序列，当 K P时， ?K,K 可能并不为零，这样就难于做直观判断，必须进行统计推断。统计推断方法是：

12、取显著水平 0.05 ，若 ?K,K 1.96 ，则可接受 ?K,K异于 0 的假 n设。 例如，某河流年径流偏相关系数如图 8-10 所示，该图表明， 只有 ?1,1超过 95%容许极限，即 ?1,1故该模型阶数应为1，即 AR（ 1）模型。(二) 单站月径流随机模型的建立对已有 n 年实测月径流资料的单站通常有两条途径建立其月径流随机模型。一是先建立 年径流模型，再通过建立解集模型把年径流分解成各月月径流。解集模型公式为：Y AQ B( 8-32 )公式中 Y - 各月月径流流量， Y (y1,y2, y12)T；Q- 年径流量；A - 模型参数， A (a1,a2, a12)T ，反映各

13、月月径流流量平均分配水平；B -12 × 12 的参数矩阵，反映各月之间的相关关系程度；- 模型残差项， ( 1, 2, 12)T 相互独立，可以是正态或偏态分布。以上参数 n 由年实测资料估算。本模型结构简单，概念清晰，但因参数多，故所需实测 资料较长。二是直接建立月径流随机模型，通常采用季节性一阶自回归模型，即假定可用 12 个一阶自回归模型来描述各月月径流流量及相关系数。各月月径流模型如下：(Qi,j 1 Qj 1)j 1 rj2,j 1 j(8-33 )jQi, j Qjrj,j 1j1公式中i - 年份， i 1,2,. ；j - 月份， j 1,2,.12 ；Qi,j-

14、第 i年第 j 月的月径流量；Qj ， j - 第 j 月的月径流均值和均方差， Q0 Q12 ， 012 ；rj,j 1- 第 j 月和第 j -1 月月径流之间的相关系数， r1,0表示第一月和上一年第十二月月径流 量相关系数；j - 第 j 月纯随机变量， 是模型残差想， 可以是标准正态分布或标准 P- 分布， 各月之间 jj 采用标准 P- 分布，则月径流量相互独立，且 j 与 Qi,j 1 相互独立。 以上各有关参数可由 n 年实测月径流资料用矩法估算。 若 j 采用标准正态分布，则月径流量也是正态分布。若为近似 P- 分布，生成P-分布 j 时，还需估算 j 的偏态系数CsCsQj

15、rj ,j 1CsQj 1(12 3/2 r j, j 1)8-34 )公式中 CsQj - 第 j月径流偏态系数， j 0时，表示的是第 12 月的偏态系数，该参数可用实 测 n 年第 j 月月径流序列估算。五、年月径流序列的生成(模拟)年月径流序列模拟涉及到纯随机变量t的生成及由 t 生成年月径流序列的问题。(一) 纯随机变量随机数的产生纯随机变量 t的分布可以是正态，也可以是偏态。它的生成一般先生成0，1 均匀分布随机数 u，再通过变换生成指定分布的随机数。1 均匀分布随机数的生成 生成方法有随机数表法、物理方法及数学方法。由于前两种方法存在严重缺陷，故常用 数学方法生成，其中应用最广的

16、是乘同余法。乘同余法生成随机数递推公式是：xn 1 MOD( xn,M) n 0,1,2,. (8 35)un 1 xn 1 /M (8 36)式(8 35)、式(836)中 x0为初值， 为乘子，为模，它们均为非负整数，而且<M。xn 1是 xn被 M整除后的余数，于是 xn 1<M，故 un 1即为0 ， 1 上的随机数。这种方法生成的随机数存在着循环周期， 因此， un 1不是真正意义上的随机数， 俗称“伪 随机数”。但由于 M往往取值很大，周期也很长，目前微机上周期可达109 以上，实用上完全能满足需要。正因如此，实际大都使用该法生成0，1 均匀分布随机数。不过使用前要对生

17、成的随机数作均匀性、独立性等检验。2 正态分布随机数的生成通常用 Box Muller 变换生成，即N12lnu1 cos(2 u2) (8 37)N12lnu1 cos(2 u2)中(8 37)中 u1， u2为0 ，1上的随机数， N1， N 2为相互独立标准化正态随机数。对于任意正态分布 N( , 2 ) ，随机数 N ' 生成公式 N'N ，N 为标准化正态分布随机数。3 P型分布随机数生利用舍选法生成 P型分布随机数，详见图 8 一 11。其中 ui (i 1,2, ,n 3)为为 0 ，1均匀分布随机数， z 为所需生成 P型分布随 机数。该分布 3 个参数为：均值

18、 Ex ，离势系数 Cv及偏态系数 Cs。在 Ex，Cv，Cs已知的 情况下，框图中 n、 p、 a、 的计算公式如下：2n INT (4 / Cs2)(838)p 4/Cs2 n(8 39)2/(ExCvCs) 或2/( Cs )( 8-40 )a Ex 2ExCv /Cs( 8-41 )(二)年月径流的生成AR(1) ，且参数下面仅介绍年径流生成的方法，月径流生成的方法类同。设所建立的模型为 已估计出。其模型为Qt Q1,1(Qt 1 Q) t(8-42 )式(8 42) 中 t 分布参数：均值为 0，均方差为 t 、，偏态系数为 Cs t 。下面分两种情 tt 况介绍年径流生成步骤：1.

19、 考虑年径流为正态分布这种情况下， t 为正态分布，即 Cs 可认为等于 0。t生成步骤：(1) 以 Q 或 Qt (t l 、2、 n)为 Q0 ；2(2) 生成一个符合 N(0, 2) 的正态随机数 1； t(3) 以Q0及 1代入式 (8 42)生成一个年径流 Q1；(4) 同步骤 (2) ，生成一个 2 ；(5) 以 Q1和 2代人式 (8 42) ，计算出 Q2；(6) 重复上述步骤，可得到一个很长 Qt 生成序列，如容量为 NN 50 的序列 Qt ；(7) 考虑到前 50 项可能受初值影响，应舍去，故剩下 NN 年生成的年径流系列。 2考虑年径流为偏态分布这种情况一般考虑 t为

20、P型分布， t 三个参数是：均值(均值为 0)，方差 2，偏态系 数 Cs t 。生成年径流 Qt 序列的方法与考虑年径流为正态分布时几乎一样，唯一不同的是上述第二步 t改用了 P型分布随机数生成。这样生成的Qt 序列可近似认为是 P型分布。六、模型及生成系列的检验模型检验是指所选用的模型残差 t 为独立随机变量是否成立， t 分布是否为假定分布。 生成序列检验是指所生成年月径流序列是否能反映实测径流的统计特性。1 残差独立性检验在模型及参数确定后，根据实测样本 Qt ，用式( 8-15 )可推出残差序列 t(t=P+1， n)，由 t序列可计算其各阶自相关系数 rK ，再对 rK 作独立性假

21、设检验。 当检验通过， 即 t 满 足独立性，说明建模时对 t 独立性假定是成立的，否则要分析产生的原因。若t序列存在密切自相关，应考虑使用其他模型。2 生成系列检验10一般要求生成序列与实测序列统计特征相近。如差异很大，要分析原因，确实是模型结 构问题，应考虑改变模型。七、实例红水河龙滩站年径流序列生成红水河龙滩水库为一多年调节水库，坝址处有自1946 1979年共 34 年的资料。为了分析工作需要，要求生成年径流系列。1径流组成分析 经过分析，末发现有趋势、突变、周期等确定成分，故可直接对实测序列建立平稳模型。 2模型选择和参数估计为了便于选择模型，已估算了年径流序列统计参数Q 、 Cv

22、、 Cs及自相关系数 r1、 r2、r15 ，见表 83 及图 812(a) 。从图 8 12(a) 可看出，自相关系数呈指数衰减趋势，故选用常用的 AR(P) 模型作为年径流模型。11为了确定模型阶数，计算了 ?K,K(K 1,2,.,15 )，互 1，2， 15)见图 812(b) 。从 图中可知， K3时， ?K ,K超过置信限范围， 而K>3时， ?K ,K均在置信限范围内， 说明 ?K,K 在 K> 3 久可以认为是 0，因此，应选定阶数为 3。模型参数利用前面介绍的矩法估算， 结果为： 3,1一 0021， 3,2 0143， 3,3 0333，t 0931。此外，利用

23、实测系列仅推出4、 5、 34 、计算出 t的Cs t，结果是 Cstt t t 接近于 0。3 年径流序列生成由于 Cst 很小，故把 t当作正态分布，即用正态分布生成 t ，利用三阶自回归模型生成 1000F 年径流系列。4 模型及生成系列的检验(1) 对于反推序列 4 、 5 、 34 ，求自相会系数 rK ，经检验可认为是独立的，因此，模 型的假定是成立的。(2) 生成系列的检验。对生成 1000F 年径流系列计算 Q 、 、 Cv 、 Cs 及自相关系数 r1 、r2 、，见表 83。通过对比发现两序列各项统计参数相对误差较小。 以上检验结果表明，所建模型是可接受的。总之，径流随机模

24、型是随着数学中随机过程理论和电子计算机技术在水文学中应用而逐渐 发展起来的，所有各种随机模型都是建立在水文现象的统计特性基础上，模型中参须基于实 测资料所提供信息加以估计，因此，只有深入了解水文现象的特性以及获得尽能多的可靠信 息，才能更有效地使随机模型应用于工程水文中。图 8 12 龙滩站自、偏相关函数示意图第七节 设计枯水流量分析计算(未吃透)枯水流量亦称最小流量， 是河川径流的一种特殊形态。 枯水流量往往制约着城市的发展规模、 灌溉面积、通航的容量和时间，同时，也是决定水电站保证出力的重要因素。 按设计时段的长短，枯水流量又可分为瞬时、日、旬、最小流量。其中又以日、旬、月最 小流量对水资源利用工程的规划设计关系最大。一、有实测水文资料时的枯水流量计算当设计代表站有长系列实测径流资料时，可按年最小选样原则选取一年中最小的时 段径流量，组成样本系列。枯水流量常采用不足概率q，即以小于和等于该径流的概率来表示，它和年最大选择的概率 P有 ql 一 P的关系。因