GEE集成了Landsat、MODIS、Sentinel等30余种卫星数据资源，在生态系统动态分析、农业生产力预测及生物多样性保护等跨学科领域展现显著应用价值^[13]。本研究中利用谷歌地球搜索引擎进行石漠化等级划分、植被覆盖度的泰尔-森(Theil-Sen，T-S)趋势分析和曼-肯德尔(Mann-Kendall，M-K)显著性检验。

本研究涉及数据包括：数字高程(DEM)数据、归一化植被指数(NDVI)数据、自然因子数据、社会经济因子数据、研究区行政边界(表1)。气温、降水和归一化植被指数均来源于逐月栅格数据。气温和降水采用全年数据的单日平均作为年值，归一化植被指数则取每年12个月的月最大值。在此基础上，提取以上数据在重庆市各个区(县)范围内的栅格平均值，得到重庆市38个区(县)21年间的逐年数据。其他数据则直接以年为单位按区(县)行政边界求平均值处理后汇总，路网密度根据《重庆市统计年鉴》各区(县)数据直接得到。

通过GEE在线编写代码调用MODIS卫星数据集，参考文献[14]，使用T-S斜率分析计算时间序列中所有数据对之间斜率的中位数，得到β，分析研究区归一化植被指数变化趋势。β＞0表示时间序列存在上升趋势，β＜0则相反。

M-K显著性检验^[15]则用于判断变化趋势的统计显著性。该方法通过计算标准化检验统计量Z来评估趋势的置信水平。在本研究中，Z值用于判断β所指示的归一化植被指数增减趋势是否具有统计学意义，并根据|Z|与特定阈值(1.65、1.96、2.58)的比较，将变化趋势划分为不显著、微显著、显著和极显著等不同等级^[16]。通过结合β值与Z值(表2)，可实现对研究区植被变化趋势及其空间分异特征的全面评估。

在生态系统的非线性响应机制研究中，由于环境因子与生物学变量常呈现非单调响应、阈值效应等复杂关系，传统参数化线性模型往往难以准确解析此类非线性过程。GAM作为半参数化统计工具^[17−19]，通过引入光滑函数建立响应变量与预测因子的柔性关系式，可以用来解释变量之间的非线性关系。

本研究基于R4.3.3平台，利用GAM探索植被变化趋势驱动机制的非线性关系，采用mgcv包^[20]实现。将研究区38个区(县)的逐年归一化植被指数作为响应变量，土层厚度、年平均气温、年平均降水量3个自然因子和国内生产总值、路网密度、人口密度3个社会经济因子作为解释变量，构建GAM模型，绘制响应曲线。

GAM模型依靠参考自由度、决定系数(R²)、方差解释率、F统计值和P值来判断解释变量对归一化植被指数的影响程度^[21]。其中参考自由度用来判断是否为非线性关系，参考自由度为1时说明解释变量与响应变量为线性关系，当参考自由度＞1时解释变量与响应变量为非线性关系，且参考自由度的值越大，解释变量与响应变量的非线性关系越显著。R²的取值范围为0~1，R²越接近于1则表明模型越精确。方差解释率越高，表明拟合效果越优。F值用来判断影响因子的重要性程度，F值越大，表明其相对重要性越大。P值用来判断显著性，P值越小则结果越显著。

通过GEE在线编写代码，调用Landsat数据，参考侯自航^[22]、LIAO等^[23]的研究，结合重庆市地理特点，选择植被覆盖度、岩石裸露率、坡度作为指标计算石漠化指数(quaternary rocky desertification index，QRI)，并以此进行石漠化分级(表3)。植被覆盖度是反映地表覆盖情况的重要指标，与石漠化有着较显著的相关性^[22]；岩石裸露率表征基岩出露状况，是石漠化的直接视觉标志；坡度则影响水土保持能力和土壤侵蚀速率，是石漠化发生发展的重要地形条件^[24]。三者综合能够较全面地刻画喀斯特地区的石漠化程度。

在研究期间，研究区植被整体呈现良好恢复态势，其中西部和东北部部分区域植被恢复最为显著，东南部地区植被恢复程度相对较弱。植被退化区域呈零星分布，中西部存在小范围连片退化区。由图1可知：研究区内93.00%的面积区域内归一化植被指数呈增加趋势，呈现下降趋势的仅为7.00%。根据显著性划分结果判断，69.60%的面积区域内归一化植被指数具有极显著上升或显著上升趋势。这一结果表明，过去21 a间重庆市植被恢复成效显著。

图  1  重庆市植被变化趋势及显著性面积占比

Figure 1.  Trends of vegetation changes and proportion of significant areas in Chongqing

