近百年来全球气温持续上升, 根据IPCC发布的《气候变化2021：公众摘要》报告显示：与19世纪末(工业革命前)相比, 2011年至2020年地球表面的平均气温升高了1.1℃, 并且比过去12.5万年的任何时候都高^[1]。历年的IPCC发布的报告均认为：近百年来, 人类过度的开荒及大规模活动, 导致所排放的温室气体有增无减, 这是造成全球变暖最主要的原因^[2]。人类活动不仅导致了全球气候的变暖, 还诱发了如极端天气、植被退化、水土流失等生态安全问题。生态环境是人民生存的基础条件, 人类可持续发展的保障。而植被作为生态系统主要的组成部分, 对区域生态系统环境变化有着重要指示, 对全球碳循环和气候调节具有非常重要的作用。植被覆盖度(FVC)是指植被在地面的垂直投影面积占总面积的百分比, 它量化了植被的茂密程度, 反映了植被的生长态势, 是描述生态系统变化的重要基础数据^[3]。怒江-萨尔温江是东南亚最长的自由流动河流和最重要的跨界河流之一, 在中国-中南半岛经济走廊建设和中缅经济走廊建设和生态保护中发挥着不可替代的作用。在过去的50年里, 该流域冬季和春季气温呈上升趋势, 气温上升的幅度随着海拔高度的升高而增加^[4]。然而, 目前对该流域关注不够, 尚没看到其植被覆盖研究报道。因此, 研究分析其植被覆盖时空演变, 并探讨其驱动力, 可为其生态环境保护和社会经济发展提供重要指导作用^[5]。

估算植被覆盖度(FVC)是植被覆盖相关研究的重要基础, 有多种估算方法, 如地表实测法、经验模型法、植被指数法、决策树分类法等, 部分方法存在对实测数据依赖性强、计算繁琐、需要大量样本数据等缺点, 特别对于偏远地区来说很难实现^[6], 而遥感技术的快速发展以及遥感大数据时代的到来, 使得在大时空尺度上精确监测偏远地区的植被行为成为可能^[7], 其中MODIS NDVI产品被大量学者用于植被覆盖度的计算及其动态变化的定量分析, 而像元二分模型法是常用的FVC估算方法, 其操作较为简单, 可直接利用NDVI数据来进行计算^[8—10]。目前针对植被覆盖演变规律的研究主要是运用一些常规的线性趋势分析方法, 如线性回归分析方法^[11—12]、相关性分析法^[13—14]、Theil-Sen斜率与Mann-Kendall检验结合的趋势分析法^[15], 很少有学者关注其非线性变化特征。然而, 植被覆盖的时空分布与变化受多方面的因素影响, 使得其演变规律也具有复杂性, 并不仅仅呈现为简单的线性变化, 而BFAST法可用于分析不同类型的时间序列, 对长时间序列过程中趋势的突变检测具有很好的应用效果, 已有学者将其用于植被的非线性特征^[15—17], 以及气温变化^[18]、气溶胶变化^[19]等领域研究中, 并显示了其较好的适用性^[20—21]。对于南北跨度极大的怒江-萨尔温江流域, 其不同纬度区域的植被变化特征及原因也会存在差异, 因此有必要运用BFAST法来更加准确地探究其植被变化的非线性特征及其未来发展趋势。

对于植被未来可持续性的研究, Hurst指数法是广大学者常用的方法^[22—23]。耦合BFAST法与Hurst指数可以实现植被未来趋势的更加精确表达。驱动植被变化的因素, 有人为因素和自然因素, 以往学者大多采用偏相关分析、回归分析等方法研究植被覆盖度变化的影响因素^[24—25], 这些方法假设了植被与各环境因子存在线性关系。然而植被在环境中的生长过程影响因素复杂, 可能不存在严格的线性关系^[26]。而由王劲峰等学者提出的地理探测器模型(Geographical Detector, 简称Geodetector), 不仅可用于探测空间分异性并识别其主要的影响因子, 还能实现不同因子间的交互作用探测^[27]。此方法无线性假设, 可以避免多重共线性问题, 在计算驱动因子的贡献率时结果较为准确^[28]。其改进的模型OPGD模型(Optimal Parameters-Based Geographical Detector, OPGD)加入了参数优化模块, 可根据自变量的特征进行多种方式和多种断点数目的离散化处理, 选取出最优的离散化方式和最优的断点数目, 提高各类别间的显著差异性^[29]。因此, 研究采用该方法进行植被非线性特征的驱动因子探测。

