基于GF-2 PMS影像和随机森林的甘肃临夏花椒树种植监测

柳明星， 刘建红， 马敏飞，蒋娅， 曾靖超

(1.西北大学城市与环境学院，西安 710127; 2.西北大学陕西地表系统与环境承载力重点实验室，西安 710127)

0 引言

甘肃省临夏回族自治州(以下简称临夏州)位于甘肃省西南部，不仅是我国农业农村部重点扶持的“三州三区”深度贫困区之一，也是我国主要的花椒生产地。临夏州自然条件非常适合花椒种植，距今已有1 300多年的栽培历史。自1979年国家启动天然林保护计划后，花椒种植面积迅速扩大，并在1999年实施退耕还林工程后得到进一步扩大[1-2]。州内大部分乡镇80%以上的耕地和适宜的荒山、荒坡均栽植了花椒树，如莲花镇、南塬乡、银川乡、安集乡、河滩乡等[2-3]。花椒树核心种植带分布在黄河沿线以及环刘家峡库区周围等光热条件较好的河谷、川台塬及浅山坡。据统计，截至2020年底，全州共有57个乡镇267个行政村种植花椒树，全州花椒树种植总面积达到1万km2。

花椒树不仅是当地最重要的生态树种，也为当地农户带来了持续的经济收益。临夏州依托各项生态工程，借助国家精准扶贫和农业供给侧改革等政策机遇，因地制宜大力种植花椒树，积极发展地方特色经济，形成了以花椒为主导的生态扶贫产业[4-5]。2020年，临夏州实现花椒年产量2 100多万kg，总产量1.66万t，总产值6.6亿元。然而，花椒树常因极端环境变化发生病变、死亡，如自然灾害(干旱、泥石流等)、极端天气(持续高温、寒潮等)、花椒病虫害(黑胫病、枝枯病、蚜虫、铜色花椒跳甲等)等，这严重制约了当地生态保护和社会经济的协调发展[3，6]。因此，及时、准确地获取花椒树的种植面积和空间布局信息，对于及时、合理地调整生态工程实施方案，保障生态扶贫成效，促进经济、生态和社会可持续发展具有重要的意义。

在近几十年间，随着遥感技术和图像分析技术不断发展，植被遥感监测受到人们的广泛关注[7]，如植被类型识别[8-10]、作物产量监测[11-15]、植被动态变化分析[16-18]、植被生长及病虫害监测[19-21]、碳氮和生物物理参数机制研究等[22-26]，涉及到人工神经网络[27-30]、支持向量机[31-33]、随机森林[34-35]、深度学习等[36-37]多种算法。然而，在目前植被遥感监测中，大类作物往往容易获得更多的关注，如榆树[38]、小麦[11，39]、甘蔗[40]、玉米和大豆[41]、水稻[11，42]等，而对于具有特殊生态功能和经济效益的小类树种却明显关注不足，特别是以花椒树为代表的特殊或不常见树种[43]。虽然此类少数作物种植面积相对较少，但也是生态系统中不可或缺的组成部分。其不仅能够有效增加生物多样性，进而提高生态复原力和养分利用率，而且对于全球粮食和营养安全也有着重要的作用[44-46]。对特殊或不常见树种进行及时准确的监测是保障生态保护、环境治理及规划决策分析的重要内容。

与其他生态工程常见树种(如松树、杨树和桦树等)相比，种植花椒树每年都可获得较高的经济效益和持续的生态效益。综合全面地调查临夏州花椒树的种植状况，有助于提高农户经济收入，促进当地生态、经济、社会可持续发展。本文基于GF-2 PMS影像和随机森林算法，综合使用光谱波段、归一化植被指数(normalized difference vegetation index,NDVI)、纹理以及数字高程模型(digital elevation model,DEM)特征，设计了3种分类方案，深入探讨了花椒树遥感监测的可行性，以期为后续生态工程布局、其他生态树种或小类植被遥感监测提供指导意义。

1 研究区及其数据源

1.1 研究区概况

