高原气象  2018, Vol. 37 Issue (5): 1241-1253  DOI: 10.7522/j.issn.1000-0534.2018.00026
0

引用本文 [复制中英文]

任梅芳, 庞博, 徐宗学, 等. 2018. 基于随机森林模型的雅鲁藏布江流域气温降尺度研究[J]. 高原气象, 37(5): 1241-1253. DOI: 10.7522/j.issn.1000-0534.2018.00026
[复制中文]
Ren Meifang, Pang Bo, Xu Zongxue, et al. 2018. Downscaling of Air Temperature in the Yarlung Zangbo River Basin Based on the Random Forest Model[J]. Plateau Meteorology, 37(5): 1241-1253. DOI: 10.7522/j.issn.1000-0534.2018.00026.
[复制英文]

资助项目

国家自然科学基金项目(211700044,210200002)

通信作者

徐宗学(1962-), 男, 山东人, 教授, 主要从事气候变化对水循环的影响研究.E-mail:zxxu@bnu.edu.cn

作者简介

任梅芳(1987-), 女, 青海人, 博士研究生, 主要从事气候变化及城市洪涝灾害研究.E-mail:renmeifang@mail.bnu.edu.cn

文章历史

收稿日期: 2017-11-11
定稿日期: 2018-02-06
基于随机森林模型的雅鲁藏布江流域气温降尺度研究
任梅芳, 庞博, 徐宗学, 赵彦军     
北京师范大学 水科学研究院 水沙科学教育部重点实验室/城市水循环与海绵城市技术北京市重点实验室, 北京 100875
摘要: 采用随机森林RF(Random Forest)模型对雅鲁藏布江流域22个站点的日平均气温进行降尺度研究,为了探求在雅鲁藏布江流域更适宜的气温降尺度方法,采用多元线性回归MLR、人工神经网络ANN和支持向量机SVM三种方法作为对比模型,并且采用主成分分析PCA和偏相关分析PAR两种分析方法,进行特征变量筛选。采用纳西效率系数NASH、均方根误差RMSE系数、绝对误差MAE和相关系数r值四种标准来评价模型的模拟效果。结果表明,RF模型的模拟效果要明显优于其他几种方法的模拟结果;采用PAR筛选特征变量的模型计算结果,不仅优于采用PCA筛选特征变量模型的模拟结果,且较稳定,另外,各种模型验证期的NASH效率系数都在0.86以上,相关系数都在0.93以上,所用几种模型都能较好地模拟雅江流域平均气温。选取MPI-ESM-LR模式在未来(2016-2050年)两种极端典型浓度路径RCP(Representative Concentration Pathway)排放情景RCP2.6和RCP8.5下的试验数据,研究雅鲁藏布江流域未来气温变化趋势表明,雅鲁藏布江流域未来2016-2050年在RCP2.6和RCP8.5两种排放情景下,平均气温都呈现出持续上升的趋势,在RCP2.6排放情景下日平均气温平均上升0.14℃,在RCP8.5排放情景下日平均气温平均上升0.30℃。
关键词: 统计降尺度    随机森林    雅鲁藏布江    气温    
1 引言

目前, 气候模式GCMs(General Circulation Models)和降尺度技术是未来气候特征研究的主要手段。GCMs通过数值模拟技术输出全球气候状态, 是进行气候模拟和未来气候变化情景预估的主要工具(IPCC, 2007)。然而, 尽管GCMs的模拟结果能在全球尺度上进行宏观分析, 但其输出结果都是大网格尺度的空间数据, 空间分辨率较低, 很难满足区域尺度上水文过程的研究, 因此, 降尺度技术就成为了连接GCMs模式大尺度数据与区域小尺度之间气候特性的纽带(赵娜等, 2017; 刘文丰等, 2014)。降尺度技术一般分为统计降尺度、动力降尺度和动力与统计结合的混合降尺度研究方法, 而统计降尺度以其简单、计算量少、效果明显等特点成为了低分辨率资料到流域尺度资料转换的有效手段。

统计降尺度方法有转换函数法、天气分类法和天气发生器法。常用的统计降尺度方法有很多, 例如多元线性回归MLR(Multiple Linear Regression)、人工神经网络ANN(Artificial Neural Networks)、支持向量机SVM(Support Vector Machines)及典型相关分析CCA(Canonical Correlation Analysis)等。国内外关于不同统计方法在降尺度中效果对比的研究也比较广泛, Schoof et al(2001)指出在日均气温的降尺度研究中ANN的模拟效果要优于MLR; Duhan et al(2015)研究指出最小二乘法支持向量机LS-SVM(Least Square Support Vector Machine)对日最高与最低极端气温的降尺度模拟效果要优于MLR和ANN的模拟效果; Kostopoulou et al(2007)指出MLR和CCA在极端气温降尺度研究中的模拟效果要优于ANN。综上所述, 基于不同统计方法的降尺度研究中, 不同研究背景下不同方法的表现能力都是不同的。

大量研究表明, 筛选合适的特征变量是降尺度研究中最重要的步骤之一(Hewitson et al, 1996), 多种方法可用来筛选特征变量, 例如主成分分析PCA(Principal Component Analysis)、偏相关分析PAR(Partial Correlation Analysis)、最大协方差分析MCA(Maximum Covariance Analysis)、交互式模型拟合方法(Interactive Model Fitting Approaches)等。然而, Sharma(2000)筛选特征变量的方法在应用过程中也存在一些缺陷, 例如传统相关性分析方法在非线性关系及稳定性方面的缺陷, 如果特征变量与模型预测值之间存在一种线性的关系(特征变量与模型预测值之间存在某种线性函数), 通常情况下采用相关性系数, 但如果特征变量与模型预测值之间的关系相对较复杂, 而大多数气象因素之间的关系属于这一种类型, 传统的线性相关性分析很可能会导致误导性预测。因此, 选用合适的特征变量筛选方法在降尺度研究中是至关重要的。

