基于地形位置指数的赤水河流域植被时空变化研究

马伟波， 赵立君， 田佳榕,2， 王 楠， 高媛赟， 李海东*

1.生态环境部南京环境科学研究所， 江苏 南京 210042

2.南京林业大学南方现代林业协同创新中心， 江苏 南京 210037

地形和灾害评价是资源环境承载力和国土空间开发适宜性评价的重要内容，也是村镇建设生态安全和产业规划的约束条件之一. 在山区地形条件下，地形因子对灾害危险性评估[1]、聚落分布优化[2]、“三生空间”协调发展[3]、生态功能提升等村镇建设问题将产生重要影响. 研究山区植被时空变化及地形分异特征对于村镇建设具有重要指导意义. 地形位置指数(Topographic Position Index, TPI)是一种测量地形位置和地形分类的方法[4]，坡度位置方法在其基础上综合了坡度信息，通过地形地貌划分消除了绝对高程特征在区域尺度范围的影响，这两种分析方法可在微地形、中尺度及景观地貌等多个尺度表达地形特征，被应用于热带森林生物量估算[5]、森林生态系统中植被土壤属性和地形对微生物的影响分析[6]、野生动物保护活动范围探索[7]、林场管理[8]、景观异质性分析[9]以及科尔沁沙地植被生长分析[10]等研究.

赤水河流域地处云贵高原向四川盆地倾斜的斜坡面，是大娄山区水源涵养与生物多样性保护重要区(简称“大娄山区”)的重要组成部分. 赤水河流域地形以山地丘陵为主，长期以来由于过度垦殖、滥砍乱伐等，上游地区植被破坏，水土流失严重，生态系统退化；中下游地区小煤窑、酒作坊和城镇扩张对赤水河水环境威胁较大[11]. 地形影响着区域水、热、养分等的再分配过程[12]，进而影响植被生长和分布格局，对区域生态系统服务价值有一定作用[13-14]. 同时，地形因子约束人类生活与生产活动，影响区域土地利用的方式和程度[15]，是村镇建设可持续发展的重要影响因素. 有学者从地形位指数等指标出发，分析得到不同地形梯度下土地利用变化与植被生长的关系，揭示了地形因素对于区域生态功能和生产建设的影响规律[16-18].

为阐释影响村镇建设的生态约束条件，笔者基于SPOT_VGT NDVI产品和高程等数据，通过Sen′s slope、M-K检验、TPI和坡度位置方法，分析了赤水河流域植被生长季NDVI时空变化及地形分异特征，以期为赤水河流域村镇建设生态安全格局优化提供支撑.

1 材料与方法

1.1 研究区概况

赤水河流域位于川、滇、黔三省接壤地带(见图1)，地理位置为26°49′N～28°54′N、104°09′E～107°10′E，干流总长度433.32 km；流域总面积1.83×104km2，占大娄山区面积(32 854 km2)的55.7%；高程在190～2 219 m之间，落差达 2 029 m. 赤水河发源于云南省镇雄县，流经四川省、贵州省，最终汇入长江. 赤水河流域上游属于云贵高原，中上游为喀斯特地区，下游为丹霞地貌地区；流域降水主要集中在6—9月，年降水量为700～1 100 mm；流域土地利用类型以林地为主，耕地和草地次之，2018年林地总面积为 10 556.77 km2. 赤水河主要支流包括上游的二道河，中游的桐梓河、古蔺河，下游的大同河和习水河.

图1 赤水河流域地形特征和大娄山区位置

赤水河流域分布有4个国家级自然保护区，分别为长江上游珍稀特有鱼类国家级自然保护区、赤水桫椤国家级自然保护区、四川画稿溪国家级自然保护区和遵义习水中亚热带常绿阔叶林国家级自然保护区，保护区总面积达 1 450 km2，占赤水河流域总面积的7.9%，流域水源涵养与生物多样性维护功能突出[19]. 赤水河上游地区分布有国家级贫困县和国家级开发扶贫区；中下游地区经济发展水平明显高于上游地区，仁怀市地区生产总值最高，2019年达到1.30×1011元；流域城镇化水平较低，产业发展以农业生产为主，酿造业尤为发达.

1.2 数据来源

