基于地理探测器模型的中亚NDVI时空变化特征及其驱动因子分析

王 伟，阿里木·赛买提，吉力力·阿不都外力

(1.中国科学院新疆生态与地理研究所荒漠与绿洲生态国家重点实验室，乌鲁木齐 830011；2.中国科学院大学，北京 100049；3.中国科学院中亚生态与环境研究中心，乌鲁木齐 830011)

0 引言

在全球气候变化的大背景下，中亚地区从20世纪70年代起气温开始迅速上升，平均增温速率达到0.4℃/10 a，远高于同期全球变暖速率[1-2]。与此同时，不合理的土地开发利用方式[3]、粗犷的农田灌溉模式[4]等人类活动的不断增强更加剧了区域生态与资源的竞争局面。于是，中亚地区开始出现了土壤盐渍化[5]、湿地退化[6]和湖泊萎缩[7]等一系列生态环境问题。植被作为联结大气圈、土壤圈、水圈和生物圈的重要纽带，对生态环境变化具有极高的敏感性[8]。因此研究中亚地区的植被时空变化特征及其对气候变化和人类活动的响应，对进一步研究中亚地区陆表环境变化过程及区域生态环境保护具有重要意义。

近几年不少学者通过经验正交模型(empirical orthogonal function，EOF)[9]、变异系数(coefficient of variation，CV)[10]、Hurst指数[10]和趋势分析[11]等方法研究了中亚地区的植被空间分布特征和变化趋势，通过相关分析[9]、偏相关分析[12]和奇异值分解(singular value decomposition，SVD)[12]等方法分析不同空间和时间尺度上中亚地区植被覆盖变化对气候变化的响应特征。但前人研究大多集中于气温和降水这2个气候因子对植被变化的影响，忽略了其他气候因素以及地形、土壤等多种因素的作用。由于中亚地区地形较为复杂，不同地形区域的气候因素对植被生长的影响可能完全不同，比如在中亚东南部山区，年均气温与植被变化呈显著的正相关[13]，而在中亚西南部地区气候转暖趋势则加剧了当地的干旱[9]。近些年来，人类活动对土地利用与覆被变化的影响加剧，尤其是1991年苏联解体以来，更是发生了显著变化[14]。这些因素无一不影响着中亚地区植被覆盖的时空变化特征，除此之外，这些因素之间的交互作用往往也会增强其对植被覆盖变化的影响。

有关自然和人为因子对中亚植被覆被特征的影响方面，前人已经做了大量的工作，然而还有一些问题仍然没有得到有效解决。例如，对中亚地区植被覆盖变化与驱动因子的研究应该不仅是从两者相关系数空间分布的角度出发，而且更要理清驱动因子与植被覆盖变化之间的规律与内在机理，这样才有助于全面了解中亚乃至整个亚洲中部干旱区植被覆盖时空变化的趋势和未来生态环境演化的方向。应用王劲峰等[15]提出的地理探测器模型不仅可以弥补这一不足，也可为其他地区植被时空变化特征影响因素研究提供可行的方法。本文利用分区统计和趋势分析的方法研究中亚地区1991—2015年间植被覆盖的时空变化特征，结合地理探测器模型分析气候、地形、土壤和土地覆被等多种变量的空间分异特征和变化趋势对中亚地区植被覆盖时空变化的驱动作用，并研究了不同因子之间的交互作用对植被覆盖时空变化的影响。进一步揭示气候变化和人类活动对中亚植被覆盖变化的作用机理，旨在为中亚地区生态环境建设和制定社会经济可持续发展模式提供理论基础和科学依据。

1 研究区概况及数据源

1.1 研究区概况

本文所研究的中亚范围西达里海与伏尔加河，东与中国接壤，北至额尔齐斯河的分水岭，并延伸至西西伯利亚平原南部，南临伊朗与阿富汗的边界，地理位置在N35°5′2.24″～52°33′30.49″，E46°45′28.13″～87°21′47.81″[16]之间。行政区域涵盖中亚五国，即哈萨克斯坦、吉尔吉斯斯坦、乌兹别克斯坦、塔吉克斯坦和土库曼斯坦，总面积近4×105km2，人口接近6.6×106，人口密度约17人/km2，是典型的地广人稀区域[17]。整体上中亚地势呈东南高，西北低，由于受东南部高山阻挡，来自印度洋和太平洋的水汽无法深入，这就使得中亚五国成为北半球温带面积最大的干旱区。同样受地形和气候等多种因素影响，中亚的植被分布具有明显的空间分异性。因此本文在高程数据的基础上，结合已有地形图将中亚地区划分为8大地形区(图1)，来研究不同地形区植被覆盖的变化特征。