基于以上分析可知, 尽管已经开展了许多关于植被时空演变与驱动因素的研究, 但目前还没有学者采用BFAST法研究其历史非线性趋势的同时结合Hurst指数表达其未来更加精确的可持续性, 也就是考虑了植被非线性特征的未来趋势预测。另外, 怒江-萨尔温江是东南亚一条重要的国际河流, 对南亚、东南亚的气候调节具有重要作用, 但由于其南北跨度大, 涉及中国、缅甸和泰国, 对其进行相关研究尤其数据收集较为困难, 目前还没有学者对该流域进行植被覆盖的相关研究。因此, 研究以该流域为研究对象, 基于MODIS NDVI数据, 在遥感与GIS技术的支持下, 采用像元二分模型估算其近22年的植被覆盖度, 并借助BFAST01模型及Hurst指数等, 对区域内植被覆盖的时空分布、非线性变化趋势和未来可持续性进行分析, 并结合坡度、海拔、降雨、气温、人口、GDP、土地类型等自然和人类活动因子, 引入OPGD对影响研究区植被覆盖度的主要因素进行定量探究, 以期为研究区的生态环境保护和可持续经济发展提供科学依据。

1 研究区概况及数据来源 1.1 研究区概况

怒江-萨尔温江发源于青藏高原上的唐古拉山脉, 经中国云南(保山、临沧)流入缅甸, 最后注入印度洋的安达曼海, 经过中国、缅甸和泰国(图 1)。该流域形状极为狭长, 河流全长约3200 km, 流域总面积为32.5万km², 它在中国的部分被称为怒江, 其余的部分称为萨尔温江，是东南亚最长的自由流动河流, 也是最重要的跨境河流之一。杨帆等的研究中, 将嘉玉桥以上划分为上游, 嘉玉桥到中缅边界为中游, 中国边界以南为下游^[30]。由于该流域南北跨度和海拔高差都很大, 横跨多个气候带和生态区, 其气候多样, 地形复杂, 不同区域的地形、气候、土地利用类型都存在很大的空间差异性, 本研究根据流域海拔特征并结合行政区划界限, 将云南和西藏交界线以北区域划分为上游, 属于青藏高原, 为典型的高原山地气候区, 气候变化较为敏感；云南西藏交界线到中缅边界划分为中游, 主要为高山峡谷区, 即狭窄又陡峭, 气候垂直差异明显, 流经三江并流保护区；中国边界以南为下游, 海拔相对北部较低, 地势较为平坦, 属于热带季风气候, 以此分区作为后期相关研究和分析的基础。

1.2 数据源及预处理

研究数据主要包括MODIS NDVI数据, 地形数据, 土地利用数据, 气象数据以及人类活动数据, 研究区流域边界和行政区划数据等。详述如下：

1.2.1 NDVI数据

