联合国政府间气候变化专门委员会(IPCC)第六次报告指出^[1], 在全球气候变暖加剧的大背景下, 干旱频率和强度呈现增加趋势, 水资源长期短缺引起的干旱可能导致土壤水分不足和地下水下降^[2], 影响植物生长和物种更替^[3—5], 从而改变生态系统的结构和功能^[6]。近年来, 愈来愈多的学者聚焦于评估干旱的生态影响, 生态干旱的概念也由此提出。Crausbay等^[7]定义21世纪的生态干旱为由于间歇性的供水不足, 导致生态系统超出其脆弱性阈值, 影响生态系统服务, 并触发人类和自然系统的反馈。我国学者粟晓玲等^[8]在此基础上提出生态干旱是由自然与人类活动引发的生态系统可获得的水量低于其需水阈值, 影响生态系统服务, 并在自然和/或人类系统中触发反馈。目前国内外学者针对生态干旱监测和评估的研究已经取得了一定的进展。如Bradford等^[9]基于SOILWAT2生态系统水平衡模型, 结合年平均温度、年平均降水以及年均降水与潜在蒸散发的比值三个指标, 来监测生态干旱的地理格局和变化趋势。Park等^[10]基于河流鱼类生态系统生态流量, 建立生态干旱监测指标。然而, 基于水文、气象要素构建的生态干旱监测指标无法定量估计生态干旱。一些学者通过构建植被指数如植被状况指数^[11]、温度植被干旱指数^[12—13]等量化生态干旱。Chang等^[14]基于降水量和潜在蒸散量之差得出生态系统的可用水量, 构建归一化生态干旱指数来量化干旱的严重程度。Zhang等^[15]基于蒸散发、气候和生态系统的直接联系, 构建标准化蒸散差干旱指数来评估干旱对生态系统的影响, 对比发现该指数比降水和温度指数更能突出气候干旱演变过程中的植被效应以及植被对干旱强度变化的响应。然而, 由于植被对干旱的响应具有滞后性, 只能反映生态系统的长期状态, 忽略了生态系统对短期干旱的响应, 例如, 当存在温度适宜的短期降水时, 土壤水分并不会立刻恢复, 植被缺水和死亡可能仍在继续^[16]。基于植被生态用水量的供需平衡关系, Jiang等^[17]提出的标准化生态缺水指数能够实时识别生态干旱事件, 确定干旱特征以制定应对措施。研究表明^[18—20], 该指数能够较好地描述中国区域生态干旱特征。

20世纪50年代以来, 嫩江流域源头区、嫩江平原等地区遭受大面积垦殖、毁林开荒, 水文调蓄和水土保持功能显著降低, 生态系统供水量与需水量间存在不平衡关系, 极易造成生态缺水, 从而引发生态干旱。目前, 已有研究从气象、水文等角度分析嫩江流域的干旱分布特征及演变规律^[21—22], 但无法揭示干旱影响生态系统的机制。近年来, 有研究聚焦于干旱影响下嫩江流域植被的动态变化^[23], 但无法实时监测生态系统用水量的供需平衡关系, 反映生态干旱的演变规律。此外, 干旱指数一般反映月、季或年尺度下的短期干湿交替变化, 难以反映长期干湿期对单次干湿事件的影响^[24]。基于此, 本研究以嫩江流域作为研究区域, 构建1980—2017年标准化生态缺水指数序列, 并结合旋转经验正交函数和游程理论分析该指数序列及其干旱特征的时空演变规律。此外, 本研究利用Pearson相关系数分析生态干旱与气象水文影响因子间的相关关系, 并利用累积和曲线识别嫩江流域生态干湿期。研究结果旨在为生态干旱和生态脆弱性的长期监测提供参考, 为嫩江流域应对气候变化和生态环境保护提供科学理论依据。

1 材料与方法 1.1 研究区域概况

嫩江发源于内蒙古境内的伊勒呼里山, 主要流经内蒙古、吉林及黑龙江边界, 全长约1370km, 是松花江的最大支流。嫩江流域位于我国松嫩平原的农牧交错带, 面积约29.7万km², 占松花江流域总面积的52%。流域地处全球气温变暖增幅最大的北方高纬度地带, 属大陆性季风气候, 夏季高温多雨, 降水量可占全年降水量70%以上, 冬季寒冷干燥。根据地形地貌特征, 嫩江流域可分为上、中、下游三段。嫩江流域上游森林密布, 湿地众多, 具有多个我国重点森林湿地自然保护区。嫩江中下游河网密布, 有广阔的平原地带, 是良好的农牧业产区。研究区区域概况如图 1所示。

