土壤光谱数据来自数据库, 为ICRAF-ISRIC^[6]提供的可见-近红外漫反射光谱, 所用245份中国样本来自不同的土层深度, 参照FAO74标准^[7]分类如表2。 测量过程为标准化流程: 每份样品大约20 g, 经风干、 研磨后, 筛选2 mm细颗粒, 放在7.4 cm直径的玻璃器皿上, 样品的堆放厚度大约1 cm; 再由功率4.5 W、 色温3 000 K的卤素灯照明, 用FieldSpec FR(ASD)光谱仪记录350~2 500 nm的漫反射光谱; 为消除样品差异和仪器噪声, 探头每旋转90° 记录25条光谱的平均值, 每份样品读取前参考白板对光谱反射率进行校正。 测量结束后间隔10 nm对光谱反射率进行采样, 得到含216个数据点的反射率光谱曲线。 而土壤有机碳含量(SOC)通过测定使用的氧化剂重铬酸钾消耗量计算, 其测量与反演精度的单位为%, 探究有机碳含量与光谱反射率之间的关系。

首先, 对光谱反射率进行包络去除处理, 突出光谱特征, 使得测量、 组成的差异减轻, 不同土壤样品具有归一化光谱背景, 为多元回归做准备; 再者, 利用单位方差法(unit variance), 对模型输入的不同量纲的变量, 做均值中心化和归一化处理。 同时, 因土壤漫反射成分吸收的重叠严重, 不符合朗伯-比尔定律的使用条件, 故直接使用反射率光谱及其变换作为输入, 不进行吸光度变换。

另外, SOC参考值存在偏态分布, 采用幂函数或对数变换加以优化。 从试验设计的角度出发, 建模回归宜采用测量值范围较宽、 分布具有代表性的样本, 但碍于数据库所测样本有限, 频数分布不均衡, 故变换偏斜的数据反演。

1.3.1 双包络去除法

包络去除法^[8]是对光谱诊断特征的一种计算表达, 将反射率数值根据连结肩峰的“ 外壳(Hull)” 做归一化处理, 能够有效突出光谱特征。 利用连续折线段近似光谱曲线的包络, 可逐点计算对应波长λ 处的“ 包络去除值”

CRλ=RλCλ(1)

