广西壮族自治区短轮伐期人工林时空分布信息提取

段文胜，陈元鹏，王力，黄妮，贺原惠子，4，张昌赛，4，张阳坚，周泉，4，牛铮，4

1.中国科学院空天信息创新研究院 遥感科学国家重点实验室,北京 100101;

2.中国科学院大学 电子电气与通信工程学院,北京 100049;

3.自然资源部国土整治中心,北京 100035;

4.中国科学院大学,北京 100049

1 引言

森林生态系统在全球碳循环、能量平衡以及物质交换中扮演着重要角色（Brockerhoff等，2008；Williams，2015；葛全胜 等，2008）。人类对自然森林资源的过度掠夺造成了全球森林覆盖面积的锐减以及严重的生态环境问题（Hansen 等，2013；张帅帅 等，2020），短轮伐期人工林SRPs（Short-Rotation Plantations）是指轮伐期短、材质好、造林成活率高的人工经济林，是目前国内主要的木材供应原料。其轮作形成的重复性森林砍伐和造林是区域碳汇变化的重要影响因素（步巧利 等，2020）。SRPs 种植的快速推广虽然缓解了社会发展对木材的需求，但是粗放式的快速发展也带来了一系列生态环境问题（Brockerhoff 等，2008；Williams，2015；黄国勤和赵其国，2014）。因此，SRPs客观科学的监测研究逐渐受到广泛的关注。

目前，虽然有诸多学者利用遥感技术对森林覆盖进行了大量的研究（吴雪琼 等，2010），但针对SRPs 进行遥感信息提取及变化监测的研究十分匮乏。以往人工林信息提取的研究多基于Landsat 系列中高分辨率影像，仅局限在单景影像的小范围内，且精度有限（Qiao 等，2016；蒙良莉 等，2019）。由于SRPs 在两个轮伐期间存在休耕期，仅利用单时相的影像进行SRPs 监测会将休耕期的裸土地块识别为非短轮伐期人工林（Non-SRPs）用地，增加SRPs 用地信息提取的漏分误差（Qiao 等，2016）。时序遥感数据能较好地解决这一问题，例如，Maire 等（2014）基于MODIS（Moderate Resolution Imaging Spectroradiometer）影像重构归一化植被指数NDVI（Normalized Difference Vegetation Index）时序，提高了巴西SRPs提取精度（Le Maire等，2014）。虽然利用高时间分辨率的遥感影像提高了SRPs 信息的提取精度，但中国南方地区受粗放式发展以及地形地貌影响导致SRPs 林分破碎，低分辨率MODIS 影像无法满足SRPs 提取的要求。因此，具有30 多年历史存档数据的中分辨率Landsat 系列卫星影像数据成为更合适的选择。

由于中国南方地区多云多雨，选取每年相同儒略日附近的影像重构长时序遥感数据较为困难（沈文娟和李明诗，2017），而逐像素重构的年度合成数据能够有效的解决此问题。SRPs 的轮伐时序特征，本质上是人为造成的周期性森林干扰和恢复事件。LT（LandTrendr）算法是基于Landsat长时间序列监测森林干扰的主流算法，能有效的识别连续变化的森林干扰事件（Kennedy 等，2010；Shen等，2017）。

本研究以广西壮族自治区（以下简称广西）为研究区，基于GEE（Google Earth Engine）云平台以及近34年的Landsat系列影像，首先采用逐像素合成技术重构34年的归一化燃烧指数NBR（Normalized Burn Ratio）长时间序列数据，然后运用LT 算法对NBR 时序数据进行分割拟合，提取森林区域SRPs的时空分布信息，最后借助Google Earth 高分辨率影像选取样本验证分类提取的精度，分析SRPs 种植面积变化的时空特征与相关因子。本研究的开展对于林业管理、人工林经济估算、生态环境保护以及碳循环研究都有着重要的意义。

2 研究区域和数据源

2.1 研究区域

广西壮族自治区（简称“广西”）位于中国华南地区，20°54′N—26°24′N，104°28′E—112°04′E，陆地面积为2.38×107ha。广西属亚热带季风气候和热带季风气候，气候条件非常适合SRPs的种植。据统计，2018 年广西全区森林面积达到1.48×107ha，森林覆盖率达到62.31%（图1）（Hansen等，2013；兰秀 等，2019）。桉树作为广西的主要SRPs 类型，全区种植面积达2.33×106ha，桉树木材产量约占全国70%（广西壮族自治区人民政府门户网站：http：//www.gxzf.gov.cn［/2021-02-09］）。

