0 引言

植被是指地球表面某一地区所覆盖的植物群落。作为生态系统的核心组成部分，植被是联结大气、水、生物、岩石、土壤等自然环境要素的纽带，与物候变化、水循环、净化空气和能量流动等自然生态现象密切相关^[1-3]。因此，掌握植被时空演变规律及驱动机制是科学认识生态系统的有效途径^[4]。

归一化植被指数(Normalized Differential Vegetation Index, NDVI)是表征植被覆盖和生长状况的最佳指标，可表征植被的种类、长势和密度等特点，并已被广泛应用于植被变化及驱动力特性的研究^[5-7]。近年来，为应对物候变化，提高生态系统适应力，许多学者对植被NDVI演变及其驱动力进行了研究。Zoungrana等^[8]以基纳法索、西非为研究区，分析了植被覆盖度的影响因子，发现可达性、气候和地形等是植被NDVI的主要影响因子。Lamchin等^[9]以亚洲为研究区，借助相关分析法探究植被覆盖度变化与气候的相关性，发现气温是影响植被绿度的重要因子。高江波等^[10]以中国为研究对象，利用加权回归法分析植被NDVI与气候变化的响应规律，发现植被NDVI与气温呈现非平稳关系，降水影响显著区为青藏高原。李双成等^[11]以青藏高原为研究区，借助小波分析法分析植被NDVI与地形因子的相关性，指出人类活动是植被NDVI的制约因子。孔冬冬等^[12]也以青藏高原为研究区，利用回归分析探究植被NDVI的物候变化，认为温度是影响植被NDVI变化的重要因子。彭文甫等^[13]以四川为研究区，利用地理探测器探测影响植被覆盖度的自然因子，认为土壤、高程和温度对植被可变性有较强解释力。祝聪等^[14]以岷江上游为研究区，利用趋势分析法和地理探测器分析植被覆盖度的时空演变规律及驱动力特性，发现该区植被覆盖度良好、稳定，而海拔、气温、土壤类型和降水为主要解释因子。尽管已有研究对认识植被NDVI演变特性及驱动力有重要意义，但仍然存在以下问题：(1)对于研究手段，此前研究多选取回归、相关分析等线性计量手段分析植被NDVI的驱动力，但实际中植被的响应过程十分复杂，很难存在严格的线性关系，因此线性分析难以量化自然因子对植被NDVI的影响；(2)对于驱动因子，此前多选取气候、人类活动和地形等单因子作为植被NDVI的驱动因子，但植被受复杂的自然环境影响，驱动因子并非独立作用，因此单因子分析难以解释复杂的因子影响机制；(3)对于研究区域，此前研究区域多选大陆内部、流域、高原等，而针对海陆交互带的研究较少。

广西北部湾位于广西壮族自治区的东南部，是我国西部陆海新通道的重要关口。该地季风气候显著，地形独特，植被类型多样。近年来，西江经济带、北部湾经济区等地区的经济建设加剧了植被砍伐、生态环境恶化的状况。广西北部湾地处复杂的陆海交错地貌区，加上人为干扰日渐增强，植被覆盖率易受生态环境恶化的影响而下降。因此本研究选取广西北部湾作为研究区域，利用一元线性回归法与变异系数法对植被NDVI的时空演变格局进行分析，借助地理探测器探测其多因子影响下的影响特性，以便为该地的生态建设提供科学有效的参考依据。

1 数据来源与方法 1.1 研究区域概况

广西北部湾地处广西东南部，地理位置介于21°24′-22°43′N，107°27′-109°52′ E之间(图 1)。气候类型属亚热带季风气候，水热充足，年平均温度约22℃，年均降水量约1 600 mm。地势北高南低，地形以平原为主，东北、西北方分别横贯六万大山山脉及十万大山山脉，平均海拔为800-1 000 m。在气候和地形等要素的影响下，植被类型以亚热带季风常绿阔叶林、热带雨林为主，四季常绿，无明显的换季落叶现象。

