MATLAB数字图像处理在根系固土中的应用研究

阎梦晴，于司杭，王培亮，王逸飞，李明宝

(东北林业大学土木工程学院，黑龙江哈尔滨150040)

我国是世界上水土流失最严重的国家之一［1］。水土流失会破坏土壤的肥力，造成土壤干层，淤积水库、阻塞河道、抬高河床，恶化生态环境。虽然自新中国成立以来，我国的水土流失治理已经逐渐成为国建的生态建设重点工程，但是由于我国国土面积大，自然地理条件特殊，以及长期以来对水土资源的过度利用，我国水土流失仍然面积大、分布广，防治任务艰巨［2］。防止水土流失，采用科学、环保的植物固土措施，引入数字图像处理技术研究植物根系固土机制，能够为我国农业水土保持增添新的技术手段，促进我国生态文明建设的发展。

水土流失的防治措施有工程措施和植物措施，其中植物措施不仅可以防治水土流失，还可改善局部的生态环境，节约经济成本，具有生态可持续性。近年来，由于生态环境建设的需要及各学科科学技术的发展，植物根系固定土体机理的研究成为根系研究的焦点［3］。根系与土壤形成复合结构体，起到固定植物地上部分以及固定土壤、防止水土流失的作用。对植物根系固定土体的机理进行综合定量研究，对于搭建优良的人工水土保持生态系统不仅具有广阔的应用前景，而且具有重要的学术意义［1］。对于根系固土力学机理的研究，已经得到越来越多学者的重视，植物根系防治水土流失已经成为水土保持领域的一个研究热点。

近年来，国内外科研专家运用各种科学技术方法对根系固定土体的力学机理进行了研究，测定植物根系本身的抗拉力［4－5］，采用原位剪切试验得到的有根土和无根土的抗挤强度值作为衡量根系固土性质的指标［6］。随着科学技术手段的不断进步，土壤中根系的几何形态和分布情况能够通过各种试验方法得到，在研究植物根系在土壤中的分布以及根系几何形态的提取与仿真方面，也不可避免地需要对根土断层图像、土壤裂缝图像进行数字图像处理，为科研人员提供根系的形态参数数据。

运用数字图象处理技术，对根系固土抗裂性试验中的试验图像进行特征提取，将裂缝的图像特征(如裂缝宽度、深度等)用以衡量根系固土的抗裂有效性，补充完善了宏观的试验指标，通过这些指标可以直观地反映根系固土的程度。通过对图像特征的提取，建立更加全面的参数，这些参数可作为研究抗裂性及其影响因素之间的桥梁。运用数字图像处理，能够使根系固土抗裂性这一性质转化为更加具体且全面的数学指标，通过这些量化的指标可建立抗裂性与其影响因素之间的关系，探究它们的影响机理，这将为研究植物根系固土防治水土流失奠定很好的理论基础。数字图像处理技术的引入，将成为水土保持研究领域的一个重要环节，将为分析根系固土机制提供有效的技术手段。

1 基于MATLAB的数字图像处理

数字图像处理，把一幅图像定义为一个二维函数f(x，y)，其中x、y是空间坐标，f是任意坐标(x，y)处的幅度(称为亮度或灰度)。MATLAB(矩阵实验室)是处理矩阵和线性代数的强大的数值计算软件，将图像作为矩阵来处理，选用MATLAB软件能有效地提高图像处理的效率，准确地提取出数字图像的特征。将坐标值(x，y)数字化称为采样，将幅值8数字化称为量化，采样和量化后得到一个实数矩阵。假设，对图像f(x，y)采样后得到M×N的图像，坐标值为离散量。将基于MATLAB的数字图像处理技术应用于根系固土抗裂性试验，旨在得到土壤裂缝的图像特征。

1.1 边缘检测算法 从输入的图像中提取属性，进行图像分割。分割算法的基本原理是基于灰度的突变来分割图像，基于检测灰度值的不连续来检测边缘。灰度的不连续性的检测采用一阶导数和二阶导数。图像处理中的二阶导数由于对噪声过于敏感而很少直接用于边缘检测，仅在与其他边缘检测技术结合时使用，在此不作讨论。

图像二维函数f(x，y)的梯度为

向量∇f的幅值为

梯度的幅值通常简称为梯度，由式(1)可知梯度具有微分的性质，当灰度恒定时梯度值为零，当灰度变化时梯度值相应变化，由此可通过计算梯度值来检测灰度的变化程度，从而检测边缘。确定图像中灰度快速变化的位置，基本准则是寻求灰度的一阶导数的幅值大于某个指定阈值的位置。在MATLAB中，采用edge函数进行边缘检测。

在MATLAB中，根据土壤裂缝的图像以及边缘检测器的功能，选用 Sobel、Prewitt、Roberts、LoG、Canny 这 5 种边缘检测器，对比它们边缘检测的效果，以确定最适用于土壤裂缝图像做分割的边缘检测器。

1.2 算法比较

1.2.1 Sobel边缘检测器。Sobel边缘检测器使用图1(b)中的Sobel近似导数发现边缘，即图1(a)所示的3×3邻域的行、列间的离散差来计算梯度，每行每列的中心像素以2加权提供平滑效果［7］。

图1 边缘检测器模板及其实现的一阶导数

式中，z为灰度。在(x，y)处，当∇f≥T(T为指定阈值)，则(x，y)处的像素为边缘像素。

1.2.2 Prewitt边缘检测器。Prewitt边缘检测器使用图1(c)中的近似导数发现边缘。与Sobel检测器相比，Prewitt边缘检测器的计算更简单，但是检测结果中的噪声较大。

1.2.3 Roberts边缘检测器。Roberts边缘检测器使用图1(d)的近似导数发现边缘。Roberts边缘检测器是最古老且最简单的数字图像处理边缘检测器之一，优点在其速度快，但功能不全面。

1.2.4 LoG 检测器。LoG(Laplacian of a Gaussian)检测器使用了高斯拉普拉斯算子与图像卷积(滤波)。对于高斯函数

高斯拉普拉斯算子为

使用∇2G(x，y)卷积图像，使得图像平滑，降低了噪音，计算拉普拉斯算子得到双边缘图像，查找双边缘之间的零交叉即可定位边缘。

1.2.5 Canny边缘检测器。与LoG检测器相比，Canny边缘检测器指定了高斯滤波器中的标准差σ，以此来平滑图像，减少噪声。在(x，y)处计算局部梯度［+］1/2和边缘方向arctan(gx/gy)，通过寻找f(x，y)的梯度局部最大发现边缘。Canny边缘检测器使用了两个阈值来检测强边缘和弱边缘，此方法更有可能检测到真正的弱边缘。

2 结果与分析

为得到土壤裂缝的形态特征参数，要选用适宜的边缘检测器以及阈值处理方法。通过5种不同边缘检测器的处理结果图与原裂缝图的对比，选择出精度最高的边缘检测器。在此基础上，采用两种不同的阈值处理方法，对比所得图像的精确度，确定最适宜的阈值处理方法。

2.1 边缘检测器的选用 土壤的裂缝图像来自于牛国权［8］所做的根系固土抗裂有效性试验，见图2(a)。

图2为5种不同边缘检测器处理土壤裂缝图像后的结果图。由图2(b)、(c)可知，Prewitt边缘检测器和Sobel边缘检测器的效果均不佳，二者显然无法检测到图像中的弱边缘，无法有效分割出土壤裂缝的图像。Roberts边缘检测器原理如图1所示，与Prewitt和Sobel边缘检测器相比，Roberts边缘检测器模板是非对称的，虽然简单，但用于检测土壤裂缝的图像特征时，如图2(d)所示，其效果不好。由图2(e)可知，LoG检测器结果比前3种检测器的效果好，但不及Canny边缘检测器。由图2(f)可知，Canny边缘检测器功能强大，能检测到细微的图像特征，但是在检测土壤裂缝时，还需调整参数尽可能地减少不相关的细节。

图2 不同边缘检测器的处理结果

2.2 图像阈值处理 做分割后得到的裂缝图像是初步处理结果。在初步处理图像时，选择的是各个检测器自动选择的阈值，得到的结果不精确。需根据裂缝图像的图形特征，选择最优阈值。

阈值处理后的图像为二值图像，它定义为

式中，a像素对应物体，b像素对应背景，按照惯例，a=1(白)，b=0(黑)。

2.2.1 使用Otsu方法进行最佳全局阈值处理。将图像的直方图的成分表示为

其中，n为图像像素总和，nq为具有灰度级q的像素数量，L为图像中可能的灰度级总数(灰度级为整数)。

假设初始阈值为 k，C1 是灰度级为［0，1，2，…，k］的一组像素，C2是灰度级为［k+1，k+2，…，L －1］的一组像素，则最大类间方差为

