土壤有机碳(soil organic carbon, SOC)影响土壤肥力, 可以调控土壤生态系统的功能发挥[1, 2]、 维持结构稳定[3, 4], 对于维持农业可持续发展, 缓解温室气体增加和全球气候变化的影响具有重要意义[5]。 传统的土壤有机碳含量测定技术成本高、 周期长, 高光谱技术连续波段数量足、 光谱分辨率高, 可以快速、 高效的进行土壤有机碳的测定[6, 7]。 土壤高光谱数据在采集过程中, 易受外界测试环境及土壤样品质量、 仪器本身因素的影响, 导致光谱变量之间存在大量的无关冗余信息[8]。 连续小波变换(continuous wavelet transform, CWT)通过多尺度分解, 能够有效提取光谱信号中的近似特征和细节信息, 捕捉土壤光谱特征, 提取有效信息, 在高光谱数据的噪声去除和数据压缩方面具有独特优势[9, 10]。 Lin等[11]通过CWT放大被噪声掩盖的隐藏光谱信息, 估算了矿区的土壤有机质。 Li等[12]利用CWT进行光谱预处理, 提升了土壤光谱对有机质含量的相关性, 经CWT变换的光谱估算能力明显优于传统光谱变换。 经CWT变换的光谱数据被分解若干尺度, 建立估算模型过程繁琐、 计算量大, 相关学者利用特征波段筛选方法进行光谱数据量的简化。 赵海龙等[13]应用相关系数法(correlation coefficient, CC)结合CARS选择CWT不同分解尺度的特征波段, 可以有效地减少建模的自变量数, 模型稳定性及估算能力良好。 勾宇轩等[14]利用CWT结合稳定自适应重加权采样(stability competitive adaptive reweighted sampling, sCARS)算法进行土壤有机质含量的反演, sCARS算法能够去除冗余、 重叠的光谱信息, 提高建模效率。 蔡亮红等[15]研究表明利用竞争自适应重加权采样(competitive adaptive reweighted sampling, CARS)算法, 在建立模型时可以减少光谱细节损失, 同时还能对光谱无信息变量进行有效去除。 Wang等[16]研究发现CWT经连续投影算法(successive projections algorithm, SPA)处理后, 能够选取有效波长, 减少冗余变量。 Banskota等[17]指出遗传算法(genetic algorithm, GA)可以消除小波系数产生的噪音, 获取敏感信息。 然而, 特征波段筛选方法, 虽然在保证CWT构建模型精度的前提下一定程度降低了光谱数据量, 但是多分解尺度的小波系数仍然难以充分利用, 单一分解尺度小波系数不能有效代表不同尺度小波系数。 章涛等[18]构建小波能量特征向量(energy features vector, EFV), 结果表明EFV不仅能够代表不同尺度小波系数整体信息, 并且小波能量特征向量方法可以有效减少噪声影响, 基于EFV的模型输入的光谱数据维度可控, 估算能力提高。 目前, 应用高光谱小波能量特征向量EFV估算土壤有机碳含量的研究鲜有报道。

以博斯腾湖湖滨绿洲为研究区, 应用CWT对光谱数据进行预处理, 将土壤光谱分解为若干个尺度的小波系数, 利用CC、 sCARS、 CARS、 SPA、 GA常见5种特征波段筛选方法优选不同尺度的小波系数, 结合全波段小波系数分别计算各尺度小波系数的均方根, 将各尺度数据结合构建EFV。 选取具有较强的非线性处理能力及良好的拟合优度的BP神经网络模型(BP neural network model, BPNN)[19]。 以EFV作为建模输入量, 构建EFV-BPNN模型, 探究土壤高光谱噪声去除和信息压缩方法, 实现土壤有机碳含量快速估算。

