CN101354359B - 液体药品中运动可见异物与气泡的检测识别方法 - Google Patents

Abstract

本发明公开了一种液体药品中运动可见异物与气泡的检测识别方法，其特征在于，包括以下步骤：1)对待检测的当前液体药品获取多帧连续图像；2)图像预处理；3)运动目标提取，获得增强后图像；4)运动目标分割：对所述增强后图像运用改进型二维最大熵阈值分割算法，采用灰度和邻域平均灰度构成的二维直方图选择阈值，得到分割后图像；5)运动目标跟踪；6)图像识别与判断：通过跟踪结果得到目标运动轨迹，利用目标轨迹的连续性和方向性识别区分出异物与气泡目标，从而最终判断当前液体药品是否合格。本发明保证算法精度的条件下充分简化了计算，提高了算法的实时性，大大提高了药液中异物与气泡识别的准确率，降低了药品的误检率。

Description

液体药品中运动可见异物与气泡的检测识别方法

技术领域

本发明属于图像处理以及自动化技术领域，涉及一种液体药品中运动可见异物与气泡的检测识别方法。

背景技术

目前，在我国广泛使用的液体药品包括注射针剂、大输液、口服液、滴眼液等，通常采用玻璃瓶或塑料瓶的形式灌装，这些药品在生产过程中由于过滤不佳、容器清洗不佳、封装时碰撞等原因，会使药液中存在玻璃屑、悬浮物、毛发等可见异物，为了检测灌装后的液体药剂中是否混入可见异物，需要对每一瓶药品进行可见异物检测。我国药典自1985年起对此类药品的检测即规定进行澄明度检测和不溶性微粒检查，对药品内的异物和微粒进行严格控制。1999年8月24日，国家药监局正式印发《关于实施<药品生产质量管理规范>有关规定的通知》，明确要求药品生产企业必须在2000年底前通过GMP(good)认证。《中国药典》2005年版已于2005年7月1日起执行，将原“澄明度检查细则和判断标准”修订为“可见异物检查法”，进一步明确了注射剂、滴眼剂等药品必须在符合药品生产质量管理规范(GMP)的条件下生产，产品出厂前应采用适宜的方法逐一检查并同时剔除不合格产品，可见异物是指存在于注射剂、滴眼剂中，在规定条件下目视可以观测到的不溶性物质，其粒径或长度通常大于50um。可见异物检查法有灯检法和光散射法，一般常采用灯检法。目前国内药品可见异物检测大都采用人工灯检的方法，人工手动轻轻旋转或翻转药瓶容器，使药液中的可能存在的可见异物运动(注意不使药液产生气泡)，通过人眼目测进行判断，若存在可见异物将其手动剔除。人工检测时不可必避免的存在效率低，漏检率高，精度低等问题，检测结果的重复性和一致性差，且长时间从事检测工作也给灯检人员的健康带来危害。随着GMP认证的推广，我国医药生产装备和检测设备自动化和智能化水平得到了很大的提高，出现了用于药品可见异物检测的半自动灯检机和全自动灯检机，其中半自动灯检机采用简单的机械传动与光学***实现人工灯检的半自动化，该检测装置能在一定程度上减轻人工劳动，但仍然通过灯检工分析判断，检测精度和准确率没有得到提高。全自动灯检机只有德国、日本、意大利等少数几个国家的几家公司有能力生产，但至今国外的全自动灯检机在我国应用仍然很少，究其原因存在以下问题：国内制药环节中的过滤、包装材料等和国外存在很大差异，致使检测质量低下。但随着国内制药环境的不断改善，外国的先进医药检测设备势必改变国内制药市场的检测现状。面对国外先进成套医药设备的冲击，研发具有自主知识产权的医药检测设备和核心检测识别方法具有十分重要的价值。

在医药液体内异物视觉检测过程中，需用特制的机械装置使沉淀在药液底部的异物带动到液体上部，以便于摄像机成像，但由于机械装置的运动往往使成像后的图像中多少存在有气泡和噪声，而图像中异物目标和气泡识别中所面临的技术难题主要有：

1、异物大小不一，检测分辨率要达到50um，异物目标没有尺寸、形状和纹理等信息，传统的图像处理方法无法应用；

2、由于瓶子表面存在刻度浮雕、斑点，成像后的图像背景复杂，噪声多样；

3、图像中信噪比低，异物目标往往淹没在背景噪声中；

4、生产线检测速度要求高，对识别算法的实时性提出了很高的要求。

发明内容

本发明所要解决的技术问题是，提供一种液体药品中运动可见异物与气泡的检测识别方法，用于在自动化液体药品检测线中迅速检测出具有异物的不合格产品，与成像硬件***配套进行信号分析与处理，为医药检测设备提供一种可靠的检测方法，从而提高药品可见异物检测精度和重复性、彻底解决检测中误检率高的问题，满足国内现有灯检***的性能需求。

本发明为解决上述技术问题所采用的技术方案是：

一种液体药品中运动可见异物与气泡的检测识别方法，其特征在于，包括以下步骤：

1)对待检测的当前液体药品获取多帧连续图像；

2)图像预处理：采用7×7的圆形模板作为结构元素对获取的图像进行高帽(top-hat)形态学滤波，以去除比结构元素小的各类亮噪声和暗噪声，得到去噪图像；

3)运动目标提取：根据异物目标在时域上表现为帧间的运动特性，利用新型的三帧图像差分算法提取图像中的运动目标，其操作过程是：第一帧与第二帧做差分处理，第二帧与第三帧做差分处理，然后通过图像累积增强目标图像，获得增强后图像；

4)运动目标分割：对所述增强后图像运用改进型二维最大熵阈值分割算法，采用灰度和邻域平均灰度构成的二维直方图选择阈值，得到分割后图像；

5)运动目标跟踪：当检测到分割后图像中的可能目标连续3帧出现在一定邻域内，则认定运动轨迹起始；根据图像中异物目标和气泡目标运动的规律性，采用基于匀速直线运动模型的卡尔曼滤波更新目标状态，并由此生成椭圆跟踪门；当落入椭圆跟踪门的量测值只有一个，直接进行目标轨迹更新；当落入椭圆跟踪门的量测值多于一个，目标轨迹通过距离进一步预测值最近的测量值来更新；当目标丢失时，采用以前时刻的滤波值外推预测后几帧目标的可能位置，保持稳定的记忆跟踪能力；

6)图像识别与判断：通过跟踪结果得到目标运动轨迹，利用目标轨迹的连续性和方向性识别区分出异物与气泡目标，从而最终判断当前液体药品是否合格。

所述步骤2)为：

设图像的灰度函数为f(x，y)，结构元素为B(x，y)，(x，y)为像素坐标，则4种基本操作定义如下：

腐蚀：(fΘB)(x，y)＝min{f(x+i，y+j)-B(i，j)|(x+i，y+j)∈D_f；(i，j)∈D_B}；

膨胀： $(f\&CirclePlus;B)(x,y)=\max \{f(x-i,y-j)+B(i,j)|(x-i,y-j)\&Element;D_f;(i,j)\&Element;D_B\};$

开运算：

闭运算： $f\&CenterDot;B=(f\&CirclePlus;B)\mathrm{\&Theta;B};$

式中，D_f和D_B分别是函数f和B的定义域，

引入高帽变换滤波算子：Hat(f)＝f-(foB)；

式中(foB)为结构元素B对图像f进行灰度开运算，选用7×7的圆形模板作为结构元素对原图像进行高帽形态学滤波，得到含有异物目标、气泡的图像，即去噪图像。

所述步骤3)为：

针对去噪图像，按照图像序列对每连续三帧图像中的第一帧与第二帧做差分处理，第二帧与第三帧做差分处理，然后按下式增强目标图像，获得增强后图像：

d(x，y)＝|f₁(x，y)-f₂(x，y)|×|f₂(x，y)-f₃(x，y)|；

其中f₁(x，y)，f₂(x，y)，f₃(x，y)分别为经高帽滤波后的第一帧、第二帧、第三帧图像，|f₁(x，y)-f₂(x，y)|为第一帧与第二帧差分图像，|f₂(x，y)-f₃(x，y)|为第二帧与第三帧差分图像。。

步骤4)包括以下步骤：

(a)将图像分割成一系列子图像；

(b)二维最大熵算法计算出每个子图像的阈值；

(c)根据各子图的阈值，对图像进行分割；

其中阈值的求取方法如下：设图像的尺寸为M×N，图像共有L个灰度级，每个像素点的邻域灰度平均值也为L个灰度级，总的灰度级为L×L个；在二维平面上作出二维灰度直方图，灰度值和邻域灰度平均值形成灰度二元值，记为每个像素点的坐标值为(i，j)，i和j分别代表像素点的灰度值和其邻域的灰度平均值，该灰度组合出现的频率为f_ij，该像素点函数值为联合概率密度函数p_ij：

$p_{\mathrm{ij}}=\frac{f_{\mathrm{ij}}}{M\&times;N},i,j=\mathrm{0,1},...L-1;$

若以(s，t)为分割阈值，对物体O的概率分布p_ij(i＝1，2，...s；j＝1，2，...t)进行归一化操作；定义归一化操作后与物体O的概率分布相关的熵为：

$H\left(O\right)=-\underset{i=0}{\overset{s}{\mathrm{\&Sigma;}}}\underset{j=0}{\overset{t}{\mathrm{\&Sigma;}}}\frac{p_{\mathrm{ij}}}{P_{\mathrm{st}}}\ln \frac{p_{\mathrm{ij}}}{P_{\mathrm{st}}}=-\frac{1}{P_{\mathrm{st}}}\underset{i=0}{\overset{s}{\mathrm{\&Sigma;}}}\underset{j=0}{\overset{t}{\mathrm{\&Sigma;}}}(p_{\mathrm{ij}}\ln p_{\mathrm{ij}}-p_{\mathrm{ij}}\ln P_{\mathrm{st}})=\ln P_{\mathrm{st}}+\frac{H_{\mathrm{st}}}{P_{\mathrm{st}}};$

其中， $P_{\mathrm{st}}=\underset{i=0}{\overset{s}{\mathrm{\&Sigma;}}}\underset{j=0}{\overset{t}{\mathrm{\&Sigma;}}}{p}_{\mathrm{ij}},$ $H_{\mathrm{st}}=-\underset{i=0}{\overset{s}{\mathrm{\&Sigma;}}}\underset{j=0}{\overset{t}{\mathrm{\&Sigma;}}}{p}_{\mathrm{ij}}\ln p_{\mathrm{ij}};$

与背景B的概率分布相关的熵为

$H\left(B\right)=-\underset{i=s+1}{\overset{L-1}{\mathrm{\&Sigma;}}}\underset{j=t+1}{\overset{L-1}{\mathrm{\&Sigma;}}}\frac{p_{\mathrm{ij}}}{1-P_{\mathrm{st}}}\ln \frac{p_{\mathrm{ij}}}{1-P_{\mathrm{st}}}=-\frac{1}{1-P_{\mathrm{st}}}\underset{i=s+1}{\overset{L-1}{\mathrm{\&Sigma;}}}\underset{j=t+1}{\overset{L-1}{\mathrm{\&Sigma;}}}[p_{\mathrm{ij}}\ln p_{\mathrm{ij}}-p_{\mathrm{ij}}\ln (1-P_{\mathrm{st}})];$

$=\ln (1-P_{\mathrm{st}})+\frac{H_{L-1L-1}-H_{\mathrm{st}}}{1-P_{\mathrm{st}}}$

其中， $H_{L-1L-1}=-\underset{i=0}{\overset{L-1}{\mathrm{\&Sigma;}}}\underset{j=0}{\overset{L-1}{\mathrm{\&Sigma;}}}{p}_{\mathrm{ij}}\ln p_{\mathrm{ij}};$

准则函数Ψ(s，t)取为H(O)，H(B)之和，即

$\mathrm{\&Psi;}(s,t)=H\left(O\right)+H\left(B\right)=\ln P_{\mathrm{st}}+\frac{H_{\mathrm{st}}}{P_{\mathrm{st}}}+\ln (1-P_{\mathrm{st}}+\frac{H_{L-1L-1}-H_{\mathrm{st}}}{1-P_{\mathrm{st}}});$

$=\ln P_{\mathrm{st}}(1-P_{\mathrm{st}})+\frac{H_{\mathrm{st}}}{P_{\mathrm{st}}}+\frac{H_{L-1L-1}-H_{\mathrm{st}}}{1-P_{\mathrm{st}}}$

使Ψ(s，t)最大的灰度值(s，t)即是所求的最优阈值t^*，即

$t^{*}=\mathrm{Arg}\left\{\underset{s,t\&Element;G}{\mathrm{Max}}\mathrm{\&Psi;}(s,t)\right\};$

其中Arg{}表示求反函数。

作为改进，所述步骤4)还包括：引入遗传算法对二维最大熵法的分割阀值进行优化，以缩短图像处理时间；具体步骤如下：

初始化种群：随机产生，选择均匀分布；

种群规模：经过分析，种群规模为popsize＝20，最大进化代数为g_max＝100；

适应度函数：直接使用目标函数式Ψ(s，t)作为适应度函数；

编码：由于图像灰度值在0～255之间，采用8位灰度图像，则门限参数满足0≤s，t≤255，故将个体编码为16位二进制码，高8位表示s，低8位表示t；

解码：将16位二进制码的高8位解码为s值，低8位解码为t值；

选择：采用比例选择算子；

交叉算子：采用双点交叉，随机产生的两个交叉点分别位于编码串的前8位和后8位；以交叉概率Pc进行交叉操作；搜索前期选取Pc＝0.8，搜索后期选取Pc＝0.6；

变异算子：变异操作为按位取反；将变异概率Pm定义为类抛物线型变异算子：

$P_m=P_m\left(g\right)=P_{m\_\max }-\frac{4\&times;(P_{m\_\max }-P_{m\_\min })}{g_{\max }^2}{(g-\frac{g_{\max }}{2})}^2$

其中g为进化代数，且g∈[1，g_max]，P_{m_min}＝P_b，P_{m_max}＝10P_{m_min}

P_b为基本变异概率，按下式进行估计：

$P_b\&ap;\frac{175}{\mathrm{popsize}\&times;\sqrt{\mathrm{bits}}};$

其中bits表示个体的位数；当popsize＝20，bits＝16时，P_b≈0.02；

当P_b≈0.02，g_max＝100时，可得：

P_m＝P_m(g)＝0.2-0.00008×(g-50)²；

终止准则：当算法执行到最大代数或群体中的最高适应度值稳定时，算法停止运行，具有最高适应度值的个体即为所求的解；算法终止时，具有最高适应度值的个体作为最佳阈值。

所述步骤5)包括以下步骤：

通过基于改进型二维最大熵阀值分割算法可以确定得到单帧检测结果，将目标用匀速直线运动模型近似；其跟踪***模型为：令状态向量为 $X\left(k\right)={(x\left(k\right),\stackrel{\&CenterDot;}{x}\left(k\right),y\left(k\right),\stackrel{\&CenterDot;}{y}\left(k\right))}^T,$ 测量向量为Z(k)＝[x(k)，y(k)]^T，***方程可写成为

X(k+1)＝ΦX(k)+Gw(k)；

Z(k+1)＝HX(k)+v(k)；

其中 $\mathrm{\&Phi;}=\left[\begin{array}{cccc}1& T& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& T\\ 0& 0& 0& 1\end{array}\right],$ $G=\left[\begin{array}{cc}\frac{T^2}{2}& 0\\ T& 0\\ 0& \frac{T^2}{2}\\ 0& T\end{array}\right],$ $H=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right],$

w(k)，v(k)分别是方差为Q(k)和R(k)的零均值高斯白噪声；式中T为采样周期；轨迹起始：如果检测的可能目标连续3帧出现在3×3邻域内，则轨迹起始；

跟踪门及数据关联：假设k-1时刻的状态向量的估计为

，状态向量的一步预测为：

$\hat{X}\left(k\right|k-1)=\mathrm{\&Phi;}\hat{X}(k-1|k-1);$

则***的新息为：

$d\left(k\right)=Z\left(k\right)-H\hat{X}\left(k\right|k-1);$

新息方差阵为：

S(k)＝HP(k|k-1)H^T+R(k)；

其中，P(k|k-1)是一步预测方差；

令新息向量d(k)的范数为：

g(k)＝d^T(k)S^-1(k)d(k)；

其中g(k)服从

分布，M1为测量维数；

在测量空间中定义一个跟踪门，使得测量以一定的概率分布于跟踪门内；

$\tilde{V}\left(\mathrm{\&gamma;}\right)=[Z:g\left(k\right)\&le;\mathrm{\&gamma;}],$ 其中γ可以通过x²分布表差查得；

如果在某帧处理完后，在跟踪门内只有一个回波，则目标轨迹直接更新；如果跟踪门内有多于一个回波，则目标轨迹通过距离一步预测值最近的测量值来更新；

目标的轨迹通过标准卡尔曼滤波算法来更新

X(k|k-1)＝ΦX(k-1|k-1)

P(k|k-1)＝ΦP(k-1|k-1)Φ^T+GQ(k-1)G^T

K(k)＝P(k|k-1)H^T[HP(k|k-1)H^T+R]^-1；

X(k|k)＝X(k|k-1)+K(k)[Z(k)-H(k)X(k|k-1)]

P(k|k)＝[I-K(k)H]P(k|k-1)

在目标被检测并跟踪上之后，若目标短暂消失，可以根据目标此前的位置信息和运动状态，预测出目标下一步可能的位置，当目标再次出现时，仍可稳定跟踪而不至于丢失目标，具体过程如下：

假设k时刻的状态向量的估计为，则在k+n帧的目标的预测值为：

$\hat{X}(k+n|k)={\mathrm{\&Phi;}}^n\hat{X}\left(k\right|k);$ 选择n<6。

所述步骤6)为：对于目标跟踪得到的各条假设轨迹按以下3条原则来区分轨迹是异物形成的还是噪声或气泡形成的，判断准则如下：

(a)由于图像是在药瓶翻转180度后静止时拍摄的，可见异物是向下运动的，所以其形心纵坐标应该是依次变大的，以左上角为坐标原点；

(b)气泡是向上运动的，运动方向与异物相反；

(c)异物形成的轨迹曲线是光滑的，噪声形成的轨迹是无序的；

因此，如果检测到运动轨迹的形心纵坐标依次变小，说明该轨迹由气泡产生；如果检测到运动轨迹平滑且形心纵坐标依次变大的轨迹，说明有异物存在，判断当前液体药品不合格。

本发明的有益效果有：

本发明采用数学形态滤波获得预处理图像，其算法可通过硬件并行实现，大大提高了处理速度。在运动目标提取中，利用新型三帧图像差分算法很好的克服了瓶侧壁可能存在差异、简单序列图像差分对微小异物目标检测效果不好的缺陷，并且显著提高了输出图像的信噪比，大大简化了后续分割检测算法的难度。采用改进型二维最大熵阀值运动目标分割算法充分利用了图像中目标像素的灰度分布信息和像素间的相关信息，提高了阈值分割的抗噪能力，并通过遗传算法对二维最大熵法的参数进行优化，减少了图像处理时间，大大提高了算法的实时性和鲁棒性。根据图像中运动目标的运动特性，采用基于匀速直线运动模型的卡尔曼滤波更新目标状态，在保证跟踪精度的情况下简化了滤波的计算量。本发明保证算法精度的条件下充分简化了计算，提高了算法的实时性，大大提高了药液中异物与气泡检测识别的准确率，降低了药品漏检率和误检率，可应用于医药、饮料、酒类等各种检测***中，具有广阔的市场前景和实际应用价值。

附图说明

图1是本发明所涉及方法的总流程框图；

图2是本发明中获取的液体药品在瓶内的原始图像；

图3是本发明中目标提取后得到的增强图像；

图4是本发明中运动目标分割后得到的图像；

图5是本发明中运动目标跟踪得到的运动轨迹图。

标号说明：1-气泡轨迹，2-异物轨迹。

具体实施方式

以下将结合图1～5和具体实施例作进一步详细说明。

实施例1：

如图1所示，本发明提出了一种用于液体药品检测的运动可见异物与气泡识别方法。各部分具体实施细节如下：

1、Top-hat(高帽)形态学滤波处理

对获取的图像采用灰度形态学滤波滤除图像噪声。灰度形态学最基本的运算就是腐蚀、膨胀、开运算和闭运算。设图像的灰度函数为f(x，y)，结构元素为B(x，y)，则4种基本操作定义如下；；式中(x，y)是像素坐标；

腐蚀：(fΘB)(x，y)＝min{f(x+i，y+j)-B(i，j)|(x+i，y+j)∈D_f；(i，j)∈D_B}；

膨胀： $(f\&CirclePlus;B)(x,y)=\max \{f(x-i,y-j)+B(i,j)|(x-i,y-j)\&Element;D_f;(i,j)\&Element;D_B\};$

开运算：

闭运算： $f\&CenterDot;B=(f\&CirclePlus;B)\mathrm{\&Theta;B};$

式中，D_f和D_B分别是函数f和B的定义域。腐蚀灰度图像的结果是，比背景暗的部分得到扩张，而比背景亮的部分受到收缩；膨胀灰度图像的结果是，比背景亮的部分得到扩张，而比背景暗的部分受到收缩；开运算使图像的轮廓变得光滑，断开狭窄的间断和消除细的突出物；闭运算同样使图像的轮廓变得光滑，但与开运算相反，它能消除狭窄的间断和长细的鸿沟，消除小的孔洞，并填补轮廓线中的裂痕。

基于灰度形态学的Top-hat变换滤波算子定义为：

[0118]Hat(f)＝f-(fοB)；    (1)

[0119]式中B为结构元素，(fοB)为结构元素B对图像f进行灰度开运算，其结构即是估计出的背景，相减后得到Hat(f)，即为滤除了噪声的目标图像，当然，经过滤波后的图像还包含有一些单点脉冲噪声，它们可以在后面的处理中滤除。由于获取的图像中的噪声为小的亮斑，因此，选用7×7的圆形模板作为结构元素对原图像进行Top-hat形态学滤波，得到含有异物目标、气泡的去噪图像，如图3所示。

2、基于三帧图像差分算法的运动目标提取

针对Top-hat滤波去噪处理后的图像，按照图像序列对第一帧与第二帧做差分处理，第二帧与第三帧做差分处理，然后按下式增强目标图像：

d(x，y)＝|f₁(x，y)-f₂(x，y)|×|f₂(x，y)-f₃(x，y)|   (2)

其中f₁(x，y)，f₂(x，y)，f₃(x，y)分别为经Top-hat滤波后的第一帧、第二帧、第三帧图像，|f₁(x，y)-f₂(x，y)|为第一帧与第二帧差分图像，|f₂(x，y)-f₃(x，y)|为第二帧与第三帧差分图像；

通过上述三帧图像差分算法处理后，d(x，y)中只有运动目标对应的像素点位置不为零。但实际上，由于各种噪声的存在，使得差分图像中在运动目标外的很多像素点位置上也不为零，这就要求后续的处理能有效识别出异物信息。

3、运动目标分割

设图像的尺寸为M×N，图像共有L个灰度级，每个像素点的邻域灰度平均值也为这L个灰度级，总的灰度级为L×L个。在平面上作出二维灰度直方图，灰度值和邻域灰度平均值形成灰度二元值，记为(i，j)，i和j分别代表像素点的灰度值和其邻域的灰度平均值，该灰度组合出现的频率为f_ij，该像素点函数值为联合概率密度函数p_ij。

$p_{\mathrm{ij}}=\frac{f_{\mathrm{ij}}}{M\&times;N},i,j=\mathrm{0,1},...L-1---\left(3\right)$

若以(s，t)为分割阈值，对应于物体O的概率分布p_ij(i＝1，2，...s；j＝1，2，...t)的和应当为1，于是需进行归一化操作。

定义归一化操作后与物体O的概率分布相关的熵为：

$H\left(O\right)=-\underset{i=0}{\overset{s}{\mathrm{\&Sigma;}}}\underset{j=0}{\overset{t}{\mathrm{\&Sigma;}}}\frac{p_{\mathrm{ij}}}{P_{\mathrm{st}}}\ln \frac{p_{\mathrm{ij}}}{P_{\mathrm{st}}}=-\frac{1}{P_{\mathrm{st}}}\underset{i=0}{\overset{s}{\mathrm{\&Sigma;}}}\underset{j=0}{\overset{t}{\mathrm{\&Sigma;}}}(p_{\mathrm{ij}}\ln p_{\mathrm{ij}}-p_{\mathrm{ij}}\ln P_{\mathrm{st}})=\ln P_{\mathrm{st}}+\frac{H_{\mathrm{st}}}{P_{\mathrm{st}}}---\left(4\right)$

其中， $P_{\mathrm{st}}=\underset{i=0}{\overset{s}{\mathrm{\&Sigma;}}}\underset{j=0}{\overset{t}{\mathrm{\&Sigma;}}}{p}_{\mathrm{ij}}---\left(5\right)$

$H_{\mathrm{st}}=-\underset{i=0}{\overset{s}{\mathrm{\&Sigma;}}}\underset{j=0}{\overset{t}{\mathrm{\&Sigma;}}}{p}_{\mathrm{ij}}\ln p_{\mathrm{ij}}---\left(6\right)$

同理，与背景B的概率分布相关的熵为

$H\left(B\right)=-\underset{i=s+1}{\overset{L-1}{\mathrm{\&Sigma;}}}\underset{j=t+1}{\overset{L-1}{\mathrm{\&Sigma;}}}\frac{p_{\mathrm{ij}}}{1-P_{\mathrm{st}}}\ln \frac{p_{\mathrm{ij}}}{1-P_{\mathrm{st}}}=-\frac{1}{1-P_{\mathrm{st}}}\underset{i=s+1}{\overset{L-1}{\mathrm{\&Sigma;}}}\underset{j=t+1}{\overset{L-1}{\mathrm{\&Sigma;}}}[p_{\mathrm{ij}}\ln p_{\mathrm{ij}}-p_{\mathrm{ij}}\ln (1-P_{\mathrm{st}})]---\left(7\right);$

$=\ln (1-P_{\mathrm{st}})+\frac{H_{L-1L-1}-H_{\mathrm{st}}}{1-P_{\mathrm{st}}}$

其中， $H_{L-1L-1}=-\underset{i=0}{\overset{L-1}{\mathrm{\&Sigma;}}}\underset{j=0}{\overset{L-1}{\mathrm{\&Sigma;}}}{p}_{\mathrm{ij}}\ln p_{\mathrm{ij}};$

只要在区域(i＝s+1，...L-1且j＝1，2，..t)和区域(i＝1，2，..s且j＝t+1，...，L-1)中p_ij≈0，以上定义的H(B)就能成立，这样假设有助于节省计算时间。由于许多情况下，非对角线的概率分布可以忽略的，因此作出上述假设是合理的。

准则函数Ψ(s，t)取为H(O)，H(B)之和，即

$\mathrm{\&Psi;}(s,t)=H\left(O\right)+H\left(B\right)=\ln P_{\mathrm{st}}+\frac{H_{\mathrm{st}}}{P_{\mathrm{st}}}+\ln (1-P_{\mathrm{st}}+\frac{H_{L-1L-1}-H_{\mathrm{st}}}{1-P_{\mathrm{st}}})---\left(8\right);$

$=\ln P_{\mathrm{st}}(1-P_{\mathrm{st}})+\frac{H_{\mathrm{st}}}{P_{\mathrm{st}}}+\frac{H_{L-1L-1}-H_{\mathrm{st}}}{1-P_{\mathrm{st}}}$

使Ψ(s，t)最大的灰度值(s，t)即是所求的最优阈值t^*，即

$t^{*}=\mathrm{Arg}\left\{\underset{s,t\&Element;G}{\mathrm{Max}}\mathrm{\&Psi;}(s,t)\right\}---\left(9\right)$

其中Arg{}表示求反函数。

在此，针对(s，t)两个参数问题的最优求解采用遗传算法优化。遗传算法作为一种求解问题的高效并行的全局搜索方法，其主要特点是群体搜索策略和群体中个体之间的信息交换，能在搜索过程中自动获取和积累有关搜索空间的知识，自适应地控制搜索过程以求得最优解或近似最优解。遗传算法利用简单的编码技术和繁殖机制来表现难以用传统搜索方法解决的复杂和非线性的问题。

初始化种群：随机产生，选择均匀分布。

种群规模：经过分析，种群规模为popsize＝20，最大进化代数为g_max＝100。

适应度函数：直接使用目标函数式(6)作为适应度函数。

编码：由于图像灰度值在0～255之间，采用8位灰度图像进行仿真实验，则门限参数满足0≤s，t≤255，故将个体编码为16位二进制码，高8位表示s，低8位表示t。

解码：将16位二进制码的高8位解码为s值，低8位解码为t值。

选择：采用比例选择算子；

交叉算子：采用双点交叉，随机产生的两个交叉点分别位于编码串的前8位和后8位。以事先设定好的交叉概率Pc进行交叉操作。交叉概率过高意味着个体更新较快，可以达到更大的解空间，并能降低取得非最优解的概率，但对已有的较优模式的破坏性也随之增加。交叉概率过低，搜索会因搜索范围的减小而变得迟钝。本文在搜索前期选取Pc＝0.8，搜索后期选取Pc＝0.6。

变异算子：由于采用二进制编码方式，因此变异就是按位取反。变异操作控制了新基因引入种群的比例。变异概率太低，一些有用的基因就不能引入，变异概率太高，即随机变化太多，那么后代就可能失去双亲的优良特性。本实施例将变异概率Pm定义为类抛物线型变异算子：

$P_m=P_m\left(g\right)=P_{m\_\max }-\frac{4\&times;(P_{m\_\max }-P_{m\_\min })}{g_{\max }^2}{(g-\frac{g_{\max }}{2})}^2---\left(10\right)$

其中g为进化代数，且g∈[1，g_max]，P_{m_min}＝P_b，P_{m_max}＝10P_{m_min}；

P_b为基本变异概率，可以按下式进行估计：

$P_b\&ap;\frac{175}{\mathrm{popsize}\&times;\sqrt{\mathrm{bits}}};$

其中bits表示个体的位数。当popsize＝20，bits＝16时，P_b≈0.02；

当P_b≈0.02，g_max＝100时，可得：

P_m＝P_m(g)＝0.2-0.00008×(g-50)²；

终止准则：当算法执行到最大代数或群体中的最高适应度值稳定时，算法停止运行，具有最高适应度值的个体即为所求的解。算法终止时，具有最高适应度值的个体作为最佳阈值。

分块二维最大熵阀值分割算法的基本步骤如下：

(1)将图像分割成一系列子图像；

(2)计算出每个子图像的阀值(用二维最大熵算法来计算每一个阀值)；

(3)根据各子图的阀值，对图像进行分割。

分割后的图像如图4所示。

4、运动目标跟踪

通过基于改进型二维最大熵阀值分割算法可以确定得到单帧检测结果，由于异物和气泡在水中运动的速度较慢，因此目标可以用匀速直线运动模型近似。

跟踪***模型：令状态向量为 $X\left(k\right)={(x\left(k\right),\stackrel{\&CenterDot;}{x}\left(k\right),y\left(k\right),\stackrel{\&CenterDot;}{y}\left(k\right))}^T,$ 测量向量为Z(k)＝[x(k)，y(k)]^T，***方程可写成为：

X(k+1)＝ΦX(k)+Gw(k)         (11)

Z(k+1)＝HX(k)+v(k)          (12)

其中 $\mathrm{\&Phi;}=\left[\begin{array}{cccc}1& T& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& T\\ 0& 0& 0& 1\end{array}\right],$ $G=\left[\begin{array}{cc}\frac{T^2}{2}& 0\\ T& 0\\ 0& \frac{T^2}{2}\\ 0& T\end{array}\right],$ $H=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right],$

w(k)，v(k)分别是方差为Q(k)和R(k)的零均值高斯白噪声。式中T为采样周期。

轨迹起始：如果检测的可能目标连续3帧出现在3×3邻域内，则轨迹起始。

跟踪门及数据关联：假设k-1时刻的状态向量的估计为

状态向量的一步预测为：

$\hat{X}\left(k\right|k-1)=\mathrm{\&Phi;}\hat{X}(k-1|k-1)---\left(13\right)$

则***的新息为：

$d\left(k\right)=Z\left(k\right)-H\hat{X}\left(k\right|k-1)---\left(14\right)$

新息方差阵为：

S(k)＝HP(k|k-1)H^T+R(k)        (15)

其中，P(k|k-1)是一步预测方差。

令新息向量d(k)的范数为：

g(k)＝d^T(k)S^-1(k)d(k)；     (16)

其中g(k)服从

分布，M1为测量维数。

在测量空间中定义一个椭圆体(又称为跟踪门)，使得测量以一定的概率分布于跟踪门内

$\tilde{V}\left(\mathrm{\&gamma;}\right)=[Z:g\left(k\right)\&le;\mathrm{\&gamma;}]---\left(17\right)$

其中γ可以通过x²分布表查得。

如果在某帧处理完后，在跟踪门内只有一个回波，则目标轨迹直接更新；如果跟踪门内有多于一个回波，则目标轨迹通过距离一步预测值最近的测量值来更新。

目标的轨迹通过标准卡尔曼滤波算法来更新

X(k|k-1)＝ΦX(k-1|k-1)

P(k|k-1)＝ΦP(k-1|k-1)Φ^T+GQ(k-1)G^T

K(k)＝P(k|k-1)H^T[HP(k|k-1)H^T+R]^-1        (18)

X(k|k)＝X(k|k-1)+K(k)[Z(k)-H(k)X(k|k-1)]

P(k|k)＝[I-K(k)H]P(k|k-1)

图5给出了目标跟踪得到的运动轨迹图。

在目标被检测并跟踪上之后，若由于某些光学因素使目标短暂消失，可以根据目标此前的位置信息和运动状态，预测出目标下一步可能的位置，当目标再次出现时，仍可稳定跟踪而不至于丢失目标。这样做的目的大大降低了将噪声误判为异物的可能性。由于异物在液体中的运动速度较慢，可以通过上边所述的基于匀速直线运动模型的目标跟踪得出目标状态向量的滤波值，外推预测得到目标的可能位置。

假设k时刻的状态向量的估计为

，则在k+n帧的目标的预测值为：

$\hat{X}(k+n|k)={\mathrm{\&Phi;}}^n\hat{X}\left(k\right|k)---\left(19\right)$

在实际检测中应用时，随着时间的增加及可能物体运动的变化，过去的更总跟踪数据与未来的情况越来越不相关，随着n的增大，预测精度会下降。这里一般选择n<6。

5、目标识别与判断

这一步主要是利用一些图像识别的技术，对于目标跟踪得到的各条假设轨迹按以下3条原则来区分轨迹是异物形成的还是噪声或气泡形成的，判断准则如下：

(1)由于图像是在药瓶翻转180度后静止时拍摄的，可见异物是向下运动的，所以其形心纵坐标应该是依次变大的(以左上角为坐标原点)；

(2)气泡是向上运动的，运动方向与异物相反；

(3)异物形成的轨迹曲线是光滑的，噪声形成的轨迹是无序的。

因此，如果检测到运动轨迹的形心纵坐标依次变小，说明该轨迹由气泡产生；如果检测到运动轨迹平滑且形心纵坐标依次变大的轨迹，说明有异物存在，记录判断结果，并向***发出剔除信号。异物和气泡的轨迹及识别结果如图5所示。

Claims (7)

1.一种液体药品中运动可见异物与气泡的检测识别方法，其特征在于，包括以下步骤：

1)对待检测的当前液体药品获取多帧连续图像；

2)图像预处理：采用7×7的圆形模板作为结构元素对获取的图像进行高帽形态学滤波，以去除比结构元素小的各类亮噪声和暗噪声，得到去噪图像；

3)运动目标提取：根据异物目标在时域上表现为帧间的运动特性，利用新型的三帧图像差分算法提取图像中的运动目标，其操作过程是：第一帧与第二帧做差分处理，第二帧与第三帧做差分处理，然后通过图像累积增强目标图像，获得增强后图像；

4)运动目标分割：对所述增强后图像运用改进型二维最大熵阈值分割算法，采用灰度和邻域平均灰度构成的二维直方图选择阈值，得到分割后图像；

5)运动目标跟踪：当检测到分割后图像中的可能目标连续3帧出现在一定邻域内，则认定运动轨迹起始；根据图像中异物目标和气泡目标运动的规律性，采用基于匀速直线运动模型的卡尔曼滤波更新目标状态，并由此生成椭圆跟踪门；当落入椭圆跟踪门的量测值只有一个，直接进行目标轨迹更新；当落入椭圆跟踪门的量测值多于一个，目标轨迹通过距离进一步预测值最近的测量值来更新；当目标丢失时，采用以前时刻的滤波值外推预测后几帧目标的可能位置，保持稳定的记忆跟踪能力；

6)图像识别与判断：通过跟踪结果得到目标运动轨迹，利用目标轨迹的连续性和方向性识别区分出异物与气泡目标，从而最终判断当前液体药品是否合格。

2.如权利要求1所述的液体药品中运动可见异物与气泡的检测识别方法，其特征在于，所述步骤2)为：

设图像的灰度函数为f(x，y)，结构元素为B(x，y)，(x，y)为像素坐标，则4种基本操作定义如下：

腐蚀：

膨胀：

开运算：

闭运算：

式中，D_f和D_B分别是函数f和B的定义域，

引入高帽变换滤波算子：Hat(f)＝f-(fοB)；

式中(fοB)为结构元素B对图像f进行灰度开运算，选用7×7的圆形模板作为结构元素对原图像进行高帽形态学滤波，得到含有异物目标、气泡的图像，即去噪图像。

3.如权利要求1所述的液体药品中运动可见异物与气泡的检测识别方法，其特征在于，所述步骤3)为：

针对去噪图像，按照图像序列对每连续三帧图像中的第一帧与第二帧做差分处理，第二帧与第三帧做差分处理，然后按下式增强目标图像，获得增强后图像：

d(x，y)＝|f₁(x，y)-f₂(x，y)|×|f₂(x，y)-f₃(x，y)|

其中f₁(x，y)，f₂(x，y)，f₃(x，y)分别为经高帽滤波后的第一帧、第二帧、第三帧图像，|f₁(x，y)-f₂(x，y)|为第一帧与第二帧差分图像，|f₂(x，y)-f₃(x，y)|为第二帧与第三帧差分图像。

4.如权利要求1所述的液体药品中运动可见异物与气泡的检测识别方法，其特征在于，所述步骤4)包括以下步骤：

(a)将图像分割成一系列子图像；

(b)二维最大熵算法计算出每个子图像的阈值；

(c)根据各子图的阈值，对图像进行分割；

