在水循环中, 降水是最活跃的因素^[1-2]。作为地表水资源的补给来源, 降水与人类生产生活及生态息息相关, 对于区域资源的时空分布、生态环境形成与演变及农业生产等起着决定性的作用^[3-4]。关于降水演变特征的分析早已成为国内外学者关注的焦点, 降水的变化趋势研究也已成为水文及气候系统研究中的主要组成部分^[5-9]。这些研究有利于增强对未来气候变化认识, 也有利于剖析区域水文、自然灾害及其生态环境的影响因素。由于降水受到多方面自然因素及人类生产活动的共同影响, 其成因及受影响因素具有不确定性, 对降水的研究仍有很大的空间, 仍需不断的研究与探索^[10-12]。

“三北工程建设是同我国改革开放一起实施的重大生态工程, 是生态文明建设的一个重要标志性工程。经过40a不懈努力, 工程建设取得巨大生态、经济、社会效益, 成为全球生态治理的成功典范”^[13]。由于三北工程区大部分区域地处西北干旱半干旱地区, 降水量不足、蒸发量过大, 水资源一直以来都是这里的稀缺资源。工程早期植被建设由于忽略了水资源承载力, 致使出现了“造林成活率低、生长缓慢、停滞甚至枯死等衰退现象”^[14-15]。加上西部大开发和一带一路战略开发的陆续实施, 三北工程区水资源供需矛盾更加突出。大气降水是地表水资源的补给来源, 对降水量演变特征及其变化趋势的研究有利于增强对未来气候变化认识并采取有效应对措施。而目前关于三北工程区多年来的降水量及其变化的研究相对较为薄弱。尽管王强在分析植被覆被变化对气候变化的响应时, 采用1982—2006年气象数据, 对三北工程区25a的气温降水变化作了研究, 但时间序列过短, 对于工程建设前后降水量变化的阶段特征缺乏对比^[16]；王鹏涛等选取117个站点1960—2011年的气象数据, 采用气候倾向率、距平分析法及空间插值等方法, 分析三北工程区近52a降水量的总体变化趋势和降水量距平空间变化, 但未揭示出降水量变化趋势的空间分异及其突变现象^[17], 未能全面揭示三北工程区的降水变化特征。

本文在三北工程区、重点建设区、重点县3个空间尺度上, 基于1951—2018年长时间序列降水量观测数据, 采用Mann-Kendall非参数检验、ARIMA模型预测和ArcGIS空间分析等方法, 开展年降水量、生长季降水量和非生长季降水量多尺度变化趋势与突变研究, 并对未来30年年降水量预测分析, 可为落实“以水定林草”的发展理念, 和因地制宜发展“雨养林草植被”提供科学支撑。

1 研究区域概况

早在1978年, 国家就决定在西北、华北及东北风沙危害和水土流失严重的地区, 建设大型防护林工程, 即三北防护林体系建设工程^[18]。三北工程区大部分地处西北干旱半干旱的水资源稀缺地区, 年均降水量在400 mm以下, 形成了“十年九旱, 不旱则涝”的气候特点, 曾严重制约着区域经济和社会发展^[12]。经过40年的建设, 林草植被初步得到恢复, 沙化土地和水土流失治理成效明显。随着造林面积及工农业和生活用水量的增加, 水资源供需矛盾日趋严重。树立“以水定林草”的发展理念, 因地制宜发展“雨养林草植被”, 是三北工程未来的建设目标和方向。工程区范围及18个重点建设区和5个重点县分布如图 1所示。

2 研究方法 2.1 数据来源

本研究数据源包括:中国科学院资源环境科学数据中心(http://www.resdc.cn/Default.aspx)739个站点1951—2018年共68a的逐日降水资料, 数据格式为文本文档。从统计的角度, 序列长度为68a的降水资料用于分析区域降水特性是比较可靠的。三北工程区和重点建设区边界来源于国家林业和草原局调查规划设计院项目共享数据, 数据格式为1：100万县域边界矢量数据。

2.2 研究方法

本研究主要采用Mann-Kendall趋势检验、ARIMA模型预测和ArcGIS空间分析等方法, 开展年降水量、生长季降水量和非生长季降水量多尺度变化趋势、突变分析与预测研究。

(1) Mann-Kendall趋势检验。Mann-Kendall方法由Mann和Kendall提出, “是一种基于秩的非参数统计检验方法, 不需要样本遵从一定的分布, 也不受少数异常值的干扰, 适用性强, 是时间序列趋势分析有效方法之一, 对揭示整体时间序列演变趋势与突变情况有良好的表现”^[19-21]。由世界气象组织(WMO)推荐并已广泛应用于气温、降水等要素时间序列的变化趋势分析。

1) 非参数Mann-Kendall趋势检验定义检验统计量S为：

