基于三维光谱指数的土壤含水率反演研究

周益安，周 昊，劳聪聪

（1.杭州市南排工程建设管理服务中心，杭州310019；2.浙江省水利河口研究院（浙江省海洋规划设计研究院），杭州310020；3.河海大学水利水电学院，南京210098）

0 引 言

土壤含水率（SMC）是土壤最重要的组成部分之一，影响着土壤的物理性质和土壤养分的溶解、转移，故而对土壤含水率进行监测具有重要意义。传统的土壤含水率监测方法主要有依靠土壤介电特性、土壤电导率、电磁波、中子法等[1,2]，但这些方法易受土壤物理化学特性的影响而产生局限性，无法满足实时、在线以及高精度测量的要求[3]。而高光谱技术能够高效且精准的性能[4,5]，弥补了传统方法的局限性。

高光谱技术凭借其高效、快速的优势在土壤含水率监测的研究中得到重视。以往的研究主要是利用光谱反射率进行直接反演，如彭翔等[6]，尹业彪等[7]和姚艳敏等[8]的研究。虽然这些研究对土壤含水率也具有较好的效果，但是模型具有多个自变量，相对复杂。Wang 等[9]发现利用光谱指数在预测土壤含水率方面是可行，而且利用光谱指数建模有效降低了模型的复杂程度。但是未探索三维光谱指数在土壤含水率预测方面的效果。已有研究表明在土壤有机质含量预测方面，三维光谱指数相比于一维二维光谱指数更有优势[10]。所以有必要研究三维光谱指数在土壤含水率预测方面的效果，为土壤含水率快速估算提供参考。

本文选取浙江永康地区的土壤样本进行研究，测定土样的光谱反射率，并对原始光谱反射率（R）进行一阶微分（FD）和二阶微分（SD）等处理，再构建不同维度的光谱指数，最后采用偏最小二乘回归模型(PLSR)分析不同光谱指数与土壤含水率之间的定量关系，以期找到对土壤含水率最敏感的光谱指数。

1 材料与方法

1.1 研究区概况

研究区域为浙江永康地区(28°45′N，119°53′E)。该地区海拔96 m，为亚热带季风气候，夏热冬温，四季分明，季风发达，年降水量介于800～1 600 mm 之间。该地区主要农作物为水稻，因此，有效监测土壤含水率，对该地区的精准农业发展是必要的。

1.2 研究材料采集及土壤样本含水率的测定

本次研究以永康地区采集的159个不同含水率的土样为研究样本，采样深度为0～20 cm。收集时剔除浸入体，然后将土样放入采样袋中密封、编号、称重并带回实验室。从采样袋内取20 g 左右有代表性的土样放入铝盒内，盖上盖称重，并记录铝盒的编号和重量。将没加盖的铝盒放入干燥箱内，在105 ℃、24 h 恒温条件下用干燥法[11]测得土壤样本质量含水率ωm。其中土壤样本质量含水率ωm的公式为：

式中：M1为土样质量（含铝盒）；M2为干燥后土样质量（含铝盒）；M3为空铝盒质量。

1.3 光谱测定

利用ASD 公司生产的FieldSpec 3 光谱仪（波长350～2 500 nm），在照明控制的暗室中测量土壤的光谱反射率。该光谱仪的采样间隔：350 到1 000 nm 为1.4 nm，1 000 到2 500 nm 为2 nm，重采样间隔为1 nm。将159 份制备好的土壤样品装入黑色容器（直径10cm，深度2cm），填充后将表面刮平。光谱测定的光源为50 W 卤素灯，使用视场角5°的光纤探头，根据洪永胜等[12]的研究将测量时距土壤样品表面距离定为50 cm，光源的天顶角定为50°，探头至待测样品的表面距离定为15 cm。每次光谱测定之前，去除暗电流并校准白板。每个土壤样品分别在4 个方向测量(3 次旋转，每次90°)，每个方向保存5 条光谱曲线，总共20 条。最后利用ViewSpecProV6.0.11 软件进行算术平均计算，得到土壤样品的实际反射光谱数据。

1.4 光谱数据预处理

为消除光谱测试环境、高频随机噪声和杂散光等干扰因素引起的噪声，提高数据的信噪比。故剔除信噪比较低的边缘波长350～400 nm 和2 401～2 500 nm，并利用Savitzky-Golay（SG）滤波法对所有光谱数据（400～2 400 nm）进行平滑去噪。为了降低光谱数据的冗余，将光谱数据按10 nm 的间隔重采样，得到400～2 400 nm，采样间隔为10 nm，共201 个波长数据用以后续的分析计算。此外，利用Matlab 软件计算光谱曲线的一阶微分（FD）及二阶微分(SD)，以此降低基线漂移的影响，突出曲线特征。