博斯腾湖湖滨绿洲位于新疆焉耆盆地东南面, 典型的干旱区湖滨绿洲, 地理位置介于86° 15'— 86° 55'E, 41° 45'— 42° 10'N, 面积为1 360 km2; 年均降水量47.7~68.1 mm, 年均蒸发量1 880.0~2 785.8 mm, 年平均气温8.2~11.5 ℃, 蒸降比大于40∶ 1, 无霜期176~211 d; 地下水埋深1.0~2.5 m, 矿化度为0.1~10.0 g· L-1; 土壤有机碳平均含量为13.31 g· kg-1, 研究区主要土壤类型为草甸土、 潮土、 沼泽土等[20]。

1.2.1 土壤有机碳含量测定

依据研究区土壤类型、 土地利用类型以及植被类型状况等因素, 于2020年10月进行土壤样品的采集。 采样点需遍及研究区主要土地利用类型, 共计选取69个样点, 采用分层随机取样的方法, 每个样本点以0~10和10~20 cm采样土层进行, 共计138份0~20 cm表层土壤样品[21]。 处理土壤中杂物后以四分法混合均匀, 选取约200 g土壤样品装袋; 将土样风干、 研磨, 过2 mm筛后重新封装, 以备实验分析; 采用重铬酸钾-外加热法进行土壤有机碳含量的测定[22], 剔除3个异常值, 统计特征见表1。

1.2.2 土壤光谱测定

土壤样品的高光谱数据利用ASD FieldSpec3地物光谱仪室外采集(表2)。 考虑环境因素的影响, 高光谱数据首尾波段的信噪比低, 保留400~2 450 nm, 并去除1 301~1 449与1 801~1 949 nm水汽影响波段, 减少高频噪音对光谱数据的影响[23]。 应用Savitzky-Golay滤波方法进行光谱曲线的平滑处理, 去除干扰波段的光谱曲线。

1.3.1 连续小波变换

连续小波变换(CWT)源于傅里叶算法, 通过小波基函数将一维的土壤高光谱数据分解包含波段和不同分解尺度的二维小波系数, 具有一定程度的高分辨和适应性[12], 见式(1)

Wf(a, b)≤fΨa, b≥∫-∞+∞f(t)Ψa, b(t)dt(1)

其中

Ψa, b=1aΨt-ba(2)

式中: t为高光谱波段数; a为尺度因子; b为平移因子; f(t)为光谱反射率; Ψ a, b为小波基函数。

1.3.2 高光谱小波能量

计算第i个分解尺度的小波系数的均方根, 作为该尺度的E Fi[18], 见式(3)

EFi=1j∑j=1jWj(ai, bj)2(3)

构建EFV

EFV=[EF1, EF2, …, EFi](4)

1.3.3 模型构建

将光谱数据进行连续小波变换预处理, 选择Gaussian4作为小波变换基函数[14], 提高光谱数据的相关性, 放大有效光谱信息。 以5种常规特征波段筛选方法, CC、 sCARS、 CARS、 SPA、 GA算法筛选10个分解尺度的特征波段[18, 20]。 结合全波段数据, 分别计算不同尺度下的均方根, 用以构建高光谱小波能量特征向量, 简化建模数据量, 将复杂光谱数据化简为10维数据; 利用小波能量特征向量作为建模输入量, 构建BPNN模型, 利用R2、 RMSE、 RPD评估小波能量特征向量在土壤有机碳含量估算建模的适用性, 并探究基函数对模型精度的影响。 流程图见图2。

由图3相关分析的结果可知, 土壤有机碳含量与土壤光谱反射率R存在明显的负相关关系, 相关系数范围绝对值|0.005|< r< |0.406|; 相关曲线在400~1 800 nm区间较为平滑, 在1 950~2 450 nm区间呈现波动变化; 相关系数在400~616 nm呈现增加趋势; 相关系数在616~1 300和1 451~1 800 nm呈现降低趋势; 578~616 nm波段内相关性整体最强, 相关系数最大值r=|0.406|; 相关系数平均值r=|0.208|。

