三江平原位于中国黑龙江省东北部, 由松花江、 黑龙江和乌苏里江汇流冲积而成的沼泽化低平原^[10]。 地处湿润、 半湿润、 中温带大陆性季风气候, 夏季多雨潮湿, 冬季寒冷干燥, 年平均气温在1~4 ℃, 年平均降水量在450~650 mm。 土壤类型主要包括黑土、 草甸土和白浆土^[12]。 该地区耕地面积达667.79万hm², 占三江平原总面积的55.28%, 是中国北方主要商品粮基地^[11]。 选择了三江平原典型耕地区域— — 友谊农场, 位于三江平原中部(图1)。 该农场始建于1954年, 2013年场区总面积1 684.28 km², 其中耕地面积占比达80.67%, 为1 358.66 km^2[12]。 友谊农场作为中国最大的机械化国营农场, 除西南少数区域外, 地势相对平缓, 土地接连成片, 适合大面积机械化生产。 该地区主要种植农作物包括玉米、 大豆和水稻, 作物一年一熟^[13]。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 友谊农场土壤采样点空间分布Fig.1 Distribution of sampling points in Youyi Farm

为保证样本的空间代表性, 考虑土壤类型和土地利用类型, 并考虑道路可达性限制, 于2021年10月作物收获后, 在友谊农场农田内共采集103个表层(0~20 cm)土壤样品(图1), 其中旱地土壤样品46个, 水田土壤样品57个。 每个土壤样品由5 m范围内采集的3~5个土壤子样混合而成。 在去除作物秸秆、 根系等非土壤物质后, 将土壤样品保存在密封袋中。 土壤样品采集时, 同时记录了采样点的经纬度。 野外采集土壤样品运送至实验室后, 进行室内风干然后进行土壤样品研磨过筛, 用于土壤光谱和肥力属性测量。 SOM的测定采用重铬酸钾容量法, TN采用半微量凯式定氮法测定, TP的测定采用酸消化和比色法, TK采用酸解和原子吸收光谱法测定^[14]。

采用ASDFieldSpec4便携式光谱仪进行土壤光谱测量, 光谱仪波长范围为350~2 500 nm, 光谱分辨率为3 nm(350~1 000 nm)和8 nm(1 000~2 500 nm)。 土壤样品放置在直径为5.6 cm、 深0.8 cm的样品盒中。 光谱测量在暗室中进行, 采用75 W的卤素光作为照明光源, 光源的天顶角为45° , 其与土壤样品表面的距离为20 cm。 采用8° 视场角的光谱探头进行光谱测量, 测量时探头垂直于土壤表面, 距土壤样品表面为18 cm。 每个土壤样品通过90° 旋转测量四个方向的光谱, 取四个方向测量光谱的平均值作为该样品土壤光谱。 每隔30 min, 采用标准Spectrolon参考板对光谱仪进行校正。

由于光谱仪在波段首尾部分的信噪比相对较低, 本研究去除350~399和2 401~2 500 nm噪声较大的波段, 保留400~2 400 nm范围内的波段进行模型构建。 为进一步降低噪声的影响, 对测量光谱采用Savitzky-Golay(SG)平滑算法进行平滑处理, 并采用标准正态变换(standard normal variate, SNV)和多元散射校正(multiplicative scatter correction, MSC)进行基线校正和散射影响的消除; 通过一阶微分变换(first derivative, FD)处理突出土壤光谱吸收特征。 本研究将探讨这四种光谱处理方法在土壤肥力属性预测中的表现。

采用土壤光谱范围为400~2 400 nm, 共计2 001个波段, 土壤样品共计103个, 属于小样本高维数据, 且不同波段间存在一定共线性。 直接采用所有波段构建模型容易产生过拟合问题, 通过波段筛选可以过滤无效波段, 提升模型构建效率和精度。 采用CARS算法结合偏最小二乘回归(partial least square regression, PLSR)模型(CARS-pls)分别针对各预测土壤肥力属性进行变量筛选。 算法执行步骤: (1)采用蒙特卡罗(Monte Carlo, MC)方法随机选择一定数量的土壤样本作为建模集, 使用该建模集构建PLSR模型; (2)采用PLSR模型中波段系数的绝对值, 采用指数递减函数(exponentially decreasing function, EDF)确定的变量个数选择有效波段; (3)采用自适应重加权采样(adaptive reweighted sampling, ARS)方法进一步去除不重要波段; (4)计算所选择波段下基于交叉验证(CV)的均方根误差值(RMSE)^[15]。 CARS-pls的初始主成分数为15, MC随机选择土壤样本的次数为50, 每次均可生成一个RMSE, RMSE为最小值时的波段集合为筛选出的波段, 用于模型构建。

以CARS算法筛选出的波段, 采用PLSR算法分别构建SOM、 TN、 TP和TK的预测模型。 由于土壤样本数量相对较少, PLSR模型训练过程中, 通过留一交叉验证进行主成分个数寻优并进行模型评价。 主成分的最大个数设为20。 采用决定系数(R²)、 均方根误差(RMSE)、 相对分析误差(RPD)以及性能与四分位间隔距离的比率(RPIQ)评价模型预测效果。 与RPD相比, RPIQ不容易受到变量分布的影响。 当R²、 RPD和RPIQ越大, RMSE越小时, 模型预测效果越好。 见式(1)— 式(4)

$R^{2}=1-\frac{\sum_{i=1}^{n}\left(y_{i}-\hat{y}_{i}\right)^{2}}{\sum_{i=1}^{n}\left(y_{i}-\bar{y}\right)^{2}}$(1)

$\mathrm{RMSE}=\sqrt{\frac{\sum_{i=1}^{n}\left(\hat{y}_{i}-y_{i}\right)^{2}}{n}}$(2)

RPD=SD/RMSE(3)

RPIQ=IQRMSE(4)

上述中, n为土壤样本个数, y_i为实测值, ŷi为预测值, y̅为均值, SD为标准差, IQ是第三和第一四分位数之间的差值。

针对PLSR的最佳预测模型, 采用变量投影重要性(variable importance in projection, VIP)^[16]衡量土壤肥力属性预测中各波段相对重要性。 一般认为VIP大于1的波段为重要波段, VIP得分越高说明该波段在土壤肥力属性预测中所发挥的作用越强^[17]。

四种土壤肥力属性: SOM、 TN、 TP和TK的描述性统计结果见图2。 SOM含量范围为2.27%~6.84%, 平均值为4.18%, 采集土壤样本SOM主要分布在2.50%~5.00%之间, 属于较高水平; TN含量范围为0.12%~0.34%, 平均值为0.20%, 以TN含量在0.15~0.25%的土壤样本较多; TP含量范围为388.22~884.99 mg· kg^-1, 平均值为647.52 mg· kg^-1, 较多样本TP含量在550~700 mg· kg^-1; TK含量范围为1.95%~3.08%, 平均值为2.39%, 以含量在2.20%~2.60%的采集土壤样本居多。 尽管经过50多年的农业开垦利用, 根据第二次全国土壤普查的土壤养分分级标准^[18], SOM、 TN、 TK和TP分别为一级、 二级、 三级和二级标准, 土壤肥力整体仍处于较高水平。 这四种土壤肥力属性整体上符合正态分布(图2), 因此构建光谱预测模型时无需对土壤肥力属性值进行转换处理, 直接采用测量值进行预测模型构建。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 土壤肥力属性直方图Fig.2 Histogram of soil fertility attributes

图3为四种不同光谱预处理后的土壤光谱。 结果表明, 在整个波段范围内SG平滑光谱的反射率在0.05~0.50之间, 小于750 nm时, 由于SOM含量较高, 反射率小于0.2。 在1 400、 1 900和2 200 nm附近存在明显的吸收谷。 其中1 400和1 900 nm附近是土壤矿物中水的吸收特征, 而2 200 nm是土壤中矿物的吸收特征。 经FD变换后, 土壤光谱的吸收特征更为明显。 与SG光谱相比, FD光谱中1 400、 1 900和2 200 nm附近的吸收特征明显增强。 在可见光波段410和430 nm附近出现了2个小的吸收峰, 在2 000 nm后也可见较多小吸收峰。 因此通过FD处理可以更好突出土壤的光谱特征。 SNV和MSC处理后的土壤光谱曲线较为相似, 而且光谱曲线上的光谱特征也与SG土壤光谱相近。 与SG光谱相比, 不同土壤样品间SNV和MSC光谱的变异范围更小。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 不同光谱预处理土壤光谱(平均值和± 标准差)Fig.3 Soil spectra (mean and± standard derivation) with different spectral processing methods

