门限自回归模型及其在水文随机模拟中的应用.王文圣

门限自回归模型及其在水文随机模拟中的应用*王文圣,　袁　鹏,　丁　晶,　邓育仁(四川大学水电学院,四川成都　610065)摘　要:为了客观描述日流量变化的非线性特性,将一种非线性时序模型——门限自回归模型引入日流量随机模拟。

根据我国金沙江流域屏山站观测资料建立了日流量随机模拟的门限自回归模型。

实用性检验结果表明,该模型用于模拟日流量过程是可行的,成果实用。

这种尝试为日流量随机模拟提供了一种考虑日流量非线性变化特性的新模型。

关键词:门限自回归模型;日流量随机模拟;实用性检验中图分类号:P33;P333.6文献标识码:　B文章编号:1001-2184(2001)增-0047-041　引　言日流量随机模拟利用日流量涨落的统计特性,具体说是利用日流量在时序上的统计关系。

这种统计关系非常复杂,为简化处理常常以线性来表征前后日流量的关系。

在一般情况下,这种简化尚能反映日流量时序变化的主要特性。

所以在日流量随机模拟时,当前广泛使用线性时序模型。

但是日流量在时序上的前后流量关系是非线性的。

例如,对大流域一次洪水的日流量过程涨水段的下部、中部和上部有着明显不同的涨率,前后流量关系显然不是线性的;同样在落水段的下部、中部和上部有着明显不同的退水率,前后流量关系也不是线性的。

因此,为更全面地反映日流量时序变化的特性,最好考虑日流量在时序变化上的非线性特性。

近来,非线性时序的分析获得了迅速的发展,并且相继出现了一系列非线性时序模型,比如门限自回归模型,双线性模型,指数自回归模型,状态依赖模型等。

对双线性模型曾初步研究了在洪水模拟中应用的可能性[1]。

门限自回归模型最近尝试应用于水文预报并获得较好的效果[2]。

鉴于门限自回归模型在表征非线性特性上具有其独到之处,笔者将之引入日流量随机模拟并以某站日流量资料为基础,全面探讨了这种模型在日流量模拟中的可行性,模拟效果和优缺点等。

2　门限自回归模型的形式和基本特性2.1　模型形式门限自回归模型由汤家豪1978年提出[3],用来解决一类非线性问题。

其思路是:对研究对象按照不同区间建立若干个线性时序模型;然后将这些线性时序模型组合起来描述该对象非线性时序变化特性。

对于时间序列{Z t},门限自回归模型的一般形　收稿日期:2000-08-14*　基金项目:国家自然科学基金(49871018);高速水力学国家重点实验室开放基金资助项目(编号2008)式为:Z tU(1)0+∑p1i=1U(1)i Z t-i+E(1)t Z t-d F r1U(2)0+∑p2i=1U(2)i Z t-i+E(2)t r1<Z t-d F r2……U(L)0+∑pLi=1U(L)i Z t-i+E(L)t r L-1F Z t-d(1)式中　r1,r2……,r L-1为门限值;d为门限迟时;L为门限区间数;U(j)0,U(j)1,……U(j)p j为第j区间自回归系数;P\-j为第j区间模型阶数(j=1,2,……,L); E(j)t是独立随机变量。

从(1)式可看出,门限自回归模型是将{Z t}按{Z t-d}值的大小采用对应区间的线性自回归模型来描述整个系统。

2.2　模型特性门限自回归模型是具有广泛意义的一种非线性时序模型,实质上是以分区间线性自回归模型来描述研究对象在整个区间的非线性变化特性。

下面以简单的一阶非线性自回归模型(NLAR(1))说明。

NLAR(1)模型的形式为:Z t=f(Z t-1)+a t(2)式中　f(Z t-1)为区间[a,b]上的连续非线性函数。

{a t}为白噪声序列。

若将{Z t-1}(t=2,3……,N)的取值区间划分成L个小区间,见图1。

当某一Z落在j区间(r j-1,r j)时,f(Z)由fd(Z)任意逼近,即f(Z)≈fd(Z)=f(r j-1)+A j Z(j=1,2……,L),其中f(r j-1)和A j的取值由第j区间(r j-1,r j)决定。

因此: Z t=f(r j-1)+A(j)j Z t-1+a(j)t　(j=1,2,……,L)(3)取f(r j-1)=U(j)0,A(j)j=U(j)1,则NLAR(1)模型可由下式近似表达:Z t=U(j)0+U(j)1Z t-1+a(j)tr j-1<Z t-1F r j(j=1,2,……,L)(4)式(4)即所谓的一阶门限自回归模型。

类似地,对于第20卷增刊2001年7月四　川　水　力　发　电Sichuan　Water　Pow erVol.20,S upplemen tJul.,　2001一般高阶非线性自回归模型可用高阶门限自回归模型逼近。

图1　非线性自回归示意图3　门限自回归模型建立门限自回归模型是分区间的自回归模型(AR 模型),因此在建模中,只需沿用一般AR 模型的参数估计方法和模型检验准则。

目前多采用最小二乘法和A IC 准则。

与AR 建模有所不同的是,建立门限自回归模型的关键不在于参数估计和模型检验,而在于确定门限参数(门限区间数L ,门限值r 1,r 2,……,r L -1和门限迟时d )。

严格讲,这是一个对L ,r 1,r 2,……,r L -1和d 的多维寻优问题[5]。

通常以A IC 准则作为模型定阶(确定各门限区间AR 模型阶数P 1,P 2,……P L )和确定L ,r 1,r 2,……,r L -1,d 参数的目标函数,即:A IC =f (P 1,P 2,……,P L ,d ,r 1,r 2,……r L -1)(5)A IC 值最小时对应的模型参数、门限参数为所求。

可见,上述寻优工作量相当大,它是多维空间中超曲面的极小点的寻找。

为减少工作量,近来出现了一些改进办法,如D .D .C 方法,局部区间搜索法等[5],其目的是缩小寻优区间,但这些方法在实际计算时工作量仍很大且较繁琐,有时甚至无法实现。

鉴于此,笔者在建模过程中提出了一些改进。

为了减少寻优工作量,门限区间数L 和门限迟时d 的选取要结合研究对象的物理成因和一般变化规律,取值范围可得到有效的缩小。

比如结合研究对象L 取2,3,4等,d 取1,2,3等。

门限值可按研究对象在变幅范围内取值的经验频率分布来确定。

一般要求在各区间内样本的资料数据大致相当,以便可靠地估计模型参数。

为此,先均匀划分经验频率,然后据经验频率在频率分布曲线上找出相应的门限值。

所以,关键是确定门限值所对应的频率。

如L =2时,门限值对应的频率可取0.5;L =3,门限值对应的频率可取0.33和0.67。

上述结合实际问题初选门限参数的考虑,使寻优工作量大大降低了。

综上所述,门限自回归模型建模过程如下:首先固定一组d ,L ,r j (j =1,2,……,L -1),由最小二乘法和A I C 准则分别确定各区间AR 模型,从而得到相应的门限自回归模型;然后再分别改变d ,L ,r j (j =1,2,……,L -1),同样,建立分区间AR 模型以得到门限自回归模型。

比较各种d ,L ,r j 情况下的A I C 值,其中AI C 值最小对应的模型即为所求。

4　模型应用以我国金沙江流域屏山站日流量资料为例,将门限自回归模型尝试于日流量过程的随机模拟。

4.1　基本资料预处理选用屏山站53年(1940～1992年)日流量序列作为建模的基本资料。

为方便,每年2月份取28d ,故每年共365d 。

以X i ,j 表示第i 年j 截口日流量,并假定同一截口所有元素来自同一总体,各截口具有相同的分布型式。

由基本资料用概率权重矩法[4]计算截口统计参数(X -j ,S j ,Cv j ,Cs j ,j =1,2,……365)。

结果表明,它们随截口有明显的变化,即该站日流量过程是非平稳的,其分布是非正态的。

为便于建模需要对原始日流量作标准化和正态化预处理。

前者消除日流量过程均值和标准差的非平稳性,后者将偏态分布转换为正态分布。

标准化由(6)式进行:Y i ,j =X i ,j -X -jS j(6)式中　Y i ,j 为标准化日流量序列;X -j ,S j 分别为第j截口均值和标准差;X i ,j 同前。

将偏态序列转换为正态序列的方法较多,本文用H-W 逆变换法[6],式为:Z i ,j =6Cs jCs j2Y i ,j +113-1+Cs j 6(7)式中　Cs j 为截口j 的偏态系数;Y i ,j 同前;Z i ,j 为标准化正态序列。

经(6)和(7)式预处理后的{Z i ,j }为平稳正态序列,用它即可进行模型参数估计和建立门限自回归模型。

4.2　模型建立及其模拟4.2.1　模型建立对日流量时序过程,前后流量统计关系的差异,一般说来,在过程上、中和下三部分之间表现较为明显。

因此,门限区间数不宜定得太多,以3左右为宜,即L =2,3,4三个优选备用值。

门限区间数确定后,相应的门限值可由均匀划分的经验频率估计。

至于门限迟时由于日流量前后关系密切,取1或2即可。

根据上面选定的各种门限参数,互相组合便可以用前述的A I C准则进行优化最后确定所求的门限自回归模型。

由该站预处理日流量序列{Z i,j}通过优化得到d=1,L=3,r1=-0.520,r2=0.533和相应的门限自回归模型:Z t=-0.08+0.9674Z t-1-0.0501Z t-2+0.013Z t-3+0.0348Z t-4-0.0537Z t-5+E(1)tZ t-1F-0.520,R(1)E=0.3001.3678Z t-1-0.3524Z t-2-0.0714Z t-3+0.0365Z t-4+E(2)t-0.520<Z t-1F0.533,R(2)E=0.2000.01+1.3371Z t-1-0.4428Z t-2+0.0799Z t-3-0.0653Z t-4+0.0589Z t-5+E(3)tZ t-1>0.533,R(3)E=0.208(8)式中　E(1)t,E(2)t,E(3)t是均值为0,方差为0.3002, 0.2002,0.2082的独立正态序列。

4.2.2　模型模拟程序模型建立后,就可以进行全年日流量序列随机模拟。

模拟程序如下:¹由模型式(8)模拟标准正态日流量Z i,j(i=1,2,……,N;j=1,2,……,365);º由式(9)进行H-W变换,获得均值为0,方差1的P-Ⅲ型分布Y i,j;»由逆标准化(X i,j=X-j+S j Y i,j)求得日流量过程X i,j。

实际计算按¹～»步一年一年地模拟。

根据模拟序列估计有关参数和统计有关特征。

H-W变换式为:Y i,j=2Cs j(1+Z i,j Cs j6-Cs2j36)3-2Cs j(9)4.3　模型检验对建立的门限自回归模型是否满足建模要求,是否反映日流量过程的各种统计特性和客观变化规律,必须进行检验。