CN106228529A - 一种激光散斑图像处理分析方法 - Google Patents

Abstract

一种激光散斑图像处理分析方法，用于提取出激光散斑图像中的不同频率的信息，包括如图所示的8个步骤，该方法解决了现有通过计算衬比提取图像特征的方法由于成像对象不均匀性以及图像噪声，无法分离特征层，导致提取图像特征及测量相关参数的精度低的问题，提供了一种基于自适应噪声的完备经验模态分解的激光散斑图像处理分析方法，由于使用了自适应噪声的完备经验模态分解(Complete Ensemble Empirical Mode Decomposition With Adaptive Noise，CEEMDAN)算法，可以较好地解决图像中的噪声干扰，从而获得图像中不同频率的信息能够剔除图像中的高频噪声，分离特征层提取出图像中不同频率的信息。

Description

一种激光散斑图像处理分析方法

技术领域

本发明涉及图像处理分析方法，属于计算机图像处理领域。

背景技术

激光散斑血流成像(Laser Speckle Flowgraphy，LSFG)技术是近些年发展起来的一种新型的血流波动监测技术、具有无需扫描即可全场成像的特点。其原理是：当使用激光照射成像表面(如生物组织),入射的相干光被组织中的散射颗粒(如血流中的红细胞)所散射,而散射光经过随机干涉便形成了明暗相间的图样,这种图样被称为散斑(speckle)。组织中散射颗粒的运动会引起散射光光强信号的相移,进而使得散斑图样的光强发生变化。当被观测物体中的散射粒子保持静止，则其后向散射光形成的散斑图案也保持静止，并具有明暗相间的高对比灰度起伏；当被观测物体中的散射粒子处于运动状态，其后向散射光形成的散斑图案强度分布会不断波动，则在图像采集设备的有限曝光时间内强度起伏被平均，运动越剧烈，积分平均作用越显著。由此，激光散斑图像中包含有散斑图案波动频率的信息，也就是被观测区域散射子的运动速度的信息。

目前激光散斑图像的处理分析主要有空间衬比分析(Laser Speckle SpatialContrast Analysis，LSSCA)和时间衬比分析(Laser Speckle Temporal ContrastAnalysis，LSTCA),如：“一种激光散斑血流成像分析方法”(中国发明专利CN101485565)，“一种激光散斑血流成像处理系统及方法”(中国发明专利CN102357033)，“一种激光散斑血流成像衬比分析方法”(中国发明专利CN102429650)等，都是通过计算衬比来反映出激光散斑图像各区域血流速度的快慢。空间衬比分析方法是针对于一帧图像，取Ns×Ns大小的空间滑动窗，用滑动窗内像素灰度值的标准偏差除以像素灰度值的平均值，计算得到衬比(对比)值Ks并将其赋值于中心像素，当滑动窗滑过整幅图像，得到对应的空间衬比图。时间衬比分析方法针对于同一组图像中的多帧图像，利用多帧图像同一位置的像素灰度值计算衬比Kt,并将其赋值于相同位置像素，则得到对应的时间衬比图。

在空间衬比图和时间衬比图中，不同区域具有不同的衬比值，也就表明不同区域具有不同的速度信息。但是利用空间衬比分析方法会损失激光散斑图像的空间分辨率，同样，时间衬比分析方法会损失激光散斑图像的时间分辨率。同时，激光散斑图像由于成像对象的不均匀性以及图像噪声等导致最终图像中包含有高频的噪声信息，这两种分析方法都无法去除或削弱噪声的影响，无法分离特征层，造成最终的处理得到的图像特征或参数精度低。

因此上述的空间衬比分析(Laser Speckle Spatial Contrast Analysis，LSSCA)和时间衬比分析(Laser Speckle Temporal Contrast Analysis，LSTCA)技术方案对激光散斑图像的处理都不能很好的解决图像中噪声剔除的问题。

发明内容

