超级会员免费看
光谱转移与变换：去除土壤水分等外部影响
1. 光谱校正与土壤水分影响
在光谱分析中，不同仪器获取的光谱可能存在差异，需要进行校正。同时，土壤的红外光谱对外部环境条件（如土壤水分）非常敏感。土壤的吸光度光谱会随着水分含量的增加而增大。在实验室中，土壤通常在标准风干条件下进行扫描，但在野外很难控制土壤的含水量。尽管许多研究成功地使用野外采集的光谱来校准土壤性质（如土壤有机碳，SOC），但土壤水分含量的变化确实会对土壤性质的预测产生重要影响。
1.1 光谱校正示例
这里展示了参考仪器的光谱和通过分段直接标准化（PDS）校正后的从属仪器光谱。需要注意的是，由于对光谱应用了移动窗口操作，应去除光谱的边缘部分。
1.2 土壤水分对光谱的影响
为了研究土壤水分含量对光谱的影响，我们使用了包含 100 个土壤样本（来自更大数据集的子集）在 3 种不同水分条件下的光谱数据：
-
spectra0
：土壤在风干条件下的吸光度光谱（平均水分 5%）。
-
spectra1
：土壤湿润后的吸光度光谱（平均水分 12%）。
-
spectra2
：湿润土壤风干 1 天后的吸光度光谱（平均水分 9%）。
以下是加载和处理数据的代码：
# load the required packages
require(prospectr)
require(soilspec)
# load the package data
data("datsoilspc")
# convert reflectance to absorbance
spectraA <- log(1/datsoilspc$spc)
# embed the soil property and the spectra in one single table
datsoilspc$spcA <- spectraA
# apply a Savitzky-Golay filter
datsoilspc$spcASg <- savitzkyGolay(datsoilspc$spcA, w = 11, p = 2, m = 0)
# apply some smoothing to the spectra
old.wavs <- as.numeric(colnames(datsoilspc$spcASg))
new.wavs <- seq(500, 2450, by = 5)
datsoilspc$spcARs <- prospectr::resample(datsoilspc$spcASg,
wav = old.wavs,
new.wav = new.wavs,
interpol = "linear")
# apply a standard normal variate transformation for baseline correction
datsoilspc$spcASnv <- standardNormalVariate(datsoilspc$spcARs)
# import the data
data("datEPO")
# show names
names(datEPO)
# extract the objects from datEPO
soilC <- datEPO$soilC
spectra0 <- datEPO$spectra0
spectra1 <- datEPO$spectra1
spectra2 <- datEPO$spectra2
# apply Savitzky-Golay filter
spectra0_sg <- savitzkyGolay(spectra0, w = 11, p = 2, m = 0)
spectra1_sg <- savitzkyGolay(spectra1, w = 11, p = 2, m = 0)
spectra2_sg <- savitzkyGolay(spectra2, w = 11, p = 2, m = 0)
# make some preprocessing to the spectra
spectra0_rs <- prospectr::resample(spectra0_sg,
wav = old.wavs,
new.wav = new.wavs,
interpol = "linear")
spectra1_rs <- prospectr::resample(spectra1_sg,
wav = old.wavs,
new.wav = new.wavs,
interpol = "linear")
spectra2_rs <- prospectr::resample(spectra2_sg,
wav = old.wavs,
new.wav = new.wavs,
interpol = "linear")
# apply a standard normal variate transformation for baseline correction
spectra0Snv <- standardNormalVariate(spectra0_rs)
spectra1Snv <- standardNormalVariate(spectra1_rs)
spectra2Snv <- standardNormalVariate(spectra2_rs)
1.3 不同水分条件下光谱的可视化
我们可以绘制样本 1 在三种不同水分条件下的光谱，以直观地观察水分含量对光谱的影响。
# plot the spectra0
plot(colnames(spectra0Snv), spectra0Snv[1, ],
type = "l",
xlab = "Wavelength /nm",
ylab = "Absorbance",
col = "blue")
# add line for spectra2
lines(colnames(spectra2Snv), spectra2Snv[1, ],
col = "red")
# add line for spectra1
lines(colnames(spectra1Snv), spectra1Snv[1, ],
col = "green")
# add a legend
legend("topright",
legend = c("average moisture 5%", "average moisture 9%", "average moisture 12%"),
col = c("blue","red" ,"green"),
lty = 1,
cex = 1)
1.4 土壤水分对预测的影响
通常，我们会基于干燥样本的光谱库建立模型，然后使用该模型从野外采集的光谱（水分含量变化）中预测土壤性质。为了评估水分含量对预测的影响，我们进行了以下操作：
# load required package
require(Cubist)
# compute the integer number corresponding to the 25% of the samples
nsamples2select <- round(x = nrow(datsoilspc$spcASnv) * 0.75, digits = 0)
# separate between calibration and validation
isrow <- sample(1:nrow(datsoilspc$spcASnv), size = nsamples2select)
datsoilspcC <- datsoilspc[isrow,]
datsoilspcV <- datsoilspc[-isrow,]
# generate a Cubist model
soilcCubistModel <- cubist(x = datsoilspcC$spcASnv, y = datsoilspcC$TotalCarbon)
# predict on the calibration dataset
soilvCubistPred <- predict(soilcCubistModel, datsoilspcV$spcASnv)
# plot the predicted calibration data
plot(datsoilspcV$TotalCarbon, soilvCubistPred,
xlab = "Observed",
ylab = "Predicted",
xlim = c(0, 12),
ylim = c(0, 12),
pch = 16)
abline(0, 1)
我们使用相关函数评估了预测结果，结果如下表所示：
| 数据集 | ME | RMSE | r2 | R2 | rhoC | RPD | RPIQ |
| ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- |
| 校准集 | -0.09 | 0.56 | 0.86 | 0.85 | 0.91 | 2.6 | 2.12 |
| 低水分（5%） | -0.25 | 0.79 | 0.79 | 0.84 | 0.91 | 2.54 | 1.62 |
| 中水分（9%） | 0.31 | 1.08 | 0.64 | 0.71 | 0.84 | 1.86 | 1.18 |
| 高水分（12%） | 1.81 | 2.61 | 0.14 | -0.72 | 0.16 | 0.77 | 0.49 |
从结果可以看出，干燥样本的预测较为准确，而当土壤样本含水量较高时，碳含量会被低估。
1.5 流程总结
以下是上述操作的流程图：
graph LR
A[加载数据] --> B[数据预处理]
B --> C[建立模型]
C --> D[模型预测]
D --> E[评估预测结果]
2. 外部参数正交化（EPO）方法
为了去除光谱校准中水分的影响，我们可以使用外部参数正交化（EPO）方法。EPO 算法将所有土壤光谱投影到与不需要的变化（即水分）空间正交的方向，从而有效地去除由于土壤水分引起的变化。
2.1 EPO 原理
在 EPO 中，假设光谱中包含的信息可以分解为三个部分：
-
有用的化学光谱响应
：也称为光谱生色团。
-
由外部因素引起的部分
：与第一部分无关。
-
残差
。
用矩阵形式表示，光谱 X（大小为 n × b）可以写成：
[X = XP + XQ + R]
其中，P 是光谱有用部分的投影矩阵（大小为 b × b），(\tilde{X} = XP)；Q 是光谱中不需要部分（即受水分影响的部分）的投影矩阵（大小为 b × m），(X_T = XQ)；R 是大小为 n × b 的残差矩阵。EPO 的目的是获得光谱的有用部分：
[\tilde{X} = X(I - Q)]
矩阵 Q 可以写成 (Q = GG^T)。矩阵 G 是光谱中与有用部分正交的无信息部分；可以使用差分光谱 D 的主成分作为近似。这里，我们将 D 定义为湿润和风干（标准）条件下光谱的差异。
2.2 EPO 实现步骤
2.2.1 计算差分矩阵 D
# difference matrix (spectra difference between wet and dry)
D = as.matrix(spectra0Snv - spectra1Snv)
2.2.2 定义 EPO 函数
# define the EPO function
epo <- function(D, npc){
# npc is the number of components to use
# return: P: the projection matrix
D <- as.matrix(D)
n <- nrow(D)
p <- ncol(D)
dtd <- t(D) %*% (D)
# singular value decomposition of the D (n x n) matrix
s <- svd(dtd)
# extract the no. factors
ld <- s$v[, 1:npc]
# projection matrix
P <- diag(p) - ld %*% t(ld)
return(P)
}
2.2.3 计算投影矩阵 P
# define number of principal components
npc <- 3
# compute the projection matrix P:
P <- epo(D, npc)
2.2.4 可视化投影矩阵 P
# load required package
require(soilspec)
# plot the projection matrix P
myImagePlot(P,
zlim = c(-0.1,0.1))
2.2.5 投影光谱
# EPO projected spectra of spec0
Z0 <- as.matrix(spectra0Snv) %*% P
colnames(Z0) <- colnames(spectra0Snv)
# EPO projected spectra of spec1
Z1 <- as.matrix(spectra1Snv) %*% P
colnames(Z1) <- colnames(spectra1Snv)
# EPO projected spectra of spec2
Z2 <- as.matrix(spectra2Snv) %*% P
colnames(Z2) <- colnames(spectra2Snv)
2.2.6 绘制校正后的光谱
# plot corrected spectra0
plot(colnames(spectra0Snv), Z0[1,],
type = "l",
xlab = "Wavelength /nm",
ylab = "Absorbance",
col = "blue")
# add line for corrected spectra2
lines(colnames(spectra2Snv), Z2[1,],
col = "red")
# add line for corrected spectra1
lines(colnames(spectra1Snv), Z1[1,],
col = "green")
# add a legend
legend("topright",
legend = c("average moisture 5%", "average moisture 9%", "average moisture 12%"),
col = c("blue","red" ,"green"),
lty = 1,
cex = 1)
2.3 EPO 处理后的模型预测
# project the whole calibration spectra
specZ <- datsoilspc$spcASnv %*% P
colnames(specZ) <- colnames(datsoilspc$spcASnv)
specZC <- specZ[isrow,]
specZV <- specZ[-isrow,]
# make a Cubist model from EPO transformed dry spectra
epoCubistModel <- cubist(x = specZC, y = datsoilspcC$TotalCarbon)
# predict on the calibration dataset
soilvCubistPred <- predict(epoCubistModel, specZV)
# plot the predicted calibration data
plot(datsoilspcV$TotalCarbon, soilvCubistPred,
xlab = "Observed",
ylab = "Predicted",
xlim = c(0, 12),
ylim = c(0, 12),
pch = 16)
abline(0, 1)
我们同样评估了 EPO 处理后不同水分条件下的预测结果，结果如下表所示：
| 数据集 | ME | RMSE | r2 | R2 | rhoC | RPD | RPIQ |
| ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- |
| 低水分（5%） | 0.4 | 0.82 | 0.64 | 0.83 | 0.91 | 2.43 | 1.55 |
| 中水分（9%） | 1.29 | 1.83 | 0.5 | 0.15 | 0.66 | 1.09 | 0.69 |
| 高水分（12%） | 0.26 | 1.83 | 0.23 | 0.16 | 0.42 | 1.09 | 0.69 |
从结果可以看出，去除水分的影响显著提高了预测精度。
2.4 EPO 流程总结
以下是 EPO 方法的操作步骤列表：
1.
EPO 校准
：
- 计算差分光谱 D 作为湿润和干燥光谱之间的差异。
- 使用 EPO 函数对 D 进行主成分分析（PCA），定义 c 个因子。
- 从 epo() 函数中获得投影矩阵 P。
2.
模型校准
：
- 将光谱库中的光谱 X 转换到不受水分影响的 EPO 空间（Z），(Z = XP)。
- 使用转换后的光谱 Z 校准模型以预测土壤性质。
3.
模型预测
：
- 转换未知样本的光谱 x：(z = xP)。
- 使用 z 在已校准的模型中预测土壤性质。
2.5 EPO 方法流程图
graph LR
A[计算差分矩阵 D] --> B[定义 EPO 函数]
B --> C[计算投影矩阵 P]
C --> D[投影光谱]
D --> E[建立模型]
E --> F[模型预测]
F --> G[评估预测结果]
通过以上方法，我们可以有效地去除土壤水分等外部因素对光谱的影响，提高土壤性质预测的准确性。
3. 基于威尔克斯 Lambda 确定主成分数量
在使用 EPO 方法时，需要选择的唯一参数是主成分的维度数量 npc，它代表需要去除的不需要部分。确定 npc 的一种方法是通过对基于转换后光谱的校准模型进行交叉验证，但这可能因不同的土壤性质模型而异。另一种方法是使用威尔克斯 Lambda（Wilks’ Λ）。
3.1 威尔克斯 Lambda 原理
威尔克斯 Lambda 是一种聚类分离度量，它是组间方差与总方差的比率。在这种情况下，(\Lambda = tr(B)/tr(T))，其中 (tr(.)) 是矩阵的迹，B 是聚合转换后光谱（即所有水分水平上 EPO 转换后光谱的平均值）的组间方差 - 协方差矩阵，T 是 EPO 转换后光谱的方差 - 协方差矩阵。
在 EPO 转换之前，样本在不同水分含量下的光谱（样本内变化）比两个不同样本的光谱（样本间变化）差异更大。EPO 转换后，不同样本应在光谱空间中得到良好分离。威尔克斯 Lambda 等于 1 表示样本完全分离，值为 0 表示没有分离。
3.2 计算威尔克斯 Lambda
以下代码计算威尔克斯 Lambda 作为主成分数量的函数，我们将 npc = 0 分配为不进行转换：
D = as.matrix(spectra0Snv - spectra1Snv)
SC <- (spectra0Snv + spectra1Snv)/2
B <- cov(SC)
TrB <- sum(diag(B))
T <- cov(spectra1Snv)
nc <- 20
xc <- seq(0:nc)
wilks <- matrix(0, nrow = nc + 1, ncol = 1) # create
# no transformation
wilks[1] <- TrB/sum(diag(T))
# make a for loop to estimate the Wilks" lambda
for (i in 1:nc) {
npc <- i
P <- epo(D, npc)
Z0 <- as.matrix(spectra0Snv) %*% P
Z1 <- as.matrix(spectra1Snv) %*% P
T <- cov(Z1)
SC <- (Z0+Z1)/2
B <- cov(SC)
TrB <- sum(diag(B))
wilks[i+1] <- TrB/sum(diag(T))
}
# plot the Wilks" lambdas against number of component
plot(xc, wilks,
type = "l",
xlab = "Number of components",
ylab = expression(paste("Wilks", ~Lambda)))
3.3 结果分析
通过绘制威尔克斯 Lambda 与主成分数量的关系图，我们可以找到威尔克斯 Lambda 的第一个最大值。在这个例子中，威尔克斯 Lambda 的第一个最大值出现在 npc = 3，这表明之前使用的三个主成分是一个不错的选择。
3.4 威尔克斯 Lambda 计算流程
以下是计算威尔克斯 Lambda 的操作步骤列表：
1. 计算差分矩阵 D。
2. 计算平均光谱 SC。
3. 计算组间方差 - 协方差矩阵 B 和总方差 - 协方差矩阵 T。
4. 初始化主成分数量序列 xc 和威尔克斯 Lambda 矩阵 wilks。
5. 对于每个主成分数量 npc：
- 计算投影矩阵 P。
- 投影光谱得到 Z0 和 Z1。
- 更新 T 和 B。
- 计算威尔克斯 Lambda 并存储在 wilks 中。
6. 绘制威尔克斯 Lambda 与主成分数量的关系图。
3.5 威尔克斯 Lambda 计算流程图
graph LR
A[计算差分矩阵 D] --> B[计算平均光谱 SC]
B --> C[计算 B 和 T]
C --> D[初始化 xc 和 wilks]
D --> E[循环计算威尔克斯 Lambda]
E --> F[绘制关系图]
4. 其他校正方法说明
虽然 DS 和 PDS 方法也可用于校正影响光谱的其他外部因素，如土壤水分，但不建议使用它们来处理野外土壤水分可变的情况。因为野外土壤的水分含量是变化的，而不是固定在一个水分含量上。不过，有研究表明 DS 方法在处理田间水稻土的水分问题时是有效的（假设田间土壤始终是湿润的）。
5. 总结
5.1 主要内容回顾
- 土壤水分会显著影响光谱和土壤性质的预测，干燥样本的预测较为准确，而高水分样本的碳含量会被低估。
- 外部参数正交化（EPO）方法可以有效地去除土壤水分等外部因素对光谱的影响，提高预测精度。
- 可以通过计算威尔克斯 Lambda 来确定 EPO 方法中主成分的数量，以优化 EPO 转换效果。
5.2 方法对比
| 方法 | 适用情况 | 优点 | 缺点 |
|---|---|---|---|
| EPO | 去除土壤水分等外部因素影响 | 有效去除水分影响，提高预测精度 | 需要确定主成分数量 |
| DS 和 PDS | 处理固定水分含量的情况 | 可校正外部因素 | 不适合野外水分可变情况 |
5.3 未来展望
在实际应用中，可以进一步研究如何更准确地确定 EPO 方法中的主成分数量，以提高预测的准确性。同时，可以探索将多种校正方法结合使用，以更好地应对不同的外部因素和复杂的土壤环境。
通过本文介绍的方法和技术，我们可以更好地处理光谱数据，去除外部因素的干扰，从而更准确地预测土壤性质，为土壤科学研究和农业生产提供有力的支持。
扫一扫
871
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
