CN111008259B - 一种流域降雨相似性搜索方法 - Google Patents

Abstract

本发明公开了一种流域降雨相似性搜索方法，包括以下步骤：1）根据流域的轮廓范围生成栅格矩阵；2）对栅格矩阵进行编码；3）降雨事件划分；4）生成降雨矩阵样本集合；5）生成降雨事件相似性矩阵；6）流域降雨事件相似性搜索，最终获得相似的降雨事件。本发明能够同时考虑降雨的时间、空间分布，计算流域（区域）降雨事件相似性矩阵，利用相似性矩阵搜索相似的降雨事件，为流域的暴雨洪水预警及区域的内涝风险管理提供有力支撑。

Description

一种流域降雨相似性搜索方法

技术领域

本发明属于水利工程技术领域，尤其涉及防洪预报技术领域，具体为一种流域(区域)降雨相似性搜索方法。

背景技术

近年来我国极端暴雨事件频发，暴雨是诱发山洪及城市内涝的主要因素。对于暴雨事件，除雨量、雨强、持续时间等特征外，雨型及降雨空间分布作为对暴雨过程的描述，表现了暴雨强度在时间、空间尺度上的分配，同样是暴雨事件的主要致灾特征之一，即便具有相同的雨量及雨强，不同时程、空间分布的暴雨过程，致灾性截然不同。

对于诱发山洪的流域范围内的降雨，从降雨时程分布的角度出发，雨峰在中部或后部的雨型所产生的洪峰要远大于均匀形雨型所产生的洪峰，从降雨空间分布的角度出发，降雨中心由上游向下游移动的降雨过程所产生的洪峰要远大于降雨中心由下游向上游移动所产生的洪峰。因而，若简单的考虑降雨量及雨强去进行流域的洪水预警预报，会产生较大的误差。对于诱发城市内涝的区域范围内的降雨，从降雨时程分布的角度出发，雨峰位置不同所造成的城市内涝积水的最大范围和最大深度差异明显，从降雨空间分布的角度出发，不同降雨中心位置对内涝风险的分布有直接的影响，如何根据历史降雨过程分析频发的暴雨中心移动轨迹及停留位置，对于识别不同地区的内涝风险、加强内涝的风险管理具有重要意义。

目前在进行山洪临界雨量确定和城市内涝数值模拟时，降雨时程分布通常采用设计雨型，如芝加哥雨型、Huff雨型、Pilgrim雨型、Yen&Chow雨型等。降雨的空间分布往往被忽视，而采用面平均降雨过程。同时，设计雨型往往根据站点观测数据分析得到，利用点数据来代替面数据，难以表征降雨空间的分布特征。忽视降雨的空间分布及空间变化过程会使山洪临界雨量的确定以及城市内涝分析、排水***设计产生偏差，直接影响防洪治涝工作效果。

发明内容

本发明的目的在于针对此问题，提出一种流域(区域)降雨相似性搜索方法。本发明的目的是通过以下技术方案实现的：

一种流域降雨相似性搜索方法，包括以下步骤：

步骤1)根据流域的轮廓范围生成栅格矩阵；

步骤2)对栅格矩阵进行编码：通过对栅格进行编码获得三个编码矩阵，第一矩阵的元素为每个栅格中心点投影坐标X值；第二矩阵的元素为每个栅格中心点投影坐标Y值；第三矩阵的元素为栅格编号n，以Row代表矩阵行数，Col代表矩阵列数，n＝(i-1)*Col+j，i、j为矩阵栅格的行列索引，i＝1，...，Row；j＝1，...，Col；

步骤3)降雨事件划分：将记录降雨数据的连续的矩阵序列或一维时间序列根据降雨过程的起止时间进行分割；通过降雨事件划分，获得各降雨场次的开始时间点T_start_i与结束时间点T_end_i，并生成开始时间点序列{T_start_i，i＝1，...N}与结束时间点序列{T_end_i，i＝1，...N}，N为降雨事件个数；

步骤4)生成降雨矩阵样本集合：降雨矩阵样本集合包含N个样本元素，N为步骤3)获得的降雨事件个数，样本元素i对应第i个降雨事件，为长度为n_i的矩阵序列{Matrix_ij，j＝1，...，n_i}，j为矩阵序列中的时间索引，n_i为降雨事件i持续时长，n_i＝T_end_i-T_start_i，序列中每个矩阵的大小与步骤1)中划分的栅格矩阵相同；对于步骤3)中提取的一维时间序列，需利用站点权重矩阵序列{Matrix_k，k为流域内站点编号}将其转化为矩阵序列作为样本；

步骤5)生成降雨事件相似性矩阵：通过降雨事件相似性分析获得降雨事件相似性矩阵Matrix_P，矩阵的元素(i，j)为样本集中的第i个降雨事件与第j个降雨事件的相似度，以距离作为度量，距离越小则相似度越高；相似性矩阵Matrix_P的计算方式如下：

5-1.以步骤2)中的第三矩阵为索引，遍历栅格矩阵的每一个元素，根据流域边界和第一矩阵、第二矩阵中保存的栅格中心点投影坐标判断栅格中心点是否位于流域范围内，若位于边界内则参与计算，使用数组记录参与计算栅格的编号，记为Array_c；

5-2.通过对降雨矩阵样本集进行分析，获得嵌套权重矩阵Matrix_D，其大小为N×N，N为降雨事件个数，其元素(i，j)为降雨事件i与降雨事件j关联的权重矩阵，该权重矩阵的大小与步骤1)中的栅格矩阵相同，其计算方式为：对于第i个与第j个降雨事件，查询两场降雨中累积降雨量非零的栅格并取并集，设并集中栅格总数量为N_c，则对于第i个与第j个降雨事件，累积降雨量非零的栅格的距离权重为

其余栅格的距离权重为0，以矩阵的形式记录距离权重，依此方法完成嵌套权重矩阵Matrix_D各元素计算；

5-3.针对流域范围内的各栅格计算相似性矩阵序列{Matrix_m，m为Array_c中保存的栅格编号}，序列中每一个矩阵的大小均为N×N，N为降雨事件个数，矩阵Matrix_m的元素(i，j)为编号为m的栅格记录的第i个与第j个降雨事件的时间序列的相似度；根据5-1步中获得的编号数组Array_c遍历流域边界内的栅格，针对每一个栅格利用栅格记录的时间序列计算相似度，采用DTW距离作为相似度量，获得流域内各栅格对应的Matrix_m，生成矩阵序列{Matrix_m}；

5-4.计算降雨事件相似性矩阵Matrix_P，其大小为N×N，N为降雨事件个数，其元素(i，j)的计算方式为：根据5-1步得到的数组Array_C遍历流域边界内的栅格，根据Array_C中保存的栅格编号提取5-3步中得到的相似性矩阵序列{Matrix_m}中该栅格对应的Matrix_m乘以5-2步中获得的嵌套权重矩阵Matrix_D的元素(i，j)中记录的该栅格对应的权重，此权重值首先通过栅格编号查找步骤2)中的第三矩阵获得栅格行列号，进而通过行列号查找权重矩阵记录的数值获得，累加得到的矩阵的元素(i，j)即为矩阵Matrix_P的元素(i，j)值；依此方法计算各元素值得到降雨事件相似性矩阵Matrix_P，通过Matrix_P的元素(i，j)衡量样本集中的降雨事件i与降雨事件j的相似程度，值越小代表相似度越高；

步骤6)流域降雨事件相似性搜索：根据步骤5)中生成的降雨事件相似性矩阵Matrix_P，依次搜索矩阵中的最小值，最小值所在行列号对应的两个降雨事件即为最为相似的降雨事件。

进一步的，步骤3)中，对于雷达、卫星遥感和数值预报降雨数据，降雨事件划分可以通过限定阈值Th_p、Th_time、Th_int、Th_amount、Th_per来实现，其中阈值为Th_p降雨数值超过0的栅格占总栅格数目的百分比，阈值Th_time为降雨数值超过0的栅格的平均降雨时长，阈值Th_int为时间段，阈值Th_amount为降雨数值超过0栅格的平均降雨量，阈值Th_per为降雨数值超过0的栅格占总栅格数的百分比；具体流程为：当发生降雨的栅格百分比超过阈值Th_p时，认为降雨时间开始，降雨开始后，当超过时段Th_int未再有栅格发生降雨，则认为降雨事件结束，当占栅格总数百分比为Th_per的栅格的平均降雨量超过阈值Th_amount时，则保留此次降雨事件，否则作为微量降雨事件筛除，将平均时长低于阈值Th_time的降雨事件作为短时降雨筛除，从而筛选掉短时微量降雨。

进一步的，步骤3)中，对于站点记录的降雨数据，降雨事件划分可以通过限定阈值Th_p、Th_time、Th_int、Th_amount、Th_per来实现，其中阈值Th_p为降雨数值超过0的站点个数占总站点个数的百分比，阈值Th_time为降雨数值超过0的站点的平均降雨时长，阈值Th_int为时间段，阈值Th_amoumt为降雨数值超过0的站点的平均降雨量，阈值Th_per为降雨数值超过0的站点个数占总站点数的百分比；具体流程为：当发生降雨的站点数量百分比超过阈值Th_p时，认为降雨时间开始，降雨开始后，当超过时段Th_int未再有站点发生降雨，则认为降雨事件结束，当占站点总数百分比为Th_per的站点的平均降雨量超过阈值Th_amount时，则保留此次降雨事件，否则作为微量降雨事件筛除，将平均时长低于阈值Th_time的降雨事件作为短时降雨筛除，从而筛选掉短时微量降雨。