本发明目的是为了解决现有通过计算衬比提取图像特征的方法由于成像对象不均匀性以及图像噪声，无法分离特征层，导致提取图像特征及测量相关参数的精度低的问题，提供了一种基于自适应噪声的完备经验模态分解的激光散斑图像处理分析方法，能够剔除图像中的高频噪声，分离特征层提取出图像中不同频率的信息。

一种激光散斑图像处理分析方法，用于提取出激光散斑图像中的不同频率的信息，其特征在于，包括以下步骤：

步骤一，使用激光散斑血流成像装置采集散斑图f(x,y)，该图像的尺寸为m×n；

步骤二，将散斑图按行列分割，得到m行×n列数据，将每一行数据分别作为分解对象；

步骤三，使用图像的自适应噪声的完备经验模态分解CEEMDAN算法对第一行分解对象进行分解，得到该行的加总平均值并将其定义为目标信号x[n]的IMF1(n)，最终分解得到对应该行的各级加总平均值IMF及一个残余值R；

步骤四，分解完第一行的分解对象后，继续分解下一行的分解对象，直至分解完整幅散斑图的每一行的分解对象；

步骤五，将每一行分解对象分解得到的各级加总平均值IMF和残差值R重新组合成为新的图像；

步骤六，通过去除第一个加总平均值IMF1来剔除新的图像中的高频噪声；

步骤七，通过图像的自适应噪声的完备经验模态分解CEEMDAN算法，将各加总平均值IMF分量按照满足从高频到低频系列的要求进行不同频率加总平均值IMF层的分布，这样不同的加总平均值IMF反映了不同的频率信息；

步骤八，通过合并某一个或某几个不同频率加总平均值IMF层得到另外的一幅新图像，该新图像即为提取得到的激光散斑图像中的不同频率的信息。

本发明提供的激光散斑图像处理分析方法，还可以具有这样的特征，其特征在于：其中，步骤三中使用图像的自适应噪声的完备经验模态分解CEEMDAN算法进行分解，对于一个信号x[n]，该CEEMDAN算法步骤如下：

A，首先定义算子E_j(·)

E_j(·)＝x[n]+ε_kωⁱ[n](i＝1,…,I) (1)

ωⁱ[n](i＝1,…,I)是单位方差的零均值高斯白噪声的不同实现，ε_k是噪声标准偏差；

B，使用不同的噪声实现通过经验模态分解算法EMD重复分解过程I次，计算加总平均值IMF并将其定义为目标信号x[n]的IMF1(n)，公式为

${\mathrm{IMF}}_1\left(n\right)=\frac{1}{I}\underset{i=1}{\overset{n}{\Σ}}{E}_1\left(x\right(n)+{\mathrm{\ϵ}}_0{\mathrm{\ω}}^i(n\left)\right)---\left(2\right);$

C，对k＝1,计算一级残余值r₁[n],公式为

