CN113705877B - 基于深度学习模型的实时月径流预报方法 - Google Patents
基于深度学习模型的实时月径流预报方法 Download PDFInfo
- Publication number
- CN113705877B CN113705877B CN202110966434.7A CN202110966434A CN113705877B CN 113705877 B CN113705877 B CN 113705877B CN 202110966434 A CN202110966434 A CN 202110966434A CN 113705877 B CN113705877 B CN 113705877B
- Authority
- CN
- China
- Prior art keywords
- month
- forecasting
- runoff
- training
- data
- 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
Links
- PEDCQBHIVMGVHV-UHFFFAOYSA-N Glycerine Chemical compound OCC(O)CO PEDCQBHIVMGVHV-UHFFFAOYSA-N 0.000 title claims abstract description 32
- 238000013136 deep learning model Methods 0.000 title claims abstract description 30
- 238000013277 forecasting method Methods 0.000 title claims abstract description 20
- 238000012549 training Methods 0.000 claims abstract description 81
- 238000000034 method Methods 0.000 claims abstract description 65
- 238000012795 verification Methods 0.000 claims abstract description 14
- 239000013598 vector Substances 0.000 claims abstract description 13
- 238000013527 convolutional neural network Methods 0.000 claims abstract description 10
- 238000012545 processing Methods 0.000 claims abstract description 10
- 238000004458 analytical method Methods 0.000 claims abstract description 8
- 238000010606 normalization Methods 0.000 claims abstract description 7
- 238000013135 deep learning Methods 0.000 claims abstract description 6
- 230000006870 function Effects 0.000 claims description 24
- 239000011159 matrix material Substances 0.000 claims description 22
- 238000011176 pooling Methods 0.000 claims description 22
- 238000004364 calculation method Methods 0.000 claims description 18
- 238000013528 artificial neural network Methods 0.000 claims description 15
- 230000004913 activation Effects 0.000 claims description 10
- 230000008859 change Effects 0.000 claims description 10
- 230000008569 process Effects 0.000 claims description 10
- 238000012937 correction Methods 0.000 claims description 8
- 230000001419 dependent effect Effects 0.000 claims description 6
- 238000001704 evaporation Methods 0.000 claims description 6
- 230000008020 evaporation Effects 0.000 claims description 6
- 238000013507 mapping Methods 0.000 claims description 5
- 238000001556 precipitation Methods 0.000 claims description 5
- 241001502050 Acis Species 0.000 claims description 4
- 230000006835 compression Effects 0.000 claims description 4
- 238000007906 compression Methods 0.000 claims description 4
- 238000002790 cross-validation Methods 0.000 claims description 3
- 239000002689 soil Substances 0.000 claims description 3
- 241001123248 Arma Species 0.000 claims 5
- 238000003064 k means clustering Methods 0.000 abstract description 8
- 238000012216 screening Methods 0.000 abstract description 8
- 238000005192 partition Methods 0.000 abstract description 3
- 210000004027 cell Anatomy 0.000 description 8
- 238000010586 diagram Methods 0.000 description 5
- 238000005457 optimization Methods 0.000 description 5
- 238000004422 calculation algorithm Methods 0.000 description 3
- 238000000611 regression analysis Methods 0.000 description 3
- ORILYTVJVMAKLC-UHFFFAOYSA-N Adamantane Natural products C1C(C2)CC3CC1CC2C3 ORILYTVJVMAKLC-UHFFFAOYSA-N 0.000 description 2
- 238000003491 array Methods 0.000 description 2
- 238000010276 construction Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000010801 machine learning Methods 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 210000002569 neuron Anatomy 0.000 description 2
- 238000005096 rolling process Methods 0.000 description 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 2
- 108010014173 Factor X Proteins 0.000 description 1
- 238000007792 addition Methods 0.000 description 1
- 230000002776 aggregation Effects 0.000 description 1
- 238000004220 aggregation Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 125000004122 cyclic group Chemical group 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 238000007726 management method Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000008239 natural water Substances 0.000 description 1
- 230000000737 periodic effect Effects 0.000 description 1
- 230000002035 prolonged effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000012706 support-vector machine Methods 0.000 description 1
- 238000012731 temporal analysis Methods 0.000 description 1
- 238000000700 time series analysis Methods 0.000 description 1
- 238000010200 validation analysis Methods 0.000 description 1
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
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
- G06F18/214—Generating training patterns; Bootstrap methods, e.g. bagging or boosting
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/23—Clustering techniques
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/02—Neural networks
- G06N3/04—Architecture, e.g. interconnection topology
- G06N3/045—Combinations of networks
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/02—Neural networks
- G06N3/08—Learning methods
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A10/00—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE at coastal zones; at river basins
- Y02A10/40—Controlling or monitoring, e.g. of flood or hurricane; Forecasting, e.g. risk assessment or mapping
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- General Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Life Sciences & Earth Sciences (AREA)
- Artificial Intelligence (AREA)
- Evolutionary Computation (AREA)
- Business, Economics & Management (AREA)
- Software Systems (AREA)
- Biomedical Technology (AREA)
- Computing Systems (AREA)
- General Health & Medical Sciences (AREA)
- Computational Linguistics (AREA)
- Mathematical Physics (AREA)
- Biophysics (AREA)
- Human Resources & Organizations (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Evolutionary Biology (AREA)
- Health & Medical Sciences (AREA)
- Molecular Biology (AREA)
- Economics (AREA)
- Strategic Management (AREA)
- Development Economics (AREA)
- Game Theory and Decision Science (AREA)
- Entrepreneurship & Innovation (AREA)
- Marketing (AREA)
- Operations Research (AREA)
- Quality & Reliability (AREA)
- Tourism & Hospitality (AREA)
- General Business, Economics & Management (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明提供基于深度学习模型的实时月径流预报方法,包括:步骤1.基于历史信息和未来气象信息收集预报因子,并根据流域历史时期月径流的自相关性分析确定前期月径流对预报月影响的最长滞时;步骤2.对训练期的预报因子和月径流数据分别进行归一化处理,再采用基于嵌入式思想的LASSO回归方法自动筛选预报因子;步骤3.采用基于划分思想的K均值聚类方法对训练期样本集进行聚类,将样本分成互不重合的K类;步骤4.计算验证集预报因子向量与K个训练集的聚类中心的距离,找到距离最近的训练集然后以该数据集训练卷积神经网络与门控循环单元网络结合的组合深度学习预报模型;步骤5.采用自回归滑动平均模型对预报残差进行实时校正。
Description
技术领域
本发明属于水文预报技术领域,具体涉及基于深度学习模型的实时月径流预报方法。
技术背景
月径流预报是水文学领域的重要工程技术难题之一,不仅能为解决天然来水与人为用水不协调、指导流域水资源开发管理提供信息支持,还是决策者掌握防洪抗旱工作主动权,有效规避自然灾害的前提和依据之一。月径流过程作为弱相关且高度复杂的非线性动力系统,对预报模型构建要求较高。
通常而言,月径流预报模型可分为过程驱动和数据驱动两类。过程驱动模型又称物理成因分析法,需要借助能够反映流域产汇流特征的水文模型,以未来预报气象信息作为输入获取预报结果。然而,构建过程驱动模型较为复杂且与流域特性关系密切,工程应用性较差。数据驱动模型直接基于历史数据建立预报对象与预报因子的数学关系,进而预报未来水文变量。数据驱动模型通常包括时间序列分析(差分自回归滑动平均模型等)、回归分析(岭回归等)以及机器学习方法(支持向量机、人工神经网络等)。近年来机器学习方法在径流预报领域受到广泛关注,尤其是深度学习算法展现了良好的预报性能。例如,李文武等提出了基于变分模态分解和深度门控网络的径流预测方法(水力发电学报,2020,39(3):34-44);专利CN202010285986.7提出了基于多模型组合的中长期径流集合预报方法;岳兆新等提出了基于改进深度信念网络模型的中长期径流预测方法(水力发电学报,2020,39(10):33-46)。
然而,现有的基于深度学习的月径流预报方法存在三大问题:(1)预报因子来源单一,仅考虑了历史时期的观测数据,未考虑将数值预报产品预报的未来气象信息作为预报因子;(2)预报因子筛选方法单一,当前大部分研究均采用过滤法(如皮尔逊相关系数、互信息系数等)筛选预报因子,预报因子个数基于主观意愿,存在较大的任意性和不确定性;(3)预报模型单一,未能考虑月径流的时间异质性,尤其是对于汛期和非汛期径流,采用统一的模型欠缺考虑。此外,现有深度学习模型输出的月径流预报的精度与实际需求仍存在一定的差距。
发明内容
本发明是为了解决上述问题而进行的,目的在于提供一种基于深度学习模型的实时月径流预报方法,能够保证良好的月径流预报精度。
本发明为了实现上述目的,采用了以下方案:
本发明所提供的基于深度学习模型的实时月径流预报方法,包括以下步骤:
步骤1.基于历史信息和未来气象信息收集预报因子,并根据流域历史时期月径流的自相关性分析确定前期月径流对预报月影响的最长滞时,并以此值统一作为其它预报因子的最长影响滞时;历史信息包括:流域当地气象信息,包含大气环流指数ACIs、海温指数SSTs和其它指数OCIs的全球气候指数,包含地表温度、实际蒸散发、潜在蒸散发、归一化植被指数及土壤湿度的多源遥感数据;包含预报月以前的历史月径流数据的前期径流;未来信息包括数值气象预报的月降水和月气温数据;
步骤2.对训练期的预报因子和月径流数据分别进行归一化处理,再采用基于嵌入式思想的LASSO回归方法筛选预报因子;
步骤3.采用基于划分思想的K均值聚类方法对训练期样本集进行聚类,应用最近邻规则将样本分成互不重合的K类;包括如下子步骤:
步骤3.1采用肘部法则确定聚类数K值;核心指标是误差平方和SSE,式中,X为待归类的点,C为聚类中心点;肘部法则的计算原理是成本函数,成本函数是类别畸变程度之和,每个类的畸变程度等于每个变量点到其类别中心的位置距离平方和即误差平方和,若各类内部的成员彼此间越紧凑则类的畸变程度越小;将聚类数K从1开始依次增加,统计误差平方和与聚类类别数K的关系,找到误差平方和下降速度较快和误差平方和变化变缓的临界点,作为最佳的聚类个数;
步骤3.2假设训练期的长度为M个月,每个月对应N个预报因子;随机从M组数列中选中K个点(每个点都是一个长度为N的行向量),每个点代表每簇的初始聚类中心,完成初始化;
步骤3.3计算剩余各个点到聚类中心的欧式距离,公式为:式中,N为每个点的维数;根据最近邻规则,将其归为离它最近的簇;
步骤3.4重新计算每一个簇的平均值,将其作为新的聚类中心;
步骤3.5更新M组数列的归类结果,即计算每个点到聚类中心的欧氏距离,根据最邻近规则,将其归为离它最近的簇;
步骤3.6重复步骤3.4和步骤3.5,直至M组数列所属的聚类中心不再发生变化或满足设定的迭代次数,则终止全过程;
步骤4.计算验证集预报因子向量与K个训练集的聚类中心的距离,找到距离最近的训练集然后以该数据集训练卷积神经网络与门控循环单元网络结合的组合深度学习预报模型;
步骤5.采用自回归滑动平均模型对预报残差进行实时校正:
步骤5.1将深度学习模型输出的训练期的所有月径流值作如下处理:处理过后的x值作为预报残差自回归滑动平均模型的输入;
步骤5.2模型定阶:利用AIC定阶准则确定ARMA模型的自回归项阶数p和移动平均项阶数q,认定最小AIC值对应的模型是最好的模型,完成模型定阶;
步骤5.3采用矩估计法对定阶的ARMA模型参数进行估计;
步骤5.4根据确定的参数建立ARMA模型;
步骤5.5将预测月之前的月径流数据作为ARMA模型的自变量x,则ARMA模型对应的因变量为y,将y进行如下变化后作为修正后的月径流预报值,
优选地,本发明提供的基于深度学习模型的实时月径流预报方法,还可以具有以下特征:在步骤1中,待选预报因子乘lag_k月滞时组成超高维矩阵;在步骤2中,设训练期的自变量矩阵为预报因子矩阵X=(X1,X2,...,Xj,...,Xn),其中Xj=(x1j,x2j,...,xmj),训练期的因变量矩阵为待预报月径流矩阵Y=(y1,y2,...,ym)T;LASSO回归首先建立预报因子X与待预报月径流Y之间的线性模型Y=α+β1X1+β2X2+...+βnXn;其中,α为常数项,n为待选预报因子乘lag_k月滞时所组成的超高维矩阵中的总因子数,β为各变量系数;LASSO回归设置最小二乘形式的目标函数如下:其中λ为惩罚参数,其所在项为惩罚项;在LASSO回归中通过调整惩罚参数λ的大小,以达到产生稀疏解的目的,即让一些预报因子的系数值等于零从而达到因子筛选的目的;LASSO回归中惩罚参数λ越大,回归系数压缩幅度越大,相反的,λ越小,则回归系数压缩幅度较小,过小的λ可能引入过多的预报因子而无法解决过拟合问题;采用十折交叉验证法,来确定惩罚参数λ的大小。
优选地,本发明提供的基于深度学习模型的实时月径流预报方法,还可以具有以下特征:在步骤2中,是基于分界点将总样本集进行分界,第一次训练先以位于分界点以前的样本作为训练期的数据,预报训练期下一个月的月径流,第二次训练再以第一次训练期的数据和预报数据作为训练期的数据,预报训练期下一个月的月径流;以此类推,按照前述过程完成涵盖总样本集中首月至倒数第二个月时间范围内所有月份数据的训练。
优选地,本发明提供的基于深度学习模型的实时月径流预报方法,还可以具有以下特征:步骤4包括如下子步骤:
步骤4.1将训练期的数据输入卷积神经网络的卷积层的第一层;训练期的数据形式为((X1,Y1),(X2,Y2),...,(Xt,Yt),...,(XM,YM));式中,Xt∈RM×N(1≤t≤M),X为预报因子即模型输入,Y为待预报月径流即模型输出,M为训练期长度,N为预报因子个数;卷积层包括m个一维卷积核,用于将输入数据转化为特征映射,卷积层的计算公式为Mj=fcov(P*Wj+bj),式中,*代表卷积算子,Mj、Wj和bj分别代表第j个卷积核获得的特征映射,第j个卷积核的权重矩阵,和第j个卷积核的偏置,fcov代表神经网络激活函数,推荐使用线性整流函数(ReLU),该函数保留正数元素,同时将负数元素替换为0;
步骤4.2将卷积层输出的特征映射输入池化层,池化层的主要作用为完成下采样和减少参数量,发明人推荐使用最简单也最常用的最大池化层;
步骤4.3一个卷积层加一个池化层构成一组卷积神经网络结构,优选为放置两组卷积神经网络结构;
步骤4.4为了将池化层的输出变成循环神经网络GRU可接受的输入形式,需加入Flatten层进行数据处理,通过Flatten层将池化层的输出转化成一维数列;
步骤4.5放置两层GRU层;每个GRU神经网络单元包含一个重置门和一个更新门,Rt和Zt分别代表重置门和更新门,计算公式分别表示为Rt=σ(WR·[Ht-1,Xt]+bR)和Zt=σ(WZ·[Ht-1,Xt]+bZ),式中,Ht-1代表记忆单元在t-1时刻的输出状态,WR和WZ分别代表重置门和更新门的权重矩阵,bR和bZ分别代表重置门和更新门的偏置向量,Xt代表记忆单元在t时刻的输入,σ代表激活函数sigmoid,而记忆单元在t时刻的输出状态的计算公式为Ht=Ct,记忆单元在时刻t的新状态Ct计算公式又为Ct=(1-Zt)⊙Ct-1+Zt⊙Ct′,当前可能的单元状态Ct′的计算公式又为Ct′=tanh(WCXt+WC(Rt⊙Ht-1)+bC),式中,WC和bC分别代表记忆单元的权重矩阵和偏置向量,⊙代表向量元素的乘积,tanh代表激活函数;
步骤4.6在循环神经网络层后放置两层全连接层,并在两层全连接层中加入Dropout层,以减少过拟合现象;
步骤4.7在此步之前的操作全部都基于归一化的数值,故需要进行反归一化后将数据输出。
优选地,本发明提供的基于深度学习模型的实时月径流预报方法,还可以具有以下特征:在步骤4.3中,两组卷积神经网络结构为:卷积层1+池化层1+卷积层2+池化层2。
优选地,本发明提供的基于深度学习模型的实时月径流预报方法,还可以具有以下特征:在步骤4.5中,激活函数为双曲正切函数。
发明的作用与效果
本发明所提供的基于深度学习模型的实时月径流预报方法,首先基于历史信息和未来气象信息收集预报因子,采用基于嵌入式思想的LASSO回归方法对所有的独立变量同时进行处理、自动筛选预报因子,然后基于预报因子通过样本训练集聚类方法建立组合式深度学习模型,将训练期月径流数据离散化,将总样本训练集分为不同的训练样本子集,根据不同的子样本数据集分别建立组合式深度学习模型,将验证集的数据逐一输入,通过预报因子性质判断该月应归于哪一类训练集,从而驱动通过以该训练集为基础搭建的深度学习模型,进一步通过预报残差自回归滑动平均模型实时校正深度学习模型输出的月径流预报序列,因此能够根据实时更新的上月径流实测值与上月径流预报值之间的关系,及时有效地校正组合式深度学习模型输出的当月径流预报值,可实现预见期长达12个月甚至更长的滚动月径流预报,并保持良好的预报精度。
附图说明
图1为本发明实施例涉及的基于深度学习模型的实时月径流预报方法的流程图;
图2为本发明实施例涉及的月径流自相关性(1)与偏自相关性(2)分析的示意图;
图3为本发明实施例涉及的模型输入与输出数据的示意图;
图4为本发明实施例涉及的LASSO回归中惩罚参数λ值与均方误差(MSE)值关系的示意图;
图5为本发明实施例涉及的在K均值聚类中聚类数K值与误差平方和(SSE)关系的示意图;
图6为本发明实施例涉及的采用逐步回归分析法和LASSO回归法进行预报因子优选的预报结果对比图;
图7为本发明实施例涉及的不采用K均值聚类与采用聚类方法的预报结果对比图;
图8为本发明实施例涉及的在采用LASSO回归法进行预报因子优选、采用K均值聚类法对训练集样本进行划分的基础上、不采用自回归滑动模型进行预报结果校正与进行预报结果校正的预报结果对比图。
具体实施方式
以下结合附图对本发明涉及的基于深度学习模型的实时月径流预报方法进行详细地说明。
<实施例>
如图1所示,本实施例所提供的基于深度学习模型的实时月径流预报方法包括以下步骤:
步骤1.基于历史信息和未来气象信息收集预报因子,并根据流域历史时期月径流的自相关性分析确定前期月径流对预报月影响的最长滞时,并以此值统一作为其它预报因子的最长影响滞时。步骤1进一步包括以下子步骤:
步骤1.1本实施例收集的待选预报因子包括历史信息和未来信息,下述所有待选预报因子均为月平均值,历史信息包括:①流域当地气象信息,包括降水、最高气温、最低气温等;②全球气候指数,包括来自国家气候中心的气候系统指数,含大气环流指数(ACIs)、海温指数(SSTs)和其它指数(OCIs);③多源遥感数据,包括降水量、地表温度、实际蒸散发、潜在蒸散发、归一化植被指数及土壤湿度等;④前期径流,即预报月以前的历史月径流数据;⑤未来信息,包括数值气象预报的月降水和月气温数据等。
本实施例收集的数据序列为从1951年1月到2020年的长度为70年的月尺度数据,数据依次标为1至840号;本实施例收集了130项气候系统指数,包括88项ACIs、26项SSTs和16项OCIs。
步骤1.2考虑到月径流与预报因子之间存在滞后关系,本实施例根据流域历史时期月径流的自相关性分析确定前期月径流对预报月影响的最长滞时,并以此值统一作为其它预报因子的最长滞时。
本实施例研究流域月径流的自相关性和偏自相关性分析的示意图如图2所示。自相关系数呈现周期性变化,而偏自相关系数在滞后编号为36个月之后,均落入置信区间内,即可理解为预报月以前36个月的径流对该月径流预报结果均有一定的影响,故本实施例选取36个月作为月径流和气候系统指数的最长影响滞时。
步骤2.对训练期的预报因子和月径流数据分别进行归一化处理,再采用基于嵌入式思想的LASSO回归方法自动筛选预报因子。
步骤1中共有131项待选预报因子(即1项月径流数据与130项气候系统指数),均考虑了36个月的最长影响滞时,故待选预报因子个数共有131*36个。如图3所示,给出了模型输入与输出数据的示意图,本实施例以1954年1月至2003年12月为训练期,2004年1月至2020年12月为验证期。
本实施例未直接将总样本集分为训练期和验证期,而是完全依照实际应用需求,采取实时径流预报方法,例如以2003年12月为分界点,逐步将验证期中的数据补入训练期,即训练期的数据长度从600个月逐步延至803个月,共完成204次(即验证期的17年*12个月)训练。
例如第1次训练为:以1954年1月至2003年12月为训练期,预报2004年1月径流;第2次训练为:以1954年1月至2004年1月为训练期,预报2004年2月径流;…;以此类推。在这一步中使用LASSO回归方法进行预报因子的优选,每训练一次模型之前都需要进行一次预报因子的筛选。且在筛选之前应对训练集中的预报因子(输入)和月径流数据(输出)分别进行归一化处理,以消除指标之间的量纲影响。
本实施例以下内容均以第1次训练为例,即以1954年1月至2003年12月为训练期,预报2004年1月径流。使用LASSO回归筛选预报因子的过程如下:首先,采用十折交叉验证法,来确定惩罚参数λ的大小。
如图4所示,本实施例的λ的值为0.00067时,均方误差MSE值最小,所以此值为LASSO回归中最终采用的惩罚参数λ的值。
设训练期的自变量(即预报因子)矩阵为X=(X1,X2,...,Xj,...,Xn),其中Xj=(x1j,x2j,...,xmj),训练期的因变量(即待预报月径流)为Y=(y1,y2,...,ym)T,且如前所述,自变量与因变量均已进行归一化处理。首先建立预报因子X与待预报月径流Y之间的线性模型Y=α+β1X1+β2X2+...+βnXn;其中,α为常数项,n为待选预报因子个数(即131*36个),β为各变量系数。LASSO回归设置最小二乘形式的目标函数如下:
本实施例中惩罚参数λ的值为0.00067,经过LASSO回归后让一些自变量的系数值等于零从而达到了因子筛选的目的。
步骤3.采用基于划分思想的K均值聚类(K-means)方法对训练期样本集进行聚类,应用最近邻规则将样本分成互不重合的K类。步骤3进一步包括以下子步骤:
步骤3.1,采用肘部法则确定聚类数K值,肘部法则的核心指标是误差平方和(SSE):
式中:X为待归类的点,C为聚类中心点。肘部法则的计算原理是成本函数,成本函数是类别畸变程度之和,每个类的畸变程度等于每个变量点到其类别中心的位置距离平方和即误差平方和,若各类内部的成员彼此间越紧凑则类的畸变程度越小。将聚类数K从1开始一次增加,统计误差平方和与聚类类别数K的关系。
本实施例中统计误差平方和与聚类类别数K的关系如图5所示,当类数从1增加到6,误差平方和下降速度较快,类数超过6以后,误差平方和变化变缓。即K=6时为误差平方和的“肘部”,因此6即为本实施例最佳的聚类个数。
步骤3.2随机从600组数列中选中6个点(本实施例中每个点都是一个长度为131*36的行向量),每个点代表每簇的初始聚类中心,完成初始化。
步骤3.3计算剩余各个点到聚类中心的欧式距离,公式为:
式中:N为每个点的维数(131*36)。根据最近邻规则,将其归为离它最近的簇。
步骤3.4重新计算每一个簇的平均值,将其作为新的聚类中心。
步骤3.5更新600组数列的归类结果,即计算每个点到聚类中心的欧氏距离,根据最邻近规则,将其归为离它最近的簇。
步骤3.6重复步骤3.4和步骤3.5,直至6组数列所属的聚类中心不再发生变化或满足设定的迭代次数,则终止全过程。
经过以上步骤,本实施例训练期600组样本被分成6类,每类的样本数分别为214,100,150,50,36和50个。
步骤4.计算验证集预报因子向量与K个训练集的聚类中心的距离,找到距离最近的训练集然后以该数据集训练卷积神经网络(CNN)与门控循环单元(GRU)网络结合的组合深度学习预报模型,即CNN-GRU模型。步骤4进一步包括以下子步骤:
步骤4.1将训练期的数据输入卷积神经网络的第一层,即卷积层。训练期的数据形式为((X1,Y1),(X2,Y2),...,(Xt,Yt),...,(XM,YM)),其中,Xt∈RM×N(1≤t≤M)。式中,X为预报因子即模型输入,Y为待预报月径流即模型输出,M为训练期长度,N为预报因子个数。在这里,M为600,N为18。卷积层包括m个一维卷积核,用于将输入数据转化为特征映射。卷积层的计算公式为Mj=fcov(P*Wj+bj),式中,*代表卷积算子,Mj、Wj和bj分别代表第j个卷积核获得的特征映射,第j个卷积核的权重矩阵,和第j个卷积核的偏置,fcov代表神经网络激活函数,使用ReLU,该函数保留正数元素,同时将负数元素替换为0。
步骤4.2将卷积层输出的特征映射输入最大池化层。最大池化层的主要作用为完成下采样和减少参数量。
步骤4.3一个卷积层加一个池化层构成一组卷积神经网络结构,放置两组卷积神经网络结构,即卷积层1+池化层1+卷积层2+池化层2。本实施例中各层的参数如下:卷积层1中的卷积核个数为8个,尺寸为1×2;卷积层2中的卷积核个数为16个,尺寸为1×3。两个最大池化层的尺寸为1×2。
步骤4.4为了将池化层的输出变成GRU可接受的输入形式,需加入Flatten层进行数据处理。通过Flatten层将池化层的输出转化成一维数列。
步骤4.5放置两层GRU层。本实施例中每层GRU层中的神经元个数设置为128个。每个GRU神经网络单元包含一个重置门和一个更新门。Rt和Zt分别代表重置门和更新门,计算公式分别表示为Rt=σ(WR·[Ht-1,Xt]+bR)和Zt=σ(WZ·[Ht-1,Xt]+bZ)。式中,Ht-1代表记忆单元在t-1时刻的输出状态,WR和WZ分别代表重置门和更新门的权重矩阵,bR和bZ分别代表重置门和更新门的偏置向量,Xt代表记忆单元在t时刻的输入,σ代表激活函数sigmoid。而记忆单元在t时刻的输出状态的计算公式为Ht=Ct,记忆单元在时刻t的新状态Ct计算公式又为Ct=(1-Zt)⊙Ct-1+Zt⊙Ct′,当前可能的单元状态Ct′的计算公式又为Ct′=tanh(WCXt+WC(Rt⊙Ht-1)+bC)。式中,WC和bC分别代表记忆单元的权重矩阵和偏置向量,⊙代表向量元素的乘积,tanh代表激活函数(双曲正切函数)。
步骤4.6本实施例在循环神经网络层后放置两层全连接层,全连接层中的神经元个数设置为128个;并在两层全连接层中加入Dropout层,dropout率设置为0.5,以减少过拟合现象。
步骤4.7在此步之前的操作全部都基于归一化的数值,故需要进行反归一化后再将数据进行输出。
在模型训练过程中,本实施例采用Adam梯度优化算法更新CNN-GRU模型的权重和偏置。Adam算法的目标为使训练期的损失函数值最小,这里选取均方误差(MSE)为损失函数。
步骤5.采用自回归滑动平均(ARMA)模型对预报残差进行实时校正。步骤5进一步包括以下子步骤:
步骤5.1将深度学习模型输出的训练期的所有月径流值作如下处理:处理过后的x值作为预报残差自回归滑动平均模型的输入。
步骤5.2模型定阶。利用AIC定阶准则确定ARMA模型的自回归项阶数p和移动平均项阶数q,最小AIC值对应的模型被认为是最好的模型。本实施例中模型的定阶结果为选取ARMA(1,0)模型。
步骤5.3采用矩估计法对定阶的ARMA模型参数进行估计。本实施例中,ARMA(1,0)模型中的自回归系数为0.3951,滑动平均系数为0,常数项为-0.8566。
步骤5.4根据确定的参数建立ARMA模型。
步骤5.5将本步骤5.1处理过后的预测月之前的月径流数据作为ARMA模型的自变量x,则ARMA模型对应的因变量为y,本实施例中,y=0.3951*x-0.8566。最后,y需进行如下变化才能作为修正后的月径流预报值,
进一步,将以上提供的实时月径流预报方法与现有技术进行比较:
①运行时间。在本实施例中,使用LASSO回归法进行预报因子筛选的运行时间不超过3分钟,而使用逐步回归法的时间为10分钟左右。
②预报精度。根据《南方电网水文气象情报预报规范》,月径流预报精度的计算公式为:
式中,Q代表月径流量,下标s和o分别代表模拟值和实测值,A代表各月的预报精度。预测值大于实测值2倍以上,预报精度按等于0处理。
接下来分三步说明本发明所采用的方法对预报精度的提升效果:
(1)采用逐步回归分析法和本发明LASSO回归法进行预报因子优选的预报结果对比:
如图6所示,采用逐步回归法进行预报因子的筛选,在验证期的平均预报精度为54.9%;采用本发明方法进行预报因子的筛选,在验证期的平均精度为61.4%。
(2)在采用LASSO回归法进行预报因子优选的基础上,不采用K均值聚类与采用本发明聚类方法的预报结果对比:
如图7所示,不采用K均值聚类,在验证期的平均预报精度为61.4%;采用本发明方法K均值聚类,在验证期的平均精度为66.9%。
(3)在采用LASSO回归法进行预报因子优选,采用K均值聚类法对训练集样本进行划分的基础上,不采用自回归滑动模型进行预报结果校正与本发明进行预报结果校正的预报结果对比:
如图8所示,不进行预报结果校正,在验证期的平均预报精度为66.9%;本发明进行预报结果校正,在验证期的平均精度为70.1%。
以上数据充分证实了本发明的技术方案有助于提高月径流的预报效率和精度。
以上实施例仅仅是对本发明技术方案所做的举例说明。本发明所涉及的基于深度学习模型的实时月径流预报方法并不仅仅限定于在以上实施例中所描述的内容,而是以权利要求所限定的范围为准。本发明所属领域技术人员在该实施例的基础上所做的任何修改或补充或等效替换,都在本发明的权利要求所要求保护的范围内。
Claims (6)
1.基于深度学习模型的实时月径流预报方法,其特征在于,包括以下步骤:
步骤1.基于历史信息和未来气象信息收集预报因子,并根据流域历史时期月径流的自相关性分析确定前期月径流对预报月影响的最长滞时,并以此值统一作为其它预报因子的最长影响滞时;历史信息包括:流域当地气象信息,包含大气环流指数ACIs、海温指数SSTs和其它指数OCIs的全球气候指数,包含地表温度、实际蒸散发、潜在蒸散发、归一化植被指数及土壤湿度的多源遥感数据;包含预报月以前的历史月径流数据的前期径流;未来信息包括数值气象预报的月降水和月气温数据;
步骤2.对训练期的预报因子和月径流数据分别进行归一化处理,再采用基于嵌入式思想的LASSO回归方法筛选预报因子;
步骤3.采用基于划分思想的K均值聚类方法对训练期样本集进行聚类,应用最近邻规则将样本分成互不重合的K类;包括如下子步骤:
步骤3.1采用肘部法则确定聚类数K值;核心指标是误差平方和SSE,;式中,X为待归类的点,C为聚类中心点;成本函数是类别畸变程度之和,每个类的畸变程度等于每个变量点到其类别中心的位置距离平方和,若各类内部的成员彼此间越紧凑则类的畸变程度越小;将聚类数K从1开始依次增加,统计误差平方和与聚类类别数K的关系,找到误差平方和下降速度较快和误差平方和变化变缓的临界点,作为最佳的聚类个数;
步骤3.2假设训练期的长度为M个月,每个月对应N个预报因子;随机从M组数列中选中K个点,每个点代表每簇的初始聚类中心,完成初始化;
步骤3.3计算剩余各个点到聚类中心的欧式距离,公式为:;式中,N为每个点的维数;根据最近邻规则,将其归为离它最近的簇;
步骤3.4重新计算每一个簇的平均值,将其作为新的聚类中心;
步骤3.5计算每个点到聚类中心的欧氏距离,根据最邻近规则,将其归为离它最近的簇;
步骤3.6重复步骤3.4和步骤3.5,直至M组数列所属的聚类中心不再发生变化或满足设定的迭代次数,则终止全过程;
步骤4.计算验证集预报因子向量与K个训练集的聚类中心的距离,找到距离最近的训练集然后以该数据集训练卷积神经网络与门控循环单元网络结合的组合深度学习预报模型;
步骤5.采用自回归滑动平均模型对预报残差进行实时校正:
步骤5.1将深度学习模型输出的训练期的所有月径流值作如下处理:,处理过后的x值作为预报残差自回归滑动平均模型的输入;
步骤5.2模型定阶:利用AIC定阶准则确定ARMA模型的自回归项阶数p和移动平均项阶数q,认定最小AIC值对应的模型是最好的模型,完成模型定阶;
步骤5.3采用矩估计法对定阶的ARMA模型参数进行估计;
步骤5.4根据确定的参数建立ARMA模型;
步骤5.5将预测月之前的月径流数据作为ARMA模型的自变量x,则ARMA模型对应的因变量为y,将y进行如下变化后作为修正后的月径流预报值,。
2.根据权利要求1所述的基于深度学习模型的实时月径流预报方法,其特征在于:
其中,在步骤1中,待选预报因子乘lag_k月滞时组成超高维矩阵;
在步骤2中,设训练期的自变量矩阵为预报因子矩阵,其中,训练期的因变量矩阵为待预报月径流矩阵/>;
LASSO回归首先建立预报因子X与待预报月径流Y之间的线性模型;其中,α为常数项,n为待选预报因子乘lag_k月滞时所组成的超高维矩阵中的总因子数,β为各变量系数;
LASSO回归设置最小二乘形式的目标函数如下:;其中λ为惩罚参数,其所在项为惩罚项;在LASSO回归中通过调整惩罚参数λ的大小,以达到产生稀疏解的目的;LASSO回归中惩罚参数λ越大,回归系数压缩幅度越大;采用十折交叉验证法,来确定惩罚参数λ的大小。
3.根据权利要求2所述的基于深度学习模型的实时月径流预报方法,其特征在于:
其中,在步骤2中,是基于分界点将总样本集进行分界,第一次训练先以位于分界点以前的样本作为训练期的数据,预报训练期下一个月的月径流,第二次训练再以第一次训练期的数据和预报数据作为训练期的数据,预报训练期下一个月的月径流;以此类推,按照前述过程完成涵盖总样本集中首月至倒数第二个月时间范围内所有月份数据的训练。
4.根据权利要求1所述的基于深度学习模型的实时月径流预报方法,其特征在于:
其中,步骤4包括如下子步骤:
步骤4.1 将训练期的数据输入卷积神经网络的卷积层;训练期的数据形式为;式中,/>,/>,X为预报因子,Y为待预报月径流,M为训练期长度,N为预报因子个数;卷积层包括m个一维卷积核,用于将输入数据转化为特征映射,卷积层的计算公式为/>,式中,/>代表卷积算子,/>、/>和/>分别代表第j个卷积核获得的特征映射,第j个卷积核的权重矩阵,和第j个卷积核的偏置,代表神经网络激活函数;
步骤4.2 将卷积层输出的特征映射输入池化层;
步骤4.3 一个卷积层加一个池化层构成一组卷积神经网络结构;
步骤4.4 为了将池化层的输出变成循环神经网络GRU可接受的输入形式,需加入Flatten层进行数据处理,通过Flatten层将池化层的输出转化成一维数列;
步骤4.5 放置两层GRU层;每个GRU神经网络单元包含一个重置门和一个更新门, 和 分别代表重置门和更新门,计算公式分别表示为/>和,式中,/>代表记忆单元在t-1时刻的输出状态,/>和/>分别代表重置门和更新门的权重矩阵,/>和/>分别代表重置门和更新门的偏置向量,/>代表记忆单元在t时刻的输入,/>代表激活函数sigmoid,而记忆单元在t时刻的输出状态的计算公式为/>,记忆单元在时刻t的新状态/>计算公式又为/>,当前可能的单元状态/>的计算公式又为/>,式中,/>和/>分别代表记忆单元的权重矩阵和偏置向量, ⊙ 代表向量元素的乘积,/>代表激活函数;
步骤4.6 在循环神经网络层后放置两层全连接层,并在两层全连接层中加入Dropout层,以减少过拟合现象;
步骤4.7 进行反归一化后将数据输出。
5.根据权利要求4所述的基于深度学习模型的实时月径流预报方法,其特征在于:
其中,在步骤4.3中,两组卷积神经网络结构为:卷积层1+池化层1+卷积层2+池化层2。
6.根据权利要求4所述的基于深度学习模型的实时月径流预报方法,其特征在于:
其中,在步骤4.5中,激活函数为双曲正切函数。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110966434.7A CN113705877B (zh) | 2021-08-23 | 2021-08-23 | 基于深度学习模型的实时月径流预报方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110966434.7A CN113705877B (zh) | 2021-08-23 | 2021-08-23 | 基于深度学习模型的实时月径流预报方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN113705877A CN113705877A (zh) | 2021-11-26 |
CN113705877B true CN113705877B (zh) | 2023-09-12 |
Family
ID=78653873
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202110966434.7A Active CN113705877B (zh) | 2021-08-23 | 2021-08-23 | 基于深度学习模型的实时月径流预报方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN113705877B (zh) |
Families Citing this family (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN114282614B (zh) * | 2021-12-27 | 2023-09-26 | 淮阴工学院 | 基于随机森林和ifda优化cnn-gru的中长期径流预测方法 |
CN115115148B (zh) | 2022-08-30 | 2022-12-30 | 武汉大学 | 基于过程-数据协同驱动的长期径流预报方法及系统 |
CN115766607B (zh) * | 2022-11-08 | 2024-09-17 | 陕西通信规划设计研究院有限公司 | 一种基于5g物联网的数据压缩传输方法 |
CN115689368B (zh) * | 2022-11-10 | 2023-08-01 | 华能西藏雅鲁藏布江水电开发投资有限公司 | 一种基于全寿命周期的径流预报模型评价方法 |
CN116187501A (zh) * | 2022-11-29 | 2023-05-30 | 伊金霍洛旗那仁太能源有限公司 | 基于CatBoost模型的低温预测 |
CN116343554B (zh) * | 2023-05-26 | 2023-08-22 | 山东建筑大学 | 一种基于机器学习的智能化信息处理校射雷达仿真方法 |
CN117494862B (zh) * | 2023-08-09 | 2024-05-28 | 长江勘测规划设计研究有限责任公司 | 一种基于假设检验的有限样本下数据驱动径流预报模型优选方法 |
CN117057253B (zh) * | 2023-09-28 | 2023-12-08 | 中国水利水电科学研究院 | 基于空间离散洗牌复形进化算法的水文模型参数率定方法 |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110555561A (zh) * | 2019-09-06 | 2019-12-10 | 清华大学 | 一种中长期径流集合预报方法 |
CN111461453A (zh) * | 2020-04-13 | 2020-07-28 | 中国水利水电科学研究院 | 一种基于多模型组合的中长期径流集合预报方法 |
CN111597758A (zh) * | 2020-05-14 | 2020-08-28 | 河海大学 | 一种基于负相关学习的中小河流集成预报方法 |
AU2020104000A4 (en) * | 2020-12-10 | 2021-02-18 | Guangxi University | Short-term Load Forecasting Method Based on TCN and IPSO-LSSVM Combined Model |
WO2021120787A1 (zh) * | 2019-12-20 | 2021-06-24 | 华中科技大学 | 一种流域干支流大规模水库群模拟调度方法 |
-
2021
- 2021-08-23 CN CN202110966434.7A patent/CN113705877B/zh active Active
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110555561A (zh) * | 2019-09-06 | 2019-12-10 | 清华大学 | 一种中长期径流集合预报方法 |
WO2021120787A1 (zh) * | 2019-12-20 | 2021-06-24 | 华中科技大学 | 一种流域干支流大规模水库群模拟调度方法 |
CN111461453A (zh) * | 2020-04-13 | 2020-07-28 | 中国水利水电科学研究院 | 一种基于多模型组合的中长期径流集合预报方法 |
CN111597758A (zh) * | 2020-05-14 | 2020-08-28 | 河海大学 | 一种基于负相关学习的中小河流集成预报方法 |
AU2020104000A4 (en) * | 2020-12-10 | 2021-02-18 | Guangxi University | Short-term Load Forecasting Method Based on TCN and IPSO-LSSVM Combined Model |
Non-Patent Citations (1)
Title |
---|
水文大数据共享平台研究与设计;陈华;徐坚;肖志远;杨家伟;陈杰;郭生练;许崇育;;水资源研究(第01期);全文 * |
Also Published As
Publication number | Publication date |
---|---|
CN113705877A (zh) | 2021-11-26 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN113705877B (zh) | 基于深度学习模型的实时月径流预报方法 | |
CN108304668B (zh) | 一种结合水文过程数据和历史先验数据的洪水预测方法 | |
CN110471950B (zh) | 一种中小河流实时洪水预报智能模型预报方法 | |
CN112288164B (zh) | 一种计及空间相关性和修正数值天气预报的风功率组合预测方法 | |
CN107992976B (zh) | 热点话题早期发展趋势预测系统及预测方法 | |
CN111860982A (zh) | 一种基于vmd-fcm-gru的风电场短期风电功率预测方法 | |
CN113468803B (zh) | 一种基于改进的woa-gru洪水流量预测方法及系统 | |
CN106503867A (zh) | 一种遗传算法最小二乘风电功率预测方法 | |
CN117035201B (zh) | 平原河网水工程集群多目标调度规则制定方法及系统 | |
CN109143408B (zh) | 基于mlp的动态区域联合短时降水预报方法 | |
CN113139329B (zh) | 一种基于水文相似性和人工神经网络的新安江模型参数率定方法 | |
CN111461453A (zh) | 一种基于多模型组合的中长期径流集合预报方法 | |
CN110674965A (zh) | 基于动态特征选取的多时间步长风功率预测方法 | |
CN116596044A (zh) | 基于多源数据的发电负荷预测模型训练方法及装置 | |
Calp | A hybrid ANFIS-GA approach for estimation of regional rainfall amount | |
CN114117852B (zh) | 一种基于有限差分工作域划分的区域热负荷滚动预测方法 | |
CN113052373A (zh) | 一种基于改进elm模型的月径流变化趋势预测方法 | |
CN114897204A (zh) | 一种海上风电场短期风速预测方法和装置 | |
CN115713144A (zh) | 基于组合cgru模型的短期风速多步预测方法 | |
CN116826737A (zh) | 一种光伏功率的预测方法、装置、存储介质及设备 | |
CN115206444A (zh) | 基于fcm-anfis模型的最佳投药量预测方法 | |
CN117290673A (zh) | 一种基于多模型融合的船舶能耗高精度预测系统 | |
CN116663404A (zh) | 一种耦合人工智能和贝叶斯理论的洪水预报方法及系统 | |
CN115794805A (zh) | 一种中低压配网量测数据补齐方法 | |
Velasco et al. | Performance analysis of multilayer perceptron neural network models in week-ahead rainfall 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 |