1.5 光谱指数构建

选取了六种光谱指数，包括3 个二维光谱指数(RI，DI，NDI)，3 个三维光谱指数（TBI1，TBI2及TBI3），用以分析比较不同维度指数在土壤含水率预测方面的效果[13]。具体公式如下：

式中：Rλ1，Rλ2和Rλ3为400～2 400 nm 中任意3 个波长，且Rλ1≠Rλ2≠Rλ3。

1.6 建模策略及模型评价

偏最小二乘回归模型（PLSR）集成了主成分分析、典型相关分析和多元线性回归等的优点。在建模过程中，降维、信息集成和波长优选等方法极大地提高了系统提取主成分的能力，因而在光谱数据建模中得到了广泛的应用[14,15]。本文选取偏最小二乘回归模型进行建模分析。

利用K-S 算法将159 个样本划分为两部分，其中106 个样本用于建模，其余53 个样本用以精度验证。土壤含水率的样本特征见图1。由图1 可知，3 个样本集的最大值及最小值接近，平均土壤含水率值均在23%左右。由图中3个小提琴子图的形状可以看出，3个样本集之间的样本分布基本一致，这为后面的分析建模提供了一定的数据基础。

图1 描述样本特征的小提琴图Fig.1 Violin diagram describing the characteristics of the sample

利用偏最小二乘回归模型建立土壤含水率的预测模型，通过对比各模型的建模均方根误差（RMSEc）、建模决定系数（）、预测均方根误差（RMSEv）、预测决定系数（Rv2）、相对分析误差（RPD），筛选出最优模型并对研究区土壤含水率进行预测。Rc 2，用以判定模型的稳定程度，越接近于1 说明模型的稳定性越好；RMSEc及RMSEp用于表征模型的准确性，其值越小表明模型的精度越高。另外，利用RPD对模型性能进行评估。当RPD＜1.5 时，模型几乎不能对样本进行预测；当1.5＜RPD＜2时，模型对样本可以粗略估计；当2＜RPD＜2.5时表明模型具有较好的定量预测能力；当2.5＜RPD＜3 时模型具有很好的预测能力；当RPD＞3.0 时表示模型具有极好的预测能力[16]。其中计算R2、RMSE及RPD的公式如下：

式中：yi和分别为检验样本的观测值和预测值；为样本观测值的平均值；n为预测样本数。

式中：yi和y^i分别为验证样本的观测值和预测值；n为预测样本数。

式中：S.D为样本观测值的方差；RMSE为均方根误差。

1.7 变量投影重要性分析

变量投影重要性分析（VIP）是基于PLSR 的一种变量重要性分析方法，通过计算VIP得分来表征各自变量对因变量的解释能力，并根据解释能力实现自变量的排序[17]，计算公式为：

式中：p是自变量的个数；m是主成分的总数；SSYf是第f个主成分解释方差的平方和；SSYtotal是因变量平方和；F是主成分总数；给出了第f个主成分中第j个变量的重要性。VIPj值越大，自变量重要性越高。一般认为，当自变量的VIP值大于1时，该自变量被判定为重要的自变量。

2 结果与分析

2.1 不同土壤含水率的光谱曲线分析

图2为不同预处理的土壤光谱曲线。由图2（a）可知，土壤光谱在1 400、1 900 和2 200 nm 等3 个波长处，均有明显的吸收谷，有研究表明这些吸收谷与土壤含水率有关[18]。经过一阶微分处理之后，见图2（b），可以发现光谱曲线水分吸收谷进一步凸显，也在一定程度上消除光谱曲线的基线漂移。而二阶微分[见图2（c）]虽然有效消除了基线漂移，但是同时也凸显了许多噪声，可能会对后面的建模分析造成一定的影响。

图2 不同预处理的土壤光谱曲线Fig.2 Soil spectral curves of different pretreatments

2.2 一维相关性分析

由图3可知，原始光谱反射率与SMC之间的敏感波长主要集中在1 400～1 600 nm 以及1 900～2 300 nm 区间。这个区间与土壤水分敏感区间一致。经过一阶微分处理之后，进一步加强对SMC 敏感的波长，削弱对SMC 敏感度不高的波长，出现“强者越强，弱者越弱”的现象。而二阶微分光谱更加突出了此现象，只有几个较强波长的R2较大，其余的波长的R2都相对较小，甚至接近于零（呈现紫色）。通过提取原始光谱反射率和一阶微分光谱及二阶微分光谱的最佳波长反射率作为一维光谱指数，用以建模分析。具体的指数见表1。

图3 一维决定系数矩阵图Fig.3 One-dimensional determination coefficient matrix diagram

2.3 二维相关性分析

