长江流域生态系统恢复力评价及其空间异质性研究

刘孝富，张志苗，刘柏音，王 莹，罗 镭，邱文婷

中国环境科学研究院环境信息研究所，北京 100012

长江是中国和亚洲的第一大河、世界第三大河，是中国水量最丰富的河流.长江流域以其庞大的水系孕育了丰富的生态系统类型，具有强大的水源涵养、生物繁育、释氧固碳等功能，在防洪、供水、发电、航运等方面取得了长足的发展，在我国生态安全格局与社会经济发展中具有重要的战略地位.习近平总书记提出了“要把修复长江生态环境摆在压倒性位置，共抓大保护，不搞大开发”的重大战略思想，部署谋划了一系列长江保护与发展的重大战略事项，体现了党中央对于保护好长江生态环境、永葆长江生机活力的高度重视.长江的健康与安全关乎着中华民族的永续发展，因此开展长江流域生态系统恢复力定量研究，评估长江流域健康程度和生态系统稳定性十分必要.

恢复力表示物体在压力释放后的回弹，也称之为弹性力，从20世纪70年代开始被应用于自然系统[1-2]、社会系统[3-4]、经济系统[5-6]以及复合生态系统[7-8].在生态学领域，有两个恢复力概念：一是Holling恢复力，其指系统吸收状态变量、驱动变量和参数变化并持续存在的能力，这种“持续存在”表现为不发生状态转移和质的变化[1]；二是Pimm恢复力，其指系统在遭受扰动后恢复到原有稳定态的速率或时间[9].对恢复力的定性研究较多，定量测算研究较少[10]，主要源于直接衡量恢复力的难度较大，但是可以用系统中的某些属性来替代，间接衡量恢复力，这个替代的属性必须与系统恢复力相关且能被衡量[11-12].Holling恢复力需与特定的干扰联系起来，识别干扰会对哪些参数产生影响，构建响应关系，识别临界减速的特征，并将其与临界阈值连接起来，这就决定了Holling恢复力评估更多停留在理论层面，实践验证较少.与Holling恢复力相比，Pimm恢复力较易评价，也更具有可操作性[13-15].模糊评价法是Pimm恢复力应用最广的定量评价方法，其是通过建立指标体系采用加权叠加方式求得恢复力综合指数的粗略模拟方法，被广泛应用于植被[16]、湿地[17-18]、城市[19-20]的系统恢复力评价，但其评价指标主要依靠研究者的经验判断，即使研究的案例区、指标相同，也会因权重的不同而取得不一样的结果，存在很强的主观性.为了规避主观性的影响，Lanfredi和Simoniello提出了利用植被覆盖的维持概率评价恢复力的方法[21].概率衰减法评价相对模糊评价法具有较强的客观性，无需主观筛选评价指标和确定权重，避免了经验性误差；而且数据易获取，尺度可大可小，评价过程可操作性强；同时具备空间恢复力评估潜力，可识别空间差异性，已在意大利地中海区域[21]以及中国汉江流域[22]、汶川地震灾区[15]进行了应用，证明该方法应用于生态系统恢复力评价的可行性.目前，针对长江流域如此庞大生态系统的恢复力评价及影响因素研究仍相对薄弱，因此该研究基于长时间序列的增强型植被指数(enhanced vegetation index,EVI)数据，采用概率衰减法评估长江全流域及其小流域的生态系统恢复力，识别生态系统恢复力的空间异质性，研究生态系统恢复力的影响因素，以期为长江大保护和可持续发展提供技术决策支撑.

1 研究区域与方法

1.1 研究区域

长江流域面积约为180×104km2，自西向东流经青海省、四川省、西藏自治区、云南省、重庆市、湖北省、湖南省、江西省、安徽省、江苏省、上海市等11个省(自治区、直辖市)注入东海，流域植被类型主要以针叶林、阔叶林、栽培植被、灌丛、草甸和草丛为主.根据地形地貌、河流水系、行政区划等，将长江流域划分为通天河、金沙江、雅砻江、岷沱江、嘉陵江、长江上游(宜宾至宜昌)、长江中游(宜昌至湖口)、长江下游(湖口以下)、乌江、汉江、洞庭湖水系、鄱阳湖水系、太湖水系等13个二级流域，细化为6 767个小流域(见图1).

图1 长江二级流域空间分布Fig.1 Spatial distribution map of thesecond-level basins in the Yangtze River Basin