随机森林RF(Random Forest)是一种由Breiman(2001)最早提出的基于统计学习理论的分类与回归算法, 主要是利用Bootstrap重新抽样的方法从原始数据中抽取多个样本, 再对每个样本构建决策树, 从而通过对所有决策树的预测进行组合和投票的方式得出最终结果。RF模型被广泛应用于各个研究领域, 例如风险分析(Malekipirbazari et al, 2015), 地表水研究(Baudron et al, 2013), 遥感影像分析(Tatsumi et al, 2015)和洪水风险评估(Wang et al, 2015)等, 大量的实例研究表明, RF具有较强的数据挖掘和预测能力, 对异常值具有较好的容忍度, 且不容易出现过度拟合, 甚至被称为当前最好的算法之一(赖成光等, 2015)。

采用随机森林模型, 对雅鲁藏布江流域(简称雅江流域)22个站点的日平均气温进行降尺度研究, 结果与多元线性回归MLR、人工神经网络ANN和支持向量机SVM三种模型进行对比, 寻求在雅江流域更适宜的气温降尺度方法。另外, 采用主成分分析PCA和偏相关分析PAR两种方法, 在建立MLR, ANN和SVM模型之前进行变量筛选, 分别表示为MLR-PCA, ANN-PCA, SVM-PCA和MLR-PAR, ANN-PAR, SVM-PAR, 而RF模型内置变量重要性评价功能, 故建立模型之前不需要进行因子筛选。

雅江流域因其特殊的地形、地理位置和气候条件, 受观测资料所限, 例如流域测站稀少、遥感监测数据时序短、不同来源的数据不匹配等, 在流域气候特征等方面的研究较少(尚可政等, 2010; 建军等, 2012)。本文通过对比与分析各种模型雅江流域的适用性, 探求雅江流域相对较适宜的气温降尺度方法, 以期为雅江流域气候变化研究提供技术支撑; 并选取MPI-ESM-LR模式在未来(2016-2050年)两种极端典型浓度路径(Representative Concentration Pathway, RCP)排放情景RCP2.6和RCP8.5下的试验数据, 对雅江流域未来平均气温进行模拟与分析, 以期为雅江流域气候变化研究以及为青藏高原气候变化适应性决策和水资源管理等提供背景数据。

2 研究区域概况与资料选取 2.1 研究区域概况

雅江流域是青藏高原上最大的河流, 也是最重要的国际河流之一, 发源于西藏自治区南部, 喜马拉雅山北麓的杰玛央宗冰川。雅江流域横贯青藏高原南部, 河流总长为2 229 km, 流域面积为2.42×105 km2, 年径流量为1.39×108 m3(You et al, 2007; Li et al, 2014)。雅江流域位于81°9′E 97°1′E, 29°9′N 31°9′N(图 1), 是世界上最高的河流之一, 平均海拔为4 600 m, 流域地势西高东低, 平均坡降为2.6 ‰, 其上较大支流有拉萨河、帕隆藏布、年楚河、多雄藏布和尼洋河等(Tang et al, 1998)。雅江流域的气候条件主要受高原地理位置和地势特点支配, 上游地区属高原寒温带气候, 具有气候寒冷、雨量稀少的特点; 中部地区属高原温带气候区, 具有长冬无夏、光照充足、蒸发量大及雨季短的气候特点; 下游属热带亚热带气候, 天气变化较为复杂, 季节现象明显; 流域多年平均气温的地域分布较明显, 西部气温较低, 东南部气温较高(刘文丰等, 2014)。

图 1 雅江流域位置和气象站点位置分布 Figure 1 Geographical location of the Yarlung Zangbo River basin and the distribution of meteorological stations
2.2 数据来源 2.2.1 站点气温数据

