CN110826431B - 基于蒙特卡洛的可见光静脉显像方法 - Google Patents

Abstract

本发明公开了种基于蒙特卡洛的可见光静脉显像方法，属于信息感知与识别技术领域，本发明的方法为1）、建立皮肤模型，对皮肤进行光学模型分析；2）、基于Monte Carlo方法进行皮肤成像模拟；3）、皮肤组织成像，最终实现静脉的显像；本发明利用Monte Carlo方法模拟皮肤图像颜色值，再以该颜色值为输入，皮肤参数值为输出，训练神经网络。对于给定的皮肤图像，从获取的皮肤参数值分布图中即可实现静脉显像，解决了现有技术中存在的问题，在针对数码成像的可见光图片进行静脉显像，为后续的识别奠定重要基础。

Description

基于蒙特卡洛的可见光静脉显像方法

技术领域

本发明属于一种信息感知与识别技术领域，具体是一种基于蒙特卡洛的可见光静脉显像方法。

背景技术

近年来，身份识别与个人利益、社会稳定和国家安全息息相关。身份识别技术有基于特定物品、特定知识和生物特征三个类别。前两类如身份证、IC卡、钥匙、密码、条纹码等，往往会存在复制、遗忘、丢失、盗用等隐患，这将会给个人、社会甚至国家造成巨大的损失。与前两类身份识别技术相比，第三类提取人的生物特征信息，主动对个人身份进行认证，早年运用于司法鉴定。指纹和人脸识别技术日趋成熟，但是在司法领域往往很难获取罪犯的脸部及指纹信息，这给识别和鉴定带来了极大困难。

静脉识别，作为一种新型的生物特征识别技术，提取纹络中的特征信息来实现身份识别，识别率高，安全性高，使用便利，因而具备不可比拟的优越性。但是，传统的静脉识别仪器依赖近红外成像设备，维护困难且研究成本很高。目前已有的一些可见光图像静脉显像方法存在着模型简单、鲁棒性差的问题。

发明内容

本发明针对现有技术中存在的鲁棒性差等的问题，公开了一种基于蒙特卡洛的可见光静脉显像方法，利用MC方法模拟皮肤图像颜色值，再以该颜色值为输入，皮肤参数值为输出，训练神经网络。对于给定的皮肤图像，从获取的皮肤参数值分布图中即可实现静脉显像，解决了现有技术中存在的问题。

本发明是这样实现的：

一种基于蒙特卡洛的可见光静脉显像方法，其方法步骤如下：

步骤一、建立皮肤模型，对皮肤进行光学模型分析；生物特征具有不易遗忘、防伪性能好、不易被盗、随身“携带”和随时随地可用等优点，因此利用人自身的生理和行为特征进行身份识别的生物特征识别技术可以克服传统认证方法的缺陷，本发明建立七层皮肤模型进行分析。

步骤二、基于Monte Carlo方法进行皮肤成像模拟；

步骤三、皮肤组织成像，最终实现静脉的显像。

进一步，所述的步骤一中建立角质层、生发层、真皮层、皮下脂肪层的七层皮肤模型，所述的真皮层依次包括真皮乳头层、上层血管网真皮层、网状真皮层、深部血管网真皮层，对七层皮肤进行光学模型分析；在利用MonteCarlo方法进行皮肤组织模拟时，需要输入五个光学特性参数值：吸收系数μ_α(cm^-1)、散射系数μ_s(cm^-1)、各向异性因子g、折射率n和皮肤各层厚度d(cm)。

进一步，所述的步骤一具体如下：

1.1、分析光在皮肤组织中的吸收作用，七层皮肤的吸收系数如下：

角质层：角质层中对光的吸收主要是水和组织的固有吸收，公式如下：

其中，

是除了水和血液之外的其它组织的固有吸收系数，

是水的体积分数，

是水的吸收系数；

生发层：生发层中的吸收体主要为黑色素，其吸收系数表达式为：

其中，

是黑色素的吸收系数，C_mel是黑色素的体积分数，其取值范围是1.3％～43％；

真皮层：真皮层包含真皮乳头层、上层血管网真皮层、网状真皮层、深部血管网真皮层四层，根据血氧饱和度S和血液中血红蛋白总的体积分数γ，该层的吸收系数可由下式表达：

其中，

是水的体积分数；C_blood是血液的体积分数；γ＝F_HbF_RBCH_t，F_Hb是红细胞中血红蛋白的体积分数，F_RBC是红细胞在所有血细胞中的体积分数，H_t是红细胞比容；

是脱氧血红蛋白的吸收系数，

是有氧血红蛋白的吸收系数，表达式如下：

其中c^hb是血液中血红蛋白的浓度，值为150g/L，m^hb是血红蛋白的分子量，值为66500g/mole，ε^ohb(λ)和ε^dhb(λ)分别为有氧血红蛋白和脱氧血红蛋白的消光系数；