1.2 研究方法

1.2.1 数据源与数据处理

该研究采用MOD13Q1产品EVI数据，MOD13Q1是搭载在Terra卫星上MODIS传感器生成的植被指数专题数据，投影方式为Sinusoidal，空间分辨率250 m，每隔16 d提供一期数据，即每月两期数据.该研究数据时间为2001−2020年，月时相为1−12月，每年数据23期，每期8景，共3 680景数据.数据来自美国地质调查局网站(http://ladsweb.nascom.nasa.gov/data/search.html).

由于长江流域属亚热带湿润季风性气候，降雨量大，植被茂密，蒸腾作用强烈，常年云雾较多，单期MODIS-EVI数据可能因云量较大对计算结果造成较大误差，因此该研究采用最大值合成方法来计算年度EVI最大值〔见式(1)〕.

式中，EVIy为y年份的EVI，EVIi为第i期的EVI.

1.2.2 概率衰减法

1.2.2.1 基本原理

生态系统可以看作由若干“栅格”组成，每个栅格的EVI随时间的推移呈增大或减小，斜率表现为正数或负数.时间越长，栅格的EVI变化斜率始终维持为正或维持为负的概率越低.降低的趋势可以采用指数衰减函数进行拟合，正负趋势衰减时间的差值可衡量生态系统恢复力.当正向衰减速率小于负向衰减时，生态系统恢复力较好，反之则较差，差值越大表示生态系统恢复力越高.

1.2.2.2 技术方法

生态系统恢复力评价需包括时段划分、栅格EVI变化斜率计算、斜率栅格表面计算、维持性栅格表面计算、指数衰减拟合、生态系统恢复力计算等步骤.

a)时段划分.首先，将长时间序列的EVI划分为起始时段和拟合时段，拟合时段需有充分的样本数据来拟合衰减函数.2010年三峡水库蓄水到175 m，是长江流域生态系统变化过程中具有标志性的事件，因此该研究将2001−2010年作为起始时段，2011−2020年作为拟合时段，计算公式：

式中，T(t0,t)为 研究时段，T(t0,t′)为 起始时段，T(t′,t)为拟合时段，t0为初始年，t'为参考年.

b)栅格斜率计算.将长江流域划分为若干个1 km×1 km的栅格，计算各栅格在初始年(t0)到参考年(t′)以及初始年(t0)到评估年(ti)的EVI变化斜率.以最小二乘法拟合各栅格参考时段内EVI的变化线，斜率的计算公式：

式中，Slope为斜率，EVI 为评估时段的EVI平均值，t为评估时段的平均时间.

c)斜率栅格表面计算.根据栅格的EVI变化斜率来计算各年份斜率栅格表面，即当栅格的EVI变化斜率(Slope)为正时，将斜率栅格表面赋值为1；当栅格的EVI变化斜率为负时，斜率栅格表面赋值为−1.

式中，S(ti)为评估年的斜率栅格表面.

d)维持性栅格表面计算.根据各年份的斜率栅格表面计算维持性栅格表面.参考年的维持性栅格表面与斜率栅格表面相同，将某年的斜率栅格表面与上一年的维持性栅格表面进行比较，如果栅格值相等，则该栅格的维持性栅格表面赋值与上一年相同，反之则赋为0.

式中，P(ti)为 评估年维持性栅格表面，P(ti−1)为上一年维持性栅格表面，S(ti)为评估年的斜率栅格表面，S(t′)为参考年的斜率栅格表面.

e)指数衰减拟合.以参考年维持性栅格数量(Nt′)为基准，将其余年份维持性栅格数量(Nti)进行归一化，即计算各年维持性栅格数量占参考年维持性栅格数量的比例〔q(ti)〕.采用指数衰减函数拟合维持性栅格比例随时间下降的趋势，计算得到衰减时间(τ).计算公式：

式中：τ为趋势衰减时间，a；A为常数.

f)生态系统恢复力计算.某趋势衰减时间越长，表示维持该变化趋势的概率越高.生态系统恢复力(RS)是通过正趋势衰减时间(τp)与负趋势衰减时间(τn)的差值来衡量，当正趋势的衰减时间大于负趋势时，表明该生态系统的恢复力好，反之生态系统恢复力差，计算公式：

2 结果与讨论

2.1 长江全流域生态系统恢复力

2.1.1 EVI年际变化