进一步的，步骤4)中，对于雷达、卫星遥感和数值预报降雨数据，其原始数据即为矩阵序列；对于雨量站点采集的数据，需要通过转换矩阵序列{Matrix_k，k为流域内站点编号}建立站点数据与栅格数据的转换关系，将站点观测时间序列转化为矩阵序列。

进一步的，对于利用权重矩阵序列将站点观测时间序列转化为矩阵序列，具体流程如下：

建立权重矩阵序列{Matrix_k，k为流域内站点编号}，其中的每个矩阵大小均与步骤1)中的栅格矩阵相同。序列中的第k个权重矩阵记录站点k的数据占各栅格数据的比重。以各栅格作为目标点，计算站点权重：首先在不考虑方向的情况下，选择流域内的站点利用反向距离权重法计算每个站点对目标点的距离权重：

w_k代表站点k的降雨量占目标点降雨量的权重。其中，k为站点索引，n为站点个数；d_k为站点至网格中心点的距离；

其中，(x，y)为网格中心点坐标；(x_k，y_k)为站点坐标，然后根据流域内站点的方向对权重进行修正，修正系数为：

其中l为除本站之外其余站点的索引，θ_j为站点1与站点k与目标点连线的夹角，修正后的总距离方向权重为：

W_k＝w_k(1+α_k)

即为站点k对目标点的距离权重，依次计算各站点对每个栅格目标点的距离方向权重保存于对应的权重矩阵中，从而获得权重矩阵序列{Matrix_k}。

获得权重矩阵序列后，利用流域内站点的降水量进行加权平均，得到目标网格的降水量：

其中，n为站点个数，k为站点索引，Pgrid为目标网格的降水量，P_k，obs为站点k的观测降水量。针对每个场次降雨，根据权重矩阵序列中保存的权重值利用上式将站点观测时间序列转化为栅格点时间序列，进而成为矩阵序列。

根据上述方法，计算每个栅格点所关联雨量站点的权重，根据权重值将站点观测时间序列转化为栅格点时间序列，进而成为矩阵序列。

进一步的，步骤5)中5-3步DTW距离的计算方法为：对时间序列X＝{x₁，x₂，...，x_i，...，x_m}和Y＝{y₁，y₂，...，y_i，...，y_n}，其中m、n分别代表两个时间序列的长度，寻找一个扭曲路径W来表示时间序列X与Y间的映射关系W＝{w₁，w₂，...，w_k，...，w_K}，max(n，m)≤K≤n+m-1，W的第k个元素记为w_k＝(i，j)，表示时间序列X的第i个元素与时间序列Y的第j个元素的对应关系，点(i，j)的累积距离计算公式为：γ(i，j)＝d(x_i，y_j)+min{γ(i-1，j-1)，γ(i-1，j)，γ(i，j-1)}，给定初始条件γ(1，1)＝d(x₁，y₁)，迭代计算得到累积距离矩阵，

即为时间序列X与Y的DTW距离。

本发明的有益效果：

本发明提出一种流域(区域)降雨相似性搜索方法，通过采集流域(区域)站点观测数据并自动划分为场次降雨过程，将流域划分为栅格计算单元，采用机器学习方法，能够同时考虑降雨的时间、空间分布，计算流域(区域)降雨事件相似性矩阵，利用相似性矩阵搜索相似的降雨事件。对于降雨事件，雨强、降雨量、持续时间是最常用的描述降雨事件的三个特征，而随着防洪治涝工作的日益精细化，对降雨过程的描述提出了更高的要求，已有研究表明，暴雨的致灾性不仅区分在雨强、降雨量、持续时间上，降雨的时程以及空间分布同样对暴雨洪水、内涝产生很大的影响。本发明提出一种降雨事件相似性的搜索方法，可以同时考虑降雨时间、空间、雨强等特征，评价降雨过程的相似性。基于此方法，可以根据历史数据归纳出流域(区域)常见的降雨过程特征，也可以分析当前或未来降雨可能的发展趋势和致灾性，从而提示洪水、内涝风险。结合此方法分析获得的山洪临界雨量以及城市暴雨内涝过程具有更强的合理性。

下面结合附图及具体实施方式对本发明作进一步详细说明。

附图说明

图1本发明方法流程示意图；

图2某城区轮廓及栅格矩阵；

图3某自然流域轮廓及栅格矩阵；

图4矩阵降雨数据示意图；

图5时间序列的动态扭曲路径；