将经不同预处理的土壤光谱与土壤肥力属性进行相关分析, 结果如图4所示。 SG光谱在整个波段范围内与土壤肥力属性间的相关系数在-0.50~0.50之间, SOM和TN与光谱反射率间主要呈负相关关系, SOM的相关系数绝对值的最大值出现在775 nm, 为-0.47, TN的相关系数绝对值的最大值出现在798 nm, 相关系数为-0.41。 与SOM和TN相比, TP与光谱反射率间的相关性相对较弱, 在1 364 nm前呈负相关关系, 之后呈较弱的正相关关系, 相关系数绝对值的最大值出现在678 nm, 但仅为-0.25。 与之类似, TK与光谱反射率间的相关性也相对较弱, 多呈正相关关系, 仅在558 nm前呈负相关关系, 相关系数绝对值的最大值也在此区间, 仅为-0.36(405 nm处)。 FD处理后的土壤光谱与土壤肥力属性间的相关性均有了一定程度的提升, SOM、 TN、 TP和TK的相关系数绝对值的最大值分别为0.81(1 378 nm处)、 0.79(1 376 nm处)、 0.46(2 239 nm处)、 0.64(2 228 nm处)。 SOM和TN的相关系数曲线在600和1 400 nm附近存在负值峰, 在1 300、 1 700、 2 200和2 300 nm附近存在明显正值峰。 与SG土壤光谱相比, SNV和MSC光谱与各土壤肥力属性间的相关性均得到了明显提升, SOM和TN相关系数绝对值的最大值接近0.75, 而TP和TK的最大值接近或超过了0.50。 上述分析说明, 与SG土壤光谱反射率相比, 经FD、 SNV和MSC处理可以增强土壤光谱与土壤肥力属性间的相关性。 需要指出的是, 尽管各土壤肥力属性与不同光谱预处理间的相关性大小有差异, 而不同土壤肥力属性相关系数曲线间的关系在不同预处理光谱之间较为一致。 在四种预处理光谱中, SOM和TN的相关系数曲线十分接近, 而TP和TK的相关系数曲线呈现一定的对称性。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 不同预处理土壤光谱与土壤肥力属性间的相关性Fig.4 Correlation between soil spectra with different spectral processing methods and soil fertility attributes

使用全波段与CARS算法筛选波段、 采用不同预处理光谱的土壤肥力属性交叉验证结果如表1所示。 使用未经变量筛选全部光谱波段时(400~2 400 nm), SOM和TN的预测精度相对较高, SOM的预测R²在0.85~0.89之间、 RPD在2.60~2.98之间, TN的预测R²在0.86~0.89之间、 RPD在2.68~3.00之间。 TP的预测精度较低, R²仅在0.08~0.34之间、 RPD小于1.40。 与TP相比, TK的预测精度相对较高, 但低于SOM和TN, 预测R²在0.63~0.72之间、 RPD在1.65~1.89之间。 使用CARS算法筛选波段的四种土壤肥力属性预测精度均得到了一定程度提升。 不同光谱预处理方法SOM和TN的预测精度均较高, SOM的最佳预测精度来自FD处理光谱, R²、 RMSE、 RPD和RPIQ分别为: 0.97、 0.19%、 5.54和8.11, TN的最佳预测精度同样来自FD光谱, R²、 RMSE、 RPD和RPIQ分别为: 0.96、 0.01%、 4.75和7.56。 TP的预测精度得到了明显提升, 最佳预测精度时R²、 RMSE、 RPD和RPIQ分别为: 0.82、 44.48 mg· kg^-1、 2.39和3.49。 但不同光谱预处理方法的TP预测精度差异明显, SG光谱的预测精度最低, R²、 RMSE、 RPD和RPIQ分别为: 0.49、 75.25 mg· kg^-1、 1.41和2.06, 该预测精度已经优于全波段时SNV光谱得到的最佳预测结果。 TK的预测精度也得到了一定提升, 最佳预测精度时R²、 RMSE、 RPD和RPIQ分别为: 0.92、 0.06%、 3.66和5.52。 不同光谱预处理方法间的TK预测精度也存在一定差异, 而CARS筛选波段时得到的最低预测精度(SG光谱时)仍高于全波段时得到的最佳预测结果(SNV光谱时)。 以上结果表明, 经CARS算法进行波段筛选后各土壤肥力属性的预测精度均有一定提升, 其中TP的提升幅度最大。

表1

Table 1

表1(Table 1)

