林 帅，刘会欣，张 强，刘 纯，樊晶晶，周 斌，刘 阳，黎云云

（1.河北工程大学，河北 邯郸 056000；2.昆明市水利水电勘测设计研究院有限责任公司，云南 昆明 650224；3.水利部海委漳卫南运河汤阴河务局，河南 汤阴 456150；4.绵阳师范学院，四川 绵阳 621000）

植物作为生态系统中的生产者，是自然界中联系分解者与消费者之间的重要纽带，研究植被变化，认识到生态系统的变化，制定相应生态管理措施对维护生态平衡以及维持经济的可持续发展有着重要意义。

归一化植被指数（NDVI）是表现植被变化的一个直观且有效的指标。在早期的植被研究中，为了更加直观地研究某一区域的植被变化，最早于1974 年由Rose 等提出了归一化植被指数这一概念，且现如今归一化植被指数已被广泛的应用于现代科学研究的各个方面。

Zhou 等为研究北半球不同大洲之间的植被变化，选取1981 年—1999 年北半球的NDVI 数据进行了植被变化研究[1]；WuD 等研究了植被与气象要素的滞后性，结果发现植被与气候之间存在着相应的滞后性[2]；在世界植被研究的过程中，国内的研究者也对中国大部分区域植被变化进行了相应的研究。李军瑶等在研究东北地区黑龙江多年冻土区植被NDVI在时间和空间上的尺度变化与地表温度的响应关系的过程中指出，黑龙江多年冻土区生长季GIMMS NDVI 和MODIS NDVI 对植被具有很好的指示作用[3]。基于大量的水文站数据与1982年—2015年的GIMMS NDVI数据，张勃等对中国15 个省极端天气对植被覆盖的响应进行了分析[4]；对于西北内陆陈文静等利用2001 年—2017 年的MODIS 数据，对NDVI 的变化趋势和其与气候因素相应的滞后性进行了分析[5]。变异诊断作为现在数据序列研究的一个重要环节，为各类数据趋势的研究提供了可靠的分界点。在变异研究过程中，刘琦等采用M-K 法和滑动T 法对穆棱河流域径流进行了变异诊断，确定了穆棱河流域的变异时间[6]。孙伯才采用改进的有序聚类法对北京朝阳区降雨的变异点进行了确定[7]。

为研究张家口市进入21世纪以来植被的变化情况，利用张家口市2000年—2019年的NDVI数据，对21世纪以来张家口市的植被变化进行时间与空间上的研究。

1 研究区域概况

张家口市位于河北省西北部地区，分为坝上和坝下两部分，坝上位于内蒙古高原，包括康保、沽源、尚义和张北县，平均海拔在1 400～1 600 m；坝下包括市区、宣化县、怀来县、万全县、怀安县、崇礼县、赤城县、涿鹿县、阳原县和蔚县，平均海拔在500～900 m[8]。张家口市四季分明，总面积为3.68×104km2，东经跨过113°50′～116°30′，北纬经过39°30′～42°10′，为亚热带大陆性季风气候。张家口市为缺水城市，全市年均降水量为400 mm左右，且降水主要集中在5月—8月，平均温度较低为7.7℃。张家口市地处重要位置，南邻北京市，向东靠河北省承德市，向西及西南地区与山西省接壤，向北及西北地区则与内蒙古自治区交界，张家口市流域位置图，见图1。

图1 张家口市流域位置图Fig.1 The location map of Zhangjiakou watershed

2 研究方法

2.1 滑动T法

在进行滑动T检验时需要将待检验的一组样本分为两段进行滑动检验，检验过程中通过判断两组数据的差异来诊断变异。具体公式如下：

在计算和判断过程中首先选择一个显著性水平α，在计算出来的结果中，如果|t|＜a，则说明变化不明显，如果|t|≥a，则说明变化显著且序列在ta时刻发生变异[9-10]。

2.2 改进的有序聚类法

改进的有序聚类法在诊断过程中需要对要评估的因素进行权重设置，并对其进行相应的有序归类，在计算过程中首先进行权重计算，其权重函数的计算方程为：

