基于多源遥感数据橡胶林识别及其种植适宜区研究

蒋永泉,史正涛,高书鹏

(云南师范大学 地理学部,昆明 650500)

橡胶因其天然独特的物化特性而成为我国重要的战略基础物资之一[1]。橡胶树于20世纪初引入云南省种植[2],目前,已成为我国最大的天然橡胶生产基地。20世纪60年代云南省瑞丽市开始大面积种植橡胶林,截至2017年,全市共种植橡胶7 053.55hm2,总产值高达6 000万元余[3]。长久以来,橡胶林种植多以牺牲天然林为主,导致了当地物种多样性减少、景观连通性降低、水土流失等[4]诸多问题。因此,准确监测橡胶林的种植状况对于平衡经济发展和生态环境保护具有重要意义。目前,遥感技术在橡胶林识别中取得了长足的进展,MODIS数据是识别橡胶林使用最广的遥感数据源之一[5-6],但由于MODIS空间分辨率低,识别结果较差,很难将斑块较小的橡胶林识别出来。随着遥感技术的发展,多源遥感数据识别地物已成为遥感应用的趋势和热点。与单一数据源相比,多源数据可充分利用各种数据的高光谱、高时空分辨率构建各种特征,提高识别精度以满足应用需求。高书鹏等[7]基于时空数据融合算法获取的ETM+，OLI，Sentinel-2A时序数据提取西双版纳的橡胶林,Chen等[8]和Kou等[9]都基于雷达数据和Landsat数据分别提取海南省和西双版纳的橡胶林,结果显示,同时运用2种数据源均比单一数据源的提取精度高。因此,多源遥感数据可以弥补单一数据存在的缺陷,实现不同数据源之间的优势互补,将多源遥感数据进行融合,在一定程度上优化了分类结果,提高了识别精度[10],在橡胶林监测管理中具有较好的应用前景。

原种植于亚马逊河流域的橡胶树,要求高温、多雨、静风、光照充足的气候条件,而我国属于非传统植胶区,适宜种植区非常有限[11]。瑞丽市作为我国北缘橡胶分布区,最冷月平均气温仅有12.5℃,是世界上种植橡胶气温最低的地区,经常受到寒潮和寒害的影响[12]。1970—1971年、1975—1976年以及1982—1983年瑞丽市橡胶林都分别遭受了不同程度的寒害,使当地经济受到重创[13-14]。随着市场需求量扩大和胶价上涨,2004年后,胶农开始大面积砍伐森林种植橡胶林并向不适宜种植区扩张,导致了橡胶生长速度慢、产胶量低[15-16]。因此,掌握橡胶林的空间分布特征及其种植适宜区对规划管理具有重要的意义。本文综合利用Sentinel-2和MODIS-NDVI多源遥感数据识别橡胶林,该方法数据易于获取,可实现大范围的监测,基于识别结果,对瑞丽市橡胶林空间分布特征及其种植适宜区分析,可为相关部门决策和生态恢复提供参考。

1 研究区概况

研究区位于云南省西部瑞丽市,地处横断山脉高黎贡山余脉的向南延伸部分,地势西北高东南低,位于北纬23.38′～24.14′,东经97.31′～98.02′之间,其西北、西南、东南三面与缅甸接壤[17]。该区属于南亚热带季风性气候区,冬无严寒,夏无酷暑,全年分旱雨两季,年平均气温21℃,年降水量1 394.8mm,年平均日照2 330h[18]。瑞丽市总面积为1 020km2,其中森林面积占总面积的46.67%,主要森林类型有热带山地雨林,亚热带季风常绿阔叶林,亚热带山地落叶阔叶林,针叶林和竹林。高等植物类型丰富,共有高等植物284科,943属,3 733种,其中有国家级保护植物45种,省级保护植物24种[19]。全市植被覆盖度较高,其中橡胶林是当地种植面积最大的人工林[20]。

2 研究方法

2.1 数据准备及预处理

