基于高分5 号影像的东北典型黑土区土壤分类

孟祥添，鲍依临，刘焕军,3※，张艾明，刘云超，王丹丹

（1. 赤峰学院资源环境与建筑工程学院，赤峰 024000；2. 东北农业大学公共管理与法学院，哈尔滨 150030；3. 中国科学院东北地理与农业生态研究所，长春 130012）

0 引 言

精准监测土壤类型的空间分布有助于更好的制定土地资源利用策略和指导农业生产。目前，遥感技术已被广泛应用于土壤分类的研究中，利用遥感技术可低成本、快速、全面的进行土壤分类，避免了传统土壤分类中必须根据诊断层和土壤的理化性质进行分类的过程。土壤光谱反射率是土壤理化性质和内在结构的综合反映，已有大量研究利用土壤反射光谱数据进行土壤分类[1-3]。

依据土壤反射光谱数据进行土壤分类的研究中，数据来源可大致分为室内测量的高光谱数据、机载高光谱数据、星载多光谱数据和星载高光谱数据。基于室内测量的高光谱数据易于获取，且已经实现了高精度的土壤分类。如Wang 等[4]提取土壤反射光谱的形状特征，实现在土类、土属级别上进行高精度的土壤分类；Zhang 等[5]提取了 148 个土壤样本的光谱特征参数，实现了土壤在土类级别的划分，总精度高达90.54%。然而，基于室内测量的高光谱数据的样点采集在空间是非连续的，无法实现区域性乃至更大尺度的土壤分类及制图，且不利于土壤类型的实时更新与监测。利用机载高光谱数据已实现了对植物、地物类型的高精度分类[6-8]，然而部分机载传感器存在光谱带个数有限、低信噪比（Signal to Noise Ratio，SNR）和覆盖范围有限的问题，导致使用机载高光谱数据在土壤分类研究中存在一定的劣势。星载多光谱遥感数据是发展时间最久、使用范围最广的数据，已广泛应用于土壤分类中[9-10]，然而，多光谱数据波段数量较少，某些与土壤类型相关的细节信息被综合或忽略掉，导致利用星载多光谱数据进行土壤分类的精度不高。如Meng等[11]证明了Landsat8 OLI星载多光谱数据范围在土壤属性相关性较高的波段范围（600～800 nm）内，仅包含其中的2/5 的信息，难以实现高精度的土壤属性预测。利用星载高光谱数据进行土壤分类及制图可弥补室内测量高光谱数据进行土壤性质预测时无法获取连续的土壤性质的缺陷、解决机载高光谱数据覆盖范围有限、星载多光谱数据波段数量少的问题。在高分 5 号（GF-5）全谱段高光谱卫星正式投入使用前，仅NASA EO-1 平台上的Hyperion 传感器和欧洲航天局PROBA 上的CHRIS 光谱成像平台具有足够空间分辨率用于土壤研究[12-13]。但可进行土壤属性预测应用的Hyperion 与CHRIS 成像光谱仪在定量土壤估算中均存在一定的局限性：Hyperion 数据的短波红外区域（2 200 nm 附近）的SNR 较低，无法满足土壤属性遥感预测的要求[14]。CHRIS 数据由于光谱范围（415～1 050 nm）受限，在短波红外区域没有光谱信息[15]。基于上述原因，有必要验证 GF-5 星载高光谱数据进行高精度、大范围的土壤分类及制图的可能性，以期为基于星载高光谱数据进行土壤分类的研究提供参考。

仅考虑土壤光谱特征难以实现高精度的土壤分类及制图，而成土因素的引入可有效地提高土壤分类精度[16-17]。在五大成土因素（气候、生物、母质、地形和时间）中，地形因子可减弱光谱特征相互混淆的现象且适用于大尺度的土壤分类研究[18]。同时，与生物，母质要素相比，地形可直接通过光学影像进行观测，其显著的高程变化幅度可以缩小其他干扰信息的特征。地形通过引起物质、能量的再分配而间接地作用于土壤的形成及演变，使得不同的土类在地形上具有显著的差异，且不同土壤类型的交界位置与空间地形特征紧密关联。

目前，由于可以使用的星载高光谱数据少，且对星载高光谱数据进行包络线去除和主成分分析处理以及加入地形因子（Terrain，TA）是否能够有效的提高土壤分类精度有待验证，因此，本研究利用拜泉县、明水县星载高光谱影像数据，对高光谱数据分别进行包络线去除和主成分分析处理，明确进行包络线去除处理能否提高不同土类之间的差异性；在经过包络线去除和主成分处理后的数据中加入地形因子，并结合随机森林模型构建土类级别的土壤分类模型，以期为利用GF-5 高光谱数据实现高精度土壤分类及制图提供新的方法。

1 材料与方法

1.1 研究区概况

东北黑土区的粮食商品量、调出量均居全国首位，精准的土壤分类及制图结果可为精准农业提供依据。拜泉县、明水县位于黑土带中心，包含黑土区中三大土壤类型，以此为研究区具有代表性。拜泉县、明水县，位于松嫩平原北部、黑龙江西部（124°18′～126°31′E，46°44′～47°55′N），研究区面积约为 6 000 km2，温带大陆性季风气候，夏季多雨温暖，冬季干燥寒冷，年平均气温在2.9 ℃，年平均降水485 mm 左右，主要农作物种植类型为玉米和大豆。平均海拔为240.5 m，地势整体呈北高南低，空间差异显著（图1）。从全国第二次土壤普查数据可知，研究区内成土母质多为黄土母质和沉积物，少部分为冲积物，主要土壤类型共计3 种，分别为黑土、黑钙土和草甸土，占总面积的 95%以上，黑钙土和黑土均是黑土资源的重要组成部分，其有机质含量高，是非常适合粮食生产的土壤类型。不同土壤类型具有不同的地形特征。其中黑土主要分布在高程较高、地表起伏度较低且坡度平缓的区域；黑钙土主要分布在研究区的西南部，整体高程较低；草甸土主要分布在地势较低洼的地区。

东北地区农作物为一年一熟，每年 4 月初期，积雪融化，土壤表面完全裸露出来。当地农民将农田里的农作物残留物（如秸秆与地膜）进行焚烧处理。为了使土壤变得疏松，促进土壤中潜在养分转化为有效养分和作物根系的伸展，在耕种前通常进行翻地处理，由于重复耕作，耕层的土壤性质通常是均匀的，4 月初到5 月初期，拍摄的影像不被农作物残留物和积雪覆盖，未受到混合像元影响，这一时期被称之为“裸土窗口期”[1]。

图1 研究区位置及海拔Fig.1 Location and elevation of the study area

1.2 数据来源

GF-5 高光谱数据从黑龙江省高分卫星数据中心申请获取。本研究选取了拜泉县、明水县共 4 幅符合“裸土窗口期”的GF-5 号高光谱影像，影像的具体拍摄时间为4月 18 日，无云覆盖，获取的前一周无雨，且土壤的表面没有被大片植被和积雪所覆盖。全国第二次土壤普查数据从黑龙江省农垦科学院科技情报研究所申请获取，该土壤普查结果是根据土壤发生学理论对土壤进行划分。从http://glovis.usgs.gov/网站获取了研究区 30 m 空间分辨率数字高程模型（Digital elevation model，DEM）数据。

1.3 数据处理

1.3.1 数据预处理

利用ENVI5.1 软件对GF-5 高光谱影像进行辐射定标和大气校正处理，随后使用Arcgis10.1 按照县域耕地范围将影像进行裁剪，用于后期处理。提取各采样点所对应的反射光谱曲线，由于传感器在400～430、2 450～2 500 nm波谱范围具有较低的SNR，在900～1 050 nm 受到传感器自身缺陷的影响[12]、且在1 350～1 451、1 771～1 982 nm受到大气中水蒸气吸收的影响，造成光谱数据不连续。因此，本文选取 430～900、1 050～1 350、1 451～1 771、1 982～2 450 nm 作为本次研究的光谱区间。

1.3.2 包络线去除

包络线去除也称连续统去除法，去除土壤中特定物质化学键内电子跃迁引起的特征吸收，分离出土壤本身的吸收特征。经包络线去除处理后，反射率（Original Reflectance，OR）被归一化到0～1 之间[19]，原始光谱中“峰值点”为 1，其余点小于 1，从而得到去包络线数据（Continuum Removal，CR）。该处理可以有效地突出不同土壤的光谱曲线的吸收和反射光谱特征，最大限度地增加不同土壤类别之间的差异。

1.3.3 主成分分析

主成分分析（Principal Component Analysis，PCA）是一种常用的数据分析方法，PCA 通过线性变换将高维度的原始数据转换成低维度且各维度线性无关的向量，可用于提取数据的主要特征分量，常用于数据降维。本文分别对OR 和CR 进行PCA 处理，提取累计贡献率大于85%主成分信息，分别记作OR-PCA 和CR-PCA 作为土壤分类模型的输入量。

1.3.4 地形因子获取

地形（Terrain，TA）是五大成土因素之一，其中海拔高度是基础的地形因子，衍生因子如坡向，曲率及地表起伏度共计 5 种地形因子参与土壤分类。地形因子均提取自30 m 空间分辨率的DEM 数据。分别在OR-PCA和CR-PCA 的基础上加入这5 种地形因子，得到输入量OR-PCA-TA 和 CR-PCA-TA。

1.4 土壤可分性分析方法

土壤的可分性代表了该土壤与整个土壤集中的其他土壤的差异性，可分性越大，该土壤与其他土壤的差异性越大，原则上更易实现较高的分类精度。可分性（S）的计算公式如下：

式中Vi为土类内差异性；Vb为土类间差异性；STmn为土类m和土类n之间的标准差；Mmn为土类m和土类n之间的均值；STm为土类m的标准差；Mm为土类m的平均值。

1.5 随机森林分类模型

随机森林（Random Forest，RF）是由多个决策树组成的集成模型。该模型是美国统计学家Breiman[20]提出的一种基于树的集成学习模型，训练过程采用随机采样方式，增加了模型的泛化能力；目前被广泛地应用于解决分类和回归问题[20-22]。本文在建立RF 模型时，通过观察袋外误差选择最佳回归树的数量（ntree）和分裂节点数（mtry），从而建立最优 RF 预测模型。在训练模型时，确定了不同的参数值：ntree 为500，mtry 为输入量个数的平方根[23]。

1.6 样点选取及精度验证

根据全国第二次土壤普查图选取训练样本和验证样本（图 2），样点采集原则如下：1）同一土壤类型在不同区域的反射光谱特性具有差异，因此，同一土壤类型需要在不同的区域分别采样；2）为了避免不同土壤类型的交界处模糊而造成分类误差，尽量在不同土壤类型的中心区域选取样点；3）不同土壤类型采样点数量的比例要与研究区中不同土壤类型之间面积的比例大致相同，以保证样点分布均匀[24]。选择样本后将所有的样本进行60 m（2 个像元）的缓冲区处理，最终得到训练样本2 269个（黑土863 个、黑钙土461 个、草甸土945 个）、验证样本 1 169 个（黑土 492 个、黑钙土 294 个、草甸土383 个）。使用混淆矩阵验证分类的准确性，并使用总精度和Kappa 系数代表混淆矩阵中的结果，具体计算公式如下：

图2 研究区采样点F i g.2 S a m p l e s i t e s o f s t u d y a r e a

式中k为混淆矩阵中列的数量；Pii为混淆矩阵中第i行第i列的像元数，表示正确分类的个数；Pi+和P+i分别第i行和第i列总像元个数；N代表验证像元的总个数。

2 结果与分析

2.1 不同土壤类型反射光谱特征分析

不同类型的土壤具有不同的反射光谱特征[25]。图 3a为不同土壤类型的反射率曲线，不同土壤类型的原始反射率曲线的形状差异较小。其中，黑钙土的反射率最高，这是由于黑钙土的表面有一层钙化层。与黑钙土相比，黑土的反射率稍低，这主要由于其保水保肥能力强，有机质含量较高，颜色比较暗。由于草甸土主要分布在地势较低的地方，其水分含量较高，从而导致其反射率最低。图3b 代表不同土壤类型的去包络线，草甸土在400～900 nm 的去包络线值最低，该区间对应的面积也是最大的。黑钙土在整个光谱区间内的去包络线值都比较高，特别是在1 100～1 350 nm 区间，通过去包络线处理后，不同土壤类型的光谱特征差异增加，有助于土壤分类精度的提升。

2.2 可分性描述

包络线去除处理可明显的增加不同土壤类型之间的差异，从而提高不同土类的可分性，与OR、OR-PCA相比，CR、CR-PCA 的可分性的绝对值分别提高了0.04、0.09。PCA 处理有效的提高了 OR、CR 的可分性，与OR、CR 相比，OR-PCA、CR-PCA 可分性的绝对值分别提高了0.04、0.09（图4a）。图4b 显示了不同土壤类型中，5 种地形因子的可分性大小。在不同土壤类型中，5 种地形因子的可分性由高到低的顺序是不同的。与OR、CR 的可分性相比，地形因子整体的可分性更高。

图3 训练样本不同土类的原始反射率曲线和去包络线Fig.3 Original reflectance curves and continuum removal curves of different soil types in training samples

2.3 土壤分类结果

2.3.1 未经PCA 的土壤分类结果

表 1 为基于星载高光谱遥感影像数据不同输入量的土壤分类结果。以OR 为输入量的不同土壤类型分类精度由高到低依次是黑钙土、草甸土、黑土。在黑土中，有213 个土壤样本被错分到草甸土中，在草甸土中共有106个样本被误分到黑土中，这是由于黑土和草甸土的土壤反射率特征差异较小导致的，而黑钙土的反射率特征与黑土和草甸土的差异较大，因此只有少量的土壤样本被误分至黑土或草甸土中。

CR 作为输入量时，不同类型土壤的分类精度均高于OR 作为输入量时的分类精度。与OR 作为输入量相比，CR 作为输入量的分类总精度提高了5.48%，Kappa系数提高了0.12，黑土、黑钙土、草甸土的正确分类数分别提高了20、17 和27 个。这是由于经过包络线去除处理，有效的增加了不同土壤类型之间的差异性，特别是黑土与草甸土之间的差异。

图4 不同输入量、不同地形因子的可分性Fig.4 Separability of different inputs and terrain factors

表1 基于不同输入量的土壤分类结果Table1 Soil classification results based on different inputs

2.3.2 PCA 后的土壤分类结果

对OR 进行PCA 处理，提取了3 个主成分信息（累计贡献率为96.28%），以OR-PCA 作为输入量时，不同土壤类型分类精度由高到低的顺序是草甸土、黑钙土、黑土。与OR 作为输入量相比，土壤分类的总精度提高了1.71%，Kappa 系数提高了0.02，黑土、草甸土的分类精度提高，黑钙土的分类精度降低。

对CR 进行PCA 处理，提取了7 个主成分信息（贡献率为85.95%），在CR-PCA 中，不同土壤类型的分类精度依次是黑钙土、草甸土、黑土。与CR 作为输入量相比，土壤分类总精度提高了 3.67%，Kappa 系数提高了0.02。与OR-PCA 作为输入量相比，土壤分类总精度提高了7.54%，Kappa 系数提高了0.12，黑土、黑钙土、草甸土的正确分类数分别提高了21、39 和27 个，较大程度地提高了土壤分类精度。

2.3.3 引入地形因子后土壤分类结果

从表 1 中可以看出，加入地形因子后各土壤类型的分类精度显著提升。在OR-PCA-TA 中，各土类的分类精度都提高了10%以上，与OR-PCA 作为输入量相比，土壤分类的总精度提高了14.12%，Kappa 系数提高了0.23。在 CR-PCA-TA 中，黑钙土和草甸土的分类精度达到了80%以上，与 CR-PCA 相比，黑土、黑钙土、草甸土的正确分类数分别提高了81、30 和41 个。

2.4 土壤类型结果

在 OR、CR 作为输入量时，制图结果中（图 5）出现了明显的“椒盐”现象。经过PCA 处理后，制图结果中的“椒盐”现象减弱，土壤类型中心区域的空间分布更符合第二次土壤普查结果。以 OR-PCA-TA、CR-PCA-TA 作为输入量的制图中，不同土壤类型的边界划分的更加细致，研究区中部及东北部大部分是黑土，空间分布成条带状，黑钙土主要分布在研究区的西南部，草甸土穿插分布于各类土之间。这主要由于黑土主要分布在地势较为平坦处，本研究区大部分地区地势平坦，草甸土主要分布在低洼地势中。

图5 基于不同输入量的土壤制图结果Fig.5 Mapping results of soil based on different inputs

3 讨 论

过去利用遥感手段进行土壤分类的研究中，包络线去除方法主要对室内测量的高光谱数据进行处理[4,5,26]，而对星载高光谱数据进行处理的研究有限。本文利用GF-5 星载高光谱数据进行土壤分类，证明了CR 作为输入量可明显提高土壤分类精度，这是由于包络线去除处理可明显的增加不同土壤类型之间的差异，从而提高不同土类的可分性（图4a）。

经过PCA 处理可有效的降低制图结果中的“椒盐”现象，提高了土壤分类精度，尤其是黑土和草甸土的正确分类数得到了明显的提升，证明了PCA 可以增加黑土和草甸土之间的可分性；同时PCA 能够降低高光谱数据的维度，极大减少了分类模型的计算时间，提高模型的运算效率，且降低了遥感影像中像元的错分，减弱了不同影像获取时间对于土壤分类结果的影响，使土壤类型的变化更加圆滑。试验证实了在星载高光谱数据中，PCA 仍然是降低数据冗余性并提高土壤分类精度的有效方法。

地形作为五大成土要素之一，对母质起着重新分配的作用，同时，支配地表径流，通过改变局地小气候影响环境。以往的研究表明，地形因子的选取通常与研究区范围有关[27-30]，一般而言，研究区尺度越大，地表空间差异越明显。同时，本研究发现不同土壤类型的可分性最高的地形因子具有差异（图 4b），在黑土中，海拔高度的可分性最高，这是由于黑土主要分布于研究区北部，其中分布着大量的草甸土，而草甸土地势低。黑钙土位于研究区中部地带，地势相对于黑土较低，且主要分布在地表变化较大的地区，因此地表起伏度在黑钙土中的可分性较高。草甸土发育于地势低平区域，如洼地，沟壑等，在研究区内较为分散，土壤边界位置高程差异明显，因此海拔高度与地表起伏度对提升草甸土分类精度具有重要意义。

4 结 论

本文利用东北黑土区拜泉县、明水县高分 5 号星载高光谱影像数据，对原始反射率进行包络线去除处理得到去包络线，再对原始反射率和去包络线进行主成分分析（Principal Component Analysis，PCA），并引入地形因子，建立随机森林土壤高光谱遥感分类模型，进行土壤研究并制图，得到如下结论：

1）经过包络线去除处理可有效地增加不同土壤类型间的可分性，从而提高土壤分类的精度。与原始反射率相比，经过包络线去除处理的总精度提高了 5.48%，Kappa 系数提高了0.12。

2）利用 PCA 可以提高土壤分类精度，尤其是黑土和草甸土的正确分类数的提高，同时，降低高光谱数据的维度，减少数据冗余，提高模型的分类效率。在制图结果中，经过PCA 处理后可明显降低了像元的错分，使土壤类型的变化更趋于实际。

3）引入地形因子，有效的提高了土壤分类精度，且不同土壤类型可分性最高的地形因子是不同的，其中以经过PCA 处理后的去包络线结合地形因子作为输入量的精度最高，分类精度为81.61%，Kappa 系数为0.72，实现了区域尺度的高精度土壤分类及制图。且研究结果可为进行大范围土壤分类及制图提供参考。