应用GF-5高光谱遥感影像提取山区茶园

钱 瑞, 徐伟恒, 黄邵东, 王雷光, 3, 4, 鲁 宁, 欧光龙

1. 西南林业大学林学院, 云南 昆明 650233

2. 西南林业大学大数据与智能工程学院, 云南 昆明 650233

3. 西南林业大学大数据与智能工程研究院, 云南 昆明 650233

4. 西南林业大学生态大数据国家林业与草原局重点实验室, 云南 昆明 650233

引 言

近年来, 随着茶园不断扩张, 导致森林面积减少。 造成水土流失、 生物多样性减少及森林破碎化严重等负面影响[1-2]。 因此, 准确掌握茶园空间分布情况极其重要。 遥感技术具有宏观、 实时、 快速、 监测范围广以及持续动态监测等优势, 已经广泛用于茶园监测[3]。 如李龙伟等[1]利用Sentinel-2影像的红边与短波红外波段构建归一化茶园指数(normalized difference tea index, NDTI)和决策树(decision tree, DT)提取茶园, 其制图精度(producer’s accuracy, PA)达到92.28%, 用户精度(user’s accuracy, UA)为95.67%。 Xu等[4]利用支持向量机(support vector machine, SVM)结合地形、 植被指数及雷达数据提取了西双版纳勐海县的茶园, 其PA和UA分别达到了90.83%和91.67%。 黄邵东等[5]基于高分一号(GF-1)结合纹理和空间特征采用面向对象(object-oriented, OO)分类方法进行茶园提取, 基于OO的RF茶园提取的PA为70.54%、 UA为87.13%。 熊皓丽等[6]基于GEE平台使用Sentinel-1、 Sentinel-2结合光谱特征、 地形特征、 纹理特征利用支持向量机递归特征消除(support vector machine-recursive feature elimination, SVM-RFE)算法筛选特征, 结果表明利用SVM-RFE可以筛选出有利于茶园提取的特征集, 其茶园提取PA为91.64%、 UA为92.91%。 综上所述, 以往利用遥感技术对茶园提取的数据源主要基于中高空间分辨率遥感影像。 相较于Landsat-8、 Sentinel-2等光学数据以及Sentinel-1、 PALSAR等雷达数据, 高光谱影像具有空间分辨率高、 谱段数目多、 光谱连续、 信息量大等特点, 极大的增强了遥感对地物的鉴别能力和观斥能力[3, 7-8]。 因此, 以负载丰富高光谱信息的GF-5遥感影像更容易捕捉地物类别之间细微的光谱差异。

GF-5号在提取水体、 土壤、 植被方面也有较为广泛的应用[9-11]。 Chen等[12]采用GF-5号影像结合Sentinel-1A数据基于LPPSubFus的特征级图像融合方法, 提取光谱指数和纹理特征等因子进行茶园提取, 实验证明LPPSubFus整合高光谱和SRA数据用于茶园制图具有一定潜力。 GF-5影像在光谱细节的检测能力已经在地物类型的分类和提取中有所表现。 如任广波等[10]基于GF-5影像在黄河三角洲滨海湿地开展互花米草的分类, 发现在不同数据中SVM均能到达较好的分类效果, 总体精度(overall accuracy, OA)达到94.23%, Kappa系数0.93。 此外, Zheng等[13]基于GF-5高光谱数据对东北地区森林分类, 利用卷积神经网络(convolutional neural network, CNN)结合SVM构建基于CNN-SVM融合的高光谱森林分类模型, 其森林分类OA为88.67%, Kappa系数0.84。 然而, 目前基于GF-5号影像利用随机森林(random forest, RF)对茶园高精度提取研究还处于探索阶段。

RF具有运算速率快、 分类精度高、 稳定性强、 处理多维数据能力强等特点, 被广泛用于地表覆盖分类及植被识别等领域。 如宗影等[14]基于OO结合RF进行滨海湿地植被分类研究, 结果表明特征选择(feature selection, FS)后结合RF算法OA达到87.07%, Kappa系数0.84。 吕杰等[15]利用RF结合纹理特征识别森林类型, 在SVM、 最小距离分类(minimum distance classification, MDC)和RF三种分类器中, RF算法具有最好分类效果(OA: 81.01%, Kappa: 0.76)。 杨迎港等[16]基于GF-2影像结合RF算法进行土地利用类型分类, 结果表明, 基于RF算法的OA达到93.44%, Kappa系数0.93。 张磊等[17]利用Sentinel-2数据及RF进行黄河三角洲的植被提取, 其OA到达90.93%, Kappa系数0.90。 可见, 相较于其他分类算法, RF在土地利用类型识别领域应用较为广泛且精度较高。