r₁[n]＝x[n]-IMF₁(n) (3）；

D，经验模态分解算法EMD实现r₁(n)+ε₁E₁(ωⁱ(n)),直到满足第一个加总平均值IMF(n)条件，并定义总体平均值为IMF₂(n),公式为

${\mathrm{IMF}}_2\left(n\right)=\frac{1}{I}\underset{i=1}{\overset{n}{\Σ}}{E}_1\left(x_1\right(n)+{\mathrm{\ϵ}}_1{E}_1({\mathrm{\ω}}^i\left(n\right)\left)\right)---\left(4\right);$

E，对于k＝2，···，K,计算k阶残差rk[n],公式为

rk[n]＝rk-1[n]-IMFk(n)(5)；

F，提取r_k(n)+ε_kE_k(ωⁱ(n))的加总平均值IMF₁(n)分量，计算他们的总体平均值得目标信号的IMF_k+1(n)，公式为

${\mathrm{IMF}}_{k+1}\left(n\right)=\frac{1}{I}\underset{i=1}{\overset{n}{\Σ}}{E}_1(r_k+{\mathrm{\ϵ}}_k{E}_k({\mathrm{\ω}}^i\left)\right)---\left(6\right);$

G，重复步骤(E)(F)，直到残差不能再被分解为止，则得到最终残差R(n)为

$R\[n\]=x\[n\]-\underset{k=1}{\overset{k}{\Σ}}{\mathrm{IMF}}_k\left(n\right)---\left(7\right)$

其中K是加总平均值IMF(n)的总数.因此，目标信号x(n)可以被表达为

$x\[n\]=\underset{k=1}{\overset{k}{\Σ}}{\mathrm{IMF}}_k\left(n\right)+R\[n\]---\left(8\right);$

H，分解结束，最终分解得到各级加总平均值IMF及一个残余值R。

发明作用与效果

根据本发明提供的激光散斑图像处理分析方法，由于使用了自适应噪声的完备经验模态分解(Complete Ensemble Empirical Mode Decomposition With AdaptiveNoise，CEEMDAN)算法，可以较好地解决图像中的噪声干扰，从而获得图像中不同频率的信息。

附图说明

图1为本发明激光散斑图像处理分析方法流程示意图；

图2模拟血管激光散斑图像，a为截取前图像，b为a中黑框部分的放大图像，a、b图中的左边软管的血液流速较右边快；

图3为CEEMDAN算法处理单行激光散斑图像数据结果示意图，图中从上往下ECG分别表示处理前，IMF1至IMF7分别表示1-7级的IMF，R表示残余值R；以及

图4为CEEMDAN算法激光散斑处理结果，图中1级IMF至7级IMF分别1值7级IMF重构出的图像。

具体实施方式

为了使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，以下实施例结合附图对本发明的激光散斑图像处理分析方法的流程和具体操作方法、原理作具体阐述。

图1为本发明激光散斑图像处理分析方法流程示意图。

图2模拟血管激光散斑图像，a为截取前图像，b为a中黑框部分的放大图像，a、b图中的左边软管的血液流速较右边快。

具体操作步骤如附图1所示。

1、使用激光散斑血流成像装置采集散斑图f(x,y)，图像尺寸为m×n，例如图2所示。

2、根据图像尺寸建立多个m×n大小的空矩阵，用于存储分解结果。

3、将图像按行(列)分割，得到m行(列)数据，将每一行(列)数据分别作为分解对象。

4、使用CEEMDAN算法进行分解

A，首先定义算子E_j(·)

E_j(·)＝x[n]+ε_kωⁱ[n](i＝1,…,I) (1)

ωⁱ[n](i＝1,…,I)是单位方差的零均值高斯白噪声的不同实现，ε_k是噪声标准偏差。

B，使用不同的噪声实现通过经验模态分解算法EMD重复分解过程I次，计算加总平均值IMF并将其定义为目标信号x[n]的IMF1(n)，公式为

${\mathrm{IMF}}_1\left(n\right)=\frac{1}{I}\underset{i=1}{\overset{n}{\Σ}}{E}_1\left(x\right(n)+{\mathrm{\ϵ}}_0{\mathrm{\ω}}^i(n\left)\right)---\left(2\right).$

C，对k＝1,计算一级残余值r₁[n],公式为

r₁[n]＝x[n]-IMF₁(n) (3)。

D，经验模态分解算法EMD实现r₁(n)+ε₁E₁(ωⁱ(n)),直到满足第一个加总平均值IMF(n)条件，并定义总体平均值为IMF₂(n),公式为

${\mathrm{IMF}}_2\left(n\right)=\frac{1}{I}\underset{i=1}{\overset{n}{\Σ}}{E}_1\left(x_1\right(n)+{\mathrm{\ϵ}}_1{E}_1({\mathrm{\ω}}^i\left(n\right)\left)\right)---\left(4\right).$

E，对于k＝2，···，K,计算k阶残差r_k[n],公式为

r_k[n]＝r_k-1[n]-IMF_k(n) (5)。

F，提取rk(n)+εkEk(ωi(n))的加总平均值IMF1(n)分量，计算他们的总体平均值得目标信号的IMFk+1(n)，公式为

${\mathrm{IMF}}_{k+1}\left(n\right)=\frac{1}{I}\underset{i=1}{\overset{n}{\Σ}}{E}_1(r_k+{\mathrm{\ϵ}}_k{E}_k({\mathrm{\ω}}^i\left)\right)---\left(6\right).$

G，重复步骤(E)(F)，直到残差不能再被分解为止，则得到最终残差R(n)为

$R\[n\]=x\[n\]-\underset{k=1}{\overset{k}{\Σ}}{\mathrm{IMF}}_k\left(n\right)---\left(7\right)$

其中K是加总平均值IMF(n)的总数.因此，目标信号x(n)可以被表达为

$x\[n\]=\underset{k=1}{\overset{k}{\Σ}}{\mathrm{IMF}}_k\left(n\right)+R\[n\]---\left(8\right).$

H，分解结束，最终分解得到各级加总平均值IMF及一个残余值R。

图3为CEEMDAN算法处理单行激光散斑图像数据结果示意图，图中从上往下ECG分别表示处理前，IMF1至IMF7分别表示1-7级的IMF，R表示残余值R。

5、第一行数据分解完毕，分解结果如图3所示，得到各级IMF及残余值R，随后分解下一行数据，直至遍历整幅图像。

图4为CEEMDAN算法激光散斑处理结果，图中1级IMF至7级IMF分别1值7级IMF重构出的图像。

6、将每一行(列)数据分解得出的各阶IMF和残余值重新组合重构出整幅图片，各级IMF重构出的图像如图4所示，从IMF2后续处理可以分辨出，左边软管区域比右边软管区域频率高，判断左边流速快，该判断也和实际情况吻合。在IMF2以及IMF4中，可以分辨出血管与血管壁具有完全不同的频率信息，据此可以分辨出血管与血管壁。在IMF4中，可以清楚分辨出血管边缘，从而可以分辨出血管在图像中的具体位置。

7、去除掉IMF1后对其余IMF重新组合，可重构成去噪后图像。

8、根据运算结果进行后续处理及结果分析。

实施例的有益作用和效果

根据本实施例提供的激光散斑图像处理分析方法，由于使用了自适应噪声的完备经验模态分解(Complete Ensemble Empirical Mode Decomposition With AdaptiveNoise，CEEMDAN)算法，可以较好地解决图像中的噪声干扰，从而获得图像中不同频率的信息。

本发明相比现有的激光散斑图像处理技术，本发明的技术方案具有以下有益效果：

本发明采用CEEMDAN算法，它是在分解的每一阶段添加一个特定白噪声，并计算一个唯一残差以得到每个IMF，这种方法分解出的各级IMF能够对原始信号精确重构。有效解决了EMD的模态混叠问题，同时也保留了EMD处理非平稳信号的优势，如自适应性、二进滤波特性等。

本发明采用CEEMDAN算法，相比传统计算衬比的方法，可以剔除图像中的高频噪声，使处理结果中图像特征或参数精度高。

本发明采用CEEMDAN算法，可以将不同频率信息的分离，在激光散斑血流成像中将可以将血管与管壁分离，实现某频率范围的信息突出显示。

本发明可以将激光散斑图像不同频率信息的分离，被测物中不同速度的散射颗粒反映在散斑图像中即具有不同的频率，因此，可以使用本发明区分出激光散斑成像被测区域散射颗粒的运动速度信息。

Claims (2)

1.一种激光散斑图像处理分析方法，用于提取出所述激光散斑图像中的不同频率的信息，其特征在于，包括以下步骤：

步骤一，使用激光散斑血流成像装置采集散斑图f(x,y)，该图像的尺寸为m×n；

步骤二，将所述散斑图按行列分割，得到m行×n列数据，将每一行数据分别作为分解对象；

步骤三，使用图像的自适应噪声的完备经验模态分解CEEMDAN算法对第一行所述分解对象进行分解，得到该行的加总平均值并将其定义为目标信号x[n]的IMF1(n)，最终分解得到对应该行的各级加总平均值IMF及一个残余值R；

步骤四，分解完第一行的所述分解对象后，继续分解下一行的所述分解对象，直至分解完整幅所述散斑图的每一行的所述分解对象；

步骤五，将每一行所述分解对象分解得到的各级加总平均值IMF和残差值R重新组合成为新的图像；

步骤六，通过去除第一个加总平均值IMF1来剔除所述新的图像中的高频噪声；

步骤七，通过图像的自适应噪声的完备经验模态分解CEEMDAN算法，将各加总平均值IMF分量按照满足从高频到低频系列的要求进行不同频率加总平均值IMF层的分布，这样不同的加总平均值IMF反映了不同的频率信息；

步骤八，通过合并某一个或某几个所述不同频率加总平均值IMF层得到另外的一幅新图像，该新图像即为提取得到的所述激光散斑图像中的不同频率的信息。

2.根据权利要求1所述的激光散斑图像处理分析方法，其特征在于：

其中，步骤三中使用图像的自适应噪声的完备经验模态分解CEEMDAN算法进行分解，对于一个信号x[n]，该CEEMDAN算法步骤如下：

A，首先定义算子E_j(·)

E_j(·)＝x[n]+ε_kωⁱ[n](i＝1,…,I) (1)

ωⁱ[n](i＝1,…,I)是单位方差的零均值高斯白噪声的不同实现，ε_k是噪声标准偏差；

B，使用不同的噪声实现通过经验模态分解算法EMD重复分解过程I次，计算加总平均值IMF并将其定义为目标信号x[n]的IMF1(n)，公式为

${\mathrm{IMF}}_1\left(n\right)=\frac{1}{I}\underset{i=1}{\overset{n}{\Σ}}{E}_1\left(x\right(n)+{\mathrm{\ϵ}}_0{\mathrm{\ω}}^i(n\left)\right)---\left(2\right);$

C，对k＝1,计算一级残余值r₁[n],公式为

r₁[n]＝x[n]-IMF₁(n) (3)；

D，经验模态分解算法EMD实现r₁(n)+ε₁E₁(ωⁱ(n))，直到满足第一个加总平均值IMF(n)条件，并定义总体平均值为IMF₂(n)，公式为

${\mathrm{IMF}}_2\left(n\right)=\frac{1}{I}\underset{i=1}{\overset{n}{\Σ}}{E}_1\left(x_1\right(n)+{\mathrm{\ϵ}}_1{E}_1({\mathrm{\ω}}^i\left(n\right)\left)\right)---\left(4\right);$

E，对于k＝2，···，K，计算k阶残差r_k[n]，公式为

r_k[n]＝r_k-1[n]-IMF_k(n) (5)；

F，提取r_k(n)+ε_kE_k(ωⁱ(n))的加总平均值IMF₁(n)分量，计算他们的总体平均值得目标信号的IMF_k+1(n)，公式为

${\mathrm{IMF}}_{k+1}\left(n\right)=\frac{1}{I}\underset{i=1}{\overset{n}{\Σ}}{E}_1(r_k+{\mathrm{\ϵ}}_k{E}_k({\mathrm{\ω}}^i\left)\right)---\left(6\right);$

G，重复步骤(E)(F)，直到残差不能再被分解为止，则得到最终残差R(n)为

$R\[n\]=x\[n\]-\underset{k=1}{\overset{k}{\Σ}}{\mathrm{IMF}}_k\left(n\right)---\left(7\right)$

其中K是加总平均值IMF(n)的总数.因此，目标信号x(n)可以被表达为

$x\[n\]=\underset{k=1}{\overset{k}{\Σ}}{\mathrm{IMF}}_k\left(n\right)+R\[n\]---\left(8\right);$

H，分解结束，最终分解得到各级加总平均值IMF及一个残余值R。