长江流域EVI的年际变化趋势和空间分布特征(见图2)表明，2001−2020年长江流域EVI总体呈增加趋势，增幅为7.76%，说明长江流域植被覆盖总体逐渐变好.2001年的EVI最小，仅为0.490；EVI峰值出现在2016年，达0.536，之后稍有降低，2020年为0.532.随着时间的推移，EVI的标准差也在逐渐增加，说明长江流域植被覆盖空间分布差异在逐渐增大.长江流域植被覆盖呈西部低、中部高的空间分布规律，EVI较高的地区主要为四川省、贵州省、湖南省、湖北省等地区，长江入海口区域EVI在逐渐变小.

图2 2001−2020年长江流域EVI变化及空间分布情况Fig.2 Variation of average EVIand spatial distribution in the Yangtze River Basin from 2001 to 2020

2.1.2 维持性栅格数量变化

统计2010−2020年EVI变化维持正、负趋势的栅格数量，计算各年份维持性栅格数量占2010年维持相同趋势栅格数量的比例(见表1).由表1可见，随着时间的推移，EVI变化维持负趋势栅格数量减幅较大，降速较快，2020年EVI变化维持负趋势栅格数量仅为2010年的33.82%；而EVI变化维持正趋势栅格数量占2010年的76.58%，说明长江流域植被整体上维持着正向演替.

表1 2001—2020年长江流域EVI变化维持正、负趋势的栅格数量统计Table 1 Statistics of the number of EVI maintained grid with positive and negative trends in the Yangtze River Basin from 2001 to 2020

采用空间分析方法将2010−2020年维持性栅格进行叠加，得出每个栅格在11个年份中维持同一趋势的次数(见图3).维持正趋势的年份次数越多(绿色越深)表明维持正趋势的时间越长，维持正趋势时间长的栅格主要分布在陕西省的安康市、商洛市，湖北省的十堰市，以及重庆市的城口县、巫溪县、巫山县、奉节县等区域；维持负趋势的年份次数越多(红色越深)表明维持负趋势的时间越长，维持负趋势时间长的栅格主要分布在上海市，江苏省的苏州市、无锡市、常州市、镇江市、南京市，四川省的成都市，以及浙江省的湖州市、嘉兴市等区域.

图3 2001−2020年长江流域栅格维持趋势的时间Fig.3 Time of grid maintaining trend in the Yangtze River Basin from 2001 to 2020

2.1.3 生态系统恢复力

以2010−2020年维持正、负趋势栅格数量占参考时刻的比例为纵坐标，时间为横坐标，拟合趋势衰减曲线得到正趋势的衰减时间为38.02 a，负趋势的衰减时间为8.54 a，差值为29.48 a，表明在前期的干扰状态下，2010年后长江流域生态系统整体上在逐渐恢复.

2.2 分流域生态系统恢复力

统计长江流域6 767个小流域各年维持不同趋势的栅格数量，计算各年每个小流域维持性栅格数量占2010年维持相同趋势栅格数量的比例，拟合正负趋势衰减曲线，得到各小流域生态系统恢复力.由于小流域数量多、面积小，每个小流域包含的栅格数量较少，拟合时可能出现一些结果异常情况，采取如下处理方法：①有些小流域在2010−2020年维持正、负趋势栅格的比例保持不变，衰减时间将呈现无限大，对该小流域做出标记；②某些小流域仅存在维持正趋势或维持负趋势的栅格，即以维持该趋势栅格数量衰减时间为生态系统恢复力；③对于无统计值的小流域或正负趋势衰减时间均为无穷大的小流域进行标记.

综合考虑长江各小流域的生态系统恢复力大小范围、排序情况、数量占比等，根据自然断点法原理，结合相关研究结果[22]，将长江小流域生态系统恢复力划分为5个等级(见表2).由表2可见，长江流域生态系统恢复力大于0的小流域面积为1 415 014.73 km2，占流域总面积的81.65%，说明长江流域大部分区域生态系统为正向演替.生态系统恢复力等级为一般和较好的小流域较多，分别占流域总面积的47.58%、24.60%；生态系统恢复力极好的小流域面积为164 060.30 km2，占流域总面积的9.47%；生态系统恢复力极差的小流域面积为18 425.86 km2，仅占流域总面积的1.06%.

表2 长江小流域生态系统恢复力(RS)等级划分Table 2 Classification of ecological resilience in the small watershed of the Yangtze River Basin