图1 研究区域概况Fig.1 Overview of the study area

2.2 实验数据

Landsat 系列卫星传感器数据的波段范围、时空分辨率较为一致，因此，可以很好的重构长时序（Claverie 等，2015；Roy 等，2016；Woodcock等，2008；Zhu 等，2015；邱布布，2017；沈文娟和李明诗，2014；张志杰 等，2015）。本文选择1986年—2019年的Landsat系列卫星（Landsat 5/TM、Landsat 7/ETM+、Landsat 8/OLI）影像作为数据源（表1）。研究中使用的Landsat 影像数据共13051景，其中Landsat 5/TM 5846 景，Landsat 7/ETM+4566景，Landsat 8/OLI 2639景（图2）。在研究的时间范围内，有效观测皆能覆盖广西全境。

表1 Landsat系列数据Table 1 Landsat series data details

图2 研究区Landsat数据量统计Fig.2 Landsat data volume statistics in the study area

2.3 GEE云计算平台

GEE 是当今世界上最先进的专门处理卫星影像等地理空间观测数据的云计算平台。GEE 云端数据库中集成了近40 年的Landsat 系列卫星的历史存档数据，给个人用户提供了强大的算力和云存储空间，同时提供了方便快捷的JavaScript 语言API 接口进行数据处理、算法实现以及结果分析（Dong 等，2016；Gorelick 等，2017；Hansen 等，2013；Padarian 等，2015）。本研究应用GEE 云计算平台进行数据处理，大大减少了数据准备的前期工作，也极大的降低了数据处理与算法实现过程中对本地硬件设备的依赖。

3 研究方法

本研究首先利用GEE 平台进行数据预处理以及结合Google Earth 进行样本数据处理。然后，利用GEE JavaScript API 实现LT 时序分割以及SRPs信息提取。最后，进行SRPs 时空分布制图、精度评价以及趋势分析与相关分析。具体技术路线如下图3所示。

图3 技术路线图Fig.3 Research flowchart

3.1 样本数据处理

在Landsat 影像上SRPs 和Non-SRPs 容易产生“异物同谱”现象（图4），目视解译无法将其区分，因此无法利用Landsat影像进行样本选择。本研究基于马里兰大学M.C.Hansen团队研发的全球森林变化产品GFC（Hansen Global Forest Change v1.7（2000年—2019年））辅助进行样本选择（Hansen等，2013）。

图4 SRPs与Non-SRPs对比Fig.4 Comparison of SRPs and Non-SRPs

利用GFC 产生预选样本点的方法如下（图5）：首先利用GFC 产品中“treecover”波段对研究区进行掩膜处理，生成一个森林覆盖范围的掩膜。然后，利用“loss”波段（定义为林分置换干扰，即从森林状态到非森林状态的变化）和“gain”波段（定义为“loss”的反变化，即由非森林到森林的变化）进行交并运算，把森林覆盖区域分为两个互补的类层：“Both Loss&Gain”（轮伐将导致SRPs种植区域即发生了“loss”又发生了“gain”）和“Not Both”。最后，应用分层随机采样方法，生成两类预选样本点。

图5 样本选取流程图Fig.5 Samples selection flowchart

基于GFC生成的两类样本点并不能确认为SRPs或Non-SRPs，需要对样本点进一步验证。本文基于Google Earth Pro平台的高分辨率卫星影像对样本点进行筛选和验证。结合对SRPs 地面调查的先验知识，SRPs判别依据以下3个特征：（1）轮伐迹象特征（“SRPs—裸土—SRPs”时序特征）；（2）规律的人工种植特征；（3）人工作业痕迹（图6、图7）。通过Google Earth高清影像的验证发现，基于“Both Loss&Gain”生成的预选样本点中仅部分为真实SRPs 样本点。根据随机采样获取的真实SRPs 样本点，人工解译补充了一部分SRPs 样本点。基于“Not Both”生成的样本点中仅小部分为非Non-SRPs样本，进行简单的筛选即可获得Non-SRPs样本。

图6 SRPs样本示例（谷歌地球影像）Fig.6 SRPs sample example（Google Earth image）

图7 Non-SRPs样本示例（谷歌地球影像）Fig.7 Non-SRPs sample example（Google Earth image）