为了探究GF-2 PMS在花椒树遥感识别中的可行性，本文选择了临夏州花椒核心种植区进行分析。且受到地形、季风和水汽等原因，在花椒树生长季(3月25日—10月30日)中，研究区遥感影像含云量较高，仅获得了2景云量低于5%的GF-2 PMS遥感影像。研究区地理位置为N35.6°～35.8°，E103.1°～103.4°，总面积为1 026.59 km2(图1)。临夏州地处黄土高原与青藏高原、中原农业与西部牧区的过度地带，州内大部分地区属大陆性季风气候[47]。年平均气温为6.3℃，最高气温为32.5℃，最低气温为-27.8℃，年均降雨量为442 mm， 年平均蒸发量在1 198～1 745 mm之间，年均日照时数为2 572.3 h[48]。研究区海拔位于1 683～2 584 m之间，相对高差为901 m。

图1 研究区地理位置与范围

临夏花椒栽植区土层深厚、气候温和、光照充足，昼夜温差大，其自然条件非常适宜花椒树种植。临夏州种植的花椒品种以刺椒(伏椒)、棉椒为主，树龄普遍在20～30 a之间，郁密度普遍小于0.2。刺椒质量较优，但不耐冻，主要分布于河谷、库区周边光热条件较好、海拔于1 740 ～1 900 m的地区，占全州花椒总产量的30%左右。棉椒果实颜色艳红，产量高，较抗冻，但质量稍次，主要分布在南塬、安集等海拔在1 900 ～2 300 m之间的地区。临夏州的花椒树普遍呈顺地垄行状分布，具有明显的行距与株距(2 m×4 m)，每年的7—8月花椒成熟并进行采摘。

1.2 数据来源与预处理

1.2.1 GF-2 PMS影像

GF-2是我国首颗自主研发，且空间分辨率为亚米级的民用遥感卫星，分别搭载了2台含有全色(1 m)波段和多光谱(4 m)波段的相机(PMS)[49]。GF-2 PMS具有快速姿态机动、高定位精度、地表信息丰富等特点，在精细地物的识别和提取中发挥着重要的作用[49]。其基本卫星参数规格如表1所示。本文从中国资源卫星应用中心获取了2018年6月11日2景无云GF-2 PMS影像。使用ENVI 5.3软件对GF-2 PMS影像进行预处理。首先，对全色波段、多光谱波段分别做正射校正； 然后，对多光谱波段做辐射校正(包括辐射定标和大气校正)； 其次，对全色和多光谱波段做几何精纠正，再将全色、多光谱波段进行数据融合，得到空间分辨率为1 m的多光谱影像； 最后，将2景影像镶嵌在一起。

表1 GF-2 PMS卫星参数规格

1.2.2 地面参考数据

2019年7月开展了野外实地调查，并记录了396个样本的地理位置和作物种植信息。在调查过程中，发现不同的土地覆盖类型表现出明显的种植差异，如花椒树的种植密度(1 500～1 800棵/ hm2)较低,而玉米的种植密度(52 500～75 000株/hm2)相对较高。因此，根据实地调查所获取的信息，进一步在GF-2 PMS影像上随机选择了2 604个样本，共3 000个样本，包括训练样本和验证样本(图1)。样本数据集包含纯花椒(只种植花椒树的地块)、混合花椒(花椒树与其他作物混合种植的地块)、玉米、树林、茂密草地、稀疏草地、浑浊水体、清澈水体、人工地表、裸地共10种土地覆盖类型(表2)。

表2 本文采集的地面样本数据集

1.2.3 DEM数据

花椒树喜阳、不耐荫、耐寒，且根系耐水性较差，实际种植情况受地形条件影响较大[50]。根据当地花椒种植户反映，花椒树多分布于海拔1 740～2 300 m之间。因此，本文从地理空间数据云下载了空间分辨率为30 m的ASTGTM DEM数据作为花椒树识别的附加预测变量。为了与GF-2 PMS影像共同参与分类过程，基于ArcGIS 10.2软件计算了研究区的坡度和坡向变量，然后重采样至1 m分辨率。

2 研究方法

2.1 解译标志建立

