基于高时空分辨率可见光遥感数据的桉树种植面积提取方法研究

Extraction of Eucalyptus Planting Area Based on High Spatio-Temporal Resolution Visible Remote Sensing Data

摘要: 为了缓解国内木材市场的供求矛盾，我国高度重视速生林基地建设，而桉树作为世界三大速生树种之一，因其易于种植、适应性强且轮伐周期短的特点而具有较大的经济价值。大量、单一的速生桉代替了原有的自然林和耕地，造成生物多样性损失、局部水资源短缺等生态环境问题。对该区域桉树林的种植面积和种植结构变化的精确监测是客观评价该地区桉树林种植与生态环境变化关系的关键。本文与澜沧地区为例，分别对Sentinel-2和MOD09GA进行植被指数提取和缨帽变换，并基于时空数据融合技术，分别融合植被指数和缨帽变换指数，获取高时空分辨率NDVI数据和TCA数据，提取桉树林的植被变化特征，实现高精度的桉树林分布制图。结果表明：1) 基于时空融合数据提取的橡胶林物候变化特征能够实现桉树林的识别，识别精度可以达到80.12%，Kappa系数达到0.71；2) 用NDVI指数数据分类时，能够获取比TCA指数数据更高精度的分类结果，表明植被指数数据在高时空数据融合及植被遥感应用中有较好的应用前景。

Abstract: In order to alleviate the contradiction between supply and demand in the domestic timber market, China attaches great importance to the construction of fast-growing forest bases. As one of the world’s three fastest-growing tree species, eucalyptus has great economic value because of its easy planting, adaptability and short rotation cycle. A large number of single fast-growing rafts have replaced the original natural forests and cultivated land, causing ecological and envi-ronmental problems such as biodiversity loss and local water shortages. Accurate monitoring of the planting area and planting structure changes of the eucalyptus forest in this area is the key to objectively evaluate the relationship between the plantation of the eucalyptus forest and the ecological environment in the region. In this paper, we take sentinel-2 and MOD09GA in Lancnag region as examples to extract vegetation index and transform tassel cap, and based on the spatio-temporal data fusion technology, the vegetation index and the cap transformation index are respectively combined to obtain the high spatial-temporal resolution NDVI data and TCA data, to extract the characteristics of vegetation changes in the eucalyptus forest and to achieve high-precision eucalyptus distribution mapping. The results show that: 1) The phenological change characteristics of rubber forest based on spatio-temporal fusion data extraction can realize the identification of eucalyptus forest; the recognition accuracy can reach 80.12%, and the Kappa coefficient reaches 0.71. 2) When using NDVI index data classification, it can obtain higher precision classification results than TCA index data, indicating that vegetation index data has a good application prospect in high-temporal data fusion and vegetation remote sensing ap-plications.

1. 引言

桉树是生长在澳大利亚的土著树种，又名尤加利树，是桃金娘科、桉属植物的统称。桉树的轮伐周期短，易于种植、且适应性较强，因此具有较大的经济价值。目前桉树人工林种植面积较大的有中国、西班牙、印度、智力、葡萄牙等 [1] [2]。据2013年有关数据统计、我国桉树人工林面积达450万hm² [3]，主要分布于广西、广东、福建、海南等热带亚热带地区，其中以广西种植面积最大 [4]。近年来，不少农户为了追求经济效益，炼山种植和耕地种植桉树的现象频繁发生，原有的自然林被大量的桉树林所代替，造成生物多样性损失和局部水资源短缺等生态环境问题。因此，及时准确地掌握桉树的种植面积、种植结构的变化是客观评价桉树种植与生态环境变化关系的关键。

