黑龙江省黑河市(图1), 包括北安市西北部、五大连池市中南部及克东县东北部, 地处47° 42′ ~ 51° 03′ N、124° 45′ ~129° 18′ E, 面积68 726 km² , 春季气温低、干旱多风, 夏季温热多雨, 秋季凉爽舒适, 冬季漫长、严寒而干燥, 冬长夏短、四季分明。全市平均年降雨量500~550 mm, 有效积温1950~2300℃, 日照时数2560~2700 h, 无霜期90~120 d, 年均气温-1.3~0.4℃, 最高气温38.2℃, 最低气温-40℃。耕地面积192万公顷, 是国家重要商品粮基地和绿色食品主产区。主要耕地土壤为黑土, 其次包括草甸土、暗棕壤及沼泽土。第一产业是黑河市的主要产业, 2015年, 农林牧渔总产值达264.5亿元, 主要农作物为玉米、水稻、小麦、薯类和大豆。自2015年以来农业部在该市开展粮豆轮作试点, 优化作物种植结构, 适当调减玉米种植面积, 增加大豆种植面积。2015年共落实轮作面积1800 hm²。

图1

Fig. 1

	Figure Option ViewDownloadNew Window
	图1 研究区地理位置Fig. 1 Location of study area

1.2.1 遥感数据处理 主要使用Landsat 8卫星数据, 该卫星由NASA于2013年2月11日发射, OLI是其携带的主要传感器, 包括7个波段, 分别是海岸/气溶胶(430~450 nm)、蓝(450~510 nm)、绿(530~590 nm)、红(640~670 nm)、近红(850~880 nm)、短波红外1 (1570~1650 nm)和短波红外2 (2110~ 2290 nm), 空间分辨率均为30 m。根据研究区主要农作物玉米和大豆的生育期特征, 本文选取了覆盖整个研究区的2014年8月7日Landsat 8 OLI卫星影像, 该景影像无云覆盖, 只在东部有少量薄雾。

使用ENVI 5.0软件对获取的遥感影像进行辐射定标、大气校正和几何精校正处理。其中辐射定标公式如下:

\[{{L}_{z}}({{\lambda }_{z}})=\text{Gain}\times \text{DN}+\text{Bias}\] (1)

式中, \[{{L}_{z}}({{\lambda }_{z}})\]为传感器入瞳处的光谱辐射亮度(W m^-2 sr^-1 μ m^-1), Gain为定标斜率, DN为影像灰度值, Bias为定标截距, Gain及Bias都由卫星数据供应方提供, 可从原始Landsat影像的元数据文件中读取。Landsat 8 OLI各波段的辐射定标系数如表1所示。

表1

Table 1

表1(Table 1)

表1 Landsat 8 OLI影像各波段辐射定标系数 Table 1 Radiometric calibration coefficient of Landsat 8 OLI image 波段 Band 定标斜率 Gain 定标截距 Bias 海岸蓝波段 Coastal aerosol 0.01298 -64.89967 蓝 Blue 0.013292 -66.45805 绿 Green 0.012248 -61.24053 红 Red 0.010328 -51.64146 近红外 Near infrared 0.0063204 -31.602 短波红外1 SWIR 1 0.0015718 -7.85913 短波红外2 SWIR 2 0.00052979 -2.64895

表1 Landsat 8 OLI影像各波段辐射定标系数 Table 1 Radiometric calibration coefficient of Landsat 8 OLI image

采用ENVI/FLAASH大气校正模块进行大气校正, ENVI/OLI校正模块进行几何校正。

1.2.2 地面样方调查 地面样方是进行监督分类的重要基础, 根据研究区作物分布情况和各类作物的光谱特征, 在整个研究区随机选取均匀分布的131个样方地块, 利用同一研究区更高分辨率的RapidEye影像进行目视解译分类, Landsat原始影像及样方的分布如图2所示。样方的总面积为5961.15 hm², 其中大豆面积为1351.98 hm², 占样方总面积的22.68%, 玉米面积为1404.09 hm², 占样方总面积的23.55%, 其他地类主要包括城镇、林地、水体、湿地、裸地、水稻等, 共3205.08 hm², 占样方总面积的53.77%。

图2

Fig. 2

	Figure Option ViewDownloadNew Window
	图2 研究区Landsat 8 OLI影像及样方分布 a: Landsat 8原始影像及样方分布; b: 样方原图; c: 样方分类结果。Fig. 2 Landsat 8 OLI image and distribution of ground sample in study area a: Landsat 8 image and distribution of ground samples; b: original image of sample; c: classification of sample.

