0 引言

海洋是资源的宝库,自然界的发展、人类的生存与进步都直接受益于海洋。然而,作为海洋生态环境中最重要的初级生产者,海洋浮游微生物的种类和数量直接影响到海洋总体生产力的水平。近年来,我国海洋环境污染日益加重,赤潮频频发生^[1-2],爆发规模大,持续时间长,不断有赤潮危害升级事件报道^[3-4],因此掌握海洋浮游生物数量分布、生理特点和繁殖规律,分析引发赤潮机制,通过建立赤潮预测模型开展海洋管理,已经成为生物学家、海洋学家亟待研究的课题,其中如何实时准确判定赤潮优势藻种这一关键问题亟待解决。

以往的浮游生物分类主要依赖于专业知识丰富的藻类专家借助显微观察技术对未知种群进行识别,这也是目前最基本、最传统的鉴定方法。但许多有害藻种与无害藻种生理形态差异不大,故需要分类者具有丰富的生物学知识和分类经验,且人为干预多,精度低,远不能满足对海量样品和现场及原位实时测量的需要,迄今我国尚缺乏大尺度、准确实时的检测手段。

近年来,随着水下光学成像技术的发展,为了满足海洋勘查需要,浮游微生物野外原位成像系统不断涌现^[5],如浮游生物录像记录仪VPR^[6]、流式影像术FlowCAM ^[7]、流式细胞仪CAM ^[8]、视频荧光分析仪VFA^[9]等已被广泛应用于水下图像的获取;有关藻类鉴定与定量检测的方法也发展迅速^[10],诸如分子探针^[11]、酶联免疫分析^[12]等接触式鉴定技术可靠性较高,但只能适应较小的空间尺度;与之不同,紫外检测^[13-14]、荧光检测^[15]等高效液相色谱法属于非接触式检测方法,可实现大尺度快捷检测,但对藻种识别较为粗糙,有时只能区分到大的门类,不可细分。相比这些检测方法,以专家知识为判定标准,计算机、显微镜辅助的图像分析和模式识别技术可以形象、直观地完成藻种参数测量,可精细鉴定到种属,是藻种分类鉴定的一个重要发展方向。

基于形态差异的海洋浮游微生物图像识别研究始于20世纪90年代,随着图像采集、图像处理、计算机视觉等技术的发展,近年来取得了大量的进展^{[5, 16-17]}。但论其研究范围,大都围绕特定门类(比如硅藻门、甲藻门等)或少数几个藻种展开,涉及多种门类的多类别识别研究非常有限;论其特征提取,传统的方法均是对藻种尺寸、纹理等参数的简单测量,几乎未公开报道过有关顶刺、腹孔等生理形态方面的目标检测和特征提取方法。

本研究将结合藻种两级生物形态特征探讨大类别多形态藻种显微图像识别方法。

1 赤潮藻显微图像识别方法

本研究所用的藻类样品均采集自我国南海赤潮高发水域,通过对多类赤潮藻生物学特征分析和归纳,得出细胞轮廓和壳面花纹是赤潮藻种分类的一级特征,角毛、顶刺等突起部位、腹孔等凹入部位和横纵沟等甲板排列规律是可用于分类的二级特征,准确有效地提取以上两级特征在藻种机器识别渐进过程中占着举足轻重的地位。

据此,本研究制定了以生物形态学为科学依据的分级鉴定分类方法,如图 1所示。核心步骤主要包括:赤潮藻细胞目标精确定位、两级特征选择与提取、特征降维和分类器构建。

2 藻类图像分割

本研究实验图像是通过图像采集子系统触发CCD显微照相系统以3帧/s的曝光频率对海水样本进行采集并获取的。但在图像采集过程中,流动室中海水杂质、未知悬浮颗粒等沉积和细胞高速运动都会造成随机干扰,为了不使这些噪声影响藻类分析和识别,需要对采集后图像做预处理操作,以剔除无价值的背景信息。 传统的阈值分割方法忽略了图像中的空间信息,噪声分散在细胞体四周,且赤潮藻细胞尺寸一般大于10μm,在图像中占最大连通面积,这是区分目标和杂质的重要线索。本研究将目标的连通性引入阈值分割过程,采用多方向投影积分的方法滤除背景中不连续的杂质沉积区域,提取完整的细胞体。