式中：f n(xj)——各分类因素的权函数变量；xj——评估指标；λn——分类值；λn-1——各类别的中心权重值。

在计算过程中，有序聚类法需要得到中心聚类值，因此需要对变量中心权重值进行标准化处理，得到标准化的各指标中心聚类值，其计算方程如下：

该方法需要结合评价指标及影响因素的种类进行综合性分析，其具体分析如下：

最后需对各影响指标之间进行定量分离,其分离度计算方程如下:

式中：dk——各影响指标的综合分离度，其他变量含义同上述方程中的变量含义[11-12]。

2.3 植被覆盖划分

为了更好地对植被进行空间变化分析，研究将NDVI 值进行了等级划分，根据不同NDVI 指数的范围变化将研究区的植被分布划分了4 个等级，分别为[13]：无植被区（NDVI≤0.2）、低植被区（0.2≤NDVI＜0.4）、高植被区（0.4≤NDVI＜0.6）和密集植被区（NDVI＞0.6），并通过Arcgis 软件来将栅格数据进行重分类，对植被的空间分布情况进行了读取分析[14]。

2.4 差值法

差值法可用来衡量一个地区不同年份的变动情况[14]，其基本表达式为：

式中：NDVIi——第i年NDVI 数值，NDVIj——地j年的NDVI数值，NDVIi-j——第i年与第j年之间的差值。

3 结果分析

3.1 张家口市植被的时空变化分析

3.1.1 年均植被变化分析

首先采用Arcgis 软件对张家口市的年均NDVI 变化进行分析，将2000年—2019年20张年值的栅格图像进行合成，计算20 年间张家口市NDVI 的整体空间变化，结果见图3，从空间上可以看出张家市的植被覆盖在坝上和坝下有着明显的分界，坝上的植被覆盖率偏低，坝下的植被覆盖率偏高，即由东南到西北呈现出逐渐降低的趋势。从时间上看，张家口市的NDVI 变化整体上呈增加趋势，多年均值为0.38，但增幅较弱，每年的增长率为0.001，且在2015 年前后存在较大变幅，见图3。

图2 2000年—2019年张家口市多年平均NDVI空间分布图Fig.2 The spatial distribution of annual average NDVI in Zhangjiakou from 2000 to2019

图3 年均NDVI一元线性变化曲线Fig.3 The annual average NDVI linear variation curve

3.1.2 生长季植被变化分析

同年均植被变化相近，经过空间计算，计算出生长季20年的均值空间分布图，见图4，对图4 进行空间分析，发现张家口市生长季植被覆盖在坝上与坝下亦存在着明显的界限，表现为坝上植被覆盖率偏低，坝下植被覆盖率偏高，植被覆盖率东南到西北逐渐降低。从时间变化上看，张家口生长季NDVI 亦呈现出了整体增加的趋势，其多年均值为0.52，增长率相较于年均值较高为0.002 4，且在2015 年也出现在了明显的波动，见图5。

图4 2000-2019年张家口市生长季平均NDVI空间分布图Fig.4 The spatial distribution of average NDVI in growing season in Zhangjiakou from 2000 to 2019

图5 生长季年均NDVI的一元线性变化曲线Fig.5 The annual average NDVI linear variation curve in growing season

综合分析，张家口市植被在变化过程中，以2015 年为节点存在着明显的变化，因此需要利用变异诊断法判断张家口市在2015年左右是否存在变异，且进一步分析变异原因。

3.2 植被变异分析

引入滑动T法和改进的有序聚类法对张家口市的植被进行变异诊断。

在采用滑动T 法进行诊断的过程中选用3 步和4 步对张家口年均植被变化和生长季植被变化进行滑动检验，图6为年均NDVI 和生长季NDVI 的滑动T 法变异诊断图，在两种步长的计算结果中，年均NDVI的变异点为2009年、2014年和2015年，生长季NDVI则出现在2009年、2015年和2016年，生长季与年均NDVI在2015年均出现了变异点。使用改进的有序聚类法对变异进行诊断，结果见图7，生长季和年均值的植被变异点出现在了2015年和2016年，综合来看，变异点2015年出现的频率最高，所以2015年为张家口市植被变化的变异点。