皮下脂肪层：皮下脂肪层的吸收系数可从现有技术中已有数据中获得。分析七层皮肤的光的散射系数μ_s，光的散射系数μ_s是指每单位长度光被散射的概率，单位是cm^-1，七层皮肤的各层散射系数如下：

角质层的散射系数：

其中，λ是波长；

生发层的散射系数：

真皮层的散射系数：

皮下脂肪层的散射系数：

1.2、分析七层皮肤的折射率：

角质层的折射率：

n^str＝1.55   (13)

生发层的折射率：

n^liv＝1.684-1.8723×10⁴λ^-2+1.0964×10¹⁰λ^-4-8.6484×10¹⁴λ^-6   (14)

真皮层的折射率：

n^dermis＝1.309-4.346×10²λ^-2+1.6065×10⁹λ^-4-1.2811×10¹⁴λ^-6   (15)

皮下脂肪层的直射率：

n^fat＝1.44   (16)。

进一步，分析皮肤各层厚度的参数影响，皮肤各层的厚度变化成正相关，即同时变厚，或者同时变薄；以上层血管网真皮层的厚度(cm)为基准，构造线性映射函数如下：

其中，d^pap是真皮乳头层的厚度；d^upp是上层血管网真皮层的厚度；d^ret是网状真皮层的厚度；d^dee是深部血管网真皮层的厚度；d^fat是皮下脂肪层的厚度。

进一步，所述的步骤二中的MC方法模拟时，发射一个光子，当光子走完这一随机步长时，与生物组织发生吸收、散射作用，导致自身的权重衰减、传播方向改变；行走到组织表面时，从上表面逃逸变成漫反射光，从下表面逃逸变成透射光；当该光子自身的权值小于给定的阈值时，光子被组织吸收；最后，对所有光子重复上述过程，实现了整个光束的MC模拟。

进一步，所述的步骤二具体如下：

2.1，发射光子：

入射光束垂直组织平面xoy发射，方向沿着z轴方向，初始位置为原点(0,0,0)，则方向余弦的初始值为(0,0,1)；令光子的初始权值为1，经过反射作用，光子的权值减小至M＝1-R_sp；其中，R_sp是反射系数，R_sp＝(n₀-n₁)²/(n₀+n₁)²，n₀是组织外部介质的折射率，n₁是组织表面的折射率；

2.2，光子的步长：

设区间(a,b)上随机变量η的归一化概率密度函数为p(η)：

设随机变量ξ在区间(0,1)上服从均匀分布，并且η和ξ之间一一映射，即η＝f(ξ)；η＝f(ξ)是增函数时，经过推导一定存在ξ₁∈(0,1)，满足下式：

η＝f(ξ)是减函数时，经过推导也一定存在ξ₁∈(0,1)，满足下式：

对变量η进行抽样，可以得到随机变量η和随机变量ξ之间的关系式η＝f(ξ)；变量ξ通过随机数产生器获取；

运用式(20)对随机变量进行抽样：光子的步长s是随机的，s的计算基于其概率分布的采样，依据互相作用系数μ_t的定义(μ_t＝μ_α+μ_s)，区间(s',s'+ds')中每单位长度的光子-组织互相作用概率为：

或者写成下式：

d(ln(P{s≥s'}))＝-μ_tds'   (22)

将式(22)在(0,s₁)范围内对s'积分，得到指数分布，其中P(s≥0)＝1：

经过重新排列得到步长s的累积分布函数：

则步长s的概率密度函数为：

将式(25)代入到式(19)得：

s₁＝-ln(1-ξ)/μ_t   (26)

ξ的对称性为0.5，用ξ替代1-ξ得到光子的步长：

s₁＝-ln(ξ)/μ_t   (27)

2.3，光子的传输：

光子当前位置用坐标(x,y,z)表示，方向余弦(u_x,u_y,u_z)决定其传输方向；光子到下一个位置时，对坐标进行更新为(x',y',z')：

x'＝x+s·u_x

y'＝y+s·u_y   (28)

z'＝z+s·u_z

2.4，光子的吸收：

光子在生物组织中传输时会与内部粒子发生吸收作用，导致光子的权值衰减，衰减大小为Δw＝w·μ_α/μ_t，则当前光子的权值变为：

w_left＝w-Δw＝w-w·μ_α/μ_t＝w·μ_s/μ_t   (29)

设发射光子的初始权值为w₀，光子与组织内部粒子发生n次吸收作用后，光子的权值变为w_n：

2.5、光子的散射：

光子在发生散射时，其传播方向会跟着改变，进而方向余弦一同改变；θ和ψ一旦确定，光子的方向余弦更新为：

如果光子的方向靠近z轴，|u_z|＞0.99999，那么式(31)修改为下式：

u′_x＝sinθcosψ

u′_y＝sinθsinψ   (32)

u′_z＝sign(u_z)cosθ

用符号函数sign(·)来判断光子z轴的方向。

2.6、光子在边界上的反射和透射：

组织边界分为上、下表面，当光子抵达组织边界时，会产生逃逸出组织和反射回组织两种情况，设定光子的步长s和约化步长s₁，约化步长s₁指光子顺着当前方向与组织边界的间隔，s₁由下式表示：

在直角坐标系中，z₀和z₁分别表示组织上、下表面的z轴坐标；若光子移动的步长s＜s₁，光子处于组织内部；若光子移动的步长s≥s₁，光子将会达到边界，产生全反射或折射；到达组织边界时入射角为α_i＝arccos(|u_z|)，若产生折射，记折射角为α_t，由折射定律可得：

n_i sinα_i＝n_t sinα_t   (34)

若产生全反射，则反射率R(α_i)可由菲涅耳公式得：

若随机数ξ＞R(α_i)，光子产生折射，逃逸出组织，上表面逃逸变成漫反射光，下表面逃逸变成透射光；反之若ξ≤R(α_i)，光子产生全反射，因此方向余弦要更新为(u_x,u_y,u_z)，步长更新为s-s₁；若光子步长不大，在行走完随机步长时光子便会与组织产生吸收和散射作用；

2.7、组织交界处的反射和透射：

假设光子在第i层中传输，光学特性参数为第i层的吸收系数μ_αi、第i层的散射系数μ_si和第i层的折射率n_i，步长s的方向余弦为(u_x,u_y,u_z)；若s＞s₁，光子会行走到第i层与第i-1层或第j＝i+1层，光学特性参数为第j层的吸收系数μ_αj、第j层的散射系数μ_sj和第j层的折射率n_j的交界处，此时由式(22)判断光子的情况；若产生全反射，光子继续在第i层传输，步长更新为s-s₁，

方向余弦更新为(u_x,u_y,u_z)；若产生折射，光子传输到第j＝i+1层。

2.8：光子的终止：

光子的终止情况有：1)光子的权值衰减为0；2)光子逃逸；对于光子的权值衰减为0，设置一个阈值w_ε，当光子的权值w≤w_ε时，光子达到终止条件。

进一步，所述的步骤2.8中对于光子的终止的光子的权值衰减为0的情况时，还可以给定参数M，产生一个服从均匀分布的随机数ξ；当ξ＞1/M，光子的权值变为0，终止；当ξ≤1/M，光之权值变为Mw，继续传输。

进一步，所述的步骤三具体为：

3.1、相机拍摄的可见光范围内的图像，其RGB颜色值由光源、皮肤漫反射率和相机光谱响应函数共同决定；即：

其中，R,G,B表示图像的颜色值，I(λ)是光源，λ是波长，取值范围是400～800nm；R(λ)是皮肤的漫反射率，S(λ)是相机的光谱响应函数，可由相机模型得到；将I(λ)、R(λ)、S(λ)的数据导入，生成RGB图像，得到颜色值数据；

3.2、使用一个三层前馈神经网络模型来训练皮肤图像RGB值到皮肤生理参数C_mel，C_blood，D^dermis之间的映射关系；输入层有3个参数，以MC模拟得到的RGB为输入；以生理参数C_mel,C_blood,D^dermis为输出；隐含层有5个神经元；

3.3、隐含层和输出层使用的激活函数分别为对称型s函数和线性函数，利用BP算法进行训练；对于一张待测的皮肤图像，输入其RGB值，便可从输出的映射图像中实现静脉的显像。

本发明与现有技术相比的有益效果在于：

1、本发明利用静脉识别作为一种新型的生物特征识别技术，提取纹络中的特征信息来实现身份识别，识别率高，安全性高，使用便利，因而具备不可比拟的优越性；由于本发明静脉位于皮肤以下，因此静脉显像是识别中极为关键的一步；传统的静脉识别必须依赖近红外成像设备，成本代价高且使用不便；相对于超声波和红外设备，数码成像技术已经普遍地应用于日常生活当中，目前为大众广泛使用的智能手机都具有数码相机的功能；本发明旨在针对数码成像的可见光图片进行静脉显像，为后续的识别奠定重要基础；

2、本发明采用的蒙特卡洛(Monte Carlo，MC)方法能够有效处理光在生物组织中的传输问题，神经网络亦是近年来的研究热点，将两者结合应用于皮肤成像的模拟与静脉显像，具有极高的研究价；本发明的MC方法可以模拟出不同皮肤参数、不同光照条件和相机模型下的皮肤图像颜色数值，从而获得大量准确有效的训练数据，神经网络具备自主学习特征的能力，其中卷积神经网络经常用来对二维图像数据进行特征提取，在图像处理过程中可以保留像素间的空间位置关系；本发明利用MC方法模拟皮肤图像颜色值，再以该颜色值为输入，皮肤参数值为输出，训练神经网络；

3、本发明对于给定的皮肤图像，从获取的皮肤参数值分布图中即可实现静脉显像，该发明成果对司法、医疗领域具有重要的理论价值和现实意义。

附图说明

图1为本发明建立的七层皮肤模型示意图；

图2为本发明实施例中光子的步长s和约化步长s₁的示意图；

图3为本发明实施例中组织交界处的反射与透射的示意图；

图4为本发明MC方法模拟光子在皮肤中传输的流程图；

图5为本发明实施例中静脉显像结果图；

图6为本发明实施例中静脉显像及其提取结果对比图；

图7为本发明实施例中JAI同步相机拍摄图像的显像结果图；

图8为本发明实施例中不同相机不同身体部位的显像效果。

具体实施方式

为使本发明的目的、技术方案及效果更加清楚，明确，以下列举实例对本发明进一步详细说明；应当指出此处所描述的具体实施仅用以解释本发明，并不用于限定本发明。

步骤一、皮肤光学模型分析：

皮肤是一种非常复杂的结构，在所记载的文献中，对皮肤的研究多为三层或者五层，本发明的研究将建立七层皮肤模型，如图1所示。七层皮肤模型包括：表皮层(角质层、生发层)、真皮层(真皮乳头层、上层血管网真皮、网状真皮层、深部血管网真皮)、皮下脂肪层。

在利用MC方法进行皮肤组织模拟时，需要输入五个光学特性参数值：吸收系数μ_α(cm^-1)、散射系数μ_s(cm^-1)、各向异性因子g、折射率n和皮肤各层厚度d(cm)。

光在皮肤组织中的吸收作用主要由血红蛋白、黑色素等发色团引起。七层皮肤的吸收系数如下：

1)角质层：角质层中对光的吸收主要是水和组织的固有吸收，公式如下^[19]：

其中，

是除了水和血液之外的其它组织的固有吸收系数，

是水的体积分数，

是水的吸收系数。

2)生发层：生发层中的吸收体主要为黑色素，其吸收系数表达式为：

其中，

是黑色素的吸收系数，C_mel是黑色素的体积分数，其取值范围是1.3％～43％。

3)真皮层：真皮层包含四层：真皮乳头层、上层血管网真皮、网状真皮层、深部血管网真皮。根据血氧饱和度S和血液中血红蛋白总的体积分数γ，该层的吸收系数可由下式表达：

其中，

是水的体积分数；C_blood是血液的体积分数；γ＝F_HbF_RBCH_t，F_Hb是红细胞中血红蛋白的体积分数，F_RBC是红细胞在所有血细胞中的体积分数，H_t是红细胞比容。

是脱氧血红蛋白的吸收系数，

是有氧血红蛋白的吸收系数，表达式如下：

其中c^hb是血液中血红蛋白的浓度，值为150g/L，m^hb是血红蛋白的分子量，值为66500g/mole，ε^ohb(λ)和ε^dhb(λ)分别为有氧血红蛋白和脱氧血红蛋白的消光系数。

4)皮下脂肪层：皮下脂肪层的吸收系数如下：

光的散射系数μ_s是指每单位长度光被散射的概率，单位是cm^-1。在本发明中，七层皮肤的各层散射系数如下：

1)角质层的散射系数：

其中，λ是波长。

2)生发层的散射系数：

3)真皮层的散射系数：

4)皮下脂肪层的散射系数：

在可见光和近红外光光谱范围内，各向异性因子g的取值范围是0.8～0.95。在本发明中，七层皮肤的各层散射系数如下：

1)角质层的g：

g^str＝0.198+0.304×[(1-exp(-(λ-507.4)/2404)]   (9)

2)生发层的g：

g^liv＝0.754+0.546×[(1-exp(-(λ-500)/1806)]   (10)

3)真皮层的g：

g^dermis＝0.715+0.00038×[(1-exp(-(λ-542)/1129)]   (11)

4)皮下脂肪层的g：

g^fat＝0.75   (12)

光子在组织的边界或者层与层交界处，组织的折射率对光路的测量有着极大的影响。因此，各层皮肤组织的折射率有十分重要的地位：

1)角质层的折射率：

n^str＝1.55   (13)

2)生发层的折射率：

n^liv＝1.684-1.8723×10⁴λ^-2+1.0964×10¹⁰λ^-4-8.6484×10¹⁴λ^-6   (14)

3)真皮层的折射率：

n^dermis＝1.309-4.346×10²λ^-2+1.6065×10⁹λ^-4-1.2811×10¹⁴λ^-6   (15)

4)皮下脂肪层的直射率：

n^fat＝1.44   (16)

皮肤各层的厚度直接影响着光的漫反射率和透射率的计算，因此皮肤厚度是MC模拟方法中的一个重要的参数。角质层主要由角质细胞构成，厚度约为20μm；生发层中主要是生发细胞，厚度约为80～100μm；真皮层有四层：真皮乳头层厚度约为150～200μm；上层血管网真皮厚度约为80～100μm；网状真皮层的厚度约为1400～1600μm；深部血管网真皮层的厚度约为80～120μm；皮下脂肪层主要为脂肪细胞，厚度约为1000～6000μm。

一般来说，皮肤各层的厚度变化成正相关，即同时变厚，或者同时变薄。本发明在波长为400～800nm的光谱范围内进行MC方法模拟，为了符合变化规律，本发明以上层血管网真皮层的厚度(cm)为基准，构造线性映射函数如下：