1.2 数据来源

植被NDVI的变化与高程、坡向^[15]、气候^[16]、土地利用变化、土壤类型^[17]等因素有关，而广西北部湾地处陆海交界，有独特的植被类型、地貌特征。因此，本研究依据系统性、典型性、动态性、科学性、可量化和可获取等原则，选取高程、坡度、坡向、气温、降水、土壤类型、地貌类型、植被类型和土地利用等9个因子作为研究植被NDVI的主要驱动因子。高程、坡度、坡向数据来源于90 m分辨率的中国海拔高度(DEM)空间分布数据。气温、降水数据是基于全国主要气象站点数据，经反距离加权平均内插和DEM校正获得的栅格数据；土地利用数据是基于Landsat 8遥感影像，通过人工目视解译生成。基于SPOT/VEGETATION NDVI卫星遥感数据，通过拼接镶嵌、投影变换获得每月的植被NDVI数据，最终采用最大值合成法将月数据合成得到年植被NDVI的数据。所有数据均依照研究区矢量边界裁剪。表 1展示了主要研究数据的来源。

1.3 方法 1.3.1 变异系数法

变异系数法是衡量某对象相对变化(波动)的统计方法，通过逐像元的植被NDVI变异系数来探究其稳定性，其计算公式如下^[18]:

$ {C_v} = {\rm{ }}\frac{{\partial \left( {{\rm{NDVI}}} \right)}}{{m\left( {{\rm{NDVI}}} \right)}}, $

(1)

式中：C_v为变异系数，∂(NDVI)为植被NDVI的标准差，m (NDVI)为植被NDVI均值。

1.3.2 趋势分析法

选用最小二乘法中的一元线性回归模型，对植被NDVI的年变化进行趋势分析，计算公式为^{[19, 20]}

$ S = \frac{{n\sum\limits_{i = 1}^n {\left( {i{v_i}} \right)} - \sum\limits_{i = 1}^n i \sum\limits_{i = 1}^n {{v_i}} }}{{n\sum\limits_{i = 1}^n {{i^2}} - {{\left( {\sum\limits_{i = 1}^n i } \right)}^2}}}, $

(2)

式中：S为植被NDVI变化趋势斜率，当S为正值时，表明植被NDVI呈上升趋势，负值则为下降趋势；n为研究总年数；v_i为第i年的NDVI (i=1，2，…，15)。采用F检验对趋势结果进行显著性检验，计算公式如下^[21]:

$ \left\{ {\begin{array}{*{20}{l}} {F = U\frac{{n - 2}}{Q}, }\\ {U = {{\sum\limits_{i = 1}^n {\left( {{{\hat v}_i} - {{\bar v}_i}} \right)} }^2}, }\\ {Q = \sum\limits_{i = 1}^n {{{\left( {{v_i} - {{\hat v}_i}} \right)}^2}} , } \end{array}} \right. $

(3)

式中：F为检验统计量；U为误差平方和；Q为回归平方和；v_i为i年植被NDVI的平均值；${\hat v_i}$为i年植被NDVI回归值。根据变化趋势斜率S的大小将植被NDVI的变化趋势分为5种类型，如表 2所示。

1.3.3 地理探测器模型

地理探测器是基于空间分异角度统计分析驱动机制的方法，包括因子探测器、生态探测器、风险探测器、交互探测器4个探测内容^[22]。因子探测器是将因子分类区的方差和植被NDVI区域的总方差进行空间差异性对比分析，以探测不同因子对植被NDVI的解释力，其计算公式如下^[23]：

$ P = 1 - \frac{{\sum\limits_{h = 1}^L {{N_h}\partial _h^2} }}{{N{\partial ^2}}} = 1 - \frac{{{\rm{SSW}}}}{{{\rm{SST}}}}, $

(4)

式中：P为因子对NDVI的解释力，值越大表示该因子的解释力越强；h=1，2，…，L，表示植被NDVI或因子的分类；N_h、N分别为分类h和总区域的区域单元数；∂_h²、∂²分别是分类h和总区域的植被NDVI的方差；SSW、SST分别为分类区域的方差之和与区域总方差。

生态探测器用于比较两因子间对植被NDVI空间分布的影响是否存在显著差异，以F统计量来衡量：

$ F = \frac{{{N_{{x_i}}} \times \left( {{N_{{x_j}}} - 1} \right) \times \sum\limits_{h = 1}^{{L_i}} {{N_h}\partial _h^2} }}{{{N_{{x_j}}} \times \left( {{N_{{x_i}}} - 1} \right) \times \sum\limits_{h = 1}^{{L_j}} {{N_h}\partial _h^2} }}, $

(5)

式中：N_xi、N_xj分别表示两个不同因子(x_i和x_j)的样本数量；L_i、L_j分别表示变量x_i和x_j分类数。

交互探测器是用于比较因子交互前后解释P值的大小，用以判别因子共同作用后对植被NDVI的解释力是增强还是减弱，同时也用于判别因子交互的关系(线性/非线性)。

风险区探测器是判断两因子间的属性均值是否有显著差别，用于搜索植被NDVI最佳区域的探测。用t统计量来检验：

$ t = \frac{{{{\bar Y}_{h = 1}} - {{\bar Y}_{h = 2}}}}{{{{\left[ {\frac{{{\rm{Var}}\left( {{Y_{h = 1}}} \right)}}{{{n_h} = 1}} + \frac{{{\rm{Var}}\left( {{Y_{h = 2}}} \right)}}{{{n_h} = 2}}} \right]}^{\frac{1}{2}}}}}, $

(6)

式中：Y_h为分类h内植被NDVI值，Y_h为分类h内植被NDVI均值，n_h为子区域h内样本数量，Var表示方差。

基于文献[24]，选用自然间断法作为因子的分类方法。将气温、降水分为10级，高程、坡度、坡向分为9级。土壤类型、地貌类型、植被类型、土地利用类型则是在原来的分类基础上，经过聚类合并将二级别分类合并成一级别分类，减少分类级数。最终将土壤类型分10类，地貌类型分为5类，植被类型分为9类，土地利用分为6类。由于地理探测器的数据容纳量最多约为30 000，故借助ArcGIS软件，按照3 km×3 km的网格，在植被NDVI栅格数据中随机提取2 245个采样点。把植被NDVI采样点图层与因子分类图层叠加，获取采样数据与分类数据的空间关联属性表，最后将属性表导入地理探测器运行。

2 结果与分析 2.1 植被NDVI的时空演变

2000-2015年广西北部湾地区植被NDVI年季变化统计分析表明，16年间研究区植被生长状况整体良好(图 2)。植被NDVI年均值为0.753，整体以0.006 1的增长率缓慢上升。2000-2005年植被NDVI上升较快，2005年之后上升速度趋于平缓。从季节分配来看，研究区植被NDVI有明显季节差异，均值大小排序为秋季＞夏季＞春季＞冬季。夏季和秋季的植被NDVI值为0.68-0.75，且较平稳；冬季和春季的植被NDVI值为0.5-0.75，有波动上升的趋势。以上结果表明，研究区植被NDVI在16年间不断增加，且全年常绿的植被类型数量可能逐渐增多。

将植被NDVI类型划分为低、中低、中、中高、高5个等级(表 3)。各等级的面积排序为中高＞中＞高＞中低＞低。2000年和2015年研究区的植被覆盖类型分布表现出相同的趋势，两者均以中高等级为主，面积分别达17 499.52 km²和12 339 km²，分别占总面积的87.97%和62.20%；中、高等级的面积次之；而低、中低等级的面积占比均低于5%。同时，在16年间，中高等级的面积显著下降，减少面积为5 110.52 km²，降比达25.77%；而中、高等级的面积显著增加，增加面积为1 500-2 700 km²，增幅为8%-14%；低、中低等级的面积仅有小幅度上升，且比率小于4%。

从空间演变趋势分析(图 3a)，植被NDVI增加(S＞0)的区域占总面积的55%，减少(S＜0)的区域占45%，增加的面积比减少的面积多10%，整体以增加为主。其中，显著增加和极显著增加面积占54.32%，主要分布在中部地区以及高山地区; 显著减少和极显著减少面积占39.34%，主要分布在沿海地区，最为明显的是北海的铁山港区及合浦县沿岸开发区。不显著变化面积占总面积的6.33%，分布在早年开发的沿海海岸区域。以上结果表明, 植被NDVI呈沿海较内陆地区退化更明显的状况。

从变异角度分析(图 3b)，研究区植被NDVI整体较稳定，植被NDVI平均变异系数为3.91%。变异系数低于15%的稳定类型占比为48.55%，主要分布在高、中高等级植被覆盖地区；变异系数为15%-40%的不稳定类型占比为45.34%，主要分布在中等级植被覆盖区域；变异系数高于40%的极不稳定类型占比为6.11%，主要分布在沿海海岸的低等级植被覆盖区。以上结果表明，植被NDVI在高等级植被覆盖区较稳定，在中、中低、低等级植被覆盖区较不稳定。

2.2 地理探测器探测结果 2.2.1 影响因子探测

基于因子探测器的研究结果(表 4)表明，各因子对植被NDVI的影响力度排序为气温＞高程＞土壤类型＞坡度＞地貌类型＞土地利用＞植被类型＞降水＞坡向。其中气温的解释力最大，P值为0.381，解释力超35%。高程、土壤类型的解释力次之，P值分别为0.306，0.298，解释力在29%以上。坡向和降水的解释力最弱，P值分别为0.065，0.109，解释力均小于11%。

2000-2015年各影响因子的P值年际变化(图 4)表明，植被类型、坡度、降水及土地利用的P值呈增加趋势，气温与高程的P值呈下降趋势，土壤类型、地貌类型、坡向的P值则几乎不变。2000-2005年土地利用和降水的P值呈上升趋势，其余各因子的P值均呈下降趋势；2005-2010年除高程的P值呈下降趋势外，其余各因子的P值均呈上升趋势；2010-2015年，除坡向的P值稍有上升趋势外，其余各因子的P值均呈减少趋势。16年间，P值增加和下降最明显的因子分别是土地利用和高程，表明人类活动对植被NDVI的影响越来越显著。

2.2.2 生态探测

生态探测器分析表明, 气温与其他因子对植被NDVI的影响存在显著差异(表 5)；除坡向与降水对植被NDVI的影响无显著差异外，坡向与其他因子均对植被NDVI存在显著差异。坡度与高程、地貌类型、土壤类型对植被NDVI的影响无显著差异，与其余因子对植被NDVI的影响存在显著差异。土地利用与所有因子对植被NDVI的影响均存在显著差异。高程与坡度、地貌类型及土壤类型对植被NDVI的影响无显著差异，与其余因子对植被NDVI的影响存在显著差异。地貌类型与坡度、高程以及土壤类型对植被NDVI的影响无显著差异外，与其余因子对植被NDVI的影响均存在显著差异。植被类型除与降水对植被NDVI的影响无显著差异外，与其他因子对植被NDVI的影响均存在显著差异。土壤类型与坡度、高程及地貌类型对植被NDVI影响无显著差异，与其余因子对植被NDVI的影响均存在显著差异。降水与坡向、植被类型对植被NDVI的影响无显著差异，与其余因子对植被NDVI的影响均存在显著差异。结合因子探测的分析结果，筛选得出影响植被NDVI变化的主要解释因子是气温。

2.2.3 因子交互探测

交互探测器分析表明(表 6)，因子之间交互协同关系存在一定的差异，交互作用均能加强单因子对植被NDVI的影响。各因子之间存在线性及非线性交互关系，无独立作用因子。其中，以坡向为基准，与气温产生最大影响力(0.423)，与土壤类型产生的影响力次之(0.374)；以坡度为基准，与土壤类型产生最大影响力(0.475)，与气温产生的影响力次之(0.466)；以土地利用为基准，与气温产生最大影响力(0.475)，与高程产生的影响力次之(0.445)。排序发现因子交互作用后的解释力均比单因子独立作用的解释力大，气温与各因子交互作用后易产生较大的解释力。坡向与其他因子间交互协同作用多为非线性增强关系，其余因子间交互协同的关系多为双因子增强关系。

2.2.4 适宜性探测

风险探测器分析结果表明：同因子但不同分类对植被NDVI影响力度不同(表 7)。将气温划分为10个等级，用A1-A10表示，数字越大表明气温值越大。随着气温的变化，植被NDVI也不同。在A2等级产生了最大的植被NDVI值，为0.833，但A2等级与A1、A3等级无明显差异，表明温度在17.6-21.3℃内植被NDVI值最高，植被覆盖率较大。对A1-A10分区内的植被NDVI大小排序，A2 (0.833)>A3 (0.832)>A1 (0.831)>A4 (0.825)>A5 (0.819)>A6 (0.798)>A7 (0.778)>A8 (0.776)>A9 (0.766)>A10 (0.707)，植被NDVI值在A2等级之后随着温度的上升而不断下降，表明一定范围内的温度利于植被NDVI的增长，超出范围则致使其下降。

将高程划分为9个分区，用E1-E9表示(表 8)。不同的高程范围有不同的植被NDVI，在E9区内植被NDVI值最大，为0.838。E8区与E9区的影响无显著差异，与其余因子的影响均存在明显差异，即表明海拔在674-1 430 m能产生较大的植被NDVI值。E1-E9区产生的植被NDVI大小排序为E1 (0.719)＜E2 (0.771)＜E3 (0.787)＜E4 (0.797)＜E5 (0.822)＜E6 (0.828)＜E7 (0.831)＜E8 (0.833)＜E9 (0.838)，表明在0-1 430 m随着海拔的升高，植被NDVI值也在上升。原因可能是海拔稍高的地区，开发力度小，生态环境保护得好，人类活动干扰较少，植被乱砍滥伐的现象较少，植被覆盖度也较高。

将土壤类型划分为10个分区，用S1-S10表示(表 9)，不同类型的土壤存在不同的植被NDVI值。在S9区能产生最大的植被NDVI值，为0.833。S8区与S9区的影响无明显差异，这表明黄色赤红壤和黄红壤利于植被NDVI值的增加。原因可能是北部湾海岸带水热充足，水稻为主要粮食作物，水稻土中常年种植水稻，因而该类土壤的植被覆盖率高。黄色赤红壤、黄红壤是研究区内分布较为广泛、肥力也相对优质的土壤类型，因而也利于植被NDVI值的增长。

在其余影响因子(植被类型、坡度、地貌类型和土地利用)中，植被类型分类以亚热带季风常绿阔叶林分布区产生的植被NDVI最高，为0.832-0.835;37.53-70.34坡度产生最高的植被NDVI，为0.833;地貌类型以中/大起伏山地中产生最高植被NDVI，为0.809-0.816;土地利用则是在林地中产生最高植被NDVI，为0.785。

3 讨论

本研究利用对比分析、回归分析等方法对广西北部湾的植被NDVI演变进行探究，从时空尺度挖掘其变化特点。结果发现广西北部湾地区的植被NDVI类型以中高覆盖度为主，高山地区植被覆盖度高，沿海地区植被覆盖度低。从变异视角看该区域植被覆盖度稳定性能较好，高山地区较稳定，内陆沿海地区稳定性较差。田义超等^[25]用标准差法探究广西北部湾植被NDVI，发现其呈现出“两头高、中心地带低”的分布规律，这与本研究的结论高度一致。同时，本研究结果与广西北部湾发展现状较符合，中部高山地带，交通不便，难以开发，生态保护较好，植被NDVI亦较好。沿海港口城市的开发力度大，植被NDVI较低。在一带一路倡议下，近年来广西北部湾被大力开发，是高植被NDVI增加缓慢的一个重要原因。同时，本研究发现在一些开发建设较大的县区内部区域出现植被NDVI变异较大的情况，这也表明本研究能合理反映出研究区域的植被覆盖度的时空变化情况。

地理探测器是一种能较好地识别驱动力以及影响机制的计量方法。本研究引入地理探测器以辨识出最具解释力的影响因子，最终筛选出气温为主要影响因子。从年际变化看，土地利用解释力明显上升，表明人类因子的影响力度逐渐增加。该结论与成方妍等^[26]使用hurst指数对广西北部湾地带植被NDVI变化驱动因子分析所提出的，地形和人为因素影响是主要影响因子的结论高度契合。交互探测器探测表明，因子之间能相互协同作用，并加强对植被NDVI的解释力，但是交互增强机制存在复杂的非线性关系，因而单纯的线性计量手段难以解释影响因子之间的交互关系，难以反映因子间的交互机制。另外，风险探测器已探测出各因子中较利于植被NDVI增加的分类。以上研究均能为广西北部湾生态保护和建设提供科学参考。

4 结论

本研究采用一元线性回归法以及变异系数法对2000-2015年广西北部湾植被NDVI时空变化及稳定性进行评价；借助因子探测器、生态探测器辨识植被NDVI的主要解释因子，并分析各因子解释力的时序变化特点。利用交互探测器探测因子两两交互作用的协同规律。最后，通过风险探测器筛选出植被NDVI最佳区域，甄别最适宜植被NDVI增长的因子分类。主要结论如下：

2000-2015年研究区植被覆盖状况良好，植被NDVI年均值为0.753，并呈现夏季和秋季高，冬季和春季低的现象。在16年间，植被NDVI有缓慢上升的趋势，上升的幅度较小，上升率仅为0.006 1；植被NDVI增加面积占比55%，下降占比45%，上升下降的面积占比几乎持平。从稳定性看，平均变异系数为3.91%，变异稳定区域占总面积的48.55%，多分布在高、中高等级植被覆盖区。变异极不稳定区占比为6.11%，多分布在沿海的低等级植被覆盖区，内陆的植被NDVI较沿海的更稳定。

2000-2015年研究区中高等级的植被覆盖类型面积占比大，2000年面积高达17 499.52 km²，2015年仍有12 339 km²。但有往中、中低等级的植被类型变化的趋势，降幅较大，降比达25.77%。高、中、中低、低等级这4种覆盖类型面积占比较小，但有缓慢上升的趋势。植被NDVI空间分布状况呈高山地区植被覆盖度高，沿海港口地区植被NDVI低，中部以中高植被覆盖类型为主。

因子探测器探测发现气温是解释力最大的因子，解释力超35%。高程、土壤类型的解释力次之，解释力分别为30.6%、29.8%。基于16年间的变化趋势，气温的解释力在波动下降，而土地利用的解释力在不断上升，这表明人为因子的解释力在上升，人类活动的干扰在加剧；交互探测器探测发现因子两两交互作用后均能增强单因子解释力，无独立作用的因子，非线性协同作用强，增强机制并非单纯地叠加或线性增强。利用风险探测器探测出植被NDVI最佳区域，搜索出最适合植被NDVI增长的因子分类。此后可依据最适宜分类进行适当人为干预，以便获得较高的植被覆盖率。