选取雅江流域22个气象站点的观测日尺度平均气温数据, 其中9个站点在雅江流域内, 13个站点在流域外周边区域。22个气象站点的详细信息如表 1所示, 其地理位置如图 1所示。气象数据来源于中国气象科学数据共享服务网(http://cdc.cma.gov.cn), 选取1978年1月至2015年12月的日平均气温数据, 由于观测气温数据序列中存在缺测的现象, 对观测气温数据进行了前期处理, 由同期多年当月平均值代替序列中的缺测数据。

表 1 气象站点基本信息 Table 1 Basic information of the meteorological stations
2.2.2 再分析数据

选用再分析数据NCEP(The National Centers for Environmental Prediction)作为气温降尺度研究中的特征变量, 数据来源于: http://www.cdc.noaa.gov采用双线性插值法将NCEP数据内插到每个站点上。

表 2 NCEP特征变量代号和定义 Table 2 NCEP characteristic variable codes and definition
2.2.3 GCM数据介绍

GCM数据选用IPCC第五次评估报告AR5的第五次耦合模式比较计划CMIP5气候模式的试验数据(数据来源于: http://esgf-node.llnl.gov/search/cmip5)。CMIP5提供四种典型浓度路径排放情景, 分别为RCP2.6, RCP4.5, RCP6.0和RCP8.5四种情景(常燕等, 2016; 张蓓等, 2017; 伍清等, 2017)。RCP2.6代表辐射强迫在2100年之前达到峰值, 到2100年下降到2.6 W·m-2, RCP8.5代表辐射强迫在2100年上升至8.5 W·m-2。研究表明MPI-ESM-LR模式在青藏高原具有较好的适应性(Su et al, 2013; Zhu et al, 2013; Xu et al, 2017), 能够相对较好地模拟青藏高原气候变化趋势, 因此选取MPI-ESM-LR模式在未来(2016-2050年)RCP2.6和RCP8.5两种排放情景下的试验数据, 对雅江流域未来平均气温进行模拟与分析。MPI-ESM-LR模式是德国Max-planck气象研究所开发的大尺度全球气候模型, 该数据为1.875°×1.875°大尺度网格数据, 本研究中采用双线性插值法将GCM数据内插到每个站点上。

3 方法介绍 3.1 随机森林

随机森林方法最早是由Breiman(2001)提出的一种由多棵分类与回归决策树CART组合构成的一种新型的机器学习算法。RF模型由众多CART决策树构成, 模型的决策能力取决于每一棵CART决策树, 每一棵决策树由根节点、子节点和叶子节点组成, 每个子节点包含一个评判规则, 叶子节点则对应一个评判级别。RF算法的步骤大致分为以下几步:首先从原始数据集中抽取训练集, 每个训练集数据数量约为原始数据集的2/3, 之后为每个训练集建立决策树, 众多决策树就构成了森林。设RF模型中共有M个风险指标, 随机抽取m个(mM)作为节点指标, Breiman推荐m=$\sqrt M $, 在这m个节点指标中选择基尼最小值作为该节点的分支标准, 最后根据所有决策树的预测结果, 采用投票的方式来决定新样本的类别。每次抽样未被抽中的1/3数据构成了袋外数据(out-of-bag, OOB), 利用这些袋外数据估计内部误差, 称之为袋外误差(error of out-of bag, EOOB)。其计算原理如下:

$ {E_{{\rm{OOB}}}} = \frac{1}{n}{\sum\limits_{i = 1}^n {\left[ {\hat Y\left({{X_i}} \right) - {Y_i}} \right]} ^2} $ (1)

式中: n为OOB样本个数; $\hat Y$(Xi)为根据给定样本Xi, RF模型的输出数据; Yi为观测数据。

除了分类与回归, RF模型还可以评价特征变量的重要性, 特征变量的重要性通过RF算法中的OOB误差进行估计。首先计算每个决策树的袋外误差EOOB1, 然后对风险指标i的数据随机加入噪声并计算袋外误差EOOB2, 风险指标i的重要性就可表示为:

$ V\left(i \right) = \frac{1}{N}\sum {\left({{E_{{\rm{OOB2}}}} - {E_{{\rm{OOB1}}}}} \right)\;\;, } $ (2)

改变指标i造成的袋外误差EOOB2越大, 精度减少的越多, 说明变量i就越重要(赖成光等, 2015; 马玥等, 2016; 康有等, 2014; 张雷等, 2014)。

3.2 多元线性回归

多元线性回归MLR(Multiple Linear Regression)是一种影响因变量Y的自变量不只是一个, 通常为p≥2时采用的回归模型。本研究中建立观测气温数据与NCEP特征变量之间的多元线性回归数学模型, 其数学表达式可表示为:

$ {Y_n} = {\beta _0} + {\beta _1}{x_{n1}} + {\beta _2}{x_{n2}} + \cdots + {\beta _p}{x_{np}} + {\varepsilon _n}\;\;, $ (3)

式中: x为自变量; y为因变量; β0为常量, β1, β2, …, βp是回归系数矩阵; p为统计变量的个数; n为实测数据的个数; εn为随机误差。

3.3 人工神经网络

人工神经网络是基于人脑的一种抽象、简化和模拟, 以大量的神经元广泛地连接而成的非线性信息处理系统。人工神经网络模型能够实现输入层和输出层之间多维非线性精确映射, 而不需要在两者之间提前建立任何数学计算模型(李云良等, 2015)。从结构上, 人工神经网络通常可以分为前馈性网络(Feedforward Neural Networks)与反馈性网络(Feedbackward Neural Networks)。本研究中采用反向传播人工神经网络(Back-Propagation Neural Networks, BP算法), 该算法属于典型的前馈性网络。BP神经网络其主要特点是信息的分布式储存和并行协同处理方式。典型的BP模型为具有输入层-隐含层-输出层三层网络的模型, 一个或多个隐含层在输入层和输出层之间, 模型通过信息的正向传播和误差的反向传播来进行计算, 并且通过激活函数、权重和偏差来增强模型解决非线性复杂问题的能力(高翔等, 2017)。BP模型分为训练与检验两个部分, 模型由率定数据集进行训练, 当达到最小误差限制或者最大迭代次数时, 模型训练结束。模型的数学表达由下式所示:

$ {O_k} = {g_2}\left[ {\sum\limits_{j = 1}^M {{w_{kj}}{g_1}\left({\sum\limits_{i = 1}^N {{w_{ji}}{x_i} + {w_{j0}}} } \right) + {w_{k0}}} } \right]\;\;\;, $ (4)

式中: xi为节点i的输入值; Ok为节点k的输出值; g1为隐含层的激活函数; g2为输出层的激活函数; MN分别是输入层和隐含层中神经元的个数; wj0为隐含层中第j个神经元的偏差; wk0为输出层中第k个神经元的偏差; wji为输入节点i与隐含节点j之间的权重; wkj为隐含节点j和输出节点k之间的权重。

3.4 支持向量机

支持向量机(Support Vector Machine, SVM)是最早由Vapnik(1997)提出的分类机器学习算法, 后来逐渐衍生到回归算法。SVM算法的基本思路是通过求解最大间隔的问题将输入空间数据映射到一个高维特征向量空间, 之后在高维特征空间中构造最优分类超平面进行线性分类。最大间隔问题是通过求解二次约束方程来实现, 具体详见参考文献(Jing et al, 2016)。

3.5 模型的建立

采用的MLR、ANN、SVM和RF模型进行计算是在MATLAB(2016)环境下进行计算的。为了能够更加系统地分析ANN模型的计算效果并取得最优的模拟结果, 本次计算将ANN模型的隐含层从1取到20, 训练ANN模型直到模型效率达到最小误差限制或者最大迭代次数, 本次计算将最小误差限制取为0.005, 将最大迭代次数取为1 000次, 学习率lr取为0.1, ANN模型隐含层和输出层的传递函数借鉴Matlab中BP神经网络工具的设计原则, 隐含层和输出层的传递函数分别使用Sigmoid函数和线性函数; 在SVM模型中, 本次计算选用线性核函数; RF算法中需要定义2个参数:决策树树木Ntree和决策树节点, 划分时随机选取的候选特征变量个数m, 本文通过多次试算研究发现, 当决策树数目Ntree≥500时, OOB误差趋于稳定, 据此, 本文中将决策树的树木Ntree取为500, 决策树节点划分时随机选取的候选变量m=$\sqrt M $进行分类, 其中M为RF模型中风险指标个数。

3.6 模型评价指标

选取4种模型评价指标来评价各个模型的表现能力:纳西效率系数(Nash efficiency coefficient, NASH)、均方根误差(Root Mean Squared Error, RMSE)、平均绝对误差(Mean Absolute Error, MAE)和相关系数(Coefficient of Correlation, r), 其计算公式如下所示:

$ NASH = 1 - \frac{{\sum\limits_{i = 1}^n {{{\left({{Y_{i, {\rm{obs}}}} - {Y_{i, {\rm{sim}}}}} \right)}^2}} }}{{\sum\limits_{i = 1}^n {{{\left({{Y_{i{\rm{, obs}}}} - \overline {{Y_{{\rm{obs}}}}} } \right)}^2}} }}\;\;, $ (8)
$ RMSE = \sqrt {\frac{{\sum\limits_{i = 1}^n {{{\left({{Y_{i, {\rm{obs}}}} - {Y_{i{\rm{, sim}}}}} \right)}^2}} }}{n}} \;\;, $ (9)
$ MAE = \frac{{\sum\limits_{i = 1}^n {\left| {{Y_{i, {\rm{obs}}}} - {Y_{i, {\rm{sim}}}}} \right|} }}{n}\;\;, $ (10)
$ \begin{array}{l} r = \\ \frac{{\sum\limits_{i = 1}^n {\left[ {\left({{Y_{i, {\rm{obs}}}} - \overline {{Y_{{\rm{obs}}}}} } \right)\left({{Y_{i.{\rm{sim}}}} - \overline {{Y_{{\rm{sim}}}}} } \right)} \right]} }}{{\sqrt {\left[ {\sum\limits_{i = 1}^n {{{\left({{Y_{i, {\rm{obs}}}} - \overline {{Y_{{\rm{obs}}}}} } \right)}^2}} } \right]} \sqrt {\left[ {\sum\limits_{i = 1}^n {{{\left({{Y_{i, {\rm{sim}}}} - \overline {{Y_{{\rm{sim}}}}} } \right)}^2}} } \right]} }} \end{array} $ (11)

式中: Yi, obs为第i个站点的实测气温数据; Yi, sim为第i个站点的模拟气温数据; $\overline {{Y_{{\rm{obs}}}}} $为所有站点实测气温数据的平均值; $\overline {{Y_{{\rm{sim}}}}} $为所有站点模拟气温数据的平均值。

3.7 Mann-Kendall非参数检验法

Mann-Kendall(M-K)趋势检验用于检验数据集变化趋势的显著性, 在M-K检验分析时, 假设H0表示数据集的数据是样本独立同时分布, 不存在趋势, 假设H1是数据集中存在的一个单调趋势。M-K检验统计变量S计算如下所式(徐宗学等, 2006):

$ Z = \left\{ \begin{array}{l} \frac{{S - 1}}{{\sqrt {{\rm{var}}\left(S \right)} }}, \;\;S > 0\\ \;\;\;\;\;\;\;0, \;\;\;\;\;\;S = 0\\ \frac{{S + 1}}{{\sqrt {{\rm{var}}\left(S \right)} }}, \;\;S < 0 \end{array} \right\}\;\;\;, $ (12)
$ S = \sum\limits_{i = 1}^{n - 1} {\sum\limits_{k = i + 1}^n {{\rm{sgn}}\left({{x_k} - {x_i}} \right)\;\;\;, } } $ (13)
$ {\rm{sgn}}\left({{x_k} - {x_i}} \right) = \left\{ \begin{array}{l} \;\;1, \;\;\;\;\;{\mathit{x}_k} - {x_i} > 0\\ \;\;0, \;\;\;\;{x_k} - {x_i} = 0\\ \; - 1, \;\;\;{x_k} - {x_i} < 0 \end{array} \right\}\;\;\;, $ (14)

式中: n为数据集的长度; xkxi为样本数据集。Kendall倾斜度β可用来量化单调趋势, β表达式如下:

$ \beta = {\rm{Median}}\left({\frac{{{x_i} - {x_j}}}{{i - j}}} \right), \;\;\forall j < i, \;1 < j < i < n $ (15)

式中:当β < 0时表示数据集下降的趋势, β > 0表示数据集上升的趋势。

4 结果与讨论 4.1 因子筛选结果 4.1.1 RF模型

RF模型内置变量重要性评价功能, 模型可对每一个特征变量的相对重要性做出评估, 这将有利于RF模型在降尺度研究中特征变量的筛选。根据RF计算结果, 特征变量在每个站点的重要性排名都不相同, 这可能与每个站点的地理位置及气象与地理因素有关, 以普兰站和左贡站为流域上下游代表站, 以拉萨站和当雄站为流域中部代表站, 图 2为代表站各特征变量的因子重要性排名示意图。

图 2 代表站特征变量重要性排名 Figure 2 The significance rank of characteristic variables at representative stations

图 3为RF模型对22个站点特征变量的重要性评估, 图中最下边第一条线代表 22个站点中该变量对于平均气温重要性排名数据中的最高排名, 最上边一条线代表 22个站点中该变量对于平均气温重要性排名数据中的最低排名; “□”下边界代表 22个站点中该变量对于平均气温重要性排名数据的第一四分位点, 上边界代表 22个站点中该变量对于平均气温重要性排名数据的第三四分位点, 中间的线代表 22个站点中该变量对于平均气温重要性排名数据的中值(部分有重叠)。从图 3中可以看出, 850 hPa高度经向速度分量p8_v为影响最显著的因子; 排名第二显著的变量为850 hPa高度纬向速度分量p8_u; 地表相对湿度rhum, 500 hPa高度经向速度分量p5_v, 500 hPa高度位势高度p500, 500 hPa高度纬向速度分量p5_u和850 hPa高度位势高度p850分别排在第3~7位; 而地表气温temp和地表比湿shum为大部分站点中重要性相对较差的两个变量。

图 3 22个站点特征变量重要性排名 Figure 3 The significance rank of characteristic variables at 22 stations

由于RF模型中的袋外误差EOOB可以对RF模型的模拟效果进行偏差估计, 图 4给出了22个站点考虑1~10个特征变量的袋外误差EOOB分布, 图中最下边第一条线代表 22个站点中分别考虑1~10个变量时RF模型袋外误差数据中的最小值, 最上边一条线代表 22个站点中分别考虑1~10个变量时RF模型袋外误差数据中的最大值; “□”下边界代表 22个站点中分别考虑1~10个变量时RF模型袋外误差数据中的第一四分位点, 上边界代表 22个站点中分别考虑1~10个变量时RF模型袋外误差数据中的第三四分位点, 中间的线代表 22个站点中分别考虑1~10个变量时RF模型袋外误差数据中的中值。由图 4可以看出, 随着预报因子的增加, EOOB误差逐渐减小, 这就说明RF模型能够避免过度拟合的现象。当所考虑的特征变量达到某一固定值时, EOOB误差将趋于稳定, 从而RF模型的表现能力也就趋于相对稳定, 本次研究结果显示, 当特征变量个数达到7个时, RF模型的表现能力趋于稳定。考虑到RF模型强大的计算功能, 也为了能够考虑更多的信息, 本研究中将10个统计变量全部作为研究特征变量因子。

图 4 基于RF模型的袋外误差 Figure 4 The EOOB based on the RF model
4.1.2 主成分分析与偏相关分析变量筛选方法

(1) 主成分分析

主成分分析(Principal Component Analysis, PCA)最早是由Pearson于1901年提出的一种通过对协方差矩阵进行分析, 主要目的是在降低数据维数的条件下保持数据集对方差贡献最大的一种统计分析方法(王莺等, 2014)。其原理为通过降维的方法将原来的多个指标转化为一个或几个综合指标, 即主成分(PCs), 而各主成分之间互不相关。采用主成分分析方法PCA进行因子筛选时, 本文选取累计方差贡献率高于90 %的主成分作为统计变量。由于论文篇幅所限, 以安多站为代表站, 表 3给出代表站各主成分的贡献率及累计贡献率, 安多站所选主成分为1~5主成分。

表 3 安多站主成分累计贡献率及所选主成分 Table 3 The cumulative contribution and choices of PCs at Amdo station

(2) 偏相关分析

偏相关分析(Partial correlation analysis, PAR)是在多要素构成的自然系统中, 对其他变量的影响进行控制的条件下, 衡量多个变量中某两个变量之间的线性相关关系的密切程度的分析方法, 以偏相关系数来度量偏相关程度。大量研究表明, 多个统计变量的组合在降尺度研究中的表现能力要优于单个变量, 而在某一个给定的站点上, 降尺度研究的结果基本上是依赖于偏相关系数最高的前4~7个因子(Yang et al, 2016; Tatsumi et al, 2015)。为了能够考虑更多的信息, 本研究选用具有最高偏相关系数的前7个因子作为研究特征变量。图 5为雅江流域22个站点NCEP特征变量与各站点平均气温的偏相关系数整体分布, 图中最下边第一条线代表 22个站点中该变量对于平均气温偏相关系数的最小值; 最上边一条线代表 22个站点中该变量对于平均气温偏相关系数的最大值; “□”下边界代表 22个站点中该变量对于平均气温偏相关系数的第一四分位点, 上边界代表 22个站点中该变量对于平均气温偏相关系数的第三四分位点, “□”中间的线代表 22个站点中该变量对于平均气温偏相关系数的中值。总的来说, 500 hPa高度位势高度p500具有最高的偏相关系数, 850 hPa高度位势高度p850的偏相关系数排在第二位; 平均海平面气压mslp, 地表相对湿度rhum, 850 hPa高度经向速度分量p8_v, 850 hPa高度纬向速度分量p8_u和500 hPa高度纬向速度分量p5_u分别排在第3~7位; 而500 hPa高度经向速度分量p5_v和地表气温temp为大部分站点中偏相关系数最小的两个因子。

图 5 22个站点因子偏相关系数 Figure 5 The PAR correlation coefficients of variables at 22 stations
4.2 各种模型之间对比

将数据集中1978年1月1日至2002年12月31日25年数据作为模型的率定期, 将2003年1月1日至2015年12月31日数据作为模型的验证期。从各种模型计算结果(表 4)可以看出, RF模型的NASH效率系数为0.97, 高于MLR, ANN和SVM三种方法的NASH效率系数(0.87~0.93), RF的RMSE系数为1.20, 低于MLR, ANN和SVM三种方法的RMSE系数(1.87~2.44), RF的MAE系数为0.87, 低于MLR, ANN和SVM三种方法的MAE系数(1.39~1.86), RF的r系数为0.99, 高于MLR, ANN和SVM三种方法的r系数(0.94~0.96)。RF的NASH效率系数为0.92, 高于MLR, ANN和SVM三种方法的NASH效率系数(0.86~0.90), RF的RMSE系数为1.88, 低于MLR, ANN和SVM三种方法的RMSE系数(2.01~2.46), RF的MAE系数为1.45, 低于MLR, ANN和SVM三种方法的MAE系数(1.55~1.91), RF的r系数为0.96, 高于MLR, ANN和SVM三种方法的r系数(0.93~0.96)。

表 4 率定期与验证期模型效果 Table 4 Model performance in calibration and validation periods

图 6为各种模型率定期和验证期的NASH效率系数分布, 图中最下边第一条线代表 22个站点中不同方法NASH效率系数的最小值, 最上边一条线代表 22个站点中不同方法NASH效率系数的最大值; “□”下边界代表 22个站点中不同方法NASH效率系数中第一四分位点; 上边界代表 22个站点中不同方法NASH效率系数中的第三四分位点; 中间的线代表 22个站点中不同方法NASH效率系数的中值。从图 6中可以看出, RF模型的模拟效果要明显优于其余几种方法的模拟结果; 而采用不同的特征变量筛选方法对模拟结果有较显著的影响, 采用PAR筛选特征变量的模型计算结果, 不仅优于采用PCA筛选特征变量模型的模拟结果, 且采用PAR筛选特征变量的模型效果要比采用PCA筛选特征变量模型的模拟结果更加稳定; 另外, 各种模型验证期的NASH效率系数都在0.86以上, 相关系数都在0.93以上, 据此, 以上几种模型都能较好地模拟雅江流域平均气温。

图 6 各种模型率定期与验证期NASH效率系数分布 Figure 6 The spatial distribution of NASH in calibration periods and the validation periods with different models

图 7为RF模型与其他模型的NASH效率系数的空间分布对比, 由于采用PAR筛选特征变量的模型的效果要优于采用PCA筛选特征变量模型的模拟结果, 本次仅对RF模型和各个模型采用PAR筛选特征变量模型的效果进行对比分析。从NASH效率系数空间分布来看, 在雅江流域22个站点上, RF模型的表现效果最好, 其次为ANN模型, 而SVM模型和MLR模型的NASH效率系数相对较低。

图 7 模型NASH效率系数空间分布 Figure 7 The spatial distribution of NASH with different models
4.3 未来气候情景分析

为了探求雅江流域在21世纪中期的气候情况, 将MPI-ESM-LR模式在RCP2.6和RCP8.5两种极端排放情景下在2016-2050年的特征变量试验数据, 导入已建好的RF模型中, 从而对雅江流域未来气温进行模拟分析。图 8为雅江流域日平均气温在各站点的变化情况空间分布。

图 8 雅江流域2016-2050年未来气温升高空间分布(单位: ℃) Figure 8 Spatial distribution of future temperature increased in the Yarlung Zangbo river basin from 2016 to 2050.Unit: ℃

研究结果表明, 雅江流域未来2016-2050年在RCP2.6排放情景下日平均气温平均上升0.14 ℃, 在RCP8.5排放情景下日平均气温平均上升0.30 ℃; 从图 8中可以看出, 雅江流域虽然在RCP2.6排放情景下的气温上升趋势不是十分显著, 但在未来2016-2050年流域在RCP2.6和RCP8.5两种排放情景下, 都表现出气温持续上升的趋势, 在RCP8.5排放情景下的气温升高值明显高于在RCP2.6排放情景下的升高值, 这也符合两种排放情景的模式设置。

本文对雅江流域未来2016-2050年日平均气温的模拟值进行了M-K趋势检验分析, 并对未来时间序列进行逐月趋势检验分析, 全年和逐月趋势检验分析结果如表 5所示, 图 9为雅江流域未来2016-2050年日平均气温Kendall倾斜度β插值图。

图 9 雅江流域未来2016-2050年日平均气温Kendall倾斜度β插值图(单位: ℃·a-1) Figure 9 The Kendall Test β contour map of the daily average temperature in the Yarlung Zangbo River basin from 2016 to 2050.Unit: ℃·a-1

表 5可以看出, 雅江流域未来2016-2050年在RCP2.6和RCP8.5两种排放情景下, 全年M-K检验的Z值分别为4.362和10.85, β值分别为0.004和0.010 ℃·a-1, 雅江流域未来2016-2050年在两种排放情景下都呈现出增温趋势; RCP8.5排放情景下的增温趋势高于RCP2.6排放情景; 对于RCP2.6排放情景, 全年中的1月、2月和3月的增温趋势高于其他月份; 对于RCP8.5排放情景, 4月、5月和6月的增温趋势相对较高。从图 9中可以看出, 雅江流域未来2016-2050年在RCP2.6排放情景下, 上、下游两端的增温趋势高于流域中部区域, 在RCP8.5排放情景下, 增温趋势呈现出从北到南依次降低。

表 5 研究期全年与逐月日平均温度变化趋势 Table 5 Annual and seasonal trends for average temperature in study periods
5 结论

基于雅鲁藏布江流域内与周边共22个气象站点的观测日平均气温, 采用NCEP再分析数据, 应用随机森林RF模型对雅江流域22个站点的日平均气温进行降尺度研究, 结果与多元线性回归MLR、人工神经网络ANN和支持向量机SVM三种方法进行对比分析, 结合主成分分析PCA和偏相关分析PAR两种筛选因子方法, 探求雅江流域更适宜的气温降尺度方法。采用纳西效率系数NASH、均方根误差RMSE系数、绝对误差MAE和相关系数r值四种标准来评价模型的模拟效果。得到以下主要结论:

(1) 根据不同模型的计算结果, RF模型的模拟效果要明显优于其余几种方法的模拟结果; 采用PAR筛选特征变量的模型计算结果, 不仅优于采用PCA筛选特征变量模型的模拟结果, 且更加稳定; 另外, 各种模型验证期的NASH效率系数都在0.86以上, 相关系数都在0.93以上, 据此可以推断, 以上几种模型都能较好地模拟雅江流域平均气温。

(2) 研究结果表明, 雅江流域未来2016-2050年在RCP2.6和RCP8.5两种排放情景下, 都表现出温度持续上升的趋势, 在RCP2.6排放情景下日平均气温平均上升0.14 ℃, 在RCP8.5排放情景下日平均气温平均上升0.30 ℃。

参考文献
Baudron P, Alonso-Sarría F, García-Aróstegui J L, et al. 2013. Identifying the origin of groundwater samples in a multi-layer aquifer system with Random Forest classification[J]. J Hydrology, 499: 303–315. DOI:10.1016/j.jhydrol.2013.07.009
Breiman L. 2001. Random Forest[J]. Machine Learning, 45(1): 5–32.
Duhan D, Pandey A. 2015. Statistical downscaling of temperature using three techniques in the Tons River basin in Central India[J]. Theor Appl Climatol, 121(3-4): 605–622. DOI:10.1007/s00704-014-1253-5
Hewitson B C, Crane R G. 1996. Climate downscaling:Techniques and application[J]. Climate Res, 7(2): 85–95.
IPCC. 2007. Contribution of working groups Ⅰ, Ⅱ, and Ⅲ to the fourth assessment report of the intergovernment panel on climat change[M]. Cambridge and New York: Cambridge University Press.
IPCC. 2013. Summary for policymaker of climate change 2013:The physical science basis.Contribution of working group 1 to the fifth assessment report of the intergovernment panel on climate change[M]. Cambridge University Press: Cambridge, UK.
Jing W L, Yang Y P, Yue X F, et al. 2016. A comparison of different regression algorithms for downscaling monthly satellite-based precipitation over north China[J]. Remote Sens, 8(10): 1–17.
Kharin V V, Zwiers F W, Zhang X, et al. 2013. Changes in temperature and precipitation extremes in the CMIP5 ensemble[J]. Climate Change, 119: 345–357. DOI:10.1007/s10584-013-0705-8
Kostopoulou E, Giannakopoulos C, Anagnostopoul C, et al. 2007. Simulation maximum and minimum temperature over Greece:A comparison of three downscaling techniques[J]. Theor Appl Climatol, 90(1/2): 65–82.
Li F P, Zhang Y Q, Xu Z X, et al. 2014. Runoff predictions in ungagged catchments in southeast Tibetan Plateau[J]. J Hydrology, 511(511): 28–38.
Liu T, Chen B D. 2000. Climate warming in the Tibetan Plateau during recent decades[J]. Climate, 20: 1729–1742.
Malekipirbazari M, Aksakalli V. 2015. Risk assessment in social lending via random forest[J]. Expert Systems with Applications, 42(10): 4621–4631. DOI:10.1016/j.eswa.2015.02.001
Schoof J T, Pryor S C. 2001. Downscaling temperature and precipitation:a comparison of regression-based methods and artificial neural networks[J]. Int J Climatol, 21970: 773–790.
Screen J A. 2014. Arctic amplification decreases temperature variance in Northern mid-to high-latitudes[J]. Nature Climate Change, 4(7): 577–582. DOI:10.1038/nclimate2268
Sharma A. 2000. Seasonal to interannual rainfall probabilistic forecasts for improved water supply management:Part1-A strategy for system predictor identification[J]. J Hydrology, 239(1-4): 232–239. DOI:10.1016/S0022-1694(00)00346-2
Su F G, Duan X L, Chen D L, et al. 2013. Evaluation of the global climate models in the CMIP5 over the Tibetan Plateau[J]. J Climate, 26(10): 3187–3208. DOI:10.1175/JCLI-D-12-00321.1
Tang Q C, Xiong Y. 1998. River hydrology in China[M]. Beijing: Science Press.
Tatsumi K, Oizumi T, Yamashiki Y. 2015. Effects of climate change on daily minimum and maximum temperatures and cloudiness in the Shikoku region:a statistical downscaling model approach[J]. Theor Appl Climatol, 120: 87–98. DOI:10.1007/s00704-014-1152-9
Vapnik V N. 1997. The nature of statistical learning theory[M]. New York: Springer-Verlag.
Wang Z, Lai C, Chen X, et al. 2015. Flood hazard risk assessment model based on random forest[J]. J Hydrology, 527: 1130–1141. DOI:10.1016/j.jhydrol.2015.06.008
Xu J W, Gao Y H, Chen D L, et al. 2017. Evaluation of global climate models for downscaling applications centred over the Tibetan Plateau[J]. Int J Climatol, 37(2): 657–671. DOI:10.1002/joc.4731
Yang C L, Wang N L, Wang S J, et al. 2016. Performance comparison of three predictor selection methods for statistical downscaling of daily precipitation[J]. Theor Appl Climatol: 1–12. DOI:10.1007/s00704-016-1956-x
You Q L, Kang S C, Wu Y H, et al. 2007. Climate change over the Yarlung Zangbo River Basin during 1961-2005[J]. J Geograp Sci, 17(4): 409–420. DOI:10.1007/s11442-007-0409-y
Zhu X H, Wang W Q, Fraedrich K. 2013. Future climate in the Tibetan plateau from a statistical regional climate model[J]. J Climate, 26(4): 10125–10138.
常燕, 吕世华, 罗斯琼, 等. 2016. CMIP5耦合模式对青藏高原动土变化的模拟和预估[J]. 高原气象, 35(5): 1157–1168. Chang Y, Lü S H, Luo S Q, et al. 2016. Evaluation and projections of permafrost on the Qinghai-Xizang Plateau by CMIP5 coupled climate models[J]. Plateau Meteor, 35(5): 1157–1168. DOI:10.7522/j.issn.1000-0534.2015.00090
高翔, 毛献忠. 2017. 基于BP神经网络方法计算河网区面源负荷[J]. 给水排水, 43(1): 156–160. Gao X, Mao X Z. 2017. Study on non-point source load in river network based on the BP neural network method[J]. Water & Wastewater Engineering, 43(1): 156–160. DOI:10.3969/j.issn.1002-8471.2017.01.032
建军, 杨志刚, 卓嘎. 2012. 近30年西藏汛期强降水事件的时空变化特征[J]. 高原气象, 31(2): 380–386. Jian J, Yang Z G, Zhuo G. 2012. Temporal-spatial distributions of severe precipitation event during flood period in Tibet in last 30 years[J]. Plateau Meteor, 31(2): 380–386.
康有, 陈元芳, 顾圣华, 等. 2014. 基于随机森林的区域水资源可持续利用评价[J]. 水电能源科学, 32(3): 34–38. Kang Y, Chen Y F, Gu S H, et al. 2014. Assessment of sustainable utilization of regional water resources based on random forest[J]. Water Resource and Power, 32(3): 34–38.
赖成光, 陈晓宏, 赵仕威, 等. 2015. 基于随机森林的洪灾风险评价模型及其应用[J]. 水利学报, 46(1): 58–66. Lai C G, Chen X H, Zhao S W, et al. 2015. A flood risk assessment model based on random forest and its application[J]. Journal of Hydraulic Engineering, 46(1): 58–66.
李云良, 张奇, 李淼, 等. 2015. 基于BP神经网络的鄱阳湖水位模拟[J]. 长江流域资源与环境, 24(2): 233–240. Li Y L, Zhang Q, Li M, et al. 2015. Using BP neural networks for water level simulation in Poyang lake[J]. Resources and Environment in the Yangtze Basin, 24(2): 233–240. DOI:10.11870/cjlyzyyhj201502008
刘文丰, 徐宗学, 李发鹏, 等. 2014. 基于ASD统计降尺度的雅鲁藏布江流域未来气候变化情景[J]. 高原气象, 33(1): 26–36. Liu W F, Xu Z X, Li F P, et al. 2014. Climate change scenarios in the Yarlung Zangbo river basin based on ASD model[J]. Plateau Meteor, 33(1): 26–36. DOI:10.7522/j.issn.1000-0534.2012.00176
陆晴, 吴绍洪, 赵东升. 2017. 1982-2013年青藏高原高寒草地覆盖变化及与气候之间的关系[J]. 地理科学, 37(2): 292–300. Lu Q, Wu S H, Zhao D S. 2017. Variations in alpine grassland cover and its correlation with climate variables on the Qinghai-Tibetan in 1982-2013[J]. Sci Geograp Sinica, 37(2): 292–300.
马玥, 姜琦刚, 孟治国, 等. 2016. 基于随机森林算法的农耕区土地利用分类研究[J]. 农业机械学报, 47(1): 297–303. Ma Y, Jiang Q G, Meng Z G, et al. 2016. Classification of land use in farming area based on random forest algorithm[J]. Transactions of the Chinese Society for Agricultural Machinery, 47(1): 297–303.
潘保田, 李吉均. 1996. 青藏高原:全球气候变化的驱动机与放大器:Ⅲ.青藏高原隆起对气候变化的影响[J]. 兰州大学学报:自然科学版, 32(1): 108–115. Pan B T, Li J J. 1996. Qing-hai-Tibetan plateau:A driver and amplifier of the global climate change:Ⅲ.The effect of the uplift of Qinghai-Tibetan plateau on climatic changes[J]. Journal of Lanzhou University (Natural Sciences), 32(1): 108–115.
任国玉, 徐铭志, 初子莺, 等. 2005. 近54年中国地面气温变化[J]. 气候与环境研究, 10: 717–727. Ren G Y, Xu M Z, Chu Z Y, et al. 2005. Changes of surface air temperature in China during 1951-2004[J]. Climatic Environ Res, 10: 717–727.
尚可政, 周海, 陈录元, 等. 2010. 西藏南部河谷气候变化趋势及预测-以江孜站为例[J]. 高原气象, 29(5): 1111–1118. Shang K Z, Zhou H, Chen L Y, et al. 2010. Trend and Prediction of climate change in southern Tibet valley-taking Gyangze meteorological station for example[J]. Plateau Meteor, 29(5): 1111–1118.
王莺, 王静, 姚玉璧, 等. 2014. 基于主成分分析的中国南方干旱脆弱性评价[J]. 生态环境学报, 23(12): 1897–1904. Wang Y, Wang J, Yao Y B, et al. 2014. Evaluation of drought vulnerability in southern china based on principle component analysis[J]. Ecology Environ Sci, 23(12): 1897–1904. DOI:10.3969/j.issn.1674-5906.2014.12.003
伍清, 蒋兴文, 谢洁. 2017. CMIP5模式对西南地区气温的模拟能力评估[J]. 高原气象, 36(2): 358–370. Wu Q, Jiang X W, Xie J. 2017. Evaluation of surface air temperature in Southwestern China simulated by the CMIP5 models[J]. Plateau Meteor, 36(2): 358–370. DOI:10.7522/j.issn.1000-0534.2016.00046
徐宗学, 张楠. 2006. 黄河流域近50年降水变化趋势分析[J]. 地理研究, 25(1): 27–34. Xu Z X, Zhang N. 2006. Long term trend of precipitation in the Yellow River basin during the past 50 years[J]. Geograp Res, 25(1): 27–34.
张蓓, 戴新刚. 2017. 基于CMIP5的2006-2015年中国气温预估偏差分析及订正[J]. 高原气象, 36(6): 1619–1629. Zhang B, Dai X G. 2017. Evaluation and correction for the deviation of the surface air temperature based on 24 CMIP5 models over China for 2006-2015[J]. Plateau Meteor, 36(6): 1619–1629. DOI:10.7522/j.issn.1000-0534.2016.00136
张雷, 王琳琳, 张旭东, 等. 2014. 随机森林算法基本思想及其在生态学中的应用-以云南松分布模拟为例[J]. 生态学报, 34(3): 650–659. Zhang L, Wang L L, Zhang X D, et al. 2014. The basic principle of random forest and its applications in ecology:A case study of Pinus yunnanensis[J]. Acta Ecologica Sinica, 34(3): 650–659.
赵娜, 岳天祥, 史文娇, 等. 2017. 基于HASM方法对气候模式气温降水的降尺度研究-以黑河流域为例[J]. 中国沙漠, 37(6): 1–10. Zhao N, Yue T X, Shi W J, et al. 2017. Downscaling simulation of Annual Average Temperature and Precipitation of CMIP5 outputs by using HASM-A case study in Heihe river basin[J]. J Desert Res, 37(6): 1–10.
Downscaling of Air Temperature in the Yarlung Zangbo River Basin Based on the Random Forest Model
REN Meifang , PANG Bo , XU Zongxue , ZHAO Yanjun     
Key Laboratory for Water and Sediment Science, Ministry of Education, College of Water Sciences, Beijing Normal University/Beijing Key Laboratory of Urban Hydrological Cycle and Sponge City Technology, Beijing 100875, China
Abstract: Random forest (RF) model was used to downscale the daily air temperature at 20 meteorological stations in and around the Yarlung Zangbo River basin.For the purpose to explore the better downscaling method of air temperature in the study area, three methods were used to compare, namely, the Multiple Linear Regression (MLR), Artificial Neural Network (ANN), and Support Vector Machine (SVM).Principal Component Analysis (PCA) and Partial Correlation Analysis (PAR) were used to select the characteristic variables.The model performance was assessed using four criteria, namely, the NASH coefficient of efficiency (NASH), the root mean squared error (RMSE), the mean absolute error (MAE), and the coefficient of correlation (r).The results showed that the performance of RF model was obviously better than other models; the results obtained by PAR to select characteristic variables were not only better than those used by PCA method, but also more stable.In addition, the NASH efficiency coefficients of various models were all above 0.86 and the correlation coefficients were all above 0.93 in the validation periods.Therefore, all of models used in this study can well simulate the average temperature in the Yarlung Zangbo River basin.The experimental data of two typical extreme concentration paths(Representative Concentration Pathway, RCP)emission scenarios RCP2.6 and RCP8.5 of MPI-ESM-LR model in the future (2016-2050)were chosen to study the future trend of temperature in the Yarlung Zangbo River basin.The results showed increasing trend of daily air temperature both under RCP2.6 and RCP8.5 scenarios in the future years of 2016-2050 in the Yarlung Zangbo River basin.The daily air temperature will increase by 0.14℃ under RCP2.6 scenario, and 0.30℃ under RCP8.5 scenario from 2016 to 2050.
Key words: Statistical downscaling    random forest    Yarlung Zangbo River basin    temperature