其中，d^pap是真皮乳头层的厚度；d^upp是上层血管网真皮层的厚度；d^ret是网状真皮层的厚度；d^dee是深部血管网真皮层的厚度；d^fat是皮下脂肪层的厚度。

步骤二、基于Monte Carlo(MC)方法进行皮肤成像模拟：

MC方法存在以下四点优势：①关于生物组织的光源和边界没有要求；②无需求解传输方程；③关于生物组织的光学特性参数没有要求；④计算机模拟工作，因而该方法被广泛应用。MC方法模拟时，发射一个光子，当光子走完这一随机步长时，会与生物组织发生吸收、散射作用，导致自身的权重衰减、传播方向改变。行走到组织表面时，从上表面逃逸变成漫反射光，从下表面逃逸变成透射光。当该光子自身的权值小于某一个给定的阈值时，光子被组织吸收。最后，对所有光子重复上述过程，实现了整个光束的MC模拟，MC方法模拟光子在皮肤中传输的流程图如图4所示。具体过程如下：

1)发射光子：入射光束垂直组织平面xoy发射，方向沿着z轴方向，初始位置为原点(0,0,0)，则方向余弦的初始值为(0,0,1)。令光子的初始权值为1，经过反射作用，光子的权值减小至M＝1-R_sp。其中，R_sp是反射系数，R_sp＝(n₀-n₁)²/(n₀+n₁)²，n₀是组织外部介质的折射率，n₁是组织表面的折射率。

2)光子的步长：设区间(a,b)上随机变量η的归一化概率密度函数为p(η)：

设随机变量ξ在区间(0,1)上服从均匀分布，并且η和ξ之间一一映射，即η＝f(ξ)。η＝f(ξ)是增函数时，经过推导一定存在ξ₁∈(0,1)，满足下式：

η＝f(ξ)是减函数时，经过推导也一定存在ξ₁∈(0,1)，满足下式：

对变量η进行抽样，可以得到随机变量η和随机变量ξ之间的关系式η＝f(ξ)。变量ξ通过随机数产生器获取。由于ξ₁和1-ξ₁分布一样，且彼此间能够相互转化，因此运用式(20)对随机变量进行抽样。

光子的步长s是随机的，s的计算基于其概率分布的采样。依据互相作用系数μ_t的定义(μ_t＝μ_α+μ_s)，区间(s',s'+ds')中每单位长度的光子-组织互相作用概率为：

或者写成下式：

d(ln(P{s≥s'}))＝-μ_tds'   (22)

将式(22)在(0,s₁)范围内对s'积分，可以得到指数分布，其中P(s≥0)＝1：

经过重新排列得到步长s的累积分布函数：

则步长s的概率密度函数为：

将式(25)代入到式(19)得：

s₁＝-ln(1-ξ)/μ_t   (26)

因为ξ的对称性约为0.5，我们可以用ξ替代1-ξ得到光子的步长：

s₁＝-ln(ξ)/μ_t   (27)

3)光子的传输：光子当前位置用坐标(x,y,z)表示，方向余弦(u_x,u_y,u_z)决定其传输方向。光子到下一个位置时，需要对坐标进行更新为(x',y',z')：

x'＝x+s·u_x

y'＝y+s·u_y   (28)

z'＝z+s·u_z

4)光子的吸收：光子在生物组织中传输时会与内部粒子发生吸收作用，导致光子的权值衰减，衰减大小为Δw＝w·μ_α/μ_t，则当前光子的权值变为：

w_left＝w-Δw＝w-w·μ_α/μ_t＝w·μ_s/μ_t   (29)

设发射光子的初始权值为w₀，光子与组织内部粒子发生n次吸收作用后，光子的权值变为w_n：

可见，光子与组织产生的吸收作用次数越多，权值衰减越快。

5)光子的散射：光子在发生散射时，其传播方向会跟着改变，进而方向余弦一同改变。θ和ψ一旦确定，光子的方向余弦更新为：

如果光子的方向足够靠近z轴(如|u_z|＞0.99999)，那么式(31)修改为下式：

u′_x＝sinθcosψ

u′_y＝sinθsinψ   (32)

u′_z＝sign(u_z)cosθ

用符号函数sign(·)来判断光子z轴的方向。

6)光子在边界上的反射和透射：组织边界分为上、下表面，当光子抵达组织边界时，会产生逃逸出组织和反射回组织两种情况。约化步长s₁指光子顺着当前方向与组织边界的间隔，如图2所示。s₁由下式表示：

在直角坐标系中，z₀和z₁分别表示组织上、下表面的z轴坐标。若光子移动的步长s＜s₁，光子处于组织内部；若光子移动的步长s≥s₁，光子将会达到边界，产生全反射或折射。到达组织边界时入射角为α_i＝arccos(|u_z|)，若产生折射，记折射角为α_t，由折射定律可得：

n_i sinα_i＝n_t sinα_t   (34)

若产生全反射，则反射率R(α_i)可由菲涅耳公式得：

若随机数ξ＞R(α_i)，光子产生折射，逃逸出组织，上表面逃逸变成漫反射光，下表面逃逸变成透射光；反之若ξ≤R(α_i)，光子产生全反射，因此方向余弦要更新为(u_x,u_y,u_z)，步长更新为s-s₁。若光子步长不大，在行走完随机步长时光子便会与组织产生吸收和散射作用。

7)组织交界处的反射和透射：假设光子在第i层中传输，光学特性参数为μ_αi、μ_si和n_i，步长s的方向余弦为(u_x,u_y,u_z)。若s＞s₁，光子会行走到第i层与第i-1层或第j＝i+1层(光学特性参数为μ_αj、μ_sj和n_j)的交界处，此时由式(22)判断光子的情况。若产生全反射，光子继续在第i层传输，步长更新为s-s₁，方向余弦更新为(u_x,u_y,u_z)。若产生折射，光子传输到第j＝i+1层。由于两层的散射系数不一样，于是光子的步长更新为s←s·μ_si/μ_sj，方向余弦也做出相应的改变，如图3所示。