图6区域轮廓、站点位置分布及栅格矩阵；

图7为2011年6月7日降雨过程；

图8为2015年9月30日降雨过程；

图9为降雨事件相似性矩阵；

图10为2016年9月25日降雨过程；

图11为2018年6月16日降雨过程。

具体实施方式

实施例1

一种流域降雨事件相似性搜索方法，包括以下步骤：

1)生成流域(区域)栅格矩阵：

如图2、3所示，根据流域(区域)的轮廓范围生成栅格矩阵，栅格矩阵生成的原则为矩阵范围能够覆盖流域(区域)的轮廓范围且无过多冗余栅格。采用等边长栅格，根据流域(区域)面积和栅格总数量控制栅格大小，栅格过大会导致空间特征分辨率不足，栅格过小会导致计算量过大，可以根据计算设备的性能合理选择栅格总数量。对于雷达、卫星遥感和数值预报降雨数据，通常此类数据的原始形式即为栅格，因而可以根据流域(区域)范围直接选择参与计算的栅格，并保留其原有的栅格大小。

2)对栅格矩阵进行编码：

通过对栅格进行编码获得三个编码矩阵，第一矩阵的元素为每个栅格中心点投影坐标X值；第二矩阵的元素为每个栅格中心点投影坐标Y值；第三矩阵的元素为栅格编号n，以Row代表矩阵行数，Col代表矩阵列数，n＝(i-1)*Col+j，i、j为矩阵栅格的行列索引，i＝1，...，Row；j＝1，...，Col。

3)降雨事件划分：

将记录降雨数据的连续的矩阵序列或一维时间序列根据降雨过程的起始时间进行分割，从而获得降雨事件集合。对于雷达、卫星遥感和数值预报降雨数据，降雨事件划分可以通过限定阈值Th_p、Th_time、Th_int、Th_amount、Th_per来实现，其中阈值为Th_p降雨数值超过0的栅格占总栅格数目的百分比，阈值Th_time为降雨数值超过0的栅格的平均降雨时长，阈值Th_int为时间段，阈值Th_amount为降雨数值超过0栅格的平均降雨量，阈值Th_per为降雨数值超过0的栅格占总栅格数的百分比；具体流程为：当发生降雨的栅格百分比超过阈值Th_p时，认为降雨时间开始，降雨开始后，当超过时段Th_int未再有栅格发生降雨，则认为降雨事件结束，当占栅格总数百分比为Th_per的栅格的平均降雨量超过阈值Th_amount时，则保留此次降雨事件，否则作为微量降雨事件筛除，将平均时长低于阈值Th_time的降雨事件作为短时降雨筛除，从而筛选掉短时微量降雨。

对于站点记录的降雨数据，降雨事件划分可以通过限定阈值Th_p、Th_time、Th_int、Th_amnount、Th_per来实现，其中阈值Th_p为降雨数值超过0的站点个数占总站点个数的百分比，阈值Th_time为降雨数值超过0的站点的平均降雨时长，阈值Th_int为时间段，阈值Th_amount为降雨数值超过0的站点的平均降雨量，阈值Th_per为降雨数值超过0的站点个数占总站点数的百分比；具体流程为：当发生降雨的站点数量百分比超过阈值Th_p时，认为降雨时间开始，降雨开始后，当超过时段Th_int未再有站点发生降雨，则认为降雨事件结束，当占站点总数百分比为Th_per的站点的平均降雨量超过阈值Th_amount时，则保留此次降雨事件，否则作为微量降雨事件筛除，将平均时长低于阈值Tg_time的降雨事件作为短时降雨筛除，从而筛选掉短时微量降雨。

通过降雨事件划分，获得各降雨场次的开始时间点T_start_i与结束时间点T_end_i，并生成开始时间点序列{T_start_i，i＝1，...N}与结束时间点序列{T_end_i，i＝1，...N}，N为降雨场次个数。

4)生成降雨矩阵样本集

设N为通过步骤3)降雨事件划分得到的降雨事件个数，则降雨矩阵样本集为一个包含N个样本元素的集合。样本元素i对应第i个降雨事件，为长度为n_i的矩阵序列{Matrix_ij，j＝1，...，n_i}，n_i为降雨事件i持续时长。根据步骤3)中生成的开始时间点序列{T_start_i，i＝1，...N}与结束时间点序列{T_end_i，i＝1，...N}，n_i＝T_end_i-T_start_i。矩阵序列的元素Matrix_ij为如步骤1)中所示的Row×Col的矩阵，代表降雨事件i在j时刻的降雨空间分布。

对于雷达、卫星遥感和数值预报降雨数据，其原始数据即为矩阵序列，如图4所示，矩阵样本集可基于步骤3)中获得的开始、结束时间点序列对原始的矩阵序列进行分割获得。

