研究区位于莱州湾南岸滨海平原(图1), 在37° 30'N— 36° 37'N 和118° 43'E— 119° 42'E范围内共采集111个样点, 全区以平原为主, 80%以上为耕地, 以种植玉米和小麦为主, 北部存在轻微盐渍化主要是种植棉花, 土壤类型为潮土和盐化潮土为主。 样点的设计原则是基于交通可达性在道路两侧农田内根据区域变异性和先验知识抽取。 采集时间为2016年6月, 在农作物收割和播种的间隔时间内采样, 可以减小环境因素对实测土壤反射光谱的影响。 实验室中土壤有机质含量的测定采用高温外热重铬酸钾氧化-容量法(VOL), 其统计描述特征如表1所示。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 研究区位置及采样点示意图Fig.1 Geographical position map of study area and locations of sampling points

表1

Table 1

表1(Table 1)

表1 土壤有机质含量统计描述特征(%) Table 1 Characteristics of the soil sample sets(%) 样本布置 样点数 均值 最大值 最小值 中位值 标准差 总样本 111 0.73 3.82 0.09 0.81 0.48 校准样本 89 0.73 3.82 0.09 0.69 0.48 验证样本 22 0.77 1.36 0.15 0.85 0.31

表1 土壤有机质含量统计描述特征(%) Table 1 Characteristics of the soil sample sets(%)

111个自然状态下土壤样点反射光谱使用FieldSpec HH(Analytical Spectral Devices, USA)在2016年6月份上午11点— 下午2点之间且采集土壤样品前进行测定。 FieldSpec HH提供Vis(325~780 nm)和NIR(781~1 075 nm)两种探测方式的测定, 光谱分辨率为1 nm。 激光探头保持距地表10 cm且连续获取10条光谱。 在测定前和测定50条光谱后, 进行白板校正, 以减少噪声。 每个采样点五处共得到50条光谱, 作算术平均后得到土样实际的反射光谱数据。

对原始光谱进行平滑和一阶倒数变换处理可以突出和增强反射光谱所包含的有机质信息^[16]。 首先, 将每个样点获得的50条反射光谱平均处理, 处理软件采用ASD自带的ViewSpecPro 4.07中进行; 然后, 在Envi 5.1(Exelis Visual Information Solutions, Boulder, USA)中对原始光谱进行平滑去噪处理; 最后反射光谱经过一阶导数变换后与有机质进行相关性计算(SPSS 19.0, IBM, USA), 根据波段间隔(Δ =50 nm)和相关性的大小(r> 0.4)提取敏感波段。

选取合适的光谱特征可以有效的提高预测精度, 这在复杂的地表系统下实现是极其困难的, 但是可以基于主成分的变换结果引入一些算法组合去提高预测精度^[17]。 Å smund等^[9]在研究中发现, 土壤反射率不仅仅是土壤理化属性和过程的集中体现, 微环境也也控制实测光谱的变化特征。 对于土壤有机质的含量的高低, 土壤颜色是一个重要的特征, 引入与表征颜色指标密切相关的植被光谱特征(DI₁)和水汽光谱特征(DI₂), 在建立模型时将DI₁和DI₂作为输入以提高模型预测精度。 对获取到的13个敏感波段在SPSS软件中进行主成分分析, 一般来说主成分特征值大于1和方差贡献达到90%以上即可作为主要分量的标准。 为尽可能多挖掘潜在有效光谱信息共提取到6个主成分并在光谱分析的基础上建立光谱特征指数。 光谱特征指数构成形式通常为差值、 比值和归一化指数, 归一化指数更适合PCA变换后的引入。 在Vis-NIR范围对部分特征波段(R₄₉₁, R₆₇₀, R₇₂₅, R₈₀₀)进行差值和比值变换, 得到水分光谱调整指数(DI₁)和植被光谱调整指数(DI₂)^{[1, 18, 19]}。

DI1=λ1[maxRi-Rj(maxRi-Rj)]j∈(600, 601, …, 1075)(1)

DI2=λ2[Rj(maxRi-Rj)-maxRj(770≤j≤870)]j∈(600, 601, …, 1075)(2)