森林类型图的获取方法包括野外调查和遥感制图两种。野外调查方式的成本高、周期长，不利于大面积森林类型图的制作。遥感技术因其快速、可重复、覆盖地域广等优势，已经成为当前获取大面积森林类型空间分布信息的主要手段 [5]。目前国内外学者在森林树种遥感识别方面已经取得了一定的研究成果。伍静基于TM数据的速生桉遥感分类技术研究，对TM数据进行植被指数提取、主成分分析、缨帽变换和最小噪声分离，共生成14个特征波段，加上原始的TM影像的6个波段，作为讨论速生桉分类最优波段组合的基础，运用一般监督分类法、BP神经网络分类法、决策水分类法、混合像元线性分解法，将地类分为水域、裸地。桉树、杉木、幼林、经济林 [6]。牟智慧等运用TM影像，基于面向对象的方法，利用植被光谱信息和树种生长环境特征提取单一数据信息，完成桉树林的信息提取 [7]。雷光斌等基于多源多时相遥感影像，获取不同森林类型的季相节律信息来提高深林类型的制图精度 [5]。章辰宇等利用30年的Landsat影像数据获取NDVI时间序列，并综合考虑桉树林周期性轮伐的信息来完成桉树林分布制图 [2]。付晓等在Landsat数据的基础上，通过不同植被的波谱特征分析，选择4、5、3波段组合合成遥感图像，进行纸浆源地桉树资源的信息提取；在地面验证和训练样地的支持下，进行图像分类结果归类汇总，获得不同时相桉树资源的面积与分布 [8]。王学成等基于Raypideye高分辨率遥感影像和面向对象的方法，提取桉树林的分布信息 [9]。梁文海等基于不同的面向对象分类方法，对GF-2数据中的桉树信息信息进行提取，并对比了不同方法提取结果的可靠性 [10]。曾志康等以高分一号影像和谷歌高分辨影像为数据源，开展了面向水源地保护的地块尺度桉树遥感识别 [11]。与单一时相的中分辨率数据相比，高分辨率数据在桉树林识别中能够更加充分地利用桉树林的纹理信息，并有效提高桉树林的识别精度。多源遥感信息能综合利用多源遥感数据优势，提高分类效果 [12]。目前关于应用多源遥感数据提取树种的相关研究有很多，例如高书鹏等基于高时空分辨率可见光遥感数据提取橡胶林 [13]，李磊等基于BEMD多源遥感影像融合数据提取漓江源竹林 [14]，邬明权等利用遥感数据时空融合技术提取水稻种植面积 [15]，但用遥感提取桉树的研究很少，而且是一个难题，桉树在中低分辨率遥感影像上的光谱和自然林相近，较难区分。目前国外对桉树的研究集中在叶面积指数、种植密度、生物量等试验的研究 [15] [16]，而对于桉树提取的研究较少。

本研究分别对Sentinel-2影像(10 m)和MOD09GA影像(250 m)进行植被指数提取和缨帽变换，并基于时空数据融合算法，获取高时空分辨率的NDVI数据和TCA数据，提取桉树林的植被变化特征，实现高精度的桉树林分布制图。这不仅将桉树人工林的纹理和多光谱特征相结合，而且还有效利用了桉树林的时相特征，有效提高提取精度。

2. 研究区与数据

2.1. 研究区概况

研究区位于云南省澜沧拉祜族自治县，位于中国西南边陲(图1)。易于种植、适应性强且轮伐周期短

Figure 1. Geographical location and sample point distribution of the study area

图1. 研究区地理位置及样本点分布图

等特点使得桉树具有较大的经济价值，不少农户为了追求经济效益，炼山种植和耕地种植桉树的现象频繁发生，刺激了桉树种植面积的扩张，区域内桉树的种植面积迅速发展。因此，研究区内桉树林具备云南地区桉树林时空分布与发展的代表性。本文在此研究区内，研究不同传感器数据在包含桉树林地的复杂地表环境下的时空数据融合精度。

2.2. 数据