(1)

式中：“i≠j, 且i, j≤n, Sign()为符号函数。当D_i－D_j小于、等于或大于零时, Sign(D_i－D_j)分别为-1、0或1”。当为长时间序列时(n＞10), 统计量Z为：

(2)

“当Z>0时, 时间序列呈增加趋势；Z < 0时, 时间序列为减少趋势。当Z的绝对值≥1.28、1.64、2.32时, 表示判别结果分别通过了信度为90%、95%、99%的显著性检验”^[20]。

2) Mann-Kendall突变检测

设有时间序列：D₂, D₃, …, D_n, 构造一秩序列r_i, r_i表示D_i>D_j(1≤j≤i)的样本累积数。定义S_k：

(3)

式中当D_i>D_j时, r_i=1；当D_i≤D_j时, r_i=0 (j=1, 2, …, i)。S_k的期望值E(S_k)及其序列方差Var(S_k)分别由以下两式定义：

(4)

(5)

假定数据序列具有独立性, 定义统计量UF_k：

(6)

式中“UF_k遵从标准正态分布, 给定一个显著性水平α, 查标准正态分布表可得临界值U_α。如α取0.05时, 其临界值U_α=±1.96, 当UF_k>U_α时, 时间序列存在显著的增加或减少趋势。将历年UF_k点绘成一条曲线, 可判定其是否具有增加或减少趋势。重复上述各步反序计算并将结果乘以-1, 得到新时间序列UB_k”。分别绘出UF_k和UB_k时序图, “当UF_k大于0时, 序列为增加趋势, 反之属于减少趋势。如UF_k超过临界值, 则表示增加或减少趋势达到显著水平。当UF_k与UB_k两条曲线出现交点, 且交点在临界值之间, 则交点所对应的时间可认为就是突变开始时间”^[22]。

(2) ARIMA模型。全称“自回归积分滑动平均模型(简记ARIMA)”, 是指“将非平稳时间序列转化为平稳时间序列, 然后将因变量仅对它的滞后值以及随机误差项的现值和滞后值进行回归所建立的模型”^[23]。“模型共有3个参数, 一般形式为ARIMA(p, d, q), AR为自回归, p为自回归阶次, d为时间序列成为平稳时所做的差分次数, q为移动平均阶次。如果模型建模合适, 模型残差序列一定要为白噪声序列, 模型会自动拟合预测值, 得到有较高精度的预测模型”^[24]。

2.3 数据处理

(1) 降水量变化趋势检测与突变分析的数据处理。对原始降水资料数据进行预处理, 对缺失数据和粗差进行剔除, 采用Mann-Kendall非参数检验方法, 开展三北工程区和重点建设区的年降水量、生长季降水量和非生长季降水量的多尺度变化趋势检测分析, 并针对重点县年降水量进行突变分析。

(2) 未来30年降水预测的数据处理。根据三北工程区范围内国家地面气象观测站1951—2018年739个站点的年降水量数据, 采用在Matlab中编程实现对每个站点自动筛选ARIMA模型参数和对模型残差序列进行白噪声检验, 预测得出739个站点未来30年的年降水量数据, 结果显示残差序列的自相关图均在二倍标准差内, 说明时间序列信息提取完全, 且Ljung-Box检验P值均无统计学意义, 模型可以用来预测研究区降水量。

3 结果与分析 3.1 降水量变化趋势检测与突变分析 3.1.1 三北工程区降水量变化趋势检测

根据三北工程区范围内国家地面气象观测站1951—2018年739个站点的降水日值观测数据, 计算各个站点1980—2015年逐年降水量, 并计算各个站点年降水量的非参数Mann-Kendall检验统计量Z值, 进行ArcGIS空间插值, 并按其信度将Z值分为-2.32—-1.64, -1.64—-1.28, -1.28—0, 0—-1.28, 1.28—1.64, 1.64—2.32和Z>2.32七个等级, 得到三北工程区年降水量、生长季降水量和非生长季降水量的非参数Mann-Kendall检验统计量Z值空间分布图如图 2所示。

