拒马河上游年径流序列特性分析及其随机模拟
2015-11-24王玲
王 玲
(河北省保定水文水资源勘测局,河北 保定071000)
1 基本情况
拒马河是大清河的干流,主要为泉水,从涞源县城发源,流经易县紫荆关,涞水野三坡,北京房山十渡,在房山张坊分为南拒马河和北拒马河。
北拒马河流经北京南尚乐乡,于二合庄村东出市境,入河北涿州市境内,至东茨村以下称白沟河,在白沟村与南拒马河汇合入大清河,汇入白洋淀[1]。
长期以来拒马河上游未建大型水库,河畔两侧多山峦,遇暴雨天气,极易发生洪水,且来势凶猛。紫荆关水文站位于拒马河上游,对于拒马河的防洪有着重要作用。2012年7月21~22日拒马河发生了1963年以来最大的一场暴雨洪水,紫荆关水文站洪峰流量达2580m3/s,因此对紫荆关站径流系列的随机模拟具有重要意义。
紫荆关水文站位于河北省易县境内,流域分布在涞源县、易县境内,地势西高东低,山峦起伏为太行山脉,一般山峰海拔在2100m以下,紫荆关水文站流域面积1760km2,主河道长81.5km,河道纵坡5.5%,流域平均宽度25.4km,最大洪峰流量4490m3/s(1963年)[2]。
2 资料审查
选用紫荆关站径流序列资料 (1965~2013年)均系水文部门实测资料。多年来,位于拒马河上游的紫荆关站集水面积和气候条件并没有发生大的变化,也没有大的水利工程建设,从而能够满足资料一致性的要求。由于紫荆关水文站40多年的年径流序列有很好的稳定性,而且表现为丰枯水年组的周期性变化,可以表明所选径流资料有很好的代表性。
3 紫荆关站年径流序列自相关性分析及模型识别
根据模型识别的基本原理,径流序列究竟符合哪一类随机模拟模型,取决于该径流序列的自相关性和其偏态特性,为此首先分析紫荆关站年径流序列的自相关性[3]。
年径流序列Xt的自相关系数表示为:
式中 rk为年径流序列Xt的自相关系数;n为序列长度,当n>50时,可以取m<n/4,常取m在n/10左右;当n<50时,取m在n/4左右的数值,有的取m<n-10;k为年径流序列的滞时;为年径流序列的均值。
本文所选紫荆关站年径流序列年数为49a,取m=12,计算自相关系数并绘制自相关图,并加绘了独立序列自相关置信水平为95%的容许限,如表1和图1。
表1 紫荆关站49a径流序列自相关系数及其容许限
图1 紫荆关站年径流自相关图
表1和图1表明:紫荆关站年径流量序列的一阶和二阶自相关系数均显著异于独立序列,因此为一组相依序列。对正相依序列,可以选用不少类型的模型来描述其统计特征。考虑以下几点,选用自回归模型AR(p)来模拟紫荆关站的径流序列。
(1)AR(p)模型表征径流序列的统计特性有一定的物理基础;
(2)AR(p)模型参数的估计可以用简单的矩法,而且精度高;
(3)AR(p)模型形式简单,数学处理方法简便,为大家熟悉。
4 紫荆关站年径流序列偏态特性分析
选定AR(p)模型后,如何确定p,需分析其序列的偏态特性,主要看是否在k>p后,φp,p偏态系数在95%的容许限范围内,当假设成立,可以认为能够用AR(k)模拟该径流序列[4]。
对于自回归AR(p)模型有
式中 φp,1为p阶模型的第一系数;φp,p为p阶模型的第p个系数。
式(6)中,偏相关系数φp,p是反映消除(p-1)阶自相关系数影响后所剩余的自相关程度,计算方法上,偏相关函数和自相关函数有一定联系,利用自相关函数可以计算偏相关函数。可令φk,j表示自回归AR(k)模型中的第j个系数,这样φk,k为最后一个系数。则:
可求的φk,j(j=1,2,…,k),对于自回归AR(1)模型
对于高阶自回归模型,可应用以下述递推算法:
计算结果如表2,偏相关函数图如图2。
表2 紫荆关站49a径流序列偏相关系数及其容许值
续表2
图2 紫荆关站年径流偏相关图
图2表明:当k≥2时,φk,k落于容许限内,故可推断p=1。换言之,据偏相关系数的统计分析,AR(1)模型可以用来描述年径流量序列的统计变化。但是,在序列长度n相当大且序列为正态分布的情况下才是完全正确的,就紫荆关年径流量序列而言,其长度仅有49a,而其二阶自相关系数较大,这暗示该序列可能为AR(2)序列。有以上分析,可以认为紫荆关年径流量序列可能是AR(1)序列,也可能是AR(2)序列。下面对这两种可能的模型均作参数估计,以便进一步检验。
5 参数估计及年径流序列模拟
5.1 参数估计
通过对参数的平稳性分析,最终求得AR(1)模型和AR(2)模型。
AR(1)模型为:
式中 φt为服从均值0、方差1的偏态系数;Csφ的P-Ⅲ型分布,用公式(11),(12)可求得。
AR(2)模型为:
5.2 模型形式的进一步判断
估计出参数后,便可利用AIC准则进一步识别紫荆关站年径流序列是AR(1)还是AR(2)模型。由式(11)计算AIC值,从表3可以看出,由AIC准则来判断,AR(1)模型优于AR(2)模型。
表3 紫荆关站49a径流序列AIC值
5.3 模型检验
主要检验εt是否独立,对AR (1)模型,令εt=0.735φt,则有:
根据实测样本序列可计算得εt,由此计算出r1(ε),r2(ε),…,r12(ε),最后计算统计量Q=8.91。根据自由度和显著水平,得Q<χ2α,故为独立的假设可以接受。
φt为自回归模型的随机变量部分,用式(11)、式(12)计算:
式中 Csφ为φ分布的偏态系数;Csx为X分布的偏态系数;εt为随机数。
模拟随机数εt时,本文采用变换法,即对均匀随机数做下列变换:
则ξ1,ξ2为相互独立的标准正态分布N(0,1)变量,依次模拟出ξt。
6 模型检验和实用性分析
重点分析AR(1)能否保持实测序列的主要统计特征[5]。分别根据AR(1)模型的模拟径流序列,然后按长序列和短序列法计算各种参数,并和实测序列的相应参数做对比,结果如表4。
表4 模拟序列的实用性检验
由表4可以看出,长序列的模拟结果优于短序列的结果,长序列的AR(1)的基本参数的误差在10%以内,均能很好地保持,短序列模拟的各参数均偏小,这种偏小是由于计算Cv,Cs,r等公式在样本容量较小时皆为负偏而造成。因此使用一阶自回归模拟AR(1)能够较好地模拟拒马河上游的年径流序列,基本满足实际工程应用需要。
7 结语
通过对拒马河水系上游紫荆关站年径流变化趋势分析和随机模拟,可以清晰掌握拒马河的年径流的变化规律。通过进一步的水文模拟,可以更清楚地看出拒马河紫荆关站年径流量的变化范围,从而能够在一定程度上服务于紫荆关站的防洪和水利工程建设。
[1]郭庆宏,朱虹.河北省拒马河山区降水特性分析[J].河北工程技术高等专科学校学报,2012(2):1-4.
[2]李绍飞,余萍,等.紫荆关流域洪水径流过程变化及影响因素分析[J].武汉大学学报(工学版),2012,45(2):167-170.
[3]丁晶,邓育仁.随机水文学[M].成都:成都科技大学出版社,1998.
[4]Bras R A.Random Functions and Hydrology [M].America:Addison Wesley Publication Co.,1985.
[5]金光炎.水文统计原理与方法[M].北京:中国工业出版社,1954.