CN102901697A - 土壤孔隙度检测方法 - Google Patents

Abstract

本发明公开了一种土壤孔隙度检测方法，涉及农业及计算机图像技术领域，包括：S1：采集土壤样品的显微图像，并对所述显微图像进行预处理，以将所述显微图像转化成灰度图像并去除噪声；S2：生成去噪后灰度图像的灰度共生矩阵；S3：计算灰度共生矩阵的特征值：对比度CON、相关性COR、逆差矩IDM和熵ENT；S4：根据事先设定的拟合指标对所述对比度CON、相关性COR、逆差矩IDM和熵ENT进行拟合得到孔隙度的拟合方程；S5：将待测土壤的灰度共生矩阵的对比度CON、相关性COR、逆差矩IDM和熵ENT代入所述拟合方程得到所述待测土壤的孔隙度。本发明能够更加快速而准确的计算出土壤的孔隙度。

Description

土壤孔隙度检测方法

技术领域

本发明涉及农业及计算机图像技术领域，特别涉及一种土壤孔隙度检测方法。

背景技术

土壤孔隙度及颗粒组成是土壤重要的物理性质参数，是研究土壤最基础的数据资料，尤其是在土壤模型研究和土工试验方面。其中颗粒分布与土壤结构、土壤质地有密切关系，孔隙度直接影响水气二相在土体中的迁移与根系在土壤中的生长。

目前，测量土壤孔隙度的方法主要分为两种，一种是传统的实验法，一种是计算机图形学方法。传统的实验法以烘干法为代表，实验步骤为：首先用环刀取所需土样并称重，然后放入烘箱中于105℃下烘干24小时，取出后称重。计算过程如下：

1、土壤容重的计算：

$\mathrm{rs}=\frac{g}{V(1+w)}$

rs表示土壤容重（克/立方厘米），g表示环刀内湿土重（克），V表示环刀容积（立方厘米），w表示土壤含水率（g/kg）（如含水率为200g/kg应写成0.2）。

2、土壤含水率的计算：

$w=\frac{g}{1+w}$

3、土壤总孔隙度的计算：

$f(\%)=\left(\frac{\mathrm{ds}-\mathrm{rs}}{\mathrm{ds}}\right)\×100$

ds表示土粒密度（一般认为耕地土壤表土的土粒密度为2.65克/立方厘米，土粒密度与土壤比重的数值相等，密度单位为克/立方厘米，比重无单位）。

计算机图形学方法为，首先将图像转化为灰度图像如图1中（a）所示，然后选取阈值将图像转化为二值图像，如图1中（b）所示，阈值以上即为孔隙，最后计算二值图像中孔隙面积与总面积的比值即为孔隙度。

土壤颗粒分布的测量方法主要是吸管法，由筛分和静水沉降组成。通过2mm筛孔的土样经化学和物理方法处理成悬浮液定容后，根据司笃克斯（Stokes）定律及土粒在静水中的沉降规律，大于0.25mm的各级颗粒由一定孔径的筛子筛分，小于0.25mm的粒级颗粒则用吸管从其中吸取一定量的各级颗粒，烘干称量，计算各级颗粒含量的百分数。无相关计算机图形学方法。

以上所述方法缺陷分别为：传统方法速度慢、试验周期长，烘干法需要10小时左右，吸管法则需数天；操作复杂，对实验员操作水平有较高要求；而且实验结果相对误差较大。通过计算机图形学计算土壤孔隙度的方法最大的限制在于，阈值的选取没有公认的标准，不同的阈值选取会得到不同的孔隙度结果，对研究者的经验依赖较强。

发明内容

（一）要解决的技术问题

本发明要解决的技术问题是：如何快速而又较为准确的测量土壤孔隙度。

（二）技术方案

为解决上述技术问题，本发明提供了一种土壤孔隙度检测方法，包括以下步骤：

S1：采集土壤样品的显微图像，并对所述显微图像进行预处理，以将所述显微图像转化成灰度图像并去除噪声；