1.2 数据收集与处理

本文使用的1980—2017年2m风速、2m处温度、降水、径流、相对湿度和太阳净辐射等数据来自于欧洲中期天气预报中心(ECMWF)的ERA5再分析数据集(https://cds.climate.copernicus.eu/), 并采用双线性插值法将空间分辨率统一为0.1°×0.1°。本文采用的1980年、1990年、2000年、2010年30m精度土地利用数据、归一化植被指数NDVI数据以及中国100万植被类型空间分布数据来源于中国科学院资源环境科学与数据中心(https://www.resdc.cn/)。本文采用1980—2017年全球陆地实际蒸散发数据^[25]来源于国家青藏高原科学数据中心(http://data.tpdc.ac.cn/), 原始空间分辨率为0.25°×0.25°, 双线性插值后分辨率为0.1°×0.1°。

1.3 研究方法

本研究首先计算标准化生态缺水指数SEWDI, 并计算其与水文、气象因子间的Pearson相关系数。其次根据旋转经验正交函数计算所得载荷向量的空间分布对嫩江流域进行分区, 并分析各分区下SEWDI序列的年际变化特征及变化趋势。由于SEWDI序列只能显示了短期干旱或湿润情况, 本研究利用累积和曲线识别SEWDI的长期干旱期或湿润期及其转换, 从而更清晰的辨别单次干湿事件的影响。

1.3.1 标准化生态干旱指数

标准化生态缺水指数(standardized ecological water deficit index, SEWDI)是基于陆地生态系统中植被供需水关系构建的反映生态干旱状况的指数^[17], 采用植被的实际耗水量表征生态系统供水量, 根据生态需水量与生态系统供水量之差得到生态系统缺水量。计算公式如下：

(1)

式中, EWD代表生态缺水量, mm；EWC代表生态耗水量, mm；EWR代表生态需水量, mm。

根据联合国粮食与农业组织(FAO)的定义, 在土壤水分、盐分、病虫害对植物生长或蒸发无限制的标准条件下, 生态需水量在数值上可表示为最大植被蒸散发。本研究采用FAO-56推荐的单作物系数法计算植被生态需水量, 其中作物参考蒸散发ET₀采用FAO-56推荐的Penman-Monteith法计算。参考Jiang等^[17]的研究, 根据NDVI值最大的月份及其前后两个月作为生长中期, 生长初期为生长中期前一个月, 生长后期为生长中期后一个月。嫩江流域每年NDVI值最大的月份均为每年6月, 因此本文拟定生长初期为3月, 生长发育期为4月, 生长中期为5—7月, 生长后期为8—9月。计算公式如下：

(2)

(3)

(4)

式中, k_c为作物系数, 根据FAO-56确定k_c(tab), 并通过式4进行气象修正；Δ为饱和水汽压曲线的斜率；R_n为作物表面的净辐射, MJ m^-2 d^-1；G为土壤热通量, MJ m^-2 d^-1; T为2m处平均气温, ℃；u₂为2m处风速m/s；e_s为饱合水汽压, kPa；e_a为实际水汽压, kPa；e_s-e_a为饱和水汽压差, kPa；r为湿度计常数, kPa/℃。RH_min为相对湿度, h为植被株高。

SEWDI指数干旱划分标准参考标准化降水蒸散发指数SPEI。

1.3.2 旋转正交经验函数

本研究采用旋转正交经验函数对SEWDI序列进行特征分区, 以分析不同分区的生态干旱时空分布特征。经验正交函数(Empirical Orthogonal Function, EOF)是由Pearson等^[26]首次提出的能够分离气象变量场的时空变化, 提取变量场主要特征的方法, 其基本原理为将随时间变化的气象变量场分解为相互正交的空间模态和时间系数的组合, 计算公式参考文献^[27]。然而, EOF分离的空间分布结果无法清晰表示不同地理区域的特征。针对该问题, 旋转经验正交函数(Rotated Empirical Orthogonal Function, REOF)在EOF的基础上对载荷特征向量场进行方差极大旋转变换, 可以较好地反映不同地域的变化和相对分布状况, 克服EOF的局限性, 目前广泛应用于气候和气象^[28—30]、水文^[31]、海洋^[32—33]等学科以及干旱^[34—35]的研究中。

1.3.3 累积和曲线

为放大SEWDI序列中的微小变化, 识别生态干旱的干湿变化期及其转换, 本研究计算SEWDI序列累积和并绘制累积和曲线(Cumulative Sum, CUSUM)。累积和是一种通过信号的均值和方差来检测信号波动异常值的统计方法, 能够对突变点是否发生变化及发生时刻进行估计, 在质量控制、故障检测、经济学和医学等领域应用广泛^[36—38]。累积和通过对微小变化的累积, 放大微小波动对整个序列的影响, 计算方法如下^[39]：

对于随机时间序列X = {x(k)}, k = 1, 2..., 其中

(5)

式中, μ为序列X的均值；A_k和B_k均为随机变量；τ为发生突变的时间；D为均值变化量；β为白噪声。

(6)

(7)

式中, g为非参数化CUSUM算法中的统计函数；μ₀为稳态下的均值；g_k⁺和g_k^-分别代表在均值为0附近波动的随机变量。

2 结果与分析 2.1 生态干旱相关性分析

为比较生态干旱和气象水文因子间的变化关系, 图 2为SEWDI分别与气温(T)、降水(P)、相对湿度(RH)、径流(R)、归一化植被指数(NDVI)和生态需水量(EWD)各因子的1980—2017年变化, 图 3为Pearson相关性检验结果。由于时间序列较长, 从图 2中可以看出, 1980—2017年嫩江流域总体呈现除气温升高外, 其余因子如降水、相对湿度、归一化植被指数和生态需水量均呈减少趋势。其中, 气温和降水的趋势线斜率最大, 分别为0.04和0.017。结合气温、降水与SEWDI的变化关系来看, 年均气温高、降水量少的年份, 往往同时伴有生态干旱的发生。根据Pearson相关性检验结果, 各因子与SEWDI指数间的Pearson相关性分析排序为相对湿度>NDVI>径流>降水>气温, 相关系数分别为0.59、0.47、0.46、0.46、-0.36。相对湿度与SEWDI相关性系数大于0.5, 呈显著线性关系, 其余相关系数均小于处于0.3至0.5之间, 为低度线性关系。相对湿度与生态干旱相关最为密切, 在较大程度上能够反映生态干旱。NDVI、径流和降水的相关系数接近, 且均为正相关, 说明这三个因子一定程度上能够反映生态干旱。而气温与SEWDI为负相关关系, 说明随着气温的升高, SEWDI降低, 生态干旱呈加剧趋势。

2.2 生态干旱指数的时空特征 2.2.1 REOF分区结果

为突出流域的干湿特征区的时空分布特征, 本文基于REOF方法对嫩江流域年尺度下的SEWDI进行分解, 得出嫩江流域干湿状况分布, 表 2显示了REOF分解后前4种模态的方差贡献率和累计方差贡献率。由表 2可知, 各旋转模态的方差贡献率呈现较为均匀的分布规律, 4个模态的累积方差贡献率为60.2%, 能够基本代表嫩江流域生态干旱指数的时空变化特征^{[31, 35]}。其中, 第一旋转模态的方差贡献率最大, 为19.4%, 代表嫩江流域生态干旱指数最主要的信息。图 2展示了各旋转模态的空间模态值的空间分布图。从图中可以看出, 经过极大方差旋转后分离出的各空间模态在一块区域有高载荷区, 其余区域空间模态值接近于0。根据各模态的高载荷区分布特征, 对SEWDI进行区域划分, 结果如图 2所示。分区Ⅰ代表嫩江流域源区和上游地区, 分区Ⅱ为嫩江流域东部远离干流地区, 分区Ⅲ和Ⅳ基本可代表流域中游地区和下游地区。

2.2.2 各分区下SEWDI序列时空变化特征

根据REOF分区结果, 得出1980—2017年各分区的月尺度下SEWDI指数年际变化如图 5所示。为进一步分析流域生态干旱的年际变化特征, 图 6为年尺度下SEWDI的时间变化图。21世纪以来, 嫩江流域多地区陆续被列入自然保护区, 生态保护取得一定成效。因此, 本文以2000年为分界点, 将时间序列分为前期(1980—1999年)和后期(2000—2017年), 并分别作前期和后期的变化趋势线, 对比分析2000年前后流域生态干旱的演变趋势。

图 5显示了四个分区的月尺度SEWDI值时间变化特征, 对比发现各分区SEWDI值变化趋势较为类似, 大部分时段下几乎同时出现序列的波峰或波谷。前期(1980年—1999年)生态干旱以轻旱和中旱为主, 其中分区Ⅱ和Ⅳ在1982—1983年间、分区Ⅱ在1990—1991年间、分区Ⅳ在1997—1998年间均出现序列的波谷, SEWDI值为-1.5左右, 说明对应分区在这些时段内可能发生中旱或重旱, 且以分区Ⅳ持续时间最长。从图中可以看出, 后期(2000年—2017年)嫩江流域生态干旱频率和严重性显著增加, 尤其是发生重旱和特旱的频率明显增加。在2000年—2004年间, 四个分区均发生多次不同程度的生态干旱。2003—2004年间, 四个分区均出现SEWDI>-2, 即流域发生特旱。2007—2008年间以及2012—2013年间, 流域内各分区出现不同程度的生态干旱, 其中分区Ⅰ和Ⅱ均出现特旱, 分区Ⅲ重旱, 分区Ⅳ中旱或轻旱。2012—2016年间, SEWDI值基本在0以上, 流域出现频繁的极端湿润事件。以上结果表明, 从月尺度上看, 2000年以前, 生态干旱指数SEWDI值基本上在0线上下波动, 干旱程度较轻, 而2000年后, SEWDI值波动显著增加, 流域频繁出现极端干旱或极端湿润的情况。

从图 6中可以看出, 1980—1999年四个分区SEWDI指数呈现递减趋势, 说明前期生态干旱主要呈现加剧的变化特征。2000—2017年嫩江流域90%的分区SEWDI指数趋势线为增加趋势, 说明2000年后生态干旱总体有所缓解。从趋势线斜率上看, SEWDI序列的变化幅度存在前期分区Ⅲ>分区Ⅰ>分区Ⅳ>Ⅱ, 后期分区Ⅲ>分区Ⅳ>分区Ⅰ>分区Ⅱ的分布特征。然而, 总体趋势分析忽略了序列极值对分析结果的影响, 而SEWDI值的大小只能反映单次干旱的特征, 忽略了微小变化的累加影响, 无法反映流域整体的干湿变化, 难以准确代表干旱的影响性。例如, 当流域处在较长的湿润期时, 短期干旱的发生可能不会改变流域地表径流和土壤水分。反之, 当流域长期处于干旱期时, 短期湿润对流域干旱的影响较小, 但轻旱或中旱对流域干旱的影响将被放大。因此, 利用累积和曲线识别干旱开始前流域的干湿状况, 对于确定流域发生干湿变化的时间, 提高干旱监测的准确性具有重要意义。

从累积和曲线(图 7)上看, 1980—1995年间SEWDI累积和曲线在0线上下波动, 说明SEWDI值处于正常波动状态, 流域没有明显的偏湿或偏干特征。1995年后曲线开始下降, 表明生态干旱期开始。1996—2009年间, 虽然曲线在2004—2006年间有短暂的上升期, 但总体仍以下降为主, 说明该时段下流域生态总体出现偏干趋势。2009年—2012年流域出现短暂的湿润期后转变为干旱期, 直到2012年, 流域持续较长的时间均表现为偏湿润特征。因此, 1996—2009年是嫩江流域的生态干旱期, 2010—2012年为转折期, 2012—2017年为生态湿润期。以上结果表明, 累积和曲线能够更进一步地观察到SEWDI序列的持续变化, 以及微小变化对整个序列的影响, 从而进一步识别嫩江流域生态干旱的持续性, 以及干湿期间的转换。

2.3 生态干旱特征

为进一步研究生态干旱的干旱特征及其空间分布, 本研究采用游程理论提取嫩江流域生态干旱历时、干旱强度和干旱烈度。干旱历时代表生态干旱开始至结束的持续时间；干旱烈度和干旱强度分别代表干旱发生的严重性以及干旱的强度, 其值越大代表该区域生态干旱越严重。各分区的平均干旱历时、平均干旱强度、平均干旱烈度空间分布特征如图 8所示, 后期干旱特征较前期变化比例如图 9所示。

图 8显示了1980—2017年嫩江流域生态干旱特征。从图中可以看出, 干旱烈度与干旱历时的空间分布特征基本一致, 流域南部Ⅲ区和Ⅳ区的干旱历时和干旱烈度均大于北部Ⅰ区和Ⅱ区, 且干旱强度较大的区域位于Ⅳ区南部, 说明生态干旱的严重性呈现流域北部 < 南部的特点。从图 8显示的2000—2017年嫩江流域生态干旱特征中可以看出, 2000—2017年嫩江流域基本也呈现流域生态干旱严重性北部 < 南部的特点。与前期不同的是, 后期流域南部干旱强度明显降低, 甚至低于流域北部, 而流域北部干旱强度基本不变。对比前期和后期生态干旱特征的空间分布, 发现2000年后流域干旱历时和干旱烈度均有不同幅度的增加。嫩江流域生态干旱最长时间达近15个月, 最大干旱烈度达22.61, 均位于Ⅳ区的中部, 发生于2000—2017年。

从图 9中可以看出, Ⅰ区和Ⅱ区北部以及Ⅲ、Ⅳ区中部干旱历时和干旱烈度均出现大幅增加, 增加比例最高可达1.59倍和3.91倍、2.59倍和5.22倍。嫩江流域最南部区域增长比例为负值, 即干旱历时和干旱烈度均减少。就干旱强度而言, Ⅰ区和Ⅱ区呈现增长趋势, 而Ⅲ区和Ⅳ区呈现不变或微小减少趋势。总体来看, 嫩江流域中部偏南部区域生态干旱最为严重, 较易发生中长期干旱, 且可能存在加剧趋势。而嫩江流域北部较易发生单次强度较大的生态干旱, 但每次干旱历时较短, 且同样在2000年后加剧。

3 讨论与结论 3.1 讨论

本研究基于嫩江流域1980—2017年气象数据构建标准化生态缺水指数SEWDI, 并利用Pearson相关分析法分析SEWDI与年均降水、年均气温、相对湿度、径流、NDVI的相关性。其次, 本研究利用REOF对嫩江流域生态干旱进行分区, 并利用游程理论提取干旱特征, 分析嫩江流域生态干旱特征的空间分布。从时间上看, 根据吴燕锋等^[40]对松花江流域气象水文干旱演变特征的研究, 发现嫩江流域在1996—2010年间气象干旱频发, 这与本研究中1996—2009年为生态干旱期的结果基本一致, 说明该时段内生态干旱的发生可能受到气象干旱的影响。从空间上看, 该研究中指出嫩江流域中下游平均气象干旱持续时间最长, 再结合以往对嫩江流域水文干旱的研究中指出, 流域水文干旱呈现从上游到下游逐渐增加趋势^{[22, 41]}, 与本研究中生态干旱的空间分布特征类似。根据朱光磊等^[42]的研究, 1970—2019年嫩江流域西北部参考作物蒸散量呈显著增加趋势, 东南部松嫩平原处呈减少趋势, 与本研究中生态干旱特征变化趋势的空间分布基本一致。本研究中造成嫩江流域干湿期转换以及2000年后极端干旱和极端湿润事件增加的原因尚不明确, 可能有：(1)自然气候如降水的影响；(2)农田开垦、植被破环或退耕还林、湿地保护等人类活动对生态的影响具有滞后性等原因, 可作进一步分析。

根据游程理论对干旱烈度的定义, 干旱持续时间内干旱指数之和即为干旱烈度, 因此, 干旱历时长的区域往往烈度较大, 干旱历时与干旱烈度的空间分布基本一致。而干旱强度为干旱烈度除以持续时间, 更能代表单次干旱对植被生态系统影响的严重性。本研究中, 干旱历时较长的区域位于流域西南部, 而干旱强度大的区域集中在嫩江流域北部。因此, 嫩江流域西南部地区应加强干旱监测的长期部署, 北部地区应注意应对发生短期强度较大干旱的破坏性。

3.2 结论

(1) 从相关性分析结果来看, 生态干旱与各影响因子的相关系数大小排序为：相对湿度>NDVI>径流>降水>气温, 除相对湿度外, 其余因子相关系数均小于0.5, 呈低度线性关系, 相对湿度呈显著线性相关。其中气温与SEWDI为负相关关系, 其余因子均为正相关关系。

(2) 根据月尺度下SEWDI序列的年际变化, 发现1980—1999年序列基本处于正常波动状态, 生态干旱以轻旱和中旱为主, 而2000—2017年流域极端干旱和极端湿润事件显著增加, 且多次出现重旱和特旱。结合累积和曲线, 可划分出1980—1995年为正常波动期, 1996—2009年为生态干旱期, 2010—2012年为转折期, 2012—2017年为生态湿润期。

(3) 结合REOF四个分区分析年尺度SEWDI序列的演变规律, 发现1980—1999年生态干旱总体呈现加剧趋势, 而2000年后生态干旱总体呈现缓解趋势, 这是因为在极端干旱事件增多的同时, 流域极端湿润事件也增加。

(4) 嫩江流域生态干旱不同特征空间分布存在较大差异, 且后期干旱程度明显加剧。流域西南部干旱历时和干旱烈度均大于其他区域, 较易发生长期干旱, 而干旱强度较大的区域则集中于流域北部, 该区域需提高对单次强度较大生态干旱的预防。