嘉陵江流域作为长江流域中最大的子流域，是长江流域生态环境监测和修复重点区域。近年来，在人口急剧增长、城市快速扩张，以及全球极端气候的多重压力下，该流域的生态环境正面临着严峻的考验，比如小流域经济发展与环境保护的双重目标压力、水生态环境退化以及矿业开采导致的次生环境风险等^[1]，这些问题严重威胁着长江经济带建设战略目标的实现。因此，急需开展嘉陵江流域长时间尺度的生态环境监测与评估，进而为流域生态文明建设提供数据支撑。

目前，生态环境质量评价技术方法已从对生态环境的单因子调查发展到多因子综合指数评价^[2]，比如生态环境状况指数(Ecological Index，EI)就是一种常用的评价指标。但由于数据获取难、权重和归一化系数设定不合理等问题，在实际应用中，EI对大尺度、长时序、高频次和高精度的生态监测而言具有较大的难度和不确定性。相对于EI，遥感生态指数(Remote Sensing Ecological Index，RSEI)是一种基于遥感技术的生态评价方法，该方法可利用自然环境因子对区域生态环境状况实现快速精准的监测与评估^[3]，具有客观、准确等优势，且结果可比性强，易于空间可视化，因此现已被广泛应用于生态环境质量评价领域^[4-5]。

多年RSEI的构建需要长时间序列的遥感影像作为基本输入，但是受云层覆盖和遥感影像时间分辨率等因素的限制，较大区域尺度、长时间序列的遥感影像数据往往较难获取。因此，以往研究多基于非时序或者跳跃式时间间隔的遥感影像来构建RSEI，例如，王丽春等^[6]仅以2000、2006、2016年3期的Landsat遥感影像为数据源，借助RSEI分析了玛纳斯湖湿地生态环境状况；欧阳玲等^[7]采用1990、2000、2010、2015、2020年5期的遥感指数，对科尔沁沙地生态环境状况进行了动态评价。但上述方式仅仅反映流域生态环境质量变化的间断趋势，较难对连续变化趋势进行系统分析^[8]。另外，大尺度的RSEI建模面临着庞大的数据量、繁杂的数据预处理和指标计算等瓶颈问题，因此对海量遥感数据的处理能力提出了更高的要求^[4]，而行星级地理空间分析云平台谷歌地球引擎(Google Earth Engine，GEE)^[9]为用户提供了丰富的开源数据和强大的计算服务，为长时间序列、大尺度的RSEI快速批量时序计算提供了可能性^[10-11]。已有学者采用此方法对流域进行生态监测和评价，并取得较好的效果^[12-13]，表明GEE平台在长时间序列、大尺度的遥感生态环境监测研究中具有广阔的应用前景^{[2, 8, 14-15]}。

基于此，本研究借助GEE平台开展以下研究内容：(1)采用像元级最小云量影像合成方法，构建2001—2022年的无云遥感影像，并对嘉陵江流域目标年份逐年的长时间序列RSEI进行提取、突变点检测和趋势分析，得出流域生态环境质量时空演变特征；(2)综合对比RSEI与EI两种方法的评价结果，阐明两种方法在本流域对不同生态环境质量等级的表征能力；(3)使用地理探测器进行归因分析，探明影响流域生态环境质量的关键因子，揭示各驱动因子间的协同作用。本研究成果将为流域生态环境质量状况评价有效方法的选择，以及流域生态环境质量保护政策的制定、生态恢复工程成效的评估提供科学依据和技术支持。

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