S2：生成去噪后灰度图像的灰度共生矩阵；

S3：计算灰度共生矩阵的特征值：对比度CON、相关性COR、逆差矩IDM和熵ENT；

S4：根据事先设定的拟合指标对所述对比度CON、相关性COR、逆差矩IDM和熵ENT进行拟合得到孔隙度的拟合方程；

S5：将待测土壤的灰度共生矩阵的对比度CON、相关性COR、逆差矩IDM和熵ENT代入所述拟合方程得到所述待测土壤的孔隙度。

其中，所述步骤S1中对所述显微图像进行预处理包括：灰度转化，阈值滤波去除噪声，灰度直方归一化，其中，灰度转化公式为：

Gray=0.3Rd+0.59Grn+0.11Bl

其中，Gray为转化后灰度值，Rd为红色通道值，Grn为绿色通道值，Bl为蓝色通道值。

其中，所述步骤S2中灰度共生矩阵的生成方法如下：

对于灰度图像中任一个像素点，以该像素点为中心，定义该像素点所在的水平方向为0°，逆时针旋转依次定义为45°、90°和135°，沿该四个角度统计灰度值为i，j，距离为d的像素点对，以得到四个不同的灰度共生矩阵。

其中，所述距离d取8~32范围之内的值。

其中，所述步骤S3具体计算方式如下：

对比度CON： $f_1={\mathrm{\Σ}}_{n=0}^{N_{g-1}}{n}^2\left\{{\mathrm{\Σ}}_{i=1}^{N_g}{\mathrm{\Σ}}_j^{N_g}p(i,j)\right|\left(\right|i-j\left|\right)=0\}$

相关性COR：f₂=(∑_i∑_j(ij)p(i，j)-μ_xμ_y)/σ_xσ_y

逆差矩IDM：f₂=(∑_i∑_jp(i,j))/(1+(i-j)²)

熵ENT：f₄=-∑_i∑_jp(i,j)log(p(i,j))

其中，μ_x，μ_y是p_x，p_y的均值，σ_x，σ_y是p_x，p_y的标准差，p_x，p_y为p(i，j)的边缘概率密度，p(i，j)为灰度值分别为i、j的像素对出现的次数组成的矩阵归一化后得到的灰度共生矩阵，N_g为灰度图像的灰度等级，n取值为0~N_g-1，（ij）表示i点乘j。

其中，所述步骤S4具体包括：

S4.1：从对比度CON、相关性COR、逆差矩IDM和熵ENT中任选其一按多元回归方程进行拟合，并计算出测定系数、校正测定系数和标准差三种拟合指标；

S4.2：加入一种特征值再次拟合，计算拟合指标与前次进行比较，若校正测定系数减小说明新加的特征值与多元回归方程的因变量无关；若校正测定系数增加，说明新加值与因变量有关，改变特征值距离，直到测定系数不再增加；当校正测定系数和测定系数变化不超过5%时选择标准差最小时的一组特征值参与拟合；

S4.3：继续步S4.2直到四种特征值都经过拟合。从而得到孔隙度的拟合方程如下：

F₀=0.00191CON₃₂+159.3COR₁₆-51.055ENT₁₆-3545.13IDM₁₆+267.18

其中，F₀为孔隙度，CON₃₂表示d为32时的灰度共生矩阵的对比度，COR₁₆表示d为16时的灰度共生矩阵的相关性，ENT₁₆表示d为16时的灰度共生矩阵的熵，IDM₁₆表示d为16时的灰度共生矩阵的逆差矩。

其中，步骤S3之后还包括：

根据事先设定的拟合指标对所述对比度CON、相关性COR、逆差矩IDM和熵ENT进行拟合得到土壤颗粒含量的拟合方程；

将待测土壤的灰度共生矩阵的对比度CON、相关性COR、逆差矩IDM和熵ENT代入所述土壤颗粒含量拟合方程得到所述待测土壤的颗粒含量。

其中，所述土壤颗粒含量的拟合方程为：

