基于Copula函数的HW一致性检验修正方法

梁冀雨，林炳章

(南京信息工程大学应用水文气象研究院，南京 210044)

1 研究背景

相较于传统的单站频率分析，地区频率分析法通过整合运用某一特定地区内所有站点历史资料序列的总体平均特征，优化对该区域内各站点的频率分析结果，使之更为准确、稳健[1]。在使用地区频率分析法时，需满足假设：研究区域内所有站点除了有一个不同的尺度系数外，频率分布的线型和参数应保持一定程度的一致[2]。满足上述假设时，称该研究区域为一致区。一致区的划分是地区频率分析的基础，区域的一致性将直接影响频率估计值的准确性[3]。因此,准确有效的区域一致性判别方法具有重要的实践意义和理论价值。

在地区频率分析法的早期应用中，Dalrymple[4]于1960年提出了一种基于样本离差系数(Cv)或偏态系数(Cs)的区域一致性检验方法，被大量应用于早期的降雨以及径流的地区频率分析中。但是该法对经验阈值的依赖较大，导致在应用于某些特定类型的洪水资料序列时，对异质站点的分辨能力较差。为了解决这一问题，Wiltshire[5]在1984年提出两种改进方案，其一仍然基于样本离差系数，其二基于区域的频率分布曲线。前者因与分布无关，适用范围更广。之后，Lu等[6]于1992年利用归一化的概化极值分布(GEV)的方差特性构建了新的一致性检验方法。在这些一致性检验方法的基础上，Hosking和Wallis[7]在1993年提出了至今应用最为广泛的Hosking-Wallis一致性检验(后文简称HW检验)。HW一致性检验以样本线性矩离差系数(L-Cv)和线性矩偏态系数(L-Cs)代替传统的离差系数和偏态系数，分析研究区域内各站点资料序列线性矩系数的离散程度，能够很好地与地区线性矩法(RLMA)相结合，对异质性的鉴别能力较强。

虽然HW一致性检验至今仍被大量应用，但是却存在理论上的缺陷。Hosking和Wallis在先前的研究中发现，站点资料间的相关性会降低地区频率分析中一致性检验的效力[3],其后2人提出的HW检验也未针对该缺陷进行调整。Castellarin等[8]发现HW检验的效力也会随着资料相关性的增大而减小，即当资料存在相关性时，HW检验可能将存在异质站点的区域判定为一致。Castellarin等给出的解决方案是一个由站点平均相关系数确定的经验修正系数，但该经验修正系数在每次使用时需事先根据资料律定参数，无法直接应用。

本文根据HW检验的理论基础，重点分析资料相关性影响其效力的原因，并使用蒙特卡罗模拟法揭示资料相关性对HW检验的影响。同时利用江西省年最大降雨量序列(AMS)研究实际资料的相关性特征，通过结合资料和理论，给出一种不受资料相关性影响的一致性检验。最后再次通过蒙特卡罗模拟对新的一致性检验方法进行评估。结果显示，该法可以为地区频率分析法，特别是地区线性矩法提供更有效的理论支持。

2 Hosking-Wallis一致性检验

为了将线性矩[9]应用于地区频率分析，Hosking和Wallis在1993年提出了基于资料线性矩系数的HW检验。HW检验将研究区域内各站点资料的各阶线性矩系数的离散程度，与符合一致区定义的、与研究区域具有相同站点数以及资料长度的人工模拟资料相比较。2者的差别若在一定的范围之内，就认定研究区域为一致区。

(1)

于是,区域内各站点的平均加权标准差V为：

(2)

最后进行蒙特卡罗模拟，生成500组模拟区域数据，根据式(2)计算得到各组模拟数据的平均加权标准差Vsim。利用Vsim的均值μV和标准差δV，对研究区域的平均加权标准差V进行类正态化处理，得到HW检验统计量H：

(3)

为了得到μV、σV，需进行蒙特卡洛模拟。首先根据式(1)，使用研究区域各站点资料计算得到的前4阶地区平均线性矩系数，拟合4参数的kappa分布[10]。接着模拟生成500个人工区域，每个人工区域都拥有和研究区域相同的站点数、资料长度，同时符合先前得到的kappa分布。值得注意的是，由此得到的每个人工区域都符合一致区的假设，并且站点间不存在相关关系。最后，对每一个人工区域分别计算平均加权标准差Vsim，并得到500个Vsim的均值和标准差。