对于雨量站点采集的数据，首先需要通过权重矩阵序列建立站点数据与栅格数据的转换关系。对于权重矩阵序列数值的获得方式，此处推荐一种同时考虑距离与方向的影响的空间插值方法，但并不限于此方法。具体流程如下：

建立权重矩阵序列{Matrix_k，k为流域内站点编号}，其中的每个矩阵大小均与步骤1)中的栅格矩阵相同。序列中的第k个权重矩阵记录站点k的数据占各栅格数据的比重。以各栅格作为目标点，计算站点权重：首先在不考虑方向的情况下，选择流域内的站点利用反向距离权重法计算每个站点对目标点的距离权重：

w_k代表站点k的降雨量占目标点降雨量的权重，其中，k为站点索引，n为站点个数；d_k为站点至网格中心点的距离；

其中，(x，y)为网格中心点坐标；(x_k，y_k)为站点坐标，然后根据流域内站点的方向对权重进行修正，修正系数为：

其中为l为除本站之外其余站点的索引，θ_j为站点1与站点k与目标点连线的夹角，修正后的总距离方向权重为：

W_k＝w_k(1+α_k)

即为站点k对目标点的距离权重，依次计算各站点对每个栅格目标点的距离方向权重保存于对应的权重矩阵中，从而获得权重矩阵序列{Matrix_k，k为流域内站点编号}。

获得权重矩阵序列后，利用流域内站点的降水量进行加权平均，得到目标网格的降水量：

其中，n为站点个数，k为站点索引，P_grid为目标网格的降水量，P_k，obs为站点k的观测降水量；针对每个场次降雨，根据权重矩阵序列中保存的权重值利用上式将站点观测时间序列转化为栅格点时间序列，进而成为矩阵序列。

5)生成降雨事件相似性矩阵：

通过降雨事件相似性分析获得降雨事件相似性矩阵Matrix_P，矩阵大小为N×N，N为降雨矩阵集中的样本数量，矩阵的元素(i，j)为样本集中的第i个降雨事件与第j个降雨事件的相似度，以距离作为度量，距离越小则相似度越高，因而Matrix_P为主对角线元素为0的对称矩阵。

两个降雨事件相似度的分析方式如下：

5-1.以步骤2)中的第三矩阵为索引，遍历栅格矩阵的每一个元素，根据流域边界和第一矩阵、第二矩阵中保存的栅格中心点投影坐标判断栅格中心点是否位于流域范围内，若位于边界内则参与计算，使用数组记录参与计算栅格的编号，记为Array_c；

5-2.通过对降雨矩阵样本集进行分析，获得嵌套权重矩阵Matrix_D，其大小为N×N，N为降雨事件个数，其元素(i，j)为降雨事件i与降雨事件j关联的权重矩阵，该权重矩阵的大小与步骤1)中的栅格矩阵相同，其计算方式为：对于第i个与第j个降雨事件，查询两场降雨中累积降雨量非零的栅格并取并集，设并集中栅格总数量为N_c，则对于第i个与第j个降雨事件，累积降雨量非零的栅格的距离权重为

其余栅格的距离权重为0，以矩阵的形式记录距离权重，依此方法完成嵌套权重矩阵Matrix_D各元素计算；

5-3针对流域范围内的各栅格计算相似性矩阵序列{Matrix_m，m为Array_c中保存的栅格编号}，序列中每一个矩阵的大小均为N×N，N为降雨事件个数，矩阵Matrix_m的元素(i，j)为编号为m的栅格记录的第i个与第j个降雨事件的时间序列的相似度；根据5-1步中获得的编号数组Array_c遍历流域边界内的栅格，针对每一个栅格利用栅格记录的时间序列计算相似度，采用DTW距离作为相似度量，获得流域内各栅格对应的Matrix_m，生成矩阵序列{Matrix_m}；DTW距离计算方法如下：

对时间序列X＝{x₁，x₂，...，x_i，...，x_m}和Y＝{y₁，y₂，...，y_i，...，y_n}，寻找一个扭曲路径W来表示时间序列X与Y间的映射关系，如图5所示，W＝{w₁，w₂，...，w_k，...，w_K}，max(n，m)≤K≤n+m-1，W的第k个元素记为w_k＝(i，j)，表示时间序列X的第i个元素与时间序列Y的第j个元素的对应关系。扭曲路径的选取有三个约束条件：扭曲路径始于矩阵的起始元素，结束于对角元素，即w₁＝(1，1)，w_K＝(m，n)；扭曲路径每一步都是连续的，即对于w_k＝(a，b)，w_k-1＝(a′，b′)，要求a-a′≤1且b-b′≤1；扭曲路径在时间轴上是单调的，即对于w_k＝(a，b)，w_k-1＝(a′，b′)，要求a-a′≥0且b-b′≥0。

