基于多时相Sentinel-2A 影像的狼毒分布识别

2024-03-28房家玮胡念钊王怀玉刘咏梅

草业科学 2024年2期

房家玮，胡念钊，王怀玉，刘咏梅, 2

（1.西北大学城市与环境学院, 陕西西安 710127；2.陕西省地表系统与环境承载力重点实验室, 陕西西安 710127）

高寒草地具有气候调节、水源涵养、水土保持、生物多样性保护等多种生态功能，是青藏高原高寒生态系统价值的重要维系，也是当地畜牧业可持续发展的重要保障[1]。近年来，高寒草地呈现不同程度的退化趋势，狼毒(Stellera chamaejasme)等毒害草在退化草地的扩散日益加剧，对高寒草地生态系统的影响日益显现[2-3]。狼毒是瑞香科狼毒属多年生草本植物，其为顶生头状花序，花白色、淡红或黄色，花期为6 月下旬至7 月下旬，广泛分布于海拔2 300～4 200 m 干燥而向阳的高山山坡或河滩台地，是青藏高原危害最严重的有毒入侵植物之一[4-6]。狼毒根系发达，对土壤水分的利用效率高且对土壤养分转化速度快，在群落物种竞争中具有明显优势，已成为中、重度退化草地的优势种或建群种[7-8]。开展青藏高原高寒草地狼毒入侵调查，监测其发生发展区域，对草地资源管理和退化草甸修复至关重要。

近年来，随着影像分辨率的提高以及分类技术的不断完善，遥感技术在草地入侵植物分布识别中的应用日渐深入。依据入侵物种的物候特征并结合关键生长期的遥感影像，Li 等[9]利用IKONOS 4 m 多光谱影像对狼毒分布区域进行识别；Jana等[10]采用Rapid Eye 5 m 多光谱影像对入侵植物巨型猪草(Heracleum mantegazzianum)的生长进行监测；杨俊芳等[11]基于GF-1 和GF-2 全色及多光谱波段对黄河三角洲互花米草(Spartina alterniflora)的分布范围进行提取。上述研究使用单时相高分辨率影像取得了较理想的识别效果，但高分辨率影像在大区域尺度的广泛应用受到一定程度的限制[12]。同时，由于植物群落构成复杂多样，单时相中低分辨率遥感影像对物种的识别精度相对较低，多时相遥感影像能够反映植物的物候差异，为其精准识别提供有效数据源。卢元兵等[13]基于Landsat8 OLI的8 期影像构建混合三维和二维卷积的神经网络识别模型，提出了番茄(Solanum lycopersicum)等农作物的有效分类方案；史飞飞等[14]使用Landsat 8 OLI和GF-1WFV 4 月－11 月的作物生长期时序数据对枸杞(Lycium barbarum)种植区进行分类，取得了较好的识别精度。现有研究侧重于选取关键生长期的中高分辨率遥感影像，以大面积规律分布的农作物作为识别目标，而将多时相遥感影像与入侵植物独特的物候信息相结合，应用于毒杂草入侵调查监测，这方面的研究有待于进一步开展[15]。

本研究在青海省海北州祁连县选取狼毒广泛分布的退化草甸作为研究区，采用狼毒花期前和盛花期的Sentinel-2 多光谱影像，基于Google Earth Engine平台开展Sentinel-2 影像时间序列去云处理，采用多时相信息筛选狼毒识别敏感指数并构建最优特征组合方案，结合随机森林算法提取狼毒分布区域，评估多时相Sentinel-2 影像在区域尺度狼毒遥感识别中的适用性，为退化高寒草地毒杂草的综合防治提供技术支持。

1 研究区及数据源

1.1 研究区概况

研究区位于狼毒入侵广泛、草地退化严重的青海省海北藏族自治州祁连县(100°11′15.69″～101°10′34.21″ E，37°45′41.84″～38°16′30.49″ N)，面积约4 000 km2，海拔2 596～4 975 m。该区域为典型的高原大陆性气候，光照充足、气温日较差大、降水量较少、雨热同期。研究区草地类型主要为高寒草甸，群落优势种有狼毒、高山嵩草(Kobresia pygmaea)、披针叶黄华(Thermopsis lanceolata)、早熟禾(Poa annua)、异叶米口袋(Gueldenstaedtia diversifolia)、委陵菜(Potentilla chinensis)等。群落平均盖度约为48.9%，狼毒最大盖度达到57.7%，空间上呈现斑块状聚集分布的特征。研究区内与狼毒同期的植物群落主要呈现绿色，在盛花期大面积密集分布的狼毒花为粉白色，其他植物群落间或出现零散的黄色或紫色花，对狼毒群落光谱识别影响微弱，狼毒盛花期独特的物候特征为其遥感识别提供了基础。

