利用高分辨率森林覆盖影像实现高山林线的自动提取

江鑫，何心悦，2，王大山，邹俊宇，曾振中

1.南方科技大学环境科学与工程学院,深圳 518055;2.利兹大学地球与环境学院,英国利兹LS29JT

1 引 言

森林的生长和分布受环境和气候类型等不同因素的影响。林线位置的改变被认为是森林生态系统对气候变化最敏感的响应之一（Bharti 等，2012；杨秀云等，2020）。林线根据森林的高程可划分为低林线和高山林线两种类型。其中，低林线区主要位于森林与草原或农田的交界处，往往受到气候变化和人类活动（例如放牧、农业扩张）的干扰而发生改变（Zhang 等，2006）。高山林线周边植被受到人类活动的影响较小，该区域的森林变化在很大程度上归因于气候变化（王亚锋等，2020）。高山林线是标定全球气候变化的重要依据，是科学管理森林资源的基础，其中林线的位置、区域面积和改变模式等也是作为环境监测与建模最基本的信息（刘锬和康慕谊，2016）。因此，高山林线的变化是全球气候变化、生物多样性和森林生态系统等领域的研究热点（Theurillat和Guisan，2001；Gruber等，2009）。

尽管目前高山林线的定义大多基于树木的高度和森林覆盖率，不同研究人员对林线的定义却不尽相同。例如，高山林线可被定义为最高海拔地区树木高度低于2 m 的地带（Kullman，2001），或是郁闭林与高山矮曲林带之间的过渡带（Körner和Paulsen，2004），抑或由恶劣自然环境造成的包含高山矮曲林带的封闭森林（Tranquillini，1979）。本文高山林线定义为由树种线和郁闭林终止线等组成的林线交错带，包括密闭林地的边界以及此海拔以上的灌丛和草地等稀疏植被过渡带（王晓春等，2005）。其中，树种线是指郁闭林之上的高山低矮植被的边界。高山林线属于典型的生态交错过渡带，其位置依赖于生态系统的环境并且具有一定的动态性。因此，准确提取高山林线尚存在一定难度。

快速准确地提取高山林线的分布情况是研究植被动态对气候变化的响应关系的重要组成部分，即通过提取的林线信息分析植被变迁特征与环境条件的影响（邵佳怡等，2019；边瑞等，2020）。目前获取高山林线分布情况的方法主要包括植物群落实地测量法、地理统计法和遥感监测法等3类（Kullman，2001；Körner 和Paulsen，2004；Gehrig-Fasel等，2007）。其中，野外调查法通过实地观测树木形态特征和分布情况来获取高山林线位置。尽管能够准确地刻画林线的位置，但是高山地区存在不易到达、测绘难度大等问题，使得研究者难以获得大范围的林线分布情况，调查成本高且工作效率低。地理统计法是基于土地覆盖数据，采用滑动窗口方式推断出森林的内生和上移等情景。然而，地物类型复杂多变，且统计模型都存在由随机抽样等引起的“不确定性”问题，如何阐述和解决模型中的参数不确定性问题需要进一步分析研究。相比而言，遥感技术具有重访周期短、观测范围大和不受地理环境等条件限制的等特点，且不受天气和地理环境等条件的限制，已被广泛应用于森林资源监测和管理（史军，2002；荆雪娟，2019）。因此，具有识别树木能力的高时空分辨率遥感技术成为了提取大范围高山林线、监测林线动态的有效手段。高空间、时间分辨率遥感卫星的出现极大提高了林线调查的效率（Luo和Dai，2013；刘锬和康慕谊，2016）。