图6 滑动T变异诊断Fig.6 The diagnosis of sliding T variation

图7 改进的有序聚类变异诊断Fig.7 The improved ordered clustering variation diagnosis

3.3 变异前后植被变化分析

以2015 年为分界点，将NDVI 数据分为变异点前2000年—2014与变异点后2015年—2019年两个阶段进行分析。

3.3.1 年均NDVI变化

时间上首先对变异前后年均NDVI变化趋势进行趋势分析，然后在空间上对张家口市NDVI 在变异点前后的植被覆盖进行分析，见图8，为时间上年均NDVI 的变化曲线图，以2015年为分界点，在时间上张家口植被在变异前呈现出了快速增长的趋势，变异后则呈现出减小的趋势，且变异后植被覆盖有着明显减小的趋势。变异前后植被覆盖柱状图，见图9，空间上，变异点前高植被区的面积为12 121.96 km2，变异点后密集植被区的面积为10 708.08 km2，高植被区的面积由39.94%减少到29.10%，低植被区的面积在变异后则有所增加，由23 923.68 km2增加到了25 833.60 km2，增加了5.09%，密集植被区则几乎退化消失，因此变异点后张家口市的年均植被覆盖呈现出明显减小的趋势。

图8 变异前后年均NDVI变化趋势图Fig.8 The annual NDVI trend chart in Before and after variation

图9 年均NDVI植被面积柱状图Fig.9 The column plot of annual NDVI vegetation area

3.3.2 生长季NDVI变化

时间上生长季NDVI的变化曲线，以2015年为分界点，与年均一元变化曲线趋势相近，张家口植被在变异前呈现出了快速增长的趋势，变异后呈现出减小的趋势，且在变异后植被覆盖呈现出较为明显减小的趋势。空间上，变异前后生长季的植被覆盖以高植被区和密集植被区为主，在变异前高植被区的面积为3 627.43 km2，要低于变异后的8 078.53 km2，分别占到了整个流域面积的9.86%和21.95%，但变异后的密集植被区则有所较少，变异前密集植被区的面积为33 040.25 km2，变异后密集植被区面积由89.78%减少到了77.23%，由此对比可见，虽然变异后生长季的植被退化相较于变异后年均植被的退化程度较小，但结合一元线性图综合分析，仍然存在下降趋势。

图10 变异前后生长季NDVI变化趋势图Fig.10 The growing season NDVI trend chart in Before and after variation

图11 生长季NDVI植被面积柱状图Fig.11 The column plot of growing season NDVI vegetation area

3.4 旅游业发展对植被变化的影响

随着现代科学文化的进步与经济的快速发展，第三次工业革命不仅有科技的快速发展与进步，还伴随着文化的发展与进步，进入21 世纪以来张家口市的社会经济快速发展，且由于坝上地区海拔高，年平均温度较低，张家口市的滑雪文化也开始快速发展起来。

3.4.1 崇礼滑雪产业

首先以崇礼为例，随着滑雪旅游产业的发展，崇礼2015年—2016 的旅游人数为218.5×104人，增长了30.80%，2016年—2017 年旅游人数为267.6×104，涨幅为22.50%，2017年—2018 年涨幅较小为，6%，2018 年—2019 涨幅最小，为2.59%[15]，旅游人数趋于稳定，且2017年—2019年植被覆盖率也开始了小幅度的回升，由此可见2015年冬奥会申请成功后张家口市的滑雪产业建设和旅游产业发展进入了一个新的发展阶段，从而导致了张家口市崇礼及周边地区的下垫面发生了较大的变化，使植被发生变异。截至2020 年，张家口市滑雪场已建成9 家，占地面积约17 200×104m2，共拥有雪道169条，162 km。但在2017年后张家口市的植被又有所上升，这与滑雪场建设完成后期滑雪场内的植被生长（前期树木的生长）与人工绿化有着密切关系。