Hosking和Wallis认为[9]，当H小于1时，研究区域为可接受的一致区；当H大于1并小于2时，研究区为可能的异质区域；而当H大于2时，研究区域为确定的异质区域。

另外，在考虑线性偏态系数和线性峰度系数的情况下，V具有2种变式，即：

(4)

式中：V2、V3分别代表同时考虑线性离差系数、线性偏态系数，以及线性偏态系数、线性峰度系数的情况。

此时，为区别3者不同，式(2)中的V改称V1。而由3者得到的H对应记作H1、H2和H3。事实上，3者结构相似，计算得到的H值往往代表相同的结果，因此Hosking等[9]建议多数情况下可仅考虑H1的结果。其后，Viglione等[11]发现，H1相较于H2以及H3，具有更强的鉴别异质性(heterogeneity)能力，且表现更为稳定。综上，本文仅对H1进行分析，相应的统计量，后文简称H以及V。

3 实测与生成雨量资料的相关性

相关系数是表征2个资料序列间相关程度的统计量。其中最为常用的是统计学家K. Pearson于1895年首次提出的皮尔森积矩相关系数，简称相关系数。对于2个成对的随机变量序列X、Y，其皮尔森相关系数的表达式为：

(5)

式中：rX,Y、cov(X,Y)分别表示随机序列X、Y的皮尔森相关系数和协方差。

由式(5)可知，皮尔森相关系数可以被看作2个随机变量之间标准化的协方差，当随机序列X、Y为标准正态分布时，2者的标准差都为1，式(5)右边可简化为cov(X,Y)。因此，选择皮尔森相关系数作为相关性的度量能够在进行后文的蒙特卡罗模拟实验时，通过分解协方差矩阵的方法[12]，生成具有特定相关性的资料序列。

为了比较HW检验中生成的模拟年最大降雨量资料与实测资料相关性的区别，本文以与HW检验相同的设置[9]进行蒙特卡罗模拟实验。从15 a到50 a，每一个序列长度生成1 000个人造区域，每个人造区域都包含100个站点。为做到与实测资料一致，令每个站点的模拟资料都服从4参数的kappa分布[10]，并拥有江西全省平均的24 h年最大雨量资料线性矩系数：t=0.18，t3=0.20，t4=0.15。同时因对资料进行归一化处理，所以所有生成资料的均值为1。之后，求出每个人工区域的站点间平均相关系数，最后计算每个序列长度1 000个人工区域平均相关系数的均值，绘于图1。

图1 实测与模拟雨量序列平均相关系数随序列长度的分布

图1中，圆形标记为实测年最大雨量平均相关系数随资料长度的分布，三角标记为HW检验中使用的模拟资料的情况。根据图1，模拟相较实测资料折线图下降趋势更为明显、平滑，且始终位于实测资料的下方，表明实测资料中普遍存在的相关性大于HW检验中生成的与之比较的模拟数据，且相差达到25%以上。

4 资料相关性对Hosking-Wallis一致性检验结果的影响

4.1 蒙特卡罗模拟中年最大雨量相关性对H统计量的影响

为了定量研究年最大雨量资料的相关性对H统计量的影

图2 资料相关性对Hosking-Wallis一致性检验的影响

由图2可知，相关性对不同的地区分布的影响程度是大致相同的，资料间的平均相关系数越大，H统计量曲线的位置越靠下，其值越小。当区域平均相关系数从0.2增加到0.8时，H值减小量均达到1.5以上，已经超过了判定有效一致区的H值区间[0,1]的长度，这大大降低了HW检验对资料相关区域异质性的辨别能力。

4.2 不同站点间相关性空间分布情况对H统计量的影响

(5)

式中：R为资料的相关系数矩阵，其形式决定了生成的资料间相关系数的均一的。

(6)