基于GAM进行了各类参数的假设性检验(表4)。在6个驱动因子中，仅石漠化地区的年平均降水量对归一化植被指数无显著影响，其余驱动因子均在P＜0.001水平上对归一化植被指数影响显著。6个驱动因子中，年平均降水量的参考自由度等于1，其与归一化植被指数之间为线性关系，但因GAM是对广义线性模型(aeneralized linear models，GLM)的推广^[25]，仍能利用GAM模型结果进行解释与分析。社会经济因子中，人口密度的方差解释率较高(88.9%)，R²较大，为0.881，作为归一化植被指数的预测变量拟合结果较优。自然因子中，非石漠化地区的土层厚度和年平均气温方差解释率较高，分别为91.3%和50.0%，R²较大，分别为0.912和0.493，作为归一化植被指数的预测变量拟合结果较优。

由图2A可知：国内生产总值对归一植被化指数变化的整体影响较小，在国内生产总值较低时归一化植被指数呈先小幅度增长再缓慢下降趋势，当国内生产总值超过约600亿元后，曲线趋于平缓。由图2B可以看出：路网密度增加，归一化植被指数逐渐减少，呈现一定同步性，反映出路网密度对植被覆盖存在一定负面影响。由图2C可知：曲线呈波动起伏趋势。在人口密度相对较低时，归一化植被指数随人口增长大幅减少；当人口密度超过约1 800人·km⁻²后，归一化植被指数开始随着人口增长小幅度上升；在约2 500人·km⁻²时到达峰值，归一化植被指数约为0.68，之后逐渐下降。

图  2  归一化植被指数-社会经济因子响应曲线

Figure 2.  NDVI-socioeconomic factors response curve

由图3可知：归一化植被指数随着年平均降水量的增加缓慢增长，降水促进植被生长的作用较弱。其中，石漠化地区归一化植被指数增加速率低于非石漠化地区，可能与石漠化地区土层保水持水能力弱有关。随着年平均气温的逐渐升高，归一化植被指数呈下降趋势，但在不同气温区间下降速率不同(图4)。石漠化地区归一化植被指数下降速率先急后缓，可能是由于石漠化地区气温升高对植被产生如高温胁迫、光合效率降低等负面影响^[26]，加之缺少有效水分的补给，对植被生长产生明显负面影响，归一化植被指数先快速下降；但因石漠化地区植被长期生长在逆境下，较强的抗逆性使其能够逐渐适应变化的环境，年平均气温进一步升高后归一化植被指数下降速率变缓。非石漠化地区归一化植被指数随年平均气温的升高，下降速率先缓后急。气温升高后，非石漠化地区植被能够依靠充足的水分产生补偿作用，减少高温毒害，但超过一定范围后补偿作用失效归一化植被指数急剧下降。

图  3  归一化植被指数-年平均降水量响应曲线

Figure 3.  NDVI-mean annual precipitation response curve

图  4  归一化植被指数-年平均气温响应曲线

Figure 4.  NDVI-mean annual temperature response curve

由图5可知，归一化植被指数均随土层厚度增加呈显著的增加趋势。相较于石漠化地区，非石漠化地区曲线变化幅度更大，其归一化植被指数最大值更高，最小值更低。

图  5  归一化植被指数-土层厚度响应曲线

Figure 5.  NDVI-soil depth response curve

综合图2~5可以发现：石漠化地区和非石漠化地区归一植被指数的总体变化趋势是一致的，但波动幅度存在一定差异。相较于石漠化地区，非石漠化地区的植被对3类自然驱动因子的响应都更加显著和敏感。驱动植被变化的机制在石漠化和非石漠化地区存在差异，这解释了2.1节中的空间异质性。同时，这个结论强烈暗示了石漠化本身是一个关键的中介力量，石漠化改变了驱动因子(主要是自然因子)的作用方式和强度，因此我们在此基础上直接研究石漠化与植被的关系。

由表5可知：2000—2020年重庆市石漠化(轻度、中度、强度、极强度)面积相对较小(占比5.7%~11.0%)，且不断下降，其中轻度石漠化、中度石漠化呈波动下降趋势。强度石漠化和极强度石漠化面积虽有所上升，但始终维持在较低水平。

结合归一化植被指数变化与石漠化动态变化可以发现，21世纪初研究区西部以轻度石漠化为主，之后面积减少迅速，其归一化植被指数也出现了一定的增长，但主城区域存在归一化植被指数下降不显著的情况。研究区东北部和东南部石漠化面积最多。经过20 a的治理后，东南部的石漠化状况得到显著改善，大面积连续石漠化基本消失，其归一化植被指数呈现微显著和显著增加趋势；东北部的石漠化状况也得到了改善，同时其归一化植被指数呈现极显著增加，但依旧存在一定面积的连续石漠化区域，一直是重庆市石漠化最严重的区域。这可能与该区域石漠化治理难度有关。