秦岭山脉是中国重要的南北地理分界线。秦岭高山林线位于高海拔针叶林、灌丛和草甸的生态交错带（马新萍，2015），植被垂直带谱分布明显，为本文提供了理想的试验区域。马里兰大学Hansen 等（2013）基于Google Earth Engine 云计算平台发布了自2000年以来20 a 的全球30 m 高分辨率森林覆盖数据。该产品包含了全球2000年的森林覆盖比例以及之后逐年的各格网森林增加/减少的信息，具有精度高、时序性强等特点，适用于秦岭地区的高山林线提取。本文首次基于此高分辨率森林覆盖数据，结合美国宇航局NASA（National Aeronautics and Space Administration） 提供的高程数据以及秦岭山脉分布数据，提出一种基于遥感的快速提取高山林线的算法。首先，本文从森林覆盖数据中提取树木分布地图，通过高程数据和秦岭山脉分布范围确定研究区的最高点；然后，采用八连通域搜索算法寻找森林与非森林的边界，以确定秦岭高山林线的范围；最后，本文结合高程数据系统分析了高山林线分布与海拔、坡度和坡向的关系，探索林线分布的规律，以期为山地生态恢复与建设提供依据。本文设计的山脉最高点搜索模块可以显著提升处理高分辨率森林覆盖数据的算法效率，并且该算法不受研究区范围大小的影响，可用于全球其他山脉高山林线的研究。

2 数 据

2.1 研究区

本文的研究区域为位于陕西省境内的秦岭高山地区，经纬度范围分别为33°39′11″N—34°19′36″N 和106°27′04″E—108°38′55″E，如图1所示。其中图1（c）是来源于2000年全球森林覆盖数据集的树冠密度图，该数值越大表示树木分布越密集，即出现高山林线的概率越小，本文将冠层密度在5%以内的区域视为林线候选区。研究区东西绵延约120 km，南北宽达75 km。研究区南北坡植被类型差异明显，南坡植被类型主要包括亚热带常绿阔叶林和落叶阔叶混交林，北坡植被则主要分布暖温带落叶阔叶林（邓晨晖等，2018）。此外，秦岭山地植被分布特征呈现明显的垂直结构。在海拔高度2400—3300 m 的区间内植被以针叶林为主，包括巴山冷杉、秦岭冷杉和太白落叶松等树种林。由此区间往上，植被逐渐过渡为灌丛、草甸等低矮植被（马新萍，2015）。

图1 研究区位置、树冠密度和地形概况Fig.1 Location，canopy density and topographic overview of study area

2.2 数据

Landsat 卫星提供了30 m 空间分辨率的遥感影像。最早的数据产品可追溯到20 世纪70年代初期，使得该系列卫星具有分析长时间序列地表覆盖变化的能力（陈军等，2017；王颖洁等，2015；Liu 等，2020；宫鹏，2021）。Hansen 等（2013）基于多种类型森林样本数据和Landsat时间序列光谱曲线，采用决策树的机器学习算法制作了2000年以来逐年全球森林变化数据产品。与MODIS 提供的MCD12Q1 土地覆盖数据相比，该套数据具有更高的空间分辨率，在植被和环境监测等方面具有一定的优越性。因此，本文基于2000年的Hansen 森林覆盖数据确定研究区森林与非森林分布区域。

目前，常用的数据高程模型包括SRTM（Shuttle Radar Topography Mission）、AW3D30（ALOS World 3D）和ASTER GDEM（Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model）3种（唐新明等，2021）。上述数据产品均可免费获取、具有较高的空间分辨率，已被广泛应用于遥感、测绘和军事等相关领域。然而，不同高程模型的适用场景存在细微差异。例如，不少研究表明在地势起伏较大的区域，SRTM 和AW3D30 的测量精度受地形影响较大，易出现高程异常。考虑到秦岭地区的地形呈现差异性变化，因此本文选择ASTER GDEM 高程模型确定林线区域所在的最高点。为提高林线搜索性能，本文从GMBA（Global Mountain Biodiversity Assessment）全球主要山脉分布数据集中提取秦岭最高峰的覆盖范围以确定林线搜索区间。另外，本文收集了GPS 地面站点观测数据（Dang 等，2015；Shi等，2020）和高分辨率的Google Earth 影像用于评估基于遥感林线提取算法的有效性。Google Earth 可视化平台不仅免费提供全球高质量的时序性遥感和地图数据（Landsat，CNES 和Maxar），而且支持导入点、线和面矢量数据和服务。基于该可视化平台，易验证林线提取算法的准确性。本文采用的数据集信息如表1所示。

表1 数据集描述Table 1 Dataset description

3 研究方法

本文提出的高山林线识别算法如图2所示。算法流程主要分为以下3个阶段：数据预处理、山脉最高点搜索和高山林线轮廓提取阶段。由于不同数据集的空间分辨率存在差异，首先需要对输入数据进行重采样预处理，然后基于山脉分布数据采用分块策略确定山的最高点，最后，利用连通域算法从最高点搜索林线最大范围。

图2 技术流程图Fig.2 Technical flowchart

3.1 数据预处理

山脉栅格数据分辨率低，为获取更准确的研究区位置，需要在同一空间分辨率下对山脉数据和森林覆盖数据进行运算。最邻近内插法（Nearest neighbor interpolation）运算速度快，且不会改变原始影像值，适用于土地利用分类、植被类型统计和森林覆盖等专题地图（Zheng 等，2017）。鉴于简单高效的原则，本文最终使用该插值方法将山脉分布数据重采样到30 m空间分辨率。

3.2 高山林线迭代提取算法

3.2.1 高点搜索算法

秦岭山脉众多，海拔较高的山峰主要有：华山（2154 m），终南山（2604 m），光头山（2887 m），鳌山（3475 m），太白山（3771 m）等。秦岭高山林线位于高海拔针叶林、灌丛和草甸的生态交错带，该区域植被覆盖度低，高程数据与树冠密度是寻找林线位置的关键信息。因此，本文第一步需确定研究区最高点坐标及其对应位置的植被覆盖度信息。预处理后的ASTER GDEM 和GFC 森林分布数据提供的高程范围均为10°×10°，行列为40000×40000，远大于图1（a）所示的秦岭山脉覆盖区。为提升最高点查找效率，本文结合高程数据和30 m 分辨率的秦岭山脉分布数据用于约束研究区高山林线的搜索范围。然后，结合秦岭山脉范围内最高点坐标和GFC 森林数据确定对应位置的树冠密度。

3.2.2 分块处理策略

秦岭山脉覆盖范围大，海拔2600 m 以上的地区占比约0.78%。由于研究区高山林线主要集中在高海拔地区，可能存在多个林线区域。为加快林线搜索效率，在第一轮搜索时以山的最高点为中心，边长为r的缓冲区作为高山林线的候选区域。经测试，行列为8000，边长r约为240 km 的方形区域能覆盖全球大部分的高海拔地区。紧接着第二轮以山脉的次高点（即未被搜索部分的最高点）再做缓冲区，以此类推，直到选出的高点位置植被类型为森林（森林覆盖率≥5%），则终止林线搜索。图3为林线搜索的示意图，第一轮搜索林线的结果如图3（a）所示，第二轮的林线迭代搜索如图3（b）所示。图3 中红色线条为提取的林线边界，绿线为分块搜索区。这种分块处理的策略适用于寻找全球不同地区的高山林线分布。

图3 高山林线迭代提取算法示意图Fig.3 Alpine treeline iterative algorithm

3.2.3 八连通域搜索

分块缓冲区包含树冠低密度区域和树冠高密度区域两部分。本文将树冠密度低于5%的区域作为林线候选区，高于该密度的区域作为背景区（FAO，2020）。随后，本文采用连通域检测算法从图像中分离背景信息，进而获取林线范围。具体地，本文先根据树冠密度值将分块缓冲区分为二类，即将森林分布数据以树冠密度5%为界进行二值化。经过该处理可突出图像整体和局部特征，简化高山林线的提取流程。再基于山顶点的经纬度坐标信息和区域增长算法对二值影像进行标记，从而获取高山林线的轮廓坐标信息。连通域检测算法示意图如图4所示。

图4 基于种子填充法的连通域搜索示意图Fig.4 Connectivity domain search with seed filling method

该算法已知图像区域内部某像素的值（山顶点坐标），向外拓展并确定该图像所有与该点相连通且等值的其他像素。基于该算法能快速填充多边形内部的特征值，可用于标识高山林线区域范围（Fiorio和Gustedt，1996；Wu等，2005）。在此，本文对遥感图像I中的像素点I(x，y)的邻域定义为

像素点的八邻域表示为

在对应的八邻域像素点集合中，若I(x，y) =I(x+i，y+j)，则代表点I(x，y)与点I(x+i，y+j)存在八连通的关系。此方法生成的连通域为二值图像中相连的前景像素或背景像素区域。其中，互相连通的前景像素为感兴趣的连通域。另外，前景区域可能存在一个或多个连通域，由于林线搜索算法确定了起始搜索点，因此，一轮迭代后只标记一个连通域，即林线区域。