1.2.3 作物发育时期 春玉米从播种开始, 依次经历出苗、三叶、七叶、拔节、抽雄、乳熟、成熟等发育时期, 研究区每年4月下旬开始播种, 8月上旬成熟, 9月下旬以后收获。播种— 出苗、出苗— 三叶、三叶— 七叶、七叶— 拔节、拔节— 抽雄、抽雄— 乳熟、乳熟— 成熟等7个生长阶段历时平均分别为9、9、12、23、15、29和19 d, 全生育期所需时间约为116 d。以北安地区4月25日播种计算, 7个关键生育期起始日期分别为5月4日、5月13日、5月25日、6月17日、6月22日、7月21日和8月9日。

大豆从播种开始, 一般经历种植期、发芽期、早期生长期、出枝期、开花期、结荚期、收割期等发育时期, 种植期一般位于6月份, 通常在种植1~2周内发芽, 在发芽后25 d内为早期生长期, 作物高度可以达到15~21 cm。一般在发芽40 d后大豆可生长出1~6枝。在发芽期(6月中旬至7月下旬)和开花期(7月中旬至7月下旬)之间的大豆生长情况将决定大豆开花的数量, 并直接影响大豆的产量; 开花期一般在种植后45~50 d左右, 持续约30 d; 结荚期一般在7月下旬和8月上旬, 而收割期则一般为9月或10月份。

综合考虑研究区内作物的发育时期, 选择位于8月7日的卫星影像进行单时相作物分类识别。该时期研究区主要地物类型(大豆、玉米、林地、城镇用地)的光谱曲线如图3所示, 可以看出, 植被生长旺盛, 大豆的近红外和2个短波红外波段的反射率都处于最高水平, 玉米和林地则较为接近, 但是在海岸、蓝、绿波段林地反射率普遍高于玉米。该时期正处于春玉米抽雄至乳熟期及大豆结荚期内, 作物生长旺盛, 作物光谱类内一致性强, 而类间差异显著, 避免作物生长早期由于播种早晚导致光谱差异以及作物生长晚期不同成熟、收割时间导致的光谱差异有利于基于遥感影像的农作物识别; 同时, 根据王利民等^[30]研究, OLI影像具有的短波红外波段对于玉米-大豆识别具有较为显著的作用, 因此, 本文选用8月7日单时相的OLI数据即可识别研究区玉米-大豆。

图3

Fig. 3

	Figure Option ViewDownloadNew Window
	图3 研究区主要地物类别光谱曲线Fig. 3 Spectral curves of main ground objects in study area

1.3.1 技术思路 整体研究思路如图4所示。针对研究区的作物分布情况, 选用2014年第218天的单景Landsat-8/OLI卫星影像, 选取合适的样本数据, 分别使用最大似然分类、支持向量机分类、随机森林分类3种分类方法对研究区的主要作物玉米-大豆进行分类识别, 利用研究区更高分辨率的RapidEye影像目视解译成果作为分类真值影像, 对各种方法的分类精度进行评价, 同时分析各类方法的分类时间, 以此评价各分类方法的适用性。另外在原始影像基础上, 依据玉米-大豆-其他三大类别具有较为明显的差异这一原则, 另外增加主成分变换的前4个波段(包含原始数据的绝大部分有效信息, 并去除噪音)、归一化植被指数NDVI (作物分类识别常用指数, 可以区分植被和非植被)、归一化水体指数NDWI (区分水库以及湿地、水稻等)等额外的信息数据, 共13个波段, 再次使用3种分类方法进行作物的分类, 以此评价增加信息对分类精度的影响, 分析利用OLI数据进行大豆-玉米分类是否有必要添加额外信息。

图4

Fig. 4

	Figure Option ViewDownloadNew Window
	图4 研究技术流程Fig. 4 Technical flow chart of the study

1.3.2 最大似然算法 最大似然法又称最大概似估计或贝叶斯(Bayes)分类方法, 是监督分类方法的一种。最大似然法以统计学原理, 利用最大似然比贝叶斯判决准则建立非线性的判别函数集, 并假定各类分布函数为正态分布, 通过训练样本, 计算待分类像元对各类别的归属概率, 哪一类归属概率高即将其划分为此类。该方法是当前农情遥感监测作物分类提取业务化流程中经常使用的监督分类方法, 具有分类精度较高、分类结果稳定可靠、分类速度快的优势^[31]。

1.3.3 支持向量机 支持向量机(SVM)是Cortes和Vapnik等于1995年首先提出的机器学习分类方法^[32]。它建立在统计学VC维(Vapnik-Chervonenkis Dimension)理论和结构风险最小原理基础上, 根据有限的样本信息在模型的复杂性和学习能力之间寻求最佳折衷, 以期获得最好的推广能力。对于影像而言, 将影像的多个波段灰度值视为一个向量, 将该向量映射到一个更高维的空间里, 并构造一个最大间隔的超平面, 即在将数据分开的超平面两边构建2个互相平行的超平面, 使两者间距离最大。平行超平面间的距离或差距越大, 分类器的总误差越小, 因此可以达到最优化分类的目的。该法能在较小样本情况下自动学习样本分类知识, 获得较高精度的分类结果, 在多个领域广泛应用。

1.3.4 随机森林分类 随机森林分类(RFC)是Breiman于2001年提出来的一种较新的多决策树分类方法^[33], 该方法通过在数据上及特征变量上的随机重采样, 构建多个CART类型决策树(不剪枝), 通过多决策树投票的方式确定数据的类别归属。随机森林方法对于遥感影像分类具有很好的抗噪声性能, 分类精度较高。该法利用样方数据自动构建分类决策树, 属于监督分类的一种。

(1) 随机森林算法从原始样本数据集中抽取N个训练样本集, 每个训练样本集都是从原始样本集中随机有放回地抽取大约2/3, 剩余的约1/3作为验证样本, 称为袋外数据(out-of-bag, OOB) 进行内部误差估计, 并利用OOB数据计算各特征变量的重要性。

(2) 随机森林在构建每颗树时并不选用全部特征, 而是通过有放回随机抽取的方式, 从原始的特征集中(假设共K个特征)抽取k个(k≤ K)特征作为决策树分类依据, 构建数据特征预测变量集。一般情况下, k值设置为K的平方根即可。

(3) 根据选取的训练样本、验证样本, 以及特征预测变量集, 参照CART决策树构建方法, 通过递归建立一个分类二叉树。假设样本有k个属性特征, 对于每个属性特征, 参照Gini指数选取一个最佳划分值x, Gini指数越小即认为划分后的类别中杂质含量越低, 分类纯度越高。假设一个样本共有m类, 则二叉树节点A的Gini指数计算方式如下。

\[\text{Gini}(\text{A})\text{=}1-\sum\nolimits_{i=1}^{m}{p_{i}^{2}}\] (2)

式中, p_i代表属于i类的概率, 当Gini(A)=0时所有样本属于一类。递归的过程则是针对当前节点, 尝试样本每一个属性特征, 计算各属性变量中Gini指数最小的值作为该节点的最佳属性划分值, 构建一个最优分支子树。根据以上分裂规则, 对样本进行充分的二叉树生长, 构建一个完整的CART树, 一般情况下不对该树进行剪枝操作。

(4) 重复步骤3, 直到构建完成N颗分类树, 进而形成一个随机分类树的森林, 将影像的每一个像元使用所有的分类树进行分类, 采用多数投票方式综合分类结果, 确定该像元的最终从属类别。

对于随机森林采用样本和特征的双重随机抽样构建决策树, 因此即使不对分类树进行剪枝操作也不会出现传统CART决策树过拟合的现象。

1.3.5 精度验证方法 主要基于5 m空间分辨率Rapideye影像, 针对整个研究区的玉米、大豆及其他地物类型进行监督分类, 结合地面调查, 采用目视方法修正, 目视解译结果作为研究结果精度验证的数据, 图5给出基于Rapideye影像目视解译结果。采用混淆矩阵、Kappa系数、总体分类精度、制图精度和用户精度5种方式描述和比较分类精度, 相关定义及详细表述可参照文献[34, 35, 36]。

图5

Fig. 5

	Figure Option ViewDownloadNew Window
	图5 基于RapidEye影像的目视解译结果 a: Rapideye影像(5/4/3波段); b: Rapideye影像的目视解译结果。Fig. 5 Visual interpreting result based on the RapidEye image a: Rapideye image (5/4/3 band); b: Result of manual visual interpretation based on RapidEye image.

利用3种分类方法和相同的地面样方进行地物分类, 结果如图6所示。随机森林方法的决策树数量设定为100, 节点分裂时输入的特征变量数设置为所有特征数量的平方根; 支持向量机方法的核函数选择为径向基函数(radial basis function, RBF), 核函数的Gamma值设置为0.071; 最大似然分类方法设置分割概率阈值为单一阈值。表2表明, 最大似然分类方法、支持向量机方法、随机森林方法的总体分类精度分别为91.68%、91.49%和94.32%, Kappa系数分别为0.87、0.87和0.91, 可以看出, 随机森林方法比另两种方法有更高的分类精度。从分项看, 最大似然分类方法玉米的产品精度较高, 但是用户精度较低, 说明存在较多的错分情况(其他错分为玉米), 而其他的用户精度较高, 但是产品精度较低, 说明存在较多的漏分情况(其他错分为大豆和玉米); 而支持向量机方法情况与最大似然分类方法类似; 随机森林相比其他2类方法, 各类作物无论是制图精度还是用户精度, 都有明显的提升, 说明随机森林方法比传统监督分类方法具有更高的作物分类识别能力。

图6

Fig. 6

Figure Option ViewDownloadNew Window

图6 基于原始影像3种方法分类结果 a: 最大似然分类结果; b: 支持向量机分类结果; c: 随机森林分类结果; d: 最大似然分类结果局部; e: 支持向量机分类结果局部; f: 随机森林分类结果局部。Fig. 6 Classification results by three methods based on original image a: maximum likelihood classification result; b: support vector machine classification result; c: random forest classification result; d: part of maximum likelihood classification result; e: part of support vector machine classification result; f: part of random forest classification result.

表2

Table 2

表2(Table 2)

表2 基于原始影像3种分类方法的混淆矩阵 Table 2 Confusion matrix of three classification methods based on original image 作物 Crop 分类方法 Method 大豆 Soybean (pixel) 玉米 Maize (pixel) 其他 Other (pixel) 总计 Total (pixel) 制图精度 Mapping accuracy (%) 大豆 Soybean MLC 864849 15393 56545 936787 91.03 SVM 878082 7498 66202 951782 92.42 RFC 911841 9285 65006 986132 95.98 玉米 Maize MLC 21443 1134475 76472 1232390 92.38 SVM 35185 1140200 93000 1268385 92.85 RFC 4158 1186858 68423 1259439 96.65 其他 Other MLC 63758 78172 1436164 1578094 91.52 SVM 36783 80342 1409979 1527104 89.85 RFC 34051 31897 1435752 1501700 91.50 总计 Total (pixel) 950050 1228040 1569181 3747271 用户精度 User accuracy (%) MLC 92.32 92.05 91.01 SVM 92.26 89.89 92.33 RFC 92.47 94.24 95.61 总体精度 Overall accuracy (%) MLC 91.68 Kappa系数 Kappa coefficient MLC 0.87 SVM 91.49 SVM 0.87 RFC 94.32 RFC 0.91 MLC: 最大似然分类; SVM: 支撑向量机分类; RFC: 随机森林分类。 MLC: maximum likelihood classification; SVM: support vector machine classification; RFC: random forest classification.

表2 基于原始影像3种分类方法的混淆矩阵 Table 2 Confusion matrix of three classification methods based on original image

分别计算原始影像的归一化植被指数NDVI、归一化水体指数NDWI, 并进行主成分变换, 提取主成分波段的前4个波段, 共获得6个辅助分类波段, 与原始OLI影像的7个波段叠加组合, 形成13个波段的待分类影像。与原始影像分类相同, 分别采用3种分类方法, 使用相同的地面样方数据进行作物分类, 并与真值影像对比, 结果如图7所示。分类精度如表3所示, 增加辅助信息后, 最大似然分类和支持向量机的分类精度基本没有变化, Kappa系数也未提高; 而随机森林方法的总体分类精度则由原来的94.32%提高到95.81%, 提高了1.49个百分点, Kappa系数则由0.91提高到0.94, 表明辅助信息的加入可以在一定程度上提高作物的分类识别能力和精度, 但是整体提升的效果有限。在添加了辅助信息后, 相比最大似然分类方法, 总体识别精度从90.22%提高到95.81%, 提高了5.59个百分点, 大豆的制图精度从90.22%提高到98.32%, 提高了8.1个百分点, 玉米的用户精度从84.27%提高到94.88%, 提高了10.61个百分点。

图7

Fig. 7

Figure Option ViewDownloadNew Window

图7 增加辅助特征后的3种方法分类结果 a: 最大似然分类结果; b: 支持向量机分类结果; c: 随机森林分类结果; d: 最大似然分类结果局部; e: 支持向量机分类结果局部; f: 随机森林分类结果局部。Fig. 7 Classification results of three methods based on the stacked image a: maximum likelihood classification result; b: support vector machine classification result; c: random forest classification result; d: part of maximum likelihood classification result; e: part of support vector machine classification result; f: part of random forest classification result.