通过Matlab 编程实现二维相关矩阵的计算，计算结果如图4 所示。由图4 可知，不同形式的二维光谱指数和不同光谱曲线均会影响R2的值，R-NDI 和R-RI 具有较好的效果，其R2可达0.88。基于二阶微分光谱的指数，R2普遍偏低，均小于0.8。光谱经过二阶微分处理后，其光谱反射率值明显减小，大部分波长都集中在0附近，而且曲线边缘震荡明显，凸显了噪声，如图4（c）所示。这些因素可能导致基于二阶微分的二维光谱指数效果较差。但二维光谱指数与一维光谱指数相比仍具有较大的提升，主要是通过在二维空间中搜寻波长的最佳组合，为指数的选择提供了更多的可能。通过各个矩阵中选取R2最大的指数作为二维光谱指数用以建模分析，具体见表1。

图4 二维决定系数矩阵图Fig.4 Two-dimensional determination coefficient matrix diagram

2.4 三维相关性分析

通过Matlab 编程实现三维相关矩阵的计算，计算结果如图5～图7所示。由这3个图可知，不同形式的三维光谱指数和不同光谱曲线均会影响R2的值，R-TBI1和FD-TBI1具有较好的效果，其R2可达0.89。除SD-TBI3较低以外，三维光谱指数的R2值都在0.84 以上。三维光谱指数与二维光谱指数相比具有一定的优势，说明将组成指数的波长扩展至3个具有一定的作用。本文选取三维相关矩阵中R2最高值的指数作为三维光谱指数用以后续的建模分析，具体见表2。

图5 TBI1的三维决定系数矩阵图Fig.5 Three-dimensional determination coefficient matrix diagram of TBI1

图6 TBI2的三维决定系数矩阵图Fig.6 Three-dimensional determination coefficient matrix diagram of TBI2

图7 TBI3的三维决定系数矩阵图Fig.7 Three-dimensional determination coefficient matrix diagram of TBI3

由表1可知，不同维度的光谱指数波长组成有所不同，但主要集中在1 400、1 900及2 100 nm 这几个对水分较为敏感的波长附近。三维光谱指数具有较高的R2可能综合考虑了3个波长组成效果。在选择指数时，主要通过统计方法，对指数的机理性解释还有一定的欠缺，即使指数具有较好的效果。

表1 不同光谱指数的波长组成Tab.1 Wavelength composition of different spectral indices

2.5 光谱指数的建模及分析

利用不同维度的指数分别建立PLSR 模型，具体结果见表2。由表2 知，基于TBI2的土壤含水率模型具有最佳的效果，其Rv2可达0.92，RPD可达3.32，说明模型具有精确预测的性能。与一维光谱指数的模型相比具有明显的优势，其RPD仅为2.28。但是和二维指数模型相比只具有一定的优势。

表2 基于光谱指数的偏最小二乘回归建模与预测效果Tab.2 Partial least squares regression modeling and prediction effect based on spectral index

为了更加直观的分析各个模型的预测性能，绘制图8基于不同模型预测值与实测值的散点图。由图8可以发现，TBI2[见图8（f）]的预测带最窄，说明其预测的精度最高，推荐TBI3模型作为所研究区域的土壤含水率预测模型。

图8 基于不同光谱指数模型的预测精度对比Fig.8 Comparison of prediction accuracy based on different spectral index models

2.6 变量重要性分析

为了确定各个指数的相对重要性，计算了各个指数的VIP值并绘制成图9。由图9 可以发现三维光谱指数只有两个在1以下，而二维光谱指数有4 个在1 以下，一维光谱指数均在1以下，说明三维光谱指数重要性高于二维及一维光谱指数。且VIP 最高的为R-TBI1，这与最佳土壤含水率预测模型（TBI2）的建模指数具有一定的偏差。但两者都是三维光谱指数，说明三维光谱指数具有一定的优势，但是在以后的研究中仍需进一步探索更适用于土壤含水率的三维光谱指数形式。

图9 不同光谱指数的VIP图Fig.9 VIP diagrams of different spectral indices

3 结 论

以浙江永康地区土壤为研究对象，使用159个样品的土壤含水率数据和相应的高光谱数据，用不同维度的光谱指数分别构建土壤含水率偏最小二乘回归模型，通过模型对比分析，得到以下结论：

（1）三维光谱指数相比于一维二维光谱指数对土壤含水率更加敏感。利用三维光谱指数构建土壤含水率预测模型更有利。

（2）基于TBI2 的偏最小二乘回归模型具有最佳的预测效果，其RPD可达3.32。

（3）对于土壤含水率预测模型而言，R-TBI1 相比于其他的光谱指数更具重要性，其VIP 值达1.04。基于三维光谱指数的土壤含水率模型的构建，对遥感精准监测土壤墒情具有重要的意义。