1) 选用的Sentinel-2遥感影像数据可通过欧空局(European Space Agency,ESA)网站(https://scihub.copernicus.eu/dhus/#/home)免费下载获取。大气校正通过Sen2Cor[21]小插件完成,其核心算法是大气辐射传输模型libRadtran[22],目的是将大气上层表观反射率转化为大气下层地表反射率。根据研究需要提取出10m的红光波段、绿光波段、蓝光波段和近红外波段,以及4个20m的植被红边波段并将其重采样为10m,然后裁剪出研究区大小的区域以备后续使用。

2) MOD13Q1是由美国国家航空航天局(NASA)提供的16d合成的250m空间分辨率NDVI产品,可在NASA网站(https://search.earthdata.nasa.gov/search)免费获取。获取的原始MODIS数据格式为HDF,投影类型与Sentinel-2数据不一致,因此本文采用MRT工具做投影转换,将其地理坐标和投影类型分别转为WGS84和UTM,然后裁剪出研究区大小的区域。为了得到10m×10m像元大小的NDVI数据,将23个波段的NDVI数据用立方卷积重采样处理,使其与Sentinel-2多光谱数据具有相同的像元大小。

3) 样本数据包括训练样本和验证样本,本文基于Google Earth高清影像,对研究区进行目视解译,将土地利用类型分为5类,分别为橡胶林、自然林、耕地、建设用地和水体。根据各地类的空间分布情况、复杂程度以及比例,一共随机选取了2 563个点,其中橡胶林527个、自然林764个、耕地681个、建设用地337、水体254个。并将各地类按1∶1的比例分开分别作为训练样本和验证样本。

2.2 地类的NDVI时序变化特征和光谱特征分析

为分析瑞丽市橡胶林与其他地类在不同时段的生长变化差异情况,选用2019年的MOD13Q1数据,分析5个地类在23个时相的NDVI响应关系(图1)。从总体上看,不同地类对NDVI的响应不同,橡胶林在整个时序变化过程中与其他地类存在明显的差异。第49d(2月)后,橡胶林开始萌芽进入返青期,NDVI值上升；第129 d(5月)—177d(6月)橡胶林的NDVI值急剧降低,这是由于夏季研究区内多云雨天气,云量较大NDVI值偏低引起的;在第241—289 d(9—10月)NDVI达到最大值;第305 d(11月)后,橡胶林的NDVI值开始缓慢下降,此时橡胶林开始落叶。

不同地类对不同波段的吸收和反射光谱存在一定的差异,图2为各地类光谱反射率均值曲线图。在蓝光波段、绿光波段和红光波段橡胶林易与水体混淆；在植被红光波段1,橡胶林的光谱值快速升高,且可与其他地类较好的分离；在植被红边2和3波段,橡胶林易与建设用地混淆；在植被红边波段4,橡胶林的光谱值达到最高,与其他有较好的可分性。因此,本文选用植被红边波段1和4用于区分橡胶林和其他地类。

图1 地类NDVI时间序列曲线图

注:“B2”蓝光波段；“B3”绿光波段；“B4”红光波段；“B5”植被红边波段1；“B6”植被红边波段2；“B7”植被红边波段3；“B8A”植被红边波段4

2.3 地物分类特征提取

基于Sentinel-2多光谱遥感影像提取分类特征,包括植被指数、纹理特征和光谱特征共23个特征,以及MOD13Q1时间序列数据的23个NDVI值(每16d一个),如表1所示。 归一化植被指数(NDVI)是反映植被生长状态的重要因子[23],与植被的覆盖度有密切的相关性。NDVI的计算公式如下[24]:

(1)

式中:ρNIR和ρred分别表示近红外波段和红波段的表光反射率。

灰度共生矩阵(GLCM)是20世纪70年代初由Haralick等[25]提出的,是通过灰度的空间特性来描述纹理的方法,能用来反映地物类别在空间特征的差异[26]。灰度共生矩阵共有14种特征指标来定量描述纹理,本文选用了最常用的5种,分别为均值(Mean)、同质性(Homogeneity)、相关性(Correlation)、熵(Entropy)和对比度(Contrast)。

(2)

(3)

(4)

(5)

(6)

式中:ux,uy,σx,σy,分别为行和列的均值和标准差；n为共生矩阵的阶数；i,j为共生矩阵的坐标；P为(i,j)处的共生矩阵数值。

表1 地物分类特征表

3 结果与分析

3.1 橡胶林识别结果

本文采用支持向量机分类器进行分类,分别对Sentinel-2多光谱遥感影像、MOD13Q1-NDVI时序数据以及多源遥感数据进行分类,结果如图3所示。从总体上看,3种不同数据源得到的瑞丽市橡胶林分布在空间上趋于一致。仅用单一的MOD13Q1-NDVI分类(图3(a)),不能较准确地识别出一些面积较小且分布零散的橡胶林。仅用单一Sentinel-2数据分类(图3(b)),虽然能将面积较小的橡胶林识别出来,但得到的结果比较破碎。图3(c)为综合运用多源数据的分类结果,由于综合了两种数据源在空间和时间上的优势,总的分类结果得到了很大的改善,能在准确识别出一些破碎地块的同时,又改善了Sentinel-2数据提取橡胶林结果的破碎程度,有效地避免了单一数据源的分类错误。

