径流随机模拟(DOC)

径流随机模拟
由于实际水文资料往往比较短，难于满足实际水文程作随机模拟。

这种随机模拟的目的之一在于充分利用是用来延长资料长度。

当所建模型及参数准确时，这种年月径流随机模拟方法，对于多站及更深入的随机模拟
一、随机过程基本知识
(一)随机过程和时间序列的定义
在实际问题中，常涉及试验过程随某个参变量的变的流量、水位是随时间变化的随机变量，气温是随时间这种随机变量为随机函数，并称以时间为参数的随机函数为随机过程，记为}),({T t t ∈ξ，T 是t 变化的范围。

随机过程在一次试验或观测中所得结果，称为随机过程的一个实现。

若时间参变量T 是连续时刻的集合，则称这种随创机过程为连续参数随机过程，如水位过程、流量过程等。

若时间参变量T 是程为离散时刻的集合，则称这种随机过离散参数随机过程，也称为随机序列或时间序列。

如年、月径流程，年最大流量过程都是时间序列，也称水文时间序列。

(二)随机过程的数字特征
随机过程)(t ξ在任一固定时刻的状态是随机变量，因此可按与前述随机变量同样的方法定义随机过程的数学期望和方差。

定义如下
数学期望 )]([)(t E t ξμ= (8—11) 方差 })]()({[)(22t t E t μξσ-= (8—12)
为了规划随机变量两个不同时刻状态间关系的密切程度，可定义随机变量的自相关函数为
})
()()]
()()][()([{),(21221121t t t t t t E t t R σσμξμξ--= （8-13）
(三)随权过程基本分类
l 。

按统计性质的稳定性分类
按随机过程的统计性质是否随时间而变化，可分成平稳和非平稳过程。

若随机过程统计数字特征不随时间的平移而变化，则称为平稳过程，否则为非平稳过程。

2．按不同时刻状态间的关系分类
可分成独立过程和马尔柯夫过程。

若过程各状态相互独立，则称为独立随机过程。

在非独立随机过程中，最重要的一类是马尔柯夫过程，其特点是n t 时刻状态只与1-n t 时刻有关，而与1-n t 以前各时刻无关。

以上各个概念的严格定义，可参见有关文献。

二、径流随机模拟一船步骤
图8—7给出了单站年径流随机模拟的一般步骤：①时间序列组成分析；②模型的建立；②序列的生成；④模型及生成系列的检验。

三、水文时间序列的组成分析
水文序列t Q 一殷可按下式表示：
t t t t t S P C T Q +++= （8-14）
T t ，C t ，P t ，S t ——分别为趋势项、跳跃项、周期项和随机项。

当水文序列t Q 中不含T t ，C t ，P t 等确定性成分时，如t Q ＝S t ，即仅包括随机成分的序列。

对年径流序列而言，这种情况是出较常见的。

但月径流序列因存在明显的年周期，所以不是仅包括随机成分的序列。

趋势项指的是水文变量的统计参数(主要是均值)的长期系统性升降，如图8—8即存在明显趋势(增长)。

一般是由于气候因子或下垫面因子逐步改变而引起的缓慢变化。

对实测水文序列，可用假设检验或滑动平均的方法查明是否存在趋势。

若存在趋势，呈线性变化时，常用线性方程拟合，然后从序列中将趋势滤掉。

跳跃项是指水文序列急剧变化的一种形式，当水文序列从一种状态过渡到另一种状态时表现出来。

见图8-9。

跳跃是人为或自然原因造成的。

如建库后面积增大，蒸发量等损失增加，有可能出现跳跃，并反映在年径流序列的均值等参数。

突变可看作跳跃的一种特殊情况，如由于地震塌方，拦截江河，形成水库以后溃坝，这样引起流量的突变，随着临时水坝的冲毁，又恢复到原来状态。

跳跃是否存在于序列中，多用分割样本的方法检验。

若存在较显著的跳跃成分，应从序列中排除掉，使得剩余序列具有原始状态或一致条件。

周期项(含近似周期)是由于天体运动的周期性影响造成的，如地球公转、自转引起以年和日为周期的变化，以及太阳黑子活动引起的旱涝多年变化。

通常可用谐波分析的方法析出，再从序列中滤掉。

随机成分是由于不规则及随机振荡引起的，一般由相依成分和纯随机成分组成。

严格地讲，几乎所有水文变量应是非平稳过程。

不过除了人为影响及自然灾变外，水文环境的变化在数十年或几百年期间都相当小，因此，从实用观点，常把水文序列中的随机成分看作平稳过程。

图8—8序列趋势变化示意图
图8—9跳跃成分及包含跳跃成分的合成序列
四、单站年月径流随机模型的建立
对年月径流序列建立随机模型一般是对原始年月径流序列排除趋势跳跃等确定性成分后的随机成分而建立的模型。

设用于建模的年径流序列为1Q ，2Q ，…，n Q 。

径流序列随机模型是实测时间序列一种概化的数学式子。

目前常见的随机模型有：线性平稳模型、非线性平稳模型，马尔柯夫模型及非马尔柯夫模型等。

对年径流序列而言，常采用线性平稳模型或马尔柯夫模型。

建立随机模型的一般步骤为：①选择模型；②确定阶数；②估计模型参数
(一)单站年径流随机模型的建立
通常采用线性自回归模型，即马尔柯夫模型 1．线性自回归模型的一船形式
t p t p p t p t p t Q Q Q Q Q Q Q Q εϕϕϕ+-++-+-+=---)()()(,22,11, （8-15） t Q ——第t 年的年径流量，t=1，2，…，常称式（8-15）的t Q 为自回归系列。

Q ——t Q 序列的平均值；
1,p ϕ，…，p p ,ϕ——自回归系数或偏相关系数，反映t Q 在时间上相依性大小；
t ε ——模型残差项，纯随机成分，t ε 与1-t Q ，t Q ，…无关且是独立随机变量，其均值为
0，方差为2t εσ。

由于2t εσ与t Q 的方差2
Q σ有确定关系，因此，一般自回归模型中参数有：Q 、Q σ和1,p ϕ，
…，p p ,ϕ，共p 十2个参数。

该模型说明第t 年年径流量仅依赖于第t 一1年、第t 一2年、…、第t 一p 年的年径流
量和一个纯随机变量t ε。

若令Q Q y t t -=，则式(8—15)变为：
t p t p p t p t p t y y y y εϕϕϕ++++=---,22,11, (8—16) 式(8—16)是中心化变量表示的自回归模型。

2．模型参数的估计
∑==n
i i Q n Q 1
1
∑=--==n
i t Q Q Q Q n s 12)(11ˆσ ⎥⎥⎥⎥⎥
⎦⎤
⎢⎢⎢⎢⎢⎣
⎡⎥
⎥
⎥⎥
⎥⎦⎤
⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣
⎡------p p p p p p p p p r r r r r r r r r r r r 211
3
2
12111213,2,1,111ˆˆϕϕϕ
（8-19） 其中K 阶样本自相关系数K r 在n 较大，K 较小时，计算公式为
∑∑=-=+----=
n
t t
K
n t K t t K Q Q
Q Q Q Q
K
n n
r 1
2
1
)()
)(( （8-20）
据推导
)ˆˆˆ1(ˆˆ,22,11,2
2p p p p p Q t r r r ϕϕϕ
σσε ---= （8-21）
在数学上序列，一般假定t ε为正态分布，故不需计算其偏态系数。

但对于具有偏态的水文系列，一般把t ε当作P-Ш型分布，因此还必须估计t ε的偏态系数t Cs ε。

∑+=---=
n
P t t t
t t P N Cs 1
33)ˆ(ˆ)3(1
ε
ε
σ
εε （8-22） t ε（t=P+1，
P+2，…，n ）是根据估计出的以上P+2个参数及观测序列t Q ，利用公式（8-15）反推得到的。

平均值t ε
ˆ是根据反推序列t ε利用矩法估计的。

3．常见AR （1）及AR （2）模型参数估计公式
AR （1）模型形式为：
t t t Q Q Q Q εϕ+-+=-)(11,1 （8-23）
模型参数估计公式为：
11,1ˆr =ϕ
211ˆˆr Q t
-=σσ
ε ∑=--=
n
t t
s t t
n C 2
23)ˆ(ˆ)4(1
ε
ε
σ
εε （8-26） AR （2）模型形式为：
)1/()1(ˆ21211,2r r r --=ϕ
)1/()(ˆ212122,2r r r --=ϕ
22,211,2ˆˆ1ˆˆr r Q t
ϕϕσσ
ε--= ∑=--=
n
t t
s t t
n C 3
33)ˆ(ˆ)5(1
ε
ε
σ
εε （8-31） 4．模型阶数P 的确定
对于AR （P ）序列，可以证明：它的自相关系数随滞时增大而减小，呈拖尾状，而偏相关系数K K ,ϕ则呈截尾状，在P K =时出现一个截止点，即在P K ≤时，0,≠K K ϕ，
当P K >时，0,=K K ϕ。

因此从理论上讲，可以通过计算不同的K K ,ˆϕ
进行模式识别。

例如，当从样本序列估计K K ,ˆϕ
在3=K 时具有明显的截尾现象，那么可以推断该水文序列3=P 即适合于AR （3）模型。

但是由于实际水文样本容量较小，故统计量K K ,ˆϕ
抽样误差较大，即使是AR （P ）序列，当P K >时，K K ,ˆϕ
可能并不为零，这样就难于做直观判断，必须进行统计推断。

统计推断方法是：取显著水平05.0=α，若n
K K 96.1ˆ,>ϕ
，则可接受K K ,ˆϕ
异于0的假设。

例如，某河流年径流偏相关系数如图8-10所示，该图表明，只有1,1ˆϕ
超过95%容许极限，即n
96.1ˆ1,1>ϕ
，故该模型阶数应为1，即AR （1）模型。

（二） 单站月径流随机模型的建立
对已有n 年实测月径流资料的单站通常有两条途径建立其月径流随机模型。

一是先建立年径流模型，再通过建立解集模型把年径流分解成各月月径流。

解集模型公式为：
εB AQ Y += （8-32）
公式中 Y --各月月径流流量，T y y y Y ),,(1221 =；
Q --年径流量；
A --模型参数，T a a a A ),,(1221 =，反映各月月径流流量平均分配水平；
B --12×12的参数矩阵，反映各月之间的相关关系程度；
ε --模型残差项，T ),,(1221εεεε =相互独立，可以是正态或偏态分布。

以上参数n 由年实测资料估算。

本模型结构简单，概念清晰，但因参数多，故所需实测资料较长。

二是直接建立月径流随机模型，通常采用季节性一阶自回归模型，即假定可用12个一阶自回归模型来描述各月月径流流量及相关系数。

各月月径流模型如下：
j j j j j j i j j j j
j j i r Q Q r Q Q εσσσ21,11,1,1
,1)(------+-+
= （8-33） 公式中
i --年份，,...2,1=i ；
j --月份，12,...2,1=j ；
j i Q ,--第i 年第j 月的月径流量；
j Q ，j σ--第j 月的月径流均值和均方差，120Q Q =，120σσ=；
1,-j j r --第j 月和第j -1月月径流之间的相关系数，0,1r 表示第一月和上一年第十二月月径流
量相关系数；
j ε--第j 月纯随机变量，
是模型残差想，可以是标准正态分布或标准P-Ⅲ分布，各月之间j ε
相互独立，且j ε与1,-j i Q 相互独立。

以上各有关参数可由n 年实测月径流资料用矩法估算。

若j ε采用标准正态分布，则月径流量也是正态分布。

若j ε采用标准P-Ⅲ分布，则月径流量为近似P-Ⅲ分布，生成P-Ⅲ分布j ε时，还需估算j ε的偏态系数
2
/321,31,)1(1
----=
-j j sQ j j sQ s r C r C C j j j ε （8-34）
公式中j sQ C --第j 月径流偏态系数，0=j 时，表示的是第12月的偏态系数，该参数可用实测n 年第j 月月径流序列估算。

五、年月径流序列的生成（模拟）
年月径流序列模拟涉及到纯随机变量t ε的生成及由t ε生成年月径流序列的问题。

（一） 纯随机变量随机数的产生
纯随机变量t ε的分布可以是正态，也可以是偏态。

它的生成一般先生成[0，1]均匀分布随机数u ，再通过变换生成指定分布的随机数。

1．均匀分布随机数的生成
生成方法有随机数表法、物理方法及数学方法。

由于前两种方法存在严重缺陷，故常用数学方法生成，其中应用最广的是乘同余法。

乘同余法生成随机数递推公式是：
),(1M x MOD x n n λ=+ ,...2,1,0=n (8—35) M x u n n /11++= (8—36)
式(8—35)、式(8—36)中0x 为初值，λ为乘子，为模，它们均为非负整数，而且λ＜M 。

1+n x 是n x λ被M 整除后的余数，于是1+n x ＜M ，故1+n u 即为[0，1]上的随机数。

这种方法生成的随机数存在着循环周期，因此，1+n u 不是真正意义上的随机数，俗称“伪随机数”。

但由于M 往往取值很大，周期也很长，目前微机上周期可达109
以上，实用上完全能满足需要。

正因如此，实际大都使用该法生成[0，1]均匀分布随机数。

不过使用前要对生成的随机数作均匀性、独立性等检验。

2．正态分布随机数的生成
通常用Box —Muller 变换生成，即
)2cos(ln 2211u u N π-= (8—37)
)2cos(ln 2211u u N π-=
中(8—37)中1u ，2u 为[0，1]上的随机数，1N ，2N 为相互独立标准化正态随机数。

对于任意正态分布),(2σμN ，随机数'
N 生成公式N N σμ+='，N 为标准化正态分布随机数。

3．P —Ⅲ型分布随机数生
利用舍选法生成P —Ⅲ型分布随机数，详见图8一11。

其中i u （3,,2,1+=n i ）为为[0，1]均匀分布随机数，z 为所需生成P —Ⅲ型分布随机数。

该分布3个参数为：均值x E ，离势系数v C 及偏态系数s C 。

在x E ，v C ，s C 已知的情况下，框图中n 、p 、a 、λ的计算公式如下：
)/4(2s C INT n = (8—38)
n C p s -=2/4 (8—39)
)/(2s v x C C E =λ或)/(2s C σλ= （8-40）
s v x x C C E E a /2-= （8-41）
(二)年月径流的生成
下面仅介绍年径流生成的方法，月径流生成的方法类同。

设所建立的模型为AR(1)，且参数已估计出。

其模型为
t t t Q Q Q Q εϕ+-+=-)(11,1 （8-42）
式(8—42)中t ε分布参数：均值为0，均方差为t εσ、，偏态系数为t s C ε。

下面分两种情况介绍年径流生成步骤：
1. 考虑年径流为正态分布
这种情况下，t ε为正态分布，即t s C ε可认为等于0。

生成步骤：
(1)以Q 或t Q (t ＝l 、2、…、n)为0Q ； (2)生成一个符合),0(2
t N εσ的正态随机数1ε； (3)以0Q 及1ε代入式(8—42)生成一个年径流Q 1； (4)同步骤(2)，生成一个2ε；
(5)以Q 1和2ε代人式(8—42)，计算出Q 2；
(6)重复上述步骤，可得到一个很长t Q 生成序列，如容量为50+NN 的序列t Q ； (7)考虑到前50项可能受初值影响，应舍去，故剩下NN 年生成的年径流系列。

2．考虑年径流为偏态分布
这种情况一般考虑t ε为P —Ⅲ型分布，t ε三个参数是：均值（均值为0），方差2
t εσ，偏态系数t s C ε。

生成年径流t Q 序列的方法与考虑年径流为正态分布时几乎一样，唯一不同的是上述第二步t ε改用了P —Ⅲ型分布随机数生成。

这样生成的t Q 序列可近似认为是P —Ⅲ型分布。

六、模型及生成系列的检验
模型检验是指所选用的模型残差t ε为独立随机变量是否成立，t ε分布是否为假定分布。

生成序列检验是指所生成年月径流序列是否能反映实测径流的统计特性。

1．残差独立性检验
在模型及参数确定后，根据实测样本t Q ，用式（8-15）可推出残差序列t ε（t=P+1，…，n ），由t ε序列可计算其各阶自相关系数K r ，再对K r 作独立性假设检验。

当检验通过，即t ε满足独立性，说明建模时对t ε独立性假定是成立的，否则要分析产生的原因。

若t ε序列存在密切自相关，应考虑使用其他模型。

2．生成系列检验
一般要求生成序列与实测序列统计特征相近。

如差异很大，要分析原因，确实是模型结构问题，应考虑改变模型。

七、实例——红水河龙滩站年径流序列生成
红水河龙滩水库为一多年调节水库，坝址处有自1946—1979年共34年的资料。

为了分析工作需要，要求生成年径流系列。

1．径流组成分析
经过分析，末发现有趋势、突变、周期等确定成分，故可直接对实测序列建立平稳模型。

2．模型选择和参数估计‘ 为了便于选择模型，已估算了年径流序列统计参数Q 、v C 、s C 及自相关系数1r 、2r 、…、15r ，见表8—3及图8—12(a)。

从图8—12(a)可看出，自相关系数呈指数衰减趋势，故选用常用的AR(P)模型作为年径流模型。

为了确定模型阶数，计算了K K ,ˆϕ（15,...,2,1=K ）
，互＝1，2，…，15)见图8—12(b)。

从图中可知，K ＝3时，K K ,ˆϕ
超过置信限范围，而K ＞3时，K K ,ˆϕ均在置信限范围内，说明K K ,ˆϕ在K ＞3久可以认为是0，因此，应选定阶数为3。

模型参数利用前面介绍的矩法估算，结果为：1,3ϕ＝一0．021，2,3ϕ＝0．143，3,3ϕ＝0．333，t εσ＝0．931。

此外，利用实测系列仅推出4ε、5ε、…34ε、计算出t ε的t s C ε，结果是t s C ε接近于0。

3．年径流序列生成
由于t s C ε很小，故把t ε当作正态分布，即用正态分布生成t ε，利用三阶自回归模型生成1000F 年径流系列。

4．模型及生成系列的检验
(1)对于反推序列4ε、5ε、…34ε，求自相会系数K r ，经检验可认为是独立的，因此，模型的假定是成立的。

(2)生成系列的检验。

对生成1000F 年径流系列计算Q 、σ、v C 、s C 及自相关系数1r 、2r 、…，见表8—3。

通过对比发现两序列各项统计参数相对误差较小。

以上检验结果表明，所建模型是可接受的。

总之，径流随机模型是随着数学中随机过程理论和电子计算机技术在水文学中应用而逐渐发展起来的，所有各种随机模型都是建立在水文现象的统计特性基础上，模型中参须基于实测资料所提供信息加以估计，因此，只有深入了解水文现象的特性以及获得尽能多的可靠信息，才能更有效地使随机模型应用于工程水文中。

图8—12龙滩站自、偏相关函数示意图
第七节 设计枯水流量分析计算（未吃透）
枯水流量亦称最小流量，是河川径流的一种特殊形态。

枯水流量往往制约着城市的发展规模、灌溉面积、通航的容量和时间，同时，也是决定水电站保证出力的重要因素。

按设计时段的长短，枯水流量又可分为瞬时、日、旬、…最小流量。

其中又以日、旬、月最小流量对水资源利用工程的规划设计关系最大。

一、有实测水文资料时的枯水流量计算
当设计代表站有长系列实测径流资料时，可按年最小选样原则选取一年中最小的时段径流量，组成样本系列。

枯水流量常采用不足概率q ，即以小于和等于该径流的概率来表示，它和年最大选择的概率P 有q ＝l 一P 的关系。

因此在系列排对时按右小到大排列。

除此之外，年枯水流量频率曲线的绘制与时段径流频率曲线的绘制基本相同，也常采用P-Ⅲ型频率曲线适线。

图
8-13为某水文站不同天数的枯水流量频率曲线的事例。

年枯水流量频率由线，在某些河流上，特别是在干旱半干旱地区的中小河流上，还会出现时段径流量为零的现象，可参阅含零系列的频率分析方法。

此处只介绍一种简易的实用方法。

设系列的全部项数为n ，其中非零项数为k ，零值项数为n —k 。

首先把k 项非零资料视作一个独立系列，按一般方法求出其频率曲线。

然后通过下列转换。

即可求得全部系列的频率曲线。

其转换关系为：
非设P n
k P (8—43) 式中设P ——全系列的设计频率；
非P ——非零系列的相应频率。

在枯水流量频率曲线上，往往会出现在两端接近P ＝20％和P ＝90％处曲线转折现象。

在P=20%以下的部分是河网及潜水逐渐枯竭，径流主要靠深层地下水补给。

在P=90%，可能是某些年份有地表水补给，枯水流量偏大所致。

二、短缺水文资料时的枯水流量估算
当设计断面短缺径流资料时，设计枯水量主要借助于参证站延长系列或成果移置，与本章第三节所述方法基本相同。

但枯水流较之固定时段的径流，其时程变化更为稳定。

因此，在与参证站建立径流相关时，效果会好一些。

或者说，条件可以适当放宽。

例如，当设计站只有少数几年资料，与参证站的相似性较好时，也可建立较好的枯水流量相关关系。

在这种情况下，甚至可以不进行设计站的径流系列延长和频率分析，而直接移用参证站的频率分析成果，经上述相关关系，转化为本站的相应频率的设计枯水流量。

在设计站完全没有径流资料的情况下还可以临时进行资料的补充收集工作，以应需要。

如果能施测一个枯水季的流量过程，对于建立30天以下时段的枯水流量关系，有很大用处；如只研究日最小流量，那么在枯；水期只施测几次流量(如10次流量)，就可与参证站径流建立相关关系。