图1 研究区域示意图Fig.1 Map of research area

1.2 数据来源与预处理

本文采用的数据如表1所示，由于数据为不同来源和不同空间分辨率，本文结合模型的运算能力和研究任务的精度要求，基于最邻近距离法对数据集进行了重采样，最后去除缺失值，共得到20 569行数据。由于地理探测器模型输入的自变量必须保证是离散型数据，因此基于R语言“GD”包的“optidisc”函数，对数据进行了离散化处理。

表1 数据来源与预处理Tab.1 Data source and data preprocessing

2 研究方法

2.1 趋势分析法

基于最小二乘方法的一元线性回归模型对植被变化进行趋势分析是研究植被动态变化的常用方法[18-19]，对1991—2015年间中亚地区逐年植被覆被变化进行逐像元趋势性分析，计算公式为[19]

(1)

式中：n为研究时段的年数；Vi为第i年的NDVI值(i=1，2，...，25)；S为植被覆被变化趋势斜率，当S为正值时，则该像元NDVI值呈上升的趋势，反之呈下降的趋势。

为了更好地表达植被覆盖的变化趋势，采用F检验对NDVI的变化趋势进行显著性检验，显著性检验结果反映了趋势变化可信度的高低，F检验统计量计算公式为[20]

(2)

表2 NDVI变化趋势类型Tab.2 Types of NDVI change trend

2.2 地理探测器模型

地理探测器主要包含以下3种探测器：

1)因子探测器，该探测器表示研究中每个影响因子变化趋势对NDVI时空变化趋势的影响大小，q值越高的影响因子，对NDVI变化的影响力越大。

2)风险探测器，探测因子对NDVI变化是否具有风险性，指示因子在不同等级范围内对NDVI时空变化的影响。

3)交互作用探测器，用于判断不同影响因子对植被时空变化的交互作用。通过比较在单因子作用时的q值、2个单因子的q值之和与双因子交互作用时的q值，根据三者之间的大小关系，将交互作用类型分为5类，如表3所示。

表3 双因子交互作用结果类型Tab.3 Types of two-factor interaction result

从上面的3个探测器可以看出，在地理探测器中，q值是衡量一个因子对因变量的贡献率或者说是解释力，其模型表达式[15]为

(3)

3 结果与分析

3.1 中亚NDVI的空间分异特征

利用栅格计算器计算得到1991—2015年间中亚地区NDVI空间分布及其变化率、土地覆被类型和土壤类型如图2所示。

(a)NDVI空间分布 (b)NDVI变化率

(c)土地覆被类型 (d)土壤类型

图2 中亚地区NDVI空间分布与变化率、土地覆被类型和土壤类型

Fig.2SpatialdistributionofNDVI,NDVIchanges,landcovertypesandsoiltypes

从图2(a)可以看出，中亚地区NDVI具有明显的空间分异性，在平原区，NDVI平行纬度呈带状分布，从南向北逐渐增加，绿洲沿河流两岸呈带状分布，河口三角洲地带普遍发育绿洲；在山区，NDVI空间分布较平原区复杂，NDVI随高程增加先升高后降低。因此可以看出中亚地区植被覆盖空间分布主要受纬度地带性、垂直地带性规律共同作用。

不同地形区和土地覆被类型对应的NDVI值变化情况如图3所示，从地形区上看(图3(a))，西西伯利亚平原和哈萨克丘陵东北部地区NDVI最高，NDVI平均值为0.551，图兰平原和于斯蒂尔盖高原最低。从土地覆被类型上看(图3(b))，林地NDVI最高，平均值为0.7，旱地次之，水域与冰川NDVI最低。

(a)不同地形区 (b)不同土地覆被类型

图3 不同地形区和土地覆被类型对应的NDVI值变化

Fig.3ChangesofNDVIvaluescorrespondingtodifferentterrainareasandlandcovertypes

3.2 中亚NDVI时空变化特征

3.2.1 中亚NDVI的时间变化特征

