CN103767698A - 一种基于近红外光谱技术的静脉动态特性分析装置及其分析方法 - Google Patents

Abstract

本发明涉及一种基于近红外光谱技术的静脉动态特性分析装置及其分析方法，该装置包括静物台、近红外摄像机和位于静物台上方的至少两个反射光源，所述反射光源以静物台的纵向中轴线为轴对称放置，所述静物台正下方设近红外透射光源，所述近红外摄像机置于静物台正上方；所述近红外摄像机与计算机连接，将所检测到的图像传送至计算机中的数据处理单元，通过预设算法实现静脉动态特性分析。本发明将近红外反射和透射光源相结合，合适的光路设计并有效抑制了单独透射光源引起的噪点问题，提高了图像对比度和成像质量，且安全无侵入性，成本较低，实时性好，结合图像增强处理和帧间差分等算法可实现动态特性的准确计算。

Description

一种基于近红外光谱技术的静脉动态特性分析装置及其分析方法

技术领域

本发明属于医疗及科学实验辅助装置，具体涉及的是一种利用近红外光谱对人类及哺乳动物静脉血管动态特性的分析装置及其分析方法。 

背景技术

研究表明，血流参数能很好地反应动物个体的血管功能状态和生理状况，可为诸多疾病的预防、早期诊断和治疗起到了非常积极的作用。尤其当许多疾病初期，个体尚未发现明显症状，而此时血管阻力和血管粘性等血流参数实际已发生变化。因此，及时、便捷、准确地监测血流动态参数对动物健康水平的评估和疫病的预防有重要意义，更是将医学诊断从基于血管形态学转向基于形态血流动力学的特征分析。 

静脉血管再填充速度是静脉动态特性分析中非常重要的参数。目前国内外对于血流速度等大多采用超声多普勒方法或激光多普勒方法，利用超声或激光方法将组织血液图像呈现到屏幕上，根据测量屏幕中血流图像的变化反应血流速度。但多普勒方法一般仪器较为庞大，成本高，且对操作人员专业技能要求高。 

此外还有利用非多普勒原理的血流速度直接测量方法。专利（公开号：CN1957266A）采用光扫描显微镜测量，利用时空分析方法测量红细胞在血流中的移动速度，此方法同样成本高，不利于推广。专利（公开号：CN1358478A）公开了一种表面安装有可探测血液温度变化热敏电阻的导管，利用血液温度变化感知原理检测血流速度。专利（公开号：CN102755157A）则基于容积脉搏波传感器检测人体指端脉搏搏动信号，处理得到系统所需加速脉搏波，提取出脉搏波传导时间，确定血流速度。 

研究表明近红外光可穿透肌肉和骨骼，仅静脉中的还原血红蛋白对近红外光有吸收特性。基于此原理国内外已有诸多近红外手掌静脉静态图像获取的产品及研究，但目前尚未拓展至静脉血管的动态特性分析。 

发明内容

针对上述现有技术中存在的问题，本发明的目的在于利用近红外透射和反射光源结 合的方法，提供一种主动照明实现静脉成像及其动态特性分析的装置。该装置安全可靠，无侵入性，成本较低，实时性好，有助于深入了解静脉血流的动态特性。 

为实现上述目的，本发明所采取的技术方案：一种基于近红外光谱技术的静脉动态特性分析装置，包括静物台、近红外摄像机和位于静物台上方的至少两个反射光源，所述反射光源以静物台的纵向中轴线为轴对称放置，所述静物台正下方设近红外透射光源，所述近红外摄像机置于静物台正上方；所述近红外摄像机与计算机连接，将所检测到的图像传送至计算机中的数据处理单元，通过预设算法实现静脉动态特性分析。 

所述近红外摄像机前端放置有近红外滤光片，所述反射光源前端放置有漫射板，所述近红外透射光源前端放置有透镜。 

所述近红外透射光源采用波段在800nm-910nm之间的LED阵列近红外灯，所述近红外滤光片（6）采用波段在800nm-910nm之间的窄带吸收滤光片。 

所述静物台可使用聚碳酸酯（PC）或聚甲基丙烯酸甲酯（PMMA）等材料，PC或PMMA材料不易破碎，对近红外光穿透性好（可达90%以上），并能滤除可见光等杂散光。 

