CN111597494A - 一种基于非平稳时间序列分解的统计降尺度方法 - Google Patents

一种基于非平稳时间序列分解的统计降尺度方法 Download PDF

Info

Publication number
CN111597494A
CN111597494A CN202010376317.0A CN202010376317A CN111597494A CN 111597494 A CN111597494 A CN 111597494A CN 202010376317 A CN202010376317 A CN 202010376317A CN 111597494 A CN111597494 A CN 111597494A
Authority
CN
China
Prior art keywords
imf
component
decomposition
downscaling
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.)
Granted
Application number
CN202010376317.0A
Other languages
English (en)
Other versions
CN111597494B (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.)
Shandong University
Original Assignee
Shandong University
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 Shandong University filed Critical Shandong University
Priority to CN202010376317.0A priority Critical patent/CN111597494B/zh
Publication of CN111597494A publication Critical patent/CN111597494A/zh
Application granted granted Critical
Publication of CN111597494B publication Critical patent/CN111597494B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • 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)
  • Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • General Physics & Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Algebra (AREA)
  • Computational Mathematics (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Complex Calculations (AREA)

Abstract

本发明涉及一种基于非平稳时间序列分解的统计降尺度方法，该方法的具体步骤包括：S1.降尺度数据准备；S2.非平稳时间序列分解；S3.最优时间序列分解结果的选择；S4.随机森林模型训练；S5.训练后时间序列分解分量的合成；S6.模型评估；S7.未来情景的降尺度。该方法在稳态的分量之间寻找统计关系，进而解决统计降尺度中有关统计关系一致性假设的问题，解决了气候模型和水文模型之间的一直存在的尺度不匹配问题，也为利用水文模型进行未来气候变化等一系列研究中的流域响应提供了一种更为可靠的途径。

Description

一种基于非平稳时间序列分解的统计降尺度方法
技术领域
本发明涉及一种基于非平稳时间序列分解的统计降尺度方法，属于气候和水文交叉领域。
背景技术
气候变化是全球变化的先驱和重要组成部分，数据分析表明，随着社会的发展，人类活动，土地利用以及城市化进程的加快，气候变化产生了多尺度、全方位、多层次的影响。气候变化通过气温、降水、蒸发等因素的改变来影响陆地水循环系统,导致不同时空尺度的水资源重新分布。全球气候模式的出现及其提供的大尺度气候因素的模拟结果为利用水文模型研究区域尺度的水资源的未来和过去的时空特征变化提供了可能。通常，全球气候模型具有较粗的空间分辨率(几百千米),而水文模型的尺度往往在几十千米甚至几千米的范围内，全球气候模型的输出分辨率太低以至于水文模型和全球气候模型之间往往存在尺度差异的问题。全球气候模型的输出无法直接驱动水文模型，也无法直接用于区域水环境或气候环境相关的应用研究。另外，由于气候模型模拟的全球尺度，考虑到计算成本和时间成本，全球气候模型无法抓住全部的区域地势特征细节，比如复杂的地形、水文条件、土地利用情况等等。而某些局部地势特征对当地水文环境的影响是十分显著的，在区域水文气候研究中无法忽视。这种全球气候模型和水文模型之间尺度不匹配的问题限制了全球气候模型输出数据在水文模型中的应用，解决尺度不匹配问题对于利用全球气候模型和水文模型研究水文环境在过去和未来对气候变化的响应是十分重要的。
为此，许多研究致力于发展降尺度技术来弥补这一尺度差距的问题。降尺度技术主要分为三类：动态降尺度、统计降尺度和混合降尺度技术。动态降尺度技术主要指利用区域气候模型来实现尺度下延。区域气候模型具有与全球气候模型相似的物理机制，但是拥有比全球气候模型更细的格网化方案，所以能在模型中刻画更多的区域细节特征，但是区域气候场景仍需要较高时间成本和计算成本。尽管动态降尺度技术已经得到广泛的认可，但是对于更大的研究区域，更紧急的模拟需求，区域气候模型很难快速的产生多种气候场景响应。区域气候模型复杂的物理数学基础也使得动态气候模型的应用受到限制。另外，动态降尺度技术直接依赖全球气候模型的输出作为区域气候模型的边界条件和输入，因此全球气候模型的错误模拟结果往往会被区域气候模型继承，从而导致错误的降尺度结果。
统计降尺度模型由于较高的计算精度和计算效率，简单的建立原理被广泛应用。通常的统计降尺度技术往往是寻找一种大尺度气候序列和小尺度气候序列之间的统计关系。小尺度气候序列往往是地面真实的观测值，统计降尺度技术也可以视为一种对全球气候模型输出的校正，往往与动态降尺度技术连用发展成为一种混合降尺度技术-统计-动态降尺度技术。统计降尺度快速校正了全球气候模型输出的偏倚，防止了区域气候模型继承这种误差。建立一种统计关系的方法多种多样的，例如多个变量和降尺度目标变量之间的线性/非线性的统计关系或者是基于某种分布的统计关系。Wilby等人提出的统计降尺度模型就是一种基于多元线性回归和随机天气发生器的统计降尺度模型，而随机天气发生器是一种基于温度/降水分布的统计方法。Themeβl等人提出了一种分位数映射方法，以纠正从历史降水记录得出的区域气候模型未来降水的经验分布。此外，建立这样统计关系的方法还有典型相关分析、主成分分析、奇异值分解、广义线性模型等。随着机器学习的流行，人工神经网络、支持向量机，深度学习包括长短时记忆人工神经网络、递归神经网络也被纷纷用于统计降尺度模型中。
然而现有的统计降尺度技术往往基于一种假设：气候模型和地面真实气候时间序列之间存在一种一致性的统计关系，即大尺度和小尺度之间的统计关系在未来不会发生改变。但是由于气候自身内部的变率以及人类活动的影响，这种假设是无法成立的。气候系统的高度非线性、非稳态性使得基于过去十几年甚至几十年建立的经验统计关系未来可能不会完全重现。这种假设使得统计降尺度模型在应用到未来场景降尺度的时候备受质疑，然而统计降尺度快速的降尺度能力在目前仍然是无可替代的。
发明内容
针对现有技术的不足，本发明提供一种基于非平稳时间序列分解的统计降尺度方法，将非平稳的气候序列分解为平稳的气候分量，通过寻找平稳分量之间的统计关系克服统计降尺度方法的非一致性假设的问题。
本发明的技术方案为：
一种基于非平稳时间序列分解的统计降尺度方法，具体步骤包括：
S1.降尺度数据准备：
收集全球气候模型GCM的大尺度数据和地面观测的小尺度数据，形成数据集；所述大尺度数据包括降水时间序列、最大温度时间序列、最小温度时间序列；小尺度数据也包括降水时间序列、最大温度时间序列、最小温度时间序列；
大尺度数据作为待降尺度数据；同一变量对应的小尺度数据作为真实值用于大尺度数据的降尺度模型的训练；所述变量包括降水、最大温度、最小温度；小尺度数据是针对小尺度系统来说的，在小尺度系统中，水平尺度约为20km以下，如城市尺度，流域尺度等；大尺度数据是针对大尺度系统来说的，在大尺度系统中，水平尺度约为2000km以上、如全球洋流循环系统、碳循环、氮循环等；
并将数据集划分为训练数据集和测试数据集；训练数据集和测试数据集经过非平稳时间序列分解和最优时间序列选取处理后，训练数据集用于训练降尺度模型，测试数据集用于对降尺度模型进行评估；
S2.非平稳时间序列分解：
使用双变量经验模态分解(BEMD)算法对步骤S1中同一变量对应的大尺度数据和小尺度数据进行非平稳时间序列分解，分别得到本征模态函数(IMF)分量组和一项残差分量，即一组分解结果；大尺度数据分解得到分量的时间尺度和小尺度数据分解得到的分量的时间尺度具有同步性；
通过调整BEMD算法中的投影方向的数量Ndir和迭代停止条件，获得多组分解结果；
通过经验模态分解(Empirical mode decomposition，EMD)算法可以实现非平稳时间序列分解。EMD算法假设原始信号由一组简单的、被认为是稳态的震荡模式叠加而成，这些震荡模式被称为本征模态函数。通过步骤S2实现非平稳气候序列的分解，寻找GCM大尺度和地面观测小尺度数据不同时间尺度的稳态分量之间的对应关系，要求保持GCM和地面观测稳态分量之间时间尺度的同步性。但是，尽管EMD算法可以实现非平稳时间序列分解，但是无法实现两个变量的同时分解，具体表现为对于同一时期同一变量同一站点的GCM和地面观测序列，其分解结果可能会包含不同的分量个数，从而无法直接寻找一种时间尺度的对应关系，而保持大尺度数据和小尺度数据分解的同尺度性是本方法建立降尺度模型的关键。所以，如果在建立模型阶段使用EMD分解将会导致GCM和地面观测值在分解中失去时间尺度的同步性。
本发明中，在建立统计降尺度模型的阶段，非平稳时间序列的分解采用了双变量经验模态分解(BEMD，Bivariate empirical mode decomposition)。BEMD是一种基于EMD(Empirical mode decomposition)的双变量经验模态分解算法，使用多方向投影解决二维数据的分解问题；BEMD算法可以将两个非平稳的时间序列同时分解成两组平稳的、具有规律振荡周期的分量组，保持时间尺度的对应性。
S3.最优时间序列分解结果的选择：
利用正交性指数OI从步骤S2得到的多组分解结果中选择最优时间序列分解结果，即确定最优Ndir值，具体为：
对于降水时间序列、最大温度时间序列、最小温度时间序列的任一组分解结果，分别计算GCM分解结果的正交指数OIGCM和地面观测分解结果的正交指数OIOB，将两个正交指数相加即得到一组分解结果的正交性指标OIT，即OIT＝OIGCM+OIOB；
从多组分解结果的正交性指标OIT中，选择正交性指标OIT最低时对应的Ndir，即为最优Ndir值；
正交性指标OIT越低，意味着该组分解结果分量之间的正交性越强，同时也意味着存在较少的模态混合；从而可以确定投影方向，选择出最优的分解结果对应S2(2)中的Ndir。确定Ndir的分解结果即最优的分解结果，用于训练RF模型。
S4.随机森林(RF)模型训练：
采用最优时间序列分解结果中大尺度数据和小尺度数据的每一IMF分量对和残差分量对训练降尺度模型，降尺度模型为RF模型；GCM的IMF分量和和残差分量是训练的自变量，地面观测的IMF分量和残差分量是训练的因变量；由于经过S2非平稳时间序列分解，GCM的IMF分量与地面观测的IMF分量及GCM的残差分量与地面观测的残差分量具有残差具有时间相关性，两类数据的IMF分量的数量相同，残差分量数量也相同，GCM和地面观测的分解结果形成IMF分量对和残差分量对；
所述RF模型中树的个数为500，最大的叶子节点数为20，节点可分的最小样本数为2；采用MSE、MAE评价一次分裂的结果，最小不纯度为1e-7，叶子节点最小的样本权重和为0；
RF模型是一种树形结构的集成式机器学习模型，由多个彼此独立的决策树组成，利用集体决策的优势降低过拟合的问题。RF模型的优势之一是随机性，主要来源于两个方面：用于训练独立决策树的训练子集是随机的；用于节点分割的数据属性是随机的。这种随机性确保了决策树之间是相互独立的，从而使得误差均匀的分布在整个森林中，而这种误差最终会随着RF模型的集体投票策略而消除，通过集成多个弱决策来实现强的集成。RF模型因为参数简单，预测结果可靠等优点已经得到广泛应用。
S5.训练后时间序列分解分量的合成：
训练后的各个时间尺度的IMF分量和残差分量重新进行逆合成，以获得正常的降尺度降水时间序列、最大温度时间序列、最小温度时间序列；
所述逆合成是将所有分量重新相加，即按照r+IMFi+IMFi-1+…+IMF3+IMF2+IMF1，其中，r代表残差分量，i代表第i个IMF分量；
S6.模型评估：
利用均方根误差、相关性系数、拟合优度对降尺度结果进行精度验证和模型评估；降尺度结果包括步骤S5得到的降水时间序列、最大温度时间序列、最小温度时间序列；
S7.未来情景的降尺度：
S7-1.获得未来情景的GCM的大尺度数据；
S7-2.进行EMD分解得到与已训练RF模型相对应的不同时间尺度的分量结果；
S7-3.利用步骤S4中训练好的RF模型进行对应时间尺度分量的降尺度，获得降尺度后的各时间尺度分量；
S7-4.将RF模型输出的降尺度分量按照步骤S5中的合成方法进行重新合成，从而获得未来的降尺度结果，用于未来区域的水文研究或者气候研究。
根据本发明优选的，步骤S2中，非平稳时间序列分解采用的BEMD算法的具体步骤为：
(1)将某一变量对应的GCM的大尺度数据和地面观测的小尺度数据组成一个二维数据x(t)，t表示时间，所述变量包括最大温度、最小温度和降水；
(2)根据公式
Figure BDA0002480216050000051
确定二维数据x(t)的投影方向
Figure BDA0002480216050000052
其中1≤n≤Ndir；Ndir表示投影方向的数量，n表示第n个投影方向；
(3)对二维数据x(t)进行投影，得到未分解的原始气候序列x(t)在方向
Figure BDA0002480216050000053
上的投影
Figure BDA0002480216050000054
二维数据采用处理复数投影的方式计算投影，即二维数据x(t)的两个维度标记为a(t)，b(t)，a(t)和b(t)代表大尺度和小尺度的气候序列；对以x(t)＝a(t)+b(t)j复数形式的数组进行复数向量的投影，即采用
Figure BDA0002480216050000055
进行数据投影，其中，j表示虚数单位，Re(·)代表数据的实数部分，
Figure BDA0002480216050000056
中
Figure BDA0002480216050000057
依据欧拉公式eiθ＝cosθ+i·sinθ计算；
(4)提取投影
Figure BDA0002480216050000058
的局部极值
Figure BDA0002480216050000059
记录极值位置为
Figure BDA00024802160500000510
(5)使用三次样条插值算法插值步骤(4)中极值，获取包络线，提取的极值数据集
Figure BDA0002480216050000061
获取投影方向
Figure BDA0002480216050000062
上的极值包络线
Figure BDA0002480216050000063
(6)重复步骤(3)-(5)获取所有投影方向上的包络线；
(7)计算所有投影方向上包络线的均值m(t)，
Figure BDA0002480216050000064
(8)二维数据x(t)减去m(t)获得h(t)，IMF本质上是一组具有简单振荡周期的函数，二维数据x(t)减去均值m(t)，相当于去平均化，然后可获得震荡变化的高频成分h(t)，高频成分h(t)即潜在的IMF分量；
判断h(t)是否满足以下条件之一：
8-a、h(t)满足IMF的物理意义，满足IMF的物理意义包括IMF分量极值的个数和过零点的个数之差的绝对值小于1；IMF分量极大值插值的包络线与极小值插值的包络线的均值为0；
8-b、满足迭代停止条件；
8-c、达到最大迭代次数；迭代次数的确定根据数据量和计算资源决定，最大迭代次数是为了防止分解无法收敛，造成内存泄漏；
如果满足上述条件，则进行步骤(9)；
如果不满足上述条件，则令x(t)＝h(t)，重复步骤(2)-(7)，进行迭代运算，直到h(t)满足条件；
(9)将步骤(8)最终得到的h(t)标记为ci(t)，ci(t)即为第i个IMF分量；从x(t)中减去第i个IMF分量，得到第i个IMF分量提取后的剩余ri(t)，即ri(t)＝x(t)-ci(t)；将得到的ri(t)作为新的x(t)；
(10)重复(3)-(9)，直至满足下列条件之一：
10-a、无法再提取IMF分量，即最终得到的残差分量ri+1(t)除两端端点外无极大值和极小值，无法再插值极值成包络线执行(4)之后的步骤进行IMF分量的提取；
10-b、满足迭代停止条件；
10-c、达到最大的迭代次数；
最后提取到第i+1个IMF分量标记为ci+1(t)，提取后的残差分量标记为ri+1(t)＝ri(t)-ci+1(t)；
(11)最终原始的二维数据被表示为
Figure BDA0002480216050000065
K为IMF分量的总数，ck(t)是第k个IMF，rk+1(t)是最终的残差分量，即趋势项。
根据发明优选的，步骤S2中，投影方向的数量Ndir为大于等于2的整数。
Ndir越大，二元信息投影的次数越多，信息保留越丰富，推荐根据计算机的计算能力和应用精度需求设置一组参数值便于进行分解结果的筛选。在不同的Ndir参数的控制下，会得到多组分解结果。
根据本法发明优选的，所述S3中，正交性指数OI的计算公式为：
Figure BDA0002480216050000071
式(I)中，NIMF是本征模态函数IMF的个数；N为IMF的长度，即原始时间序列的长度；IMFik是第i个IMF分量中的第k个值,IMFjK表示第j个IMF分量中的第k个值，xk是未分解的原始气候序列，rk是分解最终得到的残差分量。
根据本发明优选的，步骤S7-2中，进行EMD分解得到与已训练RF模型相对应的不同时间尺度的分量结果，具体步骤为：
(a)提取x(t)极值，此处x(t)是由GCM输出的未来情景大尺度数据组成的一维数组；
(b)采用三次样条插值算法插值获取极值包络线；
(c)计算包络线均值m(t)；
(d)计算h(t)＝x(t)-m(t)，判断h(t)是否满足以下条件之一：
d-1、h(t)满足IMF的物理意义，满足IMF的物理意义包括IMF分量极值的个数和过零点的个数之差的绝对值小于1；IMF分量极大值插值的包络线与极小值插值的包络线的均值为0；
d-2、满足迭代停止条件；
d-3、达到最大迭代次数；
若不满足，则x(t)＝h(t)重复步骤(a)-(d)；
若满足，则执行步骤(e)；
(e)将得到的h(t)标记为ci(t)，ci(t)即为第i个IMF分量；从x(t)中减去第i个IMF分量，得到第i个IMF分量提取后的剩余ri(t)，即ri(t)＝x(t)-ci(t)；将得到的ri(t)作为新的x(t)；
(f)重复(a)-(e)，直到无法继续提取IMF分量或满足迭代停止条件或达到最大的迭代次数；
最终提取到第i+1个IMF分量标记为ci+1(t)，提取后的残差项标记为ri+1(t)＝ri(t)-ci+1(t)；
(g)最终GCM输出的气候序列被表示为
Figure BDA0002480216050000081
K为IMF分量的总数，ck(t)是第k个IMF，rk+1(t)是最终的残差分量，即趋势项。
根据发明优选的，步骤S2及S7-2中迭代停止条件为：对于(1-α)％的数据，α为0到1之间的常数，δ(t)<sd；对于剩下的α％的数据，δ(t)<sd2；δ(t)是迭代的目标函数，且
Figure BDA0002480216050000082
其中m(t)为每次迭代剩余的残差；a(t)为平均振幅，且
Figure BDA0002480216050000083
emax表示极大值插值得到的包络线；emin表示极小值插值得到的包络线；
迭代条件决定了分解算法的结果精度，迭代条件宽松会导致BEMD和EMD分解不彻底，即所有的IMF无法被完全提取出来，最终的残差项仍包含震荡信息。
进一步优选的，sd＝0.05，sd2＝10×sd＝0.5，α＝0.05。
根据本发明优选的，步骤S1中，所述数据集的前80％的数据为训练数据集，用于训练降尺度模型；所述数据集的后20％的数据为测试数据集，用于验证降尺度模型。
本发明的有益效果为：
1.本发明提供的统计降尺度方法旨在通过将非稳态的数据分解为稳态分量，在稳态的分量之间寻找统计关系，进而解决统计降尺度中有关统计关系一致性假设的问题，解决了气候模型和水文模型之间的一直存在的尺度不匹配问题，也为利用水文模型进行未来气候变化等一系列研究中的流域响应提供了一种更为可靠的途径。
2.本发明提供的统计降尺度方法能够有效提取非稳态时间序列中的稳态分量，这些稳态分量具有不同的周期，对应不同时间尺度上的气候循环。相比于无非平稳时间序列分解的普通RF模型、线性回归模型，本发明提供方法得到的温度和降水尺度下延结果具有明显的优势。
3.本发明提供的统计降尺度方法对两个完全不同地形的站点位置的模拟结果显示，该方法具有空间尺度上的强泛化能力，可实现忽略地形因素带来的外在影响。另外，仅仅利用GCM大尺度数据也降低了统计降尺度模型中多源数据带来的输入不确定性。
附图说明
图1是实施例1提供一种基于非平稳时间序列分解的统计降尺度方法的流程示意图。
图2是实施例1采用BEMD对站点1的GCM和地面观测的最小温度时间序列的分解结果示意图。
图3-1是实施例1中站点1的地面观测最小温度时间序列分解结果的显著性检验结果。
图3-2是实施例1中站点1的GCM分解结果最小温度的时间序列的显著性检验结果。
图4是实施例1提供的分解站点1得到的正交性指数OI结果。
图5-1是实施例1提供的站点1最小温度降尺度结果对比示意图。
图5-2是实施例1提供的站点2最小温度降尺度结果对比示意图。
图6-1是实施例1提供的站点1最大温度降尺度结果对比示意图。
图6-2是实施例1提供的站点2最大温度降尺度结果对比示意图。
图7-1是实施例1提供的站点1降水降尺度结果对比示意图。
图7-2是实施例1提供的站点2降水降尺度结果对比示意图。
具体实施方式
下面结合实施例和说明书附图对本发明做进一步说明，但不限于此。
实施例1
本实施例以位于平原站点1和高山站点2的月平均最大温度、最小温度、降水的降尺度为例来说明一种基于非平稳时间序列分解的统计降尺度方法。站点1为平原站点，平均海拔高度为18.1m，为某城市一农场附近的固定气象观测站；站点2为高山站点，平均海拔高度为1890.9m，为某湖泊固定气象观测站。站点1为典型的地中海式气候，冬季温和潮湿，夏季炎热干燥；由于站点2海拔较高，该地区为地中海大陆性气候，夏季温暖干燥。降尺度方法的技术路线图如图1所示，具体步骤包括：
S1.降尺度数据准备：
分别收集站点1和站点2在1950年至2005年期间的CMIP5全球气候模型GCM的大尺度数据和地面观测的小尺度月平均数据，形成数据集；所述大尺度数据包括降水时间序列、最大温度时间序列、最小温度时间序列；小尺度数据也包括降水时间序列、最大温度时间序列、最小温度时间序列；
CMIP5：国际耦合模式比较计划第5阶段，有来自全球20多个研究组，气候模式有24个，其中13个模式包含大气化学模式分量；地球系统模式有11个，其中5个包含完整的陆地、海洋碳循环以及大气化学过程，5个只包含海洋和陆地碳循环，1个只包含海洋碳循环。中国参与CMIP5的气候系统模式有5个。这里以CMIP5包括的模型之一CanESM2为例来进行说明。CanESM2是一个完整的地球系统模型，包括了耦合大气海洋的CanCM4、耦合陆地和大气的碳循环的CTEM以及海洋碳循环模块CMOC，以重现大气、海洋和陆地之间的碳循环过程。
大尺度数据作为待降尺度数据；同一变量对应的小尺度数据作为真实值用于大尺度数据的降尺度模型的训练；所述变量包括降水、最大温度、最小温度；小尺度数据是针对小尺度系统来说的，在小尺度系统中，水平尺度约为20km以下，如城市尺度，流域尺度等；大尺度数据是针对大尺度系统来说的，在大尺度系统中，水平尺度约为2000km以上、如全球洋流循环系统、碳循环、氮循环等；
并将数据集划分为训练数据集和测试数据集；训练数据集和测试数据集经过非平稳时间序列分解和最优时间序列选取处理后，训练数据集用于训练降尺度模型，测试数据集用于对降尺度模型进行评估；
步骤S1中，所述数据集的前80％的数据为训练数据集，用于训练模型；所述数据集的后20％的数据为测试数据集，用于验证模型。
S2.非平稳时间序列分解：
使用双变量经验模态分解(BEMD)算法分别对步骤S1中站点1和站点2同一变量对应的大尺度数据和小尺度数据进行非平稳时间序列分解，分别得到本征模态函数(IMF)分量组和一项残差分量，即一组分解结果；大尺度数据分解得到分量的时间尺度和小尺度数据分解得到的分量的时间尺度具有同步性；
通过调整BEMD算法中的投影方向的数量Ndir和迭代停止条件，获得多组分解结果；
通过经验模态分解(Empirical mode decomposition，EMD)算法可以实现非平稳时间序列分解。EMD算法假设原始信号由一组简单的、被认为是稳态的震荡模式叠加而成，这些震荡模式被称为本征模态函数。通过步骤S2实现非平稳气候序列的分解，寻找GCM大尺度和地面观测小尺度数据不同时间尺度的稳态分量之间的对应关系，要求保持GCM和地面观测稳态分量之间时间尺度的同步性。但是，尽管EMD算法可以实现非平稳时间序列分解，但是无法实现两个变量的同时分解，具体表现为对于同一时期同一变量同一站点的GCM和地面观测序列，其分解结果可能会包含不同的分量个数，从而无法直接寻找一种时间尺度的对应关系，而保持大尺度数据和小尺度数据分解的同尺度性是本方法建立降尺度模型的关键。所以，如果在建立模型阶段使用EMD分解将会导致GCM和地面观测值在分解中失去时间尺度的同步性。
本发明中，在建立统计降尺度模型的阶段，非平稳时间序列的分解采用了双变量经验模态分解(BEMD，Bivariate empirical mode decomposition)。BEMD是一种基于EMD(Empirical mode decomposition)的双变量经验模态分解算法，核心算法一致，区别在于使用多方向投影解决二维数据的分解问题；BEMD算法可以将两个非平稳的时间序列同时分解成两组平稳的、具有规律振荡周期的分量组，保持时间尺度的对应性。
步骤S2中，投影方向的数量Ndir为大于等于2的整数。
Ndir越大，二元信息投影的次数越多，信息保留越丰富，推荐根据计算机的计算能力和应用精度需求设置一组参数值便于进行分解结果的筛选。在不同的Ndir参数的控制下，会得到多组分解结果。
步骤S2中迭代停止条件为：对于(1-α)％的数据，δ(t)<sd；对于剩下的α％的数据，δ(t)<sd2；δ(t)是迭代的目标函数，且
Figure BDA0002480216050000111
其中m(t)为每次迭代剩余的残差；a(t)为平均振幅，且
Figure BDA0002480216050000112
emax表示极大值插值得到的包络线；emi+表示极小值插值得到的包络线；
迭代条件决定了分解算法的结果精度，迭代条件宽松会导致BEMD和EMD分解不彻底，即所有的IMF无法被完全提取出来，最终的残差项仍包含震荡信息。
本实施例中，sd＝0.05，sd2＝10×sd＝0.5，α＝0.05。
步骤S2中，非平稳时间序列分解采用的BEMD算法的具体步骤为：
(1)将某一变量对应的GCM的大尺度数据和地面观测的小尺度数据组成一个二维数据x(t)，t表示时间，所述变量包括最大温度、最小温度和降水；
(2)根据公式
Figure BDA0002480216050000113
确定二维数据x(t)的投影方向
Figure BDA0002480216050000114
其中1≤n≤Ndir；Ndir表示投影方向的数量，n表示第n个投影方向；
(3)对二维数据x(t)进行投影，得到未分解的原始气候序列x(t)在方向
Figure BDA0002480216050000115
上的投影
Figure BDA0002480216050000116
二维数据采用处理复数投影的方式计算投影，即二维数据x(t)的两个维度标记为a(t)，b(t)，a(t)和b(t)代表大尺度和小尺度的气候序列；对以x(t)＝a(t)+b(t)j复数形式的数组进行复数向量的投影，即采用
Figure BDA0002480216050000117
进行数据投影，其中，j表示虚数单位，5e(·)代表数据的实数部分，
Figure BDA0002480216050000118
中
Figure BDA0002480216050000119
依据欧拉公式eiθ＝cosθ+i·sinθ计算；
(4)提取投影
Figure BDA00024802160500001110
的局部极值
Figure BDA00024802160500001111
记录极值位置为
Figure BDA00024802160500001112
(5)使用三次样条插值算法插值步骤(4)中极值，获取包络线，提取的极值数据集
Figure BDA00024802160500001113
获取投影方向
Figure BDA00024802160500001114
上的极值包络线
Figure BDA00024802160500001115
(6)重复步骤(3)-(5)获取所有投影方向上的包络线；
(7)计算所有投影方向上包络线的均值m(t)，
Figure BDA0002480216050000121
(8)二维数据x(t)减去m(t)获得h(t)，IMF本质上是一组具有简单振荡周期的函数，二维数据x(t)减去均值m(t)，相当于去平均化，然后可获得震荡变化的高频成分h(t)，高频成分h(t)即潜在的IMF分量；
判断h(t)是否满足以下条件之一：
8-a、h(t)满足IMF的物理意义，满足IMF的物理意义包括IMF分量极值的个数和过零点的个数之差的绝对值小于1；IMF分量极大值插值的包络线与极小值插值的包络线的均值为0；
8-b、满足迭代停止条件；迭代停止条件为：对于(1-α)％的数据，δ(t)<sd；对于剩下的α％的数据，δ(t)<sd2；δ(t)是迭代的目标函数，且
Figure BDA0002480216050000122
其中m(t)为每次迭代剩余的残差；a(t)为平均振幅，且
Figure BDA0002480216050000123
emax表示极大值插值得到的包络线；emin表示极小值插值得到的包络线；
8-c、达到最大迭代次数；迭代次数的确定根据数据量和计算资源决定，最大迭代次数是为了防止分解无法收敛，造成内存泄漏；
如果满足上述条件，则进行步骤(9)；
如果不满足上述条件，则令x(t)＝h(t)，重复步骤(2)-(7)，进行迭代运算，直到h(t)满足条件；
(9)将步骤(8)最终得到的h(t)标记为ci(t)，ci(t)即为第i个IMF分量；从x(t)中减去第i个IMF分量，得到第i个IMF分量提取后的剩余ri(t)，即ri(t)＝x(t)-ci(t)；将得到的ri(t)作为新的x(t)；
(10)重复(3)-(9)，直至满足下列条件之一：
10-a、无法再提取IMF分量，即最终得到的残差分量ri+1(t)除两端端点外无极大值和极小值，无法再插值极值成包络线执行(4)之后的步骤进行IMF分量的提取；
10-b、满足迭代停止条件；
10-c、达到最大的迭代次数；
最后提取到第i+1个IMF分量标记为ci+1(t)，提取后的残差分量标记为ri+1(t)＝ri(t)-ci+1(t)；
(11)最终原始的二维数据被表示为
Figure BDA0002480216050000131
K为IMF分量的总数，ck(t)是第k个IMF，rk+1(t)是最终的残差分量，即趋势项。
经过步骤S2，原始的气候序列(非平稳时间)序列被分解成一组不定数量的IMF分量和一个残差分量，采用BEMD分解大尺度数据和小尺度数据的IMF分量有着不同的时间周期，在物理上对应着相应时间尺度(小时/日/月/季/年/十年/百年)上的气候循环信息。残差分量代表一种趋势，如果残差分量不是常量的话意味着原始时间序列是非稳态的。通常气候信息分解得到的残差项一般都是非稳态的，因为气候系统不是一成不变的。经过分解可以得到的例如最大温度/最小温度/降水的分解分量，一组分解结果由同一变量的GCM的分解和地面观测值的两组分解的结果组成，来自GCM和地面观测的每对IMF在时间尺度上保持一致，即代表着同一时间尺度的气候循环。
如图2中的分解结果所示，针对站点1的GCM和地面观测的最小温度时间序列，左图为地面观测对应的分解结果，右图为GCM对应的分解结果，经过步骤S2、S3处理后，初始非平稳的最小温度时间序列被分解成10对稳态分量对，其中前9对为IMF分量对，最后一对为残差分量对(趋势项对)，则需要分别训练10个RF模型用于实现未来相应时间尺度分量的降尺度。
表1为站点1最小温度的地面观测的小尺度数据和GCM的大尺度数据分解得到9对IMF分量对的分解周期；
表1
最小温度 IMF1 IMF2 IMF3 IMF4 IMF5 IMF6 IMF7 IMF8 IMF9
地面站点分解周期 2.6 4 6.9 12 19.2 29 56 134.4 224
GCM分解周期 2.6 4 6.7 12 20.4 32 74.7 134.4 224
由如表1可知，图2中的9对IMF分量对，左右IMF分量相互对应，具有相同的时间周期，保持了时间尺度上的同步性。
由图2和表1的结果可知，GCM的大尺度数据和地面观测的小尺度数据被同步分解为稳态分量，保持了时间尺度上的一致性。排列靠前的分量变化频率较高通常代表较短时间尺度内的变化，例如图2中的1至3分量对应表1周期2.6至6.9不等，属于季节内变化；排列靠后的分量变化频率较低，代表更长时间尺度的变化，比如图2中的4至9分量，对应周期12至224不等，属于年际范围的变化。最后一对分量对为残差对，代表着变量的变化趋势，图2中的第10个分量，该分量表明GCM模拟的大尺度最小温度和地面观测最小温度都存在一个逐渐升高的趋势，暗示最小温度的时间序列是一个非稳态的时间变化序列。
经验模态分解的分解结果通常存在着不同程度的模态混合问题。BEMD的核心算法与EMD相同，因此BEMD也存在着明显的模态混合问题，具体如图2中第3分量和第4分量在第300个月和第500个月的位置存在明显的互补部分。模态混合是一种某个IMF混合了一个或多个其他IMF部分信号的现象。由于不同时间尺度的气候信息具有明显的间歇性(比如明显的季节变化)，模态混合在GCM和地面观测的分解中是无法避免的。模态混合会使IMF失去应有物理意义。经过显著性检验，能够剔除存在显著性不明显的结果，确保每个分量的物理意义明确。根据IMF的能量密度和周期对步骤S2得到的分解结果进行显著性检验，具体步骤包括为：A、计算IMF的能量密度
Figure BDA0002480216050000141
Cn(j)表示IMF分量；
B、计算IMF的平均周期为Tn，
C、当IMF的能量密度与IMF的平均周期之间的关系满足：lnEn+lnTn＝0，则步骤S2得到的分解结果具有显著性，确保大部分分解结果的能量密度有至少90％的置信度。
显著性检验的结果如图3-1和图3-2所示，大部分的分解结果都位于90％的置信线以上，证明大部分分解结果都具有物理意义。
S3.最优时间序列分解结果的选择：
利用正交性指数OI从步骤S2得到的多组分解结果中选择最优时间序列分解结果，即确定最优Ndir值，具体为：
对于降水时间序列、最大温度时间序列、最小温度时间序列的任一组分解结果，分别计算GCM分解结果的正交指数OIGCM和地面观测分解结果的正交指数OIOB，将两个正交指数相加即得到一组分解结果的正交性指标OIT，即OIT＝OIGCM+OIOB；
从多组分解结果的正交性指标OIT中，选择正交性指标OIT最低时对应的Ndir，即为最优的最优Ndir值；
正交性指标OIT越低，意味着该组分解结果分量之间的正交性越强，同时也意味着存在较少的模态混合；从而可以确定投影方向，选择出最优的分解结果对应S2(2)中的Ndir。确定Ndir的分解结果即最优的分解结果，用于训练RF模型。
所述S3中，正交性指数OI的计算公式为：
Figure BDA0002480216050000142
式(I)中，NIMF是本征模态函数IMF的个数；N为IMF的长度，即原始时间序列的长度；IMFik是第i个IMF分量中的第k个值,IMFjk表示第j个IMF分量中的第k个值，xk是未分解的原始气候序列，rk是分解最终得到的残差分量。rk即步骤S2中(11)得到的rk+1(t)。
如图4所示，针对站点1最小温度时间序列的多组分解结果，相应的正交指数指标OI结果列在图4，根据OIT可以选择适合最小温度分解的Ndir＝40。
S4.随机森林(RF)模型训练：
采用最优时间序列分解结果中大尺度数据和小尺度数据的每一IMF分量对和残差分量对训练降尺度模型，降尺度模型为RF模型；GCM的IMF分量和和残差分量是训练的自变量，地面观测的IMF分量和残差分量是训练的因变量；由于经过S2非平稳时间序列分解，GCM的IMF分量与地面观测的IMF分量及GCM的残差分量与地面观测的残差分量具有残差具有时间相关性，两类数据的IMF分量的数量相同，残差分量数量也相同，GCM和地面观测的分解结果形成IMF分量对和残差分量对；
所述RF模型中树的个数为500，最大的叶子节点数为20，节点可分的最小样本数为2；采用MSE、MAE评价一次分裂的结果，最小不纯度为1e-7，叶子节点最小的样本权重和为0；
RF模型是一种树形结构的集成式机器学习模型，由多个彼此独立的决策树组成，利用集体决策的优势降低过拟合的问题。RF模型的优势之一是随机性，主要来源于两个方面：用于训练独立决策树的训练子集是随机的；用于节点分割的数据属性是随机的。这种随机性确保了决策树之间是相互独立的，从而使得误差均匀的分布在整个森林中，而这种误差最终会随着RF模型的集体投票策略而消除，通过集成多个弱决策来实现强的集成。RF模型因为参数简单，预测结果可靠等优点已经得到广泛应用。
S5.训练后时间序列分解分量的合成：
训练后的各个时间尺度的IMF分量组和残差分量重新进行逆合成，以获得正常的降尺度降水时间序列、最大温度时间序列、最小温度时间序列；
所述逆合成是将所有分量重新相加，即按照r+IMFi+IMFi-1+…+IMF3+IMF2+IMF1，其中，r代表残差分量；
S6.模型评估：
利用均方根误差、相关性系数、拟合优度对降尺度结果进行精度验证和模型评估；降尺度结果包括步骤S5得到的两个站点的降水时间序列、最大温度时间序列、最小温度时间序列；
本实施例中，将测试数据集中的大尺度数据中的同一变量的各时间尺度分解分量输入训练好的对应尺度的RF模型，进行逆合成后得到测试数据集的大尺度数据对应的降尺度后的最小温度/最大温度/降水时间序列，再和测试数据集中的小尺度对应的最小温度/最大温度/降水时间序列进行评估，评估不同时间尺度下两种最小温度/最大温度/降水数据的拟合优度、均方根误差以及地面观测相关性。
同时对比未经过分解的RF模型、普通的线性回归模型的降尺度结果的拟合优度、均方根误差以及地面观测相关性。
对进行模型的评价，表2、3、4分别为某平原站点1最小温度/最大温度/降水的降尺度的对比评估结果；
表2
变量-最小温度(℃) 拟合优度R<sup>2</sup> 均方根误差RMSE 与地面观测相关性r
本发明 0.94 1.35 0.93
无分解的RF模型 0.84 1.50 0.91
线性回归 0.79 1.99 0.88
表3
变量-最大温度(℃) 拟合优度R<sup>2</sup> 均方根误差RMSE 与地面观测相关性r
本发明 0.96 2.19 0.97
无分解的RF模型 0.94 2.00 0.92
线性回归 0.87 2.99 0.90
表4
变量-降水(mm) 拟合优度R<sup>2</sup> 均方根误差RMSE 与地面观测相关性r
本发明 0.53 1.82 0.68
无分解的RF模型 0.41 1.65 0.57
线性回归 0.33 2.97 0.36
由表2、3、4可知，本发明提供的方法对地面站点位置的降尺度结果进行评估：拟合优度越接近于1表示所构建模型的稳定性越高，本发明建立的基于非平稳时间序列分解的降尺度模型稳定性皆优于未进行分解的机器学习模型RF和简单的线性回归模型；相关性系数越接近1，均方根误差越小，说明拟合的结果越接近真实的地面观测值，而本方法的相关性系数和均方根误差仍优于未分解的RF模型和简单的线性回归模型。本发明提供的采用非平稳时间序列分解的RF模型在均方根误差、相关性系数、拟合优度上均优于不采用非平稳时间序列分解的RF模型以及采用线性回归对地面站点位置的降尺度结果的预测。由于GCM对降水的模拟结果较差，温度的降尺度结果明显优于降水的降尺度结果。
S7.未来情景的降尺度：
S7-1.获得未来情景的GCM的大尺度数据；
S7-2.进行EMD分解得到与已训练RF模型相对应的不同时间尺度的分量结果；
步骤S7-2中，进行EMD分解得到与已训练RF模型相对应的不同时间尺度的分量结果，具体步骤为：
(a)提取x(t)极值，此处x(t)是由GCM输出的大尺度数据组成的一维数组；
(b)采用三次样条插值算法插值获取极值包络线；
(c)计算包络线均值m(t)；
(d)计算h(t)＝x(t)-m(t)，判断h(t)是否满足以下条件之一：
d-1、h(t)满足IMF的物理意义，满足IMF的物理意义包括IMF分量极值的个数和过零点的个数之差的绝对值小于1；IMF分量极大值插值的包络线与极小值插值的包络线的均值为0；
d-2、满足迭代停止条件；
d-3、达到最大迭代次数；
若不满足，则x(t)＝h(t)重复步骤(a)-(d)；
若满足，则执行步骤(e)；
(e)将得到的h(t)标记为ci(t)，ci(t)即为第i个IMF分量；从x(t)中减去第i个IMF分量，得到第i个IMF分量提取后的剩余ri(t)，即ri(t)＝x(t)-ci(t)；将得到的ri(t)作为新的x(t)；
(f)重复(a)-(e)，直到无法继续提取IMF分量或满足迭代停止条件或达到最大的迭代次数；
最终提取到第i+1个IMF分量标记为ci+1(t)，提取后的残差项标记为ri+1(t)＝ri(t)-ci+1(t)；
(g)最终GCM输出的气候序列被表示为
Figure BDA0002480216050000171
K为IMF分量的总数，ck(t)是第k个IMF，rk+1(t)是最终的残差分量，即趋势项。
S7-3.利用步骤S4中训练好的RF模型进行对应时间尺度分量的降尺度，获得降尺度后的各时间尺度分量；
S7-4.将RF模型输出的降尺度分量按照步骤S5中的合成方法进行重新合成，从而获得未来的降尺度结果，用于未来区域的水文研究或者气候研究。
利用本发明提供的基于非平稳时间序列分解的统计降尺度方法，获取某平原站点1和高山站点2最小温度/最大温度/降水的降尺度结果：图5-1是站点1最小温度的降尺度结果、GCM大尺度数据及地面观测对比图；图5-2是站点2最小温度的降尺度结果、GCM大尺度数据及地面观测对比图；图6-1是站点1最大温度的降尺度结果、GCM大尺度数据及地面观测对比图；图6-2是站点2最大温度的降尺度结果、GCM大尺度数据及地面观测对比图；图7-1是站点1降水的降尺度结果、GCM大尺度数据及地面观测对比图；图7-2是站点2降水的降尺度结果、GCM大尺度数据及地面观测对比图。图5、6、7中，地面观测值由虚线表示，GCM大尺度数据由实线表示，降尺度结果由点状线表示。对于平原站点1的最小温度(图5-1)，GCM输出的大尺度数据对该站点冬季的最小温度存在明显的低估，同时也低估了该站夏季的最小温度，但是对夏季高温时期的低估程度要低于该站冬季的低温时期。对于平原站点1的最大温度(图6-1)，GCM的大尺度结果依旧存在低估，尤其是GCM模拟的冬季最大温度偏低。可能的原因是GCM中没有充分考虑部分下垫面信息，而该站点位于某城区附近的一个农场内(海拔高度18.1m)，站点附近受人类活动影响下垫面信息较为复杂，GCM没有对下垫面及人类活动进行充分的描述导致冬季模拟的温度偏低。降尺度后的站点1的最小温度序列和最大温度纠正了被GCM低估的冬季最小温度和GCM模拟偏低的冬季最大温度。对于高山站点2的最小温度(图5-2)和最大温度(图6-2)，大尺度的GCM的模拟结果存在明显的高估：GCM模拟结果高估了站点2夏季的最大温度，大尺度GCM与地面观测的最大温度差距将近10℃；同时对于最小温度，无论是夏季还是冬季都存在严重的高估。这种高估可能是GCM并没有充分考虑站点2所在地区的高山信息。降尺度后的高山站点2的最大温度和最小温度更接近于实际观测值。对于降水(图7)，无论是平原站点1还是高山站点2，GCM的大尺度数据与地面观测的降雨数据存在较大的偏差，经过降尺度后，两个站点的降水时间序列和地面观测值的相关性系数上升。但是，由于GCM对降水的模拟效果较差，降水的降尺度结果明显不如温度的降尺度结果。图5、6、7中竖线的左侧训练阶段，右侧为测试阶段，训练阶段和测试阶段，降尺度结果出现没有很大偏差，也说明了所建降尺度模型性能稳定可靠。站点1和站点2的最小温度/最大温度/降水的降尺度结果也显示本发明提出的基于非平稳时间序列分解的降尺度方法，可以不用额外考虑地势因素的影响，例如传统的简单线性回归往往需要考虑包括海拔高度等在内的地势信息和包括植被归一化指数在内的下垫面等土地利用信息，具有较强的区域泛化能力，同时也减少了统计降尺度方法中多源数据所带来的输入不确定性。

Claims (7)

1.一种基于非平稳时间序列分解的统计降尺度方法，其特征在于，具体步骤包括：
S1.降尺度数据准备：
收集全球气候模型GCM的大尺度数据和地面观测的小尺度数据，形成数据集；所述大尺度数据包括降水时间序列、最大温度时间序列、最小温度时间序列；小尺度数据也包括降水时间序列、最大温度时间序列、最小温度时间序列；
大尺度数据作为待降尺度数据；同一变量对应的小尺度数据作为真实值用于大尺度数据的降尺度模型的训练；所述变量包括降水、最大温度、最小温度；
并将数据集划分为训练数据集和测试数据集；
S2.非平稳时间序列分解：
使用BEMD算法对步骤S1中同一变量对应的大尺度数据和小尺度数据进行非平稳时间序列分解，分别得到IMF分量组和一项残差分量，即一组分解结果；大尺度数据分解得到分量的时间尺度和小尺度数据分解得到的分量的时间尺度具有同步性；
通过调整BEMD算法中的投影方向的数量Ndir和迭代停止条件，获得多组分解结果；
S3.最优时间序列分解结果的选择：
利用正交性指数OI从步骤S2得到的多组分解结果中选择最优时间序列分解结果，即确定最优Ndir值，具体为：
对于降水时间序列、最大温度时间序列、最小温度时间序列的任一组分解结果，分别计算GCM分解结果的正交指数OIGCM和地面观测分解结果的正交指数OIoB，将两个正交指数相加即得到一组分解结果的正交性指标OIT，即OIT＝OIGCM+OIoB；
从多组分解结果的正交性指标OIT中，选择正交性指标OIT最低时对应的Ndir，即为最优Ndir值；
S4.RF模型训练：
采用最优时间序列分解结果中大尺度数据和小尺度数据的每一IMF分量对和残差分量对训练降尺度模型，降尺度模型为RF模型；GCM的IMF分量和和残差分量是训练的自变量，地面观测的IMF分量和残差分量是训练的因变量；
所述RF模型中树的个数为500，最大的叶子节点数为20，节点可分的最小样本数为2；采用MSE、MAE评价一次分裂的结果，最小不纯度为1e-7，叶子节点最小的样本权重和为0；
S5.训练后时间序列分解分量的合成：
训练后的各个时间尺度的IMF分量和残差分量重新进行逆合成，以获得正常的降尺度降水时间序列、最大温度时间序列、最小温度时间序列；
所述逆合成是将所有分量重新相加，即按照r+IMFi+IMFi-1+…+IMF3+IMF2+IMF1，其中，r代表残差分量，i代表第i个IMF分量；
S6.模型评估：
利用均方根误差、相关性系数、拟合优度对降尺度结果进行精度验证和模型评估；降尺度结果包括步骤S5得到的降水时间序列、最大温度时间序列、最小温度时间序列；
S7.未来情景的降尺度：
S7-1.获得未来情景的GCM的大尺度数据；
S7-2.进行EMD分解得到与已训练RF模型相对应的不同时间尺度的分量结果；
S7-3.利用步骤S4中训练好的RF模型进行对应时间尺度分量的降尺度，获得降尺度后的各时间尺度分量；
S7-4.将RF模型输出的降尺度分量按照步骤S5中的合成方法进行重新合成，从而获得未来的降尺度结果，用于未来区域的水文研究或者气候研究。
2.根据权利要求1所述的一种基于非平稳时间序列分解的统计降尺度方法，其特征在于，步骤S2中，非平稳时间序列分解采用的BEMD算法的具体步骤为：
(1)将某一变量对应的GCM的大尺度数据和地面观测的小尺度数据组成一个二维数据x(t)，t表示时间，所述变量包括最大温度、最小温度和降水；
(2)根据公式
Figure FDA0002480216040000021
确定二维数据x㈦的投影方向
Figure FDA0002480216040000022
其中1≤n≤Ndir；Ndir表示投影方向的数量，n表示第n个投影方向；
(3)对二维数据x(t)进行投影，得到未分解的原始气候序列x(t)在方向
Figure FDA0002480216040000023
上的投影
Figure FDA0002480216040000024
二维数据采用处理复数投影的方式计算投影，即二维数据x(t)的两个维度标记为a(t)，b(t)，a(t)和b(t)代表大尺度和小尺度的气候序列；对以x(t)＝a(t)+b(t)j复数形式的数组进行复数向量的投影，即采用
Figure FDA0002480216040000025
进行数据投影，其中，j表示虚数单位，Re(-)代表数据的实数部分，
Figure FDA0002480216040000026
中
Figure FDA0002480216040000027
依据欧拉公式
Figure FDA0002480216040000028
计算；
(4)提取投影
Figure FDA0002480216040000029
的局部极值
Figure FDA00024802160400000210
记录极值位置为
Figure FDA00024802160400000211
(5)使用三次样条插值算法插值步骤(4)中极值，获取包络线，提取的极值数据集
Figure FDA0002480216040000031
获取投影方向
Figure FDA0002480216040000032
上的极值包络线
Figure FDA0002480216040000033
(6)重复步骤(3)-(5)获取所有投影方向上的包络线；
(7)计算所有投影方向上包络线的均值m(t)，
Figure FDA0002480216040000034
(8)二维数据x(t)减去m(t)获得h(t)，判断h(t)是否满足以下条件之一：
8-a、h(t)满足IMF的物理意义，满足IMF的物理意义包括IMF分量极值的个数和过零点的个数之差的绝对值小于1；IMF分量极大值插值的包络线与极小值插值的包络线的均值为0；
8-b、满足迭代停止条件；
8-c、达到最大迭代次数；
如果满足上述条件，则进行步骤(9)；
如果不满足上述条件，则令x(t)＝h(t)，重复步骤(2)-(7)，进行迭代运算，直到h(t)满足条件；
(9)将步骤(8)最终得到的h(t)标记为ci(t)，ci(t)即为第i个IMF分量；从x(t)中减去第i个IMF分量，得到第i个IMF分量提取后的剩余ri(t)，即ri(t)＝x(t)-ci(t)；将得到的ri(t)作为新的x(t)；
(10)重复(3)-(9)，直至满足下列条件之一：
10-a、无法再提取IMF分量，即最终得到的残差分量ri+1(t)除两端端点外无极大值和极小值，无法再插值极值成包络线执行(4)之后的步骤进行IMF分量的提取；
10-b、满足迭代停止条件；
10-c、达到最大的迭代次数；
最后提取到第i+1个IMF分量标记为ci+1(t)，提取后的残差分量标记为ri+1(t)＝ri(t)-ci+1(t)；
(11)最终原始的二维数据被表示为
Figure FDA0002480216040000035
K为IMF分量的总数，ck(t)是第k个IMF，rk+1(t)是最终的残差分量，即趋势项。
3.根据权利要求1所述的一种基于非平稳时间序列分解的统计降尺度方法，其特征在于，步骤S2中，投影方向的数量Ndir为大于等于2的整数。
4.根据权利要求1所述的一种基于非平稳时间序列分解的统计降尺度方法，其特征在于，所述S3中，正交性指数OI的计算公式为：
Figure FDA0002480216040000041
式(I)中，NIMF是本征模态函数IMF的个数；N为IMF的长度，即原始时间序列的长度；IMFik是第i个IMF分量中的第k个值，IMFjk表示第j个IMF分量中的第k个值，xk是未分解的原始气候序列，rk是分解最终得到的残差分量。
5.根据权利要求1所述的一种基于非平稳时间序列分解的统计降尺度方法，其特征在于，步骤S7-2中，进行EMD分解得到与已训练RF模型相对应的不同时间尺度的分量结果，具体步骤为：
(a)提取x(t)极值，此处x(t)是由GCM输出的未来情景大尺度数据组成的一维数组；
(b)采用三次样条插值算法插值获取极值包络线；
(c)计算包络线均值m(t)；
(d)计算h(t)＝x(t)-m(t)，判断h(t)是否满足以下条件之一：
d-1、h(t)满足IMF的物理意义，满足IMF的物理意义包括IMF分量极值的个数和过零点的个数之差的绝对值小于1；IMF分量极大值插值的包络线与极小值插值的包络线的均值为0；
d-2、满足迭代停止条件；
d-3、达到最大迭代次数；
若不满足，则x(t)＝h(t)重复步骤(a)-(d)；
若满足，则执行步骤(e)；
(e)将得到的h(t)标记为ci(t)，ci(t)即为第i个IMF分量；从x(t)中减去第i个IMF分量，得到第i个IMF分量提取后的剩余ri(t)，即ri(t)＝x(t)-ci(t)；将得到的ri(t)作为新的x(t)；
(f)重复(a)-(e)，直到无法继续提取IMF分量或满足迭代停止条件或达到最大的迭代次数；
最终提取到第i+1个IMF分量标记为ci+1(t)，提取后的残差项标记为ri+1(t)＝ri(t)-ci+1(t)；
(g)最终GCM输出的气候序列被表示为
Figure FDA0002480216040000042
K为IMF分量的总数，ck(t)是第k个IMF，rk+1(t)是最终的残差分量，即趋势项。
6.根据权利要求5所述的一种基于非平稳时间序列分解的统计降尺度方法，其特征在于，步骤S2及S7-2中迭代停止条件为：对于(1-α)％的数据，δ(t)＜sd；对于剩下的α％的数据，6(t)＜sd2；δ(t)是迭代的目标函数，且
Figure FDA0002480216040000051
其中m(t)为每次迭代剩余的残差；a(t)为平均振幅，且
Figure FDA0002480216040000052
emax表示极大值插值得到的包络线；emin表示极小值插值得到的包络线。
7.根据权利要求1-6任一项所述的一种基于非平稳时间序列分解的统计降尺度方法，其特征在于，步骤S1中，所述数据集的前80％的数据为训练数据集，用于训练降尺度模型；所述数据集的后20％的数据为测试数据集，用于验证降尺度模型。
CN202010376317.0A 2020-05-07 2020-05-07 一种基于非平稳时间序列分解的统计降尺度方法 Active CN111597494B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010376317.0A CN111597494B (zh) 2020-05-07 2020-05-07 一种基于非平稳时间序列分解的统计降尺度方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010376317.0A CN111597494B (zh) 2020-05-07 2020-05-07 一种基于非平稳时间序列分解的统计降尺度方法

Publications (2)

Publication Number Publication Date
CN111597494A true CN111597494A (zh) 2020-08-28
CN111597494B CN111597494B (zh) 2022-06-17

Family

ID=72186859

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010376317.0A Active CN111597494B (zh) 2020-05-07 2020-05-07 一种基于非平稳时间序列分解的统计降尺度方法

Country Status (1)

Country Link
CN (1) CN111597494B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113219555A (zh) * 2020-09-30 2021-08-06 国家气候中心 一种基于统计降尺度技术的短期气候预测方法
CN113705657A (zh) * 2021-08-24 2021-11-26 华北电力大学 一种基于差分法消除多重共线性的逐步聚类统计降尺度方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100092028A1 (en) * 2008-10-10 2010-04-15 National Central University Data Decomposition Method and Computer System Therefrom
CN105046044A (zh) * 2015-05-29 2015-11-11 上海大学 基于最优小波包变换的非平稳风速预测方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100092028A1 (en) * 2008-10-10 2010-04-15 National Central University Data Decomposition Method and Computer System Therefrom
CN105046044A (zh) * 2015-05-29 2015-11-11 上海大学 基于最优小波包变换的非平稳风速预测方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
万仕全等: "基于EMD方法的观测数据信息提取与预测研究", 《气象学报》 *
胡娅敏 等: "基于多时间尺度的回归集成预测模型", 《气象》 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113219555A (zh) * 2020-09-30 2021-08-06 国家气候中心 一种基于统计降尺度技术的短期气候预测方法
CN113219555B (zh) * 2020-09-30 2021-10-26 国家气候中心 一种基于统计降尺度技术的短期气候预测方法
CN113705657A (zh) * 2021-08-24 2021-11-26 华北电力大学 一种基于差分法消除多重共线性的逐步聚类统计降尺度方法
CN113705657B (zh) * 2021-08-24 2024-01-19 华北电力大学 一种基于差分法消除多重共线性的逐步聚类统计降尺度方法

Also Published As

Publication number Publication date
CN111597494B (zh) 2022-06-17

Similar Documents

Publication Publication Date Title
Huang et al. Forecasting hourly solar irradiance using hybrid wavelet transformation and Elman model in smart grid
CN107104442B (zh) 计及参数模糊性的含风电场电力系统概率潮流计算方法
CN111144644B (zh) 基于变分异方差高斯过程回归的短期风速预测方法
CN111597494B (zh) 一种基于非平稳时间序列分解的统计降尺度方法
Sekertekin et al. Short-term air temperature prediction by adaptive neuro-fuzzy inference system (ANFIS) and long short-term memory (LSTM) network
CN112651108B (zh) 一种解耦气象要素和植被动态对水文要素影响的方法
Isaksson et al. Solar power forecasting with machine learning techniques
WO2021063461A1 (en) Method for planning a layout of a renewable energy site
Aliberti et al. Forecasting Short-term Solar Radiation for Photovoltaic Energy Predictions.
Bhola et al. Estimation of solar radiation using support vector regression
CN109671019A (zh) 一种基于多目标优化算法和稀疏表达的遥感影像亚像元制图方法
CN115759389A (zh) 基于天气类型的相似日组合策略的日前光伏功率预测方法
CN106405683B (zh) 基于g-l混合噪声特性核岭回归技术的风速预报方法及装置
Masud et al. Monitoring and predicting landuse/landcover change using an integrated markov chain & multilayer perceptron models: A case study of sahiwal tehsil
Shi et al. Deep-learning-based wind speed forecasting considering spatial–temporal correlations with adjacent wind turbines
CN107844872B (zh) 一种用于风力发电的短期风速预报方法
CN114417728A (zh) 基于温度和发射率及深度学习的近地表空气温度反演方法
Kulkarni et al. Output estimation of solar Photovoltaic (PV) system
Bürger On the disaggregation of climatological means and anomalies
Bantupalli et al. Wind Speed forecasting using empirical mode decomposition with ANN and ARIMA models
CN110175639B (zh) 一种基于特征选取的短期风电功率预测方法
CN113449920A (zh) 一种风电功率预测方法、系统及计算机可读介质
Su et al. A feature importance analysis based solar irradiance mapping model using multi-channel satellite remote sensing data
CN113298304A (zh) 基于多维气象预报因子图卷积的光伏发电功率预测方法
CN112417768A (zh) 一种基于藤结构Pair-Copula的风电相关性条件采样方法

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