本研究利用单景GF-5号影像进行茶园提取, 利用原始波段作为光谱特征(spectral feature, SF), 参考前人研究构建适用于GF-5号影像山区茶园提取的植被指数特征(vegetation index feature, VIF)以及结合数字高程模型(digital elevation model, DEM)的地形特征(TRF)进行特征组合形成12种分类方案, 并对比不同分类方案下茶园的提取精度, 最终确定最优分类方案, 为基于GF-5 AHSI卫星数据的山区茶园分布制图和资源监测提供技术支撑。

1 研究区及数据源介绍

1.1 研究区概况

如图1所示, 研究区位于云南省普洱市南部, 分布着大面积山区茶园, 面积约为36 hm2。 研究区四季差异不明显, 气候温和, 属于亚热带高原性季风气候。 年平均气温为18 ℃, 年降水量为1 414～1 800 mm, 海拔为753～2 283 m。 该研究区日照充足, 雨量充沛, 土壤以赤红壤为主, 土层深厚, 较适宜茶叶生长[18]。

1.2 遥感数据介绍及处理

1.2.1 GF-5数据源介绍

为了更好的区分茶园与其他地物, 本文选择单景GF-5号影像作为数据源。 已有学者提出茶园具有人为修剪的物候期, 而本文选择GF-5遥感数据时间正好是修剪后处于萌发阶段, 这个阶段茶园与其他地物存在一定的区分性, 因此选择成像时间为2020年3月11日的GF-5遥感影像作为本研究的数据源, 数据来源于自然资源部国土卫星遥感应用中心, 其详细参数见表1。 影像包括可见光及近红外波段150个, 短波红外波段180个, 空间分辨率为30 m。 利用ENVI 5.3 SP1软件对GF-5号影像进行辐射定标、 大气校正及正射校正等预处理。 此外, 从GEE平台合成并下载一景经预处理后覆盖研究区无云Landsat-8影像, 以Landsat-8影像为基准在ArcGIS上对GF-5影像进行绝对配准。 由于传感器不同探元间存在响应差异, GF-5影像存在水汽吸收波段和受到噪音影响波段[19-20], 如390～398、 933～954、 1 122～1 156、 1 333～1 502、 1 797～1 965、 1 999～2 015、 2 369～2 513 nm波长范围; 对于可见光及近红外波段(1 006～1 028 nm)和短波红外波段(1 004～1 029 nm)之间的光谱重叠, 为保留准确的光谱信息, 去除短波红外波长范围, 保留可见光波长范围, 移除波段见表2[19]。

表1 GF-5影像信息表

表2 GF-5卫星高光谱数据移除波段

1.2.2 DEM数据介绍

前人研究表明, DEM作为辅助数据可提高茶园提取精度。 熊皓丽等[6]研究表明地形特征对提取茶园有较好的辅助作用。 Xu等[4]利用光谱特征及海拔因子(900～2 400 m)对云南三大茶叶主产区的茶园分布进行提取的效果较好。 而本研究区正好位于云南三大茶叶主产区内, 因此基于GEE平台调用了空间分辨率为30 m的DEM数据(SRTMGL1: Global 1 arc second Voo3)[21], 并将其海拔、 坡度、 地形等因子作为特征因子, 参与对茶园提取和分布制图。

1.3 野外调查数据

实地调查了研究区内的主要地物类型并拍照、 定位, 调查路线详见图1。 根据实地调查的情况, 在GEE平台上对5类地物(茶园、 森林、 农田、 不透水层和水体)进行自视解译, 并勾绘了茶园样本1 130个、 森林样本690个、 农田样本735个、 不透水层样本330个以及水体样本194个, 共3 079个样本。 然后随机划分各类样本为训练样本和验证样本, 分别用于分类器的训练和精度验证[22], 为保证训练和测试样本数量的均衡性, 随机划分各类地物的30%作为训练样本, 70%作为验证样本。

2 实验部分

利用GF-5影像的SF、 VIF及TRF, 结合RF分类器进行茶园提取。 首先, 绘制茶园、 森林、 农田在原始波段上的光谱反射率曲线图, 基于三类地物反射率区分明显且为反射峰和吸收谷的波段构建新的光谱指数, 结合原始波段特征和地形特征构建初始分类特征空间, 应用RF对6种初始分类特征方案进行特征重要性排序, 将特征重要性按从大到小逐个依次输入RF分类器进行茶园提取, 随着特征的不断输入茶园F1-Score达到稳定不再明显上升时的特征维度则为最终分类特征空间, 其PA和UA即为最优分类特征下的茶园提取精度。 然后将特征优化前的6种初始分类方案特征输入RF进行茶园提取, 并通过混淆矩阵对RF算法下12种分类方案进行茶园提取精度评价。

