基于图像处理和GBRT模型的表土层土壤容重预测

杨 玮 兰 红 李民赞 孟 超

(中国农业大学现代精细农业系统集成研究教育部重点实验室， 北京 100083)

0 引言

土壤容重又称为土壤密度，是指在自然结构状态下单位体积干土的质量，是表征土壤物理性质的重要参数[1]。表土层土壤容重一般指地表以下28 cm深度范围内土层的容重，是影响土壤水分、盐分及养分在土壤中运移的重要因素之一[2]。土壤容重的测量方法可分为直接方法和间接方法。传统的测量方法是环刀法，属于直接测量法的一种[3]。该方法是测量土壤容重的标准方法，测量过程简单，所得结果较为精确，但是对于大范围农田土壤容重数据的获取，该方法耗时、耗力，不能满足精细农业快速测量土壤容重的需求。

近年来，国内外学者在快速高效预测土壤容重方面进行了较多研究。文献[4]以圆锥指数作为输入，应用BP神经网络建立土壤传递函数PTFs，预测了粘土容重和粉质壤土容重，决定系数R2分别为0.697 3和0.686 8。文献[5]对兰陵溪流域森林土壤数据库进行研究，利用土壤有机碳、有机质建立了与土壤容重之间的回归模型。文献[6]研究了6种不同的回归方法，将几何平均粒径、土壤容重和粘土含量作为模型的输入，预测土壤水容重，结果显示，神经模糊和基因表达编程模型具有较高的预测精度，R2分别为0.874和0.995。文献[7]将预测的土壤含水率作为辅助数据，改进协同克里格算法，实现了对土壤容重的预测，预测精度R2从0.449提高到0.852。

然而，影响表土层土壤容重的因素是多方面的，在大多数预测土壤容重的研究中，影响因素的选取较为简单，且部分参数需要经过复杂的实验才能获得。另外，由于数据获取地块和建模数据不同，大部分由实验建立的预测模型的适用范围具有针对性和局限性，因此应用于大范围农田土壤容重的获取时存在困难，不利于方法的推广使用。

现有研究表明，表土层土壤容重与土壤孔隙度存在极显著负相关关系[8]，土壤表面孔隙度可以通过土壤表面粗糙度来预测[9-11]。而图像处理在粗糙度测量中有着广泛应用[12-13]。土壤阻力和土壤容重也具有着重要关系[14-15]。为此，本文从土壤阻力和基于图像处理得到的表征土壤表面粗糙度的信息融合入手，避免复杂的参数获取实验，建立梯度提升决策回归树(Gradient boosting regression tree，GBRT)模型，以期实现利用常规易获得的土壤物理参数对表土层土壤容重进行预测。

1 数据获取

1.1 数据采集

田间实验于2019年5月在北京市海淀区中国农业大学上庄实验站进行，该站位于116.189 4°E，40.141 1°N。实验田长42 m，宽40 m，总占地面积1 680 m2。采集10列样本，每列采集10个，共100个样本采集点。首先使用手机自带的相机拍摄采样点表面的土壤图像，然后使用环刀取土法[16]采集土样并装入有序号标记的铝盒中，采土设备如图1所示。同时使用GPS测量仪记录采样点的GPS数据，从第1个采集点开始，按一条龙顺序依次采集100个土样。使用实验室自主研发的车载式土壤阻力测量系统[17]测得离地表约15 cm深处的土壤阻力和GPS信息并保存，测量场景如图2所示。为了减少天气等因素对采集样本的影响，数据采集需要保证在24 h之内完成。

图1 环刀取土设备Fig.1 Equipment for digging soil

图2 车载式土壤阻力测量场景Fig.2 Vehicle-mounted soil resistance measurement

1.2 数据处理

车载式土壤阻力测量系统采用应变片电桥法测得深松犁钩在水平方向上受到的阻力并记录GPS信息。采集的数据为连续数据，所以根据土样采集点的GPS数据在土壤阻力数据中找到对应采集点的土壤阻力信息。将采集的100个土样带回实验室，连同环刀和铝盒一起称量。然后把装有环刀和土壤的铝盒放入105℃的电热恒温干燥箱中干燥24 h，取出铝盒，冷却至室温，再次称量并记录数据。表土层土壤容重计算式为

(1)

式中d——土壤容重,g/cm3

P——环刀内湿土质量,g

V——环刀容积，cm3

W——土壤含水率,%

剔除实验中的人为误差和损坏数据，一共取得90组可用于建模分析的实验数据。选取其中75%的数据四舍五入得到68个样本进行建模，剩下的22个样本作为验证集。

2 研究方法

2.1 图像预处理

