本发明公开了一种基于机器学习的局地降雨雨型分析方法，包括以下步骤：1）数据的收集、处理与保存；2）降雨事件的自动提取；3）生成局地降雨样本集合；4）基于GPU加速的降雨事件聚类分析；5）对生成的聚类树进行分析获得代表性雨型。本发明通过采集流域站点观测数据，自动提取降雨事件，进而采用机器学习方法分析最具代表性的降雨过程作为局地降雨代表性雨型，能够极大的节省人为分析的工作量，避免了主观判断带来的差异，同时分析结果更加具有区域的针对性，可以为山洪临界雨量的分析以及城市内涝数值模拟提供有力支撑。
1.一种基于机器学习的局地降雨雨型分析方法，其特征在于：包括以下步骤：
1)数据的收集、处理与保存：收集待分析流域内的水文、气象站点的降雨数据并进行等时段处理；
2)降雨事件的自动提取：依次读取数据库中各站点的连续降雨时间序列，并将其划分为独立的场次降雨，生成多个场次的降雨时间序列；
3)生成局地降雨样本集合：利用步骤2)中提取的多个场次的降雨时间序列生成样本集合，该样本集合的元素为独立的降雨事件且经过标准化处理，集合中元素的个数与降雨场次相同；根据每个降雨事件持续时间的不同将该样本集合划分为多个子集；步骤3)中的子集划分方法为：根据降雨事件集合中每个事件的持续时间，将总时长划分为若干时间区间，将降雨持续时间位于同一区间内的事件提取生成为一个子集；
4)基于GPU加速的降雨事件聚类分析：基于步骤3)生成的样本集合的各子集进行聚类分析，生成若干聚类树，聚类分析的具体步骤为：4-1.生成初始簇：将子集中的每一个元素作为一个初始簇；4-2.计算距离矩阵：矩阵大小为(N×N)，N为该子集中包含的降雨事件个数，矩阵的元素(i,j)为i簇与j簇的距离度量，表示降雨事件i与降雨事件j的相似度，使用DTW距离作为相似性度量标准，距离越小则相似性越强；采用GPU并行计算对该矩阵的计算进行加速；4-3.基于步骤4-2中的距离矩阵对簇进行合并，找出距离最近的两个簇进行合并，将聚类簇重新编号，并计算新簇与其他各簇的距离，更新距离矩阵；4-4.重复步骤4-3直至所有的聚类簇合并为一个簇，由此生成一棵聚类树；4-5.重复步骤4-2～4-4，使基于样本集合中每个子集均生成一棵对应的聚类树；DTW距离的计算方法为：对时间序列X＝{x1，x2，...，xi，...，xm}和Y＝{y1，y2，...，yi，...，yn}，其中m、n分别代表两个时间序列的长度，寻找一个扭曲路径W来表示时间序列X与Y间的映射关系W＝{w1，w2，...，wk，...，wK}，max(n，m)≤K≤n+m-1，W的第k个元素记为wk＝(i，j)，表示时间序列X的第i个元素与时间序列Y的第j个元素的对应关系，点(i，j)的累积距离计算公式为：γ(i，j)＝d(xi，yj)+min{γ(i-1，j-1)，γ(i-1，j)，γ(i，j-1)}，给定初始条件γ(1，1)＝d(x1，y1)，迭代计算得到累积距离矩阵， 即为时间序列X与Y的DTW距离；
5)对步骤4)中生成的聚类树进行分析：将根节点作为聚类树的第1层，则聚类树的第n层包含n个节点，每个节点为1个簇，包含1个聚类中心，遍历搜索给定层的节点簇，计算各节点簇所包含降雨事件的距离矩阵，矩阵大小为(mi×mi)，i＝1,...,n，n为该层节点的个数，mi为第i个节点中包含的降雨事件个数，矩阵的元素(i,j)为降雨事件i与降雨事件j的DTW距离；计算各节点簇所包含降雨事件的距离矩阵，然后计算其各行的和值，和值最小的行索引所对应的降雨事件为该节点簇的聚类中心，即为该流域的局地降雨代表性雨型。
2.根据权利要求1所述的基于机器学习的局地降雨雨型分析方法，其特征在于：步骤1)中降雨数据的时间跨度覆盖10年或10年以上。
3.根据权利要求1所述的基于机器学习的局地降雨雨型分析方法，其特征在于：步骤2)中场次降雨的划分方法为：设定时间阈值，当降雨过程的间歇时间超过阈值则视为两次降雨过程，不足阈值的则视为一次降雨过程；设定量级阈值，当一次降雨过程的总雨量低于量级阈值，则认为是微量降雨，不纳入考虑范围。
4.根据权利要求1所述的基于机器学习的局地降雨雨型分析方法，其特征在于：步骤3)中的降雨事件标准化处理方法为： 其中n为降雨事件的长度，Pi′为标准化后
5.权利要求1所述的基于机器学习的局地降雨雨型分析方法，其特征在于：步骤4)4-2步采用GPU并行计算对该矩阵的计算进行加速的具体方法为：为每个矩阵元素分配一个线程负责DTW距离的计算，首先为线程块分配线程数，对于二维矩阵运算，可为每个线程块分配线程数(tb，tb)，其中tb2需小于GPU所允许的线程块包含的最大线程数；其次，为线程格分配线程块，对于二维矩阵运算，可为每个线程格分配线程块(bg，bg)，其中 N为降雨事件样本子集中的样本数量，bg需小于GPU所允许的线程格包含的最大线程块数量；完成线程分配后，利用GPU完成矩阵各元素的DTW距离计算并返回至CPU内存，由此完成距离矩阵的计算。
一种基于机器学习的局地降雨雨型分析方法
技术领域
[0001] 本发明属于水利工程技术领域，尤其涉及防洪预报技术领域，具体为一种基于机器学习的局地降雨雨型分析方法。
背景技术
[0002] 近年来我国极端暴雨事件频发，局地暴雨突发性强、时效短，是山洪及城市内涝的主要诱发因素。对于局地暴雨，除雨量、雨强外，雨型作为对暴雨过程的描述，表现了暴雨强度在时间尺度上的分配，是暴雨事件的主要致灾特征之一，即便具有相同的雨量及雨强，不同雨型的暴雨过程，致灾性截然不同。
[0003] 由于山丘区暴雨洪水过程陡涨陡落，难以实时预报，目前主要采用动态临界雨量的方法进行山洪预警。同时，我国城镇化快速发展且大多数城市采用的防洪排涝标准较低，内涝灾害频繁，目前通常采用数值模拟的方式对暴雨内涝进行影响评估。已有研究表明，暴雨雨型对山洪灾害临界雨量的确定以及城市内涝积水的最大范围和最大深度的确定有直接影响。
[0004] 目前在进行山洪临界雨量确定和城市内涝数值模拟时，局地暴雨雨型主要采取同频率分析法计算或采用设计雨型，常用的设计雨型有芝加哥雨型、Huff雨型、Pilgrim雨型、Yen&Chow雨型等。其中同频率方法需要人为干预比较多，由于专家经验不同导致的选样不同、理解不同，结果往往带有主观性。芝加哥雨型、Huff雨型、Yen&Chow雨型等各种设计雨型均为外国学者根据某区域暴雨样本概化设计得到，与实际降雨过程存在一定差距，且目前还没有一种公认的雨型作为设计的依据，由于各地的暴雨产生机制存在差异，选取的设计雨型并不一定具有代表性。
发明内容
[0005] 本发明的目的在于针对此问题，提出一种局地降雨雨型分析方法。
[0006] 一种基于机器学习的局地降雨雨型分析方法，包括以下步骤：
[0007] 1)数据的收集、处理与保存：收集待分析流域内的水文、气象站点的降雨数据并进行等时段处理；
[0008] 2)降雨事件的自动提取：依次读取数据库中各站点的连续降雨时间序列，并将其划分为独立的场次降雨，生成多个场次的降雨时间序列；
[0009] 3)生成局地降雨样本集合：利用步骤2)中提取的多个场次的降雨时间序列生成样本集合，该样本集合的元素为独立的降雨事件且经过标准化处理，集合中元素的个数与降雨场次相同；根据每个降雨事件持续时间的不同将该样本集合划分为多个子集；
[0010] 4)基于GPU加速的降雨事件聚类分析：基于步骤3)生成的样本集合的各子集进行聚类分析，生成若干聚类树，聚类分析的具体步骤为：4-1.生成初始簇：将子集中的每一个元素作为一个初始簇；4-2.计算距离矩阵：矩阵大小为(N×N)，N为该子集中包含的降雨事件个数，矩阵的元素(i,j)为i簇与j簇的距离度量，表示降雨事件i与降雨事件j的相似度，使用DTW距离作为相似性度量标准，距离越小则相似性越强；采用GPU并行计算对该矩阵的计算进行加速；4-3.基于步骤4-2中的距离矩阵对簇进行合并，找出距离最近的两个簇进行合并，将聚类簇重新编号，并计算新簇与其他各簇的距离，更新距离矩阵；4-4.重复步骤4-3直至所有的聚类簇合并为一个簇，由此生成一棵聚类树；4-5.重复步骤4-2～4-4，使基于样本集合中每个子集均生成一棵对应的聚类树；
[0011] 5)对步骤4)中生成的聚类树进行分析：将根节点作为聚类树的第1层，则聚类树的第n层包含n个节点，每个节点为1个簇，包含1个聚类中心，遍历搜索给定层的节点簇，计算各节点簇所包含降雨事件的距离矩阵，矩阵大小为(mi×mi)，i＝1,...,n，n为该层节点的个数，mi为第i个节点中包含的降雨事件个数，矩阵的元素(i,j)为降雨事件i与降雨事件j的DTW距离；计算各节点簇所包含降雨事件的距离矩阵，然后计算其各行的和值，和值最小的行索引所对应的降雨事件为该节点簇的聚类中心，即为该流域的局地降雨代表性雨型。
[0012] 进一步的，步骤1)中降雨数据的时间跨度覆盖10年或10年以上。
[0013] 进一步的，步骤2)中场次降雨的划分方法为：设定时间阈值，当降雨过程的间歇时间超过阈值则视为两次降雨过程，不足阈值的则视为一次降雨过程；设定量级阈值，当一次降雨过程的总雨量低于量级阈值，则认为是微量降雨，不纳入考虑范围。
[0014] 进一步的，步骤3)中的降雨事件标准化处理方法为： 其中n为降雨事件的长度，Pi为标准化后的降雨序列点，Pj为原始降雨序列点，i、j分别为标准化序列以及原始序列的时间索引。
[0015] 进一步的，步骤3)中的子集划分方法为：根据降雨事件集合中每个事件的持续时间，将总时长划分为若干时间区间，将降雨持续时间位于同一区间内的事件提取生成为一个子集。
[0016] 进一步的，DTW距离的计算方法为：对时间序列X＝{x1，x2，...，xi，...，xm}和Y＝{y1，y2，...，yi，...，yn}，其中m、n分别代表两个时间序列的长度，寻找一个扭曲路径W来表示时间序列X与Y间的映射关系W＝{w1，w2，...，wk，...，wK}，max(n，m)≤K≤n+m-1，W的第k个元素记为wk＝(i，j)，表示时间序列X的第i个元素与时间序列Y的第j个元素的对应关系，点(i，j)的累积距离计算公式为：γ(i，j)＝d(xi，yj)+min{γ(i-1，j-1)，γ(i-1，j)，γ(i，j-1)}，给定初始条件γ(1，1)＝d(x1，y1)，迭代计算得到累积距离矩阵， 即为时间序列X与Y的DTW距离。
[0017] 进一步的，步骤4)4-2步采用GPU并行计算对该矩阵的计算进行加速的具体方法为：为每个矩阵元素分配一个线程负责DTW距离的计算，首先为线程块分配线程数，对于二维矩阵运算，可为每个线程块分配线程数(tb，tb)，其中tb2需小于GPU所允许的线程块包含的最大线程数；其次，为线程格分配线程块，对于二维矩阵运算，可为每个线程格分配线程块(bg，bg)，其中 N为降雨事件样本子集中的样本数量，bg需小于GPU所允许的线程格包含的最大线程块数量；完成线程分配后，利用GPU完成矩阵各元素的DTW距离计算并返回至CPU内存，由此完成距离矩阵的计算。
[0018] 本发明的有益效果：
[0019] 通过采集流域(区域)站点观测数据，自动提取降雨事件，进而采用机器学习方法分析最具代表性的降雨过程作为局地降雨代表性雨型，能够极大的节省人为分析的工作量，避免了主观判断带来的差异，与传统设计雨型相比，根据本发明中提出的方法所分析的雨型结果更加具有区域的针对性，可以为山洪临界雨量的分析以及城市内涝数值模拟提供有力支撑。
[0020] 下面结合附图及具体实施方式对本发明作进一步详细说明。
附图说明
[0021] 图1为本发明方法整体流程图；
[0022] 图2降雨数据插值示意图；
[0023] 图3时间序列的动态扭曲路径；
[0024] 图4分的降雨事件整体；
[0025] 图5雨样本示例1；
[0026] 图6雨样本示例2；
[0027] 图7雨样本示例3；
[0028] 图8雨样本示例4；
[0029] 图9雨样本示例5；
[0030] 图10降雨样本示例6；
[0031] 图11第1子集聚类树；
[0032] 图12第2子集聚类树；
[0033] 图13第3子集聚类树；
[0034] 图14第4子集聚类树；
[0035] 图15第5子集聚类树；
[0036] 图16第6子集聚类树；
[0037] 图17第1子集的第1聚类簇；
[0038] 图18第1子集的第2聚类簇；
[0039] 图19第1子集的第3聚类簇；
[0040] 图20第1子集的第4聚类簇；
[0041] 图21第1子集的第5聚类簇；
[0042] 图22第1子集的代表雨型1；
[0043] 图23第1子集的代表雨型2；
[0044] 图24第1子集的代表雨型3；
[0045] 图25第1子集的代表雨型4；
[0046] 图26第1子集的代表雨型5。
具体实施方式
[0047] 实施例1
[0048] 一种基于机器学习的局地降雨雨型分析方法，包括以下步骤：
[0049] 1)数据的收集、处理与保存
[0050] 数据的收集：收集待分析流域(区域)内的水文、气象站点的降雨观测数据，由于聚类分析对数据量需求较大，降雨数据的时间跨度需要覆盖10年或10年以上。
[0051] 数据的处理：将降雨数据处理为等时段时间序列，若原始数据为非等时段数据，则需对数据进行插值处理，优选根据降雨量累积曲线进行插值，如图2所示，首先利用原始序列获得降雨量累积曲线，进而差分获得等时段降雨时间序列{P′1，P′2，P′3，...，P′12}。
[0052] 将处理好的等时段降雨时间序列数据保存至数据库。
[0053] 2)降雨事件的自动提取
[0054] 依次读取数据库中各站点的降雨时间序列，将其划分为独立的场次降雨。以降雨时间序列{P1，P2，P3，...，Pt}及其对应的时间标识序列{T1，T2，T3，...，Tt}为例，划分方法为：设定时间阈值ThT，当降雨过程的间歇时间Tj-Ti超过阈值ThT则视为两次降水过程，不足阈值ThT则视为一次降水过程，从而实现自动连续的降雨场次划分；设定量级阈值ThA，当一次降雨过程的总雨量低于阈值ThA，则认为是微量降雨，不纳入考虑范围。依靠上述方法，遍历计算各站点降雨时间序列得到n个场次降雨序列{Pi1，Pi2，...，Pik}及其时间标识序列{Ti1，Ti2，...，Tik}，其中i＝1，...，N，N为降雨场次个数，k为该场降雨对应的时段个数。
[0055] 3)生成局地降雨样本集合
[0056] 利用步骤2)中提取的N个场次降雨序列{Pi1，Pi2，...，Pik}，i＝1，...，N，N为降雨场次个数，k为该场降雨对应的时段个数，生成包含N个元素样本集合，集合元素为独立的降雨事件。利用公式 对降雨事件进行标准化处理，其中n为降雨事件的长度，P′i为标准化后的降雨序列点，Pj为原始降雨序列点。根据每个降雨事件持续时间的不同划分该样本集合划分为多个子集。
[0057] 4)基于GPU加速的降雨事件聚类分析
[0058] 基于步骤3)中生成的样本集合的每个子集进行聚类分析，生成聚类树。聚类分析的具体步骤为：
[0059] 4-1.生成初始簇，将子集中的每一个元素作为一个初始簇，对于一个具有N个数据对象的子集D＝{x1，x2，...，xN}，设定初始簇集合C＝{C1，C2，...，CN}，其中Cj＝{xj}；
[0060] 4-2.计算第一距离矩阵MatrixF，矩阵大小为(N×N)，N为子集中包含的降雨事件个数，矩阵的元素(i，j)为簇i与簇j的相似度，对于第一距离矩阵，元素(i，j)即为降雨事件i与降雨事件j的相似度。使用DTW距离作为相似性度量标准，距离越小则相似性越强，DTW距离计算方法如下：
[0061] 对时间序列X＝{x1，x2，...，xi，...，xm}和Y＝{y1，y2，...，yi，...，yn}，寻找一个扭曲路径W来表示时间序列X与Y间的映射关系，如图3所示，W＝{w1，w2，...，wk，...，wK}，max(n，m)≤K≤n+m-1，W的第k个元素记为wk＝(i，j)，表示时间序列X的第i个元素与时间序列Y的第j个元素的对应关系。扭曲路径的选取有三个约束条件：扭曲路径始于矩阵的起始元素，结束于对角元素，即w1＝(1，1)，wK＝(m，n)；扭曲路径每一步都是连续的，即对于wk＝(a，b)，wk-1＝(a′，b′)，要求a-a′≤1且b-b′≤1；扭曲路径在时间轴上是单调的，即对于wk＝(a，b)，wk-1＝(a′，b′)，要求a-a′≥0且b-b′≥0。
[0062] 能够满足约束条件的路径有很多条，此处寻找扭曲代价最小的路径，即：
[0063]
[0064] 其中d(wk)为wk代表的两个对应元素间的距离。
[0065] 根据动态规划思想，若点(i，j)在最佳路径上，那么从点(1，1)到点(i，j)的子路径也是局部最优解，即从点(1，1)到点(m，n)的最佳路径可以由起始点(1，1)到终点(m，n)之间的局部最优解递归搜索获得，因而可以方便地找到这个最佳路径。具体步骤为：首先构建一2
个m×n阶矩阵，矩阵元素(i，j)为两个时间序列点xi和点yj之间的距离d(xi，yj)＝(xi-yj) 。
定义点(i，j)的累积距离计算公式：
[0066] γ(i，j)＝d(xi，yj)+min{γ(i-1，j-1)，γ(i-1，j)，γ(i，j-1)}
[0067] 给定初始条件γ(1，1)＝d(x1，y1)，可以迭代计算得到累积距离矩阵。即为时间序列X与Y的DTW距离，从点γ(m，n)出发反向搜索累积距离矩阵即可得到最佳匹配路径。
[0068] 由于第一距离矩阵计算的时间复杂度为 其中N为样本集合中元素个数，同时DTW距离计算的时间复杂度为O(m·n)，其中m、n为降雨事件时间序列的长度，因此第一距离矩阵的计算耗时往往很大，传统方法难以满足大数据量分析的需求。
[0069] 介于各元素的计算互相独立，适用于并行计算，因而采用GPU并行计算对第一矩阵的计算进行加速。具体方法为：为每个矩阵元素分配一个线程负责DTW距离的计算，首先为线程块分配线程数，线程块中的最大线程数根据GPU性能的不同有所区别，对于二维矩阵运算，可为每个线程块分配线程数(tb，tb)，其中tb2需小于GPU所允许的线程块包含的最大线程数；其次，为线程格分配线程块，对于二维矩阵运算，可为每个线程格分配线程块(bg，bg)，其中 且小于GPU所允许的线程格包含的最大线程块数量，N为子集中的样本数量。完成线程分配后，利用GPU完成矩阵各元素的DTW距离计算并返回至CPU内存，由此完成第一距离矩阵的计算。
[0070] 4-3.基于第一距离矩阵对簇进行合并，找出距离最近的两个簇Ci*和Cj*，合并Ci*和Cj*：Ci*＝Ci*∪Cj*，将聚类簇重新编号，删除距离矩阵M(当前距离矩阵)的第j*行和第j*列，并计算新簇与其他各簇的距离，更新距离矩阵。
[0071] 4-4重复上一步骤直至所有聚类簇合并为一个簇，由此生成一棵聚类树。
[0072] 4-5.重复步骤4-2～4-4使基于样本集合中每个子集均生成一棵对应的聚类树。
[0073] 5)聚类中心提取及局地降雨雨型分析
[0074] 对步骤4)中生成的聚类树进行分析，将根节点作为聚类树的第1层，则聚类树的第n层包含n个节点，每个节点为1个簇，包含1个聚类中心，由根节点向下遍历搜索各层节点簇，计算各节点簇所包含降雨事件的距离矩阵MatrixD，矩阵大小为(mi×mi)，i＝1，...，n，n为该层节点个数，i为节点索引，mi为节点中包含的降雨事件个数，矩阵的元素(i，j)为降雨事件i与降雨事件j的DTW距离(计算方法可以参照步骤4)中DTW距离的计算方法)。首先计算距离矩阵MatrixD，然后计算其各行的和值，和值最小的行索引所对应的降雨事件为聚类中心，即为该流域(区域)的局地降雨代表性雨型。
[0075] 本实施例中：对长江流域的某子流域进行局地降雨雨型分析，流域面积572平方公里，流域内共18个雨量站点，站点雨量数据时间跨度为15年。
[0076] 设置阈值ThA＝10mm，阈值ThT＝6h，阈值ThL＝3h，经过降雨事件的自动提取，共提取出3092个降雨事件，降雨事件持续时长为3至120小时，降雨事件如图4所示。
[0077] 根据降雨事件计算累积降雨过程并标准化，将标准化的累积降雨过程时间序列作为样本，并根据降雨时长区间划分子集，按照[3，6)、[6，12)、[12，24)、[24，48)、[48，96)、[96，192)6个区间共划分为6个子集，每个子集中的样本个数分别为：521、874、1051、542、97、7，部分降雨过程及样本如图5～10所示。
[0078] 以每个子集为基础分别进行聚类分析，得到聚类树如图11～16所示：
[0079] 根据已生成的聚类树，选择第5层提取聚类中心，以第一子集即3至5小时降雨事件为例：对于基于第一个子集生成的聚类树，可以获得该流域3至5小时降雨的5种代表雨型，5个聚类簇和聚类中心(代表雨型)分别如图17～26所示。