本发明的另一目的是提供一中基于近红外光谱的静脉动态特性分析方法，该方法包括： 

对待测部位静脉人工施加外力后迅速撤去，利用图像采集单元获取整个血管血流再填充过程的图像，传输至计算机中的静脉动态特性分析单元。 

所述静脉动态特性分析单元包含如下步骤：（1）选取每帧图像中的待检测静脉区域，即感兴趣区域选取；（2）图像增强处理；（3）动态参数计算；其中：感兴趣区域（英文缩写为ROI）可在整幅图像中手动选取，ROI的选取可有效缩短处理时间，提高精度。 

所述图像增强处理过程具体包括灰度归一化、双边滤波消噪和直方图均衡化处理，具体包括： 

（1）假设原始图像的灰度值为I(x，y)，则经过归一化处理后图像的灰度值为 

$I_1(x,y)=\left\{\begin{array}{cc}{M}_1+\sqrt{(V_1{(I(x,y)-M)}^2/V)}& (I(x,y)>M)\\ M_1-\sqrt{(V_1{(I(x,y)-M)}^2/V)}& (I(x,y)>M)\end{array}\right.$

其中M和V分别是原始图像灰度值的均值和方差值。参考前期实验结果，优选地将 变化后图像的灰度均值M₁设定为110到150之间，方差值V₁设定为255。 

（2）对灰度值归一化后的图像进行双边滤波消噪处理（Bilateral Filtering,BF），以消除图像增强处理中引入的噪声干扰，同时可解决边界的模糊问题。假设像素点p和q的灰度值分别为I(p)和I(q)，则变化方程如下： 

$I_{\mathrm{BF}}\left(p\right)=\frac{1}{W\left(p\right)}{\mathrm{\Σ}}_{q\∈S}{G}_{{\mathrm{\σ}}_d}\left(\right|\left|(p-q)\right|\left|\right)G_{{\mathrm{\σ}}_r}(I\left(p\right)-I\left(q\right))I\left(q\right),$

其中， $W\left(p\right)={\mathrm{\Σ}}_{q\∈S}{G}_{{\mathrm{\σ}}_d}\left(\right|\left|(p-q)\right|\left|\right)G_{{\mathrm{\σ}}_r}(I\left(p\right)-I\left(q\right)),$ 双边滤波器的空间临近函数和灰度相似度函数均取参数为欧氏距离的高斯函数，即：

σ_d和σ_r为基于高斯函数的标准差。 

（3）参考实验结果，优选地选取对比度受限自适应直方图均衡化处理（Contrast Limited Adaptive Histogram Equalization,CLAHE），此方法可在有效扩大图像灰度值的动态范围从而增强图像的同时，通过限制局部直方图的高度抑制局部对比度的增强及噪声方法，具体算法如下： 

将图像划分为多个M×N（优选地可设置为10*8个）连续不重叠的子区域，根据预编程序统计得到每个区域的灰度直方图Hist(i),通过以下公式计算直方图阈值H_max $H_{\max }=S_{\max }\frac{\mathrm{MN}}{\mathrm{Bins}}$

其中，Bins为有效的灰度级数量，选择较小的值（优选地可设置为128）可以加速处理速度产生较好的输出图像质量。S_max为剪切极限系数，其意义是限定的变换函数的最大斜率，手动设置其大小（优选地可设置为4.5）。 

对每个子区域进行灰度直方图剪切，将直方图中高于H_max的部分平均分配到直方图的底部，得到Hist′(i)。 

${\mathrm{Hist}}^{\&prime;}\left(i\right)=\left\{\begin{array}{cc}\mathrm{Hist}\left(i\right)& \mathrm{Hist}\left(i\right)<H_{\max }\\ H_{\max }& \mathrm{Hist}\left(i\right)\&GreaterEqual;H_{\max }\end{array}\right.$

计算子区域的累积分布函数

p=0,1……,L-1，即图像的灰度范围为[0,L-1]，从而得到变换函数T(i)： 

$T\left(i\right)=\frac{\mathrm{Bins}}{\mathrm{MN}}{H}_c\left(i\right)$

通过变化函数T(i)，计算出每个子区域中心像素的灰度值。为了提高计算速度，采 用双线性插值法计算图像内部其他像素点的灰度值。外围中心像素点之外的像素灰度值用其对应的两个临近中心像素线性插值计算，图像的四角角点处的像素的灰度值用其所在子区域的变换函数计算。 

所述动态参数计算过程包括：（1）对增强处理后的图像进行帧间差分；（2）计算血液流过的距离；（3）计算静脉血液瞬时流速和单位长度的静脉血管恢复时间。 

结合上述方法，本发明提供一套由近红外光源、漫射板、聚焦透镜、静物台、滤光片和近红外摄像机组成的图像采集单元。该装置的特征在于其结合反射和透射两种成像方式，即有两个（或以上）反射光源位于静物台上方，以静物台的垂直平分线为轴呈对称放置，一个透射光源位于静物台正下方。反射光源前端放置有漫射板，可使得光源发射的光均匀漫射；透射光源前端放置有透镜，可对近红外光起聚焦作用。近红外摄像机置于静物台正上方，在摄像头前端放置有滤光片，用于滤除可见光，提高摄像机图像成像质量。近红外摄像机与计算机连接，将所检测到的图像传送至计算机中的数据处理单元，实现静脉动态特性分析。 

相对比现有技术，本发明的主要创新点有： 

（1）将近红外光谱方法引入静脉的动态特性分析中，该方法安全无侵入性，成本较低，实时性好，结合图像增强处理和帧间差分等算法可实现动态特性的准确计算。 

（2）将近红外反射和透射光源相结合，合适的光路设计并有效抑制了单独透射光源引起的噪点问题，提高了图像对比度和成像质量。 

附图说明：

以下结合具体实施方式给出本发明的装置和算法附图，并将作进一步详细描述。 

图1：本发明装置的结构示意框图；其中，1：近红外反射光源，2：近红外透射光源，3：漫射板，4：透镜，5：静物台，6：近红外滤光片，7：近红外摄像机，8：计算机。 

图2：本发明的算法流程框图。 

图3：图像中子区域中心像素点和其他像素点示意图。 

具体实施方式：

下面结合具体实施方式并参照附图对本发明作进一步详细说明。 

如图1所示，两个反射光源位于静物台上方，呈对称性放置，透射光源位于静物台正下方，即置于对称轴中心位置。优选地，光源为850nm LED阵列近红外灯，反射光源前端放置有漫射板，透射光源前端放置有透镜。近红外摄像机置于静物台正上方，近红外摄像机前端放置有850nm窄带吸收滤光片，将所检测到的图像传送至计算机中的数据处理单元，通过预设算法实现静脉动态特性分析。 

具体分析步骤如下： 

首先打开光源预热20-30分钟，打开摄像机开始图像采集，将待测部位放置于静物台，对待测部位静脉人工施加外力后迅速撤去，获取整个血管血流再填充过程的每帧图像，传输至计算机中的数据处理单元。 

在计算机中进行静脉动态特性分析，手动选取每帧图像中的感兴趣区域ROI（即待检测静脉区域），对选取的ROI进行按如下步骤进行图像增强处理，具体操作如下： 

（1）为减少外界因素变化（光强等环境因素、实验动物肤色，表层肤质等个体差异）干扰，首先对图像进行灰度归一化处理，原始图像的灰度值为I(x，y)，则经过归一化处理后图像的灰度值为 

$I_1(x,y)=\left\{\begin{array}{cc}{M}_1+\sqrt{(V_1{(I(x,y)-M)}^2/V)}& (I(x,y)>M)\\ M_1-\sqrt{(V_1{(I(x,y)-M)}^2/V)}& (I(x,y)>M)\end{array}\right.$

其中M和V分别是原始图像灰度值的均值和方差值。参考前期实验结果，优选地将变化后图像的灰度均值M₁设定为110到150之间，方差值V₁设定为255。 

（2）对灰度值归一化后的图像进行双边滤波消噪处理（Bilateral Filtering,BF），以消除图像增强处理中引入的噪声干扰，同时可解决边界的模糊问题。双边滤波在考虑空间邻近性的同时也考虑亮度的相似性，每个像素被相邻像素的加权平均取代，假设像素点p和q的灰度值分别为I(P)和I(q)，则变化方程如下： 

$I_{\mathrm{BF}}\left(p\right)=\frac{1}{W\left(p\right)}{\mathrm{\&Sigma;}}_{q\&Element;S}{G}_{{\mathrm{\&sigma;}}_d}\left(\right|\left|(p-q)\right|\left|\right)G_{{\mathrm{\&sigma;}}_r}(I\left(p\right)-I\left(q\right))I\left(q\right),$

其中， $W\left(p\right)={\mathrm{\&Sigma;}}_{q\&Element;S}{G}_{{\mathrm{\&sigma;}}_d}\left(\right|\left|(p-q)\right|\left|\right)G_{{\mathrm{\&sigma;}}_r}(I\left(p\right)-I\left(q\right)),$ 双边滤波器的空间临近函数和灰度相似度函数均取参数为欧氏距离的高斯函数，即：

I(p)和I(q)分别为像素点p和q的灰度值，σ_d和σ_r为基于高斯函数的标准差。 

（3）最后进行直方图均衡法处理，参考实验结果，优选地选取对比度受限自适应直方图均衡化处理（Contrast Limited Adaptive Histogram Equalization,CLAHE），此方法可在有效扩大图像灰度值的动态范围从而增强图像的同时，通过限制局部直方图的高度抑制局部对比度的增强及噪声方法，具体算法如下： 

将图像划分为多个M×N（优选地，10*8个）连续不重叠的子区域，根据预编程序统计得到每个区域的灰度直方图Hist(i),通过以下公式计算直方图阈值H_maM

$H_{\max }=S_{\max }\frac{\mathrm{MN}}{\mathrm{Bins}}$

其中，Bins为有效的灰度级数量，选择较小的值（优选地，128）可以加速处理速度产生较好的输出图像质量。S_max为剪切极限系数，其意义是限定的变换函数的最大斜率，手动设置其大小（优选地，4.5）。 

对每个子区域进行灰度直方图剪切，将直方图中高于H_max的部分平均分配到直方图的底部，得到Hisr′(i)。 

${\mathrm{Hist}}^{\&prime;}\left(i\right)=\left\{\begin{array}{cc}\mathrm{Hist}\left(i\right)& \mathrm{Hist}\left(i\right)<H_{\max }\\ H_{\max }& \mathrm{Hist}\left(i\right)\&GreaterEqual;H_{\max }\end{array}\right.$

计算子区域的累积分布函数

p=0,1……,L-1,即图像的灰度范围为[0,L-1]，从而得到变换函数T(i)： 

$T\left(i\right)=\frac{\mathrm{Bins}}{\mathrm{MN}}{H}_c\left(i\right)$

通过变化函数T(i)，计算出每个子区域中心像素的灰度值。为了提高计算速度，采用双线性插值法计算图像内部其他像素点的灰度值。如图3所示，“■”为个子区域中心像素，“●”为待计算的像素点，计算方程如下： 

M(i)=p[qI_--(i)+(1-q)I_+-(i)]+(1-p)[qI_-+(i)+(1-q)I₊₊(i)] 

其中，M(i)为像素点(x,y)的灰度值，其四周临近的中心像素点的灰度值分别用I_--(i)、I_+-(i)、I_-+(i)、I₊₊(i)表示。外围中心像素点之外的像素灰度值用其对应的两个临近中心像素线性插值计算，图像的四角角点处的像素的灰度值用其所在子区域的变换函数计算。 

最后，将经过上述增强处理后的图像进行帧间差分，在差分图像中以逐行扫描的方式，根据所设定阈值寻找在同一连通域中的血管首尾坐标值（分别为x_h(t)、y_h(t)和 x_t(t)、y_t(t))，则可得血液流过的距离： $d\left(t\right)=\sqrt{{(x_h\left(t\right)-x_t\left(t\right))}^2+{(x_h\left(t\right)-x_t\left(t\right))}^2}.$ 所设阈值的目的在于消除查分图像中其他动态信息（如噪声）干扰，确保所找出流动距离的正确性。 

静脉血液瞬时流速可由每帧血液流过的距离(d)除以帧间时间差(T)计算得到：v(t)＝d(t)/T。计算从血液开始流动到血管完全恢复的帧数差(n)可得到血管再填充时间T_recovery＝n·T，再除以所测试血管的长度(L)，可得到单位长度的静脉血管恢复时间,即T_normalixed＝T_recovery/L。 

Claims (8)

1.一种基于近红外光谱技术的静脉动态特性分析装置，其特征在于：包括静物台（5）、近红外摄像机（7）和位于静物台上方的至少两个反射光源（1），所述反射光源（1）以静物台的纵向中轴线为轴对称放置，所述静物台（5）正下方设近红外透射光源，所述近红外摄像机（7）置于静物台正上方；所述近红外摄像机（7）与计算机（8）连接，将所检测到的图像传送至计算机中的数据处理单元，通过预设算法实现静脉动态特性分析。 

2.根据权利要求1所述的基于近红外光谱技术的静脉动态特性分析装置，其特征在于：所述近红外摄像机（7）前端放置有近红外滤光片（6），所述反射光源（1）前端放置有漫射板（3），所述近红外透射光源（2）前端放置有透镜（4）。 

3.根据权利要求2所述的基于近红外光谱技术的静脉动态特性分析装置，其特征在于：所述近红外透射光源采用波段在800-910之间的LED阵列近红外灯，所述近红外滤光片（6）采用波段在800-910之间的窄带吸收滤光片。 

4.根据权利要求1所述的基于近红外光谱技术的静脉动态特性分析装置，其特征在于：所述静物台采用聚碳酸酯或聚甲基丙烯酸甲酯材料。 

5.根据权利要求1-4任一项所述的基于近红外光谱技术的静脉动态特性分析装置的分析方法，其特征在于：步骤如下：（1）选取每帧图像中的待检测静脉区域，（2）图像增强处理，（3）动态参数计算。 

6.根据权利要求5所述的基于近红外光谱技术的静脉动态特性分析装置的分析方法，其特征在于：所述图像增强处理包括： 

（a）灰度归一化：假设原始图像的灰度值为I(x，y)，则经过归一化处理后图像的灰度值为 

其中M和V分别是原始图像灰度值的均值和方差值；将变化后图像的灰度均值M₁设定为110到150之间，方差值V₁设定为240到255之间； 

（b）双边滤波消噪：对灰度值归一化后的图像进行双边滤波消噪处理，假设像素点p和q的灰度值分别为I(p)和I(q)，则变化方程如下： 

其中，

双边滤波器的空间临近函数和灰度相似度函数均取参数为欧氏距离的高斯函数，即：

σ_d和σ_r为基于高斯函数的标准差； 

（c）直方图均衡化处理：算法如下： 

将图像划分为若干个M×N，优选地可设置为10*8个连续不重叠的子区域，根据预编程序统计得到每个区域的灰度直方图Hist(i),通过以下公式计算直方图阈值H_max

其中，Bins为有效的灰度级数量，设置灰度级数量为125到130之间；S_max为剪切极限系数，限定的变换函数的最大斜率，手动设置其大小，设置为4.5； 

对每个子区域进行灰度直方图剪切，将直方图中高于H_max的部分平均分配到直方图的底部，得到Hist'(i); 

计算子区域的累积分布函数

p=0,1……,L-1，即图像的灰度范围为[0,L-1]，从而得到变换函数T(i)： 

通过变化函数T(i)，计算出每个子区域中心像素的灰度值。 

7.根据权利要求6所述的基于近红外光谱技术的静脉动态特性分析装置的分析方法，其特征在于：所述每个子区域中心像素点以及图像四角角点处的像素点的灰度值仍采用此区域的变换函数计算；中心像素之外的像素点采用与此像素点临近的四个区域的变换函数双线性差值获取；外围中心像素点之外的像素灰度值用与此像素点对应的两个临近中心像素线性插值计算。 

8.根据权利要求5所述的基于近红外光谱技术的静脉动态特性分析装置的分析方法，其特征在于：所述动态参数计算过程包括：（a）对增强处理后的图像进行帧间差分；（b）计算血液流过的距离；（c）计算静脉血液瞬时流速和单位长度的静脉血管恢复时间。 

Images