8)光子的终止：以下两种情况下，光子被认为终止：①光子的权值衰减为0；②光子逃逸。对于第一种情况，设置一个阈值w_ε(如w_ε＝1e-6)，当光子的权值w≤w_ε时，光子达到终止条件。不过为了更加契合实际情况，引入“俄罗斯轮盘赌”的办法：给定参数M(如M＝100)，产生一个服从均匀分布的随机数ξ。当ξ＞1/M，光子的权值变为0，终止；当ξ≤1/M，光之权值变为Mw，继续传输。

步骤三、皮肤组织成像：相机拍摄的可见光范围内的图像，其RGB颜色值由光源、皮肤漫反射率和相机光谱响应函数共同决定。即：

其中，R,G,B表示图像的颜色值，I(λ)是光源。λ是波长，取值范围是400～800nm。R(λ)是皮肤的漫反射率，S(λ)是相机的光谱响应函数，可由相机模型得到。生发层的黑色素体积分数C_mel、真皮层的血液体积分数C_blood和真皮层的厚度D^dermis这三个生理参数直接影响血管在皮肤图像中的颜色，使其在各自取值范围内进行变化，得到对应的漫反射率。将I(λ)、R(λ)、S(λ)的数据导入，生成RGB图像，得到颜色值数据。

使用一个三层前馈神经网络模型来训练皮肤图像RGB值到皮肤生理参数C_mel，C_blood，D^dermis之间的映射关系。输入层有3个参数，以MC模拟得到的RGB为输入；输出层有3个参数，以生理参数C_mel,C_blood,D^dermis为输出；隐含层有5个神经元。隐含层和输出层使用的激活函数分别为对称型s函数和线性函数，利用BP算法进行训练。对于一张待测的皮肤图像，输入其RGB值，便可从输出的映射图像中实现静脉的显像。

以下根据具体的实例阐述本发明的方法、结果：

将单反相机NikonD70拍摄的测试皮肤图像作为输入，网络输出得到生理参数值，进而得到静脉显像图像，如图5所示，图5静脉显像结果，Cm、Cb、depth分别为生发层的黑色素体积分数、真皮层的血液体积分数和真皮层的厚度。从图中可以看出，每一个生理参数均能够使皮肤图像达到静脉显像的效果。总体来说，在这三个生理参数中，血液的体积分数C_blood拼接出的静脉显像的效果相对稳定。

为了进一步验证静脉显像的正确性，通过血管增强操作进一步去除噪声，最后提取血管位置，如图6所示。以JAI同步相机拍摄的近红外图像为真实值，对JAI同步相机拍摄的可见光图像进行静脉显像，评价算法效果。评价指标为准确率：

准确率体现的是模型对于静脉显像的精确性，其值越大越好，说明该模型得到的静脉显像图像更加接近真实结果。部分图像对及其静脉显像图像如图7所示，评价指标数值如表1所示。

表1显像准确率指标

从表1可知，该模型的平均准确率达到了97.8749％，说明MC模拟方法结合神经网络可以十分准确地达到静脉显像的效果，且非常接近真实静脉脉络结构的分布情况。

下面验证本发明算法对不同相机拍摄的不同身体部位的可见光图像依然有很好的静脉显像效果，结果如图8所示，其中给出了C_blood图的结果。由图8可知，对于不同相机拍摄的不同身体部位，本发明算法还是能够清晰地对RGB图像进行静脉检测与显像，从而验证了本发明算法的鲁棒性。从中亦可以看出，基于MC的皮肤成像模拟并进行静脉显像的方法具有很高的实用价值。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进，这些改进也应视为本发明的保护范围。

Claims (6)

1.一种基于蒙特卡洛的可见光静脉显像方法，其特征在于，所述的方法的步骤如下：

步骤一、建立皮肤模型，对皮肤进行光学模型分析；所述的步骤一具体如下：1.1、分析光在皮肤组织中的吸收作用，七层皮肤的吸收系数如下：

角质层：角质层中对光的吸收主要是水和组织的固有吸收，公式如下：