Bior1.3、 db4、 Gaussian4与mexh基函数[24]中以Gaussian4为基函数的小波变换结果较为理想, 故选择小波变换基函数为Gaussian4函数, 进行CWT处理, 变换尺度为21, 22, 23, …, 210, 将生成的连续小波变换后的光谱数据与土壤有机碳含量进行相关性分析, 土壤有机碳含量与土壤原始光谱反射率R存在负相关关系。 由图4可知, 有效光谱信息主要集中于中高尺度分解。 在分解尺度为1~4时, 光谱细节特征缺乏, 连续有效光谱信息较少, 随着分解尺度增大, 光谱对土壤有机碳的敏感性也逐步增大。

通过CWT处理后10个分解尺度都存在一定的正负相关性, 土壤光谱与土壤有机碳含量间相关性明显提高, 其中, 4~10分解尺度的相关系数平均值均高于r=|0.208|, 相关系数平均值分别为|0.254|、 |0.299|、 |0.366|、 |0.302|、 |0.345|、 |0.446|、 |0.422|; 1~3分解尺度相关性较差, 相关系数平均值分别为|0.165|、 |0.183|、 |0.206|, 第1尺度相关系数平均值最低r=|0.165|, 第10尺度相关系数平均值最高r=|0.446|。 1~10分解尺度的相关系数最大值分别是|0.623|、 |0.623|、 |0.690|、 |0.730|、 |0.681|、 |0.696|、 |0.678|、 |0.710|、 |0.710|、 |0.677|, 相较于原始光谱与土壤有机碳含量的相关系数最大值r=|0.406|有所提升。

通过CWT处理后, 10个分解尺度的小波系数难以整合利用, 构建的小波能量特征向量能够代表不同尺度小波系数整体信息, 数据维度可控, 有效降低建模输入量[12]。 以CC、 sCARS、 CARS、 SPA、 GA, 分别筛选10个分解尺度的特征波段, 计算不同尺度下的均方根, 构建高光谱小波能量特征向量。

利用CC法, 由图5(a)可知, 分别筛选10个分解尺度通过p< 0.01显著性检验, 相关系数r≥ 0.50的特征波段。 由图5(b)可知, 第1、 3、 4、 7尺度主要分布在400~760 nm可见光波段内; 第2、 8~10尺度主要分布在400~1 500 nm可见光波段及近红外短波; 第5~6尺度主要分布在400~760和1 500~2 500 nm可见光及近红外长波。 不同尺度筛选特征数量占全波段数量的1.14%、 2.00%、 2.22%、 6.73%、 12.84%、 25.84%、 12.26%、 26.58%、 51.34%、 46.04%。 逐尺度计算所筛选特征波段的均方根, 构建CC-EFV。

应用sCARS算法, 分别筛选10个分解尺度特征波段, 以第1尺度为例筛选过程为例, 由图6(a)可知, 优选出的波段变量个数随着运行次数增加逐渐减少, 经过前20次运算, 筛选的波段变量个数明显降低, 在第23次运算后趋于平稳; RMSECV值通过前20次运算, 逐渐降低, 在第23次采样达到最低值, 而后RMSECV值呈升高趋势; 第23次采样RMSECV值最小即所选择的特征波段最优。 图6(b), 第1~3、 5、 7、 9~10尺度选择的特征波长主要位于1 500~2 500 nm近红外长波; 第4尺度主要集中于760~2 500 nm近红外波段; 第6尺度主要分布在400~760 nm可见光波段及1 500~2 500 nm近红外长波; 第8尺度主要分布在400~760 nm可见光波段。 不同尺度筛选特征数量占全波段数量的4.79%、 2.74%、 4.16%、 2.80%、 3.59%、 3.59%、 0.63%、 4.79%、 1.03%、 1.37%。 逐尺度计算所筛选特征波段的均方根, 构建sCARS-EFV。

采用CARS算法, 分别筛选10个分解尺度特征波段, 与sCARS筛选原理基本相似, 以第1尺度筛选过程为例; 由图7(a)可知, 变量个数在第27次运算后趋于平稳; RMSECV在第23次采样达到最低值, RMSECV值最小即所选择的特征波段最优。 图7(b), 第1~5、 7、 9~10尺度选择的特征波长主要分布在1 500~2 500 nm近红外长波; 第6尺度主要分布于400~760和1 500~2 500 nm可见光及近红外长波; 第8尺度主要集中于400~760 nm可见光波段; 不同尺度筛选特征数量占全波段数量的2.74%、 1.83%、 4.73%、 4.79%、 1.20%、 1.83%、 0.80%、 3.14%、 0.68%、 0.91%。 逐尺度计算所筛选特征波段的均方根, 构建CARS-EFV。