4 结 果

4.1 树冠密度与高程的关系

秦岭高山植被覆盖度与海拔高度的关系如图5。其中深黄色为冠层覆盖度大于5%的统计结果，绿色为树冠覆盖度小于或等于5%的统计结果，该值越小则表明该地区树木分布越稀疏。相关文献表明（唐志尧等，1999；马新萍，2015），海拔2600—3200 m 的区间属于针叶林带，树冠覆盖度呈现下降趋势，海拔3300—3400 m 范围为太白红杉林和亚高山灌木的生态交错带（李娜，2019）。海拔高于3400 m 的地区，主要分布着高山灌丛、草甸等低矮植被类型，与相关文献描述一致（戴君虎等，2001）。因此，植被覆盖度与高程的关系可为寻找林线提供关键信息。

图5 冠层覆盖度与海拔高度的关系Fig.5 Relationship between canopy cover and altitude

4.2 基于遥感的高山林线提取与验证

4.2.1 秦岭高山林线空间制图

本文对高山林线的定义为各坡面上由密闭林地边界与此海拔以上稀疏林地组成的过渡带，即密闭林地与灌丛、草甸等稀疏林地的分界线。本文使用基于遥感的全球森林覆盖产品，采用连通域搜索算法检测稀疏林地，由此确定秦岭高山林线的分布位置（图6）。图6（a）为所有林线区域分布情况，蓝色、黄色和橙色分别为第1 轮迭代、第2轮迭代和第3轮迭代提取的林线范围，林线轮廓如图6（b）、图6（c）和图6（d）所示。五角星为每次迭代搜索的最高点位置，高程分别为3764 m，3532 m和3411 m，更详细的统计结果，例如林线区最高点坐标和连通域面积信息，如表2所示。

表2 林线区域面积和最高点信息Table 2 Treeline area and peak point information

图6 高山林线提取结果Fig.6 Alpine treeline extraction results

4.2.2 精度评价方法

为验证高山林线提取结果的准确性，本文基于高分辨率的Google Earth 影像（Landsat，CNES 和Maxar）、GPS 地面站点观测数据（Dang 等，2015；Shi 等，2020）和归一化植被指数（NDVI），采用目视解译方法对高山林线的位置进行评估。蓝线为密闭林地与稀疏林地的分界线，区域内月均气温低，林地稀少且以耐寒的高山灌丛和草甸植被为主（秦进 等，2017）。黄色标记位置为GPS 记录的林线点观测坐标（33°57′00″N，107°36′36″E）。本文根据坐标的有效位数取500 m作为该点的缓冲区以划定坐标误差范围（图7）。从图7 中可以看出，GPS 记录的坐标位置及其缓冲区均位于森林内部，明显偏离实际高山林线。

图7 GPS记录的站点坐标Fig.7 GPS recorded field coordinates

高山地区不易到达、测绘难度大，实地站点记录少。NDVI 记录植被信息，反映了不同区域的植被覆盖度，运用该指数可判定林线附近植被过渡变化。NDVI 数值越小，植被覆盖度越低。本文采用2000年研究区30 m 分辨率的NDVI 空间分布数据作为参考（图8），用于辅助评估提取的高山林线精确度。图中林线Z1、Z2和Z3区域内与林线区域外的NDVI 具有明显差异，内部植被覆盖度低于外部。而且，对比本文提取的高山林线和研究区Google Earth 影像可以发现，本文结果与Google Earth 影像中实际林线分布基本一致，进一步说明本文提出的林线搜索算法的优越性。

图8 2000年研究区最大NDVI空间分布Fig.8 Spatial distribution of maximum NDVI in 2000

4.3 高山林线分布特征

4.3.1 林线高度信息

研究区的高山林线主要分布在海拔2400—3800 m 范围，平均海拔为3061 m（图9）。3 处高山林线区均集中分布在顶峰区域。其中，林线Z1 处于海拔最高位置，垂直范围约为2389—3662 m，平均高度为3119 m；林线Z2 和林线Z3位于海拔相对较低位置，平均高度分别为2995 m 和3072 m（表3）。

表3 秦岭高山林线高程信息Table 3 The information on the elevation distribution of treeline in Qinling Mountains