采用时间序列SPOT_VGT NDVI卫星遥感数据分析赤水河流域植被生长变化情况，该数据通过中国科学院资源环境科学与数据中心(http:www.resdc.cn)下载. 该数据集时间分辨率为30 d，空间分辨率为1 km，时间序列为1998—2018年. 为了消除云、大气和月内物候的影响，采用最大值合成法得到每月NDVI最大值；同时，为减少异常噪声，对时间序列NDVI数据进行逐像元Savitzky-Golay卷积平滑.

高程数据来自SRTM(Shuttle Radar Topographic Mission)90 m的数字高程模型(Digital Elevation Data，DEM)，通过SRTM官网(http:srtm.csi.cgiar.orgsrtmdata)下载.

1.3 研究方法

1.3.1年际变化趋势分析

选择Sen′s slope方法分析赤水河流域植被生长季NDVI年际序列变化趋势，该研究中生长季为5—10月. 该方法是一种时间序列非参数检验检测方法[20]，可以估算单位时间的变化幅度(β). 如果一个时间序列xt=(x1,x2,…,xn)存在线性趋势，则当β>0时，序列呈上升趋势；当β=0时，序列趋势不明显；β<0时，序列呈下降趋势.β的计算公式：

(1)

式中，Median为取中值函数，xj、xi分别表示第j年和第i年生长季的NDVI平均值，且j>i.

1.3.2Mann-Kendall非参数检验

采用Mann-Kendall非参数检验(M-K检验)方法对赤水河流域植被生长季NDVI年际序列变化做显著性检验. 通过M-K检验的Z检验值判断赤水河流域植被变化趋势的显著性[21-22]. 对于时间序列xt={x1,x2,…,xn}，首先计算统计值S：

(2)

(3)

Mann-Kendall检验的Z检验值计算公式为

(4)

式中，Z>0表示植被生长呈增长趋势，Z<0表示植被生长呈下降趋势. |Z|越大，说明该时间序列的变化趋势越显著；|Z|≥1.96和|Z|≥2.58时，分别代表通过95%和99%的显著性检验.

1.3.3地形位置指数

(5)

(6)

式中：n表示R半径内DEM栅格的数量；zi表示第i个栅格处的高程，m. TPI越大，表示该位置处于较高区域，即可能处于山脊；TPI越低，表示该位置处于整体较低位置，即可能处于山谷. TPI的取值范围不仅取决于高程的差异，还取决于R，较大的R可使TPI表达大尺度地貌单元，而较小的R可使TPI表达较小尺度地形特征. 根据相关研究[9-10,23]，结合赤水河流域山区地形的特点，通过对3～1 200个像元单位的测试，最终选择R为20个DEM像元尺度，即半径为 1 800 m的圆形卷积窗口.

1.3.4坡度位置

地貌变化源于坡面变化[23]. Weiss[4]将TPI与坡度信息进行组合，将地貌划分为6类(见表1)，分别是山脊、上坡、中坡、平坡、下坡和山谷. 坡度位置综合了坡度(slope)和TPI所蕴含的地形特征信息[24].

表1 基于TPI的坡度位置划分

2 结果与分析

2.1 NDVI的时空变化特征

由图2(a)可知，1998—2018年赤水河流域植被生长季平均NDVI为0.64～0.79，多年平均NDVI为0.70，2017年最高(0.79)，1998年最低(0.64). 赤水河流域内，1998—2018年植被生长季平均NDVI呈缓慢上升趋势，斜率为 0.004 7；其中，1998—2009年呈上升趋势，2010—2012年出现短暂连续下降，在2012年和2016年发生较大幅度的下降，但在2017年迅速上升并达到最高值. 由图2(b)可见，赤水河流域植被生长季平均NDVI在空间上整体处于较高水平，NDVI≥0.55的区域面积为1.10×104km2，占整个流域的60.10%；其中NDVI>0.60的集中连片区域主要分布在古蔺县北部、赤水市大部和习水县西北部，占赤水河流域总面积的8.42%. NDVI<0.40的区域面积为19.47 km2，主要分布在仁怀市的中枢街道、茅台镇，习水县的东皇镇，桐梓县的娄关山镇以及赤水市的市中街道. NDVI在0.40～0.55区间的区域面积为 7 220.94 km2，占流域总面积的39.46%，主要分布在赤水河中上游地区以及下游的合江县.

图2 1998—2018年赤水河流域植被NDVI年际变化及其空间分布