通过建立线性和非线性两种模型进行分析比较得到更适合进行有机质预测的模型, 分别是基于经典统计学和最小二次拟合法的多元逐步回归(MLR)和基于非线性差分函数权值训练的误差反向传播神经网络(BPN)。 为能准确评价模型的预测精度和效果, 我们将样品按照抽样法分为校准部分和验证部分, 分别包含89个样点和22个样点。 通常选择R²(determination coefficients), RMSE(root mean square error)和RPD值来评价模型的预测精度和效果。 R²和RMSE分别用来评价解释因变量程度和精确度的能力, R²越大(0≤ R²≤ 1)和RMSE越小意味着模型拟合性和稳定性越好, 一般来说稳定可靠的模型有高的R²和低的RMSE。 RPD值是衡量预测模型能力的一个重要标准, 其意义是表征预测模型相对于仅使用均值而提高精度的能力, 一般来说, RPD越大, 预测的结果越容易接近实测值, 而低于1, 意味着模型难以对土壤有机质含量做出预测^[20]。

土壤反射率是有机物、 无机物和水等物质组合所产生的光谱的累积性质, 同时也包含了局地微地貌的特征信息, 在一定程度上, 反射率随土壤性质的改变而改变^[22]。 图2展示不同有机质含量下的土壤光谱反射率, 从中可以看出土壤光谱曲线的大致特征均为较平滑的上升趋势, 在可见光部分曲线斜率大于近红外部分, 上升更快; 在350 nm附近存在明显的下凹吸收峰, 在600 nm附近存在一个较陡的反射峰。 而反映有机质的光谱曲线陡坡应该在800 nm处, 因此可能存在着土壤有机质的轻微退化^{[4, 6]}, 与刘文全等^[23]研究中认为此研究区内存在着轻微的盐渍化的结论一致。 此外不同有机质含量土壤光谱曲线趋势大致相同, 且随着有机质含量的增加, 反射率降低, 这些明显光谱特征差异都为建立有机质预测模型提供基础。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 不同有机质含量土壤的原始光谱反射率Fig.2 Soil original spectral reflectance with different soil organic matter contents

将反射光谱与有机质进行相关分析, 发现相关系数绝对值较小且变化幅度不大, 难以区分敏感波段。 将反射光谱进行一阶导数变换^[16]后相关系数正负交替变化, 数值存在较大差异, 可以根据相关系数值的大小判断出特征波段(图3)。 图3中可以发现相关系数r> 0.4时已经属于较高层次, 因此作为指标选取的一个重要参考^[15]。 在SOM与一阶导数光谱保持高相关性前提下, 选择的敏感波段保持Δ =50 nm间距下涵盖可见光和部分近红外光谱范围内的信息, 并且不包含在噪声波段和水吸收峰波段范围内, 选取了13个相关性较高的波段, 对应的中心波长分别为345, 491, 594, 627, 718, 723, 748, 810, 873, 889, 931, 981和1 075 nm(图3)。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 有机质含量与原始光谱、 一阶导数光谱相关关系曲线Fig.3 Correlation coefficient curves between Soil spectral reflectance and soil organic matter contents

表2

Table 2

表2(Table 2)

表2 主成分各分量载荷值 Table 2 Loading value of first six principal components 特征 波段 主分量 1 2 3 4 5 6 R345 -0.759 0.377 0.425 0.019 0.273 0.087 7 R491 0.590 0.640 -0.419 0.140 0.180 -0.007 7 R594 0.882 0.295 -0.136 -0.236 -0.138 0.043 0 R627 0.951 0.272 -0.014 -0.116 -0.090 -0.007 1 R718 0.984 0.020 0.061 0.08 -0.024 -0.023 6 R723 0.967 0.149 0.100 0.098 0.013 -0.002 0 R748 0.987 0.008 0.061 0.037 -0.005 -0.000 2 R810 0.849 0.078 0.464 -0.029 -0.067 -0.067 9 R873 0.911 0.022 0.204 0.068 0.036 -0.194 1 R889 0.906 -0.132 -0.042 -0.16 0.196 -0.079 5 R931 0.774 -0.430 -0.101 0.416 0.044 0.006 0 R981 0.751 -0.503 -0.100 -0.296 0.212 0.042 7 R1 075 0.910 -0.040 0.107 0.053 -0.025 0.368 1

表2 主成分各分量载荷值 Table 2 Loading value of first six principal components