其中, R_λ 为波长λ 处的光谱反射率, C_λ 为波长λ 处上方的Hull值, 二者的比值量化了光谱吸收特征。 通常包络去除仅给出正向吸收特征, 欠缺对分段正向发射(或反向吸收)的CR_λ =1位置处的考虑, 一定程度上造成了变量信息缺失。 因此我们提出一种“ 双包络去除法” , 既保留原光谱反射率低谷处、 正向吸收的CR_λ , 可称“ 上包络去除值” (用CR或uCR表示), 也纳入光谱反射率凸起处、 反向下凹的“ 下包络去除值” (用CR'或dCR表示)。 具体地, 由原光谱得到反射率互补谱(1-R), 再根据ASD扫描间断位置或水的吸收峰将其分若干段, 对相应位置的包络进行上、 下包络去除, 捕捉CR=1处、 CR'_λ 随波长的变化, 反映反射率在上述区间、 速度变缓的过程, 实际是对斜率等一阶导数信息的补充。

如图1所示, 连结原光谱反射率的谷底位置, 以分段后下边的外壳为基础, 定义原光谱中各点对应的互补谱CR_λ 为新的“ 下包络去除值”

CR'λ=1-Rλ1-C'λ(2)

式(2)中, 下包络C'_λ 为由原光谱谷底位置确立的Hull值, CR'_λ 为反射率互补谱(1-R_λ )对应(1-C'_λ )归一化值, 也即原光谱在波长位置λ 处的下包络去除值。 建模时将得到的CR_λ 和CR'_λ 联合作为OPLS算法的输入, 加强对坡度变化的有效提取, 补充上包络去除信息的不足。 需注意的是, 下包络去除并非简单的发射包络去除, 而是结合了区间分段后的一种背景去除, 具体取值与断点位置有关。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 双边包络去除的过程图解Fig.1 Illustration of bi-continuum removal

1.3.2 Orthogonal PLS

PLSR^[9]是对高度共线性变量的多元回归方法, 能有效提取出与组分相关性最大的因子, 使之对响应矩阵具有较好的预测能力。 OPLS^[10]则基于正交信号校正(orthogonal signal correction, OSC)对PLSR进行了改进, 虽然精度无特别差异, 但其图谱可视化效果更好, 使模型变得简单、 易理解。 它将观测光谱分解为一个预测成分和多个正交成分; 前者与响应线性相关, 后者则正交无关, 也可通过交叉检验确定其数目, 有效避免PLSR模型反演结果的过度拟合。 为充分利用数据, 模型采用K折交叉验证^[11]。 最终, 在全谱建模的基础上, 再根据载荷散点图(loadings)分布、 回归系数(regression coefficients)和变量投影重要性(variable importance in projection, VIP)等进一步筛选特征波段, 针对某种土壤类型, 提出基于几个波长反射率的“ SOC指数” ^[12]。

1.3.3 模型评价指标

评价指标, 包括模型的拟合优度、 预测能力和结果的反演精度。 其中, 拟合优度由判定系数R²(coefficients of determination)反映模型能够解释的变化比例, 更进一步, R²指标又包括拟合优度R²Y、 预测优度Q², 二者越接近1、 模型的刻画能力越优; 反演精度则包括交叉验证根均方误差(root mean square error of cross validation, RMSE_CV)和估计根均方误差(root mean square error of estimation, RMSEE)较低。 兼顾Q²与R²较大时, 保证RMSE_CV与RMSEE越小, 才能说明模型反演得越优。 此外, 为评价模型参考值分布的偏斜状况, 还将引入偏度(skewness)衡量有机碳含量参考值分布的不对称性, 以此作为优化目标值分布的依据^[13], 调整样本的平均值与中位数关系, 确立参考值变量的变换。

关于模型的自变量, 图2中分别为光谱反射率(a), 经上包络去除后(b), 样品尺寸、 光照、 测量等差异减少, 吸收位置变得更加突出; 分段下包络去除(c)以反射率低谷的分段点, 将反射率凸起处的变化体现出来; (b)经UV缩放、 SOC参考值做正交信号校正后(d), 在1 230~1 280, 1 700~1 740和2 240~2 270 nm等波段具有差异, 与C— H, C=O等理论吸收峰相近^[14], 避开在1 400, 1 900和2 200 nm附近OH^-的影响, 为建立特征谱段SOC指数提供了可能性。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 (a): 土壤光谱反射率; (b): 光谱反射率的上包络去除; (c): 分段下包络去除; (d): 上包络去除的正交信号校正Fig.2 (a): Soil spectrum reflectivity; (b): Up-Continuum removal; (c): Down-continuum removal; (d): OSC of up-continuum removal

关于模型的参考值, 样本集的SOC含量并非正态分布, 而是整体向较低测量值集中、 偏斜。 偏度校正可改变分布的偏斜和不对称程度。 如图3, 原始SOC分布正偏, 偏度较高, 值为2.85; 经开方运算后, 因对小于1的数拉伸、 对大于1的数压缩, 偏度减小至1.22, 样本也更匀称地散落在中位值0.38两侧。 可见, 中位值贴近平均值, 最大值靠拢中位值, 改善了低含量样本数过多的情况, 提高了校正模型的判决能力和稳健性。 但偏度校正并非越低越好, 而应根据数据转换后参考值的区间范围、 集中程度确定: 性质类同的样本应以略微右偏为宜; 具有更复杂偏态分布的样本, 可考虑用系数待定的Log变换进行优化处理。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 (a): SOC参考值分布; (b): SOC开方变换后的参考值分布Fig.3 Histogram of (a): SOC in %; (b): Square root of SOC in 5%

关于建模性能, 表1给出了基于全谱反射率R、 包络去除量CR和双包络去除量BCR, 分别与PLS或OPLS进行组合, 反演SOC时不同模型的结构参数和性能指标。 (1)相比R-PLS模型, 进行上边包络去除, CR-PLS模型R²Y, RMSEE基本不变, Q²和RMSE_CV略有提高; (2)改用CR-OPLS模型, 各指标略微改善, 分别为0.71, 0.61, 0.43和0.50; (3)如基于双包络去除BCR, OPLS模型各项指标显著改善至0.84, 0.72, 0.33和0.42; (4)先用Log(Y+1)对数或Sqrt(Y)幂函数开方, 优化SOC含量参考值的分布, 则BCR-OPLS模型的性能进一步提高, R²和Q²达到0.90和0.79, RMSEE和RMSE_CV分别较提高65.4%和71.8%。 图4给出直接以有机碳含量Y、 光谱反射率R进行PLS拟合, 与SOC含量先开方(Y^1/2)、 再由双包络去除量BCR进行OPLS反演的结果, 用不同色块代表不同土壤类型, 比较预测值与参考值的拟合线, 显然后者的反演结果分布更密集、 拟合度更优, 表明BCR预处理、 OPLS模型和SOC分布优化具有联合作用。

表1

Table 1

表1(Table 1)

表1 不同方法、 变量集与样本分布下的模型反演性能比较 Table 1 Models’ performance under combinations of different methods, datasets and distribution Method X Variables Skewness Transform Factor R2Y Q2 RMSEE RMSECV PLS R 2.85 Y 13 0.69 0.54 0.45 0.55 PLS CR 2.85 Y 8 0.67 0.56 0.46 0.54 OPLS CR 2.85 Y 1+9+0 0.71 0.61 0.43 0.50 OPLS BCR 2.85 Y 1+14+0 0.84 0.72 0.33 0.42 OPLS BCR 1.35 Log(Y+1) 1+14+0 0.87 0.78 0.27 0.32 OPLS BCR 1.22 Sqrt(Y) 1+14+0 0.90 0.79 0.26 0.31

表1 不同方法、 变量集与样本分布下的模型反演性能比较 Table 1 Models’ performance under combinations of different methods, datasets and distribution

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 SOC含量参考值与R-PLS, BCR-OPLS模型反演结果的拟合线Fig.4 SOC content observedvs predicted by model (a): R-PLS; (b): BCR-OPLS

根据FAO74标准, 表2给出分类土壤的BCR-OPLS模型反演结果。 其中, 土壤类型由英文缩写表示, 线性和正交因子数用“ 1+N” 的形式表示, 并给出了分布偏度、 拟合优度和误差等指标。 容易看出, 对多数种类的土壤, R²Y分布在0.9~1.0, RMSEE集中于0~0.23, RMSE_CV基本落在0.1以下。 而反演不理想的个别土质类型, 例如“ 艳色淋溶土(Lc)” , 在使用对数函数优化参考集分布后, RMSE_CV则可从0.57降至0.34, R²从0.95升至0.98。 由此说明, 只要参考集的SOC含量分布合理, 对不同土壤建立的分类BCR-OPLS反演模型, 比此前建立的综合反演模型性能更优、 精度更高, 为推导基于特征波段的SOC指数提供基础。

表2

Table 2

表2(Table 2)

表2 不同土样分类OPLS模型性能汇总 Table 2 Performance of OPLS models of different soil types 分类OPLS模型 因子数 样本数 偏度 R2Y Q2 RMSEE RMSECV 艳色淋溶土(Lc) 1+2 21 2.1 0.92 0.74 0.23 0.57 潜育始成土(Bg) 1+3 5 2.03 1 1 0 0 潜育淋溶土(Lg) 1+4 6 0.96 1 1 0 0 风化粘磐土(Ne) 1+2 6 0.92 1 0.991 0.01 0.06 黄色铁铝土(Fo) 1+1 18 0.59 0.925 0.883 0.24 0.29 石灰粗骨土(Rc) 1+4 20 0.09 0.922 0.384 0.07 0.16 铁质淋溶土(Af) 1+6 43 0.89 0.954 0.882 0.13 0.2 强淋溶铁铝土(Fa) 1+3 6 1.14 1 0.991 0.01 0.04 典型强淋溶土(Ao) 1+1 9 0.63 0.92 0.743 0.04 0.06 不饱和始成土(Bd) 1+2 8 1.12 0.997 0.98 0.2 0.46 硫酸盐冲积土(Jt) 1+4 6 -0.05 1 1 0 0 典型黑土(Hh) 1+3 9 -0.71 0.998 0.979 0.03 0.06 典型碱土(So) 1+2 5 0.64 1 0.995 0.03 0.04 薄层黑钙土(Ch) 1+2 6 -0.63 0.999 0.99 0.03 0.08 石灰性始成土(Bk) 1+1 9 0.94 0.914 0.833 0.09 0.1

表2 不同土样分类OPLS模型性能汇总 Table 2 Performance of OPLS models of different soil types

为构建多项式或波段比值表示的“ SOC指数” , 本节基于BCR-OPLS模型的载荷位置关系, 利用变量投影趋势提取特征谱段。 载荷图描绘了X和Y变量关于OPLS预测项的空间。 理论上, BCR变量越靠近Y, 二者的线性关系愈加明确; 但实践中, 由于含量叠加非线性严重、 成分存在多谱段表达, 则对在预测轴投影接近Y、 投影趋势急剧反转的边缘区域, 也可看作潜在的特征谱段。 使用SIMCA通过OPLS分析模块得到载荷图, 根据投影位置和趋势反转提取特征波段, 在载荷图5(a)中, Y变量在预测成分轴投影值处在0.2附近, 图5(b)给出X变量uCR1650— uCR1750这一红色区域, 在预测方向上投影值靠近Y、 方向迅速改变且处于正半区边缘, 符合特征波段要求。 据此性质, 还可在载荷分布图里找到CR390— CR420, CR840— CR900和CR’ 610— CR’ 700等包络去除量作谱吸收的局部特征。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 (a): Y和BCR在载荷图中; (b): 对uCR1650— uCR1750局域放大Fig.5 (a): Y and BCR in loadings scatter plot; (b): Partial enlarged drawing of uCR1650— uCR1750 in (a)

考虑多种类型土壤, 建立基于特征谱段的综合反演模型。 以前述方法选择的四片红色区域的双包络去除量作X变量, 并由Lg(Y+1)优化SOC参考值分布作Y变量, 经单位方差法缩放, 进行BCR-OPLS模型反演。 如图6(a), 反演结果的拟合线R²为0.72, Q²为0.68, RMSE_CV小于0.48, 较全谱模型反演的性能有所下降。 再据图6(b)协相关系数大小和VIP准则, 还原出光谱反射率R₈₄₀~R₉₁₀, R_{1 610}~R_{1 760}和R₆₀₀~R₇₄₀, 建立基于特征谱段的R-OPLS模型, 其结果R²为0.61, Q²为0.58, RMSEE为0.56。 根据交叉判定系数Q²> 0.5, 由此推断特征谱段对多类土壤综合的SOC反演性能尚可。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 (a): 特征波段的反演结果; (b): 特征模型的协相关系数Fig.6 (a): Fit plot of feature based OPLS; (b): Coefficients of BCR-OPLS for all samples

然而, 协相关系数确立的判定系数仍显复杂, 无法利用特征区域得到简洁的SOC指数。 现针对特定土壤类型, 建立基于特征谱段的分类反演模型。 选取样本数较多、 全谱反演效果较好的黄色铁铝土(Fo), 根据投影趋势和位置关系, 确立下包络去除量CR’ 590作为特征谱段, 解析找到光谱反射率R₄₀₀, R₅₉₀和R₉₂₀作为模型的输入如图7(a), 经OPLS回归得到该类土壤的“ SOC指数”

ySOC=1.94+0.39R400+0.44R920-1.68R590(3)

其中, y_SOC为黄色铁铝土中有机碳含量。 图7(b)显示模型三个成分累积的R², Q²达0.9以上, 图7(c)中交叉验证精度达到较优的0.25, 表明针对特定类型土壤建立的BCR-OPLS模型, 可推得仅用3个波长表示的SOC指数, 为野外筛选探测以及多光谱遥感波段确立提供依据。 值得注意的是, 无论多种土壤综合反演模型推得的R₈₄₀~R₉₁₀, R_{1 610}~R_{1 760}和R₆₀₀~R₇₄₀, 还是分类土壤BCR-OPLS反演模型所用的CR’ 590特征区域, 这些谱段附近的吸收对模型的统计校正影响深刻, 与农作物、 植物残体及土壤里某些矿物质存在物理联系。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 (a): 黄色铁铝土(Fo)的SOC指数模型CR’ 590与原光谱反射率的位置关系; (b): 累积拟合优度R2与Q2演变; (c): 基于特征的SOC指数反演结果Fig.7 (a): SOC index model for Orthic Ferralsol with CR’ 590 and its spectral reflectance; (b): R2(cum) and Q2(cum) progression; (c): Fit plot of feature based SOC index