根据年降水量、生长季降水量和非生长季降水量的非参数Mann-Kendall检验统计量Z值分布图进行空间分类统计, 得到不同信度条件下降水量变化趋势的面积百分比统计表如表 1所示。

三北工程区降水量呈增长趋势的区域面积大于呈减少趋势的区域, 年降水量和生长季降水量呈增长趋势的区域面积百分比分别为73.64%和70.10%, 主要分布在西北荒漠区；非生长季降水量呈增长趋势的区域面积百分比达92.06%, 除黄土高原南部和风沙区的少部分地区而外, 均呈增长趋势。年降水量、生长季降水量和非生长季降水量呈增长趋势且置信度为90%以上的区域面积百分比分别为45.43%、37.31%和36.79%。三北工程区生长季降水量占到年降水量的85.70%以上, 生长季降水量的增加, 有利于植被的生长和生物多样性的增加。非生长季降水量呈增长趋势的区域面积百分比达到92.06%, 但由于非生长季降水量基数小, 降水量增加总量不大。

3.1.2 重点建设区降水量变化趋势检测

将三北工程重点建设区各分区年降水量、生长季降水量和非生长季降水量的非参数Mann-Kendall检验统计量Z值进行空间分析, 对其降水量变化趋势的面积百分比进行统计, 并绘制其雷达统计图如图 3所示。

在18个重点建设区中, 生长季降水量的变化趋势与年降水量的变化趋势趋于一致, 由东向西, 从松嫩平原到河西走廊, 松辽平原、海河流域、科尔沁沙地、呼伦贝尔沙地、浑善达克沙地、晋陕峡谷、泾河渭河流域以不显著减少趋势为主；松嫩平原、毛乌素沙地、河套平原、晋西北、湟水河流域、阿拉善地区、河西走廊以不显著增长趋势和为主。西部的柴达木盆地、天山北坡谷地、准噶尔盆地南缘、塔里木盆地周边年降水量的变化趋势以通过信度为95%显著性检验的增长趋势为主；柴达木盆地和塔里木盆地周边以生长季降水量的变化趋势以通过信度为95%显著性检验的增长趋势为主, 天山北坡谷地、准噶尔盆地南缘生长季降水量的变化趋势以通过信度为90%显著性检验的增长趋势为主。非生长季降水量除了晋陕峡谷、泾河渭河流域以不显著减少趋势为主而外, 其他地区均以增长趋势为主；以增长趋势通过信度为99%的显著性检验为主的包括呼伦贝尔沙地、天山北坡谷地和准噶尔盆地南缘；以增长趋势通过信度为95%的显著性检验为主的为松嫩平原；松辽平原、海河流域、科尔沁沙地、毛乌素沙地、呼伦贝尔沙地、浑善达克沙地、河套平原、晋西北、湟水河流域、河西走廊、柴达木盆地、塔里木盆地周边和阿拉善地区以不显著减少趋势为主。

3.1.3 重点县降水量突变分析

根据三北工程5个重点县自1951年以来有数据的分县年降水量统计数据, 分别绘出5个重点县的UF_k和UB_k时序图如图 4所示。

黑龙江省望奎县年降水量的统计量UF_k在0与临界值之间, 说明其增减趋势不显著, UF_k与UB_k两条曲线出现交点, 且交点在临界值之间, 分别在1961年、1964年、2013年、2015年出现交点, 说明其对应年份年降水量突变发生。新疆库尔勒市年降水量的统计量UF_k大部分时间大于0, 且在1965年、1966年大于临界值, 年降水量呈增长趋势, 在1965年、1966年达到显著水平；1970—1972年、1979—1980年、1984—1987年UF_k小于且在临界值之内, 年降水量呈不显著减少趋势。库尔勒市的UF_k与UB_k曲线多处出现交点, 且交点在临界值之间, 年降水量突变发生频繁。内蒙古自治区磴口县年降水量的统计量UF_k大部分时间大于0, 且在临界值之内, 年降水量呈不显著增长趋势。UF_k与UB_k曲线多处出现交点, 年降水量突变发生频繁。甘肃省平凉市崆峒区的年降水量的统计量UF_k1991年以后大都小于0, 且在临界值之内, 年降水量1991年以后呈不显著减少趋势。UF_k与UB_k两条曲线出现交点, 且交点在临界值之间, 分别在1953年、1971年、1976年、1981年、2017年出现交点, 说明其对应年份年降水量突变发生。内蒙古自治区科尔沁左翼后旗的年降水量的统计量UF_k1998年以后小于0, 且在临界值之内, 年降水量1998年以后呈不显著减少趋势。UF_k与UB_k曲线多处出现交点, 说明其年降水量突变发生频繁。