根据长江小流域生态系统恢复力空间分布(见图4)得出，生态系统恢复力高的小流域集中在长江流域中部地区，主要分布在嘉陵江、汉江上游、长江上游(宜宾至宜昌段)、乌江流域、洞庭湖水系、鄱阳湖水系等二级流域，尤其是陕西省的安康市、商洛市，甘肃省的陇南市，四川省的广元市，湖北省的神农架林区、恩施土家族苗族自治州，湖南省的张家界市、湘西土家族苗族自治州、永州市、怀化市，广西壮族自治区的桂林市，江西省的吉安市、赣州市、抚州市、景德镇市，以及安徽省的黄山市等区域生态系统恢复力较好，正负趋势的衰减时间差最高可达13 628 a；生态系统恢复力差的小流域分布在太湖水系、长江下游(湖口以下)、金沙江、岷沱江等流域，主要是四川省的成都市、眉山市、自贡市，青海省的玉树藏族自治州，江苏省的南通市、泰州市、常州市、镇江市、苏州市、南京市，浙江省的嘉兴市、湖州市，上海市，安徽省芜湖市、马鞍山市、合肥市，湖北省的潜江市、孝感市、武汉市、鄂州市以及湖南省长沙市等地区，正负趋势的衰减时间差最低可达−838 a.长江流域生态系统恢复力总体上呈现中部高、东部和西部低的特点，长江流域源头区和长江入海口区域生态系统恢复力较差.

图4 长江小流域生态系统恢复力空间分布情况Fig.4 Spatial distribution of ecological resilience in the small watershed of the Yangtze River Basin

2.3 讨论

2.3.1 长江流域生态系统恢复力分析

从时间尺度上看，2001−2020年长江流域整体上EVI呈增加趋势，EVI维持下降趋势的栅格减少速率较快，长江流域整体植被覆盖正向演替趋势明显，与相关研究结果[23-26]一致.从空间尺度上看，植被覆盖空间分布差异性在逐渐增大，植被覆盖维持上升趋势的栅格主要分布在汉江上游、长江上游(宜宾至宜昌)、乌江中下游、洞庭湖水系上游来水区域，集中在陕西省的安康市、商洛市，湖北省的十堰市、重庆市的城口县、巫溪县、巫山县、奉节县等地区；植被覆盖保持下降趋势的栅格主要分布在太湖水系以及岷沱江流域中下游，如上海市，江苏省苏州市、无锡市、常州市，四川省成都市等地区，与相关研究结果[15,22-23,26-30]一致.

长江流域生态系统恢复力为正的小流域面积占流域总面积的81.65%，说明长江流域大部分区域生态系统具有一定的稳定性，为正向变化趋势.生态系统恢复力较好的小流域分布与植被覆盖持续变好区域基本重叠，集中在汉江上游、长江上游(宜宾至宜昌段)等区域，表明良好的生态环境基础为生态系统稳定性提供了保障.太湖水系的南通市、嘉兴市，以及成都市、眉山市、芜湖市等地区生态系统恢复力差，这些地区主要位于成渝经济圈、长江三角洲经济圈，经济发展水平较高，城市化进程快，说明人类的高强度干扰对生态系统恢复力有较大影响，若干扰到一定程度生态系统将发生逆转[31-34].笔者研究揭示了长江流域生态系统恢复力空间分布差异，可为流域生态系统保护、环境治理与经济发展协同管理提供参考.因此，长江流域下游省市在国土开发过程中要重视生态环境保护和修复，开展综合性、针对性、差异性管控[35-36].

2.3.2 生态系统恢复力空间异质性分析

针对长江流域生态系统恢复力的空间差异性，深入分析其影响因素.参考李观凤等[7,31,37]对于生态系统脆弱性和恢复力驱动因子的分析结果，笔者研究以长江流域6 767个小流域为对象，分析小流域生态系统恢复力与降雨量、地形、植被覆盖、人类活动等因素之间的相关性，判断各因素对生态系统恢复力的影响情况.

2.3.2.1 降雨量的异质性

分析长江小流域生态系统恢复力与2013−2020年多年平均降雨量(自然条件)、降雨量变化趋势之间的关系(见图5).由图5可见：多年平均降雨量在800~1 200 mm范围内的小流域，生态系统恢复力较高；多年平均降雨量<1 200 mm时，生态系统恢复力随着降雨量的增大总体呈增加趋势；多年平均降雨量>1 200 mm后，生态系统恢复力有所减小.生态系统恢复力与降雨量变化斜率之间无明显相关特征规律.