图3展示了地区分布为GEV，区域站点数为66时，不同相关系数分布情况对H值的影响。图3的横纵坐标与图2相同，改进后的模拟结果以圆点虚线分别绘于2张图上。图3(a)中黑色实线和三角形点虚线分别表示以式(5)的形式生成的平均相关系数为0.4、0.8的10 000个模拟区域的HW检验结果的均值，与图2(e)相同；而圆点虚线代表的人造区域中，仅有1/2的站点间具有0.8的相关系数，根据式(6)，其区域平均相关系数为0.4。图3(b)与图3(a)类似，其圆点虚线代表的人造区域中 的站点间拥有0.6的相关系数，区域平均相关系数为0.2。图3(a)、(b)中，圆点虚线几乎与黑色实线重合，而与三角形点虚线相距较远。由此不难发现，在区域平均相关系数相同时，相关性对H值的影响是近似相同的，而与站点间相关性空间分布无关。

基本指令是可编程控制编程中的重要组成部分，有 LD、LD NOT、AND、AND NOT、OR 等。基本指令多用于简单控制，编程较为灵活；对于繁琐的控制，若仅使用基本指令，编程较复杂，可用功能指令配合编程使其程序简化，且所编的程序较易阅读。

图3 不同相关系数分布对H值的影响

5 对Hosking-Wallis一致性检验的修正

5.1 修正Hosking-Wallis检验的算法及步骤

由第4节知,HW检验的偏差主要来自研究区域资料间固有的相关性与检验中生成的人工资料序列间相关性的不同。而即使站点资料序列间相关系数分布情况不同，只要整个区域的平均相关系数相同，其对HW检验结果的影响就基本相同。基于以上结论，对HW检验进行修正，使检验中生成的人工资料拥有与研究区域相同的平均相关系数，改进后的HW检验步骤如下。

(7)

式中：ρij为研究区域相关系数矩阵第i行第j列的元素；N为区域站点数。

(2)由式(5)得到N维平均相关系数矩阵R。

(3)计算研究区域各站点的各阶线性系数和区域的各阶平均线性离差系数。

(4)根据式(2)计算统计量V。

(5)生成与研究区域相同站点数、资料长度的人工区域资料yik，i=1,2,…,N，k=1,2,…,ni。ni为第i个站点的资料长度。运用正态Copula函数，使人工区域中站点资料服从多元正态分布，拥有均值零，协方差矩阵R。

(6)由(3)中求得的地区平均线性矩系数拟合4参数kappa分布。将资料yik转化为符合拟合的kappa分布的资料Qik：

Qik=KAP[φ(yik)]

(8)

式中：KAP为拟合得到的kappa函数的分位数函数；Φ为标准正态分布的累积分布函数。

(7)对人工区域资料序列Qik计算Vsim。

(8)重复1 000次步骤(5)～(7)，计算得到的1 000个Vsim的均值μV和标准差σV，并根据式(3)计算得到修正后的统计量H*。

5.2 资料相关性对修正前后Hosking-Wallis检验的影响

图4 改进前后HW检验结果对比

6 结 语

(1)利用蒙特卡洛模拟实验，得出HW一致性检验中生成的人工资料序列间的平均相关系数随资料长度的分布，并将其与实测资料进行对比。结果表明人工资料间的相关性均大于实测资料25%以上。

(2)从HW检验的定义出发，定性分析资料相关性影响H统计量值的原因，并通过蒙特卡罗模拟实验，定量分析资料相关性以及不同的相关性分布对H值大小的影响。定量分析发现当资料间存在较大相关性时，HW检验对异质性的检验能力大幅下降，而当区域平均相关系数相同时，相关性导致H值的减小量相同。

(3)针对雨量资料间的平均相关系数，对HW检验进行修正。实践证实，资料相关性对修正后HW检验结果的影响大大降低，同时修正后的HW检验依旧保留了其对区域一致性的鉴别能力。这种新的一致性检验方法适用条件等同HW检验，无需额外调整，计算速度快，可直接运用于地区线性矩频率法，帮助划分水文一致区。相较传统HW检验更为符合实测雨量资料情况、检验结果更加可靠。

(4)修正后的一致性检验通过优化水文一致区划分，提升了针对各水文时间序列资料，特别是降雨、径流序列频率分析的效率及其结果的准确性。将改进后的降雨、径流频率分析结果与气象降水量预报或流域水文模型径流量预报相比较，能够更为准确地划定预报降水或洪水的量级、重现期，为洪水预警预报提供更加坚实的理论基础。

□