对特征波段进行主成分分析后提取到前6个主成分。 表2可以看出, R₆₂₇, R₇₁₈, R₇₂₃, R₇₄₈, R₈₈₉和R_{1 075}载荷值在主成分1上较大, 所包含的光谱信息与第一主成分相关性最强, 这说明主成分1所代表的信息集中在红波段和近红外波段; R₄₉₁, R₉₃₁和R₉₈₁的载荷值在主成分2上较大, 其他波段在主成分2上的载荷值较小, 表明蓝波段和近红外波段的信息主要反映此主分量上; 主成分3的信息集中在R₄₉₁, R₈₁₀和R₈₇₃, 载荷值较高的波段来源自蓝波段和近红外波段附近; 主成分4, 5和6的有效信息来自R₈₇₃, R_{1 075}和R₃₄₅。 故可以判定敏感波段中与有机质相关的光谱信息主要是可见光波段至近红外波段附近处(600~900 nm), 与林^{[6, 23]}等的研究结果基本一致。

运用主成分提取到光谱特征作为MLR模型的输入自变量, 有机质的含量为因变量, 得到PCA-MLR[图4(a)]。 将植被水汽光谱特征(DI₁与DI₂)加入自变量后, 得到PCA-DI-MLR[图4(b)], PCA-DI-MLR的预测值曲线与实测值曲线更加吻合, 预测效果较好, 而PCA- MLR曲线整体上比较平滑, 变化幅度较小, 对极值点无法进行有效反演; 同样的, 对BPN模型将光谱特征(DI₁与DI₂)作为输入, 有机质的含量作为训练目标, 训练后分别得到最优PCA- BPN[图4(c)]和PCA-DI-BPN[图4(d)], PCA-BPN总趋势呈现较小的变化幅度, 对极小值点的预测较为精确, 但在极大点上预测误差较大, 预测值曲线的走向与实测值基本一致。 PCA-DI-BPN模型中, 预测值与实测值曲线拟合的效果最好, 曲线的趋势走向也基本一致, BPN模型非线性拟合优势体现较为明显, 对极值点的预测误差也最小, 证明了本文建立的PCA-DI-BPN可以有效的对土壤有机质含量进行预测。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 有机质实测值与预测值的比较 (a): PCA-MLR; (b): PCA-DI-MLR; (c): PCA-BPN; (d): PCA-DI-BPNFig.4 Comparison of measured and predicted values of soil organic matter (a): PCA-MLR; (b): PCA-DI-MLR; (c): PCA-BPN; (d): PCA-DI-BPN

从图5不同模型拟合效果图可以看出, PCA-DI-BPN模型的预测值与与实测值的拟合效果最好, 预测值曲线与实测值曲线的走向基本一致, 预测值的误差最小。 图6为有机质预测值与实测值相关关系图, 可以看出基于MLR建立的PCA-DI-MLR模型决定系数R²(0.712)大于PCA-MLR R²(0.643), 而PCA-DI-MLR模型的均方根误差RMSE(0.963)小于PCA-MLR(1.085)的RMSE, 说明PCA-DI-MLR的预测精度要高于PCA-MLR, 相对于PCA-MLR引入水汽植被光谱特征(DI)后R²提高14%、 RMSE降低10%, 提高了模型预测精度。 同样的, 图6可以得出PCA-BPN和PCA-DI-BPN的决定系数R²和均方根误差RMSE分别为0.704, 0.764和0.166, 0.104, BPN模型在训练过程中, 可以保持数据稳健性和抑制特异质, 使预测模型的误差值RMSE降低38%, 达到较高的预测水平。 基于PCA-BPN和PCA-MLR的预测模型RPD值分别为1.419和0.471, 相对比来说, 基于BPN的预测模型可使RPD值整体提升0.9左右, 对有机质的预测估算是可靠的。 而PCA-MLR模型的精度提升效果低于PCA-BPN模型, 预测能力有待提升。 引入水汽植被光谱特征(DI)后MLR和BPN模型的预测结果均有所提升, MLR模型中的RPD值提升了0.7, 提升效果大于BPN模型的0.2, 表明光谱特征用在MLR模型中效果更

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 不同模型下有机质实测值与预测值的比较Fig.5 Comparison of measured and inversion values of soil organic matter content in different model

好。 PCA-DI-BPN模型的R², RMSE和RPD值分别为0.764, 0.161和1.590, 模型预测值与实测值曲线基本吻合, 在相关性分析中有机质含量点基本分布在趋势线的两侧, 所以本模型能实现对有机质的较准确预测。

