基于航空高光谱反演的黑土地质量综合评价研究

赵宁博， 杨佳佳， 赵英俊， 秦凯*， 杨越超， 李明

(1.核工业北京地质研究院，遥感信息与图像分析技术国家级重点实验室，北京 100029；2.中国地质调查局沈阳地质调查中心，沈阳 110034)

土壤质量通常指土壤在生态系统中保持生物生产力、维持环境质量及促进动植物健康的能力[1]。土壤质量评价主要从土壤肥力质量、环境质量和健康质量三个维度[2-3]开展，评价指标涵盖土壤物理、化学和生物学指标。土壤质量评价结果可以直接指导土地利用规划、农业生产、土壤环境保护工作，同时也是污染土壤治理修复、土壤污染生态效应评价及地球化学灾害预测研究的基础[4]。

我国东北地区的黑土地性状好、肥力高，非常适合农作物生长，是我国重要的粮食生产基地。但近年来黑土地腐殖层变薄、肥力下降、水土流失等现象已经引起了广泛关注[5-7]，黑土地的退化将直接影响农业生态系统和农业可持续发展，因此必须科学及时地开展土壤质量评价工作，为黑土地的保护与可持续利用提供有效决策依据。

目前，我国黑土地的土壤质量评价主要采用地球化学调查手段，局限性表现为地面采样耗费较多人力物力，样品分析工作周期较长，在大比例尺调查时局限性表现的更为明显。遥感技术凭借时效性强、数据覆盖面广等优势，能够对地球化学调查进行有效补充，逐渐在评价工作中发挥作用。根据遥感数据的特点，在土壤质量评价中主要有两方面应用，一是利用多光谱影像提取植被指数(normalized differential vegetation index，NDVI)[8]、土壤退化指数(ratio vegetation index，RVI)[9]、土壤水分指数(difference vegetation index，DVI)[9]、地形因子[10]、土地利用[11]等信息；二是利用高光谱数据提取土壤有机质[12-15]、养分元素[16]、重金属元素[17-18]、含盐量[19]、质地[20]等指标。

相比于多光谱遥感，高光谱遥感提取土壤理化性质的能力更强，能够获得更为全面的土壤质量指标数据。目前，土壤质量调查中主要采用的是地面高光谱技术，数据为散点形式且工作周期较长，航空高光谱技术改善了地面调查方式下数据获取周期长及覆盖面不足的劣势，且具有“图谱合一”的特征，提高了在土壤质量评价中的实际应用能力。本文以黑龙江省海伦典型黑土区为例，利用航空高光谱数据从土壤肥力、土壤环境和农作物长势几个维度提取相关信息，构建土壤质量综合评价模型，初步探讨航空高光谱技术在黑土地质量评价中的应用能力和效果，为田块尺度的黑土地质量快速评价提供技术支撑。

1 材料与方法

1.1 研究区概况

研究区位于我国黑龙江省海伦地区，坐标126°29′10″—127°4′40″E、47°31′10″—47°35′40″N。海伦地处松嫩平原东北端，小兴安岭西麓，地形为丘陵、漫岗，属中温带大陆性气候。

海伦地区位于我国东北黑土区的中心区域，主要土壤类型为黑土、暗棕壤和草甸土等，是黑土地研究领域中的代表性区域。该地区已经开展了1∶5万多目标地球化学研究，航空高光谱与地球化学的有机结合能够为黑土地关键带的综合建模提供支撑，便于开展黑土质量评价工作。

1.2 数据获取与预处理

1.2.1航空高光谱数据 设备采用CASI-1500和SASI-600高光谱成像系统，由加拿大Itres公司生产。CASI-1500传感器谱段范围为380～1 050 nm，空间分辨率1.5 m，波段数72个；SASI-600传感器谱段范围为950～2 450 nm，空间分辨率为3.75 m，波段数100个。为了保证航空高光谱数据质量，挑选晴朗无云的天气飞行，飞行时间设定为当天10：00—14：00。数据分两期获取：第一期次实际获取日期为2018年5月29日，地面处于裸土期，用于反演土壤质量参数；第二期次实际获取时间为2018年7月31日，为夏季农作物成熟期，用于反演农作物质量参数。

数据预处理主要包括辐射定标、几何校正、大气纠正和光谱重建。航空成像光谱测量系统配备了辐射定标软件RCX(radiometric calibration xpress) 9.3.5.1和几何校正软件Geocor 3.0。首先利用RCX软件进行辐射定标，然后利用Geocor软件对原始GPS数据和基站数据进行差分处理，提取GPS时间并处理传感器姿态数据，对传感器姿态数据与GPS定位数据进行时间同步与集成，最后利用几何校正信息对辐射定标后的数据进行几何校正。然后在ENVI5.3软件中继续进行几何精校正处理，利用高分2号影像作为参考数据，利用控制点校正的方法完成几何精校正，进一步提升数据的几何精度。

在此基础上，利用ENVI5.3软件中大气辐射传输模型和地空回归方法对航空高光谱遥感数据重建光谱。首先在ENVI的FLAASH模块中利用大气辐射传输模型进行校正，然后利用地面黑、白两种定标布的光谱测量数据在ENVI的empirical line模块中进行地空回归校正，进一步消除因辐射定标、波段间相对定标、波段配准、大气参数选取等误差因素而造成的光谱误差，最终获得地物光谱反射率。

为了更全面地提取研究指标的光谱特征，利用ENVI5.3软件扩展工具中的Image Derivate模块反射率数据进行一阶导数变换，Continuum Removal模块进行去连续统变换。

本次数据反演及质量评价仅限于研究区内的旱田范围，因此在高光谱影像中提取旱田范围用于后续数据处理，地物分类提取方法采用ENVI软件扩展工具中的随机森林法。首先在影像内建立感兴趣区(region of interest，ROI)，ROI均匀分布于全区范围，保证覆盖到各种地物类型，然后利用随机森林模型对选取的ROI进行训练，获得研究区的地物分类数据，最后利用ENVI软件中的Majority Analysis工具分类数据进行后处理，去除细小斑块。

1.2.2地面数据 在航空数据获取时同步开展地面调查，调查点包括38个土壤调查点及22个农作物调查点。地面获取的数据指标类型即为参与土壤质量评价的指标类型，地面数据用于后续的航空高光谱反演建模及验证。指标选取的依据包括：①以土地质量地球化学评价规范[21]为基础，选取通用的土壤养分和环境指标，保证评价结果的科学性和通用性。②突出航空高光谱的优势，结合评价目标拓展土壤理化性质和农作物长势方面的指标。

根据上述指标选取依据，地面共获取了包括土壤肥力、土壤环境和农作物长势三个类别共15项指标： ①土壤肥力指标。包括有机质(SOM)、全氮(N)、全磷(P)、全钾(K)、硒(Se)、锗(Ge)、阳离子交换量(cation exchange capacit，CEC)；②土壤环境指标。土壤环境指标包括重金属污染和土壤退化两类。重金属元素参照规范[22]选取了砷(As)、铬(Cr)、镉(Cd)、汞(Hg)、铅(Pb)五种元素。土壤退化指标选取了全盐量，代表土壤盐碱化程度。③农作物长势指标。包括叶绿素含量和叶面积指数，在第二期次农作物航空高光谱调查时同步测量。

1.2.3土壤指标化学分析 SOM采用硫酸亚铁铵容量法[23]滴定，仪器为酸式滴定管；采用凯氏定氮法[24]测定N含量，仪器为上海市欧公司生产的全自动凯氏定氮仪；采用X射线荧光光谱法[25]测量P和K含量，仪器为荷兰帕纳科公司生产的Axios-mAX 波长色散型X射线荧光光谱仪；采用电感耦合等离子体发射光谱法[26]测量Se、Ge、Hg、As、Pb、Cd、Cr含量，仪器为美国赛默飞世尔公司生产的ELEMENT XR 等离子体质谱仪；采用乙酸铵法[27]测定CEC，仪器为上海赫田公司生产的电动离心机；采用浙江托普仪器公司生产的SPAD502PLUS叶绿素仪测定叶绿素含量；采用英国Delta-T公司生产的SUNSCAN冠层分析仪测定叶面积指数。

1.3 土壤质量参数反演

本次评价所采用的各类土壤质量参数数据均为航空高光谱数据反演所得。航空数据的光谱分辨率较高，通过不同土壤质量参数的光谱特征对其含量进行反演。将航空数据与地面采样点进行空间叠加，在航空数据中以采样点位置为中心，周围3×3个像元范围的光谱平均后作为该点的光谱数据。

采用偏最小二乘法建模。在建模时，特征波段选择直接关系着模型的精度和稳定性，由于土壤成分和性质较为复杂，针对同一指标，土壤类型或区域不同时相应的特征波段也不尽相同。在选择特征波段时，首先利用全部波段参与建模，然后根据模型系数曲线选择峰值区域的波段再次建模，不断优化调整特征波段的，最终根据建模和验证精度选择最优模型。

将总样本数量的2/3用于建模，1/3用于验证，根据样本的空间分布均匀挑选建模和验证样本，以减少模型误差。

1.4 土壤质量综合评价

以土壤及农作物多项质量参数的航空高光谱反演数据为基础，通过建立评价模型最终获得研究区的土壤质量综合评价指数。选取层次分析法(analytical hierarchy process)作为评价方法将土壤质量评价系统结构化，形成递阶层次结构模型，然后对每一层次中的各因素进行较客观的逐对比较和判断，最后通过排序结果来分析和解决问题。层次分析法所采用的软件为YAAHP10.1。

1.4.1评价层次结构 根据选取的评价指标建立层次结构(图1)，分为目标层、中间层和指标层三个层次。目标层为土壤质量，即评价工作的最终目标。中间层包括土壤肥力、土壤环境和农作物长势三个大类，土壤肥力又包括大量元素和微量有益元素类别；土壤环境包括重金属污染和土壤退化两类。指标层即上文中选取的SOM、N、P、K等各项土壤养分、环境、农作物长势类别中的具体指标。

图1 评价模型层次结构

1.4.2权重计算 权重计算是建立评价模型的重要步骤之一。根据咨询专家、参考现有规范和统计资料对判断矩阵中各指标的重要性进行两两比较，按照9标度方法[28]进行打分，构建出比较判断矩阵。根据建立的判断矩阵，采用和积法分别计算出最大特征值所对应的特征向量，所求特征向量即为相对权重值。权重计算完成后需要计算一致性比率(consistent ratio，CR)，当CR=0时，表明权重分配具有完全一致性；当CR<0.1时,具有满意一致性,说明权重分配合理；否则就要调整判断矩阵，直到具有满意的一致性为止。本次判断矩阵的CR计算结果为0.017 6，通过一致性检验。最终权重计算结果见表1。

表1 模型权重计算结果

1.4.3评价单元的确定 通常根据评价数据比例尺的不同，评价单元可以按照规则网格、土地利用现状调查图斑、实际地块等不同方式划分，本研究凭借航空数据的高空间分辨率优势，评价单元采用最高精度划分方式，即研究区的实际地块。地块数据以研究区第二次全国土地调查数据为基础，并与航空高光谱影像进行叠加后检查地块空间划分方式的变化，经检查后整体上地块划分无明显变化，并对个别地块发生改变的区域进行了修改完善。

1.4.4指标评分 通过制定打分标准对每种指标进行评分，分值区间为0～1分。评分标准制定以相关规范中的标准为基础，并考虑不同指标特征进行科学划分，具体评分标准如下。

①土壤肥力指标评价标准。肥力指标中SOM、N、P、K的评价标准参照规范[21]中的规定，评级由高至低包括五类等级：一等(1.0分)、二等(0.8分)，三等(0.6分)、四等(0.4分)、五等(0.2分)。需要说明的是，P和K采用的是全量，而不是有效态含量。速效磷和速效钾是能够被植物吸收的部分，与土壤肥力的关系比全量更为密切，但是通常情况下有效态含量占全量的比重很小，过低的含量会大大提高其光谱特征的获取难度，从而影响反演精度。本研究的重点是在保证反演精度的前提下构建综合评价体系，后续将进一步开展对速效磷和速效钾等指标的高光谱反演研究。

有益元素Se和Ge的评价标准采用隶属度函数模型。隶属度函数模型针对不同类型评价目标包分为戒上型、戒下型或峰值型三类。Se和Ge均具有生物学的两面性，在一定含量范围内对人体和动植物生长有益，过量则会造成损害，基于该特性，如果研究区指标含量中存在明显过量级别时宜采用峰值型函数，如果没有过量级别存在则采用戒上型函数。Se和Ge元素的过量临界值标准分别采用3.0和5.8 mg·kg-1[29]，根据研究区统计结果，Se最大值为0.43 mg·kg-1，Ge最大值为1.69 mg·kg-1，均低于过量的临界值，因此隶属度函数模型采用戒上型，通过隶属度函数计算得到Se和Ge的评分。

目前CEC较为常用的评价标准分为三个等级，CEC≥20 cmol·kg-1代表保肥能力强；1017～20)，0.8分；三等(>13～17)，0.6分；四等(>10～13)，0.4分；五等(≤10)，0.2分。

②土壤环境指标评价标准。土壤重金属元素As、Cr、Cd、Hg、Pb的分级标准同样参照规范[21]。

土壤全盐量数值通过电导率法测定，因此参照相关电导率数值的分级标准[30]制定评分规则。

③农作物长势指标评价标准

农作物长势指标包括叶绿素含量和叶面积指数，这两种指标含量均与农作物长势呈正相关，因此采用戒上型隶属度函数模型进行计算评价。

1.4.5土壤质量分级 评价单元赋值包括两个步骤：①在ArcGIS软件中将评价所采用的矢量点文件与指标反演的栅格文件进行空间分析，将矢量点所在空间位置对应的栅格数值赋予该点；②将矢量点文件与地块图斑进行空间分析，同一地块中所有矢量点数据取平均值作为该地块的数值。

将所有评价指标的数据赋予地块图斑后，依照上文所述的相应评价标准进行单指标得分计算，然后依照如下计算方式获得综合评价得分。

式中，F综为综合评价得分，Fi为第i个指标评价得分，Wi为第i个指标权重。

综合评价得分的分值区间为0～1.0，共分为五个等级：一等(优质)土壤，>0.9～1.0；二等(良好)土壤，>0.7～0.9；三等(中等)土壤，>0.5～0.7；四等(差等)土壤，>0.3～0.5；五等(劣等)土壤，≤0.3。

2 结果与分析

2.1 模型关键参数分析

2.1.1光谱数学变换形式对建模的影响 分别采用光谱反射率、反射率一阶导数和反射率去连续统变换数据作为模型自变量参与计算，各指标对应的模型决定系数(R2)见表2。针对每一种指标，反射率一阶导数变换对应的建模和验证精度均远低于反射率相对应的精度；基于反射率去连续统建立的模型中，N、Ge、Cr和Hg的建模精度优于反射率对应的模型精度，但是所有指标的验证精度均低于反射率对应的精度。

表2 基于不同光谱形式建立的模型R2

一阶导数和去连续统变换均有提升微弱光谱特征的作用，但此次建模效果并不如原始光谱反射率，尤其是基于一阶导数变换的模型精度最低。原因是航空光谱受到大气传输等因素的影响，数据信噪比要低于地面光谱，一阶导数和去连续统变换在增强光谱特征的同时也提高了噪声的影响，因此最终选用原始反射率数据进行建模。

2.1.2特征波段分析 首先利用全部波段参与运算，然后基于回归系数曲线挑选具有代表性的峰值波段作为特征波段再次运算，根据模型精度进行波段的优化调整，直至模型达到最优。以SOM为例，根据图2挑选了570、943、1 085、1 175、1 505、1 670、1 775、2 105、2 225和2 390 nm作为特征波段，利用上述特征波段再次建模并观察波段回归系数，最终将570 nm剔除波段组合后模型精度达到最优。根据上述方法对各指标进行特征波段选择与模型优化，各指标的特征波段及最终模型精度见表3，R2均达到较高水平，说明优化后的模型精度得到提升。

图2 SOM模型回归系数

2.2 反演模型精度分析

评价指标反演的可靠性是后续土壤质量评价的基础。从表3可以看出，SOM的建模R2为0.813，验证R2为0.726，两项均为所有指标中最高；N、P、K、Se、CEC的建模R2均超过了0.7，验证R2均超过了0.6；养分指标中Ge元素的建模R2稍低，为0.667；As、Cr、Cd、Hg、Pb几种重金属元素及全盐量的建模R2处于0.6～0.7区间，验证R2处于0.5～0.6区间。叶绿素和叶面积指数两种农作物长势指标的建模R2均超过了0.7，验证R2均超过了0.6。

表3 各指标特征波段与模型最终精度

从全部指标反演精度的横向对比来看，整体上土壤养分和农作物长势指标的模型精度要稍高于重金属元素，一方面因为重金属元素不像有机质等成分具有明显的光谱特征，另一方面研究区内重金属元素含量较低，一定程度上影响了通过土壤矿物等组分进行间接反演的效果。对比各指标的建模R2和验证R2，没有出现验证R2显著降低的情况，说明反演没有明显的过拟合现象。

2.3 反演结果分析

通过反演结果(图3)观察研究区各指标的空间分布状态。可以看出，SOM高值区主要集中于爱民乡，低值区主要位于海北镇，其他区域处于中间状态；N和SOM的空间分布规律接近，这是由于土壤中的氮素绝大多数是贮藏在有机质中的有机态含氮化合物，其次是被黏土矿物吸附的铵态氮、硝态氮和亚硝态氮，因此两者关系密切；P高值区主要分布于研究区南部，向荣乡-长发镇-东林乡一带；K的分布较为均匀，空间变异程度较低。Se元素高值区主要分布于长发镇、向荣乡和爱民乡；Ge元素分布状态较为均匀，在海北镇稍高；阳离子交换量的分布形态与有机质类似。

图3 部分指标含量反演分布

重金属元素As、Cr、Cd、Hg、Pb的含量均低于规范[22]规定的风险管控值，而在空间分布的相对趋势上，Hg元素的分布较为均匀，Pb及其他三种元素的相对高值区均位于研究区南部。叶绿素和叶面积指数的分布规律较为接近，整体上分布状态较为均匀，高值区相对集中于长发镇和东林乡。

2.4 土壤质量综合评价结果分析

评价范围为研究区的旱田范围，将所研究的地块按照等级进行分级显示，其余区域(水田、居民地等)以遥感影像作为底图，最终评价结果如图4所示。结果统计显示，研究区地块等级均在二等及以上，可见该区土壤质量整体较高。其中一等(优质)地块面积占全区的58.63%，二等(良好)地块面积占全区的41.37%。在空间分布上，两个等级地块分布较为均匀，呈交错分布状态。

图4 研究区土地质量综合评价结果

通过单指标评价结果进一步分析土壤综合等级的影响因素。根据单指标评分，SOM、N、K和CEC四种养分指标评级均以一等为主，其中SOM一等地块占比98.23%，N一等地块占比96.34%，P一等地块占比100%，CEC一等地块占比98.65%；P和Se以二等地块为主，分别占比95.31%和97.45%。

五种重金属元素及全盐量的评级中一等地块均占比100%，表明研究区土壤重金属污染及盐碱化的威胁较小，土壤环境质量较高。农作物数据由于是单期影像，因此主要评价空间上不同地块之间的长势差异性，评价结果显示叶绿素一等地块面积占比53.26%，二等地块占比46.74%，叶面积指数一等地块面积占比46.35%，二等地块占比53.65%。

根据地面数据验证航空高光谱评价结果的准确性。由于评价单元达到了地块级，所以挑选与评价单元相匹配的最大比例尺的地面数据，地面验证数据比例尺达到1∶10 000，分布区域主要位于长发镇，面积占研究区总面积的21.3%，该部分地面数据并未参与前期的建模及土地质量评价，因此可以保证验证的客观性。地面数据的综合评价方法与航空数据评价方法完全一致，以地块为单元对比两类数据评价等级的一致性，二者一致率达到97.6%，表明航空高光谱评价结果可靠。

3 讨论

本研究基于航空高光谱技术构建土壤质量综合评价模型，将航空高光谱与地球化学调查进行有机结合，发挥了航空高光谱技术空间及光谱分辨率高、时效性强等优势。土地质量地球化学评价工作在大尺度(地块级)调查中的比例尺通常需达到1∶10 000及以上，而本次航空高光谱数据空间分辨率为3.75 m，具备数据“全覆盖”的特征，数据精度有较大优势。土壤质量评价中单纯依靠地面调查工作将会耗费大量的人力、物力，工期较长，而航空高光谱调查能够有效缩减工作周期，一次飞行即可获得全区数据，能够有效减弱由于不同数据获取周期引起的系统误差。同时，所处研究区此前开展的地面调查数据能够用于综合研究，提升数据的利用率。

数据反演质量是评价的基础，目前常用的反演建模方法[31-37]包括多元逐步线性回归、支持向量机、偏最小二乘法、人工神经网络、随机森林等多种方法。本研究选择建模方法的原则是方法应用效果较为成熟并且样本需求量较少，以便于后续土壤航空高光谱调查工作的工程性推广，最终选用了偏最小二乘法。偏最小二乘法结合了建模类型的预测分析方法和非模型式的数据内侧分析方法，在一个算法下同时进行数据结构简化(主成分分析)、回归建模(多元线性回归)以及两组变量之间的相关性分析(单相关分析)，简化了多维复杂数据的分析难度。此次评价指标中有机质、养分和农作物指标的建模总体上获得了较好的反演效果，各项指标的验证精度与建模精度对比没有明显的下降，表明建模没有出现明显的过拟合现象，为评价工作提供了较可靠的数据基础。

利用航空高光谱反演获取了土壤肥力、环境与农作物长势等多项土壤质量指标，基于层次分析法构建了土壤质量综合评价模型。综合评价模型考虑了影响土壤质量的多方面因素，在为地方农业生产及规划提供决策依据时兼顾了全面性与简洁性，使整个调查系统更为完善。

目前基于航空高光谱技术开展土壤质量的综合评价尚处于起步阶段，后续工作可对以下方面进行改进和完善：①土地质量参数反演方法。不同质量参数的光谱特征不尽相同，机理较为复杂，可借助深度学习方法充分研究光谱机理，进一步提升反演精度；②评价指标补充完善。本研究筛选了一系列指标进行综合建模，但还需要充分挖掘航空高光谱对土壤质量参数的反演能力，扩充及优化评价指标和模型结构，提升评价的综合性和客观性。