建立正确的解译标志是遥感影像解译的基础，本文的主要目的是将纯花椒、混合花椒与其他土地覆盖类型进行有效区分。虽然GF-2 PMS影像的光谱信息有限(仅含4个波段)，不足以直接区分花椒树和其他植被，但是高分辨率提供了丰富的空间细节，可以很好地增强地物之间的可区分性。并且，花椒树种植较规律，会在遥感影像上表现出与其他土地覆盖类型不同的纹理信息。例如，在纯花椒树种植地块中，花椒树具有明显的行距和株距。在混合花椒地块中，由于混合种植了玉米、芝麻、蔬菜等其他作物，不如纯花椒地可以明显地区别于其他土地覆盖类型，但是依然可以看到花椒树呈点状分布的特征。2种花椒树种植类型的特征是其所特有的，与其他土地覆盖类型(玉米、树林、稀疏草地等)明显不同。结合现场调查和GF-2 PMS融合影像，确定了10种土地覆盖类型的解译标志(图2)。

(a) 纯花椒(b) 混合花椒(c) 玉米(d) 树林(e) 茂密草地

(f) 稀疏草地(g) 浑浊水体(h) 清澈水体(i) 人工地表(j) 裸地

2.2 NDVI计算

NDVI是评估绿色植被生长状态及植被空间分布密度的最佳指示因子[51]。由于其出色的植被表征性能，NDVI已经成为植被监测领域最重要的指标之一[52-54]。为了提高不同植被类型之间的可分性，本文计算了研究区的NDVI，其计算公式为：

(1)

式中Red和NIR分别代表红光、近红外波段亮度值，分别对应GF-2 PMS影像的B3和B4。

2.3 纹理特征提取

花椒树与其他土地覆盖类型具有明显的纹理差异，因此本文也将纹理特征作为附加变量参与分类过程。目前，最常用的纹理分析方法是灰度共生矩阵(gray level co-occurrence matrix，GLCM)[55-56]。因此，本文基于GLCM计算了GF-2 PMS影像的8种纹理特征： 均值(mean)、方差(variance)、同质性(homogeneity)、对比度(contrast)、异质性(dissimilarity)、熵(entropy)、角二阶矩(second moment)及相关性(correlation)。

由于GF-2融合影像数据量非常大，单景影像融合后数据量就达到了6 G，普通电脑端难以处理，为了确定纹理特征的最优窗口，本文以花椒种植最核心的区域(长为8.4 km，宽为5.2 km)为试验区，设置6种滑动窗口 (7×7,9×9,11×11,13×13,15×15,17×17)来计算GLCM纹理特征。

图3表示了8种纹理特征在不同滑动窗口下的反射率情况。从图3可得，试验区内各土地覆盖类型的均值、同质性、异质性、对比度随窗口的增大几乎不变或变化很小； 各土地覆盖类型熵的差异随着计算窗口的增大而增大，而不同土地覆盖之间相似度和角二阶矩的差异却随着窗口的增大而减小。从8种纹理特征区分不同土地覆盖类型的能力上看，均值、熵和相似度可以最大程度地反映不同地物之间的差异。因此，最终选择了在11×11窗口下计算的均值、熵和相似度，并参与分类过程。

(a) 均值 (b) 熵 (c) 相似度

(d) 角二阶矩 (e) 同质性 (f) 方差

(g) 异质性 (h) 对比度

2.4 随机森林分类

随机森林是通过设置多棵随机森林树对样本进行训练并预测结果的监督分类器，常用于处理分类、回归等问题[57]。在分类模式下，多棵树分别对样本进行投票并根据投票率决定最终分类结果。首先，从样本集中通过自举重采样的方法随机有放回地选出训练样本集，其中约2/3的样本用于模型训练，剩余1/3样本用于模型的内部验证(即“袋外误差”)。然后，对每一个训练样本分别设置对应的决策树模型，持续分裂，直到该节点的所有训练样本都属于同一类型。最后，将生成的多棵决策树组成随机森林，按照投票概率决定最优分类[57]。随机森林算法已经广泛应用于林业[35]、农业[58]、土地覆盖[59]等众多领域中。

基于随机森林算法，本文设计了3种分类方案来探究不同输入特征在花椒树种植识别中的能力： ①S1,仅使用光谱波段; ②S2,组合使用光谱波段、NDVI和DEM； ③S3，组合使用光谱波段、NDVI、纹理特征和DEM。