表3

Table 3

表3(Table 3)

表3 增加辅助特征波段后的3种分类方法分类结果混淆矩阵 Table 3 Confusion matrix of three methods based on the stacked image 作物 Crop 分类方法 Method 大豆 Soybean (pixel) 玉米 Maize (pixel) 其他 Other (pixel) 总计 Total (pixel) 制图精度 Mapping accuracy (%) 大豆 Soybean MLC 857118 3824 56780 917722 90.22 SVM 885153 7611 66296 959060 93.17 RFC 934124 2276 55865 992265 98.32 玉米 Maize MLC 56343 1201206 167854 1425403 97.81 SVM 31759 1145348 87264 1264371 93.27 RFC 2053 1205962 63085 1271100 98.20 其他 Other MLC 36589 23010 1344547 1404146 85.68 SVM 33138 75081 1415621 1523840 90.21 RFC 13873 19802 1450231 1483906 92.42 总计 Total (pixel) 950050 1228040 1569181 3747271 用户精度 User accuracy (%) MLC 93.40 84.27 95.76 SVM 92.29 90.59 92.90 RFC 94.14 94.88 97.73 总体精度 Overall accuracy (%) MLC 90.81 Kappa系数 Kappa coefficient MLC 0.86 SVM 91.96 SVM 0.88 RFC 95.81 RFC 0.94 MLC: 最大似然分类; SVM: 支撑向量机分类; RFC: 随机森林分类。 MLC: maximum likelihood classification; SVM: support vector machine classification; RFC: random forest classification.

表3 增加辅助特征波段后的3种分类方法分类结果混淆矩阵 Table 3 Confusion matrix of three methods based on the stacked image

分析随机森林方法增加辅助特征前后各特征在决策分类中的重要性, 结果如图8所示。原始OLI数据中, 红光波段和2个短波红外波段对于分类的重要性较高, 而海岸、蓝、绿、红光波段的重要性则较低; 在增加了NDVI、NDWI和PCA前4个波段后, 对于分类重要性较高的依然是近红外和2个短波红外波段。可以认为, OLI原始数据中已经包含了较为丰富的作物识别信息, 而通过波段运算、主成分变化生成的新特征波段, 实际上是一种人工的信息挖掘方式, 在原有信息已能被充分利用的情况下, 这些变换并不能提高数据中对于分类有效的信息, 因此对于分类精度提升效果并不明显。可以在表3中发现, 即使是随机森林方法, 精度也没有明显的提升。

图8

Fig. 8

	Figure Option ViewDownloadNew Window
	图8 增加辅助特征前后各特征变量的重要性 a: 增加辅助特征前各波段重要性; b: 增加辅助特征后各波段重要性。Fig. 8 Variable importance before and after adding auxiliary features a: variable importance before adding auxiliary features; b: variable importance after adding auxiliary features.

利用普通工作站统计3类方法的分类时间, 工作站的CPU为i7处理器, 频率为3G Hz, 内存为8G, 统计结果如表4所示。最大似然分类方法的耗费时间最少, 仅约145 s; 支持向量机分类时间耗费最多, 大约为11 000 s; 随机森林的分类时间居中, 约为1800 s。综合考虑各分类方法的作物分类精度及时间, 可以看出, 最大似然分类方法分类时间上具有明显的优势, 但是分类精度相对较低; 而支持向量机方法耗费时间最长, 约为最大似然分类方法的76倍, 随机森林方法的6倍多, 且作物分类识别精度并未有显著提升, 整体上劣势较大; 随机森林方法的作物识别精度最高, 且分类时间相对适中, 对于农情遥感监测作物面积提取具有较大的整体优势。

表4

Table 4

表4(Table 4)

表4 3种分类方式作物分类提取耗费时间 Table 4 Classification time of the three methods 分类方式 Classification method 分类时间 Time cost (s) 最大似然分类MLC 145 支持向量机SVM 11000 随机森林分类RFC 1800

表4 3种分类方式作物分类提取耗费时间 Table 4 Classification time of the three methods

根据随机森林方法的原理, 其决策树可以很方便地并行生成, 且每棵树都不需要剪枝处理, 使其训练时间大幅减少。且随机森林分类方法继承了决策树分类方法分类速度快的优势, 应用训练完成的随机森林模型进行影像分类的耗时也很少。综合这两个因素, 使随机森林方法在获取较高分类精度的同时, 耗时也较少。