首先将目标的二值图像分别在 (0,π/8,π/4,3π/8,π/2,5π/8,3π/4,5π/8)8个方向上投影积分,获取目标像素宽度直方图^[18],通过保留各投影结果中具有最大像素宽度总和的区域来获取目标在原始图像中的位置候选区,剔除虚假目标区。最后,对投影积分结果进行形态学运算,能够获取一个平滑的目标模板,用原图像与该模板比对提取目标区域。

为了后续精确提取藻种二级形态特征,本研究进一步对生理细节目标实现了精细分割。(1)对藻种目标进行分割,以能使像素宽度面积分布图取得极大值的像素宽度驻点确定最佳结构元尺寸^[19],结合滤波、连通区域合并等形态学运算可检测二级目标——顶刺;(2)采用融合特征灰度和约束标记的分水岭变换^[20]可以检测二级目标——横沟。图 2与图 3为顶刺和横沟的分割结果。

3 特征提取

本研究对大量赤潮藻进行了分析,发现赤潮藻类的几何形状和壳面纹理规律存有较大差别,因此结合藻类本身的特点,提取藻类细胞的偏心率、体态比、矩形度、形状参数、凸率、圆形度、球状型7个不变矩和纹理特征共八类特征作为藻类的一级特征。

然而,赤潮藻是形态会随目标游动和周围环境发生改变的三维活体目标,要精确识别这种类别多、形态丰富的目标,还需选取同一类别在不同生长时期和环境中之间保持较高稳定性且符合人类视觉的主观理解的特征。本研究选择横沟和顶刺两类二级特征,通过定义新的特征形式化描述,提取局部形态特征。

3.1 形状特征

几何形状特征蕴含有大量的信息,因此外形是藻学分类学家人工鉴定的最基本、最有效的依据。

(1) 偏心率

$e=c/a,$

(1)

式中:c为与区域等效的椭圆焦距长度;a为与区域等效的椭圆主轴长度。偏心率描述了区域与矩形的偏离程度。

(2) 体态比

${{C}_{1}}=W/H,$

(2)

式中,W和H分别为目标区域最小外接矩形的宽度与高度。体态比描述了区域的紧凑性。

(3) 矩形度

$R=A/\left( W*H \right),$

(3)

式中,A为目标区域的面积。矩形度描述了目标在外接矩形中的充满程度。

(4) 形状参数

${{C}_{2}}={{L}^{2}}/4\pi A,$

(4)

式中,L为目标区域的周长。形状参数描述了目标的密集度。

(5) 凸率

$CR=\frac{A\prime }{\sum\limits_{x=1}^{m}{\sum\limits_{y=1}^{n}{k\left( x,y \right)}}},$

(5)