1) Sentinel-2数据。Sentinel-2A卫星于2015年6月23日，从法属圭亚那库鲁发射场由“织女星”(Vega)运载火箭发射升空。时隔不到两年，Sentinel-2B卫星在2017年3月7日同一处发射升空。这两颗卫星是“哥白尼”计划下的多光谱成像卫星，用于全球高分辨率和高重访能力的陆地观测、生物物理变化制图、监测海岸和内陆水域，以及风险和灾害制图等。本研究选取研究区内无云无像元缺失的数据，该数据是与MOD09GA数据同时段的10 m空间分辨率影像，源自美国地质调查局(https://glovis.usgs.gov/)。Sentinel-2数据的大气校正均采用Sen2Cor大气校正，Sen2cor是由ESA发布，并专门用于生成L2A级数据产品的插件。

2) MOD09GA数据。中分辨率成像光谱仪(MODIS)是搭载在Terra和Aqua卫星上的一个重要载荷。Terra卫星于1999年12月18日发射成功，于每天地方时上午10:30时过境，因此也把它称作地球观测第一颗上午星(EOS-AM1)。Aqua卫星于2002年5月4日发射成功，于每天地方时下午1:30过境，因此称作地球观测第一颗下午星(EOS-PM1)。Terra卫星的数据产品时间范围从2000年至今，Aqua卫星的数据产品时间范围从2002年至今。在Terra和Aqua的相互配合下，每1~2天就可对整个地球表面进行重复观测。本研究选择2017年、2018年、2019年的影像(500 m空间分辨率MODIS地表反射率产品，Daily)，源自NASA数据共享平台(https://search.earthdata.nasa.gov/)。其中2017年和2019年的MOD09GA数据用于修复2018年中存在缺失或云覆盖的MOD09GA数据，以便于生成高质量的2018年的MOD09GA daily数据。

3) 验证数据。所有验证数据均来自Google Earth高清影像数据。优于训练样本的选取直接关系到分类结果，为了保证选取样本的代表性以及随机性，本文根据选取样本的空间分布情况以及各地类的面积比例，确定各个地类所选取的样本点数量，然后结合Google Earth高清影像目视解译的矢量文件，按比例随机生成各地类的样本点(表1)。

Table 1. List of Sentiel-2 and MODIS remote sensing image data selected

表1. 所选用的Sentinel-2与MODIS遥感影像数据列表

3. 方法

3.1. 融合算法

Zhu等人提出的时空数据融合算法ESTARFM，因其应用普遍、可靠程度高，容易实现等原因，而被广泛使用 [17]。因此，本研究采用时空数据融合算法ESTARFM，融合高空间分辨率数据和高时间分辨率数据，最终获取具有高时空分辨率的NDVI和TCA数据提取桉树林。该算法需要完全无云且质量好的两对遥感数据作为输入数据，而MODIS反射率数据能够与高空间分辨率Sentinel-2A反射率数据形成同日观测数据对，具体算法如下：

$F\left(x,y,t_0,B\right)=a\times C\left(x,y,t_0,B\right)+b$ (1)

$F\left(x,y,t_p,B\right)=F\left(x,y,t_0,B\right)+a\times \left(C\left(x,y,t_p,B\right)-C\left(x,y,t_0,B\right)\right)$ (2)

$F\left(x_{\omega /2},y_{\omega /2},t_p,B\right)=F\left(x_{\omega /2},y_{\omega /2},t_0,B\right)+{\displaystyle \underset{i-1}{\overset{N}{\sum }}{W}_i\times V_i\times \left(\left(x_i,y_i,t_p,B\right)-C\left(x_i,y_i,t_0,B\right)\right)}$ (3)

$W_i=\left(1/D_i\right)/{\displaystyle \underset{i=1}{\overset{N}{\sum }}\left(1/D_i\right)}$ (4)

$D_i=\left(1-R_i\right)\times \left(1+\sqrt{{\left(x_{\omega /2}-x_i\right)}^2+{\left(y_{\omega /2}-y_i\right)}^2}/\left(\omega /2\right)\right)$ (5)

公式(1)是ESTARFM算法的高分辨率影像反射率与低分辨率影像反射率的基本关系描述，F，C分别为高分辨率和低分辨数据得到的反射率，两者为线性关系。其中a、b分别直线的斜率和截距，(x, y)为像元位置，t₀为影像获取时间，B为波段。要解方程得到a和b，必须需要两个数据对，然后高时间分辨率数据去预测同一时刻的高空间分辨率数据(公式2)，t_p，t₀为两期遥感影像的成像时间，a是任意像元(x, y)的转换系数。如果考虑到从邻近的光谱相似像元所获得的辅助信息及这些像元本身的权重，在一个

给定的移动窗口内，中心像元的预测反射率就为(公式3)， $\left(x_{\omega /2},y_{\omega /2}\right)$ 为中心像元坐标， $\left(x_i,y_i\right)$ 为第i个

光谱相似像元坐标，N为光谱相似像元个数。W_i(公式4)为第i个相似像元的权重，V_i为相似的像素转换系数，R_i (公式5)表示待融合像元(x, y)在前后两个时期的像元值的相关性。

