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

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 研究区示意图Fig.1 Schematic diagram of study area

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

Table 1

表1(Table 1)

表1 GF-5影像信息表 Table 1 GF-5 Image information table 影像 产品号 中心经纬度 成像时间 幅宽/km 空间分辨率/m 光谱分辨率/nm GF-5 AHSI 74 538 101.18° E, 22.89° N 2020-03-11 60 30 5/10

表1 GF-5影像信息表 Table 1 GF-5 Image information table

表2

Table 2

表2(Table 2)

表2 GF-5卫星高光谱数据移除波段 Table 2 Excluded bands of hyperspectral data from GF-5 satellite 数据质量 去除波段 波段个数 未校正 SWIR: 43~50, 96~112; VNIR: 1~3 25 条带污染 SWIR: 40~42, 51~60, 95, 113~115, 119~121, 163~180 41 水汽吸收 VNIR: 128~132; SWIR: 15~19, 43~57, 96~115 45 SWIR中的光谱重叠 SWIR: 1~4 4

表2 GF-5卫星高光谱数据移除波段 Table 2 Excluded bands of hyperspectral data from GF-5 satellite

1.2.2 DEM数据介绍

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

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

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

Fig.2

	Figure Option ViewDownloadNew Window
	图2 三种地物的光谱反射率Fig.2 Spectral reflectance curves of three features

所选波谷波长值分别为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

Table 3

表3(Table 3)

表3 植被指数特征构建表 Table 3 Construction of vegetation index feature 指数 公式 解释 TPI1 (ρ 1 080 nm-ρ 505 nm)/Δ λ 1 ρ 为对应波长光谱反射率值, Δ λ 1为505~1 080 nm的波长差, Δ λ 2为1 249~2 083 nm的波长差, Δ λ 3为2 083~2 327 nm的波长差 TPI2 (ρ 2 083 nm-ρ 1 249 nm)/Δ λ 2 TPI3 (ρ 2 327 nm-ρ 2 083 nm)/Δ λ 3 UPI4 (ρ 1 080 nm-ρ 505 nm)Δ ω 3-{((ρ 1 080 nm-ρ 505 nm)Δ ω 1+ (ρ 1 080 nm-ρ 2 083 nm)Δ ω 2+(ρ 2 083 nm-ρ 505 nm)Δ ω 3)}/2 Δ ω 1为505与1 080 nm的波长差, Δ ω 2为1 080与2 083 nm的波长差, Δ ω 3为505与2 083 nm的波长差 TPI5 (ρ 1 080 nm-ρ 505 nm)/(ρ 1 080 nm+ρ 505 nm) ρ 为对应波长光谱反射率值 TPI6 ρ 1 080 nm-ρ 668 nm)/(ρ 1 080 nm+ρ 668 nm) TPI7 (ρ 1 266 nm-ρ 505 nm)/(ρ 1 266 nm+ρ 505 nm) TPI8 (ρ 1 266 nm-ρ 668 nm)/(ρ 1 266 nm+ρ 668 nm) TPI9 (ρ 2 327 nm-ρ 2 083 nm)/(ρ 2 327 nm+ρ 2 083 nm)

表3 植被指数特征构建表 Table 3 Construction of vegetation index feature

表4

Table 4

表4(Table 4)

表4 特征因子表 Table 4 Characteristic factor table 类型 特征 光谱特征(SF) VNIR4-VNIR127、 VNIR133-VNIR150、 SWIR4-SWIR14、 SWIR20-SWIR38、 SWIR61-SWIR94、 SWIR116-SWIR119、 SWIR122-SWIR162 植被指数特征(VIF) TPI1、 TPI2、 TPI3、 TPI4、 TPI5、 TPI6、 TPI7、 TPI8、 TPI9、 NDVIg、 DVI650、 NDVI705、 mNDVI705、 reNDVI、 GNDVIhyper、 GNDVIhyper2、 ND800、 mSR、 mSR670、 mSR705、 mSR445、 mSR2、 SR1、 SR2、 SR3、 SR4、 SR5、 REP、 VOG1、 VOG2、 VOG3、 PRI2、 SIPI1、 CAI2、 PSRI、 CRI1、 CRI2、 ARI1、 ARI2、 WBI、 NDWI、 CLG、 NDII2、 TVI1、 TVI2 地形特征(TRF) Elevation、 Slope、 Aspect

表4 特征因子表 Table 4 Characteristic factor table

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< N)个特征, 使用信息增益、 信息熵或基尼指数来选取特征进行节点分裂; ③重复①, ②, 通过k次样本抽取和训练得到k个决策树的模型; ④最终采取集成学习理论将k个决策树进行线性组合, 每个决策树所占权重相等, 当一个待分类样本被输入时, 分类结果则由每颗决策树通过投票决定。

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

Table 5

表5(Table 5)

表5 分类方案 Table 5 Classification scheme 分类方案 属性类型 方案1 SF 方案2 VIF 方案3 TRF+VIF 方案4 SF+VIF 方案5 SF+TRF 方案6 SF+VIF+TRF 方案7 SF+FS 方案8 VIF+FS 方案9 TRF+VIF+FS 方案10 SF+VIF+FS 方案11 SF+TRF+FS 方案12 SF+VIF+TRF+FS

表5 分类方案 Table 5 Classification scheme

针对方案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

Fig.3

	Figure Option ViewDownloadNew Window
	图3 分类特征优选及特征重要性排序Fig.3 Feature selection and importance ranking of classification

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

Table 6

表6(Table 6)

表6 基于随机森林不同方案茶园提取精度表 Table 6 Tea plantation extraction accuracy table of different schemes based on RF 方案 OA/% Kappa PA/% UA/% F1-Score 方案1 88.49 0.79 87.61 79.14 0.83 方案2 88.96 0.80 87.42 80.82 0.84 方案3 90.01 0.82 89.72 81.97 0.86 方案4 88.85 0.80 87.61 80.37 0.84 方案5 89.07 0.80 88.93 79.71 0.84 方案6 89.01 0.80 88.90 80.20 0.84 方案7 88.92 0.80 88.36 80.18 0.84 方案8 89.45 0.81 88.74 81.51 0.85 方案9 90.83 0.83 90.41 83.43 0.87 方案10 89.26 0.80 88.55 81.19 0.85 方案11 90.07 0.82 90.74 81.10 0.86 方案12 90.70 0.83 90.86 84.09 0.87

表6 基于随机森林不同方案茶园提取精度表 Table 6 Tea plantation extraction accuracy table of different schemes based on RF

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

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 随机森林茶园提取分类结果(光谱特征+植被指数特征+地形特征+特征优选)Fig.4 Extraction results of random forest tea plantations (SF+VIF+TRF+FS)