3.4.2 张家口冰雪及旅游人数变化

另一方面，对张家口市整体的冰雪及旅游人数进行分析，见图12，张家口市的旅游人数由2012 年的2 118×104人增长到了2015 年的3 848×104人，但增长率在逐年下降，2015年后伴随着冬奥会申办成功，张家口市的冰雪和旅游产业又开始了快速的增长，且增长率在2015 年后为最高，在1年内由3 848×104人增长到了5 193×104人。并且张家口市预计在2025年建成滑雪场30家，滑雪道600余条。对归一化植被指数的一元线性曲线与柱状进行综合分析，可知植被覆盖在2015 年的突变与张家口市的冰雪及旅游产业的发展有着密切的关系。

图12 旅游人数变化柱状图Fig.12 The cylindrical chart of changes in tourist numbers

3.4.3 张家口冰雪及旅游生产总值变化

张家口市冰雪及旅游产业所产生的GDP 变化与增长与旅游人数变化相似，变化过程中首先由2012年的128×108元增长到了2015 年的301.67×108元增长到了，接着2015 年—2016 年的增长则呈跳跃式由301.67×108元519.24×108元，增长了72.12%。对归一化植被指数的一元变化曲线的变化趋势与柱状图进行综合分析，进一步明确了植被的变异与张家口市旅游经济的快速发展有着密切关系，见图13。

图13 旅游生产总值变化柱状图Fig.13 The column chart of tourism GDP change

3.5 NDVI植被动态变化特征

3.5.1 年均NDVI变化

将张家口市年均NDVI变化以2000 年、2004 年、2009 年、2014年和2019年每隔5年为时间节点利用差值法对NDVI的动态变化特征进行分析。将各时间节点年均NDVI 分成：0～0.2、0.2～0.4、0.4～0.6、0.6～0.8 和0.8～1.0 5 段，将各时间节点的差值分成：＜-0.4、-0.4～0.2、-0.2～0、0～0.2、0.2～0.4与＞0.4 6段，对年均NDVI的分段特征进行分析，统计发现2000年、2004年和2009 年NDVI 主要占比在0.2～0.4，2014 年NDVI 主要占比在0.4～0.6，2019 年则回到了2000 年的水平，且年均NDVI 值为0.2～0.4时占比较大。对于差值而言2004年与2000年的差值和2014 年与2009 年的差值表现为植被的正增长，占比主要集中在0～0.2，2009年与2004年和2019年与2014年的差值表现为负增长，占比主要集中在-0.2～0，2019 年与2000 年的差值则主要集中在0～0.2 与-0.2～0，且占比相近，说明到2019年张家口市的年均NDVI 恢复到20 年前的水平，且与年均NDVI 的时间变化曲线图有很好的契合性。相关数据内容，见表1。

表1 年均NDVI植被动态变化特征表Tab.1 The annual NDVI vegetation dynamic change characteristics table

3.5.2 生长季年均NDVI变化

与年均NDVI 相似，将生长季NDVI 按各时间节点分成：0～0.2、0.2～0.4、0.4～0.6、0.6～0.8和0.8～1.0五段，将各时间节点的差值分成：＜-0.4、-0.4～0.2、-0.2～0、0～0.2、0.2～0.4 和＞0.4 6个阶段，对生长季NDVI的分段特征进行分析，分析统计发现各个时间节点生长季NDVI 的主要占比均在0.6～0.8，2019 年生长季植被覆盖要略高于2000 年的植被覆盖。对于差值而言2004 年与2000 年的差值和2009 年与2004 年的差值表现为植被的正增长，占比主要集中在0～0.2，2009年与2004年和2019 年与2014 年的差值则表现为负增长，占比主要集中在-0.2～0，且2014 年与2019 年的负增长程度较强。2019 年与2000 年的差值则主要集中在0～0.2 与-0.2～0，且0～0.2 的占比较大，说明到2019 年张家口市的生长季NDVI恢复到20年前的水平且相对有所增长，与生长季NDVI 的一元线性曲线图亦有着很好的契合性。相关数据内容，见表2。