考虑到拍摄时偏光导致图像照明不均匀，图像各部分的平均亮(灰)度产生起伏变化，同时为了提高图像的对比度，起到图像增强的作用，本文使用同态滤波技术[18]对土壤表面图像进行预处理。主要步骤如下：

使用fi(x,y)表示照明函数，fγ(x,y)表示物体的反射函数，则在光照下获得景物图像的数学模型为

f(x,y)=fi(x,y)fγ(x,y)
(0≤fi(x,y)≤∞,0≤fγ(x,y)≤1)

(2)

由于照明函数fi(x,y)在空间上变化缓慢，其频谱处于低频区域；反射函数fγ(x,y)反映图像的细节内容，且纹理本身具有较多的细节和边缘，其频谱处于高频区域，故将空间域转变为频域。首先，对式(2)两端取对数，使原来相乘的2个分量变成空间域中相加的2个分量，即

z(x,y)=lnf(x,y)=lnfi(x,y)+lnfγ(x,y)

(3)

再对式(3)两端分别进行傅里叶变换，得

Z(m,n)=I(m,n)+R(m,n)

(4)

式中，Z(m,n)、I(m,n)和R(m,n)分别为z(x,y)、lnfi(x,y)和lnfγ(m,n)的傅里叶变换。

将式(4)两端乘以传递函数H(m,n)，构成同态滤波器，使低频段受压缩，高频段得以扩展。则

S(m,n)=H(m,n)Z(m,n)=
H(m,n)I(m,n)+H(m,n)R(m,n)

(5)

利用傅里叶反变换，把式(5)变换回空间域，得

s′(x,y)=f′i(x,y)+f′γ(x,y)

(6)

由于傅里叶变换前取了对数，所以同态滤波增强后的图像函数g(x,y)为

g(x,y)=exp(s′(x,y))=
exp(f′i(x,y))exp(f′γ(x,y))

(7)

亦即

g(x,y)=i0(x,y)γ0(x,y)

(8)

式中，i0(x,y)和γ0(x,y)分别为输出图像的照射分量和反射分量。

2.2 表征土壤表面粗糙度的特征参数提取

特征提取是从图像中提取一组反映图像特性的基本元素或数字值[19]。表征图像表面粗糙度轮廓的特征向量可由反映纹理粗细程度的像元特征参数和反映纹理排列结构的区域特征参数两部分组成[20]。像元特征参数主要由图像灰度概率密度估计，通常由图像的灰度直方图来反映灰度概率密度，由熵、平均值、方差、偏度和峰度来表征直方图的分布情况，从而反映纹理的细致程度；区域特征参数主要由图像的灰度共生矩阵来反映。灰度共生矩阵通过研究图像中任意两点之间的灰度空间相关性来描述纹理，反映了关于方向、间隔和变换幅度的综合信息。通常选其能量、熵、对比度和逆方差等度量指标来描述图像的纹理特征。

2.3 特征参数选择

特征参数选择是从已经抽取的特征中选择能够更好地完成图像识别任务的参数，来表示原图像。实验使用灰度关联分析法[21]计算图像灰度直方图的熵、平均值、方差、偏度和峰度，以及图像灰度共生矩阵的能量、熵、对比度和逆方差，再加上土壤阻力共10种参数分别与表土层土壤容重的关联度，选择关联度较高的参数作为梯度提升决策回归树模型的输入变量。

灰度关联分析法的基本思想是量化特征图的几何状态，再计算出参考序列和比较序列之间的关联度。关联度较大时，比较序列的几何发展情况与参考序列更为接近[22]。利用灰色关联理论求出灰色关联度，以此来描述因素间关系的强弱和次序，为合理地选择回归模型的输入因子提供依据。

首先将表土层土壤容重作为母序列Y，将灰度直方图的熵、平均值、方差、偏度、峰度以及图像灰度共生矩阵的能量、熵、对比度、逆方差和土壤阻力作为子因素X1、X2、X3、X4、X5、X6、X7、X8、X9、X10，建立分析序列为

Y=(y1,y2,…,yk,…,yn)

(9)

(10)

式中m——样本数量

n——子因素数量

k——样本组序号

xnm——第m个样本的第n个子因素

然后对所选样本序列进行无量纲化处理

(11)

(12)

最后计算关联系数，R=[ξij]m×n，其中

(13)

式中ξij——关联系数

R——关联系数的m×n矩阵

ρ——分辨系数，一般取0.5

各因素关联度计算公式为

(14)

式中ψ0j——关联度

2.4 梯度提升决策回归树