本文随机森林算法是基于RStudio (4.0.2 版本)中的“randomForest”程序包实现的。其中，通过反复实验，发现3种分类方案均在随机森林数ntree为500时表现最好，其中每个节点上使用的特征数位mtry为输入特征总数的平方根。

2.5 准确性评估

为了评估不同分类方案的分类效果，基于随机分层抽样方法，在分类结果中分别从每种土地覆盖类型中随机选了300个，共3 000个样本，采取目视解译的方法对每个样本对应的真实类型进行记录。然后，使用以下指标对S1，S2和S3分类方案进行评估： 总体精度(overall accuracy，OA)，Kappa系数(kappa coefficient)，制图精度(producer’s accuracy，PA)，用户精度(user’s accuracy，UA)以及F1得分[60]。其中，F1得分计算公式为：

(2)

3 结果与分析

3.1 花椒树与其他土地覆盖类型的光谱差异

基于GF-2 PMS融合影像和地面样本，首先分析了花椒树与其他土地覆盖类型之间的光谱差异(图4)。可以发现，10种土地覆盖类型在不同波段上的反射率比较接近，光谱曲线十分相似，从而造成了明显的光谱重叠，特别是同为植被的纯花椒、混合花椒、玉米、林地、茂密草地、稀疏草地。GF-2 PMS仅含4个光谱波段，光谱信息弱，仅使用GF-2 PMS的光谱信息难以获得较好的花椒树分类效果。

图4 基于GF-2 PMS的土地覆盖类型光谱曲线

为量化不同土地覆盖类型之间的可分离性，计算了Jefferies-Matusita (JM)距离。JM值通常位于0～2之间，值越高说明其可分离性越强[61]。本文将JM距离划分为4个等级，代表着不同等级的可分离性： 可分离性好[1.9，2.0)、可分离性较好[1.8，1.9)、可分离性较差[1.0，1.8)、可分离性差[0，1.0)(图5)。其中，C1，C2，C3，C4，C5，C6，C7，C8，C9，C10分别代表裸地、人工地表、茂密草地、树林、玉米、混合花椒、纯花椒、 稀疏草地、浑浊水体、清澈水体。红色表示可分离性满足分类要求，灰色表示可分离性不满足分类要求。 圆的大小表示可分离性的大小。结果表明，单独使用GF-2 PMS影像的光谱信息时(图5(a))， 除浑浊水体和清澈水体能够明显区分，其余土地覆盖类型之间的可分离性都比较低，JM值多低于1.8。而增加NDVI和DEM后，整体的可分离度有所提高(图5(b))，但还有部分地物类型之间的JM值低于1.0。进一步增加纹理特征后，少部分土地覆盖类型之间的JM值在[1.0，1.8)之间，但整体土地覆盖类型的可分离性明显提高(图5(c))。说明在GF-2 PMS光谱信息的基础上，增加NDVI和DEM，特别是纹理特征，可以极大地增加地物之间的可区分性，提高分类精度。

(a) 基于光谱波段的可分离性 (b) 基于光谱波段、NDVI和 DEM的可分离性 (c) 基于光谱波段、NDVI、纹理特征 和DEM的可分离性

图5 基于JM距离的各地物类型可分离性

3.2 不同方案的分类精度

为了深入探究光谱波段、NDVI、纹理特征以及DEM在花椒树分类识别中的影响程度，分别设计了3种不同的分类方案参与随机森林分类过程。基于混淆矩阵，结果显示，S1,S2和S3的OA分别为65.90%，67.67%和74.43%，Kappa系数分别为0.62，0.64和0.72。与S1相比，增加NDVI和DEM时，OA提高了1.77个百分点； 进一步增加纹理特征，OA提高了8.53个百分点。这表明增加植被指数和地形变量只能略微改善分类效果，而增加纹理特征却明显提高分类精度。3种分类方案中不同土地覆盖类型的PA,UA和F1分数结果如图6所示。

(a) PA结果 (b) UA结果 (c) F1得分结果

图6 S1，S2和S3分类方案的准确性评估结果

