CN113191561B - 一种基于高斯混合模型的径流随机模拟方法及系统 - Google Patents

一种基于高斯混合模型的径流随机模拟方法及系统 Download PDF

Info

Publication number
CN113191561B
CN113191561B CN202110510281.5A CN202110510281A CN113191561B CN 113191561 B CN113191561 B CN 113191561B CN 202110510281 A CN202110510281 A CN 202110510281A CN 113191561 B CN113191561 B CN 113191561B
Authority
CN
China
Prior art keywords
flow
period
gaussian mixture
mixture model
distribution
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110510281.5A
Other languages
English (en)
Other versions
CN113191561A (zh
Inventor
周建中
贾本军
徐茂林
方威
陈潇
杨鑫
戴领
仇红亚
田梦琦
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Huazhong University of Science and Technology
Original Assignee
Huazhong University of Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Huazhong University of Science and Technology filed Critical Huazhong University of Science and Technology
Priority to CN202110510281.5A priority Critical patent/CN113191561B/zh
Publication of CN113191561A publication Critical patent/CN113191561A/zh
Application granted granted Critical
Publication of CN113191561B publication Critical patent/CN113191561B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
    • G06Q10/067—Enterprise or organisation modelling

Landscapes

  • Business, Economics & Management (AREA)
  • Engineering & Computer Science (AREA)
  • Human Resources & Organizations (AREA)
  • Strategic Management (AREA)
  • Economics (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Game Theory and Decision Science (AREA)
  • Development Economics (AREA)
  • Marketing (AREA)
  • Operations Research (AREA)
  • Quality & Reliability (AREA)
  • Tourism & Hospitality (AREA)
  • Physics & Mathematics (AREA)
  • General Business, Economics & Management (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Educational Administration (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明提供了一种基于高斯混合模型的径流随机模拟方法及系统，属于随机水文领域，具体为：根据已知t‑1时段流量条件下t时段流量的条件分布以及已知t‑1时段流量条件下t时段流量的条件分布值，计算t时段流量的模拟值；其中，已知t‑1时段流量条件下t时段流量的条件分布的获取方法为：基于历年逐日实测流量序列数据，采用算术平均方法获取某一时间尺度的径流序列数据；基于径流序列数据，构建高斯混合模型；根据高斯混合模型，推求已知t‑1时段流量条件下t时段流量的条件分布。已知t‑1时段流量条件下t时段流量的条件分布值的获取方法为：通过对0～1均匀分布进行随机抽样获得。本发明提高了径流随机模拟的准确性和可靠性。

Description

一种基于高斯混合模型的径流随机模拟方法及系统
技术领域
本发明属于随机水文领域，更具体地，涉及一种基于高斯混合模型的径流随机模拟方法及系统。
背景技术
实测水文径流序列在水资源管理与优化配置中起着至关重要的作用，是确定水资源开发方案、制定水库调度规则、开展流域水文预报、布设流域水文站网、评估水资源系统风险等水资源管理工作的基本依据。由于实测径流序列较短，通常仅涵盖几十年甚至十几年，难以有效满足水资源管理与优化的实际需要。因此，提供径流随机模拟方法，模拟生成长序列径流数据，以延长径流数据长度非常必要。
目前，已有的随机模拟方法主要包括自回归方法和基于Copula函数的径流随机模拟方法。自回归方法简单易用、概念清晰和结构简单，在径流随机模拟中应用较为广泛。然而，实践表明水文径流过程具有较强的偏态特性和非线性特性，基于径流过程线性相关与正态分布假设的自回归方法已难以有效刻画径流过程的滞时相关关系和准确反映径流过程的统计特性。尽管基于Copula函数的径流随机模拟方法一定程度上克服了自回归方法的不足，但依然存在需要单独构建边缘分布函数、滞时相关关系刻画不准确、径流序列统计特征反映不足、模型参数求解困难等问题，无法有效模拟多时间尺度径流序列，同时在反映实际径流序列的均值、方差、偏度、峰度、Pearson相关系数、Kendall相关系数和Spearman相关系数等统计特征方面也存在不足。因此，有必要提出一种更为准确、可靠、通用的径流随机模拟新方法，进一步提高径流随机模拟的精度。
发明内容
针对现有技术的缺陷，本发明的目的在于提供一种基于高斯混合模型的径流随机模拟方法及系统，旨在解决现有径流随机模拟方法模拟精度不高的问题。
为实现上述目的，一方面，本发明提供了一种基于高斯混合模型的径流随机模拟方法为：根据已知t-1时段流量条件下t时段流量的条件分布以及已知t-1时段流量条件下t时段流量的条件分布值，计算t时段流量的模拟值；
其中，已知t-1时段流量条件下t时段流量的条件分布的获取方法，包括以下步骤：
基于历年实测逐日流量序列数据，根据对时间尺度的实际需求，采用算术平均方法获取目标时间尺度的径流序列数据；
基于目标时间尺度的径流序列数据，构建t-1时段和t时段流量的二维联合分布；
根据t-1时段和t时段流量的二维联合分布，推求已知t-1时段流量条件下t时段流量的条件分布。
已知t-1时段流量条件下t时段流量的条件分布值的获取方法为：
通过对0～1均匀分布进行随机抽样获取已知t-1时段流量条件下t时段流量的条件分布值。
优选地，目标时间尺度的径流序列的计算公式为：
Figure BDA0003060049490000021
其中，
Figure BDA0003060049490000022
为目标时间尺度的径流序列的第i年第j时段的平均流量值；
Figure BDA0003060049490000023
为第i年第j时段的天数；
Figure BDA0003060049490000024
为第i年第j时段第k日的日平均流量值。
优选地，t-1时段和t时段流量的二维联合分布的密度值公式为：
Figure BDA0003060049490000025
其中，p为多维联合分布的密度值；x＝[Rt-1,Rt]T为某一年t-1和t时段平均流量构成的二维向量；Rt-1和Rt分别为某一年t-1和t时段的平均流量；N(x|μk,Σk)为二维向量x的二维联合正态分布密度函数，是高斯混合模型的第k个分量；θ＝{π,μ,Σ}为高斯混合模型的待优化参数集，π＝{π1,π2,…,πK}，μ＝{μ1,μ2,…,μK}，Σ＝{Σ1,Σ2,…,ΣK}；πk为第k个分量所占比重，μk为第k个分量的均值向量；Σk为第k个分量的协方差矩阵，K为高斯混合模型的分量个数，是一个超参数；
优选地，构建t-1时段和t时段流量的二维联合分布的方法，包括以下步骤：
S1：基于目标时间尺度的径流序列数据，构建参数集未确定的t-1时段和t时段流量的二维联合分布；其中，参数集未确定的t-1时段和t时段流量的二维联合分布为高斯混合模型；
S2：设置高斯混合模型的超参数K为若干不同的固定值，采用K均值聚类方法确定不同超参数K下高斯混合模型的参数集的初始值；
S3：以各参数集的初始值为优化起点，采用期望最大化方法优化各参数集；
S4：基于AIC信息准则，确定超参数K的最优取值，并获取最优K值对应的优化后的参数集。
优选地，获取参数集的初始值的方法，包括以下步骤：
S2.1：在高斯混合模型的超参数K取某一固定值的基础上，从样本数据中随机选取K个样本点生成初始聚类中心集合；其中，样本数据为
Figure BDA0003060049490000031
N为总样本数；
S2.2：计算每个数据点与聚类中心集中每个聚类中心的欧氏距离，获取距离矩阵；
S2.3：根据距离矩阵，将非中心数据点归入与其距离最近的聚类中心所在的聚类簇，聚类完成后形成聚类簇集合；
S2.4：将聚类簇集合中各聚类簇数据点的均值作为更新后的聚类中心集合；
S2.5：重复S2.2至S2.4，直至迭代次数达到上限或聚类结果收敛，停止重复，执行S2.6；
S2.6：基于S2.5获取的聚类中心集合，计算收敛的聚类簇集合；
S2.7：基于收敛的聚类簇集合，确定高斯混合模型的初始参数。
优选地，优化各参数集的方法，包括以下步骤：
S3.1：根据当前高斯混合模型的参数集的取值情况，计算各个分量对各样本的响应度；
S3.2：基于各个分量对各样本的响应度，更新高斯混合模型的参数集；
S3.3：重复S3.1和S3.2，直至收敛，确定高斯混合模型的最优参数集。
优选地，AIC信息准则的计算公式为：
Figure BDA0003060049490000041
其中，p为分量个数为K的高斯混合模型的参数总个数；使AIC取值最小的K值为高斯混合模型的最优分量个数。
优选地，已知t-1时段流量条件下t时段流量的条件分布p(Rt|Rt-1)的计算公式为：
Figure BDA0003060049490000042
其中，条件均值
Figure BDA0003060049490000043
条件协方差
Figure BDA0003060049490000044
正态分布
Figure BDA0003060049490000045
和正态分布
Figure BDA0003060049490000046
分别为：
Figure BDA0003060049490000047
Figure BDA0003060049490000048
Figure BDA0003060049490000049
Figure BDA00030600494900000410
优选地，t时段流量的模拟值为：
Figure BDA0003060049490000051
其中，
Figure BDA0003060049490000052
是已知t-1时段流量条件下t时段平均流量的条件分布值，通过对0～1均匀分布进行随机抽样获得；F-1()是已知t-1时段流量条件下t时段流量的累积分布函数F()的反函数，F()为：
Figure BDA0003060049490000053
另一方面，本发明提供了一种基于高斯混合模型的径流随机模拟系统，包括：模拟值计算模块、径流序列数据获取模块、高斯混合模型构建模块和条件分布计算模块；
模拟值计算模块用于根据已知t-1时段流量条件下t时段流量的条件分布以及已知t-1时段流量条件下t时段流量的条件分布值，计算t时段流量的模拟值；
径流序列数据获取模块用于基于历年逐日实测流量序列数据，根据对时间尺度的实际需求，采用算术平均方法获取目标时间尺度的径流序列数据；
高斯混合模型构建模块用于基于所述目标时间尺度的径流序列数据，构建t-1时段和t时段流量的二维联合分布；
条件分布计算模块用于根据所述t-1时段和t时段流量的二维联合分布，推求已知t-1时段流量条件下t时段流量的条件分布；
其中，已知t-1时段流量条件下t时段流量的条件分布值的获取方法为：
通过对0～1均匀分布进行随机抽样获取已知t-1时段流量条件下t时段流量的条件分布值。
优选地，高斯混合模型构建模块包括二维联合分布构建单元、参数集初始值获取单元、参数优化单元和最优超参数获取单元；
二维联合分布构建单元用于基于目标时间尺度的径流序列数据，构建参数集未确定的t-1时段和t时段流量的二维联合分布；其中，t-1时段和t时段流量的二维联合分布为高斯混合模型；
参数集初始值获取单元用于设置所述高斯混合模型的超参数K为若干不同的固定值，采用K均值聚类方法确定不同超参数K下高斯混合模型的参数集的初始值；
参数优化单元用于以各参数集的初始值为优化起点，采用期望最大化方法优化各参数集；
最优超参数获取单元用于基于AIC信息准则，确定超参数K的最优取值，并获取最优K值对应的优化后的参数集。
优选地，参数集初始值获取单元包括聚类中心生成器、距离矩阵计算器、聚类簇集合划分器、第一迭代器和初始参数计算器；
聚类中心生成器用于在高斯混合模型的超参数K取某一固定值的基础上，从样本数据中随机选取K个样本点生成初始聚类中心集合；并将聚类簇集合中各聚类簇数据点的均值作为更新后的聚类中心集合；
距离矩阵计算器用于计算每个数据点与聚类中心集中每个聚类中心的欧氏距离，获取距离矩阵；
聚类簇集合划分器用于根据距离矩阵，将非中心数据点归入与其距离最近的聚类中心所在的聚类簇，聚类完成后形成聚类簇集合；
第一迭代器用于判断迭代次数是否达到上限或聚类结果是否收敛；若迭代次数达到上限或聚类结果收敛，则驱动聚类簇集合划分器运行；否则驱动距离矩阵计算器运行；
初始参数计算器用于基于收敛的聚类簇集合，确定高斯混合模型的初始参数。
优选地，参数优化单元包括响应度计算器、参数集更新器和第二迭代器；
响应度计算器用于根据高斯混合模型的参数集的当前取值，计算各个分量对各样本的响应度；
参数集更新器用于基于各个分量对各样本的响应度，更新高斯混合模型的参数集；
第二迭代器用于判断参数集是否收敛，若收敛，则输出高斯混合模型的最优参数集；否则，驱动响应度计算器持续工作。
优选地，目标时间尺度的径流序列的计算公式为：
Figure BDA0003060049490000071
其中，
Figure BDA0003060049490000072
为目标时间尺度的径流序列的第i年第j时段的平均流量值；
Figure BDA0003060049490000076
为第i年第j时段的天数；
Figure BDA0003060049490000073
为第i年第j时段第k日的日平均流量值。
优选地，t-1时段和t时段流量的二维联合分布的密度值公式为：
Figure BDA0003060049490000074
其中，p为多维联合分布的密度值；x＝[Rt-1,Rt]T为某一年t-1和t时段平均流量构成的二维向量；Rt-1和Rt分别为某一年t-1和t时段的平均流量；N(x|μk,Σk)为二维向量x的二维联合正态分布密度函数，是高斯混合模型的第k个分量；θ＝{π,μ,Σ}为高斯混合模型的待优化参数集，π＝{π1,π2,…,πK}，μ＝{μ1,μ2,…,μK}，Σ＝{Σ1,Σ2,…,ΣK}；πk为第k个分量所占比重，μk为第k个分量的均值向量；Σk为第k个分量的协方差矩阵，K为高斯混合模型的分量个数，是一个超参数。
优选地，AIC信息准则的计算公式为：
Figure BDA0003060049490000075
其中，p为分量个数为K的高斯混合模型的参数总个数；使AIC取值最小的K值为高斯混合模型的最优分量个数。
优选地，已知t-1时段流量条件下t时段流量的条件分布p(Rt|Rt-1)的计算公式为：
Figure BDA0003060049490000081
其中，条件均值
Figure BDA0003060049490000082
条件协方差
Figure BDA0003060049490000083
正态分布
Figure BDA0003060049490000084
和正态分布
Figure BDA0003060049490000085
分别为：
Figure BDA0003060049490000086
Figure BDA0003060049490000087
Figure BDA0003060049490000088
Figure BDA0003060049490000089
优选地，t时段流量的模拟值为：
Figure BDA00030600494900000810
其中，
Figure BDA00030600494900000811
是已知t-1时段流量条件下t时段平均流量的条件分布值，通过对0～1均匀分布进行随机抽样获得；F-1()是已知t-1时段流量条件下t时段流量的累积分布函数F()的反函数，F()为：
Figure BDA00030600494900000812
总体而言，通过本发明所构思的以上技术方案与现有技术相比，具有以下有益效果：
本发明采用的高斯混合模型能够在保证模型易于构建与优化的条件下，以更多的参数描述t-1时段和t时段流量之间的联合分布规律，能够很好地模拟实际径流序列的均值、方差、偏度、峰度、Pearson相关系数、Kendall相关系数和Spearman相关系数等统计特征，且模拟精度较现有径流随机模拟方法相比更高。此外，K均值聚类算法和EM优化算法以及AIC信息准则在模型参数优化中的综合应用，进一步增强了本发明径流随机模拟方法的准确性和可靠性。
本发明所采用的高斯混合模型不需要依赖时段流量边缘分布函数的构建，能够直接描述t-1时段和t时段流量之间的联合分布规律，有效简化了模型构建及参数求解工作。
本发明在径流随机模拟应用中，无需拟合时段流量边缘分布函数，不存在哪一类分布曲线更适用于描述时段流量边缘分布规律的问题，增强了径流随机模拟方法的通用性。
附图说明
图1是本发明实施例提供的基于高斯混合模型的径流随机模拟方法流程图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
如图1所示，一方面，本发明提供了一种基于高斯混合模型的径流随机模拟方法，包括以下步骤：
步骤1：基于历年逐日实测流量序列数据，根据对时间尺度的实际需求，采用算术平均方法获得某一种时间尺度的径流序列数据；
为充分揭示本发明并论证其优越性，本实施例对日流量序列、旬流量序列和月流量序列三种时间尺度的径流序列进行随机模拟。日流量序列是指多个日流量组成的流量序列，旬和月流量序列同理。由于实施不同时间尺度径流序列随机模拟的步骤是相似的，本发明采用“时段”通用地代表日、旬和月等不同时间尺度；
步骤2：基于某一时间尺度的径流序列数据，构建t-1时段和t时段流量的高斯混合模型，即两个时段流量的二维联合分布；
步骤3：根据t-1时段和t时段流量的二维联合分布，推求已知t-1时段流量条件下t时段流量的条件分布；
步骤4：根据已知t-1时段流量条件下t时段流量的条件分布和已知t-1时段流量条件下t时段流量的条件分布值，计算t时段流量的模拟值；
其中，已知t-1时段流量条件下t时段流量的条件分布值的获取方法为：
通过对0～1均匀分布进行随机抽样获取已知t-1时段流量条件下t时段流量的条件分布值。
具体地，步骤1中某一种时间尺度的径流序列的计算公式如下：
Figure BDA0003060049490000101
其中，
Figure BDA0003060049490000102
为某一时间尺度的径流序列的第i年第j时段的平均流量值；
Figure BDA0003060049490000103
为第i年第j时段的天数；
Figure BDA0003060049490000104
为第i年第j时段第k日的日平均流量值。
具体地，步骤2中，高斯混合模型的数学公式如下：
Figure BDA0003060049490000105
其中，p为高斯混合模型的输出，即多维联合分布的密度值；x＝[Rt-1,Rt]T为某一年t-1和t时段平均流量构成的二维向量；Rt-1和Rt分别为某一年t-1和t时段的平均流量；N(x|μk,Σk)为二维向量x的二维联合正态分布密度函数，是高斯混合模型的第k个分量；θ＝{π,μ,Σ}为高斯混合模型的待优化参数集，π＝{π1,π2,…,πK}，μ＝{μ1,μ2,…,μK}，Σ＝{Σ1,Σ2,…,ΣK}；πk为第k个分量所占比重，μk为第k个分量的均值向量；Σk为第k个分量的协方差矩阵，K为高斯混合模型的分量个数，是一个超参数；
具体地，步骤2具体包括以下步骤：
步骤2.1建立t-1和t时段流量的高斯混合模型，即相邻时段流量的二维联合分布；
对于日尺度径流序列随机模拟，应有365个二维联合分布；对于旬尺度径流序列随机模拟，应有36个二维联合分布；对于月尺度径流序列随机模拟，应用12个二维联合分布；在此步骤中，高斯混合模型的参数集θ是未知的；
步骤2.2设置高斯混合模型的超参数K为某一个固定值后，采用K均值聚类方法确定高斯混合模型的参数集θ的初始值；
步骤2.3基于步骤2.2获取的初始参数集θ，采用期望最大化方法优化高斯混合模型的参数集θ；
步骤2.4对超参数K设置不同取值，重复步骤2.2至步骤2.3，直至达到预设K的个数；
步骤2.5基于AIC信息准则(Akaike Information Criterion，AIC)确定超参数K的最优取值，同时获取最优K值对应的高斯混合模型的最优参数集θ；
步骤2.2具体包括以下步骤：
步骤2.2.1在高斯混合模型的超参数K取某一固定值的前提下，从
Figure BDA0003060049490000111
样本数据中随机选取K个样本点生成初始聚类中心集合
Figure BDA0003060049490000112
其中，N为总样本数；
步骤2.2.2计算每个数据点xi与聚类中心集Ω中每个聚类中心
Figure BDA0003060049490000113
的欧氏距离dk,i，获取距离矩阵D；
Figure BDA0003060049490000114
D＝{dk,i},(k＝1,2,…,K,i＝1,2,…,N)
步骤2.2.3根据距离矩阵D，将非中心数据点归入与其距离最近的聚类中心所在的聚类簇Ck，聚类完成后形成聚类簇集合C＝{Ck}(k＝1,2,…,K)；
步骤2.2.4根据聚类簇集合C＝{Ck}(k＝1,2,…,K)，重新确定聚类中心集合Ω；
Figure BDA0003060049490000115
Figure BDA0003060049490000116
其中，Nk为聚类簇Ck中样本点个数；
步骤2.2.5重复步骤2.2.2至步骤2.2.4，直至迭代次数达到上限或聚类结果收敛；
步骤2.2.6执行步骤2.2.3，获取已收敛的聚类簇集合C＝{Ck}(k＝1,2,…,K)；
步骤2.2.7基于已收敛的聚类簇集合，确定高斯混合模型第k个分量的参数的初始取值，分别是分量占比πk、均值向量
Figure BDA0003060049490000121
和协方差矩阵
Figure BDA0003060049490000122
Figure BDA0003060049490000123
Figure BDA0003060049490000124
Figure BDA0003060049490000125
具体地，步骤2.3具体包括以下步骤：
基于高斯混合模型的参数集的初始取值，采用EM方法开始迭代优化高斯混合模型的参数集θ；具体如下：
步骤2.3.1根据当前高斯混合模型的参数集的取值情况，计算第k个分量对xi的响应度γi,k(E步)；
Figure BDA0003060049490000126
步骤2.3.2基于第k个分量对xi的响应度γi,k，更新高斯混合模型的参数集θ＝{π,μ,Σ}；
Figure BDA0003060049490000131
Figure BDA0003060049490000132
Figure BDA0003060049490000133
步骤2.3.3重复步骤2.3.1和步骤2.3.2，直至收敛，最终确定高斯混合模型的最优参数集θ；
步骤2.5中，AIC信息准则的计算公式如下：
Figure BDA0003060049490000134
其中，p为分量个数为K的高斯混合模型的参数总个数；使AIC取值最小的K值为高斯混合模型的最优分量个数。最优分量个数对应的通过EM方法求取的参数集θ为高斯混合模型的最终的优化参数集；
具体地，步骤3中，已知t-1时段流量条件下t时段流量的条件分布p(Rt|Rt-1)的计算公式如下：
Figure BDA0003060049490000135
其中，条件均值
Figure BDA0003060049490000136
条件协方差
Figure BDA0003060049490000137
正态分布
Figure BDA0003060049490000138
和正态分布
Figure BDA0003060049490000139
分别为：
Figure BDA00030600494900001310
Figure BDA00030600494900001311
Figure BDA00030600494900001312
Figure BDA0003060049490000141
具体地，步骤4中，计算t时段流量模拟值的计算公式如下：
Figure BDA0003060049490000142
其中，
Figure BDA0003060049490000143
是已知t-1时段流量条件下t时段平均流量的条件分布值，通过对0～1均匀分布进行随机抽样获得；F-1()是已知t-1时段流量条件下t时段流量的累积分布函数F()的反函数，F()如下所示：
Figure BDA0003060049490000144
另一方面，本发明提供了一种基于高斯混合模型的径流随机模拟系统，包括：模拟值计算模块、径流序列数据获取模块、高斯混合模型构建模块和条件分布计算模块；
模拟值计算模块用于根据已知t-1时段流量条件下t时段流量的条件分布以及已知t-1时段流量条件下t时段流量的条件分布值，计算t时段流量的模拟值；
径流序列数据获取模块用于基于历年逐日实测流量序列数据，根据对时间尺度的实际需求，采用算术平均方法获取某一时间尺度的径流序列数据；
高斯混合模型构建模块用于基于所述某一时间尺度的径流序列数据，获取t-1时段和t时段流量的二维联合分布；
条件分布计算模块用于根据所述t-1时段和t时段流量的二维联合分布，推求已知t-1时段流量条件下t时段流量的条件分布。
优选地，高斯混合模型构建模块包括二维联合分布构建单元、参数集初始值获取单元、参数优化单元和最优超参数获取单元；
二维联合分布构建单元用于基于某一时间尺度的径流序列数据，构建参数集未确定的t-1时段和t时段流量的二维联合分布；其中，t-1时段和t时段流量的二维联合分布为高斯混合模型；
参数集初始值获取单元用于设置所述高斯混合模型的超参数K为若干不同的固定值，采用K均值聚类方法确定不同超参数K下高斯混合模型的参数集的初始值；
参数优化单元用于以各参数集的初始值为优化起点，采用期望最大化方法优化各参数集；
最优超参数获取单元用于基于AIC信息准则，确定超参数K的最优取值，并获取最优K值对应的优化后的参数集。
优选地，参数集初始值获取单元包括聚类中心生成器、距离矩阵计算器、聚类簇集合划分器、第一迭代器和初始参数计算器；
聚类中心生成器用于在高斯混合模型的超参数K取某一固定值的基础上，从样本数据中随机选取K个样本点生成初始聚类中心集合；并将聚类簇集合中各聚类簇数据点的均值作为更新后的聚类中心集合；
距离矩阵计算器用于计算每个数据点与聚类中心集中每个聚类中心的欧氏距离，获取距离矩阵；
聚类簇集合划分器用于根据距离矩阵，将非中心数据点归入与其距离最近的聚类中心所在的聚类簇，聚类完成后形成聚类簇集合；
第一迭代器用于判断迭代次数是否达到上限或聚类结果是否收敛；若迭代次数达到上限或聚类结果收敛，则驱动聚类簇集合划分器运行；否则驱动距离矩阵计算器运行；
初始参数计算器用于基于收敛的聚类簇集合，确定高斯混合模型的初始参数。
优选地，参数优化单元包括响应度计算器、参数集更新器和第二迭代器；
响应度计算器用于根据高斯混合模型当前的参数集取值，计算各个分量对各样本的响应度；
参数集更新器用于基于各个分量对各样本的响应度，更新高斯混合模型的参数集；
第二迭代器用于判断参数集是否收敛，若收敛，则输出高斯混合模型的最优参数集；否则，驱动响应度计算器持续工作。
优选地，某一时间尺度的径流序列的计算公式为：
Figure BDA0003060049490000161
其中，
Figure BDA0003060049490000162
为某一时间尺度的径流序列的第i年第j时段的平均流量值；
Figure BDA0003060049490000163
为第i年第j时段的天数；
Figure BDA0003060049490000164
为第i年第j时段第k日的日平均流量值。
优选地，t-1时段和t时段流量的二维联合分布的密度值公式为：
Figure BDA0003060049490000165
其中，p为多维联合分布的密度值；x＝[Rt-1,Rt]T为某一年t-1和t时段平均流量构成的二维向量；Rt-1和Rt分别为某一年t-1和t时段的平均流量；N(x|μk,Σk)为二维向量x的二维联合正态分布密度函数，是高斯混合模型的第k个分量；θ＝{π,μ,Σ}为高斯混合模型的待优化参数集，π＝{π1,π2,…,πK}，μ＝{μ1,μ2,…,μK}，Σ＝{Σ1,Σ2,…,ΣK}；πk为第k个分量所占比重，μk为第k个分量的均值向量；Σk为第k个分量的协方差矩阵，K为高斯混合模型的分量个数，是一个超参数。
优选地，AIC信息准则的计算公式为：
Figure BDA0003060049490000166
其中，p为分量个数为K的高斯混合模型的参数总个数；使AIC取值最小的K值为高斯混合模型的最优分量个数。
优选地，已知t-1时段流量条件下t时段流量的条件分布p(Rt|Rt-1)的计算公式为：
Figure BDA0003060049490000171
其中，条件均值
Figure BDA0003060049490000172
条件协方差
Figure BDA0003060049490000173
正态分布
Figure BDA0003060049490000174
和正态分布
Figure BDA0003060049490000175
分别为：
Figure BDA0003060049490000176
Figure BDA0003060049490000177
Figure BDA0003060049490000178
Figure BDA0003060049490000179
优选地，t时段流量的模拟值为：
Figure BDA00030600494900001710
其中，
Figure BDA00030600494900001711
是已知t-1时段流量条件下t时段平均流量的条件分布值，通过对0～1均匀分布进行随机抽样获得；F-1()是已知t-1时段流量条件下t时段流量的累积分布函数F()的反函数，F()为：
Figure BDA00030600494900001712
实施例：选取长江流域宜昌水文站的逐日实测径流数据进行案例研究；实测径流序列的长度为136年，涵盖1882-2017年，检验基于高斯混合模型的径流随机模拟方法的模拟效果。
分别采用本发明方法、二维Copula方法和三维Copula方法随机模拟生成长度为2000年的日径流、旬径流和月径流序列，在此基础上，针对不同时间尺度的径流序列，分析对比模拟径流序列和实测径流序列的统计特征，检验本发明的模拟效果。
表1
Figure BDA00030600494900001713
Figure BDA0003060049490000181
表2
Figure BDA0003060049490000182
表3
Figure BDA0003060049490000191
表1给出了本发明方法、二维Copula方法和三维Copula方法模拟生成的日径流序列的统计特征与实测日径流序列统计特征之间的相对误差；表2给出了本发明方法、二维Copula方法和三维Copula方法模拟生成的旬径流序列的统计特征与实测旬径流序列统计特征之间的相对误差；表3给出了本发明方法、二维Copula方法和三维Copula方法模拟生成的月径流序列的统计特征与实测月径流序列统计特征之间的相对误差。表1、表2和表3中数字加黑表明该方法的误差最小。
由表1、表2和表3可以看出，无论模拟哪种时间尺度的径流序列，本发明都能够更为全面地准确模拟实测径流序列的统计特性，模拟效果优于其他两种方法。尤其在模拟实测径流序列的高阶统计特性和线性相关特性方面显著优于其他两种方法。
综上所述，本发明与现有技术相比，存在以下优势：
本发明采用的高斯混合模型能够在保证模型易于构建与优化的条件下，以更多的参数描述t-1时段和t时段流量之间的联合分布规律，能够很好地模拟实际径流序列的均值、方差、偏度、峰度、Pearson相关系数、Kendall相关系数和Spearman相关系数等统计特征，且模拟精度较现有径流随机模拟方法相比更高。此外，K均值聚类算法和EM优化算法以及AIC信息准则在模型参数优化中的综合应用，进一步增强了本发明径流随机模拟方法的准确性和可靠性。
本发明所采用的高斯混合模型不需要依赖时段流量边缘分布函数的构建，能够直接描述t-1时段和t时段流量之间的联合分布规律，有效简化了模型构建及参数求解工作。
本发明在径流随机模拟应用中，无需拟合时段流量边缘分布函数，不存在哪一类分布曲线更适用于描述时段流量边缘分布规律的问题，增强了径流随机模拟方法的通用性。
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种基于高斯混合模型的径流随机模拟方法，其特征在于，根据已知t-1时段流量条件下t时段流量的条件分布以及已知t-1时段流量条件下t时段流量的条件分布值，计算t时段流量的模拟值；
其中，已知t-1时段流量条件下t时段流量的条件分布的获取方法，包括以下步骤：
基于历年逐日实测流量序列数据，根据对时间尺度的实际需求，采用算术平均方法获取目标时间尺度的径流序列数据；
基于所述目标时间尺度的径流序列数据，构建t-1时段和t时段流量的二维联合分布；
根据所述t-1时段和t时段流量的二维联合分布，推求已知t-1时段流量条件下t时段流量的条件分布；
已知t-1时段流量条件下t时段流量的条件分布值的获取方法为：
通过对0～1均匀分布进行随机抽样获取已知t-1时段流量条件下t时段流量的条件分布值；
所述目标时间尺度的径流序列的计算公式为：
Figure FDA0003675304250000011
其中，
Figure FDA0003675304250000012
为目标时间尺度的径流序列的第i年第j时段的平均流量值；
Figure FDA0003675304250000013
为第i年第j时段的天数；
Figure FDA0003675304250000014
为第i年第j时段第k日的日平均流量值；
所述t-1时段和t时段流量的二维联合分布的密度值公式为：
Figure FDA0003675304250000021
其中，p为多维联合分布的密度值；x＝[Rt-1,Rt]T为某一年t-1和t时段平均流量构成的二维向量；Rt-1和Rt分别为某一年t-1和t时段的平均流量；N(x|μk,Σk)为二维向量x的二维联合正态分布密度函数，是高斯混合模型的第k个分量；θ＝{π,μ,Σ}为高斯混合模型的待优化参数集，π＝{π1,π2,…,πK}，μ＝{μ1,μ2,…,μK}，Σ＝{Σ1,Σ2,…,ΣK}；πk为第k个分量所占比重，μk为第k个分量的均值向量；Σk为第k个分量的协方差矩阵，K为高斯混合模型的分量个数，是一个超参数；
构建t-1时段和t时段流量的二维联合分布的方法，包括以下步骤：
S1：基于目标时间尺度的径流序列数据，构建参数集未确定的t-1时段和t时段流量的二维联合分布；其中，t-1时段和t时段流量的二维联合分布为高斯混合模型；
S2：设置所述高斯混合模型的超参数K为若干不同的固定值，采用K均值聚类方法确定不同超参数K下高斯混合模型的参数集的初始值；
S3：以各参数集的初始值为优化起点，采用期望最大化方法优化各参数集；
S4：基于AIC信息准则，确定超参数K的最优取值，并获取最优K值对应的优化后的参数集。
2.根据权利要求1所述的径流随机模拟方法，其特征在于，获取参数集的初始值的方法，包括以下步骤：
S2.1：在高斯混合模型的超参数K取某一固定值的基础上，从样本数据中随机选取K个样本点生成初始聚类中心集合；其中，样本数据为
Figure FDA0003675304250000022
N为总样本数；
S2.2：计算每个数据点与聚类中心集中每个聚类中心的欧氏距离，获取距离矩阵；
S2.3：根据距离矩阵，将非中心数据点归入与其距离最近的聚类中心所在的聚类簇，聚类完成后形成聚类簇集合；
S2.4：将聚类簇集合中各聚类簇数据点的均值作为更新后的聚类中心集合；
S2.5：重复S2.2至S2.4，直至迭代次数达到上限或聚类结果收敛，停止重复，执行S2.6；
S2.6：基于S2.5获取的聚类中心集合，计算收敛的聚类簇集合；
S2.7：基于收敛的聚类簇集合，确定高斯混合模型的初始参数。
3.根据权利要求1所述的径流随机模拟方法，其特征在于，优化各参数集的方法，包括以下步骤：
S3.1：根据高斯混合模型的参数集的当前取值，计算各个分量对各样本的响应度；
S3.2：基于各个分量对各样本的响应度，更新高斯混合模型的参数集；
S3.3：重复S3.1和S3.2，直至参数集收敛，确定高斯混合模型的最优参数集。
4.一种基于高斯混合模型的径流随机模拟系统，其特征在于，包括模拟值计算模块、径流序列数据获取模块、高斯混合模型构建模块和条件分布计算模块；
模拟值计算模块用于根据已知t-1时段流量条件下t时段流量的条件分布以及已知t-1时段流量条件下t时段流量的条件分布值，计算t时段流量的模拟值；
径流序列数据获取模块用于基于历年逐日实测流量序列数据，根据对时间尺度的实际需求，采用算术平均方法获取目标时间尺度的径流序列数据；
高斯混合模型构建模块用于基于所述目标时间尺度的径流序列数据，构建t-1时段和t时段流量的二维联合分布；
条件分布计算模块用于根据所述t-1时段和t时段流量的二维联合分布，推求已知t-1时段流量条件下t时段流量的条件分布；
其中，已知t-1时段流量条件下t时段流量的条件分布值的获取方法为：
通过对0～1均匀分布进行随机抽样获取已知t-1时段流量条件下t时段流量的条件分布值；
所述高斯混合模型构建模块包括二维联合分布构建单元、参数集初始值获取单元、参数优化单元和最优超参数获取单元；
所述二维联合分布构建单元用于基于目标时间尺度的径流序列数据，构建参数集未确定的t-1时段和t时段流量的二维联合分布；其中，t-1时段和t时段流量的二维联合分布为高斯混合模型；
所述参数集初始值获取单元用于设置所述高斯混合模型的超参数K为若干不同的固定值，采用K均值聚类方法确定不同超参数K下高斯混合模型的参数集的初始值；
所述参数优化单元用于以各参数集的初始值为优化起点，采用期望最大化方法优化各参数集；
所述最优超参数获取单元用于基于AIC信息准则，确定超参数K的最优取值，并获取最优K值对应的优化后的参数集。
5.根据权利要求4所述的径流随机模拟系统，其特征在于，所述参数集初始值获取单元包括聚类中心生成器、距离矩阵计算器、聚类簇集合划分器、第一迭代器和初始参数计算器；
所述聚类中心生成器用于在高斯混合模型的超参数K取某一固定值的基础上，从样本数据中随机选取K个样本点生成初始聚类中心集合；并将聚类簇集合中各聚类簇数据点的均值作为更新后的聚类中心集合；
所述距离矩阵计算器用于计算每个数据点与聚类中心集中每个聚类中心的欧氏距离，获取距离矩阵；
所述聚类簇集合划分器用于根据距离矩阵，将非中心数据点归入与其距离最近的聚类中心所在的聚类簇，聚类完成后形成聚类簇集合；
所述第一迭代器用于判断迭代次数是否达到上限或聚类结果是否收敛；若迭代次数达到上限或聚类结果收敛，则驱动聚类簇集合划分器运行；否则驱动距离矩阵计算器运行；
所述初始参数计算器用于基于收敛的聚类簇集合，确定高斯混合模型的初始参数。
6.根据权利要求4所述的径流随机模拟系统，其特征在于，所述参数优化单元包括响应度计算器、参数集更新器和第二迭代器；
所述响应度计算器用于根据高斯混合模型的参数集的当前取值，计算各个分量对各样本的响应度；
所述参数集更新器用于基于各个分量对各样本的响应度，更新高斯混合模型的参数集；
所述第二迭代器用于判断参数集是否收敛，若收敛，则输出高斯混合模型的最优参数集；否则，驱动响应度计算器持续工作。
CN202110510281.5A 2021-05-11 2021-05-11 一种基于高斯混合模型的径流随机模拟方法及系统 Active CN113191561B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110510281.5A CN113191561B (zh) 2021-05-11 2021-05-11 一种基于高斯混合模型的径流随机模拟方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110510281.5A CN113191561B (zh) 2021-05-11 2021-05-11 一种基于高斯混合模型的径流随机模拟方法及系统

Publications (2)

Publication Number Publication Date
CN113191561A CN113191561A (zh) 2021-07-30
CN113191561B true CN113191561B (zh) 2022-07-15

Family

ID=76981027

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110510281.5A Active CN113191561B (zh) 2021-05-11 2021-05-11 一种基于高斯混合模型的径流随机模拟方法及系统

Country Status (1)

Country Link
CN (1) CN113191561B (zh)

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20020129038A1 (en) * 2000-12-18 2002-09-12 Cunningham Scott Woodroofe Gaussian mixture models in a data mining system
US20170293856A1 (en) * 2016-04-07 2017-10-12 Xerox Corporation Clustering high dimensional data using gaussian mixture copula model with lasso based regularization
CN107341318B (zh) * 2017-07-17 2020-07-24 河海大学 一种基于全河的月径流时间位移二维矩阵的模拟方法
CN108898250B (zh) * 2018-06-29 2021-10-29 河海大学 一种基于D藤copula函数的月径流模拟方法
CN109344999A (zh) * 2018-09-07 2019-02-15 华中科技大学 一种径流概率预报方法
CN110334314B (zh) * 2019-06-24 2020-09-18 华中科技大学 一种基于条件降维重构的日径流季节性随机模拟方法

Also Published As

Publication number Publication date
CN113191561A (zh) 2021-07-30

Similar Documents

Publication Publication Date Title
Sun et al. Using Bayesian deep learning to capture uncertainty for residential net load forecasting
Serinaldi Multifractality, imperfect scaling and hydrological properties of rainfall time series simulated by continuous universal multifractal and discrete random cascade models
CN110619432B (zh) 一种基于深度学习的特征提取水文预报的方法
CN105719023A (zh) 一种基于混合高斯分布的风电功率实时预测误差分析方法
CN110197020B (zh) 一种环境变化对水文干旱影响的分析方法
CN107610021A (zh) 环境变量时空分布的综合分析方法
Coles et al. Extreme value modelling of hurricane wind speeds
CN113406558B (zh) 基于线性回归的电表失准检测方法、装置及电子设备
Cheng et al. Observation data compression for variational assimilation of dynamical systems
CN111859249A (zh) 一种基于解析四维集合变分的海洋数值预报方法
CN114357737B (zh) 针对大尺度水文模型时变参数的代理优化率定方法
Aljeddani et al. An extensive mathematical approach for wind speed evaluation using inverse Weibull distribution
CN111311026A (zh) 一种顾及数据特征、模型和校正的径流非线性预测方法
CN116796799A (zh) 无水文资料地区中小流域洪水降雨量阈值模型创建方法
CN113191561B (zh) 一种基于高斯混合模型的径流随机模拟方法及系统
CN112287299A (zh) 河流健康变化定量归因方法、装置及系统
CN108108860A (zh) 一种四步耦合中长期水文预报方法
Jurado et al. Implications of modeling seasonal differences in the extremal dependence of rainfall maxima
Zorzetto et al. Bayesian non-asymptotic extreme value models for environmental data
Xing et al. Nash model parameter uncertainty analysis by AM-MCMC based on BFS and probabilistic flood forecasting
CN115825387B (zh) 一种基于气象要素的月尺度土壤湿度预测方法
Cheng Error covariance specification and localization in data assimilation with industrial application
Lemos et al. Spatially varying temperature trends in a Central California estuary
WO2022217568A1 (zh) 一种耦合伯努利-伽马-高斯分布的日尺度降水预报校正方法
Efendi et al. Fuzzy random auto-regression time series model in enrollment university forecasting

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant