-
天然橡胶是国防和经济发展不可或缺的战略物资[1]，在“金三角”贫穷国家发展替代罂粟的橡胶种植业[2-4]，也是我国促进边疆繁荣稳定的国家战略。地处广义“金三角”老挝北部5省(南塔省、丰沙里省、乌多姆塞省、波乔省和琅勃拉邦省)充足的水热条件为该区橡胶种植提供了良好的环境。随着橡胶林不断扩张，原生植被的破坏与不合理的橡胶种植将损害中老边境地区生物多样性与生态安全，因此，掌握橡胶林分布和时空扩张监测机制对橡胶林种植规划、橡胶产品的宏观调控和进出口贸易等都具有重要作用。
橡胶林提取早期多以中低分辨率遥感影像通过监督分类进行，分类精度可达80%左右[5-7]。近年来，以中高分辨率遥感影像为主要数据源的面向对象方法逐渐成熟，Ekadinata等[8]、刘少军等[9]、刘晓娜等等[10]用面向对象方法提取的橡胶林总精度在85%左右。决策树方法需要对比多时相数据分析橡胶林提取时间窗口，基于光谱特征和其他特征提高分类精度[11]，目前尚处于探索阶段。综合而言，目前难以形成高效、稳定、可靠的橡胶林提取算法，快速了解大尺度橡胶林分布面积及时空特征仍是瓶颈。面向对象方法适用高分辨率影像，在保证提取精度的前提下较难适用于大尺度范围，但可作为其他方法分类的验证；而基于物候特征、像元和光谱特征、纹理特征等的优化决策树分类方法受影像数据源限制较小，分类精度较高[12]，所以，在大尺度范围内具有更好的可操作性和普适性。研究基于MODIS影像，根据老挝北部橡胶林物候特征构建橡胶林的空间分布遥感提取算法，进行橡胶林时空扩张监测研究。掌握老挝北部五省橡胶种植状况，为老挝北部国土资源开发、生态环境监测、橡胶产业可持续发展及社会经济发展提供参考。
HTML
-
老挝北部5省位于19°～22°30′N，100°～103°30′E，总面积638.41万hm2(图 1)。处于亚洲大陆和南太平洋群岛之间的陆桥区域，是中南半岛唯一的内陆国家，主要地貌类型为山地和高原。全年平均气温约14℃，旱季为11月到次年4月，雨季为每年的5月至10月，平均年降水量为1 250～3 750 mm。适宜的气候与地理位置使之成为老挝的橡胶宜植区。
Figure 1. Location of study area
-
遥感数据包括8天合成的2015全年和2010年1至3月的MODIS地表反射率产品MOD09A1、MOD09Q1，用于决策树橡胶林提取方法的构建；2015年3月2日的高分一号宽幅(GF-1 WFV)、2016年1月16日高分辨率(GF-1 PMS)影像用来验证MODIS影像的提取效果；谷歌高清影像辅助决策地类样本的选择。此外，还包括来源于地理空间数据云的DEM和2次野外调查数据。
研究区涉及到MODIS影像图幅分别为h27v06，h27v07，采用MRT软件基于命令行批处理实现影像拼接裁剪和投影转换。GF-1 WFV、GF-1 PMS影像景序列号分别为961759，2001118，利用野外调查控制点进行正射校正，并在大气校正后融合GF-1 PMS影像生成2 m多光谱遥感影像。
-
老挝北部热带季风气候使得区域内形成丰富的植被覆盖地类。本研究在国际生物圈地圈计划(IGBP)分类体系基础上，结合研究区实际情况，将植被覆盖地类分为橡胶林地、有林地、灌木林地、旱地、水田，其中橡胶林地根据不同林龄的光谱差异，分为橡胶幼林(＜6 a)和橡胶成林(≥6 a)。研究借助野外调查数据和谷歌高清影像建立5个省各地类的遥感解译标志，共300个样本。不同植被覆盖地类特征参数的选择是构建决策树的关键，研究区橡胶林独特的落叶现象易于在特定时间区别其他地类，因此，研究基于物候特征、归一化指数、主成分分析、纹理信息、原始光谱特征等采用决策树分类方法，建立老挝北部5省橡胶幼林、橡胶成林决策树提取模型。不同植被覆盖地类具体光谱差异及橡胶林提取将在下述研究中阐述。
-
研究在2010年、2015年橡胶林空间分布叠加基础上，监测老挝北部5省橡胶林面积时空扩张过程、扩张的地形因子变化特点和研究区橡胶林扩张的驱动因素。橡胶林面积扩张分析采用单一土地利用类型相对变化率模型[13]，该模型可定量表达一类土地利用类型数量变化在不同区域间的差异性。表达式如下：
式中：Ka和Kb分别表示某一区域内部某种土地利用类型在研究初期和末期的总面积，Ca和Cb分别表示整个研究区该土地利用类型在研究初期和末期的总面积。研究结合地形因素的3个主要因子(海拔、坡度、坡向)构建2010年到2015年橡胶林时空动态变化分析模型，基于该模型掌握5年间老挝北部5省地区橡胶林时空变化特征，并从土地利用/土地覆被变化(LUCC)角度探索橡胶林时空变化驱动机制。
2.1. 数据源及预处理
2.2. 研究方法
2.2.1. 橡胶林提取方法
2.2.2. 橡胶林时空扩张分析方法
-
归一化植被指数(NDVI)可间接反映植被在不同季节、生育期的不同生理特征，多时相NDVI时间序列曲线能够表现植被在不同生长时期的差异性，以及同一时期不同植被在曲线上的区别，即植被的物候特征。研究通过IDL语言批处理生成6个地类样本的MODIS-NDVI时间序列，并基于时间序列谐波分析法(Hants)对时间序列曲线平滑去噪(图 2)。
Figure 2. Time series NDVI curves of different vegetation coverage land types by Hants
6类植被NDVI总体值在8月进入低谷，结合研究区8月份降水值可知：老挝北部8月份以来频频遭受洪涝灾害，使得大气引起的NDVI值总体偏低，目前各种滤波方法对NDVI时间序列的单个大气噪点平滑效果较好[14]，而对多期数据出现的大气噪点普遍效果不理想[15]。依据该区水田生长的时间特点，Hants滤波曲线除8月份外都客观反映了区域植被生长随时间变化的特点。滤波后的NDVI时间序列结果显示：水田和旱地NDVI时间序列曲线除去8月份的低谷噪点后为单峰单谷，水田NDVI曲线也表现为5、6月份上升，10月份后迅速下降，与实地调查中北部5省水田5月份种植，10月份收割的实际状况一致。有林地和橡胶成林的曲线最为接近，有林地NDVI值从11月至次年3月高于橡胶成林且为区域内NDVI值最大的植被地类，而在4月下旬至9月下旬的雨季间，其一直低于橡胶成林，表明水分条件是影响橡胶林生长的重要因子。橡胶幼林和灌木林地的NDVI曲线趋势非常相近，且有部分重叠，但在1月上旬至4月上旬，二者的NDVI值差异明显。综合橡胶林的落叶物候特征，12月至1月下旬橡胶林NDVI下降至最低后在2月上旬开始增长至4月上旬，表明此间橡胶林正经历着落叶及生长期，其NDVI值与其他植被区分明显。橡胶林物候特征和NDVI曲线地类差异表明：1月到3月是橡胶林提取的时间窗口。基于此，本研究经过对比选择8 d合成的MODIS-NDVI时间序列的第5、6、7、8、9、10波段数据，即DOY33(DOY表示年积日，从当前1月1日开始计天数，即DOY33表示第33天。)、DOY41、DOY49、DOY57、DOY65、DOY73六期影像平均值，并通过六期影像相对应的MOD09A1、MOD09Q1第5、6、7、8、9、10波段的数据计算其他指数平均值，分类中同时采用其他季节的影像辅助识别以提高分类精度。
-
研究在橡胶林提取的最佳时间窗口基础上采用所有植被类NDVI值大于0.35的特点，以0.35作为阈值提取了2015年植被区。将对裸地信息敏感NDBI作为高、低密度植被覆盖区的区分标准，图 3表明：有林地和橡胶成林的NDBI指数小于0，灌木林地、橡胶幼林、旱地和水田NDBI指数均大于0，即NDBI值大于0即为低密度覆盖植被，小于0为高密度覆盖植被。
Figure 3. NDBI of different land types
主成分分析(PCA)是多光谱数据的压缩和去噪工具[16]。6类植被第一主成分分量(PCA1)值(图 4)表明：橡胶成林、有林地PCA1值大于0，其余地类均小于0，且有林地在PCA1上的值远大于0.04，因此，PCA1既可以用来辅助分离高密度和低密度覆盖植被，同时也是区分橡胶成林与有林地的有效方法。
Figure 4. PCA1 of different land types
纹理特征是影像信息提取的重要方法之一。4种主要的纹理特征中，角二阶矩反映了图像灰度分布的均匀程度和纹理的粗细度，熵是与角二阶矩相关的一个测度值，对比度对于检测图像边缘特征和图像反差具有良好的效果，同质性可度量图像纹理的局部变化大小[17]。影像4种纹理特征(图 5)对比发现，有林地熵值和对比度值皆最小；不同低密度植被在角二阶矩和同质性上的特征相近，表明其在图像灰度分布上均匀度一致；对比度和熵上的变化特征明显，尤其是对比度值大于1时，水田分离效果最好。
Figure 5. Texture features of different land types
植被原始光谱(图 6)中第1、3、4、7波段上各地类平均光谱值较为接近，而第2波段(B2)上灌木林平均反射率最高且大于0.25，旱地在第5波段(B5)上的平均反射率最大且大于0.3，因此，第2波段和第5波段可作为唯一区分旱地和灌木林的光谱特征。基于上述研究中逐层提取的橡胶成林和橡胶幼林，研究最终基于橡胶林落叶期的关键物候特征得到老挝北部橡胶成林、橡胶幼林决策树提取方法(图 7)。
Figure 6. Original spectral characteristics of different land types
Figure 7. Decision tree classification modelbased on MOIDS
-
采用2次野外调查得到的实地橡胶林点分布数据，结合最新Google earth高清影像获取不同地类样本(ROI)。根据ROI生成的混淆矩阵评价老挝北部5省MODIS影像橡胶林分类总精度达到80.67%，Kappa系数为0.8，橡胶成林和橡胶幼林分类总精度分别达86%和78%(表 1)。
地类 参考样本数 分类数 正确分类数 生产者精度/% 用户精度/% 有林地 50 54 45 90 83.33 橡胶成林 50 49 43 86 87.86 灌木林 50 56 38 76 67.86 橡胶幼林 50 50 39 78 78.00 旱地 50 44 35 70 79.55 水田 50 47 42 84 89.36 总计 300 300 242 总精度=80.67%，Kappa=0.805 3 Table 1. Classification accuracy evaluation
-
进而通过GF-1 WFV影像采用面向对象方法提取老挝北部南塔省橡胶林，在影像分割尺度80的基础上，以NDVI、RVI、DEM和纹理特征建立隶属度函数进行模糊分类，提取老挝北部南塔省橡胶林，提取结果采用混淆矩阵验证得到的总精度达88.76%，kappa系数为0.84，橡胶成林的分类精度91.65%，橡胶幼林为81.28%。选取若干像元用GF-1 WFV影像对MODIS影像提取的2015年南塔省橡胶林分布进行对比验证，结果表明：基于MODIS数据的橡胶林空间分布提取精度为GF-1 WFV影像的91%，同时, 利用GF-1 PMS影像通过面向对象和人机交互的解译方法提取南塔市周边橡胶林，作为GF-1 WFV影像验证结果的补充，混淆矩阵计算的总体精度为96.76%，kappa系数为0.93。研究在GF-1 PMS高精度的分类提取基础上，分别从像元和总面积2个角度验证了MODIS影像提取橡胶林的精度：(1)当GF-1 PMS影像上确定为橡胶林分布的区域时，在对应的MODIS影像上随机抽取若干像元验证基于MODIS影像的橡胶林空间分布提取精度为82.49%。(2)GF-1 PMS影像覆盖范围内橡胶林分布面积为1.37万hm2，对应范围内的MODIS提取橡胶林面积为1.63万hm2，GF-1 PMS影像提取的橡胶林面积为MODIS影像的83.95%。像元和面积2个角度的验证结果均与混淆矩阵和基于GF-1 WFV影像验证结果近似，因此，对于大面积橡胶林提取，MODIS影像可以在保证时效性的同时保证提取精度。基于验证的橡胶林提取算法，研究分别提取了2010年和2015年的橡胶林空间分布(图 8)。
Figure 8. Distribution of rubber plantations in the five provinces of northern Laos based on MODIS
-
研究区橡胶林主要分布于老挝北部南塔省、丰沙里省西部和琅勃拉邦省北部，同时在波乔省和乌多姆赛省西部有较小范围的零星散落分布。2015年橡胶成林面积和橡胶幼林面积分别为14.17、2.79万hm2。2010年橡胶林成林、橡胶幼林面积分别为1.47、11.13万hm2。5年间橡胶林面积增加了4.36万hm2(图 9)，其中，橡胶成林面积从2010年的1.47万hm2增加到2015年的14.17万hm2，而橡胶幼林面积大幅减少了8.34 hm2，这表明老挝北部5省2010年的橡胶幼林在2015年大多转化为橡胶成林，而这5年间橡胶幼林面积新增较少表明橡胶种植热度减退。
Figure 9. Newly distributed rubber forest between 2010 and 2015
根据单一土地利用类型相对变化率模型得到老挝北部5省2010到2015年期间5省橡胶林土地利用的相对变化率(图 10)。5年间橡胶林土地利用在数量方面的区域差异性非常明显，丰沙里省扩张速度最慢，相对变化率只有58.88%；乌多姆赛省橡胶种植面积扩张速度最快，相对变化率达到328.85%，其次是波乔省，这说明乌多姆赛和波乔成为橡胶种植发展的新区域；南塔省和琅勃拉邦省橡胶林相对变化率分别为70.16%和74.30%，表明2省橡胶林持续缓慢增长。
Figure 10. Rubber forest land use relative change rate between 2010 and 2015
-
地形分析中原始坡向图只代表坡面的朝向，值域为[0, 360][18]，而由于坡向不同所造成的太阳辐射能量不同是影响橡胶林分布的主要原因，因此，研究对坡向转换后生成Trasp图[19]，Trasp计算公式如下：
式(2)中：Trasp值域为[0, 1]，0表示接收到的太阳辐射能最少，1表示接收的太阳辐射能最大，aspect为原始坡向图。
研究将2010年和2015年橡胶林分布图分别与海拔、坡度和Trasp图叠加，提取老挝北部5省橡胶林扩张的地形特征。海拔以300 m为单位，Trasp值以0.1为单位，坡度依据《水土保持综合治理规划通则GB-T15772-1995》坡度分级体系划分为6级[20]，分别为微坡(0~5°)、较缓坡(5°~8°)、缓坡(8°~15°)、较陡坡(15°~25°)、陡坡(25°~35°)和急陡坡(>35°)。统计单位海拔、坡度、转后的坡向范围内的栅格数占橡胶林分布区总栅格数的比例(表 2~4)，分析研究区2010年和2015年橡胶林分布的地形扩张特征。
年份year 海拔Altitude/m 0~300 300~600 600~900 900~1200 1200~1500 1500~2000 2010 0.00 20.09 59.63 16.62 2.79 0.87 2015 0.00 16.07 50.48 26.19 6.24 1.02 Table 2. The spatial distribution of rubber forests at different altitudes between 2010 and 2015
% 年份year 坡度Slope/(°) 0~5 5~8 8~15 15~25 25~35 >35 2010 13.58 13.59 34.89 32.21 5.45 0.28 2015 10.47 11.13 32.47 37.37 8.04 0.52 Table 3. The spatial distribution of rubber forests at different slope between 2010 and 2015
% 年份year Trasp 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 2010 9.13 8.81 7.27 6.74 6.39 6.41 6.85 8.53 19.18 20.69 2015 9.50 8.80 7.42 6.94 6.62 6.74 7.02 10.63 16.18 20.15 Table 4. The spatial distribution of rubber forests at different Trasp between 2010 and 2015
% 老挝北部5省橡胶林大多分布在海拔500~1 200 m范围内，1 500 m以上只有零星分布，2010年在此梯度分布的橡胶林占研究区橡胶林总面积的96%，2015年降为93%。2010年和2015年橡胶林均集中分布于海拔600~900 m；600~1 200 m范围是橡胶林分布的主要区间。2010年到2015年的分布变化表明，1 200 m海拔以上分布比例增加，600~900 m橡胶林分布没有明显变化，但在900~1 200 m范围橡胶林面积不断増加，分布特征在海拔上总体呈现出不断升高的趋势。
2010年到2015年间，橡胶林分布的主要坡度均在5°~25°范围内。大于35°的急陡坡和0~5°的微坡范围内少有分布。5年间橡胶林扩张的主要坡度在15°以上，2015年橡胶林面积在该坡度区间大幅增加，导致15°以下区域橡胶林面积分布比例相比2010年下降，表明橡胶林在15°坡度以上的分布逐年上升，而在15°以下区域基本保持原有水平，橡胶林在坡度上的扩张总体呈现种植区域坡度不断上升的趋势，但35°是橡胶林分布的上限坡度。
研究区橡胶林在上Trasp图主要集中在值为0.8到1.0的范围，即分布在太阳辐射的高值区，同时有大于50%的橡胶林分布在值为0.1~0.7的区间。2010年到2015年，橡胶林在各Trasp值上大多呈现增加趋势，但在Trasp值为0.9和1.0时2015年橡胶林分布比例下降，结合其在海拔、坡度上的变化，5年间橡胶林的扩张受到地域因素限制，坡向造成的太阳辐射强度在较高海拔地区被削弱。
-
2015年橡胶成林源于2010年橡胶幼林，2015年橡胶幼林面积大幅减少表明2010年后老挝北部5省橡胶林定植面积很少。橡胶林在面积和地形上的扩张特征显示其生长的立地条件不断变化。橡胶林扩张变化是LUCC的一部分，自然环境和人类活动的干扰[21]是LUCC变化的驱动因素，其中，人类活动是LUCC变化的主要驱动因素，核心是土地利用政策和管理主体的变化[22]，总结老挝北部5省橡胶林面积变化的核心包括社会政策和社会经济两方面。
(1) 为从源头上根除老挝北部地区罂粟对中国及国际社会的威胁，中国政府积极在老挝北部地区发展“罂粟替代种植”，在经济利益和工业发展的驱动下，橡胶林从传统种植区不断向高纬度地区扩展，因此，2005年之后，老挝北部橡胶种植势如破竹；但琅勃拉邦省、南塔省分别于2008年2月和11月宣布暂停审批橡胶种植项目，老挝政府首先认为橡胶规模种植所带来的经济和社会效益并不明显，需要认真评价橡胶种植的经济、社会效益；其次由于2008年2月份橡胶价格降幅超过30%，政府担心橡胶开割期依然低迷的橡胶市场，鼓励当地种植见效快的农作物；再者，老挝国家土地署表明，老挝劳动力严重不足，前期的橡胶种植对经济、社会、环境影响没有做充分的可研，由此老挝北部南塔省和琅勃拉邦省2015年橡胶幼林面积非常少。此外，2012年老挝政府暂停全国橡胶、桉树种植、矿产开采的土地许可，导致2012年之后老挝北部5省橡胶幼林少有增加，2015年橡胶幼林面积相比2010年大幅减少。
(2) 根据天然橡胶生产国组织发布的数据得到2009年至2016年平均天然标准橡胶价格(图 11)。从2009年开始，天然标准橡胶价格经过2年多的持续上升到2011年3月突破历史最高纪录，达到每吨32 000元·t-1，涨幅高达158%，此时部分胶农受价格利益驱动有少量橡胶定植，而后天然橡胶价格一路大幅跌落，老挝北部橡胶种植不再快速蔓延。因此，基于社会政策和社会经济因素，老挝北部5省2015年橡胶幼林面积仅有2.79万hm2，占2015年橡胶林总面积的16.46%。
Figure 11. Standard rubber price between 2009 and 2016
3.1. NDVI时间序列曲线
3.2. 橡胶林提取决策树模型构建
3.3. 精度验证
3.3.1. 混淆矩阵验证
3.3.2. 高分一号影像验证
3.4. 橡胶林时空扩张分析
3.4.1. 面积扩张分析
3.4.2. 地形扩张分析
3.4.3. 扩张驱动因素分析
-
精度在遥感信息提取中是重要的描述参数。研究采用250 m分辨率结合500 m分辨率的MODIS遥感影像提取橡胶林，精度达到80.67%，在同类影像中相对较高；而老挝北部统计数据显示，2010年和2015年橡胶林面积分别为10.61、14.68万hm2，较研究提取的2010年和2015年橡胶林面积12.6、16.96万hm2相比偏小，可能原因有两方面：一是MODIS影像属低分辨率遥感影像，一个像元的大小为250 m×250 m或500 m×500 m，因此，提取过程中容易将橡胶林林间路、林间地和较小的山体沟壑部分都误分成橡胶林；二是统计数据存在欠缺状况，政府部门的统计数据大多为大规模橡胶园，农民零散种植的面积较小加之上报数据存在偏差，导致往往存在漏统现象，因此研究中提取的橡胶林面积偏大。
研究从面积和地形两方面提取了老挝北部5省橡胶林空间分布及时空变化特征。下一步的研究将考虑更多影响橡胶林空间分布的自然因子，如地貌因素中的土壤因子，气候因素中的温度、降水、湿度因子等，从而分析橡胶林空间分布的多自然因子变化特征，进而定量分析自然因子对橡胶林空间分布的影响，得到橡胶林空间分布驱动模型。
随着高分辨率遥感影像的易获取性，橡胶林空间分布提取在决策树分类模型基础上，可融入专家知识，非遥感决策函数等建立更高精度区域性强普适性的遥感橡胶林提取模型，提高橡胶林分类精度，以满足橡胶林分布区土地资源管理和环境保护的需要。
-
研究采用遥感数据结合橡胶林物候特征实现了橡胶林空间分布快速可重复提取的新方法，并基于社会和经济数据建立了橡胶林时空扩张模型。在此基础上，探索了橡胶林时空扩张特征及其驱动机制。
(1) 1月至3月是橡胶林提取的时间窗口，结合NDVI、NDBI、PCA、影像的纹理特征和原始光谱值等构建的橡胶林提取决策树算法分类总精度较高，达到80.67%，Kappa系数为0.8。
(2) 基于高分一号GF-1WFV和GF-1PMS影像对MODIS影像提取的2015年南塔省橡胶林分布从像元和面积2个角度进行验证，验证结果在80.67%左右浮动，不同验证方法充分表明对大面积橡胶林提取MODIS数据可以在保证时效性的同时兼顾精度。
(3) 2010年到2015年橡胶林时空扩张特点表明，老挝北部5省橡胶林空间分布在海拔和坡度上总体呈现不断升高的趋势，海拔1 500 m以上只有零星分布，坡度分布上限为35°，坡向上分布无明显规律。从土地利用变化角度分析得到，老挝北部5省的社会政策和经济因素是影响橡胶林时空扩张的主要驱动因素。
DownLoad: