技术摘要：
本发明涉及一种基于非平稳时间序列分解的统计降尺度方法，该方法的具体步骤包括：S1.降尺度数据准备；S2.非平稳时间序列分解；S3.最优时间序列分解结果的选择；S4.随机森林模型训练；S5.训练后时间序列分解分量的合成；S6.模型评估；S7.未来情景的降尺度。该方法在稳态  全部
背景技术：
气候变化是全球变化的先驱和重要组成部分，数据分析表明，随着社会的发展，人 类活动，土地利用以及城市化进程的加快，气候变化产生了多尺度、全方位、多层次的影响。 气候变化通过气温、降水、蒸发等因素的改变来影响陆地水循环系统,导致不同时空尺度的 水资源重新分布。全球气候模式的出现及其提供的大尺度气候因素的模拟结果为利用水文 模型研究区域尺度的水资源的未来和过去的时空特征变化提供了可能。通常，全球气候模 型具有较粗的空间分辨率(几百千米),而水文模型的尺度往往在几十千米甚至几千米的范 围内，全球气候模型的输出分辨率太低以至于水文模型和全球气候模型之间往往存在尺度 差异的问题。全球气候模型的输出无法直接驱动水文模型，也无法直接用于区域水环境或 气候环境相关的应用研究。另外，由于气候模型模拟的全球尺度，考虑到计算成本和时间成 本，全球气候模型无法抓住全部的区域地势特征细节，比如复杂的地形、水文条件、土地利 用情况等等。而某些局部地势特征对当地水文环境的影响是十分显著的，在区域水文气候 研究中无法忽视。这种全球气候模型和水文模型之间尺度不匹配的问题限制了全球气候模 型输出数据在水文模型中的应用，解决尺度不匹配问题对于利用全球气候模型和水文模型 研究水文环境在过去和未来对气候变化的响应是十分重要的。 为此，许多研究致力于发展降尺度技术来弥补这一尺度差距的问题。降尺度技术 主要分为三类：动态降尺度、统计降尺度和混合降尺度技术。动态降尺度技术主要指利用区 域气候模型来实现尺度下延。区域气候模型具有与全球气候模型相似的物理机制，但是拥 有比全球气候模型更细的格网化方案，所以能在模型中刻画更多的区域细节特征，但是区 域气候场景仍需要较高时间成本和计算成本。尽管动态降尺度技术已经得到广泛的认可， 但是对于更大的研究区域，更紧急的模拟需求，区域气候模型很难快速的产生多种气候场 景响应。区域气候模型复杂的物理数学基础也使得动态气候模型的应用受到限制。另外，动 态降尺度技术直接依赖全球气候模型的输出作为区域气候模型的边界条件和输入，因此全 球气候模型的错误模拟结果往往会被区域气候模型继承，从而导致错误的降尺度结果。 统计降尺度模型由于较高的计算精度和计算效率，简单的建立原理被广泛应用。 通常的统计降尺度技术往往是寻找一种大尺度气候序列和小尺度气候序列之间的统计关 系。小尺度气候序列往往是地面真实的观测值，统计降尺度技术也可以视为一种对全球气 候模型输出的校正，往往与动态降尺度技术连用发展成为一种混合降尺度技术-统计-动态 降尺度技术。统计降尺度快速校正了全球气候模型输出的偏倚，防止了区域气候模型继承 这种误差。建立一种统计关系的方法多种多样的，例如多个变量和降尺度目标变量之间的 线性/非线性的统计关系或者是基于某种分布的统计关系。Wilby等人提出的统计降尺度模 型就是一种基于多元线性回归和随机天气发生器的统计降尺度模型，而随机天气发生器是 6 CN 111597494 A 说　明　书 2/14 页 一种基于温度/降水分布的统计方法。Themeβl等人提出了一种分位数映射方法，以纠正从 历史降水记录得出的区域气候模型未来降水的经验分布。此外，建立这样统计关系的方法 还有典型相关分析、主成分分析、奇异值分解、广义线性模型等。随着机器学习的流行，人工 神经网络、支持向量机，深度学习包括长短时记忆人工神经网络、递归神经网络也被纷纷用 于统计降尺度模型中。 然而现有的统计降尺度技术往往基于一种假设：气候模型和地面真实气候时间序 列之间存在一种一致性的统计关系，即大尺度和小尺度之间的统计关系在未来不会发生改 变。但是由于气候自身内部的变率以及人类活动的影响，这种假设是无法成立的。气候系统 的高度非线性、非稳态性使得基于过去十几年甚至几十年建立的经验统计关系未来可能不 会完全重现。这种假设使得统计降尺度模型在应用到未来场景降尺度的时候备受质疑，然 而统计降尺度快速的降尺度能力在目前仍然是无可替代的。
技术实现要素：
针对现有技术的不足，本发明提供一种基于非平稳时间序列分解的统计降尺度方 法，将非平稳的气候序列分解为平稳的气候分量，通过寻找平稳分量之间的统计关系克服 统计降尺度方法的非一致性假设的问题。 本发明的技术方案为： 一种基于非平稳时间序列分解的统计降尺度方法，具体步骤包括： S1.降尺度数据准备： 收集全球气候模型GCM的大尺度数据和地面观测的小尺度数据，形成数据集；所述 大尺度数据包括降水时间序列、最大温度时间序列、最小温度时间序列；小尺度数据也包括 降水时间序列、最大温度时间序列、最小温度时间序列； 大尺度数据作为待降尺度数据；同一变量对应的小尺度数据作为真实值用于大尺 度数据的降尺度模型的训练；所述变量包括降水、最大温度、最小温度；小尺度数据是针对 小尺度系统来说的，在小尺度系统中，水平尺度约为20km以下，如城市尺度，流域尺度等；大 尺度数据是针对大尺度系统来说的，在大尺度系统中，水平尺度约为2000km以上、如全球洋 流循环系统、碳循环、氮循环等； 并将数据集划分为训练数据集和测试数据集；训练数据集和测试数据集经过非平 稳时间序列分解和最优时间序列选取处理后，训练数据集用于训练降尺度模型，测试数据 集用于对降尺度模型进行评估； S2.非平稳时间序列分解： 使用双变量经验模态分解(BEMD)算法对步骤S1中同一变量对应的大尺度数据和 小尺度数据进行非平稳时间序列分解，分别得到本征模态函数(IMF)分量组和一项残差分 量，即一组分解结果；大尺度数据分解得到分量的时间尺度和小尺度数据分解得到的分量 的时间尺度具有同步性； 通过调整BEMD算法中的投影方向的数量Ndir和迭代停止条件，获得多组分解结果； 通过经验模态分解(Empirical  mode  decomposition，EMD)算法可以实现非平稳 时间序列分解。EMD算法假设原始信号由一组简单的、被认为是稳态的震荡模式叠加而成， 这些震荡模式被称为本征模态函数。通过步骤S2实现非平稳气候序列的分解，寻找GCM大尺 7 CN 111597494 A 说　明　书 3/14 页 度和地面观测小尺度数据不同时间尺度的稳态分量之间的对应关系，要求保持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.训练后时间序列分解分量的合成： 8 CN 111597494 A 说　明　书 4/14 页 训练后的各个时间尺度的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)根据公式 确定二维数据x(t)的投影方向 其中1≤n≤Ndir；Ndir表示 投影方向的数量，n表示第n个投影方向； (3)对二维数据x(t)进行投影，得到未分解的原始气候序列x(t)在方向 上的投 影 二维数据采用处理复数投影的方式计算投影，即二维数据x(t)的两个维度标记为 a(t)，b(t)，a(t)和b(t)代表大尺度和小尺度的气候序列；对以x(t)＝a(t) b(t)j复数形式 的数组进行复数向量的投影，即采用 进行数据投影，其中，j表示 虚数单位，Re(·)代表数据的实数部分， 中 依据欧拉公式eiθ＝cosθ i·sinθ计 算； (4)提取投影 的局部极值 记录极值位置为 (5)使用三次样条插值算法插值步骤(4)中极值，获取包络线，提取的极值数据集 获取投影方向 上的极值包络线 (6)重复步骤(3)-(5)获取所有投影方向上的包络线； (7)计算所有投影方向上包络线的均值m(t)， (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分量极值的个数和过 9 CN 111597494 A 说　明　书 5/14 页 零点的个数之差的绝对值小于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)最终原始的二维数据被表示为 K为IMF分量的总 数，ck(t)是第k个IMF，rk 1(t)是最终的残差分量，即趋势项。 根据发明优选的，步骤S2中，投影方向的数量Ndir为大于等于2的整数。 Ndir越大，二元信息投影的次数越多，信息保留越丰富，推荐根据计算机的计算能 力和应用精度需求设置一组参数值便于进行分解结果的筛选。在不同的Ndir参数的控制下， 会得到多组分解结果。 根据本法发明优选的，所述S3中，正交性指数OI的计算公式为： 式(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分量极值的个数和过 10 CN 111597494 A 说　明　书 6/14 页 零点的个数之差的绝对值小于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输出的气候序列被表示为 K为IMF分量的总 数，ck(t)是第k个IMF，rk 1(t)是最终的残差分量，即趋势项。 根据发明优选的，步骤S2及S7-2中迭代停止条件为：对于(1-α)％的数据，α为0到1 之间的常数，δ(t)