D₁=2.1232CON₃₂+22265COR₃₂+675.45ENT₁₆+31427IDM₃₂－26041

D₂=-2.4851CON₃₂－26126COR₃₂－748.06ENT₁₆－34964IDM₃₂+30386

D₃=0.36186CON₃₂+3680.7COR₃₂+72.612ENT₁₆+3537.2IDM₃₂－4245，

其中，D₁表示2~0.05mm颗粒含量，D₂表示0.05~0.002mm颗粒含量，D₃表示小于0.002mm颗粒含量。

（三）有益效果

本发明所提出的土壤孔隙度检测方法对显微图像进行上述步骤的处理，相比其他方法可以更加快速而准确的计算出土壤的孔隙度。同时可以大幅缩短相关研究的试验周期，并提高相关研究结果的准确率。

附图说明

图1是现有技术中的计算机图形学方法检测土壤孔隙度的结果图，（a）为灰度图像，（b）为通过阈值处理后的二值图像；

图2是本发明实施例的一种土壤孔隙度检测方法流程图；

图3是图1的步骤S 100中图像预处理过程及结果示意图，其中，（a）为转化后的灰度图像，（b）为灰度图像滤波后的图像，（c）为归一化后的图像；

图4是图1的步骤S200中建立灰度图像共生矩阵的示意图；

图5是不同距离d下的对比度曲线图，横坐标表示距离d，纵坐标无量纲；

图6是不同距离d下的相关性曲线图，横坐标表示距离d，纵坐标无量纲；

图7是不同距离d下的熵曲线图，横坐标表示距离d，纵坐标无量纲；

图8是不同距离d下的逆差矩曲线图，横坐标表示距离d，纵坐标无量纲。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

本发明的土壤孔隙度检测方法流程如图2所示，包括：

步骤S100，采集土壤样品的显微图像，并对显微图像进行预处理，以转化成灰度图像并去除噪声。具体步骤如下：

用环刀取土壤样品后密封，使土壤保持原样，并尽量减少土壤中水分流失。利用显微镜获取显微图像。并进行灰度转化，阈值滤波去除噪声，灰度直方归一化。灰度转化公式为：

Gray=0.3Rd+0.59Grn+0.11Bl    (1)

其中Gray为转化后灰度值，Rd为红色通道值，Grn为绿色通道值，Bl为蓝色通道值，预处理过程中的图像如图3所示，（a）为转化后的灰度图像，（b）为灰度图像滤波后的图像，（c）为归一化后的图像。

步骤S200，生成灰度共生矩阵，包括步骤：

灰度共生矩阵的生成方法为：设Lx={1,2,…,Nx}和Ly={1,2,…,Ny}分别为图像在X、Y方向上的空间域，则Lx×Ly即为图像像素的集合，G∈{1,2,…,Ng}是图像灰度等级，每个像素的灰度等级由方程I:Lx×Ly→G表示。如图4所示，为像素*的与周围像素角度及距离关系示意图，区域1和5为像素*的0°最近邻像素；区域2和6为135°最近邻像素；区域3和7为90°最近邻像素；区域4和8为45°最近邻像素；区域1、3、5和7与*距离为1；区域2、4、6和8与*距离为2。

首先统计图像中两个灰度值分别为i、j，且两者之间角度为θ，距离为d的像素对数。其中，i和j为公式（1）计算得到的灰度值，角度θ的确定如图2中所示有0°，45°，90°，135°四个方向；d为街区距离（CityBlock Distance）等于两像素横坐标的坐标差的绝对值与纵坐标的坐标差的绝对值之和，如图2中像素*和像素1、3、5、7的距离为1，和像素2、4、6、8的距离为2。各方向上灰度值为i、j的像素成对出现的次数定义为：

P(i，j,d,0°)=#{[(k,l),(m,n)]|k-m=0,|l-n|=d，I(k,l)=i,I(m,n)＝j}

P(i，j,d,45°)=#{[(k,l),(m,n)]|k-m=d,l-n=d，I(k,l)=i,I(m,n)＝j}