嘉陵江流域地理范围为29° 17′30″—34° 28′11″N，102°35′36″—109°01′08″E，横跨青藏高原气候、温带季风气候和亚热带季风气候3个气候分区，地形包括高原区、盆地区和丘陵区^[16]。流域跨越陕、甘、川3省，于重庆汇入长江主河道，流域包括三大水系(嘉陵江、渠江、涪江)和8个支流域，总面积约15.95万km²，是长江水系中集水面积最大的流域(图 1)。流域大部分地区属亚热带季风气候，其上游段为高海拔山区，冬季严寒且多霜；而中下游段为盆地，冬季温暖多雾^[17]。流域年均温约10 ℃，年降水量在1 000 mm以上。受季风影响，流域降雨年内分配不均，再加上人类活动的频繁干扰，流域内水土流失严重，滑坡、泥石流等地质灾害频发，使得嘉陵江成为长江水系支流中含沙量最高的河流^[17-18]。

1.2 数据来源

考虑数据的可获得性和可操作性，选取2020年嘉陵江流域的EI与同期RSEI进行对比。其中，土地利用数据源为2020版30 m全球地表覆盖数据GlobeLand30 (https://www.webmap.cn/commres.do?method=globeIndex)。水系数据来自全国地理信息资源目录服务系统的1∶25万全国基础地理数据库(https://www.webmap.cn/commres.do?method=result25W)。土壤胁迫指数提取自资源环境科学数据平台的中国土壤侵蚀类型和程度空间分布数据(https://www.resdc.cn/)。污染物负荷指数采用国家地球系统科学数据中心的中国1 km分辨率逐年近地表NO₂、PM_2.5、SO₂数据集(http://www.geodata.cn/)进行计算。以上数据均用于EI的计算。

地理探测器自变量选取自然因子(气候因子、地形因子)和社会因子(人口密度、土地利用类型、夜间灯光指数)，这些因子具有可获得性、代表性，常用于空间分异性探测和定量归因分析。其中，气候因子包括年降水量和年均温，年降水量和年均温数据集均来源于国家地球系统科学数据中心(http://www.geodata.cn/)。地形因子海拔、坡度由数字高程模型(Digital Elevation Mode，DEM)提取获得。人口密度数据来源于WorldPop (www.worldpop.org)，空间分辨率为1 km×1 km，数据单位为人/km²。土地利用类型数据源与EI土地利用数据源相同。夜间灯光遥感数据来自全球变化科学研究数据出版系统(https://www.geodoi.ac.cn/)的中国长时间序列夜间灯光数据集。另外，考虑到数据的可获得性，自变量和因变量均选取2020年数据作为地理探测器的数据源进行归因探测分析。

1.3 基于遥感的生态环境质量指标模型

RSEI的数据源提取和指数构建均在GEE平台上完成，其具体方法和步骤如下。

(1) 最小云量影像合成。嘉陵江流域常为多云天气，使用传统的遥感数据获取和处理方法难以得到质量较好的长时间序列数据，而使用GEE平台可以优化重构年度无云影像^[2]，特别适用于常年多云覆盖区域的影像提取。本研究借助GEE平台，对输入的符合时空范围的MOD09A1_500 m地表反射率产品、MOD11A2_1000 m地表温度产品、MOD13A1_500 m植被指数产品进行计算，去除有云像元，进行无云筛选和逐年拼接裁剪等操作，再以无云像元重构目标年份夏季(6—9月)最小云量合成影像，得到嘉陵江流域逐年无云拼接影像。

(2) 分量指标计算与归一化。基于上述合成的多年遥感影像，提取植被指数(NDVI)、土壤指数(NDSI)、湿度分量(WET)、地表温度(LST)来分别代表绿度、干度、湿度和热度这4个重要的生态要素，并作为RSEI的构建因子^[3]，即：

$\mathrm{RSEI}=f(\mathrm{NDVI}, \mathrm{NDSI}, \mathrm{WET}, \mathrm{LST}), $

其中，各分量计算公式详见徐涵秋^[3]的研究。RSEI值为0—1，数值越高代表等级越高，生态状况越好。

1.4 长时间序列趋势分析

长时间序列RSEI的变异性采用变异系数(Coefficient of Variation，CV)来表征，该系数可以很好地反映空间数据在时间序列上变化的差异程度，评价数据时间序列的稳定性，其计算公式如下：

$\mathrm{CV}=\frac{\sigma}{\mu}, $

式中：σ为标准差，μ为算术平均值。CV越大，数据波动越大。

结合Mann-Kendall检验(简称M-K检验)^[18-19]和Theil-Sen Median趋势分析进行突变检测与趋势显著性分析。M-K检验是世界气象组织(WMO)推荐的一种非参数检验方法，已广泛应用于长时间序列的趋势检验与突变分析^[20-22]。Theil-Sen Median趋势分析又被称为Sen斜率估计，是一种稳健的非参数统计的趋势计算方法。本研究采用R语言“terra”^[23]、“trend”^[24]包运算M-K检验和Theil-Sen Median趋势分析。统计量Z的取值范围为[-∞, +∞]。在给定显著性水平α下，当|Z|>μ_1-α/2时，表示研究序列在水平上存在显著变化。一般置信水平设置为0.05^[22]，则|Z|>1.96表示研究序列变化显著，|Z|≤1.96表示研究序列变化不显著。将M-K检验分级结果和Theil-Sen Median趋势分析结果进行叠加，得到像元尺度上的长时序RSEI趋势变化，将其划分如下：“-2”表示严重退化、“-1”表示轻微退化、“0”表示稳定不变、“1”表示轻微改善、“2”表示明显改善。

Hurst指数(H值)用来定量描述时间序列长期依赖性^[25]，可以反映RSEI未来变化趋势相对于过去的依赖性，即时间序列是完全随机或是存在持续性。若0.5<H值<1.0，表示时间序列是一个持续性序列，即未来变化与过去的变化趋势一致，此时H值越接近1.0，持续性越强；若H值=0.5，则说明时间序列为随机序列，不存在长期相关性；若0<H值<0.5，则表示时间序列具有反持续性，即未来的变化趋势与过去的变化趋势相反，H值越接近0，反持续性越强^[22]。

1.5 驱动因子分析

地理探测器是探测空间分异性、揭示其背后驱动力的一种新统计学方法^[26]。相较于传统方法，地理探测器对多自变量共线性免疫，因而不必对变量进行相关性分析；但是其自变量必须为类型变量，因此需对自变量进行离散化处理^[26]。本研究使用R语言“GD”包^[27]，先采用多种离散化方法处理连续变量，再从中筛选出各变量的最优空间离散化方法，最后利用因子探测器和交互作用探测器进行数据分析。另外，为有效表征各生态环境质量因子的最佳适宜范围或类型^[28]，进一步运行风险探测器进行探测。

(1) 因子探测器

因子探测器可探测因变量Y(即RSEI)的空间分异性，以及不同影响因子Xi对因变量的空间分异性及影响程度，并检验解释程度的可信程度。解释程度用q表示，q值表示Xi解释了q×100%的Y，其值介于0与1之间。q值越大，表示空间分层异质性越强，Xi对Y的解释程度越强，反之越弱。计算公式详见文献[26]。

(2) 交互作用探测器

交互作用探测器用于判断Xi的两两交互作用对Y的空间分异性解释程度的增强或减弱作用^[29]。通过比较单个因子X1和X2分别对Y的解释程度q(X1)和q(X2)，以及X1、X2两个影响因子交互作用时对Y的解释程度q(X1∩X2)，确定影响因子交互作用的类型，或者这些因子对Y的影响是否相互独立。根据表 1中判别方法进行探测。

(3) 风险探测器

风险探测器是地理探测器模型中的一个核心子模块，其本质是一种通过统计量t来检验不同子区域间因变量均值是否存在显著差异的空间统计工具。该工具主要用来精准识别研究区域内风险的高值区与低值区，并通过统计验证判断这些区域间的差异是否具有显著性，从而为决策者提供科学依据。

1.6 EI的计算

根据《生态环境状况评价技术规范》(HJ 192—2015)(简称《规范》)^[30]，采用生物丰度指标、植被覆盖指数、水网密度指数、土地退化指数和污染负荷指数，按照不同权重计算得到EI值，具体方法和参数见《规范》。本研究采用相对易于获取的2020年数据，将该年份嘉陵江流域EI值进行提取和空间可视化后，与RSEI分级结果进行对比。按照《规范》中生态环境状况分级对结果进行划分，将EI分级结果划分为5级，即EI值≥75表示生态环境质量为“优”，55≤EI值<75表示生态环境质量为“良”，35≤EI值<55表示生态环境质量为“一般”，20≤EI值<35表示生态环境质量为“较差”，EI值＜20表示生态环境质量为“差”。RSEI值为0—1，按照等分法将分级结果也划分为5个等级，便于与EI分级结果进行对比。

2 结果与分析 2.1 生态环境质量时空变化特征

使用RSEI方法获取的2001—2022年嘉陵江流域生态环境质量有明显的空间差异和时序变化特征(图 2)。例如，西北部的白水江、北部嘉陵江干流流域大部分年份的RSEI值较高，生态环境质量好；涪江流域东南部、东河/西河流域、渠江流域则相对较差。2001、2010年嘉陵江流域整体生态环境质量好，2002—2005年则明显下降。RSEI年均值在2001年达到峰值，在2003、2005、2019年出现低值，整体呈缓慢增加趋势；整体上，2001—2022年研究区域的RSEI年均值为0.22—0.69，多数年份生态环境质量处于一般水平(图 3)。

通过CV、Hurst指数和M-K检验分析发现，流域多年的生态环境质量波动明显，除了白河流域的西北部生态环境质量出现较高波动，其他绝大部分区域出现了高波动[图 4(a)]。严重退化区域面积极少，轻微退化区域占比为8.11%(表 2)，主要分布在涪江流域、东河/西河流域的西南部，以及白河流域西南部、渠江流域部分区域；轻微改善区域面积占比最大(74.42%，表 2)，其次是明显改善区域(13.50%，表 2)，二者主要分布在嘉陵江流域海拔较高的北部和东部。整体上，嘉陵江流域的生态环境质量呈现向好的趋势[图 4(b)]，且多年RSEI变化在大部分区域表现为持续性序列(占79.55%，表 2)，即未来变化与过去的变化趋势一致，但在涪江流域、白水江流域、渠江流域部分地区以及嘉陵江流域南部的部分地区，未来变化可能与过去趋势不一致[图 4(c)]。

2.2 驱动因素分析 2.2.1 主要驱动因素

因子探测器分析结果显示7个因子均通过显著性检验(P<0.05)。其中，年降水量、年均温、海拔位列前三，其q值均大于0.3(图 5)，表明这些因子最大程度地解释了RSEI值的空间分异，是研究区生态环境质量的主要驱动因素。人口密度、土地利用类型的q值也较高，分别为0.231 4和0.141 9。坡度和夜间灯光指数的解释度较低。上述结果说明嘉陵江流域生态环境质量的主要驱动因素为自然因素，但人为因素也是不可忽视的重要驱动因素，尤其是人口密度较高的区域。

对比分析2.1节嘉陵江流域生态环境质量时空变化特征可知，流域生态环境质量长时间序列结果出现了明显波动的时间点，如2001、2003—2005、2010年。已有研究报道，年降水量在此时间段也有相似的波动，2006、2011年前后嘉陵江流域夏季发生了极大的干旱，降雨量较其他年份显著下降，对区域植被生长产生严重影响^[31-32]，且近30年来嘉陵江流域气温和降水均显著上升^[17]，这与上述驱动因素分析结果基本相符。

2.2.2 交互作用结果

表 3显示了所有因子在交互作用下均具有协同增强作用。其中，年降水量在其他变量因子的协同作用下，对RSEI空间分异特征的解释力最大，影响最显著，且与年均温、海拔的协同作用最强，其与人口密度、夜间灯光指数、土地利用类型也存在较强的协同作用。该结果表明自然因子和社会因子相互影响，嘉陵江流域生态环境质量变化是多种因子共同作用的结果，虽然自然因子解释力更大，但是社会因子也起到不可忽视的作用。上述结果有助于进一步理解嘉陵江流域生态环境质量变化的驱动机制。

2.2.3 风险探测分析结果

表征嘉陵江流域生态环境质量的7种因子最优区间如图 6所示。各个因子的生态适宜范围如下：年降水量为688—796 mm，年均温为12.3—14.7 ℃，海拔在1 340—1 770 m，人口密度在7.07—21.8人/km²，夜间灯光最少，坡度为25.7°—31.8°且为林地的区域。

2.3 RSEI和EI分级结果比对

将RSEI和EI对生态环境质量的表征能力进行对比，结果发现二者整体上的空间变化基本一致。其中，RSEI分级结果中“优”“良”比例高，“较差”和“差”级别区分明显，分别占有一定比例。而EI值为20.78—100，因此分级结果中无“差”级别，且“优”和“较差”级别均很少，大部分区域是“良”和“一般”等级(图 7)。RSEI分级结果中的“优”“良”等级在EI分级结果中多表现为“良”，“较差”和“差”级别在EI分级结果中则多表现为“一般”等级。相比之下，EI在空间上区分度不高，结果多集中在中间区间值。这在一定程度上表明，RSEI对于“优”和“差”两个极端更敏感，能更好地区分空间上的变化情况。

此外，对比EI和RSEI的两种指标，发现EI方法虽然在时空连续监测上有局限性，但是其指标中包括了RSEI方法中没有的污染负荷指数、土地退化指数等。因此，EI分级结果可与RSEI分级结果对比矫正，这样更有助于客观全面地了解区域特征，进而弥补RSEI方法不能完全覆盖和解释生态状况优劣方面的不足。例如，灌木、外来物种入侵草场可能伴随着绿度的提高，而RSEI方法无法有效区分生态退化带来的变绿与生态改善的变绿^[15]。在全球气候急剧变化和人类活动频繁干扰的背景下，流域宏观生态环境质量急需长期且连续的监测。这不仅需要以GEE平台等为代表的新兴技术的支持，同时也需要不同尺度的生态环境监测和评估方法的不断创新，以便及时发现问题并做出科学准确的应对。

3 结论

本研究基于GEE平台提供的高分辨率、长时间序列的遥感产品构建RSEI，分析嘉陵江流域生态环境质量的动态变化，并通过因子探测器阐明了不同因子对生态环境质量演变的影响。研究结果发现，2001—2022年嘉陵江流域生态环境质量有明显空间差异和时序变化。西北部的白水江、北部嘉陵江干流流域大部分年份的生态环境质量好；涪江流域东南部、东河/西河流域、渠江流域则相对较差。整体上，研究区域的RSEI年均值为0.22—0.69，呈缓慢增加趋势。流域多年的生态环境质量呈现出波动变化。轻微改善区域面积占比最大(74.42%)，其次是明显改善区域(13.50%)，二者主要分布在嘉陵江海拔较高的北部和东部。严重退化区域和轻微退化区域较少，主要分布在涪江流域、东河/西河流域的西南部，以及白河流域西南部、渠江流域部分区域。整体上，嘉陵江流域的生态环境质量呈现向好的趋势，且流域大部分区域未来变化与过去的变化趋势一致，但在涪江流域、白水江流域、渠江流域部分区域以及嘉陵江流域南部的部分地区，未来变化可能与过去趋势不一致。嘉陵江流域生态环境质量的主要驱动因子为自然因子，即年降水量、年均温、海拔，但社会因子也是不可忽视的重要驱动力，尤其是人口密度较高的区域。这些自然因子和社会因子之间具有协同增强作用。