3.2. 验证数据

1) 融合数据的精度验证。本文用ESTARFM融合得到的数据与获取的当日Sentinel-2数据进行相关性分析，例如验证ESTARFM数据为2018年第55天的数据(DOY055)，融合输入的数据为2018年第25天的数据(DOY025)与2018年第55天的数据(DOY089)的MODIS和Sentinel-2反射率数据。评价过程中将融合得到的DOY055 NDVI数据和TCA数据，与Sentinel-2数据计算得到的NDVI和TCA进行对比分析。评价指标采用决定系数R² (R为相关系数)和均方根误差(root mean square error, RMSE)，决定系数R²越接近于1，证明融合数据的精度越高，反之，融合精度越低。

由于融合数据的时间间隔越长，融合精度越低。因此。本研究分别对Sentinel-2和MODIS数据的NDVI和TCA数据做两次融合，其中第一次融合的输入数据为2017年第305天的数据(DOY305)与2018年第25天的数据(DOY025)的NDVI或TCA数据，验证ESTARFM数据为2017年第355天的数据(DOY355)，时间间隔为85天；第二次融合的输入数据为DOY025 (2018年)与DOY089 (2018年)的NDVI或TCA数据，验证ESTARFM数据为2018年第55天的数据(DOY055)，时间间隔为64天，融合精度评价结果如图2所示。

2) 橡胶林分类结果精度验证。验证数据为2018年的Google Earth高清影像。选取了研究区不同地类的总共967个样本点，其中的1/3的样本(约317个，各类地物训练样本不少于20个)作为分类器的训练样本，分类结果如图3，其余的650个样本点用来验证分类结果，图1为样本点的空间分布。精度评价采用误差矩阵评价方法计算各类别的生产者精度(Producer’s Accuracy, PA)、用户精度(User’s Accuracy, UA)以及分类的整体精度(Overall Accuracy, OA)及精度可信度参数Kappa系数，评价结果如表2。

Figure 2. Evaluation of fusion accuracy

图2. 融合精度评价图

Figure 3. Classification results

图3. 分类结果图

Table 2. System resulting data of standard experiment

表2. 标准试验系统结果数据

4. 结果

1) 图2融合精度评价结果显示，当输入数据的时间间隔长度为85天时，验证ESTARFM数据与Sentinel-2提取的NDVI数据的相关系数R²为0.8334，RMSE为0.0597；验证ESTARFM数据与Sentinel-2提取的TCA数据的相关系数R²为0.8103，RMSE为0.0573。当输入数据的时间间隔长度为64天时，验证ESTARFM数据与Sentinel-2提取的NDVI数据的相关系数R²为0.9076，RMSE为0.0595；验证ESTARFM数据与Sentinel-2提取的TCA数据的R²为08834，RMSE为0.0631。这表明：

在使用时空数据融合算法融合高时间分辨率数据和高空间分辨率数据时，输入数据的时间间隔长度会影响数据融合的精度，即输入数据的时间间隔越长，融合精度越低。此外，当输入数据的时间间隔相同时，输入数据的不同，得到的融合精度也会不同。本研究分别用时空数据融合算法将30 m空间分辨率与500 m的NDVI和30 m空间分辨率与500 m的TCA进行融合，最终的到30 m空间分辨的NDVI Daily数据产品和30 m空间分辨率的TCA Daily数据产品，并且融合精度评价结果显示：当输入数据的时间间隔相同时，NDVI的融合精度要高于TCA。因此，在使用融合算法融合时空数据时，要尽可能充分利用所有可用数据，以便于缩短输入数据的时间间隔，提高融合的精度。

2) 表2分类精度评价结果显示：用NDVI时间序列进行桉树林信息提取的总体精度、Kapp系数以及各类别的生产者精度和用户精度均高于TCA，其中桉树林的分类精度，其次是非植被。桉树林分类精度低的原因主要是桉树林与自然林相似程度极高，难以区分，且桉树没有特别明显的物候特征(例如，橡胶冬季落叶)。而非植被在本研究区内非常少，小面积的水体依然能够较好的识别，建物和道路由于过于破碎而难以被识别。