能够满足约束条件的路径有很多条，此处寻找扭曲代价最小的路径，即：

其中d(w_k)为w_k代表的两个对应元素间的距离。

根据动态规划思想，若点(i，j)在最佳路径上，那么从点(1，1)到点(i，j)的子路径也是局部最优解，即从点(1，1)到点(m，n)的最佳路径可以由起始点(1，1)到终点(m，n)之间的局部最优解递归搜索获得，因而可以方便地找到这个最佳路径。具体步骤为：首先构建一个m×n阶矩阵，矩阵元素(i，j)为两个时间序列点x_i和点y_j之间的距离d(x_i，y_j)＝(x_i-y_j)²。定义点(i，j)的累积距离计算公式：

γ(i，j)＝d(x_i，y_j)+min{γ(i-1，j-1)，γ(i-1，j)，γ(i，j-1)}

给定初始条件γ(1，1)＝d(x₁，y₁)，可以迭代计算得到累积距离矩阵。

即为时间序列X与Y的DTW距离，从点γ(m，n)出发反向搜索累积距离矩阵即可得到最佳匹配路径。

依照上述方法完成相似性矩阵序列{Matrix_m}的计算。

5-4.计算降雨事件相似性矩阵Matrix_P，其大小为N×N，N为降雨事件个数，其元素(i，j)的计算方式为：根据5-1步得到的数组Array_c遍历流域边界内的栅格，根据Array_c中保存的栅格编号提取5-3步中得到的相似性矩阵序列{Matrix_m}中该栅格对应的Matrix_m乘以5-2步中获得的嵌套权重矩阵Matrix_D的元素(i，j)中记录的该栅格对应的权重，此权重值首先通过栅格编号查找步骤2)中的第三矩阵获得栅格行列号，进而通过行列号查找权重矩阵记录的数值获得，累加得到的矩阵的元素(i，j)即为矩阵Matrix_P的元素(i，j)值；依此方法计算各元素值得到降雨事件相似性矩阵Matrix_P，通过Matrix_P的元素(i，j)衡量样本集中的降雨事件i与降雨事件j的相似程度，值越小代表相似度越高。

6)流域(区域)降雨事件相似性搜索

根据步骤5)中生成的降雨事件相似性矩阵Matrix_P，依次搜索矩阵中的最小值，最小值所在行列号对应的两个降雨事件即为最为相似的降雨事件。

针对我国北方某城市近年发生的降雨事件进行相似性搜索，共收集到城区49个雨量站点近10年的降雨数据，49个站点相对均匀地分布于城区，拟采用空间插值的方法将站点数据转化为栅格数据进行分析。

首先根据城市区域的轮廓生成计算栅格矩阵，基于区域范围的大小，选择栅格大小为5000m×5000m，得到栅格矩阵的行列数分别为36、33。区域轮廓、站点位置分布及计算栅格矩阵如图6所示：

对降雨事件进行划分：限定阈值Th_p为0.1，即降雨数值超过0的站点占总站点的百分比高于10％即认为降雨事件开始；Th_int为6，即若6小时内无站点监测到降雨，认为一次降雨事件结束；Th_time为3，即降雨数值超过0的站点的平均降雨时长大于3小时才考虑此次降雨事件；Th_amount分别设置为0、10与40，Th_per分别为0.5、0.2与0.01，即50％的站点监测到降雨或20％的站点监测降雨量超过10mm或1％的站点监测降雨量超过40mm即考虑此次降雨事件。依照上述阈值，共划分出181场城区降雨事件。

生成降雨矩阵样本集：采用同时考虑距离与方向的影响的空间插值方法将站点降雨数据转换为矩阵数据，生成181个降雨矩阵样本，如图7、8以2011年6月7日降雨事件、2015年9月30日降雨事件为例所示。

生成降雨事件相似性矩阵如图9所示。

依据降雨事件相似性矩阵，搜索矩阵中的最小值，判断得到2016年9月25日发生的降雨事件与2018年6月16日发生的降雨事件最为相似，两次降雨过程降雨中心均位于城区南部，且由西向东移动，降雨过程分别如图10、图11所示。

Claims (6)

1.一种流域降雨相似性搜索方法，其特征在于：包括以下步骤：

步骤1)根据流域的轮廓范围生成栅格矩阵；

步骤2)对栅格矩阵进行编码：通过对栅格进行编码获得三个编码矩阵，第一矩阵的元素为每个栅格中心点投影坐标X值；第二矩阵的元素为每个栅格中心点投影坐标Y值；第三矩阵的元素为栅格编号n，以Row代表矩阵行数，Col代表矩阵列数，n＝(i-1)*Col+j，i、j为矩阵栅格的行列索引，i＝1，...，Row；j＝1，...，Col；

