一种基于雷达降雨数据挖掘的雨量站观测网络设计方法

技术领域

本发明涉及雨量站网络优化计技术，尤其涉及一种基于雷达降雨数据挖掘的雨量站观测网络设计方法。

背景技术

雨量站测得的降雨量是校准遥感降雨量和验证数值天气模型降雨产品的重要信息来源，然而由于维护和使用雨量站的成本过高，导致世界上许多地方的雨量站数量显著减少。因此，需要一种有效可行的雨量站网络设计方案，以最少的雨量站获取最大的降雨信息。

目前，世界上大多数地区仍然没有用于雨量站网络设计的标准程序。而传统的雨量站网络设计一般分为两种：偶然性方式和定量方式。偶然性的方式通常基于许多技术准则或考虑因素，包括确定雨量站数量及其位置，流域的性质，地形影响，排水方式，拟议地点的可及性和适用性，仪表的安装和维护成本，网络的应用目的和区域气候。此外，为了代表高变异性和间歇性的降雨模式，雨量站网络设计一般都十分密集。定量分析的雨量站网络设计通常基于分析可用的有限雨量站信息或复制成熟的网络。使用的方法有许多，比如空间相关性，变异函数分析和熵论；或使用统计技术，例如方差减少算法，状态空间随机模型和广义最小二乘法。但由于这些方法的性质和各个地区的条件差异，这些方法难以实施且经常需要主观调整参数，同时它们需要用所有可能的候选网络来探索最佳方法，成本较高。

发明内容

发明目的：针对以上问题，本发明提出一种基于雷达降雨数据挖掘的雨量站观测网络设计方法。

技术方案：为实现本发明的目的，本发明所采用的技术方案是：一种基于雷达降雨数据挖掘的雨量站观测网络设计方法，具体方法如下：

(1)获取研究区域在一段时间内N个降雨事件的降雨数据；具体方法如下：

