林地资源是森林资源的重要组成部分，是森林和野生动植物生存与发展的根基，是国家发展的基础性自然资源和战略性经济资源。但是不可持续地使用和管理方式造成地球上数量惊人的天然林地严重破坏、退化和消失，并已对人类生活构成威胁。

近年来利用遥感技术提取林地信息、了解林地变化、调查植被资源日益普遍。如缨帽变换的遥感影像融合^{[1, 2]}、特征融合的SVM(support vector machine)分类算法^[3]、新型植被指数^[4]、面向对象^{[5, 6]}、决策树^[7-11]、随机森林^[12]、利用相关指数和影像确定最佳波段^[13]等成为近年来遥感信息提取的重要方法，所用遥感数据主要包括高分辨率和中分辨率影像。许多学者深入研究土地利用变化，如采用一般土地利用动态度模型、地类的变化量或者土地利用转移矩阵等方法对不同地区的土地利用变化进行分析^{[14, 15]}，时间跨度多为10年左右，很少有进行长时间序列的研究。

菲律宾植被资源丰富，在20世纪90年代初期有570万公顷的森林面积，其中有大片热带雨林区域，在沿海地带有较多红树林。本文以菲律宾Landsat 5和Landsat 8为数据源，首先对TM(thematic mapper)和OLI(operational land imager)数据进行预处理，然后进行相关植被指数的计算和缨帽变换，根据不同植被在不同指数上的阈值不同进行提取，再结合最大似然监督分类法对研究区林地进行提取，最后根据Google Earth影像进行目视解译修订获得1990年和2017年林地信息，并对各种林地变化进行空间分析，为林地利用规划和森林资源可持续管理提供丰富的信息，影响政府对现有林地资源的评估，帮助决策者为菲律宾林地制定可持续的发展政策。

1 数据和方法 1.1 研究区数据

菲律宾位于亚洲东南部，植被资源丰富，约570万hm²的森林面积，属季风型热带雨林气候，高温多雨、湿度大、台风多，年平均气温约27 ℃，大部分地区年平均降水量在2 000~3 000 mm之间，海岸线长约18 533 km，在沿海地带有较多红树林。地形多以山地为主，占总面积3/4以上，可分吕宋岛群、米沙鄢岛群、棉兰老岛群三大岛群。