P(i，j,d,90°)=#{[(k,l),(m,n)]||k-m|=d,l-n=0,I(k,l)=i，I(m,n)＝j}

P(i，j,d,135°)=#{[(k,l),(m,n)]|k-m=-d,l-n=-d，I(k,l)=i,I(m,n)＝j}(2)

式中(k,l)，(m,n)为图像中的像素点，k和l分别为像素(k,l)的横纵坐标，m和n分别为像素(m,n)的横纵坐标，I(k,l)，I(m,n)为像素灰度值；#{}即为集合Lx×Ly中的符合条件的像素对数目。

P(i，j)组成的矩阵就是图像的灰度空间依赖矩阵，将矩阵归一化既可得到图像的灰度共生矩阵。归一化过程如下式所示：

$p(i,j)=P(i,j)/{\mathrm{\Σ}}_{i=1}^{N_g}{\mathrm{\Σ}}_{j=1}^{N_g}P(i,j)---\left(3\right)$

生成灰度共生矩阵的参数选取方式如下：

生成灰度共生矩阵共需确定三个参数：灰度等级N_g，角度θ，距离d。

在灰度共生矩阵算法应用早期，由于计算能力的限制，研究者通常采用减少灰度等级的方法减小计算量，但是这样会损失一定的灰度空间依赖信息。现在计算能力已经不是研究需要考虑的重点，所以研究中灰度等级使用原图像的灰度等级256不作调整。

由公式2可知，沿4个角度统计灰度值为i，j，距离为d的像素点对时，可以得到4个不同的灰度共生矩阵，根据矩阵计算出的特征值也因此不同。对于各向同性的纹理，角度对特征值的影响不大，但对各向异性的纹理，角度的影响则相对明显。例如相关性的数值就被用来分析纹理的走向，根据各角度得到的不同相关性值，数值最大的角度即为纹理走向的角度。为了消除角度对结果的影响，实验取各角度特征值的算数平均值作为函数拟合的输入。

公式（2）中最后一个参数是像素点对之间的距离d，大部分研究者在应用是一般取1，有研究者指出，当像素距离足够大时，特征值趋于一致。但是特征值分辨力最大点不一定在稳定值，所以实验中对每个特征值分别进行讨论，选取各样本分辨能力最好的距离d。一般选取8~32范围之内的值。

步骤S300，特征值的计算

选取不相关且分辨力较强的4个特征值：对比度（Contrast），相关性（Correlation），逆差矩（Inverse Difference Moment），熵（Entropy）。计算公式分别如下所示：

对比度（CON）： $f_1={\mathrm{\Σ}}_{n=0}^{N_{g-1}}{n}^2\left\{{\mathrm{\Σ}}_{i=1}^{N_g}{\mathrm{\Σ}}_j^{N_g}p(i,j)\right|\left(\right|i-j\left|\right)=0\}$

相关性（COR）：f₂=(∑_i∑_j(ij)p(i，j)-μ_xμ_y)/σ_xσ_y

逆差矩（IDM）：f₃=(∑_i∑_jp(i,j))/(1+(i-j)²)

熵（ENT）：f₄=-∑_i∑_jp(i,j)log(p(i,j))（4）

式中μ_x，μ_y是p_x，p_y的均值，σ_x，σ_y是p_x，p_y的标准差，p_x，p_y为p(i，j)的边缘概率密度，p(i，j)为灰度值分别为i、j的像素对出现的次数组成的灰度共生矩阵，n取值为0~N_g-1，（ij）就是i点乘j。

图5~图8所示为各特征值与距离的关系，横坐标为特征值距离d的值，纵坐标依次为对比度、相关性、逆差矩和熵等纹理特征值。

步骤S400，多元线性回归

多元回归方程的一般形式如式5所示：

Y=β₀+β₁X₁+β₂X₂+…+β_kX_k+ε(5)