就1998—2018年赤水河流域植被生长趋势而言〔见图3(a)〕，Sen′s slope>0的面积为1.81×104km2，占流域总面积的99.05%，Sen′s slope<0的区域面积为172.78 km2，占流域总面积的0.95%；可见，1998—2018年赤水河流域植被生长整体呈增强趋势. 在Sen′s slope的6个梯度区间中，-0.014～-0.001区间植被生长呈快速下降趋势的地区主要分布在仁怀市、习水县和赤水市；而Sen′s slope在0.009～0.015区间时，赤水河流域植被生长增强趋势最明显，主要集中分布在赤水河中上游、二道河以及下游的大同河干流地区. 结合图2(a)和图3(a)可知，与简单的线性拟合相比，Sen′s slope分析可在空间层面展示植被生长变化趋势，且具有更强的统计学意义.

由M-K检验结果〔见图3(b)〕可知，1998—2018年赤水河流域植被生长极显著增加的地区面积为1.71×104km2，占比为93.81%，分布在赤水河流域绝大部分地区；植被生长显著增加的地区面积为580.18 km2，占比为3.19%，零散分布在未通过显著性增加检验地区的外围；1998—2018年赤水河流域植被生长极显著减少的地区面积为34.71 km2，占比为0.19%，主要分布在赤水市中山街道、习水县东皇镇以及仁怀市茅台镇、盐津街道和中枢街道；植被生长显著下降的地区面积为36.69 km2，占比为0.20%，集中分布在未通过显著性降低检验地区的外围. 总体而言，1998—2018年赤水河流域植被生长季NDVI时间序列通过极显著性检验的地区面积占比超过94%，通过显著性检验(含极显著性)的面积占比超过97.28%，表明植被生长时间序列趋势显著.

图3 1998—2018年赤水河流域植被Sen′s slope 和趋势显著性分级

2.2 NDVI的地形分异特征

由图4(a)可知，TPI负值较小的区域表征信息多为河谷特征，而TPI正值较大的地区表征信息为山脊特征. 其中，TPI在-39.4～34.3区间的面积最大，为 6 221.63 km2，占赤水河流域总面积的34.05%，表明TPI起伏较小的区域占据流域主体；而TPI在-407.8～-127.1和118.5～486.9区间的面积占比较小，分别为8.24%和9.46%. 由图4(b)可知，6种坡度位置类别中，以中坡面积(7 792.02 km2)最大，占流域总面积的42.64%，表明TPI数值较小且坡度大于5°的区域是赤水河流域地形主体；其次为上坡和下坡，面积分别为 3 063.06 和 2 928.32 km2，占比分别为16.76%和16.03%；再次为山谷和山脊，面积分别为 2 006.31 和 1 932.32 km2；平坡的面积最小，为551.12 km2，占比为3.02%. 结合表1所示提取规则可知，赤水河流域内平地分布极少.

图4 赤水河流域TPI和坡度位置空间分布

就1998—2018年赤水河流域植被生长季NDVI的TPI梯度特征而言(见表2)，通过自然断点法划分为5个区间，TPI在118.5～486.9区间的平均NDVI最高(0.762)，而该区间的平均Sen′s slope也达到最大(0.007 3)，即TPI在该区间内时，植被生长状况整体较好，且增长趋势最强. TPI在-127.1～-39.4区间时平均NDVI最低(0.715)，同时该区间平均 Sen′s slope 最低，为 0.007 1.

表2 基于TPI的赤水河流域植被生长季NDVI地形分异特征

由表3可知，赤水河流域植被在山脊的平均NDVI最高，为0.747，且山脊平均 Sen′s slope最高，为 0.007 2；而上坡和中坡的平均NDVI次之(分别为0.737和0.725)，对应的平均 Sen′s slope均为0.007 1；1998—2018年，赤水河流域山脊处植被生长状况整体最好，上坡和中坡次之；同时，山脊处植被生长增加趋势最强，上坡和中坡次之. 山谷平均NDVI最低，为0.709；平坡处平均 Sen′s slope最低，为0.006 6. 赤水河流域植被在山脊、上坡和中坡处植被NDVI较高，同时该地区植被生长增加趋势高于平坡、下坡和山谷地区.

表3 基于坡度位置的赤水河流域植被生长季NDVI地形分异特征

3 讨论