本研究共选取2569 个样本点，包括1154 个SRPs样本点和1415 个Non-SRPs样本点（图6、图7）。将其中100 个SRPs 样本点用于算法参数优化，其他的样本点皆用于精度评价。用于算法参数优化的SRPs样本点可以根据Google Earth 多时相影像判断轮伐年份RP（Rotation Point）。如图6 样本示例，根据前后历史影像，可以判断RP 对应于2015 年。用于精度评价的样本点则是可以确定2016年该像元点是SRPs或Non-SRPs。

3.2 数据预处理与时序重构

3.2.1 掩膜处理与数据标准化

GEE 数据库包含了经过辐射校正后的Landsat系列数据：USGS Landsat 8 Surface Reflectance Tier 1、USGS Landsat 7 Surface Reflectance Tier 1、USGS Landsat 5 Surface Reflectance Tier 1。在调用所需数据后，首先利用GFC 产品的“treecover”波段对森林覆盖区域进行掩膜处理，以筛选出本研究所关注的森林覆盖区的数据。

ETM+与TM 传感器的波段设置是高度一致的，但与OLI 传感器的波段设置存在一定的差别。因此，本文以ETM+的B1、B2、B3、B4、B5、B7 的6 个波段为标准波段（对应OLI 的B2、B3、B4、B5、B6、B7波段），利用Roy等（2016）提出的回归融合模型，对OLI数据进行标准化融合处理。通过对不同传感器的波段进行匹配，确保数据波段一致性。

为了消除噪声像元的干扰，本研究利用GEE数据库中Landsat Surface Reflectance Tier 1 系列产品中的云、阴影、水、积雪标识波段（QA 波段）进行掩膜处理。通过掩膜和标准化，获取了所有的有效观测数据。

3.2.2 构建年度合成数据

为了重构以年为间隔的长时序，本文先采用逐像素合成技术构建年度合成数据。对于多波段的Landsat 数据，研究表明多维中位数合成方法（Medoid）（式（1））有着更好的鲁棒性和更快的执行效率，合成数据具有更好的代表性（Flood，2013）。故本研究选用Medoid 合成法对每年的数据进行处理。

式中，χ={x1，x2，…，xn}，d代表n维欧氏空间中欧氏距离计算函数。

对于植被监测的研究，往往采用的是植物生长季的影像数据代表植物的生长状态（沈文娟和李明诗，2017；Kennedy等，2010，2018）。为了验证选择生长季时间范围内6 月10 日—9 月20 日作为年度合成数据的合成时间范围CDR（Composite Date Range），对100 个RP 样本点对应的不同光谱指数时序应用LT算法，得到每个SRPs样本点时间序列轨迹拟合结果的RMSE（Root Mean Square Error）（图8）。选择的光谱指数（×1000）不论是NBR 还是NDVI，以6 月10 日—9 月20 日作为CDR的RMSE 均小于以全年为CDR 的分割拟合结果，这表明6 月10 日—9 月20 日作为CDR 的分割拟合效果优于全年。

图8 不同CDR应用LT算法RMSE统计Fig.8 RMSE statistics of applying LT algorithm with different CDR

3.2.3 NBR时序重构

本研究利用年度合成数据计算NBR 指数，重构了近34年的NBR时序。本文分析了NBR、NDVI、EVI（增强植被指数，Enhanced Vegetation Index）、NDMI（归一化湿度指数，Normalized Difference Moisture Index）、TCB（缨帽变换亮度分量，Tasseled Cap Brightness）、TCG（缨帽变换绿度分量，Tasseled Cap Greenness）、TCW（缨帽变换湿度分量，Tasseled Cap Wetness）、TCA（缨帽变换角度分量，Tasseled Cap Angle）以及合成数据的波段（B1，B2，B3，B4，B5，B7）共14 个指数或波段在应用LT 时序分割算法后对于RP 识别提取的识别率（图9）。结果显示，不论是生长季合成还是全年合成的时序数据，NBR 对RP具有最高的识别率，表明NBR 指数对于RP提取具有最优的效果。

图9 不同指数/波段应用LT算法对RP的识别率Fig.9 The recognition rate of RP with different indices/bands applying LT algorithm

3.3 LandTrendr算法

LandTrendr时序光谱轨迹分割算法是由美国林务局和美国俄勒冈州立大学共同提出（Kennedy等，2010），可通过该算法生成基于时序轨迹分割拟合的时序数据（图10）。本研究首先利用LT 算法对以年为间隔的NBR 时序进行分割、拟合以及平滑处理，然后基于拟合时序提取SRPs 特有的时序特征，最后进行SRPs 的时空分布信息提取。算法的主要步骤如下：