研究中用到Landsat 5影像和landsat 8影像数据均来自美国地质调查局(United States Geological Survey, USGS)官方网站(https://earthexplorer.usgs.gov/)，成像时间在1988-1992年和2016-2018年两个时期内，两个时间段所用影像覆盖全部菲律宾。由于遥感影像清晰度受云量影响较大，因此选取云量较低影像，以便地物能够更好显现，更有利于提取林地。30 m分辨率数字高程数据来自地理空间数据云(http://www.gscloud.cn/search)。辅助数据是用于裁剪菲律宾边界和岛群边界的国家级和区级矢量数据。

本研究选择TM影像和OLI影像的多光谱波段。对原始数据的预处理主要包括辐射定标，并将处理后的多光谱影像与DEM(digital elevation model)数据进行波段组合。

1.2 研究方法

通过对菲律宾两个时期的遥感数据进行预处理，在计算植被指数和缨帽变化分量的基础上，根据不同林地在不同指数上的阈值不同，并结合最大似然法进行林地信息提取，最后结合Google Earth实际影像目视解译修订得到最终林地信息结果。统计不同时期3个岛群的林地面积，分析菲律宾林地信息空间变化特征和变化原因。其具体技术路线如图 1所示。

1.2.1 植被指数法

归一化植被指数(normalized difference vegetation index, NDVI)是近红外NIR波段与可见红RED波段反射率数值之差和数值之和的比值，表示总体植被绿化率或光合活性植被，应用于检测植被生长状态、植被覆盖度和消除部分辐射误差，能反映出植物冠层的背景影响，如土壤、潮湿地面、雪、枯叶、粗糙度等，是区分植被与非植被的重要判别条件^[16]。归一化建筑指数(normalized difference building index, NDBI)是在仿归一化植被指数基础上提出的^[17]，单一的归一化建筑指数用来提取建筑面积效果不佳^[18]，本文通过结合其他指数来提取再种林。归一化湿度指数^[19](normalized difference moisture index, NDMI)广泛应用于植被液态水含量反演。各指数计算公式分别为：

$ {\rm{NDVI}} = \frac{{{\rho _{{\rm{NIR}}}} - {\rho _{{\rm{RED}}}}}}{{{\rho _{{\rm{NIR}}}} + {\rho _{{\rm{RED}}}}}} $

(1)

$ NDBI = \frac{{{\rho _{{\rm{MIR}}}} - {\rho _{{\rm{NIR}}}}}}{{{\rho _{{\rm{MIR}}}} + {\rho _{{\rm{NIR}}}}}} $

(2)

$ {\rm{NDMI}} = \frac{{{\rho _{{\rm{NIR}}}} - {\rho _{{\rm{MIR}}}}}}{{{\rho _{{\rm{NIR}}}} + {\rho _{{\rm{MIR}}}}}} $

(3)

式中，ρ_NIR表示近红外波段的反射率；ρ_RED表示红波段的反射率；ρ_MIR表示中红外波段的反射率。

缨帽变换是一种针对植被信息提取的影像增强方法, 其使坐标空间发生旋转, 旋转后的坐标轴指向与植物生长有密切关系的方向^[20]。缨帽变换后的分量与植物生长过程和土壤有关，因此缨帽变换主要用于特征提取。对预处理后的多波段数据用相应公式进行缨帽变换, 得到亮度指数Light、绿度指数Green以及湿度指数Wet 3个分量。

Landsat-5 TM数据亮度指数Light、绿度指数Green以及湿度指数Wet 3个分量计算式为：

$ {\rm{Light}} = {b_1} \times 0.303\;7 + {b_2} \times 0.279\;3 + {b_3} \times \;0.474\;3 + {b_4} \times 0.558\;5 +\\ {b_5} \times 0.508\;2 + \;{b_7} \times 0.186\;3 $

(4)

$ {\rm{Green}} = - {b_1} \times 0.284\;8 - {b_2} \times 0.243\;5 - {b_3} \times \;0.543\;6 + {b_4} \times 0.724\;3 +\\ {b_5} \times 0.084\;0 - \;{b_7} \times 0.18 $

(5)

$ {\rm{Wet}} = {b_1} \times 0.150\;9 + {b_2} \times 0.197\;3 + {b_3} \times \;0.327\;9 + {b_4} \times 0.340\;6 -\\ {b_5} \times 0.711\;2 - \;{b_7} \times 0.457\;2 $

(6)

对于Landsat-8 OLI数据，本文只用到亮度指数Light和湿度指数Wet两个分量，计算式为：

$ {\rm{Light}} = {b_2} \times 0.302\;9 + {b_3} \times 0.278\;6 + {b_4} \times \;0.473\;3 +\\ {b_5} \times 0.559\;9 + {b_6} \times 0.508\;0 + \;{b_7} \times 0.187\;2 $

(7)

$ {\rm{Wet}} = {b_2} \times 0.151\;1 + {b_3} \times 0.197\;3 + {b_4} \times \;0.328\;3 +\\ {b_5} \times 0.340\;7 - {b_6} \times 0.711\;7 - \;{b_7} \times 0.455\;9 $

(8)

式中，b₁、b₂、b₃、b₄、b₅、b₆、b₇对应数据中相应波段的辐射亮度值。

研究发现由于“同谱异物”和“同物异谱”现象, NDVI、Wet在人工林和次生林处差异很小，本文利用Light、Green扩大这种差异，使用指数LNDVI和GWet，计算公式为：

$ {\rm{LNDVI}} = {\rm{Light}} \times {\rm{NDVI}} $

(9)

$ {\rm{GWet}} = {\rm{Green}} \times {\rm{Wet}} $

(10)

式(9)、式(10)可以用于区分Landsat 5 TM影像人工林和次生林。

1.2.2 最大似然法

最大似然法^[21]是假设每一个波段的每一类统计都呈正态分布，计算给定像元属于某一训练样本的似然度，像元最终被归并到似然度最大的那一类当中，可以同时定量地考虑两个以上的波段和类别，可用性很好，而且不需要扩展的训练过程，是一种广泛应用的分类器。对于指数法提取之外的林地信息，利用ENVI5.3软件选择短波红外波段、近红外波段、绿波段为RGB组合，每种类别选择40个左右样本作为训练样本，然后选择最大似然分类其进行监督分类，经目视解译得到最终林地信息。

1.2.3 决策树法

决策树方法主要是决策树学习和决策树分类两个过程。决策树学习过程是通过对训练样本进行归纳学习(inductive learning), 生成以决策树形式表示的分类规则的机器学习(machine learning)过程^[22]。

由于影像受云层的影响，使得在分类的时候不得不分出云层这一类别，再通过决策树的方法并结合自然界发展的一般规律去除部分云层信息以增加分类精度，例如1990年的云层相应部分在2017年是天然林，就可以用2017年的天然林去覆盖1990年云层信息。

1.2.4 年林地覆盖变化百分比

Puyravaud^[23]在2003年提出的公式建立年林地覆盖变化百分比(q)可以描述林地变化的速度，数值越大表明在研究期内某种林地变化迅速，计算公式为：

$ q = \left( {1/\left( {{t_2} - {t_1}} \right)} \right) \times (\ln \left( {{A_2}/{A_1}} \right) \times 100 $

(11)

式中，t₁和t₂为测量年份；A₁和A₂为测量年份的林地覆盖面积。

2 菲律宾林地遥感分类 2.1 林地分类体系

为了确定林地的土地利用类型，本次研究将林地覆盖分为以下几类: ①人工林；②再种林；③红树林；④天然林；⑤次生林；⑥灌丛；⑦草地。人工林是采用人工播种、栽植或扦插等方法和技术措施营造培育而成的森林，包括棕榈、橡胶等经济作物；再种林包括人工林被砍伐后还未长成人工林的林地和已经被规划建造人工林的林地；天然林包括自然形成与人工促进天然更新或萌生所形成的森林，具有很重要的生态服务功能；次生林经人为采伐和破坏后，天然恢复起来的森林，生态稳定性和生态功能较差；灌丛是指矮树或丛林，是一种以散布的耐旱灌木为主的地理景观；草地是指人工草地和天然草地的总称，具有特有的生态系统，是一种可更新的自然资源。

2.2 分类过程

人为观察各类林地在不同指数处的值，通过多次提取实验，最终确定指数规则：对预处理后的TM影像利用相应指数依次提取林地建立的规则为：天然林(NDVI>0.3, Wet>10, DEM>100), 人工林(LNDVI>70, GWET>100, NDVI>0.4), 次生林(LNDVI>0.3, Wet>10, DEM < 100), 再种林(NDBI>0, NDVI>0.3)；对预处理后的OLI影像利用相应指数依次提取林地建立的规则为：天然林(NDVI>0.5, Wet>-0.07, Light < 0.35, DEM>100)，次生林(NDVI>0.7, Wet>0, DEM < 100), 红树林(NDMI>0.5), DEM < 50)。实验中所用指数阈值适用于热带雨林气候的东南亚地区，由于获取到的每幅影像时间不同，云量不同，在确认阈值时每个指数阈值可根据影像进行微调，阈值的微小差异会使得提取结果有很大差异性。利用IDL语言将已经提出的林地信息归零，对剩下的林地利用最大似然法进行监督分类，最终得到7类林地信息。

2.3 精度验证和评价

通过精度验证，分类者和使用者都可以了解分类的精度以及分类结果的可用性。本文在每幅遥感影像上面每种林地选择20个左右感兴趣区，利用ENVI5.3软件中Confusion Matrix Using Ground Truth ROIs工具，得到各林地生产者精度、用户精度(见表 1)、总精度和Kappa系数。

在提取过程中会遇见整幅影像阈值难以确定的情况，可沿着建设用地边缘分成两部分进行林地提取后再合成一幅林地信息图。分析Landsat影像中所选地物特征，通过对于光谱、形状等特征分析探究和多次实验，依据每种林地确定阈值范围来构造分类规则，探究出各林地提取的最佳结果。在热带雨林地区利用植被指数、最大似然法以及决策树方法对遥感影像林地分类提取，采用混淆矩阵评价分类结果，可得到Kappa系数是0.88和0.91，分类的精度也较高。

3 菲律宾林地空间变化特征分析 3.1 林地分布特征

经过统计，菲律宾1990年、2017年人工林的在林地中所占比重分别为2.7%、4.7%;再种林的比重分别为0.3%、0.1%；红树林的比重为1.2%、1.1%；次生林的比重为43.1%、49.8%；天然林的比重为44.9%、40.1%；灌丛的比重为0.7%、0.4%；草地的比重为7.0%、3.7%；可知在总林地中所占比例和占地面积都增加的除了人工林和次生林之外，其他林地在总林地中所占比例和面积在减小；同时也得出1990年和2017年3大岛群各林地所占比例，次生林和天然林在每个岛群都占很重的比例(见表 2)，可知菲律宾大部分林地属于次生林和天然林。

根据式(11)可计算菲律宾3个岛群年林地覆盖变化百分比如表 3所示。

由表 3可知：吕宋岛群变化率百分比差距不明显，其中再种林面积变化最多，减少了2.4%；米沙鄢岛群林地变化百分比差距较大，其中次再种林面积变化最多，增加了135.2%；棉兰老岛群变化百分比和差距不明显，其中草地面积变化最多，减少了2.8%。

3.2 林地变化分析

一直以来研究学者更加注重天然林和红树林的变化情况，通过转移矩阵计算表明1990-2017年期间菲律宾整体上红树林主要向次生林(36.7%)、人工林(5.9%)、非林地(15.4%)转化，天然林主要向次生林(15.5%)、非林地(4.2%)转化。其中在吕宋岛群红树林有37.8%转化成次生林，1.2%转化成草地，9.2%转化为非林地，少量向其他林地转化；天然林有11%转化为次生林，1.1%转化为草地，3.8%转化为非林地，少量向其他林地转化。在米沙鄢岛群红树林有46.8%转化成次生林，17.7%转化为非林地, 少量向其他林地转化；天然林有15.3%转化为次生林，2.8%转化为非林地，少量向其他林地转化。在棉兰老岛群红树林有34.7%转化成次生林，10.2%转化成人工林，17%转化为非林地，少量向其他林地转化；天然林有22.6%转化为次生林，1.3%转化为人工林，5.5%转化为非林地，少量向其他林地转化。

已有研究表明，森林砍伐是由各种因素造成的，包括经济、人口和政治因素^[24]。菲律宾是个人口增长率较高的发展中国家，1990年菲律宾总人口由6 195万人增长至2017年的10 492万人，年均增长率为2.6%。人口的增加势必增加对粮食和其他生活必需品的需求量，故而更多的林地被开发。

菲律宾在2011年依然是世界毁林最严重的十个国家之一，菲律宾缺乏防止毁林的长期战略和措施，政府部门能力不强，地方社区也没有能力借助有关REDD(减少毁林和森林退化排放)的国际讨论来有效地保护和恢复天然林以及保护生物多样性，导致菲律宾在恢复天然林方面成效甚微，故菲律宾政府在2011年开始修订国家林业政策以抑制森林退化和毁林。

4 结束语

本研究基于阈值的植被指数法，在菲律宾较大范围内实现了对中低分辨率的Landsat 5和Landsat 8数据进行了林地信息提取，分类精度分别是94.54%和95.69%，Kappa系数分别是0.88和0.91，分类精度较好。菲律宾林地主要为天然林和次生林，1990-2017年期间林地信息和空间分布发生了变化，人工林和次生林的面积在增加，而天然林减少面积最大，天然林向人工林和次生林转化；再种林、红树林和灌丛面积总量减少量稍稍明显，草地面积总量比例由7%缩减为3.7%，减少量明显，但是在米沙鄢岛群再种林呈现迅速增长趋势，年林地覆盖变化百分比高达135.2%，远远大于其他林地类型的年林地覆盖变化百分比；利用遥感影像分析林地变化有助于菲律宾对本国林地信息有准确的认识，为菲律宾地方政府对于林地资源的保护与利用提供帮助，但菲律宾境内的天然林修复应该引起政府高度重视。