基于核估计的多变量非参数随机模型初步研究 王文圣1，丁晶1 (1.四川大学 水利水电学院，四川 成都 610065)
摘要：本文基于核估计理论构造了多变量非参数模型。该模型是数据驱动的、不需识别和假定序列相依形式和概率分布形式的一类随机模型，克服了多变量参数模型的不足。实例统计试验表明，建议的多变量非参数模型是有成效的，为随机水文学发展提供了一些新思路。 关键词：核估计；多变量非参数模型；随机模拟；实用性检验 中图分类号：P333.9 文献标识码：A
流域水资源的开发利用，不仅需要单站水文信息，而且需要流域内各站的水文信息。进行多站水文序列模拟的一个重要手段就是建立多站(变量)随机模型。目前，多变量随机模型[1]比较成熟的有自回归模型和解集模型。这两类模型的共同点是用有限个参数的线性函数关系描述水文现象。因此简便实用，能表征水文序列的统计特性和一般变化规律，但缺点也明显：①水文序列是一时间不可逆过程，而参数模型描述的是可逆过程，因此大多数参数模型难以反映其涨落不对称性；②水文现象受流域下垫面、人类活动、气候等多因素影响而变化错综，是一个高度复杂的非线性系统，而多数参数模型仅能表征变量及变量之间的线性相依结构，忽略了占据重要位置的非线性性；③水文变量概率密度函数复杂且未知，某一指定概率分布与真实分布存在着差异。如图1、2所示，正态分布、P-Ⅲ型分布都与直方图相差甚远，但χ2检验并不拒绝P-Ⅲ型分布和正态分布；而核估计和k最近邻估计与直方图比较接近。即概率分布具有不确定性；④模型参数由于抽样误差和估计方法不同具有不确定性。
为克服参数模型之不足，文献[2]提出了单变量非参数模型，径流模拟表明是满意的。在此基础上，本文基于核估计理论构造了多变量非参数模型。该模型避开了序列相依形式和模型结构的假设，不涉及模型参数估计，能反映各种复杂关系，较参数模型优越。以中国金沙江流域屏山站和宜宾—屏山区间两站日流量过程随机模拟为例，对建议模型进行了应用研究。
1 核估计理论[3] 1.1 多维核估计定义 设X为d维随机变量，X1，X2，……Xn为X的一样本。X的概率密度函数f(X)的核估计定义如下：
(1) 式中：X=(x1,x2,…,xd)T,Xi=(xi1,xi2,…，xid)T(i=1,2,…,n)；K(·)为核函数，是一给定概率密度函数；h为带宽系数；n为样本容量；S是X的d×d维对称样本协方差矩阵。
1.2 核估计精度评价 核估计既同样本有关，又与K(·)和h的选取有关。在给定样本后，核估计的精度取决于K(·)及h的选取是否适当。常采用积分均方误差准
则 进行度量。MISE由偏差和方差组成。当K(·)固定时，若h选得过大，偏差较大，但降低了方差，故(X)对f(X)有较大的平滑，使得f(X)的某些特征被掩盖起来；若h选得过小，偏差减小了，但增大了方差，则(X)有较大的波动。显然同时减少偏差和方差是不可能的。
图1 屏山站1月31日流量概率密度估计 图2 宜-屏区间7月径流量概率密度估计(图例见图1) 1.3 K(·)和h的确定 依潘涅契科夫[4]和Scott[5]通过统计试验发现，当给定带宽系数，不同核函数对MISE的影响是很小的。实际工作中，选择满足一定条件的核函数即可[2]。本文采用标准高斯函数。 h随n增大而减小。h的确定还要考虑数据的密集程度，在数据密集区，h选小一点；在数据稀疏区，h选大一点。h的具体计算方法很多。这里使用LSCV法。LSCV是基于积分平方误差(Integrated Square Error(ISE))最小准则的一种计算方法。对多维随机变量X，ISE为：
(2) 式(2)中最后一项与h无关。LSCV就是取式(2)中前两项进行最小化，即 (3)
式中：当LSCV(h)最小时，h即为所求。 2 基于核估计的多变量非参数随机模型 2.1 模型建立 所谓多变量，可以是同一测站的几个水文变量，也可以是不同测站上的一种或几种变量。设｛Xt｝n(t可以是年，也可以是季、月、日)为多变量相依时间序列，Xt依赖于前P个值Xt-1,Xt-2,…，Xt-p,取Vt=(Xt-1,Xt-2,…，Xt-p)T，Xt的条件概率密度函数为：
(4) (4)式中： ,m为变量个数， 为j变量t时刻水文变量值；f(Xt,Vt)为m(P+1)维联合密度函数，fV(Vt)为mP维边缘密度函数。 由多维核估计知 (5) (6) 其中 (7) 式中：S为(Xt,Vt)的样本协方差矩阵；SX为Xt的m×m阶对称样本方差阵；SXV为Xt
与Vt的m×mP阶样本协方差阵；SV为Vt的mP×mP阶对称样本方差阵。Vi=(Xi-1,Xi-2,…，
Xi-p)T，Xi=(x(1)i,x(2)i,…，x(m)i)T，Vi和Xi(i=P+1,P+2,…n)来自实测样本。 当det(SV)≠0时(一般满足)，S可分解为：
设SX-SXVS-1VSTXV=A，则S可作如下变化： (8)
(9) 将式(8)、(9)代入式(5)整理得： (10)
再将式(6)、(10)代入式(4)可得： (11) 其中， (12)
由式(11)知，条件密度函数(Xt／Vt)是n-P个m维高斯函数(均值向量Bi，方差矩阵C)的加权(权重为Wi)平均和。用式(11)可随机模拟Xt，其模拟式为： Xt=Bi+DEt (13) 式中：D为m×m阶标准差矩阵，C=DDT；Et是均值0，方差1的m维独立高斯随机变量。在条件Vt下，模拟序列Xt是来自条件概率密度函数式(11)的一个样本。到此，多变量P阶非参数模型建立完毕。
2.2 模型阶数P的确定 在K(·)给定后，建立非参数模型的关键是寻求最优h和确定P。P可由AIC准则确定。
2.3 模型算法 算法：①从实测资料中构造Xi和Vi(i=P+1,P+2,…,n)；②计算协方差矩阵S；③给Vt赋初值；④给定Vt，由式(12)计算抽样概率Wi；⑤以概率Wi抽样Xi；⑥按式(13)模拟Xt；⑦给Vt重新赋值，转向第④步，继续模拟；满足模拟数时停止。
3 多变量非参数随机模型在日流量随机模拟中的应用 3.1 流域及资料情况 溪洛渡水电站的修建不仅要涉及自身的防洪安全，而且还要兼顾下游城市宜宾的安危，因此研究溪洛渡洪水(屏山站)和宜宾—屏山区间(简称宜-屏区间)洪水将至关重要。屏山站有48年(1940～1987)日流量过程，宜-屏区间日流量由岷江高场站实测日流量(1940～1987)按面积比获得。以上述两站日流量过程为例尝试将多变量非参数随机模型用于径流随机模拟研究，以检验该类模型的可行性和有效性。
3.2 模型建立 (i=1,2，…，48；j=1,2,…，365)表示屏山站、宜-屏区间日流量。计算各截口(天)斜方差矩阵S，建立各截口两变量非参数模型。经计算，P=1，汛期(5月1日～10月31日)h=0.532，非汛期(11月1日～次年4月30日)h=0379。构造各截口对应的Xi,Vi。当j≥2时，Xi=(),Vi=( )T(i=1,2,…,48)；当j=1时，Xi=()
Vi=( )T(i=2,3,…，48)。由模型算法便可对屏山站日流量和宜-屏区间日流量同时进行随机模拟。
3.3 模型检验 建立各截口两变量非参数模型，随机模拟100个模拟样本(模拟样本容量与实测样本相同)，用短序列法[1]对模型进行实用性检验。
3.3.1 日流量过程截口统计参数检验 截口统计参数有截口均值、S、Cv、Cs、滞时1、2的自相关系数R1,R2和最大值(Max)、最小值(Min)。统计了各统计量在不同站(区)不同模型情况下，在两个均方差检验标准下的通过率，结果载于表1。表中显示截口各统计参数的通过率是很理想的，表明各模型能很好地反映实测日流量过程的截口统计特性。
表1 截口各统计参数通过率(%) 站名 均值 S Cv Cs R1 R2 Max Min 屏山 宜-屏 100 97.3 100 100 97.8 100 87.5 98.4 100 100 100 100 100 96.7 98.9 91.8
3.3.2 时段量统计参数检验 检验时段有1、3、7、15、30d共5种，检验参数有时段量均值、Cv、Cs。成果载于表2中，可以看出：均值、Cv、Cs几乎都控制在一个均方差检验标准下。可见模型对各站(区)时段量统计参数都保持得很好。
3.3.3 月径流统计参数检验 为进一步验证模型的实用性，本文对5～10月径流统计参数进行了检验。成果载于表3中。该表显示：除极个别Cs控制在两个均方差检验标准下外，其余都控制在一个均方差检验标准下。说明各模型能保持好月径流统计参数。
3.3.4 年最大日流量季节性变化检验 分别模拟4800年日流量过程，统计年最大日流量在各月出现的百分比，载于表4。由表4看出，各模型能反映各站(区)年最大日流量季节性变化特性。
表2 时段洪量统计参数实用性检验 站名 时段洪量 W1日/108m3 W3日/108m3 W7日/108m3 W15日/108m3 W30日/108m3 参数 均值 Cv Cs 均值 Cv Cs 均值 Cv Cs 均值 Cv Cs 均值 Cv Cs
屏
山 样本 14.8 0.24 0.79 42.9 0.23 0.78 93.4 0.23 0.89 182 0.22 0.87 321 0.22 0.79
模 均值 14.9 0.23 0.50 43.0 0.23 0.53 93.8 0.23 0.58 183 0.22 0.42 327 0.22 0.30