步骤3)降雨事件划分：将记录降雨数据的连续的矩阵序列或一维时间序列根据降雨过程的起止时间进行分割；通过降雨事件划分，获得各降雨场次的开始时间点T_start_i与结束时间点T_end_i，并生成开始时间点序列{T_start_i，i＝1，...N}与结束时间点序列{T_end_i，i＝1，...N}，N为降雨事件个数；

步骤4)生成降雨矩阵样本集合：降雨矩阵样本集合包含N个样本元素，N为步骤3)获得的降雨事件个数，样本元素i对应第i个降雨事件，为长度为n_i的矩阵序列{Matrix_ij，j＝1，...，n_i}，j为矩阵序列中的时间索引，n_i为降雨事件i持续时长，n_i＝T_end_i-T_start_i，序列中每个矩阵的大小与步骤1)中划分的栅格矩阵相同；对于步骤3)中提取的一维时间序列，需利用站点权重矩阵序列{Matrix_k,k为流域内站点编号}将其转化为矩阵序列作为样本；

步骤5)生成降雨事件相似性矩阵：通过降雨事件相似性分析获得降雨事件相似性矩阵Matrix_P，矩阵的元素(i，j)为样本集中的第i个降雨事件与第j个降雨事件的相似度，以距离作为度量，距离越小则相似度越高；相似性矩阵Matrix_P的计算方式如下：

5-1.以步骤2)中的第三矩阵为索引，遍历栅格矩阵的每一个元素，根据流域边界和第一矩阵、第二矩阵中保存的栅格中心点投影坐标判断栅格中心点是否位于流域范围内，若位于边界内则参与计算，使用数组记录参与计算栅格的编号，记为Array_c；

5-2.通过对降雨矩阵样本集进行分析，获得嵌套权重矩阵Matrix_D，其大小为N×N，N为降雨事件个数，其元素(i，j)为降雨事件i与降雨事件j关联的权重矩阵，该权重矩阵的大小与步骤1)中的栅格矩阵相同，其计算方式为：对于第i个与第j个降雨事件，查询两场降雨中累积降雨量非零的栅格并取并集，设并集中栅格总数量为N_c，则对于第i个与第j个降雨事件，累积降雨量非零的栅格的距离权重为

其余栅格的距离权重为0，以矩阵的形式记录距离权重，依此方法完成嵌套权重矩阵Matrix_D各元素计算；

5-3.针对流域范围内的各栅格计算相似性矩阵序列{Matrix_m,m为Array_c中保存的栅格编号}，序列中每一个矩阵的大小均为N×N，N为降雨事件个数，矩阵Matrix_m的元素(i，j)为编号为m的栅格记录的第i个与第j个降雨事件的时间序列的相似度；根据5-1步中获得的编号数组Array_c遍历流域边界内的栅格，针对每一个栅格利用栅格记录的时间序列计算相似度，采用DTW距离作为相似度量，获得流域内各栅格对应的Matrix_m，生成矩阵序列{Matrix_m}；

5-4.计算降雨事件相似性矩阵Matrix_P，其大小为N×N，N为降雨事件个数，其元素(i，j)的计算方式为：根据5-1步得到的数组Array_C遍历流域边界内的栅格，根据Array_C中保存的栅格编号提取5-3步中得到的相似性矩阵序列{Matrix_m}中该栅格对应的Matrix_m乘以5-2步中获得的嵌套权重矩阵Matrix_D的元素(i，j)中记录的该栅格对应的权重，此权重值首先通过栅格编号查找步骤2)中的第三矩阵获得栅格行列号，进而通过行列号查找权重矩阵记录的数值获得，累加得到的矩阵的元素(i，j)即为矩阵Matrix_P的元素(i，j)值；依此方法计算各元素值得到降雨事件相似性矩阵Matrix_P，通过Matrix_P的元素(i，j)衡量样本集中的降雨事件i与降雨事件j的相似程度，值越小代表相似度越高；

步骤6)流域降雨事件相似性搜索：根据步骤5)中生成的降雨事件相似性矩阵Matrix_P，依次搜索矩阵中的最小值，最小值所在行列号对应的两个降雨事件即为最为相似的降雨事件。