1.2 遥感影像及预处理

本研究主要采用Sentinel-2A 卫星多光谱影像，下载自欧洲航天局(European Space Agency，ESA)数据中心(https://scihub.copernicus.eu/dhus)，产品级别为L2A 级，选取空间分辨率为10 m 的蓝波段B2(458～523 nm)、绿波段B3 (453～578 nm)、红波段B4 (650～680 nm)和近红外波段B8 (785～900 nm)。根据该区域狼毒生长的具体物候期[8]，选用2019年、2020 年和2021 年狼毒花期前(6 月20 日至30日)与盛花期(7 月5 日至20 日)内云量少于20%的影像作为狼毒识别的基础数据。

研究区在7 月－8 月多云雨天气存在大量的云雾遮挡区域，导致影像质量下降。本研究利用Google Earth Engine 平台对所用影像进行去云处理，获取每景影像的可用像元。首先利用Sentinel-2A 质量评估(quality assessment，QA) 波段中的云掩膜波段QA60 产品对影像进行云检测并掩膜去云；其次对去云后产生的空值区采用无云影像取中值进行补充，最终得到花期前和盛花期两个时相的影像。

1.3 野外采样数据

本研究于2020 年7 月上中旬－2021 年7 月上中旬在研究区开展了植物群落调查。根据区域内植物群落的典型性，沿青羊沟和八宝河谷地布设两条样带，共布设37 个样地。为了和Sentinel-2A 多光谱影像的空间分辨率相匹配，每个样地选在群落分布均匀且面积大于10 m × 10 m 的区域，并于区域中心布设1 m × 1 m 样方 2～3 个，共布设93 个样方(图1)。其中，草地样方(无狼毒生长样方) 39 个，狼毒生长样方54 个。样方调查内容主要包括植物种类、每个物种的数量和高度，利用数码相机获取样方盖度照片，并采用GNSS RTK 测定样方坐标、坡度和高程。对样方照片进行解译后得到样方群落盖度和物种分盖度，样方群落盖度为29.80%～70.00%。研究中依据狼毒群落盖度变化梯度设置狼毒样方，狼毒群落盖度为1.88%～57.68%，不同覆盖度狼毒群落样方如图2 所示。因所选Sentinel-2A 影像波段的分辨率为10 m，当样方的狼毒盖度小于10%时识别效果不佳，故进行剔除，最终研究中采用46 个狼毒样方和39 个草地样方作为分类样本数据，同期使用GNSS RTK 在研究区内采集139 个狼毒验证点，147个非狼毒验证点(包含其他草地群落、道路、河流、居民点及农田)作为分类结果验证数据。地面验证点主要记录经纬度坐标及类别信息。

图1 研究区植物样方分布Figure 1 Distribution of the sampling plots in the study area

1.4 其他数据

本研究同时使用了地形数据和植被数据。地形数据采用了美国国家航空航天局(National Aeronautics and Space Administration, NASA)和日本经济产业省(Ministry of Economy, Trade andIndustry, METI)联合发布的V3 版本30 m 分辨率ASTER GDEM 数据(https://earthdata.nasa.gov/)。研究区植被数据为《中华人民共和国植被图(1 ꞉ 1 000 000)》[16]，下载自中国科学院植物科学数据中心(https://www.plantplus.cn/cn)。将植被类型数据转为栅格数据，并与DEM 一起重采样至10 m，使之与遥感影像数据分辨率相统一。

2 研究方法

2.1 影像逐层掩膜

研究区内狼毒广泛生长于海拔2 600～4 200 m、年均温约0 ℃的高山及亚高山草地，分布地区植被类型主要有温带丛生禾草草原、温带丛生矮禾木、矮半灌木荒漠草原、禾草、薹草高寒草甸、蒿草和杂类草高寒草甸等。根据狼毒分布的高程范围和植被类型对研究区Sentinel-2A 影像进行地形及植被掩膜。

结合谷歌影像特征，对掩膜后研究区影像的耕地、道路、水体及居民地进行采样和监督分类，利用分类结果进一步剔除影像中耕地、道路、水体、居民地等非植被区域。通过逐层掩膜形成狼毒潜在发生区影像，主要包含狼毒群落与其他非狼毒草地群落，以减少其他地类对狼毒识别的干扰，提高分类精度。

2.2 分类特征优选

根据狼毒与其他草地植被在生长期明显的光谱差异，参考前人研究基于Sentinel-2A 影像计算相应的狼毒敏感指数，分析花期前和盛花期狼毒敏感指数的变化规律选取最优敏感指数[9]。初步选取的狼毒敏感指数有归一化指数NDVIb、NDVIb/r、NDVIr/b和乘积指数MIbg、MIbr，其计算公式[9]：

式中：ρred、ρgreen、ρblue、ρnir分别为红、绿、蓝、近红外波段的地表反射率。

依据公式(1)～(5)计算狼毒敏感指数，结合红、绿、蓝、近红外波段作为狼毒分类的初选特征，并采用相关性分析和随机森林重要性排序组合的二次降维方法对分类特征进行优选。首先，对所有特征两两之间进行Spearman 秩相关分析，根据相关系数大小筛选出信息重复率低的特征及特征组。其次，计算每个特征在随机森林模型中的贡献值进行随机森林重要性排序，依据重要性大小对特征及特征组做进一步筛选，获取有效的分类特征。最终，通过特征逐步添加的方式构建分类特征组合方案作为狼毒分类的信息源[17]。

2.3 遥感分类方法

本研究使用Python 语言基于PyCharm 平台建立随机森林分类模型，将分类特征组合方案作为信息源，46 个狼毒样方和39 个草地样方作为分类样本对研究区狼毒分布进行识别。其中80% 的样本作为训练集，20%的样本作为测试集。优化后的随机森林参数设置为，分类器个数为20 个，最小叶子数为4，最大深度为5。

2.4 分类精度评价

利用地面验证点数据，在ENVI 中计算混淆矩阵对分类结果进行精度评价，采用的参数有分类总精度(overall accuracy，OA)、Kappa 系数、制图精度(producer accuracy，PA) 和用户精度(user accuracy，UA)。分类总精度指正确分类样本数与总样本数之比，Kappa 系数是用于一致性检验的指标；制图精度指某类正确分类样本数与该类总样本数之比；用户精度指某类正确分类样本数与分为该类的样本数之比。采用以上参数对各特征组合方案的分类结果进行精度评价[18-19]。

3 结果与分析

3.1 狼毒敏感指数适用性分析

基于花期前及盛花期的Sentinel-2A 影像计算5 种狼毒敏感指数，并对MIbg、MIbr的指数值进行归一化处理，在此基础上计算狼毒样方及非狼毒样方的指数均值(表1)。花期前与盛花期狼毒与非狼毒样方的5 种敏感指数差值绝对值范围分别为0.020～0.060、0.062～0.299。其中MIbr、MIbg指数在花期前及盛花期狼毒和非狼毒之间的差值绝对值分别为0.060、0.299 和0.039、0.248，NDVIb、NDVIb/r、NDVIr/b指数在花期前与盛花期狼毒和非狼毒之间的差值绝对值分别为0.020～0.031 及0.062～0.070。在盛花期5 种植被指数更明显地体系了狼毒与非狼毒的差异，特别是MIbg、MIbr体现了两者之间的最大差异。总体来看，5 种指数均不同程度体现了花期前与盛花期狼毒与非狼毒样方的影像光谱差异，因而初步选择这5 种植被指数作为狼毒识别的敏感指数。

表1 花期前与盛花期敏感指数均值对比Table 1 Comparison of the sensitivity index mean in the periods before Stellera chamaejasme flowering and at full flowering

3.2 分类特征组合方案构建

分别基于花期前及盛花期的Sentinel-2A 影像计算5 个狼毒敏感指数，结合4 个多光谱波段，共提取18 项分类特征(表2)。KS (Kolmogorow-Smironov)检验显示A1、A5、A6 等12 个波段不符合正态分布(P< 0.05)，因此选用Spearman 秩相关分析检验各分类特征之间的相关性(图3)。A4、A8 与其他各特征之间的相关性均较弱；而部分波段或植被指数间的相关性较高，如A1 与A2、A3、A9、A12、A13，A5与A6、A7，A10 与A11，A14 与A15、A16，A17 与A18，则进一步依据随机森林重要性排序结果从这几组波段中各筛选出一个代表性波段。

表2 分类特征信息Table 2 Classification feature information

图3 分类特征相关性热力图Figure 3 Thermodynamic chart of correlation coefficients between

随机森林算法中决策树每个节点的变量选择具有随机性，在特征变量重要性统计中会出现一定程度的不确定性[17]。为降低误差将随机森林算法运行10 次并计算其均值，得到各特征重要性的最终排序结果(图4)。将各特征之间的相关性分析与其随机森林重要性排序结果相结合，最终选择A6、A14、A18、A8、A1、A11、A4 作为狼毒识别的最优特征。其中，A6、A14、A18、A8 为盛花期特征，A1、A11、A4 为花期前特征。

图4 分类特征重要性排序Figure 4 Feature importance ordering

按随机森林重要性排序依次递增添加特征来构建特征组合方案，分别采用4、5、6、7 个特征构成4 种组合方案作为分类信息源，4 个特征组合方案仅包含盛花期提取特征，5、6、7 个特征组合方案同时包含花期前与盛花期提取的特征。各方案包含的特征波段：A6、A14、A18、A8 (4 个特征)，A6、A14、A18、A8、A1 (5 个特征)，A6、A14、A18、A8、A1、A11 (6个特征)，A6、A14、A18、A8、A1、A11、A4 (7 个特征)。

3.3 狼毒分类与精度评价

分别基于4 种特征组合方案采用随机森林算法对掩膜的影像进行分类，生成研究区狼毒分布图(图5)。4 种方案提取的狼毒群落分布趋势基本一致，在地势平坦的谷地附近及河流两侧狼毒群落分布较密集，总体呈现大面积、斑块状、随机分散的分布特征。图5 细节图进一步显示，4 个特征组合方案提取的狼毒群落在整个区域的密集度及广泛度相对较高，5、6、7 个特征组合方案提取的狼毒空间分布模式相近。

图5 基于4 种特征组合方案的狼毒分类结果(a－d)及其细节图(e－h)Figure 5 Classification results of four feature schemes (a－d) and detailed maps (e－h)

研究区总面积约4 000 km2，经逐层掩膜形成的狼毒潜在发生区面积为2 761.51 km2。4 种组合方案分类结果的面积统计显示，狼毒占适生区草地总面积的比例介于17.60%～21.86% (表3)，其中，4 个特征组合提取的狼毒面积最大，占比为21.86%，5、6、7 个特征组合方案提取的狼毒面积相近，反映出4 种方案提取的狼毒面积相对较为稳定。

表3 4 种特征组合方案提取的狼毒面积Table 3 Area information of the four feature combination schemes

4 种组合方案的分类总精度介于78.67%～84.62%(表4)，Kappa 系数介于0.57～0.69，狼毒识别的生产者精度范围为81.29%～87.77%，用户精度为74.69%～81.88%。其中，6 个特征组合和7 个特征组合的分类总精度均大于80%，Kappa 均大于0.60，取得了较好的分类效果。6 个特征组合方案的狼毒分类生产者精度和用户精度分别为87.77%和81.88%，取得了最优的狼毒识别精度。通过分类特征的逐步增加，狼毒识别精度有相当程度的提升。多时相的6 个特征组合比单时相的4 个特征组合总体分类精度提高5.25百分点，Kappa 系数增加了0.10，狼毒生产者精度与用户精度分别提高了0.72 百分点和7.19 百分点，反映出多时相特征组合在狼毒群落识别中的优势。

表4 4 种组合方案狼毒分类精度评价Table 4 Accuracy evaluation of four feature schemes of Stellera chamaejasme classification

4 讨论

前人研究显示，去云及掩膜等预处理能够有效提高影像质量，减少分类干扰信息[20-21]。本研究中采用时间序列去云方法去除Sentinel-2A 影像的云污染，并根据狼毒生长的地形、植被等特征逐层掩膜狼毒适生区影像，这些处理能够有效地减少狼毒分类的干扰因素。在此基础上，根据狼毒的物候特征构建并筛选狼毒敏感指数，增强狼毒与其他植物群落的影像光谱差异，进一步为狼毒识别提供了有效的信息源[22]。Rapinel 等[23]及Noujdina 和Ustin[24]的研究显示，与单一时相影像相比，多时相影像对草地植物群落及入侵植物旱雀麦(Bromus tectorum)的识别精度更高，体现出多时相遥感在物种识别研究中的优势。高寒草地狼毒群落的空间分布分散、覆盖度差异大，单一时相的4 个特征组合分类总精度及狼毒用户精度均低于80%，而包含花期前与盛花期特征的6 个特征组合各项精度参数均大于80%，分类总精度达到84.62%，表明多时相Sentinel-2A 影像在大尺度狼毒群落遥感识别中具有较优的应用潜力。当狼毒呈稀疏分布(盖度小于10%)时与其他群落光谱差异较小，利用Sentinel-2A 的10 m多光谱影像对于低覆盖度狼毒群落的分布提取具有一定的局限性。

随机森林重要性排序方法为芦苇（Phragmites australis）、大蒜（Alium sativum）等植物的识别提取了有效分类特征[25]，不同特征组合的构建对比则能够进一步筛选出最优的分类特征方案[26]。本研究显示，将相关性分析与重要性排序二次降维方法相结合，可增强、提取和优化狼毒识别特征。随着特征数目的增加，特征组合方案的狼毒识别精度有不同程度的提高，多时相6、7 个特征组合的分类效果高于4 个特征组合。随机森林算法在高寒草地狼毒遥感识别中具有较好的应用精度，6 个特征组合结合随机森林取得了最好的狼毒识别结果，本研究与前人对于入侵植物毛蕊花(Verbascum thapsus)的识别研究[27]具有一致性，进一步验证了该方法在生态学分类中的准确性。