其中，P1(k)为集合C1发生的概率，P2(k)为集合C2发生的概率，mG为全局均值(平均灰度)。

令P2(k)=1－P1(k)，则式(8)可写为

通过类间方差最大化思想即可确定最佳阈值，在MATLAB中，采用工具箱函数graythresh()可计算Otsu阈值T=0.431 4。据此阈值再采用Canny边缘检测器即可得到如图3所示图像。

2.2.2 使用基于梯度的边缘信息改进全局阈值处理。在MATLAB中，将土壤裂缝图像作为输入图像，采用imhist()函数，用plot方法绘制输入图像的直方图。如图4所示，该图像直方图是单峰的，通常，当物体远小于背景时，物体对直方图的贡献可忽略不计，为改善这种情况，往往使用边缘信息以得到梯度图像。

图3 使用Otsu方法进行阈值处理后的结果

图4 输入图像的直方图

自定义函数percentile2i()，用以设定在边缘图像中将较少的像素用于计算阈值，即使用高百分位(99.9)来估计梯度的阈值，因为梯度图像中较大的值是出现在物体与背景的边界处。处理结果如图5所示，阈值T=0.368 6。

图5 基于梯度的边缘信息改进全局阈值处理结果

比较图3和图5，图像阈值处理的结果显然前者较好，这说明在处理土壤裂缝这类图片时，使用Otsu方法进行阈值处理能得到最优阈值。结合Canny边缘检测器和Otsu阈值处理方法，最后得到的裂缝图像如图6所示。

图6 最终结果图

土壤裂缝图像不同于其他图像，为了获取裂缝的形态参数，通过对比研究了5种边缘检测器以及2种阈值处理方法，确定了最优边缘检测器和阈值处理方法，即结合Canny边缘检测器和Otsu阈值处理的方法，能够得到最精确的裂缝图像。图6所示的图像能够精确地反映原图中土壤的裂缝形态。将二维图像数字化，可得到裂缝边缘的二值化图像，在MATLAB中以矩阵形式存储了图像像素，再后期运用裂缝形态特征参数进行根系固土机制研究时，即可直接调用。

3 结论

该研究以根系固土试验中土壤裂缝图像为例，运用MATLAB软件，对图像进行数字图像处理，对比Sobel、Prewitt、Roberts、LoG、Canny这5 种边缘检测器的效果，选择最适合用于检测土壤裂缝的Canny边缘检测器。比较基于梯度的边缘信息改进全局阈值处理方法和Otsu方法进行最佳全局阈值处理的效果，确定使用Otsu方法进行阈值处理能得到最优阈值。结合Canny边缘检测器和Otsu阈值处理方法，是最适用于土壤裂缝图像的数字图像处理方法。

为防止水土流失，对植物根系固定土体机制的研究中，引入基于MATLAB的数字图像处理方法，能够得到宏观图像中如根系、土壤裂缝等图像的形态特征参数。处理后的根系、裂缝等边缘像素的特征以及根系图像像素尺寸与实际尺寸存在线性关系，根据各个参数的几何特性可分析出根系的长度、表面积、平均直径、体积以及根系间的夹角，或是裂缝的长度、宽度、面积等。通过该研究结果，可以为关于根系固土相关领域的研究人员提供可靠分析手段以及根系的形态参数数据。

［1］刘震.水土保持60年:成就·经验·发展对策［J］.中国水土保持科学，2009(4):1 －6.

［2］田卫堂，胡维银，李军，等.我国水土流失现状和防治对策分析［J］.水土保持研究，2008，15(4):204 －209.

［3］陈丽华，余新晓，宋维峰，等.林木根系固土力学机制［M］.北京:科学出版社，2008:87－115.

［4］张瑜.微根窗根系的图像处理方法研究［D］.哈尔滨:东北林业大学，2012.

［5］史敏华，王棣，李任敏.石灰岩区主要水保灌木根系分布特征与根抗拉力研究初报［J］.山西林业科技，1994(1):17－19.

［6］赵丽兵，张宝贵，苏志珠.草本植物根系增强土壤抗剪切强度的量化研究［J］.中国生态农业学报，2008，16(3):718 －722.

［7］GONZALEZ R C，WOODS R E.Digital Image Processing［M］.3rd ed.Prentice Hall，2007.

［8］牛国权.基于力学特性的根系固土抗裂性有效性研究［D］.呼和浩特:内蒙古农业大学，2010.