其中，是除了水和血液之外的其它组织的固有吸收系数，是水的体积分数，是水的吸收系数；

生发层：生发层中的吸收体主要为黑色素，其吸收系数表达式为：

其中，是黑色素的吸收系数，C_mel是黑色素的体积分数，其取值范围是1.3％～43％；

真皮层：真皮层包含真皮乳头层、上层血管网真皮层、网状真皮层、深部血管网真皮层四层，根据血氧饱和度S和血液中血红蛋白总的体积分数γ，该层的吸收系数可由下式表达：

其中，是水的体积分数；C_blood是血液的体积分数；γ＝F_HbF_RBCH_t，

F_Hb是红细胞中血红蛋白的体积分数，F_RBC是红细胞在所有血细胞中的体积分数，H_t是红细胞比容；是脱氧血红蛋白的吸收系数，是有氧血红蛋白的吸收系数，表达式如下：

其中c^hb是血液中血红蛋白的浓度，值为150g/L，m^hb是血红蛋白的分子量，值为66500g/mole，ε^ohb(λ)和ε^dhb(λ)分别为有氧血红蛋白和脱氧血红蛋白的消光系数；

分析七层皮肤的光的散射系数μ_s，光的散射系数μ_s是指每单位长度光被散射的概率，单位是cm^-1，七层皮肤的各层散射系数如下：

角质层的散射系数：

其中，λ是波长；

生发层的散射系数：

真皮层的散射系数：

皮下脂肪层的散射系数：

1.2、分析七层皮肤的折射率：

角质层的折射率：

n^str＝1.55                      (13)

生发层的折射率：

n^liv＝1.684-1.8723×10⁴λ^-2+1.0964×10¹⁰λ^-4-8.6484×10¹⁴λ^-6  (14)

真皮层的折射率：

n^dermis＝1.309-4.346×10²λ^-2+1.6065×10⁹λ^-4-1.2811×10¹⁴λ^-6 (15)

皮下脂肪层的直射率：

n^fat＝1.44                      (16)；

步骤二、基于Monte Carlo方法进行皮肤成像模拟；

步骤三、皮肤组织成像，最终实现静脉的显像。

2.根据权利要求1所述的一种基于蒙特卡洛的可见光静脉显像方法，其特征在于，所述的步骤一中建立角质层、生发层、真皮层、皮下脂肪层的七层皮肤模型，所述的真皮层依次包括真皮乳头层、上层血管网真皮层、网状真皮层、深部血管网真皮层，对七层皮肤进行光学模型分析；在利用Monte Carlo方法进行皮肤组织模拟时，需要输入五个光学特性参数值：吸收系数μ_α(cm^-1)、散射系数μ_s(cm^-1)、各向异性因子g、折射率n和皮肤各层厚度d(cm)。

3.根据权利要求2所述的一种基于蒙特卡洛的可见光静脉显像方法，其特征在于，分析皮肤各层厚度的参数影响，皮肤各层的厚度变化成正相关，即同时变厚，或者同时变薄；以上层血管网真皮层的厚度(cm)为基准，构造线性映射函数如下：

其中，d^pap是真皮乳头层的厚度；d^upp是上层血管网真皮层的厚度；d^ret是网状真皮层的厚度；d^dee是深部血管网真皮层的厚度；d^fat是皮下脂肪层的厚度。

4.根据权利要求1所述的一种基于蒙特卡洛的可见光静脉显像方法，其特征在于，所述的步骤二具体如下：

2.1，发射光子：

入射光束垂直组织平面xoy发射，方向沿着z轴方向，初始位置为原点(0,0,0)，则方向余弦的初始值为(0,0,1)；令光子的初始权值为1，经过反射作用，光子的权值减小至M＝1-R_sp；其中，R_sp是反射系数，R_sp＝(n₀-n₁)²/(n₀+n₁)²，n₀是组织外部介质的折射率，n₁是组织表面的折射率；

2.2，光子的步长：

设区间(a,b)上随机变量η的归一化概率密度函数为p(η)：

设随机变量ξ在区间(0,1)上服从均匀分布，并且η和ξ之间一一映射，即

η＝f(ξ)；η＝f(ξ)是增函数时，经过推导一定存在ξ₁∈(0,1)，满足下式：

η＝f(ξ)是减函数时，经过推导也一定存在ξ₁∈(0,1)，满足下式：

对变量η进行抽样，可以得到随机变量η和随机变量ξ之间的关系式η＝f(ξ)；变量ξ通过随机数产生器获取；

