天然橡胶不仅是重要的工业原料,更是国防和工业建设不可或缺的战略能源。橡胶树作为重要的经济作物在热带地区农民增收和经济发展中发挥着越来越重要的作用。东南亚地区产出的橡胶占全球总产量90%以上;其中泰国、印度尼西亚、马来西亚三国占全球割胶面积的72%[1]。掌握该地区橡胶种植、长势等情况,对我国政府部门及时制定或调整产业发展政策,保障天然橡胶资源供应安全,维护橡胶产业的健康发展具有重要意义。利用遥感技术提取天然橡胶林种植的空间分布信息,获取橡胶林种植的时空变化特征,是开展橡胶林长势、灾害、产量遥感监测的重要前置条件[2,3,4]。
在遥感提取国内橡胶林种植分布方面,张京红等[5]利用Landsat-TM卫星影像采用监督分类的方法提取了海南岛2008年天然橡胶种植面积信息。田光辉等[6]利用MODIS EVI数据构建橡胶树物候特征参数提取了海南天然橡胶林分布情况。杨红卫等[7]利用高分辨率遥感影像纹理和多光谱特征提取了海南岛某农场的精细化分布。在我国另一橡胶主产区云南西双版纳地区,廖谌婳等[8]、余凌翔等[9]、Senf等[10]也分别通过不同方法开展了橡胶林分布提取,并分析了橡胶林的扩张情况[11]。针对我国境外橡胶林分布情况,李阳阳等[12]利用MODIS数据及橡胶林的物候特征提取了老挝北部地区橡胶林分布及扩张情况。李宇宸等[13]应用决策树方法对Landsat OLI多时相遥感影像数据表征的橡胶树物候特征提取了中老缅交界区橡胶林分布情况。从提取采用的数据和分类依据来看,多数研究采用MODIS中分辨率影像多时相植被指数作为区分橡胶林与天然林的差异特征[12,13,14,15,16,17],这种方法适用于大范围橡胶林分布提取;还有采用雷达卫星[18]、高分辨率遥感影像纹理[4, 8-9]作为分类依据,提取小范围、精细化分布情况的报道。以上研究主要针对我国海南、云南和周边国家橡胶种植区,而对橡胶种植最多的泰国、印度尼西亚、马来西亚地区分布提取的研究较少。
东南亚橡胶主产区（泰国、印度尼西亚、马来西亚）空间范围大、遥感影像产品数据量巨大,传统的数据下载、单机运算成为提取该地区橡胶林分布情况的瓶颈。因此本研究利用Google Earth Engine（GEE）遥感数据处理云平台[19],以Landsat OLI及MODIS NDVI数据作为数据源,融合光谱、物候2种特征作为分类依据,通过云计算技术解决大尺度、长时序的海量遥感数据处理问题,快速、准确提取东南亚地区橡胶林分布。
1 材料与方法
1.1 研究区概况
研究区包括橡胶种植面积最广、产量最高的泰国、印度尼西亚、马来西亚三国,位于北纬20~南纬10°、东经96~140°之间。泰国地处中南半岛,为热带季风气候,年均气温24~30 ℃,常年温度不低于18 ℃,平均年降水量约1000 mm;11月至次年2月受较凉的东北季风影响比较干燥,3—5月气温最高,可达40~42 ℃;7—10月受西南季风影响,是雨季;农作物一般在雨季播种,旱季收获。印度尼西亚和马来西亚地处马来群岛,属热带雨林气候,终年高温多雨,年平均温度25~ 27 ℃,无四季分别,北部受北半球季风影响,7~9月降水量丰富,南部受南半球季风影响,12月、1月、2月降水量丰富,年降水量1600~2200 mm。
1.2 数据源及预处理
研究区地处热带,云量较多,遥感影像数据质量较差。因此选用GEE平台提供的Landsat 7 2012—2014年大气层顶影像产品融合数据集[20]。该数据集为GEE对经NASA辐射定标、几何校正、云雪阴影掩膜处理的大气层顶反射率产品（TOA）用Simple compose方法融合而成。
植被指数选用MODIS MOD13Q1数据集（2011—2016年）。该数据集是由美国国家海洋和大气管理局先进的中分辨率成像光谱仪（MODIS）生成的归一化植被指数（NDVI）,空间分辨率为250 m,时间分辨率为16 d[21]。
纯净、典型的样本是高精度分类的关键。为保证具有足够数量和质量的样本点用于分类以及精度验证,本文利用Google Earth高分辨率影像数据目视解译具有代表性、典型性的纯净像元作为样本点。
1.3 方法
1.3.1 光谱指数计算 热带地区遥感分类提取的难点在于获取无云数据以及从众多植被覆盖中探寻橡胶林的特征差异。遥感分类中通常利用光谱波段特征和光谱指数对土地进行分类,划分为森林类、橡胶林、农田类、水体类、人工建筑（裸地）类。光谱波段特征光为Landsat 7三年大气层顶光谱融合数据,其中1~5、7波段为大气层顶反射率,波段6为大气层顶亮温。为增加光谱区分能力,采用了2种光谱指数归一化植被指数（normalized difference vegetation index, NDVI）和归一化水体指数（normalized difference water index, NDWI）。其中,归一化植被指数（NDVI）是由红光波段和近红外波段构成的波段组合,该指数可以反映作物长势、茂密程度,以及植被分布情况。具体计算公式[12]为：
$\text{NDVI}=\frac{{{B}_{nir}}-{{B}_{vis}}}{{{B}_{nir}}+{{B}_{vis}}}$
Bnir和Bvis分别为近红外波段（B5）和可见光红波段（B4）的反射率。
归一化水体指数（NDWI）,通过绿光波段和近红外波段间组合可以有效抑制其他类型的地表覆被而有效的凸显水体信息,利于区分水体信息。具体NDWI计算公式[12]为：
$\text{NDWI}=\frac{{{B}_{green}}-{{B}_{nir}}}{{{B}_{green}}+{{B}_{nir}}}$
式中Bgreen为可见光绿波段（B3）的反射率。
1.3.2 植被物候 研究区域的橡胶与天然森林都属于常绿植被,但橡胶林区别于其他天然森林的特征是橡胶林具有典型的物候变化。典型物候变化直接体现在叶的凋落和抽发上。研究区橡胶林12—2月开始落叶,第1蓬叶抽发期为3—4月,第2蓬叶抽发期开始时期为5月,此时橡胶林进入夏花期,第3蓬叶抽发期为7—8月。所以在2月落叶后,3—4月第1蓬叶抽发时,植被指数有一低值时间段。为反映橡胶林物候的变化情况,对2011—2016年1—6月（12期16天合成数据）的MODIS MOD13Q1中NDVI通过simple compose方法求中间值,利用该序列反映橡胶林最典型的落叶-新叶抽生物候特征。
1.3.3 分类算法 本文采用的分类回归树（CART）是一种决策树分类器。该分类器考虑地物的多重属性,综合考虑各属性的重要程度,是一种分类精度较高的常用遥感分类模型。CART分类器的流程为先从已知样本中归纳总体规律,将训练样本属性分成多个训练元组,计算这些元组分裂前后基尼指数,找到最好的分裂准则及分裂值,然后将根节点一分为二,依此进行递归运算,最终拟合出一个适合样本数据的最优二叉树[22]。本文采用的CART分类器由GEE内置的smileCart函数实现。
1.3.4 精度评价 为掌握Cart分类模型的准确性,将样本点70%用于建立分类模型,30%用于精度评价。利用混淆矩阵检验分类精度,计算的总体分类精度（overall accuracy,OA）与Kappa系数有关。Kappa系数的大小能反映出提取的分布与真实地表覆盖物的空间一致性。当Kappa系数小于0.4时,表明一致性较不理想;当Kappa系数为0.4~0.6时,说明二者一致性效果较一般,当Kappa系数大于0.6时,说明分类结果与真实分布有较强的一致性[12]。二者计算公式如下：
$\text{OA}=\frac{\mathop{\sum }_{i=1}^{k}{{x}_{ii}}}{x}$
$\text{Kappa}=\frac{x\mathop{\sum }_{i=1}^{k}{{x}_{ii}}-\mathop{\sum }_{i=1}^{k}{{x}_{i*}}{{x}_{*i}}}{{{x}^{2}}-\mathop{\sum }_{i}^{k}{{x}_{i*}}{{x}_{*i}}}$
式中, xii为混淆矩阵中的第i行i列中的数, x为验证数据集总数,xi*、x*i分别为混淆矩阵中的第i行和第i列总样本数量,k为分类数。
1.3.5 技术路线图 本文采用如图1技术路线对源数据处理、构建分类模型、评价经度、提取橡胶林分布。
2 结果与分析
2.1 典型橡胶林影像
在高分辨率的遥感影像下,可以清晰发现橡胶林、大型农作物农田、天然森林具有显著的差异（图2）。天然森林植被密度高,亮区连片、阴影呈散落斑点状;橡胶林及大型作物农田在种植时被人工排列成整齐的行列,具有明显的行列纹理特征。橡胶林与大型农作物相比,都成行成列,但大型农作物的农田一般行列间距更大;在冠层形态上也有显著不同,橡胶林在一行中更密集。若对橡胶林与大型农作物存在不确定,查阅Google Earth历史影像即可发现不同：橡胶林历史影像变化不大,均为行列种植的大型树木;农田历史影像变化差异相对更大,这是由于农作物生长更快、种植的作物更换更加频繁。其他水体、城镇用地（裸地）在影像上具有更加清晰明显的差异。通过目视解译识别这些差异可以选定典型样本区,共选取橡胶样本像素点107个,天然森林113个,水体30个,人工建筑用地（裸地）240个,农田223个。
2.2 不同地表覆盖类型的光谱信息对比
遥感影像分类过程中首先考虑研究区目标地表覆盖类型,本文将地表覆盖类型划分为农田类、天然森林类、橡胶林类、水体类及城镇用地类（图3）。通过统计典型地表覆盖纯净像元的光谱信息不同波段的反射率可以发现：水体和城镇用地在波段1~3显著高于植被覆盖的区域,而在波段4又低于植被覆盖区域。这个差异可以准确区分植被覆盖区和水体、城镇用地。波段6为亮温数据,几种地表覆盖类型差距较小,普遍在300 K附近。农田在波段4的反射率明显小于橡胶林、天然森林,易与天然森林和橡胶林区分。橡胶林与天然森林二者曲线走势及数值都非常相近,单纯依靠不同波段的光谱特征区分橡胶林与天然森林较为困难。因此,需进一步选择能反映物候差异的NDVI时间序列才能较好识别橡胶。
2.3 森林与橡胶林NDVI时间序列分析
不同的植被在不同季节或生育期表现出不同的生理特征,比如长叶、落叶,这些变化能够通过多时相植被指数时间序列的变化曲线来表示。本文提取了上半年12期16 d合成的NDVI数据中间值（图4）,对比发现农田在所有时相的植被指数均小于天然森林及橡胶林;天然森林NDVI值在不同时间普遍稳定在0.7~0.8附近;橡胶林在第4~6期有一个低值时段。橡胶林在春季一般有落叶、第1蓬叶抽发的物候现象。在中国区域橡胶林一般12月份开始落叶,至2月落叶过程完成,3—4月第1蓬叶抽发完成。橡胶林落叶、新叶抽发时间主要受气温影响,气温越高落叶时间越晚,新叶抽发速率越快。分析东南亚地区典型橡胶林NDVI曲线可以发现,该地区的橡胶林与中国区域具有相似的物候特征,在2月下旬至3月间完成落叶、新叶抽发过程。因此可以通过NDVI时间序列反映出的橡胶林与天然森林不同物候特征区分二者。
2.4 分类精度
基于目视解译的典型样本数据,采用分类回归树CART方法对典型样本进行分类,并根据分类结果建立混淆矩阵（表1）。依据精度评价计算公式（3）、公式（4）,该分类结果的总体分类精度为95.8%,Kappa系数为0.94。从分类精度评价指标来看,这个分类结果精度能够满足空间分析与实际应用需求。从橡胶林样本生产者精度来看,达到93.8%;用户精度达到88.2%,也达到较高的水平。从橡胶林分类误差的来源看,主要发生在农田、森林和橡胶林之间,特别是森林与橡胶林之间。总体而言,这个分类模型的精度满足大范围提取橡胶林分布的要求。
表1 东南亚土地覆盖分类混淆矩阵Tab. 1 Land cover classification confusion matrix in Southeast Asia |
| 样本Sample | 农田Farmland | 森林Forest | 橡胶林Rubber forest | 水体Water | 建筑用地Building land |
|---|---|---|---|---|---|
| 农田 | 64 | 0 | 1 | 0 | 2 |
| 森林 | 1 | 30 | 3 | 0 | 0 |
| 橡胶林 | 0 | 2 | 30 | 0 | 0 |
| 水体 | 0 | 0 | 0 | 10 | 0 |
| 建筑用地 | 0 | 0 | 0 | 0 | 72 |
2.5 东南亚橡胶林分布状况
图5 东南亚主产区橡胶林分布示意图基于自然资源部标准地图服务网站审图号为2020（4399）号的标准地图制作,底图无修改。 Fig. 5 Rubber plantation distribution in Southeast Asia Based on the standard map with the approval number of 2020 (4399) on the standard map service website of the Ministry of Natural Resources, the base map has not been modified. |
3 讨论
本研究利用GEE云计算平台,通过目视解译高分辨率遥感影像选取典型样本区,分析Landsat影像和MODIS NDVI时间序列差异,建立CART分类回归树分类模型,提取橡胶林分布信息。目视解译发现在高分辨率遥感影像下,橡胶林具有独特的行列纹理特征,以区别于其他植被覆盖区。天然森林与橡胶林在Landsat多波段光谱曲线较为相似,与其他地物特征显著不同。MODIS NDVI时间序列反映的植被物候特征表明,东南亚地区的橡胶林与我国境内的类似,在2月下旬至3月间完成落叶、新叶抽发,相应时段的NDVI有一低值时段。利用以上影像特征建立的CART分类回归树模型,分类精度达95.8%。模型提取的橡胶林在泰国中部、南部半岛,马来半岛的东部和南部地区,苏门答腊岛分布较为集中,而泰国北部、加里曼丹岛及其他岛屿橡胶林相对稀疏。提取的橡胶林分布信息与文献[23,24]中橡胶林分布及相关文字描述相吻合。
研究区域地处热带,全年云量大、大气中水汽含量高,严重影响光学遥感影像及产品的质量,很难获取特定时间段内高质量的大范围无云影像。为克服缺少光学影像问题,尝试利用雷达等主动遥感影像来研究土地分类[25,26]。雷达遥感对土壤水分变化较为敏感,用于大范围分类时,会因土壤水分差异而影响精度,常用于小范围、高精度提取。橡胶作为一种多年生常绿植物,在空间分布上相对稳定。本研究所采用的Landsat 7 3年大气层顶影像产品融合数据集和MOD13Q1 NDVI时间序列数据集,均为多时次遥感影像通过simple compose方法合成的数据集,反映的是地物光谱特征的中间值,具有一定的稳定性,在本分类中具有较好表现。另一个影响橡胶林分布提取精度的原因是农田、天然森林、橡胶林之间的混淆[26]。当森林上空悬浮的薄云未达到云识别阈值时,会造成NDVI值的下降,可能会误识为橡胶林。东南亚三国森林覆盖率高,在55%~75%以上,即使有小比例森林辨识为橡胶林,也会造成用户精度较大下降。在这3个国家中橡胶林占国土面积的比例较小,而且种植的集中程度低于中国的海南及西双版纳,这也给橡胶林分布的提取带来了困难。未来需要进一步实地调研获取第一手典型样方,结合地形、橡胶林的年龄等信息来构建更高精度的分类模型,进一步对东南亚橡胶主产区开展长势和产量遥感监测,以满足政策制定、贸易判断等决策需求。
责任编辑：谢龙莲