CN106442526A - 一种基于matlab的活性污泥絮体分析方法 - Google Patents

Abstract

一种基于MATLAB的活性污泥絮体分析方法，包括以下步骤：(一)采集图像：(二)图像处理：(1)读入采集图像；(2)将RGB图像转化为灰度图像；(3)去除噪声处理；(4)图像二值化；(5)开运算和闭运算，光滑边缘处理；(6)确定目标菌胶团并进行标记；(7)填充；(8)目标区域灰度图；(9)计算目标区域灰度图的灰度值之和；(10)计算目标区域面积；(11)计算目标区域密实度；(12)获取目标区域的边界；(三)提取特征信息：(1)计算菌胶团尺寸指标；(2)计算菌胶团形态指标；(3)计算丝状菌指标。该方法以于MATLAB分析软件为平台，应用图像分析技术进行活性污泥絮体的定量分析，保证了分析结果的客观性、可靠性与高效性。

Description

一种基于MATLAB的活性污泥絮体分析方法

技术领域

本发明涉及一种利用图像对活性污泥絮体的密实程度及絮体边缘规则度进行分析的方法，属于污水生物处理过程中活性污泥絮体的分析技术领域。

背景技术

传统针对活性污泥絮体的分析，主要通过在显微镜下人工观测计数的方法研究表征絮体形态特征的絮体面积、圆度、菌丝长度等直观指标，但是无法分析表征絮体内细菌的生长状态、絮体密实度以及絮体边界不规则度，且整个分析过程枯燥、耗时，分析结果受主观影响较大。

发明内容

本发明针对传统活性污泥絮体分析技术存在的不足，提供一种效率高，保证分析结果客观、可靠的基于MATLAB的活性污泥絮体图像分析方法。

本发明的对活性污泥絮体图像分析方法，包括以下步骤：

(一)采集图像：

取新鲜活性污泥制成观测玻片，在显微镜下观察，以获取完整活性污泥菌胶团为标准，拍照并采集图像；

(二)图像处理：

(1)读入采集图像；

(2)将RGB图像转化为灰度图像；

(3)去除噪声处理；

采用中值滤波法进行去除噪声干扰，其不仅可以去除孤点噪声，而且可以保持图像的边缘特性，可避免图像产生显著的模糊以及对后期周长计算的影响。

(4)图像二值化；为了确定目标菌胶团，对图像进行二值化；

(5)开运算和闭运算，光滑边缘处理；

污水生物处理系统活性污泥菌胶团内往往存在一定量的丝状菌，尤其是冬季大量丝状菌伸出菌胶团，影响对菌胶团周长的计算。开运算和闭运算取一定范围的圆盘结构元素，通过调用imopen()和imclose()函数消除不能完全包含在该圆盘内的目标对象，起到了去除碎丝，光滑边缘的作用。

(6)确定目标菌胶团并进行标记；

选取最大联通区域作为最大的菌胶团，利用八邻域标记算法标记目标菌胶团；

(7)填充；

填充二值图像前景区域内出现的黑色孔洞，以降低对计算结果准确性的影响；

(8)目标区域灰度图；

以“密实度”来表示污泥采集图像中菌胶团密实程度，“密实度”为目标区域内灰度值 之和与目标区域面积的比值，图像灰度值的范围为0～255，对应图像中的颜色为从黑到白，越密实的菌胶团图像颜色越深，灰度越小；越松散的菌胶团图像颜色越浅，灰度越大。

(9)计算目标区域灰度图的灰度值之和；

灰度值之和(Sum gray)表征了密实度，将目标区域灰度图的灰度值相加，黑色为0，白色为255，灰度值之和越大表明密实度越小，灰度值之和越小表明密实度越大；

(10)计算目标区域面积；

(11)计算目标区域密实度；

在获得目标区域的面积和灰度值总和的前提下，根据公式：密实度＝灰度值之和/菌胶团面积，计算目标区域密实度；

(12)获取目标区域的边界；

为了计算目标区域的周长，进行圆度计算，通过边缘检测以获取目标区域的轮廓边界；

(三)提取特征信息：

(1)计算菌胶团尺寸指标：