土壤有机质含量的快速无损估测对耕地质量评价、 精准农业可持续发展、 高分辨率土壤数字化制图等方面具有重要意义^[24]。 通过波段间隔(Δ =50 nm)和相关性大小(r> 0.4)筛选得到敏感波段进行主成分分析, 消除共线性和进一步分析判识光谱信息的来源, 波段间进行算法运算构建植被和水汽光谱特征, 既减少了变量数, 又在BPN模型中得到最精确的预测效果。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 土壤有机质含量实测值与预测值之间的相关性 (a): PCA-MLR; (b): PCA-DI-MLR; (c): PCA-BPN; (d): PCA-DI-BPNFig.6 Correlation between predicted soil organic matter (SOM) and the measured SOM field validation (a): PCA-MLR; (b): PCA-DI-MLR; (c): PCA-BPN; (d): PCA-DI-BPN

通过提取到主成分光谱特征后进行分类判别, 引入光谱特征(DI₁和DI₂)后建立的预测模型优于原模型。 特征波段实质是有机质组分(纤维素、 木质素和果胶等)的光谱吸收特征的一个叠加表现。 特征波长的间隔过小可导致信息冗余和共线性, 间隔过大则会导致光谱信息不足。 主成分分析和典型自相关分析经常应用在环境污染物来源识别上^[25], 同样也有助于判断选择出的光谱特征来源于土壤系统功能集合^[4]。 在6个主要分量上的特征信息与R₃₄₅, R₄₉₁, R₇₁₈, R₇₂₃, R₇₄₈, R₈₇₃, R₈₁₀, R₉₃₁, R₉₈₁和R_{1 075}体现出较高的相关性, 这些信息90%以上来自近红外波段, 这与Brown等^[7]在对SOM反演预测时所得到的敏感波段区域(700~1 100 nm)结论一致, 表明了提取到的光谱特征是有效的。 Rossel等^[26]发现了在近红外波段内的光谱特征是由C— O, C=O, N— H, O— H和C— H bonds组合反应的结果。 所以说SOM在NIR附近的光谱特征, 实际上是这些化学键组分相互作用的体现。 627~748 nm之间的敏感波段94%以上的信息集中在第一主成分, 所以说在此波段内的光谱信息可以作为PC1的典型代表, Ben-Dor^[4]研究指出, 在664 nm附近土壤光谱曲线的吸收特征主要来自叶绿素等颜色残留物, 光谱特征仅仅在前两周比较明显, 随着时间变化不断分解光谱特征不断减弱。 本实验收集光谱的时间也为冬小麦收割后一周, 土壤有机质中的色素残留物含量较高, 可以判断第一主成分代表着有机质中的色素残留物。 R₈₈₉, R₉₈₁和R_{1 075}波段在PC2上的信息贡献量为负值, 光谱的反射峰在此主成分上可能占据主导作用, 本区域内存在着轻微的土壤盐渍化, 盐分的光谱信息一定程度上掩盖SOM的信息影响到有机质信息准确提取^[27], PC2所构建的光谱特征, 在MLR和BPN模型中起到的作用最小; PC3的信息在345~600和810 nm处相关性较高, 可见光波段的信息决定PC3, Ladoni^[28]研究中得出有机质估算中可见光区域内的反射率主要由发色团和腐殖酸的的浓度决定, PC4和PC6上各波段的信息集中在近红外波段, Shi等^[15]的研究表明700~1 100 nm存在着C— H和O— H的敏感波段, 因此主成分4和 6表征SOM中水汽和矿渣相关物质。 在采集样品时表层土壤干燥, 故水汽在本实验中未占据主要信息。 PC5所包含的信息分布特征不明显, 极有可能与灰分物质有关。 详细的推断较为困难, 将其归为土壤受到环境人影响为的信息集合。

MLR和BPN的预测精度受到研究区一些环境属性(大气、 盐分、 粒度粘度和土壤生物化学过程等)和人为因素(样点布设等)的影响。 不同地区土壤的成土条件、 成土母质、 成土过程具有区域独特性, 区域内土壤属性和微地貌属性无法都被考虑进去, 这些变量在异质区域会显著影响到MLR和BPN的性能。 所以在SOM的预测中光谱信息的优化选择和模型匹配可能是提高现场估测精度的一项重要的工作; 此外新引入的植被和水汽光谱特征作为表征土壤颜色的自变量, 是否与其他主成分产生共线性, 进而影响BPN模型的精度, 未进行深入的验证。