Y表示因变量；X₁，X₂，…，X_k表示自变量，k是变量个数；β₁，β₂，…，β_k为回归系数，是未知数，通过观测数据来估计；ε为随机误差。回归系数使用最小二乘法求解。本实验中Y可表示各样品的孔隙度和各颗粒含量，本实施例中k＝4，X₁至X₄分别为4种纹理特征参数：对比度、相关性、逆差矩和熵。

介绍拟合过程前，首先引入三种评判方程拟合优度的指标，测定系数R²(COD)，校正测定系数Adj.R²，标准差SD。测定系数表示方程与实际值的密和程度，值越接近1密和度越高；校正测定系数则考虑了是否引入无统计意义的自变量，若引入则下降；标准差指的是实际值与方程之间误差的标准差。

具体拟合过程如下：

1）在步骤S300选出的特征值中任选其一进行拟合，并计算出三种拟合指标；

2）加入一种特征值再次拟合，计算拟合指标与前次进行比较，若Adj.R²减小说明新加值与因变量无关，即与土壤孔隙度无关；若Adj.R²增加，说明新加值与因变量有关，改变特征值距离，即上述的街区距离d，直到R²(COD)不再增加；当Adj.R²，R²(COD)变化均不大时，如不超过5%，则选择SD最小时的一组特征值参与拟合；

3）继续步骤2直到4种特征值都经过拟合。

所得方程如下所示：

F₀=0.00191CON₃₂+159.3COR₁₆-51.055ENT₁₆-3545.13IDM₁₆+267.18

其中F₀为孔隙度，CON₃₂表示以d为32的灰度共生矩阵（GLCM）求出的对比度，COR₁₆表示d为16时的灰度共生矩阵的相关性，ENT₁₆表示d为16时的灰度共生矩阵的熵，IDM₁₆表示d为16时的灰度共生矩阵的逆差矩。

2~0.05mm颗粒含量(%)D₁，0.05~0.002mm颗粒含量(%)D₂，小于0.002mm颗粒含量(%)D₃的拟合公式分别如下：

D₁=2.1232CON₃₂+22265COR₃₂+675.45ENT₁₆+31427IDM₃₂－26041

D₂=-2.4851CON₃₂－26126COR₃₂－748.06ENT₁₆－34964IDM₃₂+30386

D₃=0.36186CON₃₂+3680.7COR₃₂+72.612ENT₁₆+3537.2IDM₃₂－4245

步骤S500，孔隙度的计算

根据第五部分拟合出的方程，及该部分确定的方程中的矩阵参数建立待检测样品的灰度共生矩阵，计算所需特征值并代入步骤S400得到的回归方程F₀求解待测土壤的孔隙度，将特征值代入D₁~D₃求解得到待测土壤的颗粒含量。

以上实施方式仅用于说明本发明，而并非对本发明的限制，有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型，因此所有等同的技术方案也属于本发明的范畴，本发明的专利保护范围应由权利要求限定。

Claims (8)

1.一种土壤孔隙度检测方法，其特征在于，包括以下步骤：

S1：采集土壤样品的显微图像，并对所述显微图像进行预处理，以将所述显微图像转化成灰度图像并去除噪声；

S2：生成去噪后灰度图像的灰度共生矩阵；

S3：计算灰度共生矩阵的特征值：对比度CON、相关性COR、逆差矩IDM和熵ENT；

S4：根据事先设定的拟合指标对所述对比度CON、相关性COR、逆差矩IDM和熵ENT进行拟合得到孔隙度的拟合方程；

S5：将待测土壤的灰度共生矩阵的对比度CON、相关性COR、逆差矩IDM和熵ENT代入所述拟合方程得到所述待测土壤的孔隙度。

2.如权利要求1所述的土壤孔隙度检测方法，其特征在于，所述步骤S1中对所述显微图像进行预处理包括：灰度转化，阈值滤波去除噪声，灰度直方归一化，其中，灰度转化公式为：

Gray=0.3Rd+0.59Grn+0.11Bl

其中，Gray为转化后灰度值，Rd为红色通道值，Grn为绿色通道值，Bl为蓝色通道值。