图9 林线高程分布Fig.9 Treeline elevation

研究区林线的上边界均高于3400 m，这与前人对秦岭高山林线高程的研究结果一致（崔海亭，1983）。林线下边界高程接近1780 m，与前人描述的2400 m（戴君虎等，2001）高差达600 m。但本文结果显示绝大部分林线高程都高于2600 m，像元占比高达96.21%，与相关文献描述一致（马新萍，2015）。通过叠加Google Earth 影像可以看出，林线下边界主要位于陡峭山脊处，这意味着以往研究可能受限于地理和气候等要素的影响，进而忽视了边缘陡峭的林线区域。基于遥感观测的高山林线提取方法弥补了上述不足。

4.3.2 林线高程随坡度分布特征

由于坡面倾斜程度不同，研究区各山坡获得的太阳辐射量不同，进而使得不同坡面的林线分布特征呈现差异。本文基于International Geographical Union 提出的7级划分标准统计各坡度范围内的林线分布（表4）。结果显示，林线集中分布在15°—55°范围内的陡坡、急坡和急陡坡，并且林线比重随着坡度的上升逐渐增加，区间累计占比高达89.88%。其中，林线Z1、Z2 和Z3 的最高点位于急陂和急陡坡上。然而，各林线在不同坡面区间的平均高程呈现明显分布差异。例如，林线Z1 平均高程的最大值位于垂直坡上，林线Z3的最大平均高程则位于陡坡上。高山林线极少分布在平原和缓斜坡区域，该范围内林线占比仅0.62%。当坡度超过55°时，林线占比显著下降，仅为2.32%，处于急陡坡与垂直坡高坡度范围内的林线最高点高差明显扩大。

表4 不同坡度类型的林线高程分布情况Table 4 Distribution of treeline elevation with different slope types

4.3.3 林线高程随坡向分布特征

山地的不同坡向可以造成辐射强度和日照时数的不同，因此高山林线分布也呈现出坡向分布差异。本文分别对Z1、Z2 和Z3 共3 处林线的8 个坡向高程最大值和平均值进行统计（图10）。图10（a）、图10（b）分别为林线各方向高程的最大值和林线各方向的高程平均值。容易发现各坡向的林线高程分布差异显著，总体呈现南高北低，东高西低的特征。而且，林线分布具有坡向趋向性。例如，林线Z1 和林线Z2 的高程最大值和平均值均位于东南方向，最大值分别为3662 m、3476 m，平均值分别为3297 m 和3164 m。林线Z3 的高程最大值和平均值位于正南方向，分别为3402 m 和3310 m。尽管林线Z3的坡向分布与其他林线不同，但东南方向上高程值接近正南方向。

图10 各个坡向林线高程分布图Fig.10 Treeline elevation by slope direction

4.3.4 林线上、下边界分布特征

秦岭山脉呈东西走向，南、北两侧高山林线空间分布特征差异明显（表5）。例如，北侧林线上边界高度比南侧高66 m，下边界高度比南侧低606 m，南侧和北侧的上边界均位于东南方向，下边界分别位于正西、正东方向。造成这种分布差异可能是由于南、北两侧高山林线的海拔、坡度和坡向存在差异，辐射强度和日照时数的不同，也有可能南、北坡气温时空变化存在差异，北坡的整体升温高于南坡（马新萍，2015）。

表5 南、北林线上、下边界分布特征Table 5 The upper and lower boundary distribution features of the southern and northern alpine treelines

山脊南侧的正东、东北、正北和西北各坡向林线平均高度高于北侧，在其他坡向的平均高度均低于北侧。标准差结果显示南侧这4个方向的地形起伏程度要低于北侧，在西南、正南和东南方向的起伏程度均高于北侧，山脊南、北坡向的林线分布差异明显（表6）。

表6 山脊南、北不同坡向林线高度分析Table 6 Elevation analysis of alpine treeline on different slopes to the south and north of ridge

4.4 讨论