区域尺度的归一化植被指数是植被覆盖度与植被长势共同作用的结果。一般而言，人类社会活动通过改变土地利用方式，如毁林开荒、退耕还林、城市建设等，直接且剧烈地改变植被覆盖度^[27−28]，而自然因子通过调节植物的生理过程(光合作用、蒸腾作用等)影响植被的长势和生产力^[29]。本研究以归一化植被指数作为衡量植被变化趋势的指标。

本研究选取了3个社会经济因子：国内生产总值、路网密度、人口密度。国内生产总值与经济发展、城镇化率存在着较高的关联度，路网密度在一定程度上表征着城市活力^[30]，人口密度能够直接代表人类活动的强度及其作用于其他因素的强度。

重庆市统计局公布的数据表明：2000—2020年重庆市第一产业占比不断下降，重庆市经济发展对自然资源的消耗与依赖较小，因此其国内生产总值与归一化植被指数的关联较弱。重庆市是中国中西部的交通枢纽，在经济发展过程中建设了各种交通设施工程^[31]，不可避免地需要占据一部分森林植被面积，对归一化植被指数生产了一定的负面影响，但未起到决定性作用。人口密度关系着社会发展和资源环境承载力之间的平衡，对植被覆盖变化具有重要意义。人口密度增长主要集中于城市建成区。随着人口密度的增长，城市化进程不断加快，住房、医疗、教育、工作等居民生存和发展所需要的资源不断增加，农业扩张、城市建设等占据了生态用地^[32−33]。这些会直接导致植被覆盖度大幅降低。本研究中，在人口高密度城市建成区，归一化植被指数维持在较低但稳定的水平。这可以解释为：①植被覆盖度已降至城市环境下的“底线”，可被侵占的自然植被已很少；②此时的归一化植被指数主要反映的是城市绿地(人工植被)的长势，其受到人为养护的支持，受其他因素影响较小。

城市建成区内植被动态主要受城市规划、人口流动与人为管护等社会活动主导^[34]，对自然因子的变化不敏感，故下面对自然因子的讨论将不适用城市建成区。

本研究中，重庆市气温常年较高，植物生长所需要的热量较为充足，但高温胁迫也会抑制植被长势^[35]。石漠化地区风化层浅薄，失水率增加^[36]，即使降雨量增加，能够被植被吸收利用的水分也较少，增加的水分难以持续有效地供给植物生长所需，导致降水促进作用有限。

土层厚度是决定植被生产力的重要因素，它影响着土壤水分、养分的储量和植被根系分布的空间范围^[37]，决定了植被的潜在长势上限。非石漠化地区响应更敏感。在山地条件下，土层厚度通常与地形、坡度、坡向和坡位有关^[38]。因坡度原因，石漠化地区土层厚度常常分布不均匀，造成植被覆盖类型和面积等在不同坡位的分异，导致土层厚度对归一化植被指数变化的解释效果较差。

非石漠化地区归一化植被指数基数较高。归一化植被指数的饱和效应可能导致对植被变化趋势的评估与实际状况有所偏离。未来研究可结合对饱和效应不敏感的指数，如增强植被指数，进行交叉验证。

在石漠化区，驱动因子首先影响水土流失等石漠化过程，再通过改变生境条件(土壤、水分等)，最终影响植被^[38]。在非石漠化地区，驱动因子更直接地作用于植被。石漠化状况就是一个“滤镜”或“开关”，它重塑了自然驱动因子对植被的作用格局。生态工程的成功，其本质是通过人为干预打破了石漠化的恶性循环，使系统进入“植被恢复—改善石漠化—生境优化—进一步促植被恢复”的良性循环^[39−42]。植被变化和石漠化两者实现协同演化。

在促进植被恢复的过程中，应在宏观政策指引下，对石漠化地区采取“顺应自然，弥补短板”的策略，如集水、采肥等，优先做好石漠化治理，改善立地条件，提升植被覆盖度；对非石漠化地区，则可采取“优化结构，提升功能”的策略，提高植被生产力与质量。

研究区内93%的区域归一化植被指数呈增加趋势，植被状况总体改善，但空间异质性明显，东部地区植被恢复状况整体比西部好。社会经济因子中，国内生产总值对植被覆盖的影响较弱；路网密度增加与植被覆盖减少呈现同步性，但不起决定作用；人口密度增长初期对植被覆盖呈现显著负面影响，超过某一阈值负面影响降低，对植被覆盖变化具有重要意义。自然因子中，降水增加促进植被缓慢恢复，气温升高对植被生长呈负面影响，土层厚度的增加明显改善植被生长状况。相较石漠化地区，非石漠化地区的植被对3类自然驱动因子的响应更加显著和敏感，且气温和土层厚度能够较好地解释植被覆盖变化。研究区石漠化整体呈下降趋势，石漠化面积减少明显的区域，植被覆盖增加趋势也较明显，两者具有一定的协同演化关系。在进行生态恢复与管理时，须考虑石漠化的中介作用，实施分区、分类的精准施策。