2.根据权利要求1所述的一种流域降雨相似性搜索方法，其特征在于：步骤3)中，对于雷达、卫星遥感和数值预报降雨数据，降雨事件划分可以通过限定阈值Th_p、Th_time、Th_int、Th_amount、Th_per来实现，其中阈值为TH_p降雨数值超过0的栅格占总栅格数目的百分比，阈值Th_time为降雨数值超过0的栅格的平均降雨时长，阈值Th_int为时间段，阈值Th_amount为降雨数值超过0栅格的平均降雨量，阈值Th_per为降雨数值超过0的栅格占总栅格数的百分比；具体流程为：当发生降雨的栅格百分比超过阈值Th_p时，认为降雨时间开始，降雨开始后，当超过时段Th_int未再有栅格发生降雨，则认为降雨事件结束，当占栅格总数百分比为Th_per的栅格的平均降雨量超过阈值Th_amount时，则保留此次降雨事件，否则作为微量降雨事件筛除，将平均时长低于阈值Th_time的降雨事件作为短时降雨筛除，从而筛选掉短时微量降雨。

3.根据权利要求1所述的一种流域降雨相似性搜索方法，其特征在于：步骤3)中，对于站点记录的降雨数据，降雨事件划分可以通过限定阈值Th_p、Th_time、Th_int、Th_amount、Th_per来实现，其中阈值Th_p为降雨数值超过0的站点个数占总站点个数的百分比，阈值Th_time为降雨数值超过0的站点的平均降雨时长，阈值Th_int为时间段，阈值Th_amount为降雨数值超过0的站点的平均降雨量，阈值Th_per为降雨数值超过0的站点个数占总站点数的百分比；具体流程为：当发生降雨的站点数量百分比超过阈值Th_p时，认为降雨时间开始，降雨开始后，当超过时段Th_int未再有站点发生降雨，则认为降雨事件结束，当占站点总数百分比为Th_per的站点的平均降雨量超过阈值Th_amount时，则保留此次降雨事件，否则作为微量降雨事件筛除，将平均时长低于阈值Th_time的降雨事件作为短时降雨筛除，从而筛选掉短时微量降雨。

4.根据权利要求1所述的一种流域降雨相似性搜索方法，其特征在于：步骤4)中，对于雷达、卫星遥感和数值预报降雨数据，其原始数据即为矩阵序列；对于雨量站点采集的数据，需要通过转换矩阵序列{Matrix_k,k为流域内站点编号}建立站点数据与栅格数据的转换关系，将站点观测时间序列转化为矩阵序列。

5.根据权利要求4所述的一种流域降雨相似性搜索方法，其特征在于：对于利用权重矩阵序列将站点观测时间序列转化为矩阵序列，具体流程如下：

建立权重矩阵序列{Matrix_k,k为流域内站点编号}，其中的每个矩阵大小均与步骤1)中的栅格矩阵相同， 序列中的第k个权重矩阵记录站点k的数据占各栅格数据的比重， 以各栅格作为目标点，计算站点权重：首先在不考虑方向的情况下，选择流域内的站点利用反向距离权重法计算每个站点对目标点的距离权重：

w_k代表站点k的降雨量占目标点降雨量的权重，其中，k为站点索引，n为站点个数；d_k为站点至网格中心点的距离；

其中，(x，y)为网格中心点坐标；(x_k，y_k)为站点坐标，然后根据流域内站点的方向对权重进行修正，修正系数为：

其中为l为除本站之外其余站点的索引，θ_j为站点l与站点k与目标点连线的夹角，修正后的总距离方向权重为：

W_k＝w_k(1+α_k)

即为站点k对目标点的距离权重，依次计算各站点对每个栅格目标点的距离方向权重保存于对应的权重矩阵中，从而获得权重矩阵序列{Matrix_k}；

获得权重矩阵序列后，利用流域内站点的降水量进行加权平均，得到目标网格的降水量：

其中，n为站点个数，k为站点索引，P_grid为目标网格的降水量，P_k，obs为站点k的观测降水量；针对每个场次降雨，根据权重矩阵序列中保存的权重值利用上述方法将站点观测时间序列转化为栅格点时间序列，进而成为矩阵序列。

6.根据权利要求1所述的一种流域降雨相似性搜索方法，其特征在于：其特征在于：步骤5)中5-3步DTW距离的计算方法为：对时间序列X＝{x₁，x₂，...，x_i，...，x_m}和Y＝{y₁，y₂，...，y_i，...，y_n}，其中m、n分别代表两个时间序列的长度，寻找一个扭曲路径W来表示时间序列X与Y间的映射关系W＝{w₁，w₂，...，w_k，...，w_K}，max(n，m)≤K≤n+m-1，W的第k个元素记为w_k＝(i，j)，表示时间序列X的第i个元素与时间序列Y的第j个元素的对应关系，点(i，j)的累积距离计算公式为：γ(i，j)＝d(x_i，y_j)+min{γ(i-1，j-1)，γ(i-1，j)，γ(i，j-1)}，给定初始条件γ(1，1)＝d(x₁，y₁)，迭代计算得到累积距离矩阵，

即为时间序列X与Y的DTW距离。

Images