图10 LandTrendr算法示意Fig.10 LandTrendr algorithm schematic

（1）去除噪声。通过设定阈值筛除原始NBR时序轨迹存在的异常值。

（2）确定潜在分割段顶点。以分段递增的方式对时序轨迹进行分段线性拟合，剔除角度变化过小的分割段，使分割段数达到设置的最大值。

（3）分割段顶点拟合。从第一个分割段开始，在保证分割轨迹连续的前提下以均方误差MSE（Mean Square Error）最小原则确定整个时序轨迹的分割段顶点。

（4）分割拟合模型优化。以MSE 增加最少为原则逐个移除分割点，直到分割段数量为1，最终以置信度为标准选择最优分割拟合模型。

3.4 SRPs制图

利用LT 算法对NBR 长时序进行轨迹分割拟合之后，得到了去除干扰信息并能够突出SRPs轮伐时序特征信息的NBR拟合时序轨迹（图11）。图11中a、b、c这3个阶段构成一个完整的轮伐周期。a阶段对应“SRPs Gain”阶段，为SRPs种植之后的快速生长阶段。c阶段对应“SRPs Loss”阶段，为SRPs成熟林分发生皆伐的阶段。“RP”指两个轮伐期之间的光谱值低谷点。本研究用“SRPs Loss+SRPs Gain”作为时序特征进行SRPs 二分类。为了优化SRPs 的时序特征，本文统计了100 个RP 样本点对应的时序特征。根据统计分析，本研究选定“SRPs Loss”阶段的特征为：持续时间小于4 a，光谱指数变化幅度小于-200。“SRPs Gain”阶段的特征为：持续时间小于4 a，光谱指数变化幅度大于240。

图11 LT拟合的SRPs时序轨迹示例Fig.11 Example of LT fitted SRPs time series trajectory

3.5 精度评价方法

本研究利用基于GFC 和Google Earth 选择并矢量化处理之后的样本数据进行精度评价。采用了混淆矩阵中的总体精度、用户精度、生产者精度以及Kappa系数对算法进行评价验证。

4 结果与分析

4.1 SRPs时空分布制图与精度验证

4.1.1 精度评价结果

SRPs 二分类的用户精度为75.93%，制图精度为79.6%（表2）。相比于以往研究中对SRPs（以往研究中皆是仅对单一SRPs 树种进行研究）进行二分类信息提取的精度都有着不同程度的提升（Thomas 等，2021；Qiao 等，2016；Le Maire 等，2014）。研究结果中存在的漏分误差和错分误差，可能是对RP 的捕捉存一定的识别误差以及非人工林类似时序特征导致的。精度评价结果表明了基于LT算法对SRPs进行二分类是可行的。对比单独使用“SRPs Gain”或“SRPs Loss”作为时序特征进行二分类的评价结果（表2），利用“SRPs Loss+SRPs Gain”作为时序特征进行SRPs 二分类提取的精度是最佳的。

表2 SRPs基于不同时序特征（SRPs Loss+SRPs Gain/ SRPs Gain/ SRPs Loss）的二分类结果混淆矩阵Table 2 Confusion matrix of SRPs binary classification results（SRPs Loss+SRPs Gain/ SRPs Gain/ SRPs Loss）

4.1.2 SRPs时空分布

利用SRPs 独特的轮伐时序特征，本研究进行了广西近30年的SRPs种植分布制图（图12）。近30年，广西SRPs种植从零星分布到大面积分布。广西东部南部SRPs 分布最为集中，且种植面积增长最为迅速。广西西部、北部地区受喀斯特地质地貌的影响，SRPs 的增长较南部和东部缓慢。基于广西地级市行政范围进行区域统计分析（图13），河池市SRPs 种植面积最大，其次百色市（其行政面积最大）和崇左市，皆超过广西SRPs 种植总面积的10%。南部的玉林市行政面积远小于桂林市，但SRPs 种植面积却相当，也放映了南部地区比北部地区的SRPs分布更为密集。

图13 广西地级市SRPs种植占比Fig.13 The proportion of SRPs planted in prefecture-level cities in Guangxi