获取研究区域在一段时间内N个降雨事件的降雨数据；其中，每一个降雨事件均得到已有雷达降雨观测格网的p次观测降雨数据；每一个降雨事件的降雨数据均用数据集X＝[X

(2)对步骤(1)所述N个降雨事件的降雨数据均进行主成分分析，由此得到每一个降雨事件所对应的主成分方差，进而确定已有雷达降雨观测格网的冗余度是否符合要求；若该冗余度符合要求，执行步骤(4)；否则，执行步骤(3)；

(3)分别对N个降雨事件迭代计算主成分方差，由此得到每一个降雨事件的最佳主成分个数，并在N个降雨事件的最佳主成分个数中选择主成分个数最多的值作为雨量站格网的最佳数量M；

(4)使用k-means聚类分析，将雨量站格网分配到聚类子集中；

(5)在每个聚类中均保留一个最佳网格位置，并设置网格中心为雨量站最佳位置，得到雨量站观测网络；

(6)对步骤(5)所述得到的雨量站观测网络进行判断；若判断结果为合格，则表示雨量站观测网络符合要求；否则，表示雨量站观测网络不符合要求，返回步骤(3)重新确定雨量站格网的最佳数量。

进一步地，步骤(2)中任意一个降雨事件通过主成分分析，得到相应的主成分方差的方法具体如下：

(2.1)将数据集X进行归一化，使数据集X的每一行减去其对应的均值，得到归一化后的数据集X'；计算公式：

Xi′＝X

式中，X'表示数据集X归一化后的数据集；X

(2.2)对降雨数据集X'计算得到协方差矩阵C；其中协方差矩阵C的第(i，j)项的计算公式如下：

∑i,j＝cov(X

式中，X

(2.3)对协方差矩阵C进行特征值分解，并计算其相应的特征向量，将特征值按照降序排列，同时特征值对应的特征向量作为列向量组成特征向量矩阵；

其中，特征值分解公式为：

[C-λ

式中，λ

特征向量计算公式为：

F

式中，F

更进一步地，步骤(2)确定已有雷达降雨观测格网的冗余度是否符合要求，具体方法如下：

(2.4)分别对N个降雨事件计算方差贡献率以及累计方差贡献率；

其中，降雨事件的第i主成分方差贡献率通过该降雨事件的第i个特征值除以该降雨事件的所有特征值的和得到，即

累计方差贡献率是通过该降雨事件的前i个主成分方差贡献率之和得到；

(2.5)预设累计方差贡献率阈值范围，当达到累计方差阈值范围时，分别计算N个降雨事件的主成分个数；

(2.6)将步骤(2.5)所述N个降雨事件的主成分个数分别与已有雷达降雨观测格网的数量进行比较，统计主成分个数小于已有雷达降雨观测格网数量的降雨事件数m；

若

优选地，为了更加方便快捷的判断出已有雷达降雨观测格网的冗余度是否符合要求，可以从N个降雨事件中随机选择L个降雨事件，并分别将L个降雨事件根据步骤(2)进行主成分分析，得到每一个降雨事件所对应的主成分方差，进而确定出已有雷达降雨观测格网的冗余度是否符合要求；

进一步地，步骤(3)对任意一个降雨事件迭代计算主成分方差，得到该降雨事件对应的最佳主成分个数的方法具体如下：

(3.1)迭代初始化，设置k＝1；

其中，k表示最佳主成分个数；

(3.2)计算前k个主成分的累计方差贡献率；

(3.3)判断步骤(3.2)所述前k个主成分的累计方差贡献率是否小于方差阈值；所述方差阈值是根据主成分方差曲线设置得到的；

若前k个主成分的累计方差贡献率小于方差阈值，k＝k+1，返回步骤(3.2)重新计算前k个主成分的累计方差贡献率；否则，得到最佳主成分个数k；

进一步地，步骤(4)的具体方法如下：

其中，聚类数量为M，即雨量站观测格网的数量；

(4.1)将研究区域内雷达观测网络随机分为M组，并将每一组的中心位置设置为该组的质心，即聚类质心，完成初始化；

(4.2)通过欧式距离计算研究区域内已有雷达观测网络中心与聚类质心之间的距离，并按照距离大小将雷达观测网络中心分配到最近的聚类中；

(4.3)按照距离将研究区域所有雷达观测网络中心均归属到最近的聚类之后，对聚类质心进行更新，将新的聚类的中心位置设为该聚类的质心；

(4.4)重复步骤(4.2)-(4.3)，直至达到收敛，即新的质心与原质心距离小于阈值或计算时间达到最大迭代，由此将研究区域的雷达降雨观测格网全部分配到聚类子集中；

进一步地，步骤(5)所述在每个聚类中均保留一个最佳网格位置的方法具体包括：

第一方法：分别对聚类中的每一个雷达网格计算p次观测降雨数据的平均值，获得各雷达网格的平均值，从中选择平均值最大的雷达网格作为最佳网格位置；

第二方法：分别对聚类中的每一个雷达网格计算p次观测降雨数据的平均值，获得各雷达网格的平均值；根据各雷达网格的平均值对各雷达网格进行降序排列，将排在中间位置的雷达网格作为最佳网格位置；若中间位置的雷达网格有两个，则在其中任意选择一个雷达网格作为最佳网格位置；其中，所述最佳网格位置的中心即为最佳雨量站位置。

进一步地，步骤(6)的具体方法如下：

通过皮尔逊相关系数和纳什-苏克利夫系数判断步骤(5)得到的最佳雨量站位置是否符合要求；

预先设置系数阈值，若计算得到的皮尔逊相关系数和纳什-苏克利夫系数均大于其对应的阈值，则表示步骤(5)所述雨量站观测网络符合要求；

否则，表示步骤(5)所述雨量站观测网络不符合要求，返回步骤(3)重新确定雨量站格网的最佳数量；

更进一步地，所述皮尔逊相关系数关系式：

式中，Rm为利用步骤(5)所述雨量站观测网络观测得到的降雨量；Ro为已有雷达降雨观测网络观测得到的降雨量；

皮尔逊相关系数越大，表示步骤(5)所述雨量站观测网络与已有雷达降雨观测网络的系统偏差越小，则设计得到的雨量站观测网络越优；

进一步地，所述纳什-苏克利夫系数表达式：

式中，R

其中，纳什-苏克利夫系数越接近1，表示设计得到的雨量站观测网络越优；

进一步地，利用信息分量曲线保证雷达降雨观测格网的数量的增加不会大幅增加降雨信息量。

有益效果：与现有技术相比，本发明技术方案具有以下有益技术效果：

与传统方法相比，采用雷达降雨测量方法进行雨量站网络设计更为直接和有效。本发明的数据来源为雷达降雨数据集，分析长期大范围的雷达降雨数据以此揭示区域降雨特征，可以在较大尺度上指导设计雨量站的布设网络；本发明充分考虑实际需求，设计得到的区域雨量站观测网络能够保证使用最少的雨量站得到相对最多的降雨信息，为雨量站的减少带来的降雨信息的缺失提供了一种新的解决方法；同时本发明所述评估雨量站观测网络的方法可以广泛应用于其他观测网络的设计评估。

附图说明

图1是本发明方法的流程图；

图2是随机抽取的10个降雨事件的主成分和其方差之间的关系图；

图3是通过平均降雨量最大的雷达格网得到的雨量站位置；

图4是通过平均降雨量中值的雷达格网得到的雨量站位置；

图5是信息方差曲线图。

具体实施方式

下面将结合附图对本发明技术方案的实施例进行详细的描述。以下实施例仅用于更加清楚地说明本发明的技术方案，因此只是作为示例，而不能以此来限制本发明的保护范围。需要注意的是，除非另有说明，本申请使用的技术术语或者科学术语应当为本发明所属领域技术人员所理解的通常意义。下面结合附图与具体实施例对本发明的做进一步详细说明：

本发明所述的一种基于雷达降雨数据挖掘的雨量站观测网络设计方法，具体方法如下：

(1)获取研究区域在一段时间内N个降雨事件的降雨数据；具体方法如下：

获取研究区域在一段时间内N个降雨事件的降雨数据，并对每一个降雨事件得到已有雷达降雨观测格网的p次观测降雨数据；每一个降雨事件的降雨数据均用数据集X＝[X

在本实施例中，获取研究区域一段时间内N＝373个降雨事件的降雨数据，并对每一个降雨事件的已有雷达降雨观测格网均获得p次观测数据，得到所有降雨事件所对应的数据集；

(2)对步骤(1)所述N个降雨事件的降雨数据均进行主成分分析，由此得到每一个降雨事件所对应的主成分方差，确定已有雷达降雨观测格网的冗余度是否符合要求；若该冗余度符合要求，执行步骤(4)；否则，执行步骤(3)；

其中，步骤(2)中任意一个降雨事件通过主成分分析，得到相应的主成分方差的方法具体如下：

(2.1)将数据集X进行归一化，使数据集X的每一行减去其对应的均值，得到归一化后的数据集X'；计算公式：

Xi′＝X

式中，X'表示数据集X归一化后的数据集；Xi′表示数据集X'的第i行；X

(2.2)对降雨数据集X'计算得到协方差矩阵C；其中协方差矩阵C的第(i，j)项的计算公式如下：

∑i,j＝cov(X

式中，X

(2.3)对协方差矩阵C进行特征值分解，并计算其相应的特征向量，将特征值按照降序排列，同时特征值对应的特征向量作为列向量组成特征向量矩阵；

其中，特征值分解公式为：

[C-λ

式中，λ

特征向量计算公式为：

F

式中，F

其中，步骤(2)确定已有雷达降雨观测格网的冗余度是否符合要求，具体方法如下：

(2.4)分别对N个降雨事件计算方差贡献率以及累计方差贡献率；

其中，降雨事件的第i主成分方差贡献率通过该降雨事件的第i个特征值除以该降雨事件的所有特征值的和得到，即

累计方差贡献率是通过该降雨事件的前i个主成分方差贡献率之和得到；

(2.5)预设累计方差贡献率阈值范围，当达到累计方差阈值范围时，分别计算N个降雨事件的主成分个数；

(2.6)将步骤(2.5)所述N个降雨事件的主成分个数分别与已有雷达降雨观测格网的数量进行比较，统计主成分个数小于已有雷达降雨观测格网数量的降雨事件数m；

若

在本实施例中，本应对步骤(1)所述研究区域一段时间内的N＝373个降雨事件均进行主成分分析，为方便快捷的判断已有雷达降雨观测格网的冗余度是否符合要求，从研究区域内的N＝373个降雨事件中随机选择10个降雨事件，根据步骤(2)所述方法判断已有雷达降雨观测格网数的冗余度是否符合要求；参考图2，对于大多数事件而言，第一个主成分已经携带了接近90％的总方差，第二个主成分方差和超过95％；这说明该网络的冗余度非常高，不符合要求，执行步骤(3)；

(3)分别对N个降雨事件迭代计算主成分方差，由此得到每一个降雨事件的最佳主成分个数，并在N个降雨事件的最佳主成分个数中选择主成分个数最多的值作为雨量站格网的最佳数量M；

其中，步骤(3)对任意一个降雨事件迭代计算主成分方差，得到该降雨事件对应的最佳主成分个数的方法具体如下：

(3.1)迭代初始化，设置k＝1；

其中，k表示最佳主成分个数；

(3.2)计算前k个主成分的累计方差贡献率；

(3.3)判断步骤(3.2)所述前k个主成分的累计方差贡献率是否小于方差阈值；所述方差阈值是根据主成分方差曲线设置得到的；

若前k个主成分的累计方差贡献率小于方差阈值，k＝k+1，返回步骤(3.2)重新计算前k个主成分的累计方差贡献率；否则，得到最佳主成分个数k；

在本实施例中，参考表1，为使所有降雨事件达到相应方差阈值，取最佳主成分个数最多的值，即为获得至少99％的信息，需要的主成分个数是12；

表1达到相应方差阈值时的主成分个数

(4)使用k-means聚类分析，将研究区域的雷达降雨观测格网分配到聚类子集中，具体方法如下：

其中聚类数量为M，即雨量站观测格网的数量；

(4.1)基于确定的聚类数量M，将研究区域内雷达观测网络随机分为M组，并将每一组的中心位置设置为该组的质心，即聚类质心，完成初始化；

(4.2)通过欧式距离计算研究区域内雷达观测网络中心与聚类质心之间的距离，并按照距离大小将雷达观测网络中心分配到最近的聚类中；

(4.3)按照距离将研究区域所有雷达观测网络中心均归属到最近的聚类之后，对聚类质心进行更新，将新的聚类的中心位置设为该聚类的质心；

(4.4)重复步骤(4.2)-(4.3)，直至达到收敛，即新的质心与原质心距离小于阈值或计算时间达到最大迭代，由此将研究区域的雷达降雨观测格网均分配到聚类子集中；

(5)在每一个聚类中均保留一个最佳网格位置，并设置网格中心为雨量站最佳位置，得到雨量站观测网络；

其中，在每一个聚类中保留一个最佳格网位置的方法包括两种，参考图3和图4，

第一方法：分别对聚类中的每一个雷达网格计算p次观测降雨数据的平均值，获得各雷达网格的平均值，从中选择平均值最大的雷达网格作为最佳网格位置；

第二方法：分别对聚类中的每一个雷达网格计算p次观测降雨数据的平均值，获得各雷达网格的平均值；根据各雷达网格的平均值对各雷达网格进行降序排列，将排在中间位置的雷达网格作为最佳网格位置；若中间位置的雷达网格有两个，则在其中任意选择一个雷达网格作为最佳网格位置；

其中，所述最佳网格位置的中心即为最佳雨量站位置；

(6)对步骤(5)两种方法下得到的雨量站观测网络进行判断是否符合要求，具体方法如下：

通过皮尔逊相关系数和纳什-苏克利夫系数判断步骤(5)得到的最佳雨量站位置是否符合要求，即是否能够保留原始雷达网络的主要降雨信息；

预先设置系数阈值，若计算得到的皮尔逊相关系数和纳什-苏克利夫系数均大于其对应的阈值，则表示步骤(5)所述雨量站观测网络符合要求；

否则，表示步骤(5)所述雨量站观测网络不符合要求，返回步骤(3)重新确定雨量站格网的最佳数量；

其中，皮尔逊相关系数关系式：

式中，Rm为利用步骤(5)所述雨量站观测网络观测得到的降雨量；Ro为已有雷达降雨观测网络观测得到的降雨量；

纳什-苏克利夫系数表达式：

式中，R

表2雨量站观测网络的皮尔逊相关系数和纳什-苏克利夫系数评价结果表

在本实施例中，参考表2，第一方法和第二方法下得到的皮尔逊相关系数值均较大，同时纳什-苏克利夫系数均接近于1；表示通过两种方法设计得到的雨量站观测格网都能够较好的保留原始雷达网络的主要降雨信息，均符合要求；但是两种方法下的皮尔逊相关系数相差较小，而第二方法下的纳什-苏克利夫系数相较于第一方法下的纳什-苏克利夫系数更接近于1，表示在该实施例中第二方法相较于第一方法更优。

进一步地，利用信息分量曲线保证雷达降雨观测格网的数量的增加不会大幅增加降雨信息量，参考图5，当观测格网个数大于12时并不能大幅增加降雨信息量，进而得知最终设计得到的雨量站观测网络验证通过。