图5 长江小流域生态系统恢复力与多年平均降雨量、年降雨变化斜率的关系Fig.5 Relationship between ecological resilience and perennial mean rainfall and annual rainfall variation slope in the small watershed of the Yangtze River Basin

2.3.2.2 地形异质性

分析小流域生态系统恢复力与平均海拔、平均坡度的关系(见图6)发现，随着海拔升高，生态系统恢复力呈先增大后减小的趋势，海拔在1 000~1 500 m范围内的小流域生态系统恢复力最好. 海拔低于1 000 m时，随着海拔的升高，生态系统恢复力趋向于升高；海拔高于1 000 m后，随着海拔的升高，生态系统恢复力趋于降低.坡度≤30°时，随着坡度增大，生态系统恢复力呈增加趋势；坡度高于30°时，生态系统恢复力有所下降；坡度在25°~30°范围内的小流域生态系统恢复力最大.

2.3.2.3 植被覆盖异质性

以林草覆盖率表征小流域植被覆盖特征，分析长江小流域生态系统恢复力与植被覆盖之间的关系(见图7)，发现生态系统恢复力随着林草覆盖率的增大而升高.

图7 长江小流域生态系统恢复力与林草覆盖率的关系Fig.7 Relationship between ecological resilience and perennial mean vegetation index in the small watershed of the Yangtze River Basin

2.3.2.4 人类活动异质性

以建设用地占比表征小流域人类活动强度，分析长江小流域生态系统恢复力与人类活动之间的相关性(见图8).由图8可见，建设用地占比在10%以上时，小流域生态系统恢复力随着建设用地占比的增大而减小，说明生态系统恢复力会受人类活动的影响.

图8 长江小流域生态系统恢复力与建设用地占比关系Fig.8 Relationship between ecological resilience and proportion of construction land in the small watershed of the Yangtze River Basin

通过以上影响因素分析结果发现，生态系统恢复力与降雨量、地形(海拔、坡度)、植被覆盖、人类活动均有一定的相关性，生态系统恢复力随坡度、林草覆盖率的升高而增大，随人类活动强度的增大而减小，随海拔的增加呈先增大后减小的趋势.笔者研究表明，生态系统恢复力的影响因素存在明显的时空差异，海拔在1 000~1 500 m、坡度在25°~30°、年降雨量在1 000~1 200 m的小流域生态系统恢复力较高，证明生态系统恢复力是受多种驱动因素综合影响，不同区域内影响因素的作用强度和波动方向各不相同.空间异质性分析也从侧面验证了小流域生态系统恢复力的空间分布特征，在长江入海口上海市、江苏省等地区，海拔较低平、降水充沛、植被较少，土地多样化程度高，城镇化建设快速的地区生态系统恢复力相对较低.

3 结论

a)2001−2020年，长江流域整体EVI呈上升趋势，生态环境呈正向变化，植被覆盖的空间差异性在逐渐增大.长江流域上游植被覆盖一直较少，且改善程度不明显；长江入海口部分区域植被覆盖率总体呈减小趋势.

b)长江流域内81.65%的区域生态系统表现为稳定向好的状态.长江流域生态系统恢复力呈现空间分布不均匀的特征，总体上呈中部高、东部和西部低，长江入海口部分地区生态系统恢复力较差的特征.生态系统恢复力受地形(海拔、坡度)、植被覆盖、人类活动、降雨量的影响，生态系统恢复力与坡度、林草覆盖率均呈正相关，与人类活动呈负相关，与海拔呈抛物型相关.长江流域生态系统恢复力受多种驱动因素的综合影响，其影响在上、中、下游既有相似性又存在显著差异，其空间不一致性与生态系统的自然条件有关.

c)长江流域生态系统整体上呈正向演替，特别是汉江流域，这得益于长江大保护以及丹江口水源地的保护，无论从国家政策还是相关人员的研究都可以得到印证.生态系统恢复力最低值出现在长江下游及入海口区，这些区域人类干扰强度大，生态系统稳定性差，在开发过程中要特别重视生态保护与修复，避免产生生态倒退的不可逆现象.上游地区如金沙江、大渡河、岷江一带生态恢复力较低，是否与生态系统结构、水电开发、地震及地质灾害相关，还需进一步研究.