Landsat8时间序列影像支持下的农田防护林提取方法研究

雷思君，孙 华，刘 华，蒋馥根，陈 松，吴 童，谢 勇

(1.中南林业科技大学，湖南 长沙 410004；2.林业遥感大数据与生态安全湖南省重点实验室，湖南 长沙 410004；3.南方森林资源经营与监测国家林业与草原局重点实验室，湖南 长沙 410004；4.中国林业科学研究院资源信息研究所，北京 100091）

农田防护林是防护林体系中的一种重要类型，是一种由树木组成的农田廊道网络系统[1]。该系统能够改善农田小气候环境，调节农田生态系统的结构与功能，增加农作物产量，保护农田生态系统的生物多样性[2-3]，世界各国都将农田防护林体系摆在现代化农业的重要位置[4]。自1978年我国三北防护林体系工程实施以来，作为三北防护林体系中重要的组成部分，农田防护林为三北地区生态环境保护与粮食生产做出了重要贡献。因此，快速准确地获取农田防护林空间分布可以为农、林业的现代化管理和服务提供及时有效的信息支持，具有重要的意义。

目前，农田防护林信息的获取方法包括野外调查与遥感监测两种。野外调查方式周期长、成本高，不适用于信息化发展迅速的当代。遥感技术大范围、多尺度的特点为农田防护林信息的获取提供了便利[5]。国内外研究学者开展了大量基于遥感技术的防护林监测研究，其中基于中高分辨率影像的研究较多，尤其是Landsat 系列数据和ZY-3 影像数据[6-7]。遥感信息提取方法多依靠基于光谱信息的决策树方法[8]，其提取结果精度不高且斑块破碎，不利于对农田防护林进行更深层次的信息挖掘。一些学者为了提高农田防护林的提取精度，引入了农田防护林的形状因子及线性特征开展信息提取[9-10]，提升了在农田防护林信息提取研究中的适用性。但对于低分辨率的遥感影像，难以提取较短的林带，容易造成信息遗漏[11]。为了克服碎斑及林带断带等问题，近年来有学者研究了多尺度遥感影像在农田防护林信息提取中的应用，结合了光谱特征以及农田防护林的几何特征，综合利用数学形态学和面向对象的方法提取农田防护林[6-7]，由于特征信息有限，提取精度难以从根本上得到提高。与此同时，也有部分研究将纹理特征应用于农田防护林信息的提取中，在一定程度上改善了分类精度[12-14]。然而鉴于地物的多尺度特点，单尺度纹理特征难以满足不同地物的分类需求。

如今多时相遥感影像是提高遥感监测精度的重要信息源之一，被大量地应用于基于遥感的生产应用中[15-17]。同时，大量覆盖面积广、回访周期短的遥感卫星为多时相遥感影像的研究与应用积累了海量的数据源。在地物信息的提取研究中，多时相或时间序列遥感数据在农作物信息提取中应用比较广泛[18]，能达到识别不同作物的目的。多时相遥感数据是土地覆被遥感制图的重要信息源[19]，其提供的季相节律信息是地物分类的重要参数。物候学为多时相遥感数据的研究与应用提供了基本理论[20]。在提取植被信息的研究中，利用多时相遥感数据取得了一系列成果，NOAA/AVHRR、MODIS、Landsat等多种数据源的多时相数据在对森林植被进行分类的研究中，均取得较好的分类效果[21-23]。

本研究基于分层分类的思想，利用多时序遥感影像以及多分类指标开展农田防护林信息的提取工作。利用2017年4个季度9 期Landsat8 OLI 遥感数据，建立2017年甘州区地物NDVI时序曲线，通过曲线分析实现甘州区农田防护林的提取。

1 材料与方法

1.1 研究区概况

选择甘州区为研究区（图1），甘州区位于张掖市，处于甘肃省西北部，河西走廊中段的黑河流域，100°06′～100°52′E，38°32′～39°24′N。甘州区属于温带大陆性气候，境内大部分农作物于3、4月份播种，7、8月份收割。该区防护林建设属于三北防护林重点工程，防护林类型绝大部分属于农田防护林，以阔叶林为主，其中杨树面积最大[24]。研究区北部以及西南部为戈壁，主要植被类型为草本。

图1 研究区位置Fig.1 Location of the study area

1.2 数据来源与处理

1.2.1 遥感数据

选取2017年4个季度9 期Landsat8 OLI时间序列数据为农田防护林信息提取的数据源，并分别进行辐射定标与大气校正处理。由于云量的影响，缺少5、8、9 三个月份的遥感数据，研究使用6、7月份数据作为夏季部分的关键多时序数据，不影响多时序数据曲线的构建与关键物候参数的提取。

1.2.2 样地设计与调查

考虑样地大小与遥感影像像元的匹配问题，样地大小设计为30 m×30 m 并采用分层随机抽样的方法进行样地布设。为了保证抽样的可靠性与合理性，选用甘州区2017年7月GF-1 影像对研究区进行地类分层预抽样研究，设置抽样精度为95%，并且增加了20%的保险系数，在研究区内布设462个30 m×30 m 样地，并于2018年7月17日至8月25日进行野外样地调查。最终将 462个30 m×30 m 样地的实测值作为检验精度的验证样本。

1.3 研究方法

为了实现研究区内农田防护林的准确提取，研究提出一套基于光谱信息、形状指数、植被物候特征相结合的分层分类方法，具体研究思路如下：

1）在遥感影像预处理的基础上，首先利用特定物候期的一景影像提取水体指数，将研究区分为水体（包括研究区内的黑河流域、公园湿地、水库、灌溉水塘等）和非水体区域。

2）根据研究区现有地物状况，南北部无人工痕迹的戈壁滩在遥感影像上十分明显，利用特定物候期影像的NDVI 训练阈值，将明显植被和非植被区域进行区分。植被类型包括耕地、有林地和草地，非植被类型包括建筑用地、未利用地和部分草地（戈壁滩上的草本植被）。

3）植被类型和非植被类型均可利用Landsat8 OLI时间序列数据提取物候参数实现细分，为了提取农田防护林，非植被信息将不再作为研究对象。通过构建NDVI时间序列曲线，提取植被类型中的有林地、耕地和草地的关键物候信息，如植被生长季长度、生长变化幅度等，实现有林地、耕地和草地的区分提取。为了提取农田防护林信息，需将研究区内少数成片的公园林地等其他林地类型从有林地中分离，余下的则为整个研究区的农田防护林。

4）根据2018年样地野外调查结果及研究区已有土地利用类型图和高分辨率影像数据，对提取结果计算混淆矩阵，进行精度评价与分析。

2 结果与分析

2.1 水体信息提取结果与分析

为了选取提取水体信息的最佳特征指数，研究比较了NDWI和MNDWI 两个特征指数在水体信息提取中的差异和精度。利用野外样点调查数据和谷歌影像，选取67个水体样点和112个位于水体边界的非水体样点，对研究区水体提取结果计算基于样点个数的混淆矩阵。精度表明：MNDWI 指数在水体信息的提取中比NDWI 指数稍居优势。NDWI 在提取水体过程中，会混有土壤、阴影信息，改进的水体指数MNDWI 则能很好地区分阴影和水体。

2.2 植被与非植被信息提取结果与分析

经过水体信息的掩膜，对研究区内剩余地物提取NDVI 信息，利用该指数可以进行植被与非植被信息区分提取。2017年7月16日的遥感影像除去水体信息的其他地物NDVI 范围分布见图2。

图2 各地物NDVI 值分布范围Fig.2 Distribution of NDVI values for each surface

从图2可以看出：影像的NDVI 范围值贴合实际情况，7月16日处于植被最茂盛的时期，提取的耕地和林地NDVI 值均较高。从提取的NDVI范围结果可以看出：除去水体，其他地类具有明显的NDVI 值界线—耕地与林地处于最茂盛的时期，NDVI 最小值均大于0.2；建筑用地和未利用地的NDVI 最大值不超过0.2；草地的NDVI 值相对偏向于非植被类型，同时与植被（耕地、有林地）存在信息混淆现象。这是因为研究区内的草地大致分为两类——戈壁滩上的草本植被和人工种植的成片牧草，其中牧草在此季节处于茂盛时期，其NDVI值相对偏大，后期将此归类于耕地类型中。根据NDVI 阈值训练，研究利用归一化植被指数最大值NDVImax将剩下的5个地类分为植被与非植被。其中NDVImax＜0.2 定义为非植被，其余定义为植被。由此，最终将从植被（有林地、耕地）中提取农田防护林。

2.3 多时序数据支持下的物候特征提取及地物分类

2.3.1 NDVI时序曲线重构

利用NDVI时间序列数据提取地物的物候信息时，需对数据进行平滑处理，以减少噪声及缺失值的影响。常用的平滑方法有：Savitzky-Golay[25]、分段高斯函数拟合[26]等。本研究利用多项式函数对NDVI时间序列进行重构（见图3），耕地与有林地的拟合精度均达到85%以上。对比平滑前后的NDVI时间序列曲线，重构的光滑曲线保持了原有曲线的基本走向和形状，去除噪声后能更好地描述一年间地物随季节的微小变化。

图3 两种地物的NDVI时间序列曲线Fig.3 Time-series NDVI curves of four classes

从两种地物的NDVI时间序列曲线中可以看出：经过多项式函数拟合的曲线噪声明显减少，曲线更加平滑，并且保持了原曲线的形状特征。不同地物的NDVI时间序列曲线形状特征不同：耕地与有林地都具有明显的生长季，每年的5月到7月之间生长速度最快，并且在之后一段时间内保持较高的NDVI 值；其中耕地的生长季变化幅度较大，NDVI 值从0.2 以下变化至接近0.8，而有林地的生长季变化幅度不如耕地大，最大值在0.55 左右。每年生长季过后，耕地的NDVI 值会逐渐降低至接近0，而林地则保持有植被特征，NDVI 值始终不会变为0。

2.3.2 物候参数提取

利用平滑去噪后的NDVI 序列曲线可以提取地物的物候参数。本研究提取了5个物候参数，以a至e表示，如表1和图4（以耕地为例），其中：a、b为生长季始期和生长季末期，分别定义为NDVI在拟合曲线上增加、减少的速率明显提升或减慢的时刻；c为生长顶峰，定义为NDVI 在拟合曲线上的峰值；d为生长季长度，是生长季末期与生长季始期的差值；e是生长幅度，是生长顶峰与生长季始期（生长季末期）的差值。

图4 物候参数示意（以耕地为例）Fig.4 Seasonal parameters marked by cultivated land

表1 物候参数Table1 Seasonal parameters table

2.3.3 地物识别

对比有林地与耕地的时间序列曲线进行分析（图5）：耕地的曲线相对对称，耕地相比林地具有较明显的生长幅度（e），生长顶峰（c）比林地大，且耕地很容易确定生长季长度（d）；另一方面，耕地除了生长季便表现为非植被特征，NDVI 值低于0.2，这与实际情况吻合。定义规则（c＞0.55 a和e＞0.4）将耕地进行分离提取有林地。利用时间序列曲线分析结果能有效降低耕地和林地之间存在的同谱异物和同物异谱现象，减少对农田防护林信息提取的干扰。

2.4 甘州区农田防护林提取结果与分析

将甘州区的有林地成功提取后，需将其他林地进行分离得到农田防护林。经野外调查获取到：其他林地主要是甘州农场和甘州区林科所等区域的成片林地，农田防护林包含农田防护林带、林网和部分片状的林地。甘州区有林地提取效果见图6，利用462个野外样点调查结果与该地区土地利用类型图对分类结果进行精度评价。结果表明：甘州区有林地提取总体精度为85.93%，kappa 系数为0.79。

图5 耕地与有林地NDVI时序曲线对比Fig.5 Comparison of NDVI sequence curves between cultivated land and forest land

将提取的有林地结果进行人机交互去除其他有林地得到农田防护林信息。甘州区的农田防护林信息提取效果见图7。甘州区农田防护林提取结果呈现林带、林网以及部分片状林地的交错分布。利用野外样点调查数据进行精度验证：33个农田防护林样点中29个被正确提取，精度达到87.8%。

3 结论与讨论

3.1 结 论

1）NDVI时间序列数据能有效提取地物物候参数。利用生长季长度、生长幅度等物候参数能将有林地与耕地信息进行分离达到提取农田防护林的目的，这种基于植被物候信息的地物区分方法在对农作物进行细分的时候应用比较广泛。本研究借鉴有林地与农作物之间的物候差异性对甘州区农田防护林进行提取，地物之间的区分条件不如农作物之间进行细分的条件苛刻，并且分类精度较传统的依靠光谱特征差异的分类方法高。

2）多时序、多特征分层分类方法能有效提取农田防护林信息。结合野外调查的33个农田防护林样点，其中29个被正确提取，分类精度为87.8%。现有研究多注重利用光谱信息、形状指数等进行决策树分类，很少有研究将多时序与多特征方法相结合来提高地物的分类精度。

3.2 讨 论

1）本研究利用多时序遥感影像提取植被物候信息，这是非常关键的一步。虽然单期遥感影像可以较好地提取植被信息，但是不能很好地解决影像中一些地类存在的同谱异物和同物异谱现象。而论文提出的利用时间序列数据提取植被物候信息，可以获取同一位置不同时间阶段的影像信息，通过分析不同时间阶段影像的特征，确定其是否属于植被。通过这种方法可以有效地解决单期影像所不能完全解决的同谱异物和同物异谱现象，有效地提高了农田防护林提取的精度。

2）农田防护林相对于其他土地利用类型具有明显的线性、网状等形状特征，利用遥感技术获取农田防护林可以挖掘丰富的特征信息。在下一步研究中，可以在时间序列曲线中提取更丰富的农田防护林特征信息，也可以考虑温度因子在多时序影像中的应用，将不同物候期温度的差异应用到地物分类的研究中。