图3 橡胶林识别结果

本文基于混淆矩阵分别计算分类的制图精度(Producer Accuracy,PA)、用户精度(User Accuracy,UA)、总体精度(Overall Accuracy,OA)以及Kappa系数来评价结果的可行性。从表2可以看出,综合应用多源遥感数据能够实现瑞丽市橡胶林的高精度识别,优于仅用单一遥感数据源分类方法,联合多源遥感数据能够有效提高分类精度。综合运用多源数据,橡胶林的制图精度PA为89.73%、用户精度UA为73.98%、总体精度OA为91.41%、Kappa系数为0.76；与仅用MODIS数据相比,PA,UA,OA和Kappa系数分别提高了18.25%,11.52%,6.10%和0.19；与仅用Sentinel-2数据相比,PA,UA,OA和Kappa系数分别提高了16.73%,3.39%,3.21%和0.12。

表2 橡胶林识别精度评价对比表

3.2 橡胶林空间分布特征

从分类结果可以看出,瑞丽市橡胶林主要分布在的西南和东北方向上,对结果进行聚类分析处理去除破碎的错分对象,基于30m的DEM数据分析橡胶林空间分布特征。瑞丽市DEM数据的海拔高度范围为682～1 997m,将其划分为6个等级:682～900m,900～1 100m,1 100～1 200m,1 200～1 300m,1 300～1 400m,≥1 400m；根据我国水利部颁发《土壤侵蚀分类分级标准》[27]坡度范围,将瑞丽市坡度划分为6个等级:0～5°,5～8°,8～15°,15～25°,25～35°,≥35°。得到表3、表4的统计结果。

表3 瑞丽市橡胶林在海拔高度上的分布特征

表4 瑞丽市橡胶林在坡度上的分布特征

截至2019年底,瑞丽市橡胶总种植面积为7 693.83hm2,占土地总面积的7.54%,其中有6 433.11hm2种植在海拔为682～900m高度范围内,其次是种植在900～1 000m范围内,分别占橡胶林总种植面积的83.61%和15.26%；主要种植的坡度为5～15°,占橡胶林总种植面积的83.61%。

3.3 橡胶种植适宜区

中华人民共和国农业部发布的《橡胶树种植基地建设标准》(NY/T221-2016)[28]根据橡胶树对温度、水分、风速和光照等胶树生长要求以及产胶潜力,将云南省西部海拔≥900m高度以上的地区划分为不宜作为常规胶园宜林地。依据此标准的海拔高度变化将瑞丽市划分为橡胶适宜种植和不适宜种植2个区域,然后对2019年瑞丽市植橡胶林适宜区进行统计分析(图4),瑞丽市有6 433.11hm2橡胶种植在适宜种植区,占橡胶林总种植面积的83.61%,已有1 260.72hm2橡胶林向不适宜种植区扩张,占橡胶林总种植面积的16.39%。

图4 瑞丽市橡胶林适宜种植区

4 结论

1) 综合Sentinel-2多光谱影像和MOD13Q1 NDVI时间序列数据的分类方法,可充分发挥2种遥感数据源的不同特征优势,与仅用单一数据的结果相比,有效提高了橡胶林的识别精度,实现了瑞丽市橡胶林高精度识别,综合运用多源遥数据橡胶林识别结果的制图精度PA为89.73%、用户精度UA为73.98%、总体精度OA为91.41%,Kappa系数为0.76。

2) 截至2019年底,瑞丽市橡胶种植面积为7 693.83hm2,有83.61%的橡胶林种植在海拔为684～900m高度范围内,占比最大；其次是种植在900～1 100m高度范围内,占橡胶林总种植面积大的15.26%。主要种植的坡度为5～15°,占橡胶林总种植面积的83.61%。

3) 根据农业部发布的《橡胶树种植基地建设标准》将瑞丽市的土地按海拔高度划分为适宜种植区和不适宜种植区,截至2019年底,瑞丽市有6 433.11hm2橡胶种植在适宜区,占橡胶林总面积的83.61%,已有1 260.72hm2橡胶林向不适宜种植区扩张,占橡胶林总面积的16.39%。