利用SPA算法, 分别筛选10个分解尺度特征波段, 以第1尺度筛选过程为例, 由图8(a)可知, 变量个数为8之后RMSE值趋于平稳, RMSE=7.82变量个数为8时, 所选择的特征波段最优。 图8(b), 第1~3、 5~6尺度选择的特征波长主要在1 500~2 500 nm近红外长波; 第4尺度主要分布在400~760和1 500~2 500 nm可见光与近红外长波; 第7尺度主要分布于760~1 500 nm在近红外短波; 第8~10尺度主要在760~2 500 nm近红外波段; 不同尺度筛选特征数量占全波段数量的0.46%、 0.40%、 0.17%、 0.11%、 0.23%、 0.29%、 0.29%、 0.29%、 0.74%、 0.29%。 逐尺度计算所筛选特征波段的均方根, 构建SPA-EFV。

采用GA算法, 分别筛选10个分解尺度特征波段, 由图9(a)可知, GA算法据预定的目标函数对样本进行评估, 给出一个适应度值, 基于此适应度值, 选择最佳样本。 图9(b), 第1~5、 9~10尺度选择的特征波长在400~2 500 nm分布较为均匀; 第6~8尺度选择的特征波长主要分布于760~2 500 nm近红外波段。 不同尺度筛选特征数量占全波段数量的2.45%、 9.87%、 11.12%、 3.59%、 3.99%、 3.88%、 1.94%、 2.05%、 2.91%、 3.94%。 逐尺度计算所筛选特征波段的均方根, 构建GA-EFV。

依据CC、 sCARS、 CARS、 SPA、 GA算法及全波段构建的EFV, 作为模型输入量, 构建土壤有机碳BPNN估算模型。 由表3可知, 利用EFV可以较好估算土壤有机碳含量, 6种模型的训练集与验证集R2平均值分别为0.73、 0.74, RMSE平均值分别为7.64、 7.28, RPD平均值为1.95。 对比模型验证集精度可知, CC-EFV-BPNN> sCARS-EFV-BPNN> Full-spectrum-EFV-BPNN> CARS-EFV-BPNN> GA-EFV-BPNN> SPA-EFV-BPNN。 以CC-EFV-BPNN构建的模型效果最好, R2分别为0.78、 0.77, RMSE分别为6.38、 6.92, RPD为2.07; 其次sCARS-EFV-BPNN构建的模型效果较好, R2分别为0.71、 0.76, RMSE分别为7.66、 6.72, RPD为2.03。 SPA-EFV-BPNN构建模型精度略差R2分别为0.70、 0.70, RMSE分别为8.44、 6.43, RPD为1.82。

由图10可知, 对比训练集样本点Full-spectrum-EFV-BPNN、 CC-EFV-BPNN、 sCARS-EFV-BPNN、 GA-EFV-BPNN基本均匀的分布在1∶ 1线两侧, 在土壤有机碳含量高于40.0 g· kg-1时, 存在少量样本偏移点; CARS-EFV-BPNN、 SPA-EFV-BPNN, 在土壤有机碳含量介于0~20.0 g· kg-1时, 存在样本偏移现象, 估算值与实测值间相差约为10.0%, 样本点在1∶ 1线两侧呈分散分布。 验证集样本点Full-spectrum-EFV-BPNN、 CC-EFV-BPNN、 sCARS-EFV-BPNN基本均匀的分布在1∶ 1线两侧, Full-spectrum-EFV-BPNN在40.0~60.0 g· kg-1间存在两个高估样本, 高估约50.26%; CARS-EFV-BPNN、 GA-EFV-BPNN、 SPA-EFV-BPNN, 与1∶ 1线间距离明显高于其余模型。