其中阈值的求取方法如下：设图像的尺寸为M×N，图像共有L个灰度级，每个像素点的邻域灰度平均值也为L个灰度级，总的灰度级为L×L个；在二维平面上作出二维灰度直方图，灰度值和邻域灰度平均值形成灰度二元值，记为每个像素点的坐标值为(i，j)，i和j分别代表像素点的灰度值和其邻域的灰度平均值，该灰度组合出现的频率为f_ij，该像素点函数值为联合概率密度函数p_ij：

若以(s，t)为分割阈值，对物体O的概率分布p_ij(i＝1，2，...s；j＝1，2，...t)进行归一化操作；定义归一化操作后与物体O的概率分布相关的熵为：

其中，

与背景B的概率分布相关的熵为

其中，

准则函数Ψ(s，t)取为H(O)，H(B)之和，即

使Ψ(s，t)最大的灰度值(s，t)即是所求的最优阈值t^*，即

其中Arg{}表示求反函数。

5.如权利要求4所述的液体药品中运动可见异物与气泡的检测识别方法，其特征在于，所述步骤4)还包括：引入遗传算法对二维最大熵法的分割阀值进行优化，以缩短图像处理时间；具体步骤如下：

初始化种群：随机产生，选择均匀分布；

种群规模：经过分析，种群规模为popsize＝20，最大进化代数为g_max＝100；

适应度函数：直接使用目标函数式Ψ(s，t)作为适应度函数；

编码：由于图像灰度值在0～255之间，采用8位灰度图像，则门限参数满足0≤s，t≤255，故将个体编码为16位二进制码，高8位表示s，低8位表示t；

解码：将16位二进制码的高8位解码为s值，低8位解码为t值；

选择：采用比例选择算子；

交叉算子：采用双点交叉，随机产生的两个交叉点分别位于编码串的前8位和后8位；以交叉概率Pc进行交叉操作；搜索前期选取Pc＝0.8，搜索后期选取Pc＝0.6；

变异算子：变异操作为按位取反；将变异概率Pm定义为类抛物线型变异算子：

其中g为进化代数，且g∈[1，g_max]，P_{m_min}＝P_b，P_{m_max}＝10_{Pm_mm}

P_b为基本变异概率，按下式进行估计：

其中bits表示个体的位数；当popsize＝20，bits＝16时，P_b≈0.02；

当P_b≈0.02，g_max＝100时，可得：

P_m＝P_m(g)＝0.2-0.00008×(g-50)²；

终止准则：当算法执行到最大代数或群体中的最高适应度值稳定时，算法停止运行，具有最高适应度值的个体即为所求的解；算法终止时，具有最高适应度值的个体作为最佳阈值。

6.如权利要求1所述的液体药品中运动可见异物与气泡的检测识别方法，其特征在于，所述步骤5)包括以下步骤：

通过基于改进型二维最大熵阀值分割算法可以确定得到单帧检测结果，将目标用匀速直线运动模型近似；其跟踪***模型为：令状态向量为

测量向量为Z(k)＝[x(k)，y(k)]^T，***方程可写成为

X(k+1)＝ΦX(k)+Gw(k)；

Z(k+1)＝HX(k)+v(k)；

其中

w(k)，v(k)分别是方差为Q(k)和R(k)的零均值高斯白噪声；式中T为采样周期；

轨迹起始：如果检测的可能目标连续3帧出现在3×3邻域内，则轨迹起始；

跟踪门及数据关联：假设k-1时刻的状态向量的估计为

状态向量的一步预测为：

则***的新息为：

新息方差阵为：

S(k)＝HP(k|k-1)H^T+R(k)；

其中，P(k|k-1)是一步预测方差；

令新息向量d(k)的范数为：

g(k)＝d^T(k)S^-1(k)d(k)；

其中g(k)服从

分布，M1为测量维数；

在测量空间中定义一个跟踪门，使得测量以一定的概率分布于跟踪门内；

其中γ可以通过X²分布表差查得；

如果在某帧处理完后，在跟踪门内只有一个回波，则目标轨迹直接更新；如果跟踪门内有多于一个回波，则目标轨迹通过距离一步预测值最近的测量值来更新；

目标的轨迹通过标准卡尔曼滤波算法来更新：

X(k|k-1)＝ΦX(k-1|k-1)

P(k|k-1)＝ΦP(k-1|k-1)Φ^T+GQ(k-1)G^T

K(k)＝P(k|k-1)H^T[HP(k|k-1)H^T+R]^-1；

X(k|k)＝X(k|k-1)+K(k)[Z(k)-H(k)X(k|k-1)]

P(k|k)＝[I-K(k)H]P(k|k-1)

在目标被检测并跟踪上之后，若目标短暂消失，可以根据目标此前的位置信息和运动状态，预测出目标下一步可能的位置，当目标再次出现时，仍可稳定跟踪而不至于丢失目标，具体过程如下：

假设k时刻的状态向量的估计为

则在k+n帧的目标的预测值为：

选择n<6。

7.如权利要求1～6任一项所述的液体药品中运动可见异物与气泡的检测识别方法，其特征在于，所述步骤6)为：对于目标跟踪得到的各条假设轨迹按以下3条原则来区分轨迹是异物形成的还是噪声或气泡形成的，判断准则如下：

(a)由于图像是在药瓶翻转180度后静止时拍摄的，可见异物是向下运动的，所以其形心纵坐标应该是依次变大的，以左上角为坐标原点；

(b)气泡是向上运动的，运动方向与异物相反；

(c)异物形成的轨迹曲线是光滑的，噪声形成的轨迹是无序的；

因此，如果检测到运动轨迹的形心纵坐标依次变小，说明该轨迹由气泡产生；如果检测到运动轨迹平滑且形心纵坐标依次变大的轨迹，说明有异物存在，判断当前液体药品不合格。

Images