广西SRPs 种植分布存在较强的空间分布特征（图14，图15）。从图14 可以看出，随着海拔的升高，SRPs 种植分布呈现明显的下降趋势。其中，SRPs种植分布大部分在海拔500 m以下，海拔200—300 m 分布面积最大。主要原因是海拔越高，气候温度越低，不利于SRPs 的快速生长从而无法进行短轮伐期轮作。从图15 的坡度分析看出，SRPs 都是种植在具有坡度的上坡地，基本没有种植在适合作为耕地的平地。其中，SRPs 种植主要集中分布在坡度为20°左右坡地，当坡度大于40°时也只有极少量的分布。主要是由于坡度较大的山地不适合人工作业。

图14 广西SRPs种植面积占比随海拔变化分布Fig.14 The distribution of SRPs planting area proportion at different altitude in Guangxi

图15 广西SRPs种植面积占比随地表坡度变化分布Fig.15 The distribution of SRPs planting area proportion on different surface slope in Guangxi

4.2 SRPs演变趋势与驱动机制

广西SRPs 的种植面积近30 年一直呈现稳定快速增长的趋势（图16），其增长趋势与人工造林面积增长趋势一致（国家林业和草原局，2019）。1990 年SRPs 种植面积仅1.93×105ha，到2019 年达到了4.04×106ha，年平均增长速度达到1.33×105ha。根据中国林业统计年鉴数据，广西林业生产总值呈快速增长趋势（图17）。林业生产总值与SRPs种植面积之间的相关性为r=0.83，p＜0.001。广西林业生产总值与SRPs 种植面积之间高度相关，这也表明了SRPs对于林业经济发展的重要影响。

图16 广西SRPs种植面积与人工造林面积变化趋势Fig.16 The change trend of SRPs planting area and artificial afforestation area in Guangxi

图17 1992年—2018年广西林业生产总值变化趋势Fig.17 The change trend of Guangxi’s forestry production value from 1992 to 2018

从图17 可以发现，广西的林业生产总值在2007 年左右出现了明显的拐点，2010 年以后林业生产总值开始快速增长。其主要原因是：（1）广西于2007年1月1日开始实施《中华人民共和国森林法》。《中华人民共和国森林法》进一步强调了森林保护的规定，对森林经营管理、植树造林、森林采伐有了更加明确的规定（广西壮族自治区人民政府门户网站：http：//www.gxzf.gov.cn［/2021-02-09］）。（2）广西的林业“十一五”规划提出从林业大省区要变强省区的目标，让林业成为广西的支柱产业。同时强调了经济林、森林生态效益补偿林、速丰林的重要性，并大力引进与发展森工企业（国家林业与草原局：http：//www.forestry.gov.cn［/2021-02-09］）。

从近30 年广西各地级市SRPs 种植面积变化趋势（图18）来看。其中河池市在广西的SRPs 种植面积始终保持最大占比，年均增长面积也最大；北海市始终保持最小占比，且年均增长面积最小。河池市的SRPs 种植面积有两个突出的快速增长时间段（1992 年—1994 年与1998 年—2000 年），百色市、南宁市、桂林市、梧州市、来宾市、贵港市以及贺州市的增长趋势都较为平稳。河池市、崇左市以及玉林市在2005 年左右，SRPs 种植面积增速开始明显放缓，而柳州市种植面积增速有一定程度的增加。

图18 广西各地级市SRPs种植面积变化趋势Fig.18 The change trend of SRPs planting area in prefecturelevel cities in Guangxi

5 结论

本文以中国SRPs种植面积最大的广西壮族自治区为研究区，探讨了基于时序分割算法LandTrendr进行SRPs 时空分布信息提取的方法，得出以下结论：

（1）利用LandTrendr 时序分割算法，能有效地对时序分割拟合，并突出SRPs特有的RP时序特征（SRPs Loss+SRPs Gain），证明了LandTrendr 算法对于SRPs 种植导致的森林干扰和快速恢复事件有着很好的监测能力。

（2）基于时序特征的SRPs 信息提取中，应用逐像素合成技术重构的NBR光谱指数效果最优。

（3）广西SRPs 种植面积在近30 年一直呈快速稳定增长的趋势，与林业产值变化趋势高度相关，也说明了经济增长对SRPs种植的需求。

（4）从空间分布来看，广西东部、南部SRPs分布较为集中，其主要分布在海拔500m 以下的低海拔地区以及地表坡度在20°左右的坡地。其中河池市是广西SRPs种植分布面积最大的地级市。

志 谢本文实验数据获取及数据处理得到谷歌公司免费平台（Google Earth Pro，Google Earth Engine）的支持，在此表示衷心的感谢！