表1 不同光谱预处理土壤肥力属性预测结果 Table 1 Estimation results of soil fertility attributes based on different spectral processing methods 土壤属性 预处理方法 全波段 CARS波段筛选 R2 RMSE RPD RPIQ R2 RMSE RPD RPIQ SOM/% SG 0.86 0.40 2.67 3.91 0.89 0.36 2.97 4.35 FD 0.85 0.41 2.60 3.81 0.97 0.19 5.54 8.11 SNV 0.89 0.36 2.98 4.36 0.93 0.28 3.89 5.69 MSC 0.87 0.38 2.81 4.11 0.91 0.32 3.35 4.90 TN/% SG 0.86 0.02 2.68 4.26 0.91 0.02 3.41 5.43 FD 0.86 0.02 2.69 4.28 0.96 0.01 4.75 7.56 SNV 0.88 0.02 2.94 4.68 0.92 0.01 3.51 5.59 MSC 0.89 0.02 3.00 4.78 0.93 0.01 3.93 6.26 TP/(mg· kg-1) SG 0.08 101.37 1.05 1.53 0.49 75.25 1.41 2.06 FD 0.27 90.56 1.17 1.71 0.82 44.48 2.39 3.49 SNV 0.34 85.89 1.24 1.81 0.58 68.17 1.56 2.28 MSC 0.34 86.19 1.23 1.80 0.58 68.52 1.55 2.26 TK/% SG 0.63 0.14 1.65 2.49 0.74 0.11 1.98 2.98 FD 0.65 0.13 1.71 2.58 0.92 0.06 3.66 5.52 SNV 0.72 0.12 1.89 2.85 0.79 0.10 2.20 3.31 MSC 0.67 0.13 1.75 2.63 0.85 0.09 2.61 3.93

表1 不同光谱预处理土壤肥力属性预测结果 Table 1 Estimation results of soil fertility attributes based on different spectral processing methods

各土壤肥力属性基于全波段和CARS筛选波段最优预测模型的预测散点图(图5)表明, 对于SOM、 TN和TK, 使用全波段时的预测值和实测值已经较为接近, 而经CARS变量筛选后的预测值和实测值更为接近, 基本都分布于1∶ 1线附近, 拟合线与1∶ 1线偏差更小, 因此预测精度更高。 使用全波段时, TP的预测值与实测值间的差异较大, 经CARS变量筛选后预测值与实测值间的关系明显改善, 大部分预测值和实测值都紧密分布于1∶ 1线附近, 而且拟合线与1∶ 1线间的偏差明显降低, 低值区预测值偏高和高值区预测值偏低的情况也得到了显著改善。 CARS算法能够显著改善土壤肥力属性预测精度, 特别是TP的预测精度得到了明显提高。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 基于全波段和CARS筛选波段最优预测模型的土壤肥力属性预测散点图Fig.5 Scatter plot of soil fertility attributes based on full bands and selected bands by CARS algorithm, respectively

针对各土壤肥力属性最佳预测模型, 采用CARS算法筛选波段计算的VIP值如图6所示。 结果表明, 经CARS筛选后SOM、 TN、 TP、 TK的特征波段数分别为: 137、 278、 90、 90, 占全部光谱数的比例分别为: 6.85%、 13.89%、 4.50%、 4.50%, 有效减少了输入波段数量。 SOM与TN使用的波段在整个光谱范围内均有分布, 而更多集中在近红外光谱区域内(大于1 000 nm), 其中610、 1 390、 1 890和2 211 nm附近波段在SOM的预测中发挥了更大作用; 775 nm附近波段对于TN的预测也很重要。 与SOM和TN相比, TP和TK筛选出的波段相对较少, 其中750~1 000 nm范围内的波段几乎没有被选择。 与SOM和TN在1 000~2 500 nm范围的选择波段类似, TP的主要波段也集中在1 390、 1 890和2 240 nm附近, 560、 2 366和2 390 nm附近波段也具有重要作用; 与SOM、 TN和TP不同, 910 nm附近波段在TK的预测中具有较大作用, 大于2 200 nm的波段中, TK的预测中具有较大作用的波段集中在2 350 nm附近。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 基于CARS算法筛选的土壤肥力属性最优预测模型特征波段的VIP值 VIP: 变量投影重要性Fig.6 VIP values in characteristic bands selected by CARS algorithm for the optimal prediction models of the four soil fertility attributes VIP: Variable importance in projection