运用式(20)对随机变量进行抽样：光子的步长s是随机的，s的计算基于其概率分布的采样，依据互相作用系数μ_t的定义(μ_t＝μ_α+μ_s)，区间(s',s'+ds')中每单位长度的光子-组织互相作用概率为：

或者写成下式：

d(ln(P{s≥s'}))＝-μ_tds'           (22)

将式(22)在(0,s₁)范围内对s'积分，得到指数分布，其中P(s≥0)＝1：

经过重新排列得到步长s的累积分布函数：

则步长s的概率密度函数为：

将式(25)代入到式(19)得：

s₁＝-ln(1-ξ)/μ_t                 (26)

ξ的对称性为0.5，用ξ替代1-ξ得到光子的步长：

s₁＝-ln(ξ)/μ_t                  (27)

2.3，光子的传输：

光子当前位置用坐标(x,y,z)表示，方向余弦(u_x,u_y,u_z)决定其传输方向；光子到下一个位置时，对坐标进行更新为(x',y',z')：

2.4，光子的吸收：

光子在生物组织中传输时会与内部粒子发生吸收作用，导致光子的权值衰减，衰减大小为Δw＝w·μ_α/μ_t，则当前光子的权值变为：

w_left＝w-Δw＝w-w·μ_α/μ_t＝w·μ_s/μ_t          (29)

设发射光子的初始权值为w₀，光子与组织内部粒子发生n次吸收作用后，光子的权值变为w_n：

2.5、光子的散射：

光子在发生散射时，其传播方向会跟着改变，进而方向余弦一同改变；θ和ψ一旦确定，光子的方向余弦更新为：

如果光子的方向靠近z轴，|u_z|＞0.99999，那么式(31)修改为下式：

u′_x＝sinθcosψ

u′_y＝sinθsinψ   (32)

u′_z＝sign(u_z)cosθ

用符号函数sign(·)来判断光子z轴的方向；

2.6、光子在边界上的反射和透射：

组织边界分为上、下表面，当光子抵达组织边界时，会产生逃逸出组织和反射回组织两种情况，设定光子的步长s和约化步长s₁，约化步长s₁指光子顺着当前方向与组织边界的间隔，s₁由下式表示：

在直角坐标系中，z₀和z₁分别表示组织上、下表面的z轴坐标；若光子移动的步长s＜s₁，光子处于组织内部；若光子移动的步长s≥s₁，光子将会达到边界，产生全反射或折射；到达组织边界时入射角为α_i＝arccos(|u_z|)，若产生折射，记折射角为α_t，由折射定律可得：

n_i sinα_i＝n_t sinα_t   (34)

若产生全反射，则反射率R(α_i)可由菲涅耳公式得：

若随机数ξ＞R(α_i)，光子产生折射，逃逸出组织，上表面逃逸变成漫反射光，下表面逃逸变成透射光；反之若ξ≤R(α_i)，光子产生全反射，因此方向余弦要更新为(u_x,u_y,u_z)，步长更新为s-s₁；若光子步长不大，在行走完随机步长时光子便会与组织产生吸收和散射作用；

2.7、组织交界处的反射和透射：

假设光子在第i层中传输，光学特性参数为第i层的吸收系数μ_αi、第i层的散射系数μ_si和第i层的折射率n_i，步长s的方向余弦为(u_x,u_y,u_z)；若s＞s₁，光子会行走到第i层与第i-1层或第j＝i+1层，光学特性参数为第j层的吸收系数μ_αj、第j层的散射系数μ_sj和第j层的折射率n_j的交界处，此时由式(22)判断光子的情况；若产生全反射，光子继续在第i层传输，步长更新为s-s₁，方向余弦更新为(u_x,u_y,u_z)；若产生折射，光子传输到第j＝i+1层；

2.8：光子的终止：

光子的终止情况有：1)光子的权值衰减为0；2)光子逃逸；对于光子的权值衰减为0，设置一个阈值w_ε，当光子的权值w≤w_ε时，光子达到终止条件。

5.根据权利要求4所述的一种基于蒙特卡洛的可见光静脉显像方法，其特征在于，所述的步骤2.8中对于光子的终止的光子的权值衰减为0的情况时，还可以给定参数M，产生一个服从均匀分布的随机数ξ；当ξ＞1/M，光子的权值变为0，终止；当ξ≤1/M，光之权值变为Mw，继续传输。

6.根据权利要求1所述的一种基于蒙特卡洛的可见光静脉显像方法，其特征在于，所述的步骤三具体为：

3.1、相机拍摄的可见光范围内的图像，其RGB颜色值由光源、皮肤漫反射率和相机光谱响应函数共同决定；即：

其中，R,G,B表示图像的颜色值，I(λ)是光源，λ是波长，取值范围是400～800nm；R(λ)是皮肤的漫反射率，S(λ)是相机的光谱响应函数，可由相机模型得到；将I(λ)、R(λ)、S(λ)的数据导入，生成RGB图像，得到颜色值数据；

3.2、使用一个三层前馈神经网络模型来训练皮肤图像RGB值到皮肤生理参数C_mel，C_blood，D^dermis之间的映射关系；输入层有3个参数，以MC模拟得到的RGB为输入；输出层有3个参数，以生理参数C_mel,C_blood,D^dermis为输出；隐含层有5个神经元；

3.3、隐含层和输出层使用的激活函数分别为对称型s函数和线性函数，利用BP算法进行训练；对于一张待测的皮肤图像，输入其RGB值，便可从输出的映射图像中实现静脉的显像。

Images