对于单一的土地覆盖类型来说，根据PA，UA和F1得分，分类效果为： S3>S2>S1。水体在3种分类方案中都表现出较高的准确性，而裸地、树林、纯花椒的识别效果则低于其他土地覆盖类型。对于花椒树而言，纯花椒在S3方案中，UA和F1得分最高，分别是53.92%和54.73%； 对于混合花椒而言，S3的PA，UA和F1得分均高于60.18%。对于花椒树而言，除PA外，混合花椒的UA和F1得分都高于纯花椒，说明混合花椒的分类效果优于纯花椒。综上所述，S3方案的分类效果是最好的，这也进一步说明了纹理特征在花椒树分类识别中的重要性。

3.3 花椒树种植的规模和空间格局

基于最优分类方案(S3)，获得了研究区2018年花椒树的空间分布情况(图7)。研究结果显示，研究区内的花椒树广泛种植于刘家峡库区岸边及东南部，总种植面积为231.59 km2，占研究区总面积的22.56%。其中，纯花椒种植区主要分布在研究区东南部，种植面积为189.06 km2。而混合花椒种植区多集中在靠近农村居民点刘家峡水库东南侧，有利于农户在管理花椒的同时种植其他作物，其种植面积为42.53 km2。

图7 研究区花椒树空间分布

为了更好地可视化随机森林在花椒树种植识别中的效果，对纯花椒和混合花椒的典型种植区域进行放大分析(图8)。其中，区域A为纯花椒种植地块，区域B为混合花椒种植地块。研究结果显示，在区域A中，存在小部分区域被错分为混合花椒，但是分类结果还是以纯花椒为主； 在区域B中，分类结果基本上为混合花椒。整体来说，虽然存在一定的错分和误分现象，但随机森林分类结果是比较准确的。

(a) 区域A纯花椒GF-2影像(b) 区域A纯花椒分类结果

(c) 区域B混合花椒GF-2影像(d) 区域B混合花椒分类结果

图8 随机森林分类结果局部放大

3.4 花椒树分布的地理特征

进一步，基于海拔和坡度信息分析了花椒树的地理分布特征(图9)。从图9(a)中可知，研究区[2 100,2 200) m海拔范围内的花椒树分布最为密集，种植面积占比为20%。其次是[1 683,1 800) m和[2 200，2 300) m海拔区域内的花椒树，其面积占比为18%。总体上，随着海拔的增加呈现出“先减少-再增加-再减少”的趋势。这种趋势变化应该与花椒树种植的品种有关，不同品种的花椒树对地形的适应性不同。在实地调查中发现，研究区内种植的花椒品种主要为刺椒和绵椒。刺椒不耐冻，主要分布于河谷、库区周边光热条件好、海拔在[1 740，1 900) m的地区； 而绵椒较为抗冻，主要分布于海拔在[1 900,2 300) m之间的地区。总的来说，90%的花椒树集中分布在[1 683,2 300) m海拔范围内。

(a) 花椒海拔分布情况(b) 花椒坡度分布情况

此外，按照临夏州的退耕还林政策，本文将花椒树分布的坡度划分为6个等级： [0，5)°，[5，8)°，[8，15)°，[15，20)°，[20，25)°和≥25°。从图9(b)可知，[8，15)°坡度范围内种植的花椒树最多，分布频率达到了24%。而50%的花椒树都种植在坡度≥15°的地区，这与临夏州的生态恢复工程的实际情况相符合。

4 讨论

作为重要的生态树种，与其他常见生态树种(如松树、杨树和桦树等)相比，花椒树每年都能获得客观的经济收益。但是，在目前林业生态遥感监测中，多集中于基于NDVI的植被动态变化、植被覆盖度变化监测等的宏观监测[62-64]，少有针对各类生态树种的微观分类研究，特别是以花椒树为代表的特殊或不常见树种。

本文基于GF-2 PMS影像和随机森林算法对花椒树种植进行监测识别，尽管可使用的遥感影像数量有限，但研究结果表明了GF-2 PMS在花椒树遥感监测中具有较高潜力。花椒树种植识别的最高精度达到了74.43%。当仅利用GF-2 PMS的光谱信息时，花椒树与其他土地覆盖类型难以完全分开，但是通过增加NDVI、纹理特征和DEM则可以进一步提高可分离性，获得更好的分类结果，这也与前人的研究相一致[65-67]。

