本发明公开了一种基于水文气象一致性地区分析的雨量频率估计方法，属于水文气象技术领域。本发明方法首先根据气侯相似性和水文相似性的分析将研究区域划分为多个一致区，属于同一一致区的各雨量站点的气候统计特性相似、水文数据的统计特征值同一；然后利用一致区内所有雨量站点的历史数据分析本一致区最优的无量纲的雨量概率分布曲线，进而推求各个站点的雨量频率估计值。本发明能够有效提高雨量频率计算的精确性和准确性。
1.基于水文气象一致性地区分析的雨量频率估计方法，其特征在于，包括以下步骤：步骤A、根据气侯相似性和水文相似性的分析将研究区域划分为多个一致区，属于同一一致区的各雨量站点的气候统计特性相似、水文数据的统计特征值同一；具体包括：步骤A1、根据气侯特征相似性对研究区域进行初步分区，具体为：将满足水汽入流和水文气象成因背景一致的雨量站点划入同一个一级子区；
步骤A2、根据水文数据统计特征同一性将步骤A1得到的每个一级子区再细分为多个二级子区，二级子区即为一致区，具体按照以下方法：步骤A201、将每个一级子区划分为多个二级子区，使得每个二级子区的非一致性度量指标H1<1，非一致性度量指标H1的计算公式如下： ，
区中雨量站点的个数， 为该二级子区内第i 个雨量站点的历史雨量数据序列长度， 为R该二级子区中第i 个雨量站点的历史雨量数据序列的线性矩的离差系数；t 为该二级子区的平均的线性矩的离差系数；V1为该二级子区中所有雨量站点的历史雨量数据序列的线性矩的离差系数的标准差； 和 分别为由蒙特卡罗统计试验方法计算出的模拟V值的均值和均方差；
步骤A202、对每个二级子区，计算其中每个雨量站点的历史雨量数据序列的线性矩的偏态系数，从中选出特大值和特小值所对应的雨量站点，判断当去除该雨量站点时，应用“小概率事件在一次试验中不可能出现的统计特性”，结合比较其余各站点实测最大值与百年一遇估计值的大小，查看该二级子区内所有雨量站点的百年一遇估计值的合理性能否有改善，如有较大的改善就去掉该雨量站点、并将其归并至相邻二级区分析，反之保留；
步骤A203、对非一致性度量指标H1小于0的二级子区进行去相关操作，最后得到的二级子区即为一致区；
步骤B、首先将一致区内每个雨量站点的历史数据分解成共性和个性分量，然后对共性分量分析本一致区最优的无量纲的概率分布函数，进而结合各站的个性分量推求各个站点的雨量频率估计值，具体包括以下子步骤：步骤B1、按照以下方法确定各一致区的区域雨量频率因子：
步骤B101、将一致区内各雨量站点的历史雨量数据分解成共性分量和个性分量，即推求各个雨量数据的模比系数，此模比系数就是共性分量，而各站的均值就是个性分量；
步骤B102、对于一致区中每一雨量站点，计算该雨量站点历史雨量共性分量数据序列的线性矩的离差系数、偏态系数和峰度系数；
步骤B103、根据以下公式计算各一致区的共性分量的加权的区域线性矩的离差系数、偏态系数 和峰度系数 ， （r=3,4）
式中， 、 、 分别为该一致区中第i 个雨量站点的历史雨量共性分量数据序列的线性矩的离差系数、偏态系数和峰度系数，N 为该一致区内雨量站点的总数， 为该一致区内第i 个雨量站点的历史雨量共性分量数据序列长度；
步骤B104、根据区域线性矩系数和雨量概率分布函数参数之间的关系，利用样本数据推求出雨量概率分布函数的参数估计值，从而确定各一致区的区域雨量概率分布曲线；
步骤B105、根据适线准则，推选出最优的一种概率分布函数的参数估计值，从而确定各一致区的配合数据最佳的无量纲的概率分布曲线，此区域无量纲的分布曲线在不同频率下的估计值就是该一致区的雨量频率因子 ，它是代表一致区内共性因子的；
步骤B2、根据下式确定第j 个一致区中第i 个雨量站点在重现期为T时的雨量频率估计值 ：式中， 为第j个一致区在重现期为T时的地区雨量频率因子； 为第j个一致区中第i 个雨量站点年最大日降雨量的历史平均值，即该站的个性分量。
2.如权利要求1所述基于水文气象一致性地区分析的雨量频率估计方法，其特征在于，所述去相关操作具体按照以下方法：（1）选出所有站点间的距离小于30km的站点组；对于每个站点组中两个站点内的历史雨量数据序列，发生最大降雨量时间差前后大于一天的数据予以舍弃；对于历史雨量数据序列小于20年的站点组予以舍弃；
（2）计算剩余的每个站点组中两个雨量站点间的Pearson相关系数r，并选出|r|>0.7的站点组；
（3）对于选出的站点组，通过灵敏度检验来判断站点是否舍弃。
3.如权利要求1所述基于水文气象一致性地区分析的雨量频率估计方法，其特征在于，所述概率分布函数采用以下任一函数：概化的逻辑斯蒂分布函数、概化的极值分布函数、概化的正态分布函数、概化的帕累托分布函数、皮尔森Ⅲ型分布函数。
4.如权利要求1所述基于水文气象一致性地区分析的雨量频率估计方法，其特征在于，步骤A中所述研究区域包括实际研究区域以及实际研究区域周边的缓冲区。
基于水文气象一致性地区分析的雨量频率估计方法
技术领域
[0001] 本发明涉及一种雨量频率估计方法，属于水文气象技术领域。
背景技术
[0002] 我国是洪涝灾害最频发的国家之一，防洪始终是一个事关人民生命财产安全和社会稳定的非常重要的课题。无论是工程防洪设计，还是地区和都市的防洪设计，都需要一个科学可靠的防洪设计标准。我国的防洪设计标准估算的一个重要依据是洪水频率计算成果，其理论基础是水文频率计算。频率计算要解决的两大问题：精确性和准确性，即参数估算方法的精确性和频率估计值的准确性。而水文频率计算面临的两大困难是：无法通过理论分析导出概率分布曲线和计算的频率估计值的真值未知。如何来有效解决或者缓解这两个困难是水文统计学家所要面临的重大任务。
[0003] 在长期的频率计算的研究中，水文统计学家们发现，如果能找到一种参数估计方法，结合其他可操作的指标，在利用样本估计总体理论概率曲线的过程中，能够提供不偏、稳健的参数估计值，就有助于第一个问题精确性的解决。自20世纪90年代以来，Hosking在Greenwood定义的概率权重矩之后提出了线性矩，利用线性矩良好的无偏性和对特大值的稳健性，有效的解决或者缓解了第一个困难问题。Hosking在线性矩的基础上，提出了区域频率分析的方法，即“空间换时间”的思想来提高频率估计值的准确性，试图解决或缓解水文频率中的第二个困难。目前，一致区划分方法通常采用数理统计的聚类分析方法。这种方法的主导思想和准则是：欧氏距离越小，相似性越好，划分在一个区的可能性越大。但是，在水文上对一场暴雨而言，距离越小，站点间的极值暴雨的相关性就好，而水文频率分析的样本首先要满足独立同分布的条件，因此，数理统计就与实际水文事件产生相互矛盾。同时，进行聚类分析的因子变量一般都选取站点的经度、纬度、高程、多年平均降雨量等自然属性，而国内外许多学者的研究都表明：这些属性和极值暴雨很难有很好的相关关系，因此用这个四个因子变量来划分一致区，从物理意义上论据不足，利用分区计算的频率估计值的准确性虽比单站分析有所改善，但结果不是很令人满意，特别是各个分区间的频率估计值难于合理协调。我国的频率分析计算一直停留在“一点（单点、单时段）、一线（P-III型曲线）加双眼（目估适线）”的水平上，已远远不能满足国内快速发展的国民经济和生产建设的需求。因此，提出一种有效方法，来同时解决或者缓解水文频率分析计算中的两个难题，是防洪设计标准估算所要面临和解决的首要任务。
发明内容
[0004] 本发明所要解决的技术问题在于克服现有水文频率估计方法的不足，提供一种基于水文气象一致性地区分析的雨量频率估计方法，能够有效提高雨量频率估计值的准确性，为防洪设计提供准确的依据。
[0005] 本发明的基于水文气象一致性地区分析的雨量频率估计方法，包括以下步骤：
[0006] 步骤A、根据气侯相似性和水文相似性的分析将研究区域划分为多个一致区，属于同一一致区的各雨量站点的气候统计特性相似、水文数据的统计特征值同一；
[0007] 步骤B、首先将一致区内每个雨量站点的历史数据分解成共性和个性分量，然后对共性分量分析本一致区最优的无量纲的概率分布函数，进而结合各站的个性分量推求各个站点的雨量频率估计值，具体包括以下子步骤：
[0008] 步骤B1、按照以下方法确定各一致区的区域雨量频率因子：
[0009] 步骤B101、将一致区内各雨量站点的历史雨量数据分解成共性分量和个性分量，即推求各个雨量数据的模比系数，此模比系数就是共性分量，而各站的均值就是个性分量；
[0010] 步骤B102、对于一致区中每一雨量站点，计算该雨量站点历史雨量共性分量数据序列的线性矩的离差系数、偏态系数和峰度系数；
[0011] 步骤B103、根据以下公式计算各一致区的共性分量的加权的区域线性矩的离差系R R R数t、偏态系数t3 和峰度系数t4，
[0012]
[0013](i) (i) (i)
[0014] 式中，t 、t3 、t4 分别为该一致区中第i个雨量站点的历史雨量共性分量数据序列的线性矩的离差系数、偏态系数和峰度系数，N为该一致区内雨量站点的总数，ni为该一致区内第i个雨量站点的历史雨量共性分量数据序列长度；
[0015] 步骤B104、根据区域线性矩系数和雨量概率分布函数参数之间的关系，利用样本数据推求出雨量概率分布函数的参数估计值，从而确定各一致区的区域雨量概率分布曲线；
[0016] 步骤B105、根据适线准则，推选出最优的一种概率分布函数的参数估计值，从而确定各一致区的配合数据最佳的无量纲的概率分布曲线，此区域无量纲的分布曲线在不同频率下的估计值就是该一致区的雨量频率因子qT，j，它是代表一致区内共性因子的；
[0017] 步骤B2、根据下式确定第j个一致区中第i个雨量站点在重现期为T时的雨量频率估计值QT，i，j：
[0018]
[0019] 式中，qT,j为第j个一致区在重现期为T时的地区雨量频率因子； 为第j个一致区中第i个雨量站点年最大日降雨量的历史平均值，即该站的个性分量。
[0020] 上述技术方案中，一致区的划分可采用现有的各种方法，但为了进一步提高雨量估计的准确性，本发明进一步利用水文统计和气象成因相结合的途径进行一致区的划分，具体地，所述步骤A具体包括：
[0021] 步骤A1、根据气侯特征相似性对研究区域进行初步分区，具体为：将满足水汽入流和水文气象成因背景一致的雨量站点划入同一个一级子区；
[0022] 步骤A2、根据水文数据统计特征同一性将步骤A1得到的每个一级子区再细分为多个二级子区，二级子区即为一致区，具体按照以下方法：
[0023] 步骤A201、将每个一级子区划分为多个二级子区，使得每个二级子区的非一致性度量指标H1<1，非一致性度量指标H1的计算公式如下：
[0024]
[0025] 式中， N为该二级子(i)
区中雨量站点的个数，ni为该二级子区内第i个雨量站点的历史雨量数据序列长度，t 为R
该二级子区中第i个雨量站点的历史雨量数据序列的线性矩的离差系数；t 为该二级子区的平均的线性矩的离差系数；V1为该二级子区中所有雨量站点的历史雨量数据序列的线性矩的离差系数的标准差；μV和σV分别为由蒙特卡罗统计试验方法计算出的模拟V值的均值和均方差；
[0026] 步骤A202、对每个二级子区，计算其中每个雨量站点的历史雨量数据序列的线性矩的偏态系数，从中选出特大值和特小值所对应的雨量站点，判断当去除该雨量站点时，应用“小概率事件在一次试验中不可能出现的统计特性”，结合比较该站实测最大值与百年一遇估计值的大小，查看该二级子区内所有雨量站点的百年一遇估计值的合理性能否有改善，如有较大的改善就去掉该雨量站点、并将其归并至相邻二级区分析，反之保留；
[0027] 步骤A203、对非一致性度量指标H1小于0的二级子区进行去相关操作，最后得到的二级子区即为一致区。
[0028] 相比现有技术，本发明具有以下有益效果：
[0029] （1）本发明方法在进行概率分布曲线参数估计时，采用了线性矩法。和以往常用的常规矩法相比，该方法具有良好的无偏性和对特大值的稳健性等优点，尤其是当资料中出现特大值时，表现出其方法的强烈优势，推求的频率估计值相当稳定。基于线性矩法推求统计参数进行频率估计在工程设计中使得工程的安全得到有效保障，使得在这个基础上进行风险分析更加合理可靠。
[0030] （2）本发明方法在进行一致区划分时，采用了水文统计和气象成因相结合途径的地区线性矩法。和目前常用的聚类分析法和神经网络模型相比，该方法具有明确的物理意义，即在分区前，先把所有站点的降雨资料分成共性分量和个性分量两部分，或称去均值化，然后用共性分量进行地区分区，把个性分量留在本地，与实际水文现象符合；同时，考虑的影响因子比较全面，基于“小概率事件在一次试验中不可能出现的统计特性”的原理、应用“实际资料检验”（即经验频率与理论概率应该趋于一致）的思想，分析比较样本实测最大值与稀遇频率估计值关系的合理性，从水文统计和气象成因两个方面以及二者相结合的途径来考虑，给分区提供了牢固的基础；同时，能够充分利用研究区的其他站点的信息和有效利用缓冲区的信息，使得分区结果更合理、可靠。
[0031] （3）本发明将反映该一致区特点的无量纲频率估计值地区分量与反映各站点特有的降雨特性的个性分量进行“叠加”求出一致区内各站点的频率估计值，充分利用了区域内其它站点的历史资料信息，使得各个站点的雨量频率估计值都具有较高的精度和准确度。
具体实施方式
[0032] 在对本发明技术方案进行详细描述之前，先对本发明所涉及的线性矩的基本内容进行简要介绍。
[0033] Hosking将线性矩定义为次序统计量线性组合的期望值，其r阶线性矩λr的通用表达式为：
[0034]
[0035] 其中，EXr-k:r是对于一个样本容量为r的(r-k)阶次序统计量期望，其表达式为：
[0036]
[0037] 在参数估计中，一般仅用到前四阶线性矩，以下为样本线性矩前四阶表达式：
[0038]
[0039]
[0040]
[0041]
[0042] 其中n为样本容量，xm∶n为次序统计量(m＝i，j,k，l)。
[0043] 根据样本的前四阶线性矩（l2、l2、l3、l4）来求出线性矩系数（均值l1、比例系数l2、离差系数L-Cv=t=l2/l1、偏态系数L-Cs=t3=l3/l2、峰度系数L-Ck=t4=l4/l2）。
[0044] 本发明的雨量频率估计方法，利用了线性矩在样本参数估计中的良好的无偏性和对特大值的稳健性等优点，并将雨量分解为共性分量和个性分量，用共性分量进行水文气象一致区划分，在每个一致区内推求一条无量纲的、配合历史降雨资料点据最优的概率分布曲线，然后再与各个站点的个性分量结合求得各个站点的频率估计值。实践证明，该方法推求的频率估计值的稳定性（即准确度）有显著的提高。本发明方法具体包括以下步骤：
[0045] 步骤A、根据气侯相似性和水文相似性的分析将研究区域划分为多个一致区，属于同一一致区的各雨量站点的气候统计特性相似、水文数据的统计特征值同一。
[0046] 上述研究区域包括实际研究区域以及实际研究区域周边的缓冲区。缓冲区可根据实际研究区域的大小合理确定，缓冲区的距离一般取子区长度的一半即可。
[0047] 本步骤可采用现有的一致区划分方法，为了提高估计结果的准确性，本发明进一步利用水文统计和气象成因相结合的途径进行一致区的划分，具体包括以下子步骤：
[0048] 步骤A1、根据气侯特征相似性对研究区域进行初步分区，具体为：将满足水汽入流和水文气象成因背景一致的雨量站点划入同一个一级子区；
[0049] 步骤A2、根据水文数据统计特征同一性将步骤A1得到的每个一级子区再细分为多个二级子区，二级子区即为一致区，具体按照以下方法：
[0050] 步骤A201、将每个一级子区划分为多个二级子区，使得每个二级子区的非一致性度量指标H1<1，非一致性度量指标H1按Hosking建议的计算公式如下：
[0051]
[0052] 式中，
[0053]
[0054]
[0055] N为该二级子区中雨量站点的个数，ni为该二级子区内第i个雨量站点的历史雨(i)量数据序列长度，t 为该二级子区中第i个雨量站点的历史雨量数据序列的线性矩的离R
差系数；t 为该二级子区的平均的线性矩的离差系数；V1为该二级子区中所有雨量站点的历史雨量数据序列的线性矩的离差系数的标准差；μV和μV分别为由蒙特卡罗统计试验方法计算出的模拟V值的均值和均方差。
[0056] 蒙特卡罗统计试验方法为现有技术，其详细过程如下：假设二级子区N个站点中的每个站点都符合比较通用的、弹性较大的4参数的Kappa分布，以避免过早先入为主地设定某种概率分布，采用蒙特卡罗统计试验方法生成与站点实测资料相同长度的降雨资料，求得各站点每次模拟生成资料的离差系数，根据上式得到模拟值V，重复模拟Nsim（Nsim取500）次得到Nsim个模拟值V，根据Nsim个模拟样本计算模拟V的均值μV和均方差σV；即μV和μV分别为由蒙特卡罗统计试验方法计算出的V值的均值和方差。
[0057] 步骤A202、对每个二级子区，计算其中每个雨量站点的历史雨量数据序列的线性矩的偏态系数，从中选出特大值和特小值所对应的雨量站点，判断当去除该雨量站点时，应用“小概率事件在一次试验中不可能出现的统计特性”，结合比较该站实测最大值与百年一遇估计值的大小，查看该二级子区内所有雨量站点的百年一遇估计值的合理性能否有改善，如有较大的改善就去掉该雨量站点、并将其归并至相邻二级区分析，反之保留。
[0058] 线性偏态系数的判定也可以采用某种非一致性度量指标H来判定，但效果并不理想。本研究依据小概率事件在一次实验中出现的不可能性来判断线性偏态系数为特大值和特小值所对应的站点对整个子区的百年一遇的频率估计值的影响程度（即在某个仅具有几十年资料的站点的最大值居然大于百年一遇的稀遇事件），再决定该站点的去留或归属。由于中国所有站点的实测资料序列长度基本都小于100年，平均长度40～50年左右，如果百年一遇的频率估计值小于实测资料序列中的最大降雨值，这个站点就要引起注意。比如一个区有10个站点，当保留线性偏态系数为最大值的站点时，算出有5个站点的百年一遇的频率估计值小于相应站点实测序列的最大降雨值；当去掉偏态系数为最大值的站点时，算出有2个站点出现上述情况，那么就认为该站点对整个子区的的频率估计值的影响程度较大（约50%及以上），因此，综合考虑去掉偏态系数最大值对应的该站点，把它与相邻二级区重新组合分析。其它站点也是同样判断，反复进行，直到整个子区达到最优。这个过程体现了水文频率计算中分析重于计算的思想，因此这个分析过程对研究人员的专业素养和经验的要求比较高。
[0059] 步骤A203、对非一致性度量指标H1小于0的二级子区进行去相关操作，最后得到的二级子区即为一致区。
[0060] 本具体实施方式中，去相关操作具体如下:（1）制定三个筛选判定准则：第一，研究范围的确定，选出所有站点间的距离小于30km的站点组；第二，每组两个站点内的雨量资料，发生最大降雨量时间差前后大于一天的数据予以舍弃；第三，在第二个准则的基础上，对于资料序列小于20年的站点组予以舍弃。（2）相关系数的计算。选取水文计算中常用的Pearson相关系数r进行相关分析，当|r|>0.7时，认为此两个站资料间存在相关性。（3）对降雨资料存在较强相关性的站点进行灵敏度检验，来判断该站点是否舍弃。灵敏度检验是检验存在相关性的站点资料取舍与否对该一致区百年一遇的雨量估计值的影响程度，若某两个站点A、B之间的相关系数|r|≥0.7，计算其取舍与否对该一致区百年一遇的雨量估计值的影响程度，即采用相对误差RE来评估，
[0061]
[0062] 其中：Q0表示未去掉站点A、B时一致区百年一遇的降雨量估计值；Q1表示去掉站点A、B后一致区百年一遇的降雨量估计值。若相对误差RE<5%，则认为其对整个一致区的降雨估计值影响不大，可以考虑保留站点A、B资料；若相对误差RE>5%，则认为其资料对整个一致区的降雨估计值影响较大，需要去掉A、B两站中序列较短的站点资料，保留序列较长的站点资料。
[0063] 步骤B、利用一致区内所有雨量站点的历史数据分析本一致区内每一个雨量站点各自的雨量频率分布曲线，进而推求各个站点的雨量频率估计值，具体包括以下子步骤：
[0064] 步骤B1、按照以下方法确定各一致区的区域雨量频率因子：
[0065] 步骤B101、将一致区内各雨量站点的历史雨量数据（一般是年最大值）分解成共性分量和个性分量，即推求各个雨量数据的模比系数（即去均值化），此模比系数就是共性分量，而各站的均值就是个性分量；
[0066] 步骤B102、对于一致区中每一雨量站点，计算该雨量站点历史雨量共性分量数据序列的线性矩的离差系数、偏态系数和峰度系数；
[0067] 步骤B103、根据以下公式计算各一致区的共性分量的加权的区域线性矩的离差系R R R数t、偏态系数t3 和峰度系数t4，
[0068]
[0069]
[0070] 式中，t(i)、t3(i)、t4(i)分别为该一致区中第i个雨量站点的历史雨量共性分量数据序列的线性矩的离差系数、偏态系数和峰度系数，N为该一致区内雨量站点的总数，ni为该一致区内第i个雨量站点的历史雨量共性分量数据序列长度；
[0071] 步骤B104、根据区域线性矩系数和雨量概率分布函数参数之间的关系，利用样本数据推求出雨量概率分布函数的参数估计值，从而确定各一致区的区域雨量概率分布曲线；
[0072] 步骤B105、根据适线准则（goodness-of-fitting），推选出最优的一种概率分布（此分布实际上是无量纲的）函数的参数估计值，从而确定各一致区的配合数据最佳的无量纲的概率分布曲线，此区域无量纲的分布曲线在不同频率下的估计值就是该一致区的雨量频率因子qT，j，它是代表一致区内共性因子的；
[0073] 所述分布函数通常从与水文频率计算分析数据配合较好、呈正偏铃型的以下5种