表2 生长季NDVI植被动态变化特征表Tab.2 The table of NDVI vegetation dynamic changes in growing season

4 结论

首先利用Arcgis软件对张家口市的栅格图像进行相应的统计处理，得到张家口市年均NDVI 数据和生长季NDVI 数据，然后对数据进行处理，了解21 世纪初张家口市20 年间植被的时空变化；又采用2 种变异诊断方法对张家口归一化植被指数序列的变异进行诊断，整理分析变异点前后植被的变化程度；在经济方面对植被的变异进行综合分析，最后则用差值法对张家口市均值NDVI 和生长季NDVI 植被的动态变化特征进行分析，得到结论如下。

（1）对张家口市均值NDVI和生长季的NDVI进行时空分析，结果发现时间上张家口市在21 世纪初的20 年间，均值NDVI和生长季的NDVI变化均呈现上升趋势，但由于2015年后的突然变小，导致整体序列的植被增长并不是很明显，年均NDVI 的增长率为0.001，生长季NDVI 的增长率为0.0024。空间上张家口市的植被覆盖由东南到西北呈现出逐渐减小的趋势，即坝下植被覆盖率要高于坝上植被覆盖率。

（2）采用滑动T 法和改进的有序聚类法对张家口市2000年—2019 年的归一化植被指数（NDVI）序列进行变异诊断，结果发现在变异过程中2015年出现变异点的频率较高，所以可以确定2015年为张家口市植被变化的变异点。

（3）以2015 年为分界点对年均NDVI 和生长季NDVI 变化进行分析，对植被覆盖率进行统计，发现年均NDVI变异前在时间上呈现出增加的趋势，空间上变异点前低植被区小于变异后的低植被区，变异前高植被区则要大于变异后的高植被区，变异后密集植被区几乎退化殆尽，所以在时空上看，年均NDVI在变异后呈现出了明显减小的趋势；对生长季NDVI的分析发现，时间上与年均NDVI变化相似，同样是变异点前植被覆盖呈现出增加的趋势，变异后呈现出减小的趋势，空间上生长季的植被分布主要为高植被区和密集植被区，且变异点前的高植被区要低于变异后的高植被区，但变异前密集植被区要高于变异后密集植被区，总体来看生长季的植被覆盖亦有所减小。

（4）以崇礼区旅游业2014 年—2019 年的变化为例，在社会经济方面对植被变异进行分析，发现2015年来随着旅游产业的扩展，客流量的增加，导致了植被覆盖变异的发生，还得知了2017 年后植被覆盖的回升与景区植被的生长和人工绿化有着一定的关系。对张家口市的冰雪及旅游人数变化与植被的一元线性曲线进行综合分析，发现张家口植被的突变与冰雪及旅游产业有着密切关系。在经济增长方面，将张家口市冰雪及旅游产业的GDP与植被变化进行综合分析，更近一步确认了张家口植被变异与冰雪及旅游产业之间的密切关系。

（5）最后利用差值法对年均NDVI和生长季NDVI的动态变化特征进行分析，发现年均NDVI 只有2014 年的主要占比在0.4～0.6，其余各时间节点均在0.2～0.4，2004年与2000年之间的差值和2014 年与2009 年之间的差值表现为被覆盖有所增加，2009 年与2004 年和2019 年与2014 年之间的差值表现为植被覆盖有所减少，2019 年的植被覆盖则回到了2000 年的水平；生长季各时间节点NDVI主要占比则均处于0.6～0.8。增减趋势与年均NDVI 相 同，2004 年 与2000 年 和2014 年与2019年之间的差值表现为植被覆盖有增长的趋势，2009年与2004年和2019年与2014年的差值表现为植被覆盖有减少的趋势。且2019年的植被覆盖程度较2000年有轻微的上升。