研究使用2000—2021年MODIS NDVI数据, 通过GEE(Google Earth Engine, GEE)遥感大数据云平台进行调用、预处理及下载, 数据原网站为USGS (https://lpdaac.usgs.gov/products/), 其时间分辨率为16d, 空间分辨率为250m。在GEE平台中对其进行镶嵌、投影、转换和裁剪等预处理, 为了消除云、大气效应、扫描角度和太阳天顶角等对影像质量的影响, 研究采用了最大值合成法MVC(Maximum Value Composite)来进行NDVI数据的合成^[15], 该方法被广泛用于植被动态研究的预处理中^[31]。借助ENVI IDL平台, 选用像元二分模型(DPM)来计算研究区月度及年度FVC, 用于后续植被覆盖的相关分析。

1.2.2 地形数据

研究使用的地形数据为2000年采集的SRTM Digital Elevation Data Version 4, 分辨率为90m, 数据生产于NASA (https://srtm.csi.cgiar.org/), 通过GEE平台对其进行拼接、裁剪和下载, 并通过ArcGIS 10.2软件提取如海拔、坡度、坡向等地形因素作为后续植被覆盖驱动力分析的指标。

1.2.3 LULC数据

研究使用的土地利用数据为ESA WorldCover 10m v100, 来自esa (https://esa-worldcover.org/en), 也是通过GEE平台对其进行投影转换和裁剪等预处理和下载。该产品将土地利用类型分为了11类, 比较详细的表示了不同地表覆盖类型, 我们将其作为植被覆盖的人类活动驱动因子之一。

1.2.4 气象数据

研究使用的气象数据主要包括2000—2021年的气温、降水、蒸散发以及气候类型。其中气温数据由ERA5提供的第五代ECMWF全球气候大气再分析数据, 通过GEE平台进行预处理与下载。我们选用的是mean_2m_air_temperature, 单位为K；降水量数据选用的是Climate Hazards Group InfraRed Precipitation With Station Data 2.0 (https://chc.ucsb.edu/data/chirps), 单位为mm, 分辨率为0.05°。蒸散发数据为GLEAM ET产品, 空间分辨率为0.25°, 下载于gleam (https://www.gleam.eu/#downloads)。蒸散发可以反映植被的干旱情况, 对植被生长状况有一定影响, 因此研究选取该指标作为驱动因子之一。气候类型数据为全球柯本气候类型空间分布数据集, 空间分辨率为0.1°, 下载于国家地球系统科学数据中心(http://www.geodata.cn/), 将其作为植被覆盖驱动力分析因子之一。

1.2.5 人类活动等其他数据

人口格网数据来源于world pop (https://hub.worldpop.org/geodata/listing?id=64)。单位为people/km², 空间分辨率为1km。全球GDP数据来源于GitHub (https://github.com/Nowosad/global_population_and_gdp), 分辨率为0.5° × 0.5°。中国区域人口、GDP数据来源于中国科学院资源环境科学数据中心(http://www.resdc.cn)。其他数据如研究区行政区划, 流域边界等分别来源于中国科学院资源环境科学数据中心(https://www.resdc.cn/data.aspx?DATAID=259) 和HydroSHEDS网站(https://www.hydrosheds.org/)。数据来源如表 1, 驱动因子指标如图 2。

2 研究方法 2.1 像元二分模型

像元二分模型是一种简单实用的遥感估算模型, 对绿度信息具有很好的敏感性, 且对不同气候区、植被类型等的植被信息具有很好的识别效果, 被广泛用于计算植被覆盖度。其计算植被覆盖度的公式如下：

(1)

式中, NDVI_soil和NDVI_veg表示全部为裸土和全部为植被的像元NDVI值。理论上, NDVI_soil为0, NDVI_veg为1。然而, 实际中植被覆盖的影响因素复杂, NDVI_soil和NDVI_veg的值并不等于0和1, 而是受图像本身质量的影响。参考相关研究^{[23, 32]}, 根据MODIS NDVI时间序列图像, 对NDVI图像进行频率直方图统计, 去除异常值以及归一化处理, 并取累计频率在2%和98%时的NDVI值分别作为NDVI_soil和NDVI_veg参数, 然后按照公式(1)计算研究区FVC。

2.2 Hurst指数

Hurst指数是定量描述时间序列长程依赖性的有效方法^[33]。近年来在植被覆盖变化的未来可持续性研究中得到广泛应用^[34—35]。研究利用matlab计算出H值即Hurst指数, 来判断植被覆盖时间序列是否存在持续性。H值介于0—1间, 参考相关研究成果^{[23, 36]}, 其分类有三种：0＜H＜0.5(反持续性)、H=0.5(无明显变化)、0.5＜H＜1(持续性), H值越大表明持续性越强, 反之越弱。

2.3 BFAST模型

BFAST是一种迭代算法, 它利用分段线性趋势和季节性模型将时间序列分解为趋势、季节性和残差组分, 并能检测趋势和季节性成分的突变^{[7, 37]}。BFAST假设非线性可以用分段线性模型来近似, 因此, 当一个序列可以更好地用多个线性段表示, 而不是单一的单调趋势时, 就可以使用结构变化测试确定一个断点并确定日期^[26]。相比其他突变检测方法, 在时间序列内, BFAST受季节差异和噪声的影响较小, 可以更快地检测到变化^[19]。Verbesselt等对该方法进行了比较完整的描述和介绍^[38—39]。在数学上, BFAST的趋势分量和季节分量通过以下分解模型得到：

(2)

式中, Y_t为t时观测到的值, T_t为趋势组分, S_t为季节性组分, e_t为残差组分。其中趋势组分可用下式表示：

(3)

季节性成分可用下式表示：

(4)

式中, a_i和b_i为趋势项系数, γ_i为振幅, δ_i为分段数量, f为频率, 其中振幅、分段数量为未知项, 而频率为已知项。

季节组分和趋势组分中突变点的识别需要确定突变点的数量和突变点在时间序列中的位置。此处利用最小二乘移动求和从季节组分和趋势组分中检验判断是否存在突变点^[40]。首先, 采用普通最小二乘方法(Ordinary Least Squares, OLS)在一个不断增加的移动窗口内进行检验判断是否存在断点, 如果根据统计检验, 窗口内的普通最小二乘残差的移动和很大, 则被描述为一个断点。接着, 通过最小化所有窗口上断点之间的最小二乘移动总和(OLS-MOSUM), 以此来确定断点的位置^{[19, 38]}。

根据BFAST原理, 本研究采用了BFAST的改进模型BFAST01。BFAST01模型相比BFAST模型的优势在于, 它同时考虑了季节模型和趋势模型, 只检测时间序列中最具影响的趋势变化, 而不是多个较小的变化, 也就是说BFAST01只检测0个或者1个断点, 其为评估长期季节性和年际变化的主要转折点提供了一种有效的、全面的变化检测方法^{[26, 41]}。本研究利用R环境, 从GitHub加载最新的bfast包(https://github.com/), 完成所有BFAST01的统计分析^[42], 提取出植被覆盖的趋势显著性、变化强度、趋势断点发生的时间和次数、突变类型、突变年份等。参考相关研究^{[15, 19, 26]}, 将BFAST01检测到的植被覆盖变化趋势类别分为8类, 如表 2所示。

2.4 地理探测器

地理探测器Geodetector, 是一种统计工具, 可以探测地理现象的空间分异性并揭示其背后驱动力, 既适用于点数据, 也适用于面数据, 对于面数据一般要先进行重分类和离散化处理^[27]。其包括4个探测器, 我们选用因子探测器和交互探测器来探测植被覆盖变化的影响因素。

因子探测器, 用于探测因变量Y的空间分异性, 及某个因素X对因变量Y空间分布的解释力, 其大小可用q值来度量(0≤q≤1), 本研究中, 若某影响因子的q值越接近于1, 说明该因子对FVC的解释力越强, 反之越弱^{[27, 43—44]}。

交互探测器可以用于识别不同因子之间交互对植被覆盖度变化的影响是增强还是减弱, 或是独立^[27]。其交互方式见表 3。

3 结果与分析 3.1 怒江-萨尔温江流域植被覆盖时空特征分析 3.1.1 空间分布特征

通过像元二分模型计算得到FVC, 并利用GIS空间分析可视化工具得到研究区2000—2021年平均FVC空间分布图(图 3), 对其进行分析可知, 怒江—萨尔温江流域低植被覆盖度(0—0.2)、较低植被覆盖度(0.2—0.4)、中植被覆盖度(0.4—0.6)、较高植被覆盖度(0.6—0.8)、高植被覆盖度(0.8—1)占比依次为：5.41%, 5.66%, 10.41%, 21.46%, 57.06%, FVC分布具有明显的空间异质性, 从图 3可以看出, 下游和中游地区植被覆盖情况明显优于上游地区, 总体表现为南高北低。整个流域平均FVC值为0.73, 高植被覆盖度和较高植被覆盖度地区面积占比达到78.52%, 说明该流域植被覆盖情况较好。结合研究区土地利用类型分析发现, 低植被覆盖区域主要分布在海拔较高的青藏高原, 主要土地利用类型为苔藓、裸地、灌木等；高植被覆盖区域主要在中游和下游地区, 土地利用类型以林地、草地和耕地为主。

3.1.2 时间变化特征

通过统计研究区年、月、季节的FVC均值, 并制作FVC变化折线图(图 4)。从图 4的年际变化图可以看出, 怒江-萨尔温江流域2000—2021年植被覆盖整体呈波动上升趋势, FVC最大值出现在2021年为0.750, 最小值出现在2015年为0.722, 在2007年和2015年出现了明显的低值拐点, 但2016年以来, 植被生长状况出现明显上升趋势。从研究区连续22年植被覆盖季节变化图可以看出怒江-萨尔温江流域植被变化具有明显的季节规律性, 这与研究区气候有关。而FVC月际变化图表明, 3月份FVC值最低, 3—9月份FVC值随时间逐渐上升, 在9月份的时候FVC值达到最高, 随后呈下降趋势；整体上夏秋季植被覆盖状况较春冬季植被覆盖状况好, 不同季节FVC随着时间变化呈现较小的波动, 春季和夏季波动较秋季和冬季大, 但整体较为稳定。

3.2 怒江-萨尔温江流域植被覆盖非线性趋势特征

运用BFAST01算法对怒江-萨尔温江流域2000—2021年FVC进行突变检测, 得到FVC的非线性趋势特征(图 5)。从图中可知, 研究区分为了8种非线性植被覆盖趋势类型, 其中单调递增占34.63%, 是占比最多的, 在流域上、中、下游各区域均有分布；单调递减类型占比12.22%, 主要分布在流域下游区域；“单调递增+”占比2.3%, 主要分布在中缅边界；“单调递减-”占比最少为0.63%；“中断-+”占比22.91%, 排名第二, 主要分布在流域最北端和云南以及缅甸北部；“中断+-”占比8.19%, 主要分布在流域最南角, 云南和北部也有些许分布；“反转+-”占比7.73%, “反转-+”占比11.39%；通过分析发现, 总体呈增加趋势类型的占比(71.24%)大于总体呈减少趋势类型的占比(28.76%), 表明研究区植被发展趋势向好。从显著性检测结果来看, 变化呈现显著的占比为44.29%, 不显著的占比为55.71%；突变时间检测表明, 发生突变的时间主要集中在2003—2018年, 2003年和2018年发生突变的数量较多, 其他年份较为均匀, “中断-+”发生的突变最多, 发生突变占比较少的类型是：“单调递增+”和“单调递减-”。需要提醒的是, 单调递增和单调递减没有发生断点, 因此图 5(突变时间)中空白区域表示没有发生突变的区域。

为了更加详细的描述植被覆盖突变情况, 研究还对FVC变化的最大突变强度、最大突变强度发生的时间、首次季节性断点发生的时间以及断点发生次数等信息进行了检测(图 6)。通过分析发现整个流域突变发生最大强度的地区主要在云南以及青藏高原北部, 最大突变强度范围介于-1.2078—1.0287之间, 最大突变强度时间主要在2002—2019年, 2006年发生最大强度突变的比例是最大的, 2013年、2015年和2018年最大突变强度数量也较多, 上游地区最大突变强度断点发生时间变化趋势与此相似, 说明上游地区在全流域FVC的变化趋势中具有重要作用；而中游和下游区域最大突变强度断点发生时间变化趋势较上游地区波动平稳, 中游地区突变数量较多的年份为2006年、2010年和2015年, 下游地区在2011年最大强度突变数量最多, 占下游地区突变总面积的9.90%。各区域最大强度突变面积从2018年开始大面积减少。此外, 根据研究区FVC时间变化可知, 其植被变化具有很强的季节性, 从图 6可以看出, 首次发生季节性断点的时间主要集中在2010—2014年。

此外, 怒江-萨尔温江流域最大突变强度断点发生突变的次数分别有1次、2次和3次(图 7), 表示某点在研究时间段内植被发生非线性突变次数的空间分布, 1次的占比最多为77.02%, 3次的最少为4.18%, 2次的区域主要分布在中缅边界和最北部。上游地区发生1次断点的比重最大, 占上游地区发生断点次数总面积的85.72%, 而发生3次断点的占比仅占本区断点总面积的0.91%；下游和中游地区各断点次数占本区域断点总次数面积比重相差不大, 也是以1次为主, 占比略小于上游地区, 而发生3次断点的次数占比均在7%左右, 略大于上游地区。

3.3 怒江-萨尔温江流域植被未来可持续性

为了表示植被覆盖的未来可持续性, 我们计算了怒江-萨尔温江流域FVC的Hurst指数, 通过统计分析发现, Hurst指数小于0.5的区域占总面积的2.76%, 这些区域具有反持续性, 表示过去增加的趋势将来总体上减少, 反之亦然；Hurst指数大于0.5的区域占比为94.89%, 表明研究区植被覆盖的正向持续性较强；Hurst指数为0.5的区域占2.35%, 表明未来植被覆盖发展趋势不明确。

为了进一步揭示植被覆盖的非线性变化趋势及其持续性, 将BFAST变化趋势结果与Hurst指数结果进行叠加, 得到非线性变化趋势与持续性的耦合信息(图 8和表 4)。耦合结果总共为17种情形, 其中“持续性&单调递增”耦合类型含义表示未来植被与过去发展趋势一样, 将持续增加；“反持续性&单调递增”表示未来植被与过去发展趋势相反, 将从增加变为减少, 其它情形解释相似, 比较详细的展示了植被覆盖的未来发展趋势。未来将呈持续改善情景的面积占比达到68.99%, 持续改善的趋势类型分别为：“持续性&反转-+”、“持续性&中断-+”、“持续性&单调增加+”、“持续性&单调递增”、“反持续性&单调递减”、“反持续性&单调递减-”、“反持续性&中断+-”和“反持续性&反转+-”, 这些区域占据了流域大部分地区；未来发展将呈持续退化的情景所占比例为29.09%, 这些趋势类型分别为：“反持续性&反转-+”、“反持续性&中断-+”、“反持续性&单调递增+”、“反持续性&单调递增”、“持续性&单调递减”、“持续性&单调递减-”、“持续性&中断+-”和“持续性&反转+-”, 这些区域主要分布在流域中游和下游；剩下1.92%的区域未来变化趋势无法确定, 主要分布在流域上游区域。持续退化和未来变化趋势无法确定的区域的植被变化状况需要研究人员继续关注。

3.4 怒江-萨尔温江流域植被覆盖驱动因子探测

研究选取了海拔、坡度、坡向、蒸散发、气温、降水、气候区划等自然因素和GDP、人口、土地利用等人为因素作为植被覆盖的影响因子, 分别为X₁—X₁₀, 植被覆盖度FVC为Y, 作为地理探测器的数据输入, 空间分辨率统一重采样为250m。由于这些因子及FVC都是连续变量, 而地理探测器的输入数据要求为类别数据, 则根据王劲峰等提出的数据离散化方法^[27], 利用OPGD在R语言环境中进行参数优化, 为了使结果更加符合实际, 研究同时对离散化进行了人为调试, 最终得到比较理想的结果, 具体离散化方法和类别见表 5。并利用ArcGIS创建渔网, 通过多次调试设置采样间隔为3km, 运用Spatial Analyst-提取分析-采样工具, 输入栅格为所有重分类的X和Y, 采样点为研究区渔网点, 将自变量和因变量的值提取到点, 作为地理探测器输入数据。

利用因子探测和交互探测器, 得到全流域以及各子流域植被覆盖驱动力探测结果(表 6和图 9)。除中游地区的坡向因子p值为0.0791外, 其余因子p值均为0, 通过显著性检验。根据表 6可知, 全流域中海拔的q值是最大的为0.5761, 接着影响较大的是气温、降水、蒸散发等气象因子, 其次是土地利用类型、GDP、人口等人为因子, 人为因子中以土地利用影响最大。各子流域因子探测结果显示同一个影响因子在不同区域对植被覆盖的影响具有差异性, 纵向上分析, 如海拔和土地利用因子对上游和中游区域的影响大于下游区域；而人口、GDP、蒸散发等因子对下游区域的影响则大于上游和中游区域。横向上分析, 上游区域影响最大的是土地利用因子, 其次是海拔, 而人口、GDP等人类活动因子对于较为偏远的上游地区影响较小。对于中游区域, 对植被影响最大的是海拔, 其次是土地利用, 气温、降水、蒸散发等也是影响中游区域植被覆盖变化的重要因素, 自然因素对中游区域植被的影响大于人为因素。与中游区域相反, 下游区域的人口和GDP是对植被影响最大的因素, 其次是蒸散发, 土地利用和降水影响也较大, 而影响最低的是海拔。

交互探测结果显示(图 9), 任何两个因子交互影响都会增强其影响力, 对于全流域来说, 发现其它每个因子和海拔交互都达到最大影响值, 说明海拔对植被覆盖有着很强的影响力, 这与整个研究区海拔差异显著具有关联性。上游区域和中游区域与全区相似, 各因子与海拔的交互作用对植被覆盖的影响增强作用是最显著的。对于下游区域来说, 各因子交互作用具有差异性, 其中气温、降水、蒸散发与坡度的交互作用对植被的影响大于它们与其它因子的交互作用；气候类型、土地利用、人口与蒸散发的交互作用对植被的影响大于它们与其它因子的交互作用；说明坡度和蒸散发对下游区域的植被覆盖影响具有很大的可利用空间；而GDP与人口的交互作用也大于GDP与其它任何因子的交互作用, 成为所有交互作用中影响最大的, 这与因子探测的结果具有一致性。

4 讨论 4.1 怒江-萨尔温江流域植被覆盖时空总体特征

怒江-萨尔温江流域区位特殊, 上、中、下游地区地形海拔和纬度等差异明显, 植被覆盖空间分布特征也形成了鲜明的对比。从全区域来看, 高植被覆盖区域占比大于50%, 而其中上游地区高植被覆盖区域占比仅占7%左右, 平均FVC值为0.53；中游和下游区域高植被覆盖区域占比达到85%和90%以上, 平均FVC值分别为0.86和0.88, 说明中游和下游地区拉高了整个区域的植被覆盖水平, 这与研究区地形地貌有关^[45—46], 这将在驱动力分析中进行探讨。

从FVC时间演变特征来看, 怒江-萨尔温江流域植被覆盖在2015年以前增长趋势不明显, 但2016年开始发展趋势向好。这与李杰、王春雅等的研究结论相似^{[23, 47]}。但在2007年、2015年和2019年出现了明显的低值拐点, 结合图 4与图 10, 发现上游地区FVC年际变化也在2007年和2015年发生了低值拐点, 变化趋势与全区FVC非常相似；而下游和中游区域FVC年际变化增长趋势明显, 其中下游地区在2004、2012、2018等年份出现低值拐点, 查询资料发现, 2004年中缅边境发生了一场大森林火灾, 2006年的飓风“马拉”和2007年的台风“利奇马”, 2012年和2013年缅甸发生了严重干旱, 2014年、2015年的厄尔尼诺现象、极端干旱事件等, 此外, 据联合国人类事务协调办公室称, 缅甸在2015年经历了几十年来最严重的洪水, 同时2015年“4·25”西藏也发生了地震灾害, 以及2018、2019年的极端气候自然灾害如印度西南部的特大洪灾等^[48], 这些事件均严重威胁到当地植被的生长状态, 可能是导致FVC出现低值拐点的原因。对于下游地区来说, 有研究表明森林火灾、农业转型、公路建设增加等是导致其植被减少的主要原因^[49]。而研究区整体植被变化趋势向好与生态工程的实施和围栏禁牧政策以及其他人类活动等有关^[47]。

4.2 怒江-萨尔温江流域植被覆盖非线性特征

从研究区FVC年、月、季节的变化分析以及相关研究表明, 植被变化往往受多种因素的影响使其很难呈现严格的线性变化^[15]。因此, 我们在进行研究时应该考虑其趋势的非线性特征及其影响因素, 特别是大尺度长时间序列的演变分析。研究运用BFAST01方法来分析FVC的演变趋势以及非线性突变, BFAST方法可以很好地适用于突变检测^[50]。在3.2小节中, 我们详细分析了怒江-萨尔温江流域FVC变化的突变类型、突变时间、突变显著性与突变次数等特征, 将研究区植被变化详细分为了8种突变类型, 每种突变类型的含义见表 2。为了进一步分析各突变类型空间分布特征及验证分类的准确性, 我们结合Google Earth time-series images来进行突变前后年份影像的比对(图 11)。

从上、中、下游区域分别选择了某局部区域进行制图比对, 上游区域选择的代表像元突变类型为“单递减-”, 突变时间为2017年, 从Google Earth time-series images上我们可以明显看到该区域范围突变时间前后土地利用从草地和稀疏植被逐渐转变为建设用地, FVC图像也可以看出其颜色逐渐变深, 也就是FVC值在减小, 这就说明了该区域植被覆盖在减少, 这和突变类型“单调递减-”的特征是相符的。中游区域选择的代表像元突变类型为“中断+-”, 该突变类型有1个明显突变点, 且断点处值突然增大, 即显著正中断, 然后显著减小, 从Google Earth time-series images和FVC图像比对来看, 也符合其突变特征。下游区域选择的两个突变类型代表分别为“中断-+”和“反转-+”, “中断-+”特征与“中断+- 相反, 从图像特征上可以看出该局部区域植被覆盖在断点2012年处突然减小, 即负中断, 然后增加；而“反转-+”表示趋势从显著减小转换为显著增加, 即在突变点处FVC显著减小, 也符合Google Earth time-series images和FVC图像的比对特征。因此, 通过图像比对验证, 发现用BFAST模型进行长时间序列植被非线性趋势检测是适用的, 结果也具有一定可靠性, 研究可为相关生态环境保护措施的制定提供科学参考依据。

4.3 怒江-萨尔温江流域植被覆盖驱动力

为了更有针对性地分析植被覆盖的影响因素, 考虑研究区地形地貌的差异性, 研究选择了适用于具有非线性特征分析的基于最优参数的地理探测器(OPGD), 对研究区植被覆盖影响因素进行分区探讨, 根据探测结果进行更有针对性的原因分析。探测结果表明, 对整个区域来说, 海拔的影响是最大的, 原因在于研究区海拔南北差异极大, 对整个区域的地形地貌起着宏观控制作用, 进而影响局部地区的水热条件, 对植被空间分异产生重要影响。此外, 影响较大的是气温、降水、蒸散发等气象因子, 这与许多研究以及现实情况也是相符的^[51—53], 植被受气候影响较大, 特别是怒江-萨尔温江流域跨越了多个不同气候区, 气候差异大, 因此气候因子对其影响较显著；人为因子对全区植被的影响程度比自然因子较低, 其中以土地利用影响最大, 原因在于研究区地形偏远, 人口经济发展较缓慢, 特别是高海拔地区和高山峡谷区, 人类活动很难影响到, 而土地利用作为能较好反映人类活动的因素, 能很好的表现植被覆盖的变化, 同时也有很大的可利用空间来增强或者减弱植被覆盖。

对各子流域来说, 植被覆盖的影响因素具有差异性。上游区域土地利用对植被覆盖的影响大于海拔, 原因在于区域内部大部分属于高原, 高原内部地势起伏不大, 因此海拔以及坡度、坡向等对植被的影响小于土地利用, 人类活动通过改变土地利用结构, 特别是耕地面积、耕作类型等从而影响植被覆盖；在剩下的因子中降水是影响较大的, 原因在于上游地区季风难以到达, 降水又是植被生长的重要因素, 因此, 降水对植被影响较大。由于其较为偏远, 人为因素影响小于自然因素。中游区域, 由于其大部分属于高山峡谷区, 海拔差异显著, 海拔成为了植被覆盖影响最大的因素, 而气温、降水、蒸散发等受海拔差异的影响, 对植被覆盖影响也较大。下游区域, 由于以平原为主, 地形较为平坦, 人类活动频繁, 因此人口和GDP是对植被影响最大的因素, 由于纬度低, 气温高, 蒸散发也成为了植被变化的重要影响因素, 与上游和中游地区不一样的是海拔对植被影响微弱。

综上可知, 由于怒江-萨尔温江流域地理特征的复杂性, 在不同环境条件下, 各环境因子对植被覆盖的影响也存在显著的差异性。因此, 建议综合考虑不同环境条件下FVC的空间分布与各影响因子的空间关联性, 明确不同区域的环境制约因子, 因地制宜制定生态修复措施, 规划流域的经济建设。对于上游地区来说, 植被覆盖较低, 未来植被呈退化趋势的占比较多为9.03%, 应防止植被覆盖率进一步下降加强生态恢复, 如优化土地利用结构、发展现代生态农业、更改耕作类型、进行轮作和休耕管理、退耕还林、建设自然保护区和实施生态移民、以及争取国家政策的扶持等, 注重社会经济和生态环境的共同发展；对于中游地区, 自然环境因子对植被覆盖影响更大, 其经过三江源等重要生态源地区, 应该确保生态环境不被破坏；对于下游地区, 植被覆盖度较高, 但影响植被覆盖的主要是社会经济因素, 且未来植被呈退化趋势的占比是最多的为17.39%, 应切实保护好现有资源, 维持生态系统的稳定性, 并发展高效、多元化、可持续的生态经济。植被动态是一个复杂的过程, 还应加强对极端事件的预防和控制, 以防止气候变化造成植被的进一步退化。总之, 应充分发挥人类活动对于生态环境的积极作用, 找准生态保护与经济发展的平衡点, 促进生态与经济的协调可持续发展。

4.4 局限和展望

研究参考杨帆等^[30]的成果重新将怒江-萨尔温江流域划分为了上游、中游和下游区域进行植被趋势分析与驱动力探讨, 我们主要考虑的是海拔和行政区划, 区域划分较为主观, 这是未来研究需要改善的地方, 比如利用具体的气候区、土壤区、地貌区、植被区等进行分析。此外, 研究在进行驱动力指标选取的时候, 由于跨境各个国家的数据统计标准不一, 精度也有差异, 使得部分数据获取较为困难, 特别是人为因素的数量和质量, 这可能会影响结果的精度, 因此, 未来应该考虑在分区探讨的同时对不同区域建立不同的指标体系, 以提高数据的精度和结果的准确性。

5 结论

研究基于怒江-萨尔温江流域2000—2021年MODIS NDVI数据, 计算出其植被覆盖度FVC, 探讨了研究区近22年植被覆盖的时空趋势变化、非线性特征、未来可持续性以及驱动力。主要结论为：

(1) 近22年, 怒江-萨尔温江流域植被覆盖总体发展趋势向好, 2007年和2015年出现较大突变, 多年平均FVC值为0.73。受水热条件、地形地貌等影响, 植被分布具有明显的空间异质性, 下游和中游区域植被覆盖明显优于上游区域, 总体表现为南高北低, 植被覆盖状况较好。(2)BFAST趋势表明, 近22年怒江-萨尔温江流域植被覆盖改善和退化的区域面积占比分别为71.24%、28.76%, 改善的区域远大于退化的区域, 说明研究区植被得到较好的保护。Hurst指数显示, 未来植被相比以前将得到改善和将有所退化的区域占比分别为94.89%、2.76%。BFAST与Hurst二者的耦合表达了植被覆盖的未来非线性趋势, 也是改善的区域大于退化的区域, 占比分别为68.99%和29.09%。应特别关注未来发展趋势退化的区域。(3)地理探测器结果表明各区域植被覆盖影响因素具有差异性, 针对全区域海拔的影响是最大的, 进一步说明了地形与研究区植被具有密切关系。而针对各子流域, 土地利用对上游地区植被覆盖的影响比中游和下游区域显著, 中游地区海拔的影响力最大, 下游地区以人口、GDP等人为因素影响为主。因此, 要合理根据研究区不同区域的环境制约因子等制定合理的植被保护措施。