3.如权利要求1所述的土壤孔隙度检测方法，其特征在于，所述步骤S2中灰度共生矩阵的生成方法如下：

对于灰度图像中任一个像素点，以该像素点为中心，定义该像素点所在的水平方向为0°，逆时针旋转依次定义为45°、90°和135°，沿该四个角度统计灰度值为i，j，距离为d的像素点对，以得到四个不同的灰度共生矩阵。

4.如权利要求3所述的土壤孔隙度检测方法，其特征在于，所述距离d取8～32范围之内的值。

5.如权利要求1所述的土壤孔隙度检测方法，其特征在于，所述步骤S3具体计算方式如下：

对比度CON： $f_1={\mathrm{\Σ}}_{n=0}^{N_{g-1}}{n}^2\left\{{\mathrm{\Σ}}_{i=1}^{N_g}{\mathrm{\Σ}}_j^{N_g}p(i,j)\right|\left(\right|i-j|=0)\}$

相关性COR：f₂=(∑_i∑_j(ij)p(i，j)μ_x-μ_y)/σ_xσ_y

逆差矩IDM：f₃=(∑_i∑_jp(i,j))/(1+(i-j)²)

熵ENT：f₄=-∑_i∑_jp(i,j)log(p(i,j))

其中，μ_x，μ_y是p_x，p_y的均值，σ_x，σ_y是p_x，p_y的标准差，p_x，p_y为p(i,j)的边缘概率密度，p(i,j)为灰度值分别为i、j的像素对出现的次数组成的矩阵归一化后得到的灰度共生矩阵，N_g为灰度图像的灰度等级，n取值为0～N_g-1，（ij）表示i点乘j。

6.如权利要求1所述的土壤孔隙度检测方法，其特征在于，所述步骤S4具体包括：

S4.1：从对比度CON、相关性COR、逆差矩IDM和熵ENT中任选其一按多元回归方程进行拟合，并计算出测定系数、校正测定系数和标准差三种拟合指标；

S4.2：加入一种特征值再次拟合，计算拟合指标与前次进行比较，若校正测定系数减小说明新加的特征值与多元回归方程的因变量无关；若校正测定系数增加，说明新加值与因变量有关，改变特征值距离，直到测定系数不再增加；当校正测定系数和测定系数变化不超过5%时选择标准差最小时的一组特征值参与拟合；

S4.3：继续步S4.2直到四种特征值都经过拟合，从而得到孔隙度的拟合方程如下：

F₀=0.00191CON₃₂+159.3COR₁₆-51.055ENT₁₆-3545.13IDM₁₆+267.18其中，F₀为孔隙度，CON₃₂表示d为32时的灰度共生矩阵的对比度，COR₁₆表示d为16时的灰度共生矩阵的相关性，ENT₁₆表示d为16时的灰度共生矩阵的熵，IDM₁₆表示d为16时的灰度共生矩阵的逆差矩。

7.如权利要求1～6中任一项所述的土壤孔隙度检测方法，其特征在于，步骤S3之后还包括：

根据事先设定的拟合指标对所述对比度CON、相关性COR、逆差矩IDM和熵ENT进行拟合得到土壤颗粒含量的拟合方程；

将待测土壤的灰度共生矩阵的对比度CON、相关性COR、逆差矩IDM和熵ENT代入所述土壤颗粒含量拟合方程得到所述待测土壤的颗粒含量。

8.如权利要求7所述的土壤孔隙度检测方法，其特征在于，所述土壤颗粒含量的拟合方程为：

D₁=2.1232CON₃₂+22265COR₃₂+675.45ENT₁₆+31427IDM₃₂－26041

D₂=-2.4851CON₃₂－26126COR₃₂－748.06ENT₁₆－34964IDM₃₂+30386

D₃=0.36186CON₃₂+3680.7COR₃₂+72.612ENT₁₆+3537.2IDM₃₂－4245，

其中，D₁表示2～0.05mm颗粒含量，D₂表示0.05～0.002mm颗粒含量，D₃表示小于0.002mm颗粒含量。

Images