①菌胶团面积(Area)；

②菌胶团周长(Perimeter)

③等效圆直径

④菌胶团长度(length)；

⑤菌胶团宽度(width)；

(2)计算菌胶团形态指标：

①圆度Ro(表征图形与圆之间的偏差，有效描述不规则边界，圆的Ro为1)：

②三维纵横比AR(表征图像的延伸，延伸越大，AR越大，圆的AR为1)：

③密实度Avergray(表征菌胶团的密室程度)：

④分形维数FD(fractal dimension，表征菌胶团边界的粗糙程度)。

分形维数采用盒子维数(即容量维数(f))分形方法，取边长为ε的小盒子对分形曲线进行覆盖，计数整个覆盖范围内非空实体盒子数量，其数量记为N(ε)，然后缩小盒子的尺寸ε，所计数非空实体盒子的数量N(ε)随之增大，当ε趋于0时，得到分形维数(FD)。

(3)计算丝状菌指标：菌丝长度FL(fialment length)。

本发明以于MATLAB分析软件为平台，应用图像分析技术进行活性污泥絮体的定量分析， 分别采用密实度Avergray以及分形维数FD来描述活性污泥絮体形态特征，其中活性污泥絮体的密实度表征了活性污泥絮体内细菌的生长状况与絮体结构松散密实状态；絮体边界分形维数则综合表征了絮体边界的不规则度。在提高工作效率的同时克服了人工分析的主观性，从图像中提取定量数据，保证了分析结果的客观性、可靠性与高效性。

附图说明

图1是两种种污泥絮体镜检图像的灰度图像。图1(a)为第一种；图1(b)为第二种。

图2是两种污泥絮体镜检图像的二值图像。图2(a)为第一种；图2(b)为第二种。

图3是两种污泥絮体镜检图像的二值图像光滑边缘处理后的图像。图3(a)为第一种；图3(b)为第二种。

图4是两种光滑边缘处理后的取反图像。图4(a)为第一种；图4(b)为第二种。

图5是两种取反图像的标记图像。图5(a)为第一种；图5(b)为第二种。

图6是两种填充后的图像。图6(a)为第一种；图6(b)为第二种。

图7是两种目标区域灰度图。图7(a)为第一种；图7(b)为第二种。

图8是目标区域面积计算示意图。

图9是两种目标区域轮廓图。图9(a)为第一种；图9(b)为第二种。

图10是最小外接矩形图。图10(a)为第一种，图10(b)为第二种。

图11是两图像指标信息。图11(a)为第一种；图11(b)为第二种。

图12是应用fractalfox软件得到的两种分形维数值fractal，图12(a)为第一种，图12(b)为第二种。

图13是菌丝长度计算图。图13(a)为第一种，图13(b)为第二种。

具体实施方式

本发明的活性污泥絮体分析方法，以MATLAB分析软件为平台。以下以a(第一种)和吧(第二种)两种污泥絮体图像为实例对本发明的分析方法详细说明。

一.采集图像

取a和b两种新鲜活性污泥制成观测玻片，在显微镜下观察，以获取完整活性污泥菌胶团为标准，选取10个视野观察、拍照并采集图像。

二.图像处理

(1)读入采集图像

进入MATLAB软件，调用imread()函数，读入采集图像。

I＝imread('F:\MATLAB图像\稀释染色2001.GIF')；％读入图像

figure,imshow(I)；％显示图像。

(2)将RGB图像转化为灰度图像

“灰度”可视为图像色彩亮度的深浅。图像灰度级越高，意味着图像可以表现更强的色彩层次。如果将黑-灰-白连续变化的灰度值量化为256个灰度级，灰度值的范围为0～255，对应图像中的颜色为从黑到白。

数码相机拍摄JPG文件格式的RGB图像，将其转化为灰度图像，可大幅度减少图像的数据量，便于图像处理。采用MATLAB软件，调用rgb2hsv函数，转化RGB色彩空间为HSV色彩空 间，并提取HSV空间中的色彩分量(即H分量)，得到的图像即为灰度图像。

调用rgb2hsv()函数，对图像进行灰度化，并保存灰度图。

I1＝rgb2gray(I)；％灰度化

figure,imshow(I1)；

imwrite(I1,'F:\MATLAB图像\灰度图.GIF')；％保存灰度图。

得到的A和B两种新鲜活性污泥的灰度图分别如图1(a)和图1(b)所示。

(3)去除噪声处理

图像采集过程中可能出现椒盐噪声，为避免图像生成和传输过程中各种噪声干扰造成图像质量降低，可采用一种有效滤除脉冲干扰噪声的方法—中值滤波法进行去除噪声干扰，其不仅可以去除孤点噪声，而且可以保持图像的边缘特性，可避免图像产生显著的模糊以及对后期周长计算的影响。调用图像处理工具箱中的medfilt2函数实现。

调用medfilt2()函数，对图像进行中值滤波去噪。

h＝medfilt2(I1,[4,4])；％中值滤波。

(4)图像二值化

为了确定目标菌胶团，对图像进行二值化，即将256个亮度等级的灰度图像通过适当的阈值重置为0或255[大于阈值置为255,小于阈值置为0]，获得可以反映图像整体和局部特征的二值化图像，整个图像呈现出明显的黑白效果。如此有利于在对图像做进一步处理时，图像的集合性质只与像素值为0或255的点的位置有关，不再涉及像素的多级值，使得处理变得简单，且数据的处理和压缩量小。

通过分析图像的灰度直方图，确定分割图像所需最佳阈值，利用该阈值进行二值化处理，调用MATLAB软件中的graythresh()函数确定最佳阈值，并调用im2bw()函数将灰度图像转换为二值图像。

调用graythresh()函数、im2bw()函数

bw＝im2bw(h,graythresh(h))；％二值化，graythresh(h)

figure,imshow(bw)；

imwrite(bw,'F:\MATLAB图像\二值图.GIF')；％保存二值图。

图2(a)给出了第一种污泥絮体镜检图像的二值图像。图2(b)给出了第二种污泥絮体镜检图像的二值图像。

(5)开运算和闭运算-光滑边缘处理

污水生物处理系统活性污泥菌胶团内往往存在一定量的丝状菌，尤其是冬季大量丝状菌伸出菌胶团，影响对菌胶团周长的计算。开运算和闭运算取一定范围的圆盘结构元素，通过调用imopen()和imclose()函数消除不能完全包含在该圆盘内的目标对象，起到了去除碎丝，光滑边缘的作用。

调用imopen()和imclose()函数

bw2＝imopen(bw1,strel('disk',10))；％开运算，圆盘直径为10

bw2＝imclose(bw1,strel('disk',10))；％闭运算，圆盘直径为10，去掉杂丝

figure,imshow(bw2)；

imwrite(bw2,'F:\MATLAB图像\开闭运算.GIF')；％此时为白底黑图。

图3(a)给出了第一种污泥絮体二值图像光滑边缘处理后的图像。图3(b)给出了第二种污泥絮体二值图像光滑边缘处理后的图像。

(6)确定目标菌胶团并进行标记

图像中可能会存在多个分散菌胶团，为了计算的科学性与合理性，需选取最大的菌胶团即最大连通区域。为了标记目标菌胶团，利用八邻域标记算法，bwlabel()函数标记联通区域后，选取最大联通区域，舍弃其他区域。得到目标菌胶团的二值图。

调用bwlabel()函数

bw3＝imcomplement(bw2)；

imwrite(bw3,'F:\MATLAB图像\反黑白.GIF')；％为了标记，取反为黑底。

L＝bwlabel(bw3)；％标记图像联通区域，默认为8邻域

stats＝regionprops(L)；

Ar＝cat(1,stats.Area)；

ind＝find(Ar＝＝max(Ar))；％仅保留取面积最大联通区域

J＝ismember(L,ind)；

figure,imshow(J)；

imwrite(bw3,'F:\MATLAB图像\面积最大区域.GIF')。

图5(a)给出了第一种标记图像；图5(b)给出了第二种标记图像。

(7)填充

为填充二值图像前景区域内出现的黑色孔洞，以降低对计算结果准确性的影响，利用Matlab图像处理工具箱中的imfill()函数进行区域填充操作。

