本发明属于流域水文分析技术领域，特别涉及了一种基于大数据统计相似分析的流域水文区划方法。
背景技术：
目前，对于水文分区的相关研究，如：罗开富等根据内外流域的分水线、水流形态和含沙量为指标将全国划分为外流区和内流区，汤奇成按气候带将全国划分为6个河流水文区，M.J.Hall和张静怡应用人工神经网络和模糊聚类对英国威尔士的西南地区及中国的江西省和福建省进行水文分区，熊怡等以径流的年内分配、径流深和径流动态为主要指标将全国划分为56个水文分区，而这些分区主要关注生态或水文单方面的因素，水文分区的需求和目的存在一定差别。
对于复杂的地区，受到多因素的影响，所以分区不能用一般的显示单方面因素的水文分区方法进行划分，需要考虑影响分区的所有的因子，而影响分区的因子众多，需要精简一些因子而又不失因子所代表的信息。
技术实现要素：
为了解决上述背景技术提出的技术问题，本发明旨在提供一种基于大数据统计相似分析的流域水文区划方法，从多元统计的角度出发，采用主成分分析生成多个综合因子来代表原来的指标，使其含有较多的信息，适用于多因素影响下的水文分区。
为了实现上述技术目的，本发明的技术方案为：
一种基于大数据统计相似分析的流域水文区划方法，包括以下步骤：
(1)将流域按照站点控制面积或者设定的流域面积标准，划分为n个子流域；
(2)选取气候因素、下垫面因素作为分区的p个相似因子；
(3)对气象站点的测量数据进行克里金插值，从而将测量数据覆盖整个流域；
(4)将p个相似因子取平均值，分别平铺在n个子流域，组成n行p列的相似因子矩阵；
(5)对相似因子矩阵进行降维和正交变换；
(6)对相似因子进行主成分分析，计算相似因子的相关系数矩阵，并得到相关系数矩阵的p个非负特征值；选出大于预设阈值的特征值，这些特征值所对应的特征向量即为主成分；
(7)将步骤(6)提取出的主成分作为K-Means聚类分析的聚类因子，根据最小类内方差准则，得到最优聚类数；
(8)将最优聚类数载入到划分的子流域中，生成水文相似性分区图；
(9)采用斯米尔诺夫一致性检验的方法，对分区结果进行验证。
进一步地，在步骤(2)中，所述气候因素包括降水、蒸发、气温和日照时数。
进一步地，在步骤(2)中，所述下垫面因素包括土地；利用类型、地形因子、土壤质地。
进一步地，在步骤(3)中，对相似因子矩阵X进行标准化，得到标准化后的相似因子矩阵Y＝(yij)n×p：
上式中，xij为相似因子矩阵X第i行第j列的元素，即第i个子流域的第j个相似因子。
进一步地，在步骤(6)中，所述预设阈值为1。
进一步地，在步骤(9)中，首先，选取分布在不同水文分区的气象站点进行斯米尔诺夫一致性检验，若检验结果表明，各气象站点间的上确界均大于临界值，则说明这些气象站点之间不相似，与分区结果相符；其次，选取分布在同一水文分区的气象站点进行斯米尔诺夫一致性检验，若检验结果表明，各气象站点间的上确界均小于临界值，则说明这些气象站点之间相似，与分区结果相符。
采用上述技术方案带来的有益效果：
本发明从多元统计的角度出发，采用主成分分析生成多个综合因子来代表原来的指标，使其含有较多的信息，精简了计算；利用K-Means聚类可以将因子间的相互作用以数值的形式体现出来，定量的进行分区，适合多因素影响下的水文分区。无资料流域水文站点稀疏，缺乏长序列的观测资料，通过水文分区则能实现有资料流域的成果向无资料流域的转化。
附图说明
图1是本发明的方法流程图。
具体实施方式
以下将结合附图，对本发明的技术方案进行详细说明。
本实施例选择的流域为雅鲁藏布江流域，雅江流域站点稀疏，属于典型的资料匮乏区。雅江发源于西藏西南部位于喜马拉雅山脉中段，我国境内流域干流总长约2104km，整体地势高亢，从西北向东南海拔逐渐降低，东南最低处为145m，河流总落差为5435m，是高海拔高落差的大河之一。流域属高原气候区，干湿季分明，气候垂直变异大。降水量由东南向西北逐渐减少，流域径流由降水、地下水和融雪(冰)水组成。
对雅江进行子流域划分时，划分面积太大，会忽略雅江复杂的地形和垂直的气候过程。本实施例以流域结构相对稳定的羊湖作为划分时流域面积的参考标准，将雅江流域划分为577个子流域，平均面积为525km2，相近于羊湖面积628km2。相似因子选择气候条件(如降水、蒸发、气温和日照时数)，下垫面因素(土地利用类型、地形因子、土壤质地)共计15个因子，并将其平均到子流域中，组成577行15列的因子矩阵。为了剔除相关性强的重复因子，需要对因子矩阵进行降维和分析，以发现关键因子。表1给出了各相似因子间的相关系数矩阵。总体上，相关系数介于-0.9～1之间，相关程度最大的是降水和日照时数，砂粒含量和粘粒含量，其次为降水和蒸发。
表1
根据计算得到的因子相关系数，进而求得因子的特征值，将特征值从大到小排列，将特征值大于1的作为主成分个数，可以确定主成分个数。表2给出了各因子的特征值。
表2
第1个主成分特征值为4.1，代表原来信息的27.33％的信息，前5个主成分特征值都大于1，可以保留原来信息的68.14％的信息。因此选取了5个主成分。表3为各个因子在该5个主成分上的综合因子值，当因子在这五个主成分上的荷载绝对值最大时，认为该主成分具有命名解释性。
表3
可以看出第一主成分解释了降水、蒸发和日照时数；第二主成份解释了气温、砂粒含量和粘粒含量；第三主成分解释了高山草甸、高山草地、荒漠草地和冰川；第四主成份解释了高程、针叶林和阔叶林；第五主成份解释了坡度和粉粒含量。
进行K-Means聚类分析时，将上一步主成分分析得到的5个综合因子作为K-Means聚类分析的聚类因子，根据最小类内方差准则，得到最优的聚类数。将聚类数载入到划分的子流域中，生成4个水文类型分区的空间分布图。
选取7个水文站进行斯米尔诺夫一致性检验。选取资料序列较长的拉萨站为代表站，由表4可以看出这7个水文站点的日平均雨量、多年平均径流深、多年平均径流模数和多年平均径流系数的上确界都小于临界值，说明代表站之间相似。
表4给出了斯米尔诺夫一致性检验的过程。
(1)设H代表两个事件相等，即Fx事件和Gx事件相等(H：Fx＝Gx)；
(2)这里以5个样本数据为例，即A1-A5和B1-B5是两个站点的数据；
(3)A(x)和B(x)为站点数据的累加，Fx和Gx为事件也叫经验分布函数；
(4)找出Fx-Gx这一列中最大的值λ；
(5)与最大值λ相对应的A、B、Ax和Bx，求得n值；
(6)在一定的显著性水平下(自己设定显著性水平值，查表得到该水平上的系数)，
(7)判断最大值和临界值之间的大小，决定接受H，还是拒绝H。如果最大值＜临界值，说明出自一个样本，一致性好，接受H；若最大值＞临界值，一致性差，拒绝H。
表4
为了使验证结果比较客观，对气象站点的资料也进行了检验。首先，选取分布在不同水文分区的气象站点进行检验，如分布在不同区的拉萨站、波密站和浪卡子站等。检验结果表明，各站点间的上确界均大于临界值，说明这些站点之间不相似，与分区结果相符。其次，选取分布在同一水文分区的拉萨、江孜、泽当和林芝站进行检验，分析结果显示拉萨与江孜和泽当的一致性检验良好。拉萨站与林芝站尽管处于同一分区，但上确界略大于临界值，表明其未通过一致性检验。这是由于林芝站位于两个分区的边界上，处于草地向林地过渡带，分区结果受子流域分辨率的影响，较粗的分辨率对边界单元内的水土和植被特征产生均化，从而影响分区结果。因此，除个别站点外，分区结果总体上通过了统计检验，其结果是可以接受的。表5给出了一致性检验的结果。
表5
图1即为整个实施例的主要流程图。
实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。