3.2 未来30年年降水量预测

将各站点预测数据经ArcGIS空间插值绘制三北工程区未来30年年降水量空间分布图, 并在此基础上计算未来30年间年降水量变化数据, 绘制其分布图如图 5所示。

未来30年“三北工程区”年降水量最小值31.76 mm, 最大值1071.90 mm, 变化范围1040.14 mm, 与现状年2018年(最小值18.04 mm, 最大值984.93 mm, 变化范围966.88 mm)相比较年降水量变化范围增加73.26 mm；均值现状年301.80 mm, 预测年均值321.58 mm, 均值增加19.78 mm；由于三北工程区东西跨度大, 标准差也较大, 现状年为225.74, 预测年为215.72, 年降水量空间差异减少。未来30年间年降水量变化最多增加244.62 mm, 最多减少168.54 mm。基于未来30年年降水量及其变化空间分布, 在三北工程区、重点建设区、重点县3个空间尺度上, 可为模拟基于水资源约束的林草植被理论与现实分布格局, 开展基于水资源承载能力的林草资源优化配置提供基础数据, 为三北工程区科学恢复林草资源、发展“雨养林草植被”提供科学支撑。

4 结论与讨论

三北工程区造林面积增加及工农业和生活用水量增加, 水资源供需矛盾日益突出, 迫切需要树立“以水定林草”的发展理念, 因地制宜发展“雨养林草植被”。本文在三北工程区、重点建设区、重点县3个空间尺度上, 基于1951—2018年长时间序列降水量, 采用Mann-Kendall非参数检验、ARIMA模型预测和ArcGIS空间分析等方法, 开展年降水量、生长季降水量和非生长季降水量多尺度变化趋势与突变分析, 并对未来30年的年降水量进行预测, 主要研究结论如下：

(1) 三北工程区降水量呈增长趋势的区域面积大于呈减少趋势的区域, 年降水量、生长季降水量和非生长季降水量呈增长趋势的区域面积百分比分别为73.64%、70.10%和92.06%, 全年、生长季和非生长季降水量呈增长趋势且置信度为90%以上的区域面积百分比分别为45.43%、37.31%和36.79%。生长季降水量占到年降水量的85.70%以上, 生长季降水量的增加, 有利于植被的生长和生物多样性的增加。非生长季降水量呈增长趋势的区域面积比例达到92.06%, 但降水量增加总量不大。

(2) 18个重点建设区生长季降水量与年降水量的变化趋势趋于一致, 西部的柴达木盆地、天山北坡谷地、准噶尔盆地南缘、塔里木盆地周边年降水量与生长季降水量以通过信度为90%或95%显著性检验的增长趋势为主, 其他区域以不显著减少趋势或不显著增加趋势为主。非生长季降水量除了晋陕峡谷、泾河渭河流域以不显著减少趋势为主而外, 其他地区均以增长趋势为主。

(3) 5个重点县的UF_k与UB_k两条曲线出现交点, 说明其年降水量有突变发生, 库尔勒市、磴口县、科尔沁左翼后旗UF_k与UB_k曲线多处出现交点, 说明其年降水量突变发生频繁。

(4) 采用ARIMA模型预测得出未来30年的年降水量, 并在此基础上计算30年间的年降水量变化数据, 绘制出空间分布图, 可为开展基于水资源承载能力的林草资源优化配置提供基础数据, 为三北工程区科学恢复林草资源、发展“雨养林草植被”提供科学支撑。

综上, 降水是三北工程区水资源的根本来源, 本研究从多个方面多个尺度对该区域的降水量变化趋势的空间分布及其未来预测进行了研究, 为开展基于水资源承载能力的林草资源优化配置提供了可靠的数据基础。同时也可为三北工程区水文及气候系统研究提供参考, 有利于增强对该区域未来气候变化的认识, 也有利于加深该区域水文、自然灾害及其生态环境的影响因素的研究。