应用CWT提取EFV, 在有效降低数据维度与模型复杂度的前提下, 将多分解尺度小波系数整合, 充分利用小波系数整体信息。 对于快速估算土壤有机碳含量具有一定的参考价值。 研究结果表明, CWT可以有效提高光谱反射率与土壤有机碳含量间的相关性, 相关系数平均值提升43.66%, 这与勾宇轩等[14]研究结论基本一致。 高光谱数据在去除水汽及噪声过大的全波段数量为1753, 而构建的EFV仅为10维的光谱特征向量, 可以有效降低数据维度, 提升建模速率, 构建的模型训练集与验证集R2平均值分别为0.73、 0.74, RMSE平均值分别为7.64、 7.28, RPD平均值为1.95, 表明EFV进行土壤有机碳含量的估算结果较好, 这与章涛等[18]研究结果一致。 本工作利用CC、 sCARS、 CARS、 GA、 SPA筛选不同尺度的特征波段进行EFV构建, 相比只利用代表部分敏感信息的特征波段, 特征波段方法结合EFV既可以减少不同尺度的光谱数据量、 突出重要光谱信息, 又能够综合表征不同尺度的光谱信息。 对比全波段数据CC、 sCARS算法可以有效提高模型估算能力, CARS、 GA、 SPA算法虽然未能显著提高模型估算能力, 但是特征波段筛选方法一定程度上降低模型的RMSE值, 提升模型精度, 对比全波段建模结果可知验证集RMSE值平均降低22.80%。 本工作估算模型RMSE值总体偏高, 李冠稳等[24]研究土壤有机质所构建的模型RMSE值也较高, 可能是由于光谱测量实验存在一定误差, 未能实现野外原位光谱测量。 利用BPNN构建土壤有机碳估算模型, 模型能力表现为CC-EFV-BPNN> sCARS-EFV-BPNN> Full-spectrum-EFV-BPNN> CARS-EFV-BPNN> GA-EFV-BPNN> SPA-EFV-BPNN, 说明不同特征波段算法所构建的EFV不同, 最终影响BPNN模型精度。 因此, 基于高光谱小波能量特征向量波段筛选算法的机理仍有待于进一步研究。

连续小波变换可以有效提高光谱反射率与土壤有机碳含量间的相关性, 1~10分解尺度相关系数分别为|0.165|、 |0.183|、 |0.206|、 |0.254|、 |0.299|、 |0.366|、 |0.302|、 |0.345|、 |0.446|、 |0.422|, 相关系数平均值提升43.65%; 相关系数最大值分别为, |0.623|、 |0.623|、 |0.690|、 |0.730|、 |0.681|、 |0.696|、 |0.678|、 |0.710|、 |0.710|、 |0.677|, 相关系数最大值平均提升67.93%。

CC、 sCARS、 CARS、 GA、 SPA算法筛选10个分解尺度的特征波段数量平均值分别占全波段数量的18.70%、 2.95%、 2.27%、 4.57%、 0.33%。 筛选的特征波段, CC算法主要分布于在400~1 500 nm可见光及近红外短波; sCARS、 CARS算法集中于1 500~2 500 nm近红外长波; GA算法筛选的特征波段在400~2 500 nm呈现均匀分布; SPA算法集中于760~2 500 nm近红外波段。

高光谱小波能量特征向量EFV可以较好估算表层土壤有机碳含量, 6种模型的训练集与验证集R2平均值分别为0.73、 0.74, RMSE平均值分别为7.64、 7.28, RPD平均值为1.95。 以CC-EFV-BPNN构建的模型效果最好, R2分别为0.78、 0.77, RMSE分别为6.38、 6.92, RPD为2.07; SPA-EFV-BPNN构建模型精度较差R2分别为0.70、 0.70, RMSE分别为8.44、 6.43, RPD为1.82。 对比模型验证集精度, CC-EFV-BPNN> sCARS-EFV-BPNN> Full-spectrum-EFV-BPNN> CARS-EFV-BPNN> GA-EFV-BPNN> SPA-EFV-BPNN。