2.1 特征优选

2.1.1 光谱分析

结合野外调查和基于GEE平台勾绘的各地类样本绘制典型地类的光谱反射率曲线, 分析各种地物的光谱差异, 为地物类型的可分性波段选择提供依据[23-24]。 统计了茶园、 森林、 农田的波段均值并绘制出三种地物类型的光谱特征曲线图, 如图2所示。 茶园与农田、 森林的光谱值在400～700、 1 974～1 990、 2 024～2 361 nm波长范围区分比较明显; 在800～1 029、 1 038～1 325、 1 511～1 788 nm波长范围茶园与森林、 农田之间的光谱值可分度较明显。 因此, 选择上述可分波段范围内的反射峰及吸收谷波段结合斜率及归一化函数构建茶园指数(tea plantation indices, TPIs)。

图2 三种地物的光谱反射率

所选波谷波长值分别为550、 668和2 083 nm, 波峰波长值分别为1 080、 1 266、 1 249和2 327 nm, 如图2所示。 基于峰谷之间的斜率构建的指数为TPI1—TPI3, TPI4则是根据505、 1 080和2 083 nm所对应的光谱反射率构建的三角形面积指数[12]。 TPI5—TPI9指数则是仿照NDVI指数[25]基于植被反射和植被吸收波段进行归一化构建。 植被指数特征构建公式见表3。 除了基于峰、 谷构建的9维植被指数外, 另外参考前人研究构建了36维植被指数, 详见文献[26], 因此, 共构建了45维VIF特征, 见表4。

表3 植被指数特征构建表

表4 特征因子表

2.1.2 特征集构建

为了比较不同特征集茶园提取精度, 本研究构建了三种特征集, SF 250维、 VIF 45维以及TRF 3维共298维初始分类特征, 如表4所示。

2.1.3 随机森林分类