然而，本文的研究还存在一定的局限性。首先，使用基于像元的分类方法(随机森林)时，对于高分辨率的遥感影像来说，会存在一定程度的错分、漏分现象。因此，基于最优分类方案(S3)与随机森林，本文进一步分析各土地覆盖类型的分类结果一致性 (图10)。分类结果与地面真实类型一样，则定义为分类结果一致，若出现偏差，则定义分类结果不一致。从图6可知，浑浊水体和清澈水体的分类精度较高，因此仅对纯花椒、混合花椒、玉米、树林等其他8种土地覆盖类型进行分析。

图10 基于最优分类方案(S3)的土地覆盖类型结果一致性分析

基于分层抽样方法，从地面参考数据集中随机抽取了约1/3共754个样本进行分析。对于纯花椒而言，主要错分为稀疏草地、裸地及人工地表，其错分率分别为29%，11%，4%。对于混合花椒而言，主要错分为玉米、茂密草地及稀疏草地，其错分率分别为14%，4%，4%。在其他6种土地覆盖类型中，纯花椒在稀疏草地中漏分情况最严重，其漏分率为18%，其次是人工地表和裸地，其漏分率分别为10%和7%。而混合花椒在树林、茂密草地、玉米、稀疏草地中的漏分率分别为19%，13%，9%，8%。究其原因，在于花椒与其他土地覆盖类型相似的光谱和纹理特征，特别是纯花椒与混合花椒、纯花椒与裸地、纯花椒与稀疏草地、混合花椒与玉米、混合花椒与树林、混合花椒与茂密草地(图4)，这很大程度上导致了分类精度相对偏低。未来工作中应该考虑面向对象的分类方法，以期能够更好地提高分类精度。

此外，各种生态工程已经实施多年，如三北防护林计划(1978年)、天然林保护计划(1979年)、退耕还林工程(1999年)等，无论是从时间序列上进行长期动态监测，还是大面积规模监测都是很有必要的。GF-2 PMS卫星于2014年发射，在轨运行时间较短，可用遥感数据并不多，且难以获取2014年之前的遥感数据，暂时无法实现生态树种长期动态监测的目标。另一方面，由于GF-2 PMS卫星的亚米级空间分辨率，数据量巨大，对于普通电脑客户端，难以实现大区域地区监测。现有的遥感卫星中，如Landsat8 OLI或 Sentinel-2 MSI等，虽然空间分辨率不足以计算纹理信息，但是都具有丰富的光谱信息[68-69]，且Landsat 系列卫星可以实现长期工程效益监测。因此，进一步改善生态工程树种动态监测的可能方法是将高分辨率影像(如GF-2)与多光谱影像(如Landsat8 OLI)等组合探索。

5 结论

本文以GF-2 PMS融合影像为基础，通过光谱波段、NDVI、纹理特征(均值、熵和相似度)以及DEM信息共4种分类特征，基于随机森林分类算法分别设计了3种不同的分类方案，以此探索甘肃省临夏州花椒树种植遥感监测的可行性。主要结论如下：

1)4种分类特征在花椒种植识别中都表现出不同的作用效果，其中纹理特征最重要。仅使用光谱波段时，总体精度为65.90%； 增加NDVI和DEM信息，总体精度为67.67%； 进一步增加纹理特征，总体精度达到了74.43%，共提高了8.53个百分点，有效改善了分类效果。

2)2018年研究区内的花椒树主要种植在黄河沿岸和刘家峡库区周边，其总种植面积为231.59 km2，占研究区总面积的22.56%。其中，只种植花椒树的面积为189.06 km2，混合种植花椒树的面积为42.53 km2。

3)花椒树分布具有特定的地理分布特征。90%以上的花椒树种植在[1 683,2 300) m海拔范围内，且随海拔升高呈现出 “先减少-再增加-再减少”的分布特征； 58%的花椒树分布在[8，15)°坡度范围内，这与当地退耕还林政策相符合。

总的来说，GF-2 PMS影像在花椒树种植监测中具有较高的潜力。开发对花椒树的遥感识别方法，不仅有助于当地生态产业调控和后续生态工程布局，并且对其他地区的生态树种或小类植被物种开展相关遥感监测工作也具有较强的指导意义。