高山林线是植被类型及其生长环境呈现过渡变化的生态界限，准确获取高山林线的分布特征是研究山地植被动态与气候变化关系的重要基础（李迈和等，2005）。秦岭山脉是中国重要的南北地理分界线，气候、植被垂直带谱分布特征明显。为此，本文基于高质量的卫星遥感森林覆盖数据和连通域搜索算法提出了一种快速提取高山林线的方法，在此基础上分析了林线分布与高程、坡度和坡向的关系。研究结果显示，秦岭地区高山林线的高程分布范围为2400—3800 m，3处主要的林线上边界均位于3400 m 以上，这与前人对秦岭高山林线高程的研究结果接近（崔海亭，1983；戴君虎等，2001）。此外，高山林线集中分布在15°—55°的陡坡、急坡和急陡坡区，很少位于坡度低于5°的平原、缓斜坡或大于55°的垂直坡区，表明植被群落交错带较少出现在坡度过缓或过陡的区域。各个坡度的林线高程平均值存在差异，研究区林线Z1 的最大平均高程值位于垂直坡区，而林线Z2 和林线Z3 的最大平均值则位于陡坡区。高山林线高程分布特征呈现显著的坡向差异，呈现南坡林线高于北坡东坡林线高于西坡的特征。

本文的评估结果显示，通过基于遥感的快速提取高山林线的算法获得的林线与研究区Google Earth 影像地图中的实际林线分布基本一致，在林线的完整性和边界连通性等方面表现出众，证明了文中高山林线制图研究的可行性。然而，对比Dang 等（2015）记录的GPS 地面站点观测数据，本文发现GPS坐标位置及其500 m 缓冲区均位于森林内部，明显偏离了高山林线区。这可能是由于高山、树木遮蔽了GPS 信号影响了定位精度所致，也有可能是因为人工操作不当而造成的数据采集错误。

尽管本文提出的林线自动提取算法具有较大的应用潜力，但受森林覆盖数据解译精度和DEM分辨率的限制，高山林线的边界点与实地调查获取的森林边界点位仍存在差异。未来，需要融合更多高分辨的影像数据（比如，高分和PlanetLab）和森林覆盖产品（陈军 等，2017；Chen 等，2019；Liu 等，2021）以提高林线边界识别的精度。此外，现阶段基于遥感影像提取高山林线的研究缺乏统一标准。目前大部分研究主要以山地针叶林的稀疏分布区域作为高山林线（Solaimani等，2011；刘锬等，2016），或直接采用海拔最高的森林斑块连线或最高海拔位置森林点位表示林线（Král，2009），上述方法很难区分林线区不同高度的植被类型及具体的边界。本文结合森林覆盖度信息和高程数据定量描述了植被过渡变化特征，进而使用连通域搜索算法准确地提取秦岭高山林线的空间分布。

林线海拔高度的波动是植被对气候变化的重要响应。林线是当地环境因素（干旱、冰雪等）干扰与调节的结果，反映了植被对气候的适应过程（王襄平等，2004）。一般而言，纬度越低，林线越高。例如，从沿海到内陆方向高山林线逐渐升高（Malyshev，1993）。一些研究也表明山地植被随着全球气候变暖向上迁移，林线呈上升趋势（白红英等，2011；马新萍，2015），向上迁移速率为0.35 m/a（Lu 等，2021）。尽管本文采用的数据空间分辨率（30 m）远大于林线每年向上迁移的高度，但遥感卫星提供的全球近50 a 的高分辨率影像（Landsat），可为揭示过去、现在与未来高山林线对植被生态系统与气候变化的响应机制提供依据。而且，IPCC 第五次评估报告指出全球气候系统变暖持续，将研究高山林线变化格局及其与温度、降水等气候因素之间的关联，以期全面认识高山林线空间分布规律及其驱动机制。

5 结 论

本研究的主要结论如下：（1）与传统野外调查的方法相比，基于森林遥感数据的林线自动提取方法突破了时间、地域等条件的限制，在节约了成本的同时提高了工作效率，有助于促进对人类活动受限地区的林线研究；（2）本研究采用的植被覆盖数据结合了地表时间序列光谱和植被高度信息，进一步提升了林线识别的精度；（3）鉴于遥感技术的大范围对地观测的能力，以及卫星影像数据较高的数据质量和可获取性，本文提出的分块迭代搜索策略适用于寻找不同级别的高山林线分布，可进一步推广至全球高山林线制图研究。

志 谢此次实验的全球森林覆盖数据由马里兰大学Matthew Hansen 课题组提供，数字高程模型由美国宇航局NASA提供。