梯度提升决策回归树(GBRT)是由FRIEDMAN[23]于2001年提出的一种函数空间优化算法。该算法有2个重要思想，分别为提升和梯度。提升，即对一个任务而言，先由一棵决策树对任务进行判断，再用下一棵决策树对先前过程中的错误进行判断，并重复循环该过程直到得到满意的结果为止。梯度，即利用梯度作为前一步中的预测误差，梯度容易获取，而一般的损失函数在计算每一次迭代过程中的增益时难以获取，将梯度代替每一步迭代过程中的损失下降，并作为下一步拟合的目标值，易于操作和实现。主要步骤为：

(1)确定学习目标

在GBRT模型构建过程中，通常采用平方误差作为模型的学习目标，计算公式为

(15)

式中yi——真实值

f(x)——模型预测值

N——样本的数量

(2)构建弱决策回归树

GBRT是以分类回归树作为基函数的前向分布算法与加法模型。构建弱决策回归树的过程即为构建一棵分类回归树的过程。构建过程分为决策树生成和决策树剪枝。决策树生成是指基于训练集数据来生成决策树，生成的决策树越大越好。决策树剪枝是指利用测试集数据，以最小化评估指标对已生成的树进行剪枝。

(3)GBRT的集成

模型最终的预测结果是将所有弱决策树的预测结果进行一定比例的求和。

3 结果与分析

分别对90组土壤表面图像利用同态滤波技术进行预处理。原始RGB图像颜色本身容易受到光照的影响，不能很好地反映图像的本质信息。图像所携带的信息主要表现在灰度形式上，颜色对提取信息并不是必需的。而且3通道处理转换为单通道处理以后，运算量可以大大减少，所以需要先将原始图像转换为灰度图像。部分处理后的图像与原始灰度图像对比如图3所示。可以看到，与图3a原始灰度图像相比，图3b经过同态滤波处理后的图像，其凹凸处的明暗对比更加明显，达到了图像增强的效果。

图3 原始灰度图像与同态滤波处理后的图像Fig.3 Original grayscale image and image after homomorphic filtering

对每幅经同态滤波技术处理后的图像进行表征土壤表面粗糙度的特征参数提取。计算图像灰度直方图的熵(设为熵1)、平均值、方差、偏度、峰度和图像灰度共生矩阵的能量、熵(设为熵2)、对比度、逆方差的结果如表1所示。

使用灰度关联分析进行特征参数的提取。计算90组实验数据的图像灰度直方图的熵、平均值、方差、偏度和峰度，图像灰度共生矩阵的能量、熵、对比度和逆方差，以及土壤阻力共10个子因素与表土层土壤容重之间的关联度，结果如表2所示。关联度从大到小排序为X10、X3、X4、X8、X5、X7、X2、X9、X6、X1。

一般关联度大于等于0.8时，子序列与母序列关联度很好；处于0.6～0.8之间关联度较好；小于0.6时，表示该子序列与母序列基本不相关。从表2可知，除图像灰度直方图的熵以外，其他的子序列与表土层土壤容重的关联度都大于0.6。考虑到模型输入变量的选择对构建出的模型精度的影响，最终选定关联度大于0.65的变量，即土壤阻力，灰度直方图的方差、偏度、峰度，图像灰度共生矩阵的熵和对比度，共6个子因素作为梯度提升决策回归树的输入变量。计算6个子因素之间的相关系数，结果如表3所示。可以看到变量之间的相关性都比较低，因此可以作为模型的输入参数。

表2 表土层土壤容重与各子因素关联度计算结果Tab.2 Calculation results of correlation between soil bulk density of top soil layer and each sub-factor

GBRT模型主要有7个参数需要调整：弱决策树的数量n_estimators，树深度max_depth、内部节点再划分所需最小样本数min_samples_split，叶子节点最少样本数min_samples_leaf、学习率learning_rate，最大特征数max_features和子采样比例subsample。模型的构建和测试均在Python 3.7的环境中完成。首先选取一个固定的学习率，利用训练样本来训练模型，根据训练样本的损失来确定弱决策树的数量。通常对训练集进行交叉验证[24-25]获取，持续增大决策树的数量直到交叉验证的结果不再降低为止，此时获取的数量就是在确定学习率下能够获取的最优决策树数量。本实验以学习率0.1进行5折交叉验证，验证集的平均绝对误差(MAE)随决策树的数量变化趋势如图4所示。

图4 交叉验证下平均绝对误差随决策树数量的变化趋势Fig.4 Changing trend of average absolute error with number of trees under cross-validation

由图4可以看出，最佳的弱决策树数量为380棵。即在380棵后继续增长弱决策树数量，会出现平均绝对误差负增长情况，模型停止树的增长，同时不再计算平均绝对误差。