RF是Leo Breiman提出的具有迭代特性, 以决策树为基础分类器的一种集成学习方法[27]。 构建RF的基本过程为: ①通过boot-strap方式从训练样本中有放回的随机的抽取样本, 假定训练样本集中共有M个样本, 每个样本都具有N个特征, 每次都有放回的抽M个样本, 那其中某个未被抽中样本的概率则为(1-1/M)M, 当M很大时, 这个值趋近于1/e≈1/3, 即表征为抽取时大约1/3的样本未被抽到, 这部分样本就被称为带外数据(out of bag, OBB); ②然后对M个样本进行训练得到决策树模型, 在决策树的各个节点处随机选择n(n

2.1.4 特征优选方法

针对高维特征空间, 为了评估各个特征对茶园提取的重要性, 优化特征组合, 最终选择RF算法进行特征重要性评价[22]。 基于GEE平台利用RF算法结合训练样本计算各个特征对茶园提取的重要性。 以RF算法作为分类器, 按照特征重要性由高到低依次增加特征, 观察茶园提取F1-Score的变化, 随着特征维数的增加F1-Score上升不明显甚至开始下降时所参与分类的特征维度为最佳分类特征空间, 此时的茶园提取PA和UA即为FS后各方案茶园提取最优精度。

2.1.5 分类方案

为了分析光谱、 植被指数及地形特征在茶园识别中的效果, 按照SF、 VIF及TRF构建12种分类方案, 见表5。 方案1—6为初始方案, 方案7—12为特征优化后方案, 其中单一特征因子分类为方案1, 2, 7, 8, 多特征组合分别为方案3, 4, 5, 6, 9, 10, 11, 12, 由于地形特征仅有三维, 因此不单独列为分类方案。

表5 分类方案

3 结果与讨论

3.1 特征优选分析

针对方案1—6共6种方案进行FS后分别得到方案7—12(见图3), 可见不同特征组合进行特征优选直到F1-Score达到稳定且不再明显上升的特征维度存在差异其中维度最低的是方案11(SF+TRF+FS), 维度为18维; 维度最高的是方案10(SF+VIF+FS), 维度为32维。 但不论何种优选方案, F1-Score均较高且差距不大, 在0.84～0.87之间。

图3 分类特征优选及特征重要性排序

3.2 不同分类方案茶园提取精度对比分析

12种方案茶园提取精度见表6, 可见不同方案茶园提取精度差异不明显, PA浮动在0.87～0.91之间, UA浮动在0.79～0.84之间, 即不论是SF、 VIF以及TRF三种特征单独分类还是组合分类, 有较少的茶园样本错分到其他地类中, 而其他地类样本则会较多被误分为茶园。 针对茶园提取精度, 不同方案在PA和UA上表现不一致, 故此, 引进F1-Score调和PA和UA, 以确定12种方案中针对茶园提取的最佳方案。 从F1角度看, 方案9(VIF+TRF+FS)和方案12(VIF+SF+TRF+FS)茶园提取精度最高, PA、 UA、 F1-Score分别为(90.41%, 83.43%, 0.87; 90.86%, 84.09%, 0.87), 两种方案中TRF全部被选中, 因为茶园生长对海拔等地形因子具有一定要求, 云南省茶园分布主要为900～2 400 m之间[4]。 另外, 在SF参与的4种分类方案(方案7, 10, 11, 12)中, b150被选中3次, b4, b5, b6, b27, b133, b281均被选中2次, 上述选中波段均位于茶园明显区别与森林与农田的光谱范畴内; 在VIF(方案8, 9, 10, 12)参与的4种分类方案中, REP、 VOG2、 SR2、 SR3、 WBI、 TPI3、 TPI9均被选中4次, 证明红边相关指数在对于茶园同其他地类区分时起到重要作用。 此外, 相较于6种初始分类方案, 经过FS之后的6种分类方案茶园提取精度, F1-Score均高于对应初始分类方案的F1-Score, 可见对高维特征进行FS, 既可避免“维数灾难”, 减少计算量的同时又可提高茶园的识别精度。

表6 基于随机森林不同方案茶园提取精度表

3.3 分类结果

选取分类精度最高的分类结果进行展示。 基于RF算法结合方案12(SF+VIF+TRF+FS)的分类结果见图4(a)。 由于不透水层及水体易于识别, 因此本文特征筛选时并未筛选不透水层和水体区别于茶园的光谱特征, 但是识别结果显示, 茶园、 不透水层和水体都能较好的区分开。 而茶园主要集中于研究区的西南部, 主要分布于倚象镇营盘山一带, 这与实际调查情况相符。 此外, 选取三块茶园连片分布区域放大进行识别结果展示, 图4(a)中区域1、 2、 3的识别细节分别见图4(b—d); 以基于91卫图下载谷歌影像为底图, ArcGIS将茶园识别结果叠加在底图之上, 并将茶园识别结果图层以绿色表示并设置60%透明度, 可见茶园识别效果显著。

图4 随机森林茶园提取分类结果(光谱特征+植被指数特征+地形特征+特征优选)

4 结 论

利用GF-5影像光谱信息丰富、 数据连续性强的特征, 提出一种适合于GF-5号影像山区茶园的提取方法。 仅用单景GF-5影像, 利用茶园同其他2类地物光谱差异明显且具有反射峰和吸收谷波段结合归一化函数及峰谷斜率构建相关光谱指数。 联合GF-5原始波段、 光谱指数及数字高程模型(DEM)等信息构建复杂山区地形下的茶园提取方法。 利用随机森林(RF)进行各类特征及特征组合的重要性排序, 按照特征重要性从大到小依次将特征输入RF分类器进行茶园提取, 随着特征的不断输入茶园的F1-Score达到饱和不再明显增加时的特征维度即为最优分类特征。 最终确定12种分类方案, 再结合RF分类算法进行茶园提取, 12种方案茶园提取的PA和UA分别高于87%和79%。 FS后6种分类方案中, SF参与分类的4种方案中被选中2次的波段为b4、 b5、 b6、 b27、 b133、 b150和b281; 在VIF参与分类的4种方案中被选中4次的指数分别为REP、 VOG2、 SR2、 SR3、 WBI、 TIP3和TIP9; TRF参与分类的4种方案中坡度、 坡向、 高程均被选中。 证明了SF、 VIF、 TRF以及三种特征的组合结合RF能够满足复杂地形下山区茶园的高精度提取, 且高维特征在经过FS后降低了运算量的同时进一步提高茶园提取的精度。

尽管本研究能够实现复杂地形下山区茶园的高精度提取, 但是仍然存在一些不足之处, 首先所使用的分类特征仅考虑光谱特征, 而忽略纹理等因素, 根据徐伟恒等研究, 茶园具有明显的“人为修剪”物候[28], 本研究在影像选择上尚未考虑这一点。 为提高茶园提取精度, 后期尝试结合光谱+纹理+地形特征进行茶园提取; 此外, 结合茶园自然生长状态的萌芽期、 茂盛期、 修剪期等物候生长差异, 选择最佳时相的遥感影像进行茶园提取具有很大潜力。 此实验中茶园提取仅使用单景高光谱影像, 未来将考虑使用高光谱结合高空间分辨率时间序列影像完成大范围的茶园提取。 但是如何克服传感器差异将异源数据进行融合, 结合多类型特征, 选择合适的分类方法, 将其应用于大范围山区茶园的提取有待进一步探索。