赤水河流域植被生长变化趋势与文献[29]所得结果相吻合，即赤水河流域植被生长整体呈增加趋势，且在赤水河中上游地区的增加趋势较为显著. 蔡宏等[18]对赤水河流域不同植被覆盖度与地形因子的相关强度研究指出，赤水河流域在高程超过 2 000 m地区植被的覆盖度较高；而笔者在其基础上增强了时间序列和显著性检验分析. 同时，就整个流域尺度而言，高程为 2 000 m的地区主要集中分布在赤水河上游地区，对于其他地区的植被覆盖分析存在盲区；而笔者研究的优势在于通过TPI区间划分和坡度位置分析，消除了绝对高程的影响，并且从坡度位置角度分析流域植被生长情况，得到的分析结果在整个流域尺度具有一定的横向对比功能，即整个流域范围的山脊和山谷可横向比较分析. 相较于地形位指数[16-17]和地形起伏度方法[14]，TPI和坡度位置分析方法省去了对数换算，简化了地形分析复杂度.

与在低山和草原地区不同[10]，在山区地形中山脊与山谷处植被光热条件有较大差异，其生长情况和增长趋势也体现出较大差异. 但受城镇扩张影响，2000—2018年赤水河流域建设用地从24 km2增至569 km2，面积增加了约23倍，主要是中坡、下坡和山谷的耕地大面积转变为建设用地，而山谷建设用地从6.32 km2增至116.34 km2，占建设用地面积的20.45%. 实地调查发现，赤水河流域下游河谷地区粗放型小企业、小作坊无序发展问题突出，在仁怀市、习水县、桐梓县和赤水市河谷地区植被覆盖度下降较明显，城镇建设扩张和国土空间开发对区域资源环境承载力能力形成较大压力，增加了区域生态风险[24-25].

赤水河流域耕地以坡耕地为主，其中山脊处的耕地面积达390 km2，而上坡和下坡处的耕地面积占耕地总面积的32.64%. 坡耕地本身土地肥力低，长期在坡耕地活动不仅加剧喀斯特地区的石漠化[26]，对该地区农业生产造成负面影响，而且山脊处地势陡峭，极易发生水土流失甚至造成山洪灾害[27-29]. 因此，可结合地形分异特征对赤水河流域耕地资源进行优化调整，通过对其加以引导和提升来加强产业生态化和生态产业化村镇建设，改变生产经营方式，发展生态农业、生态旅游及相关产业，优化村镇建设承载力条件[30]，提升耕地质量，维护流域生态安全.

相较于草地和耕地，林地对流域生态服务价值的贡献更大，而赤水河流域林业类型中针叶林和阔叶林面积相对较小，且主要分布在古蔺县北部、赤水市大部和习水县西北部. 虽然近年来赤水河流域林地面积整体上从1.01×104km2增至1.09×104km2，但在中上游地区人工造林存在林分结构单一等问题. 后期应注重在赤水河流域中上游的中坡、下坡和山谷地区垂直生境多样性的规划和建设，以促进流域生物多样性、提升流域生态服务功能；同时，通过优化植被生境和土地利用空间格局，流域水资源供给和调节能力可进一步得到提升.

植被生长变化的影响因素极其复杂，笔者仅对赤水河流域植被生长季NDVI取平均值，忽视了其他时间尺度下的植被生长特征. 后期可结合全年NDVI数据设计更加具有代表性的植被指数，同时利用土地利用格局变化和产业结构变化揭示地形因素对于村镇建设生态安全约束的影响因素. 另外，对于不同植被类型的物候特征、地带性植物、流域水源涵养能力、生境质量的地形分异特征后期可做进一步探索.

4 结论

a) 赤水河流域植被生长整体较好且呈增强趋势. 1998—2018年植被生长季平均NDVI呈缓慢上升趋势，斜率为0.004 7；NDVI >0.60的集中连片区域主要分布在古蔺县北部、赤水市大部和习水县西北部，占赤水河流域总面积的8.42%；植被生长增强趋势最明显区域对应的Sen′s slope区间为0.009～0.015，主要集中分布在赤水河中上游沿赤水河干流、二道河干流地区以及下游的大同河干流地区.

b) 赤水河流域植被分布在118.5～486.9的TPI区间或分布在山脊处时其整体生长较好，且生长增强趋势最明显.

c) 建议对山谷坡度位置的村镇生态建设加强约束和引导，以维护村镇建设生态安全.