在获取了最佳弱决策树的数量后，即确定了一个集成器模型的基本规模。接下来确定树的深度和内部节点再划分最小样本数。使用网格搜索与交叉验证函数GridSearchCV()实现自动调参，只要把参数输进去，就能给出最优化的结果和参数。每一次调参之前都把上一步调出的参数代入模型。设n_estimators为380，max_depth的范围为(3,14)，搜索步长为2，内部节点再划分所需最小样本数min_samples_split范围为(100,801)，搜索步长为200。最佳搜索结果为：max_depth为5,min_samples_split为100。由于决策树深度5是一个比较合理的值，所以可以确定。而min_samples_split和决策树其他的参数存在关联，所以需要继续调参。

使用min_samples_split和叶子节点最少样本数min_samples_leaf一起调参。设置min_samples_split的范围为(600，1 900)，搜索步长为200。min_samples_leaf的范围为(60,101)，搜索步长为10。得到的搜索结果为：min_samples_leaf最佳为60, min_samples_split最佳为800。

再对最大特征数max_features进行网格搜索。将以上调出的参数代入。设搜索范围为(1,8)，步长为2，得到的搜索结果为5。对子采样比例subsample进行调参。设置搜索数组为(0.6,0.7,0.75,0.8,0.85,0.9)，调参最佳结果为0.8。对学习率设置搜索数组为(0.005，0.01，0.05，0.1)，调参最佳结果为0.01。模型参数的最终结果如表4所示。

表4 模型参数调节结果Tab.4 Model parameters adjustment result

使用75%的数据进行模型的构建，使用25%的数据进行模型的验证。模型精度的评价标准选择决定系数R2和平均绝对误差MAE。将土壤阻力，灰度直方图的方差、偏度、峰度，图像灰度共生矩阵的熵和对比度共6个子因素作为输入，构建梯度提升决策回归树模型。使用验证集对模型进行验证，验证集的预测值和环刀取土法得到的表土层土壤容重进行相关性分析，结果如图5所示。

图5 模型验证集的预测值与环刀法测得值的相关性Fig.5 Correlation between predicted value of model verification set and value measured by ring knife method

为了验证在相同输入参数条件下模型的性能和时间复杂度，选择BPNN模型(训练目标为0.001，学习率为0.05，激活函数为sigmoid)和支持向量机回归模型SVR(C=0.5，γ为0.2，kernel为linear)的精度和计算时间加以对比，使用Python中的time模块计算模型的运算时间。BPNN和SVR模型验证集的预测值与环刀法测得值的相关性分别如图6、7所示。

图6 BPNN模型验证集的预测值与环刀法测得值的相关性Fig.6 Correlation between predicted value of BPNN model verification set and value measured by ring knife method

图7 SVR模型验证集的预测值与环刀法测得值的相关性Fig.7 Correlation between predicted value of SVR model verification set and value measured by ring knife method

3个算法在相同的运算环境和输入参数条件下，模型预测精度和运算时间如表5所示。

表5 3种模型预测结果Tab.5 Prediction results of three models

由表5可以看出，GBRT模型的预测结果最优，且运行时间最少，所以选择GBRT模型预测表土层土壤容重具有更好的性能。

GBRT相比于其他回归树模型，变量的相对重要性是所有回归树变量相对重要性的均值，在回归的过程中，可以对特征进行进一步筛选，相关度低的特征可以不用。同时可以很好地处理缺失特征的数据，因为每个节点只依赖一个特征，如果某个特征不存在，这棵树依然可以用来决策，只是少一些路径，而BPNN和SVR没有这个特点。GBRT对于特征空间的异常也有很好的鲁棒性，在低维度情况下，数据的规模对模型的影响也不大，因为对弱决策回归树的要求不高。GBRT在建模过程中数据也不需要进行归一化，特征的作用只是用来分裂节点，叶节点的值与特征值无关，是否归一化并不影响叶节点值和梯度提升的进程。

4 结论

(1)利用同态滤波技术对土壤表面图像进行预处理，提取反映图像纹理细致程度的像元特征参数和反映图像纹理排列结构的区域特征参数来表征表土层土壤表面粗糙度，从土壤阻力和土壤表面粗糙度的信息融合入手，利用灰度关联算法筛选出与表土层土壤容重关联度高的特征参数，使用梯度提升决策回归树模型对表土层土壤容重进行预测，并与BPNN和SVR模型的预测结果进行了对比，验证得到GBRT模型的预测结果具有更高的预测精度和更短的运行时间。

(2)图像处理和梯度提升决策回归树在表土层土壤容重预测上具有良好的应用前景。在低维数据下，GBRT模型可以取得良好的预测效果，避免了传统容重测量的繁琐步骤，为表土层土壤容重的预测研究提供了参考。但是其调参过程较为繁琐，而且在处理高维度数据时会加大计算的复杂度。