式中：$k\left( x,y \right)=\left\{ \begin{align} & 1,\left( x,y \right)\in 目标区域凸包部分 \\ & 0,\left( x,y \right)\notin 目标区域凸包部分 \\ \end{align} \right.;$ A′为目标区域凸包面积。凸率描述了边界不规则特性。

(6) 圆度

$T=4\pi A/{{L}^{2}},$

(6)

圆度描述了区域与圆形的相似程度。

(7) 球状型

$S={{r}_{i}}/{{r}_{c}},$

(7)

式中,r_i,r_c

分别表示区域内、外切圆的半径。球状型描述了区域与球形的符合程度。

(8) 不变矩

HU M K最早使用图像的几何矩的非线性组合推导出下列7项不变矩^[21],图像发生尺度缩放、旋转和平移时,其值具有不变性:

${{\varphi }_{1}}={{\eta }_{20}}+{{\eta }_{02}},$

(8)

${{\varphi }_{2}}={{({{\eta }_{20}}-{{\eta }_{02}})}^{2}}+4\eta _{11}^{2},$

(9)

${{\phi }_{3}}={{\left( {{\eta }_{30}}-3{{\eta }_{12}} \right)}^{2}}+{{\left( 3{{\eta }_{21}}-{{\eta }_{30}} \right)}^{2}},$

(10)

${{\varphi }_{4}}={{\left( {{\eta }_{30}}+3{{\eta }_{12}} \right)}^{2}}+{{\left( {{\eta }_{21}}+{{\eta }_{03}} \right)}^{2}},$

(11)

$\begin{align} & {{\varphi }_{5}}=\left( {{\eta }_{30}}-3{{\eta }_{12}} \right)\left( {{\eta }_{30}}+{{\eta }_{12}} \right)[{{({{\eta }_{30}}+{{\eta }_{12}})}^{2}}- \\ & 3{{\left( {{\eta }_{21}}+{{\eta }_{03}} \right)}^{2}}]+\left( 3{{\eta }_{21}}-{{\eta }_{03}} \right)\left( {{\eta }_{21}}+{{\eta }_{03}} \right) \\ & \left[ 3{{\left( {{\eta }_{30}}+{{\eta }_{12}} \right)}^{2}}-{{\left( {{\eta }_{21}}+{{\eta }_{03}} \right)}^{2}} \right], \\ \end{align}$

(12)

$\begin{align} & {{\varphi }_{6}}=\left( {{\eta }_{20}}-{{\eta }_{02}} \right)[{{\left( {{\eta }_{30}}+{{\eta }_{12}} \right)}^{2}}-{{\left( {{\eta }_{21}}+{{\eta }_{03}} \right)}^{2}}+ \\ & 4{{\eta }_{11}}\left( {{\eta }_{30}}+{{\eta }_{12}} \right)\left( {{\eta }_{21}}+{{\eta }_{03}} \right)], \\ \end{align}$

(13)

$\begin{align} & {{\varphi }_{7}}=\left( 3{{\eta }_{21}}-{{\eta }_{03}} \right)\left( {{\eta }_{30}}+{{\eta }_{12}} \right)[{{\left( {{\eta }_{30}}+{{\eta }_{12}} \right)}^{2}}- \\ & 3{{\left( {{\eta }_{21}}+{{\eta }_{03}} \right)}^{2}}]+\left( {{\eta }_{30}}-3{{\eta }_{12}} \right)({{\eta }_{21}}+ \\ & {{\eta }_{03}}\left[ 3{{\left( {{\eta }_{30}}+{{\eta }_{12}} \right)}^{2}}-{{\left( {{\eta }_{21}}+{{\eta }_{03}} \right)}^{2}} \right], \\ \end{align}$

(14)

式中,η_pq为归一化的(p+q)阶中心矩。

由于上述不变矩的幅值范围很大,同时考虑到可能取负值的情形,为了便于比较,本研究对所得数据取对数进行压缩后再取绝对值,实际上采用的不变矩为

$\phi {{\prime }_{i}}=\left| l{{g}_{_{i}}}^{2} \right|。$

(15)

3.2 纹理特征

纹理反映了一个区域中像素灰度级的空间分布的属性,是图像分析的重要指标之一。本研究基于灰度共生矩阵提取一级纹理特征。灰度共生矩阵P 描述了图像在θ 方向上,相距为d 的两个像素灰度值分别为i 和j 的概率,反映了图像灰度对于不同方向、相邻间隔和变化幅度的综合信息^[22]。其数学表达式为

$\begin{align} & p\left( i,j,d,\theta \right)=\{\left[ \left( x,y \right),\left( x+dx,y+dy \right) \right]|f\left( x,y \right)=i, \\ & f\left( x+dx,y+dy \right)=j;x=0,1,\cdots ,{{N}_{x}}-1; \\ & y=0,1,\ldots ,{{N}_{y}}-1\}, \\ \end{align}$

(16)

式中:L 表示图像的灰度级;i,j 分别表示像素的灰度;N_x,N_y 表示图像对应的行列数。

本试验选取目标图像在0,π/4,π/2,3π/4 四个方向上位置关系为1的灰度共生矩阵,由此计算以下4个最常用的二阶统计量。

(1) 角二阶矩

$ASM=\sum\limits_{i=1}^{{{N}_{g}}}{\sum\limits_{j=1}^{{{N}_{g}}}{{{P}^{2}}\left( i,j \right)}},$

(17)

用来度量纹理的粗细程度及灰度分布的均匀性。

(2) 惯性矩

$CON=\sum\limits_{i=1}^{{{N}_{g}}}{{{\left| i-j \right|}^{2}}}\sum\limits_{i=1}^{{{N}_{g}}}{\sum\limits_{j=1}^{{{N}_{g}}}{P\left( i,j \right)}},$

(18)

描述了纹理沟纹的清晰程度。

(3) 相关性

$COR=\frac{\sum\limits_{i=1}^{{{N}_{g}}}{\sum\limits_{j=1}^{{{N}_{g}}}{ijP\left( i,j \right)-{{\mu }_{x}}{{\mu }_{y}}}}}{{{\sigma }^{2}}_{x}{{\sigma }^{2}}_{y}},\text{ }$

(19)

其中,μ_x、 μ_y、σ_x和σ_y分别是${{p}_{x}}\left( i \right)=\sum\limits_{k=1}^{{{N}_{g}}}{p\left( i,k \right)}$ 和${{p}_{y}}\left( i \right)=\sum\limits_{k=1}^{{{N}_{g}}}{p\left( k,y \right)}$的均值和方差。该特征描述了灰度共生矩阵元素在行或列方向上的相似程度。

(4) 熵

$H=-\sum\limits_{i=1}^{m}{\sum\limits_{j=1}^{n}{P\left( i,j \right)\ln p\left( i,j \right)}}。$

(20)

该特征度量了目标纹理的随机性。

3.3 形态特征

不同的藻种也可能具有相同的形状和纹理特征,因此,本研究基于藻种局部细节目标分割结果,根据藻种形态特征与机器识别特征对应关系,又定义了以下4类与顶刺和横沟直接相关的二级特征,这些特征参数均为度量性比值,不难证明其具有平移、尺度缩放和旋转不变性。

(1) 顶刺与细胞面积比

$D{{C}_{1}}=\frac{{{S}_{2}}}{{{S}_{1}}},$

(21)

式中,S₁ 和S₂ 分别为细胞体与顶刺区域的面积。

(2) 顶刺分长轴比

$D{{C}_{2}}=\frac{{{L}_{2}}}{{{L}_{1}}},$

(22)

式中,L₁ 和L₂ 分别为细胞与顶刺区域长轴的长度。 (3) 横沟与细胞面积比

$H{{G}_{1}}=\frac{{{S}_{2}}}{{{S}_{1}}},$

(23)

式中,S₁ 和S₂分别为细胞体与横沟区域的面积。

(4) 细胞上、下锥部长度比

$H{{G}_{2}}=\frac{{{L}_{1}}}{{{L}_{2}}},\text{ }$

(24)

式中,L₁和L₂分别为横沟区域质心将细胞长轴分为上、下两部分长轴的长度。

4 分类器设计

试验中选取血红阿卡藻、叉状角藻、东海原甲藻、梭角藻等15种藻种共3000幅图像(训练样本2400幅,测试样本600幅),每个目标图像提取34维的特征向量(7类形状特征、7类不变矩特征、16类纹理特征、4类形态特征),采用主成分分析将特征集合降为10维,最终特征集包括3000个10维的特征向量。

本研究采用带惩罚参数C的支持向量机分类器对显微图像进行识别。支持向量机^[23-24](support vector machines,SVM)是基于统计学习^[25]和结构风险最小化原理^[26]建立的二分类方法,能较好地解决小样本、维数灾难、局部最小点、过学习和欠学习等问题,具有坚实的理论模型。本研究选用LIBSVM程序库,以一对一方法构建多类别分类模型,为K类目标中的任意两类构建两类分类器;在分类时采用投票决策机制,测试样本所属类由得票数最多的类决定;采用网格搜索与交叉验证^[27]相结合的参数选择工具,对惩罚系数C和RBF核参数γ进行了优化选择:

误差惩罚参数C=2048,

RBF径向基函数中γ=0.5 。

5 试验结果

根据以上训练参数选取2400幅训练图像对SVM支持向量机进行训练,并得到最优结构的分类模型,表 1是对600幅测试图像的分类识别结果。

分析表 1可知,对于15类藻种600幅待识别目标图像,平均识别率可达80.33％。其中,反曲原甲藻、梭角藻和叉状角藻的辨识度高,识别结果较好,拒识数均为0,正是由于这些藻种目标的凸率参数和局部形态等特征与其他藻种差别迥异,也验证了本研究定义的局部特征的有效性。另外,个别具有特殊形态的样本对分类识别影响较大,加之本试验待识别的藻种种类比以往报道要多,这是目前识别率保持在80％左右的根本原因。

6 结语

本研究介绍了一种多类别赤潮藻细胞显微图像识别方法,首先采用多方向投影积分的方法定位分割出藻类细胞,进一步对顶刺和横沟细节区域实现了精细分割;将藻种生物形态学特征与机器识别特征相对应,以形状特征、不变矩、纹理和局部形态特征构造了两级生物形态不变特征集合;使用SVM训练识别,试验验证了特征集的稳定性和有效性,结果表明该方法能克服赤潮藻细胞不同环境和生长周期造成的个体和形态差异,可对不同成像条件下拍摄的15类藻种细胞显微图像完成快速且有效的识别。识别算法对于圆形度、矩形度比较高的藻种辨识度不高,这些问题将在以后的研究中改进。由于二级形态特征在分类识别中有良好的效果,下一步研究可以继续检测和提取腹孔、纵沟等其他类型特征。