调用imfill()函数

bw1＝imfill(bw,'holes')；％填充

figure,imshow(bw1)；

图6(a)给出了第一种填充后的图像；图6(b)给出了第二种填充后的图像。

(8)目标区域灰度图

本发明提出“密实度”算法来表示污泥采集图像中菌胶团密实程度，其为目标区域内灰度值之和与目标区域面积的比值，即Avergray＝SumGray/Area。图像灰度值的范围为0～255，对应图像中的颜色为从黑到白。越密实的菌胶团图像颜色越深，灰度越小；越松散的菌胶团图像颜色越浅，灰度越大。

Avergray＝Sum Gray/Area；

X＝im2double(h)％为相减改变图像数据类型

Y＝im2double(bw1)％为相减改变图像数据类型

T＝imconplement(Y)；％注意此处再对目标区域取反为白底；

Z＝zeros(size(X))；％新建一个Z与X,Y同大小；

Z＝X-T％将灰度图与目标区域图相减，得到目标区域的灰度图；

figure,imshow(Z)；

imwrite(Z,'F:\MATLAB图像\目标区域灰度.GIF')；％得到目标区域灰度图。

图7(a)给出了第一种目标区域灰度图；图7(b)给出了第二种目标区域灰度图。

(9)计算目标区域灰度图的灰度值之和

灰度值之和(Sum)表征了密实度。将目标区域图像的灰度和相加，黑色为0，白色为255，Sum越大表明密实度越小，Sum越小表明密实度越大。

(10)计算目标区域面积

经步骤(7)填充后，图像已转化为仅显示目标区域的二值图(即黑白图，仅由0和1构成，0为黑色，1为白色)。为了计算目标区域的面积，调用bwarea()函数。如图8数字尺寸为[M,N]的图像，可将其理解为一个矩阵(M是二维矩阵的行数，N是二维矩阵的列数)，也可理解为由M*N个小方块(像素)组成的图像。bwarea()函数即计算图像矩阵中1的总数，即白色方块的个数。

如图8所示，图像总面积为5*8＝40，若其中白色方块个数为10，其他均为黑，则目标区域面积为10。

调用bwarea()函数：

Area＝bwarea(Y)；

此时计算出的Area即为目标区域的方块面积，在获取菌胶团图像时，利用显微镜可获得图片的实际尺寸。则可计算目标区域实际面积，目标区域实际面积＝目标区域方块面积*图片实际面积/图片方块面积。

调用bwarea()函数

Area＝bwarea(Y)；

(11)计算目标区域密实度

在获得目标区域的面积和灰度值总和的前提下，可根据式Avergray＝Sum Gray/Area计算目标区域密实度。

Avergray＝Sum Gray/Area；

(12)获取目标区域的边界

为了计算目标区域的周长，进行圆度计算。通过边缘检测以获取目标区域的轮廓边界。

P＝edge(J1,'log')；％使用log边缘检测算子

figure,imshow(P)；％显示轮廓图

图9(a)给出了第一种目标区域轮廓图；图7(b)给出了第二种目标区域轮廓图。

TEMP＝0；

SUM＝0；

P1＝im2double(P)；

forx＝1:m

for y＝1:n

TEMP＝P1(x,y)；

if TEMP＝＝1；

SUM＝SUM+1；％计算轮廓图内像素数，即计算轮廓长度

end

end

end

由最小外接矩形求得长宽，如图10(a)和图10(b)所示。

三.提取特征信息

(1)计算菌胶团尺寸指标

①菌胶团面积Area；

②菌胶团周长Perimeter；

③等效圆直径

④菌胶团长度length；

⑤菌胶团宽度width；

得出的a、b两种活性污泥絮体的尺度指标：

(2)计算菌胶团形态指标：

①Ro(圆度，表征图形与圆之间的偏差，有效描述不规则边界，圆的Ro为1)

②AR(三维纵横比，表征图像的延伸，延伸越大，AR越大，圆的AR为1)

③Avergray(密实度，表征菌胶团的密室程度)

④FD(fractal dimension，分形维数)表征菌胶团边界的粗糙程度