如图3(a)所示，1991—2015年间中亚地区NDVI的平均值为0.34，其年均NDVI最高值出现在2007年，达到0.365，最低值出现在2008年，为0.312。哈萨克丘陵和图尔盖高原地区在近5 a来呈持续上升的趋势，然而于斯蒂尔盖高原和图兰平原地区在近10 a间呈明显下降趋势。如图3(b)所示，旱地、灌丛地和荒地的NDVI值略有下降，水域与冰川的NDVI从2000年后开始持续上升。水域NDVI的升高可能归因于中亚地区部分水域面积减小，水位下降；水中悬浮物浓度增加、透明度减小、近岸浅水区水生植物增加，水质变差导致[21]。

3.2.2 中亚NDVI的空间变化特征

近25 a来中亚地区植被覆盖总体上呈轻微波动变化，但是部分地区植被覆盖变化也很明显，植被显著改善区域主要集中在哈萨克丘陵的中低海拔地区，虽然从图2(b)和(c)中可以看出，NDVI下降比较严重的区域主要分布哈萨克斯坦北部旱地区域。但如图4所示的中亚地区植被覆盖变化趋势(经过显著性检验)，哈萨克斯坦北部旱地NDVI下降趋势并不明显。该原因可能是因为1991年苏联解体以来，农田出现弃耕后复垦的现象，但是复垦面积仍然无法达到弃耕前的水平[22]。植被显著恶化区域主要分布在咸海盆地以南地区，以西南方向最为集中。可能是随着咸海水面退缩后，在干涸湖盆底部形成了大量的富盐细粒沉积物，这些沉积物中的盐粒及其他化学粉尘在风的作用下发生扩散。这些盐尘加重了周边区域的土壤盐碱化，从而导致土壤肥力下降、植被退化等一系列问题[23-24]。另外，在锡尔河中游和下游的植被变化情况也显著不同(图4)，其原因可能是锡尔河为中亚跨界河流，恰尔达拉水库位于哈萨克斯坦和乌兹别克斯坦的交界处，自苏联1991年解体以来，哈萨克斯坦从跨界河流中的引水能力有所下降[25]，因此锡尔河下游的植被显著降低可归因于其灌溉水资源的受限。

图4 1991—2015年中亚地区NDVI变化趋势Fig.4 Change trend of NDVI in Central Asia during the period from 1991 to 2015

3.3 NDVI空间分异与变化趋势驱动因子分析

3.3.1 不同驱动因子对NDVI空间分异的驱动分析

利用地理探测器模型分别就NDVI空间分异特征和变化趋势进行归因，因子探测结果如表4所示，不同因子按照对NDVI空间分异特征的解释能力排序：土壤湿度>湿日频率>降水量>土壤类型>潜在蒸散发>最高气温>平均气温>最低气温>地形区>云量>高程>日较差。这也说明中亚地区地处温带面积最大的干旱区，水分因素主导该地区植被的空间分布状况，其次是土壤类型和潜在蒸散发以及气温。风险探测分析和交互作用分析结果如图5所示，NDVI随着降水频率和土壤湿度的增加而增加，然而NDVI对不同等级降水量的响应呈现出先升高后降低的趋势。这可能因为在中亚东南部帕米尔高原的高海拔区域，降水量较多，土壤水分含量较高，植被的生长将不再受水分因素的限制，反而气温开始主导植被的生长状况。NDVI随着潜在蒸散发的增加而降低。随着云量增加，NDVI也呈逐渐上升的态势。NDVI随着高程先增加后降低，高程在1 500～2 400 m之间(分区5)时，达到峰值，当高程大于2 400 m后，NDVI开始急剧下降。植被生长发育的过程往往不是影响因子单独起作用，而是多种影响因子协同交互作用。从图5(b)可以看出，交互作用最强的因子组合为高程和土壤湿度，它们双因子增强交互作用的q值为0.670。这说明高程的变化显著增加了土壤湿度作为自变量对NDVI的解释。同时降水量与大多数因子结合都能产生较高的q值。高程和潜在蒸散发组合体现了较强的非线性增强交互作用(潜在蒸散发∩高程=0.647>潜在蒸散发(0.464)+高程(0.155))。除此之外，高程与多数因子的非线性增强交互作用都很强。

表4 因子探测器结果Tab.4 Results of factor detector

①地形区因子用于NDVI空间分异的驱动分析；土地覆被类型用于NDVI变化趋势的驱动分析。

(a)NDVI空间分异风险探测器结果 (b)NDVI空间分异交互探测器结果