采用CARS算法筛选波段, 并结合不同光谱预处理方法进行土壤肥力属性光谱预测。 结果表明, 与使用全波段相比, 经CARS变量筛选后SOM、 TN、 TP和TK的预测精度均有所提升, 其中TP的提升幅度最大, TK的提升幅度次之, SOM和TN的提升幅度最小。 SOM和TN的预测精度最高(R²≥ 0.95, RPD≥ 4.70), TK的预测精度次之, TP的预测精度相对较低, 而R²也大于0.80, RPD高于2.00。 与已有研究结果相比(已发表结果SOM和TN的预测统计中值R²=0.86^[4]), 本研究经CARS变量筛选后SOM和TN的预测精度更高。 TP和TK的预测精度在不同研究中存在较大差异。 李学兰等^[9]在针对贵州省土壤肥力属性的光谱预测中发现, TP的预测精度较低, R²仅为0.09, TK的精度也不高, R²为0.58。 乔璐等^[19]在黑龙江大庆市区和四县(肇源县、 肇州县、 杜蒙、 林甸县)的检测结果, TP和TK的精度较高, R²分别达0.687和0.832。 与SOM和TN相比, TP和TK缺少明确的光谱吸收特征, TP和TK的预测精度很可能受到与其他具有明确光谱特征的土壤属性之间的关系强弱影响, 体现出一定的区域性或者受到土壤类型等因素的影响^[5]。 与已有研究结果相比, 本研究TP和TK经变量筛选后的预测精度相对较高, 具有实际预测的潜力。

光谱吸收特征是土壤属性稳定精确预测的基础。 本研究针对SOM预测筛选的波段中, 1 890 nm附近波段可以归为SOM中C— H和/或C═O键的倍频与谐频振动, 2 211 nm可能为土壤粘土矿物中的Al— OH引起^[20]。 SOM中的C— H键在825 nm附近的波段并未被选择, 可能是由于CARS算法的变量筛选机制中仅考虑了被预测土壤属性的精度; 也可能是由于土壤本身组成的复杂, 导致了不同土壤物质吸收特征重叠, 进而难以识别出相应的吸收波段。 土壤中的氮主要是有机态, 有机氮中的N— H键在751、 1 000、 1 500和2 060 nm附近存在吸收特征^[5], 本研究中筛选出的775 nm附近波段可能与N— H的振动有关。 Matin等^[21]认为TN的光谱预测会受到其与SOM之间的关系影响, 本研究中SOM和TN之间的相关系数为0.96, 因此TN的高精度预测结果可能同时利用了其自身的光谱吸收特征以及其与SOM之间的关系。 与SOM和TN预测用到的波段不同, 560、 2 366和2 390 nm附近波段在TP预测中起到了较大作用, 同时考虑到TP与SOM和TN之间的相关系数不高(分别为0.37和0.40), 因此TP的较高精度预测可能由于CARS变量筛选方法筛选出了TP的有效预测波段。 与SOM和TN相比, 910和2 350 nm附近波段在TK的预测中也具有较大作用, 而与TP不同, TK与SOM和TN之间的相关系数相对较高(分别为-0.61和-0.63), 因此TK的高精度预测可能利用了其与SOM和TN之间的关系。

本研究发现光谱预处理方法能够影响土壤肥力属性的预测精度。 SOM和TN预测受到的影响相对较小, TP受到的影响最大, TK受到的影响介于两者之间。 对于TP, SG光谱的预测精度最低(R²仅为0.49), 使用FD处理光谱后预测精度有了明显提升(R²可达0.82)。 对于SOM、 TN和TK, 同样是SG光谱的预测精度最低, FD光谱的预测精度最高。 Hou等^[22]在针对黄土高原花椒林(乡镇尺度)的研究中也发现, 与MSC和对数倒数等预处理方法相比, FD处理在SOM和TN的光谱预测中具有更好的表现。 史舟等^[23]发现FD预处理方法在基于全国尺度土壤光谱的SOM预测中也具有良好表现。 本研究与其他三种光谱预处理方法相比, FD处理更能突出土壤光谱特征, 因此与CARS变量筛选方法结合可以提升土壤肥力属性预测精度。 勾宇轩等^[24]在进行SOM的光谱预测时发现, 与原始光谱、 倒数对数和连续统去除方法相比, 采用连续小波变换处理光谱并与CARS波段筛选结合可以提升不同类型土壤SOM预测精度。 在将来研究中可以探索小波变换与CARS算法结合在其他土壤肥力属性中的预测效果。 有研究表明模型构建方法在土壤属性预测中也具有重要影响^{[9, 20, 25]}, 本研究中仅使用了PLSR一种方法, 在将来的研究中将结合其他机器学习算法进一步提升土壤肥力属性的预测精度。