分形是指事物的形状、形态与组织的分解、分割、分裂与分析，其在一定程度上代表了由局部到整体的对事物的认识过程，而分形维数则是用来描述分形不规则特征的参数。

本发明采用盒子维数，即容量维数(f)分形方法。取边长为ε的小盒子对分形曲线进行覆盖，由于覆盖曲线有各种层次的空洞和缝隙，则有些小盒子为空。计数整个覆盖范围内非空实体盒子数量，其数量记为N(ε)。然后缩小盒子的尺寸ε，所计数非空实体盒子的数量N(ε)随之增大。当ε趋于0时，得到分形维数(FD)。

实际计算中只能取有限的ε，在双对数坐标纸上做logN(ε)～log(1/ε)曲线，然后采用最小二乘法进行拟合，所得直线斜率即所求分形维数FD。

在MATLAB内输入公式得图11(a)和图11(b)两种图像指标信息。

应用fractalfox软件，得到分形维数值fractal，第一种如图12(a)，第二种如图12(b)。

得出的a、b两种活性污泥絮体的形态指标：

(3)计算丝状菌指标：

如图13(a)和图13(b)所示，通过菌丝长度计算，得出的a、b两种活性污泥絮体的菌丝长度FL(fialment length，)：

图像	菌丝长度
		a	1747.94
b	1888.74

Claims (2)

1.一种基于MATLAB的活性污泥絮体分析方法，其特征是，包括以下步骤：

(一)采集图像：

取新鲜活性污泥制成观测玻片，在显微镜下观察，以获取完整活性污泥菌胶团为标准，拍照并采集图像；

(二)图像处理：

(1)读入采集图像；

(2)将RGB图像转化为灰度图像；

(3)去除噪声处理；

(4)图像二值化；为了确定目标菌胶团，对图像进行二值化；

(5)开运算和闭运算，光滑边缘处理；

(6)确定目标菌胶团并进行标记；

选取最大联通区域作为最大的菌胶团，利用八邻域标记算法标记目标菌胶团；

(7)填充；

填充二值图像前景区域内出现的黑色孔洞，以降低对计算结果准确性的影响；

(8)目标区域灰度图；

以“密实度”来表示污泥采集图像中菌胶团密实程度，“密实度”为目标区域内灰度值之和与目标区域面积的比值，图像灰度值的范围为0～255，对应图像中的颜色为从黑到白，越密实的菌胶团图像颜色越深，灰度越小；越松散的菌胶团图像颜色越浅，灰度越大；

(9)计算目标区域灰度图的灰度值之和；

灰度值之和表征了密实度，将目标区域灰度图的灰度值相加，黑色为0，白色为255，灰度值之和越大表明密实度越小，灰度值之和越小表明密实度越大；

(10)计算目标区域面积；

(11)计算目标区域密实度；

在获得目标区域的面积和灰度值总和的前提下，根据公式：密实度＝灰度值之和/菌胶团面积，计算目标区域密实度；

(12)获取目标区域的边界；

为了计算目标区域的周长，进行圆度计算，通过边缘检测以获取目标区域的轮廓边界；

(三)提取特征信息：

(1)计算菌胶团尺寸指标：

①菌胶团面积Area；

②菌胶团周长Perimeter

③等效圆直径

④菌胶团长度length；

⑤菌胶团宽度width；

(2)计算菌胶团形态指标：

①圆度Ro：

$Ro=\frac{4\mathrm{\π}Area}{{\mathrm{Perimeter}}^2},$

②三维纵横比AR：

$AR=1+\frac{4}{\mathrm{\π}}(\frac{lenth}{width}-1),$

③密实度Avergray：

$Averygray=\frac{Sumgray}{Area},$

④分形维数FD；

(3)计算丝状菌指标：菌丝长度FL。

2.根据权利要求1所述的基于MATLAB的活性污泥絮体分析方法，其特征是，所述分形维数采用盒子维数分形方法，取边长为ε的小盒子对分形曲线进行覆盖，计数整个覆盖范围内非空实体盒子数量，其数量记为N(ε)，然后缩小盒子的尺寸ε，所计数非空实体盒子的数量N(ε)随之增大，当ε趋于0时，得到分形维数FD。