(c)NDVI变化趋势风险探测器结果 (d)NDVI变化趋势交互探测器结果

图5 风险探测器结果和交互探测器结果

Fig.5Resultsofriskdetectorandinteractiondetector

3.3.2 不同驱动因子对NDVI变化趋势的驱动分析

从NDVI变化趋势的因子探测结果(表4)可以看出，不同驱动因子的q值并不高，但是依然可以反映出不同驱动因子对NDVI变化趋势的影响差异。如图5(c)所示潜在蒸散发增长最慢的区域(分区1)，NDVI增加的最快，这说明区域潜在蒸散发的增长对NDVI的增长有一定的抑制作用。降水量增加最快的区域(分区10)，NDVI增长的也最快，但是降水量变化较小和降水量下降的区域基本上NDVI变化不大。在不同高程范围内，NDVI增长较快的区域高程范围在300～1 000 m之间(分区4)，主要集中在一些低山、丘陵地区。NDVI对不同类型气温的响应基本一致，增温慢的区域NDVI增长迅速，这一点在最高温度变化方面最为明显。云量下降最快的区域，NDVI降低趋势远超其他区域，这也说明在中亚干旱区，云量一定程度上代表着区域的降水潜力。不同土地覆被类型对NDVI的变化也存在很大的空间分异特征，旱地(分区1)的NDVI下降最明显，草地(分区5)的NDVI上升最快。在土壤类型方面，黑土(分区10)和黑钙土(分区4)上生长的植被出现了明显的减少趋势，如图2(d)所示这2种土壤都是肥力较高的土壤且广泛分布于中亚北部地区，所以这些地区的土地很早就被开发为旱地。但是由于苏联解体使得大量农田弃耕，这些区域的NDVI下降趋势也较大。

在NDVI变化趋势影响因子的交互作用方面，80%以上影响因子组合方式表现出非线性增强作用，其余因子也表现出了双因子增强作用(图5(d))。交互作用解释力排在第一位的是云量和降水量组合(q值为0.208 0)，其次是潜在蒸散发和高程组合(q值为0.196 1)以及降水量和平均气温组合(q值为0.193 0)。除此之外，还可以看出潜在蒸散发和其他影响因子的交互作用表现明显，而湿日频率和其他影响因子的交互作用较弱。

4 结论

本文利用趋势分析法和地理探测器模型研究中亚地区NDVI时空分布与变化特征及其影响因子的驱动作用，并研究了双因子之间交互作用对NDVI时空分布及变化特征的影响，得到以下主要结论：

1)1991—2015年间中亚地区植被变化整体上处于平稳波动状态。从局部地区看，哈萨克丘陵中低海拔地区NDVI表现出明显的上升趋势，咸海西南部NDVI呈极显著下降，可能归因于咸海盆地盐尘的近距离扩散；同时水域NDVI的上升表明中亚湖泊整体水环境质量在下降。

2)人类活动对中亚地区植被覆盖的影响比较明显。由于苏联解体导致的弃耕后复垦现象，哈萨克斯坦北部旱地区域的NDVI下降趋势虽然很大，但是趋势并不能通过显著性检验；另外锡尔河中游和下游地区NDVI变化的显著差异也体现了中亚各个国家间围绕水资源开发利用问题上的一些矛盾。

3)根据地理探测器研究表明，水分因素主导着中亚地区植被覆盖的时空分布格局；在不同增温背景下，NDVI的变化趋势显著不同；另外不同高程范围、不同土地覆被类型、不同土壤类型下NDVI的时空变化特征也有显著差异。

4)因子间交互作用均能增强其对植被空间分异和变化趋势的解释力；潜在蒸散发和高程的非线性增强交互作用对NDVI空间分异的影响最为显著；云量和降水量的非线性增强交互作用对NDVI的变化趋势影响最为显著；另外，潜在蒸散发和高程作为影响因素可以较大地增强其他因子对NDVI空间分异以及变化趋势的解释能力。

然而，基于AVHRR卫星影像的NDVI数据虽然具有较长的时间序列和较高的时间分辨率，但空间分辨率较低。因此，下一步希望采用具有更高空间分辨率的MODIS数据，开展小尺度范围内的NDVI影响因子研究。另外，关于潜在蒸散发是如何增强其他因子对NDVI空间分异和变化趋势的解释能力，还需进行更加深入的讨论。