CN100373396C - 基于图像分割和二维小波变换的虹膜识别方法 - Google Patents

Abstract

基于图像分割和二维小波变换的虹膜识别方法，属于生物特征模式识别技术领域。先把虹膜定位分成内边缘定位和外边缘定位两个部分，并把重点放在内边缘定位上。接着采用直角坐标与极坐标映射理论把定位后的图像归一化到一个固定的灰度矩阵中。并对归一化后的图像进行两次图像分割，最终分割成18个小区域。再通过对归一化处理后虹膜图像的各个小分割区域进行二维小波变换，提取出主要小波通道的小波系数均值和方差作为特征值。在匹配识别算法中，采用方差倒数加权求和的匹配算法进行识别判决。最后用不同的系数来加权这3个识别结果，得到最终的识别结果。与已有的虹膜识别方法相比本发明在满足系统实时性要求的前提下具有较好的抗噪性能和较高的识别率。

Description

基于图像分割和二维小波变换的虹膜识别方法

技术领域

基于图像分割和二维小波变换的虹膜识别方法，属于生物特征模式识别技术领域，特别涉及虹膜特征识别方法。

背景技术

随着信息技术的发展和电子商务的广泛应用，信息安全日渐成为人们面临的一个重要而迫切的问题。可用于身份鉴别、保护信息安全的生物特征识别技术越来越受到人们的重视。所谓生物特征识别技术是指通过计算机与光学、声学和生物统计学等高科技手段密切结合，利用人体固有的生理特性和行为特征来进行个人身份的鉴定。而虹膜特征识别技术就是利用人眼虹膜纹理的不同来进行身份识别的。

人眼的虹膜中含有极丰富的信息。它的表面有一些类似于细丝、斑点、漩涡、冠状等形状的纹理。虹膜这些纹理具有唯一性，不同的人具有不同的虹膜纹理特性，即便是同一个人，它左眼和右眼的虹膜纹理特征都不一样。因此，用这些纹理特征进行身份识别是非常精确的。而且虹膜的纹理特征主要是由人还在胚胎时的环境所决定的，并且由于虹膜的外部有透明的角膜将其与外界相隔离，所以发育完全的虹膜不易受到外界的伤害而产生变化。因此虹膜识别具有绝对的可靠性。另外，由于瞳孔会根据光线的强弱发生变化，进而牵动虹膜形状也跟着发生变化，利用这一点可以鉴别用于识别的虹膜样本是否为活眼的虹膜，因此虹膜识别也具有较高的防伪性。正是由于虹膜纹理具有这些优点，使得基于虹膜特征的身份识别技术在金融、电子商务，安全保卫等各个方面都具有极大的应用前景。

基于虹膜特征的身份识别技术在国外已经获得了蓬勃的发展，并且已经逐步开展了虹膜识别的商品化过程。虹膜识别系统最初是由英国剑桥大学开发成功的。1993年，英国剑桥大学的John.G.Daugman给出了较完整的虹膜识别方法。这一方法准确性高，速度快，是当今世界上几乎所有商用虹膜识别系统的理论基础，他的开创性工作使得自动的虹膜识别成为可能。1996年，普林斯顿的Wildes等人研制成功了基于区域图像注册技术的虹膜识别系统。1998年，昆士兰大学的Boles等人提出了基于零交叉小波变换的虹膜识别方法。详见文献：J.G.Daugman.High Confidence Visual Recognition of Persons by a Test of Statistical Independence.IEEE transactions on Pattern Analysis and Machine Intelligence，1993，15(11)：1148-1161、文献：R.P.Wildes，J.C.Asmuth，G. L. Green，et al.A Machine-Vision System for Iris Recognition.Machine Vision and Applications，1996，9(1)：1-8和文献：W.Boles，B.Boashah.A HumanIdentification Technique Using Images of the Iris and Wavelet Transform.IEEE transactions onSignal Processing，1998，46(4)：1185-1188。

国内对虹膜识别技术的研究起步较晚，在最近几年内也获得了较快的发展。但与国外虹膜识别产业蓬勃兴旺的发展势头相比，还是有一定差距。目前，中国科学院自动化所已经完成了虹膜识别的实验室阶段的研究，并且申请了虹膜采集装置的专利。上海交通大学、浙江大学、华中科技大学等也在进行相关的研究，并也都取得了一定的研究成果。详见文献：王蕴红，朱勇，谭铁牛.基于虹膜识别的身份鉴别.自动化学报，2002，28(1)：1-10、文献：应忍冬，徐国治.基于小波变换过零检测的虹膜识别技术.上海交通大学学报，2002，36(3)：355-358和文献：陈良洲，叶虎年.一种新的虹膜识别算法研究.华北工学院测试技术学报，2000，14(4)：211-216。

目前，国内外已提出的虹膜识别方法中，已在实际应用中取得了较好的识别效果的方法有：

1、Daugman的基于相位分析的虹膜识别方法。它采用Gabor小波滤波的方法编码虹膜的相位特征。2D Gabor函数可以较好地达到频率域和空间域的局部化，就是说在空间局部化的同时具有很好的频率和方向选择性。通过计算2D Gabor相位系数，可以有效地从纹理中提取连续和不连续的纹理信息。详见文献：J.G.Daugman.High Confidence Visual Recognition ofPersons by a Test of Statistical Independence.IEEE transactions on Pattern Analysis and MachineIntelligence，1993，15(11)：1148-1161。

2、Boles的过零点检测的虹膜识别方法。它采用一维小波变换对沿虹膜中心同心圆的一条采样曲线进行过零点检测，通过两个自定义的相似度函数完成虹膜特征的分类。其理论基础为Mallat的信号过零点描述重建理论。详见文献：W.Boles，B.Boashah.A HumanIdentification Technique Using Images of the Iris and Wavelet Transform.IEEE transactions onSignal Processing，1998，46(4)：1185-1188。

3、谭铁牛，王蕴红等人的基于纹理分析的虹膜识别方法。该方法将虹膜的特征看作一种随机纹理，从纹理分析的角度提取虹膜的局部特征。该方法采用Gabor滤波和以Daubechies-4小波为小波基的二维小波变换提取虹膜纹理特征，并用方差倒数加权欧氏距离方法进行特征匹配，取得了较好的识别效果，详见文献：王蕴红，朱勇，谭铁牛.基于虹膜识别的身份鉴别.自动化学报，2002，28(1)：1-10。

上述3种虹膜识别方法虽然都取得了较好的识别效果，但仍有不完善的地方。方法1虽然取得了较高的识别精确度，但由于其提取的虹膜特征向量的维数较高，达到了2048维。因此对采集到的虹膜图像的清晰度有较高的要求。方法2虽然克服了以往系统受漂移、旋转和比例放缩带来的局限，也对亮度变化及噪声不敏感，并对采集图像的质量要求不是很高，但由于该算法只利用了一部分虹膜纹理特征信息，所以并没有获得较高的正确识别率。方法3采用了不同的纹理分析策略提取纹理特征，获得了较快的运行速度，但它对虹膜纹理描述比较粗糙，其正确识别率仍然不是很高。

评价一种虹膜识别方法的优劣有两个较重要的指标：识别率和运行速度。通常，这两个指标是矛盾的。一种好的虹膜识别方法是要在满足系统实时性要求的运行速度的前提下，最大限度的提高识别率。

发明内容

本发明提出的虹膜识别方法在满足实时性要求的前提下获得了较高的虹膜识别率。它通过图像分割较好地虑除了虹膜纹理图像中的噪声干扰，再通过二维小波变换较全面的提取出了二维虹膜纹理图像的特征信息，最后通过方差倒数加权特征求和的方法得到最后的识别结果。

为了方便描述本发明的内容，首先在此作一个术语定义：

1.虹膜：人眼球中介于瞳孔和巩膜中间的组织，该部分组只含有独特的，丰富的纹理信息，可用于进行身份识别。它的几何形状为环形形状。

2.虹膜内、外边缘：虹膜与瞳孔的交界处称为虹膜的内变缘，它是一个圆形。虹膜与巩膜的交界处称为虹膜的外边缘，它也是一个圆形。

3.虹膜图像采集装置：捕获虹膜图像数字信号的装置。

4.灰度图像：图像中只包含了亮度信息而没有任何颜色信息的图像。

5.中值滤波：是一种非线性的信号处理方法，也是一种典型的低通滤波器。它将领域中的像素按灰度级排序，取其中间值为输出像素值。中值滤波的效果依赖于两个要素：领域的空间范围和中值计算中涉及的像素数。

6.灰度直方图：灰度直方图(histogram)是灰度级的函数，它表示图像中具有每种灰度级的象素的个数，反映图像中每种灰度出现的频率。灰度直方图的横坐标是灰度级，纵坐标是该灰度级出现的频率，是图像的最基本的统计特征。

7.图像二值化：设定一个阈值，把图像中每一个象素点的灰度值都与这个阈值相比较，当灰度值大于阈值时即把该象素点的灰度值设为1；当灰度值小于阈值时即把该点的灰度值设为0。把图像中象素点灰度值变为只有0，1的过程称为图像二值化。

8.虹膜定位：在一幅包含瞳孔、巩膜、眼睫毛的虹膜图像中，准确的定位出环形虹膜在图像中的几何位置的过程。

9.Roberts算子：它是一种常用的边缘检测算子。Roberts算子基于用图像的垂直和水平

差分来逼近梯度算子。它由两个2×2模板共同实现：

其计算公式为： $R(i,j)=\sqrt{{(f(i,j)-f(i+1,j+1))}^2+{(f(i,j+1)-f(i+1,j))}^2}.$

10.x方向的灰度投影量P(x)：在灰度矩阵B(x，y)中，在x方向上把每个x值对应的象素灰度都加在一起。即 $P\left(x\right)=\underset{y}{\mathrm{\Σ}}B(x,y).$

11.y方向的灰度投影量P(y)：在灰度矩阵B(x，y)中，在y方向上把每个y值对应的象素灰度都加在一起。即 $P\left(y\right)=\underset{x}{\mathrm{\Σ}}B(x,y).$

12.圆形边缘检测器：它的基本数学算子为：

其基本思想就是在参数空间不断迭代(r，x₀，y₀)的值，由于每个参数值(r，x₀，y₀)都对应着一个圆，所以迭代(r，x₀，y₀)的过程，也即是迭代圆的过程。在每一个圆的圆周上，求灰度值的圆积分。随着(r，x₀，y₀)的循环变化，灰度值圆积分也跟着变化，找到圆积分变化最大的那个圆即为要检测的圆。

13.虹膜图像归一化：由于每次采集虹膜图像时虹膜的位置都会有所不同，而且由于采集系统光照的影响，会导致瞳孔的放大或收缩，这会使虹膜的大小也跟着变化。所以定位出圆环形的虹膜区域不能直接用来特征提取，必须把圆环形虹膜转换成一个大小固定的灰度矩阵图像。该过程即为图像归一化过程。

14.直方图均衡化：它是把原始图像的直方图变换为均匀的形式，这样就增加了象素灰度值的动态范围，增强了图像整体上的对比度。经过均衡化后的图像，其清晰程度有明显的提高，所需的目标信息就会被突出出来。直方图均衡化能解决由于非均匀光照对虹膜结果造成影响的问题。

15.二维小波变换：对一幅虹膜图像进行二维小波变换可以提取水平方向、垂直方向和对角方向的细节系数。它非常适合分析二维图像信号。

16.Haar小波：它是一种较简单的小波基，它的小波函数为： $\mathrm{\ψ}\left(t\right)=\left\{\begin{array}{c}1,t\∈\[0,1/2\]\\ -1,t\∈\[1/2,1\]\end{array}\right..$

17.小波通道：对一幅图像进行完全的小波分解得到一系列的小波系数。通常把这些小波分解系数构成的子图像称为小波分解通道。对一幅图像进行1阶二维小波变换可得到四种小波通道：LL，LH，HL，HH。每个通道表征了原始图像不同空间频率和方向下的信息。

18.均值：它表示小波通道中，小波系数的平均值。它表征了小波通道的能量大小。它的计算公式为： $E_n=\frac{1}{M\×N}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\left|x(i,j)\right|.$

19.样本方差：它用来衡量小波通道中，小波系数与均值的偏离程度。它的计算公式为：

$D_n=\frac{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\[\left|x(i,j)\right|-E_n\]}^2}{M\×N-1}.$

20.系统学习阶段：系统把采集到的虹膜图像，提取出完整的特征值。并把该特征值存入虹膜样本数据库，以作为匹配识别的标准样本。

21.系统识别阶段：系统读入未知身份的虹膜图像，提取一半的特征值。并把该特征值与虹膜样本数据库中的特征值按照一定的匹配识别算法进行比较，最终得出未知身份虹膜图像的识别结果。

本发明详细技术方案为：

基于图像分割和二维小波变换的虹膜识别方法，其特征是，它包含下列步骤：

步骤1：采集虹膜图像。

通过虹膜图像的采集装置采集虹膜图像，得到可用于进一步处理的虹膜图像灰度矩阵H(x，y)。

步骤2：虹膜图像中值滤波。

对步骤1所得的虹膜图像灰度矩阵H(x，y)做平滑处理得到灰度矩阵I(x，y)。在平滑处理中，滤波器为采样窗口大小为9个离散象素点的非线性中值滤波器。

步骤3：虹膜内边缘定位。具体包括以下步骤：

步骤1)：图像二值化

求出灰度矩阵I(x，y)的灰度直方图，找出灰度直方图内中灰度值在(20～125)范围内的峰值所对应的灰度值M，把灰度值M与一安全系数D相加得到可用于灰度图像二值化的阈值Y。安全系数D通常取3到7之间的数。以Y为阈值对虹膜灰度图像I(x，y)进行二值化，得到二值化图像B(x，y)。

步骤2)：找内边缘粗略圆心

计算出二值化矩阵B在x方向的灰度投影量P(x)，在y方向的灰度投影量P(y)。在一维数组P(x)中，找到其最小值，并找到该最小值对应的x₁；同理在一维数组P(y)中，找到其最小值，并找到该最小值对应的y₁。(x₁，y₁)即为虹膜内边缘的粗略圆心。

步骤3)：边缘检测

对二值化后的图像B(x，y)用Roberts算子进行图像的边缘检测，得到一个含有边缘的二值图像BW(x，y)。

步骤4)：将边缘点分成4个象限

在BW(x，y)中以步骤3)中得到的内边缘粗略圆心(x₁，y₁)为坐标原点建立坐标系。把离散边缘点分为四个象限。在每个象限内，在距离原点30个像素和50个像素的扇形范围内，随机选取3个像素相差在10个单位以上的边缘点。同理，在其他象限内也是如此选择3个像素点。这样4个象限共有12个像素点。

步骤5)：4个象限联合精确定位

在步骤4)的12个点中选出3个不在同一直线上的点，由这个3个点就能构成一个圆。并求出其余9个点到这个圆的距离d。12个点最多能构成220个圆，也即有220个距离d，在这220个距离中找到最小的d对应的那个圆即为虹膜的内边缘。设得到的虹膜内边缘的圆心为(x_a,y_a)，半径为r_a。

步骤4：虹膜外边缘定位。具体包括以下步骤：

步骤1)：限定圆边缘检测模板的迭代范围

在步骤2得到I(x，y)中，用圆形边缘检测器进行迭代求灰度积分值。在迭代过程中把内边缘的圆心(x_a，y_a)作为迭代(x₀，y₀)的初始值。并把(x₀，y₀)的搜索范围限定在以(x_a-10，y_a-10)、(x_a-10，y_a+10)、(x_a+10，y_a-10)、(x_a+10，y_a+10)为顶点的矩形内。并把r的搜索范围限定在70～110个单位像素内。在搜索的过程中，并不是在整个圆上对灰度做圆积分，而是在以(x₀，y₀)建立的直角坐标系中，角度为-45°～45°、135°～225°的圆形弧线上求灰度的积分。

步骤2)：迭代中找出外边缘

在步骤1)的范围内迭代参数空间的(r，x₀，y₀)的值，求出灰度积分变化最大的那个圆即为虹膜的外边缘。对应的参数值(r，x₀，y₀)即为虹膜外边缘的半径和圆心值(r_b，x_b，y_b)。定位结果如图1所示。

步骤5：虹膜图像归一化。具体包括以下步骤：

步骤1)：建立坐标转换模型、并将虹膜图像归一化

由步骤3和步骤4得到虹膜内、外边缘的圆周参数(r_a，x_a，y_a)、(r_b，x_b，y_b)，以内圆圆心(x_a，y_a)作为坐标系的原点建立将直角坐标转换为极坐标的数学模型(如图2所示)。在该模型中从原点开始做与水平线成θ角的射线，它与内、外边界各有一个交点，分别记作B(x_i，y_i)，A(x_o，y_o)。射线上两个交点A，B之间的任何一点的坐标(x，y)都可以用A(x_o(θ)，y_o(θ))，B(x_i(θ)，y_i(θ))的线性组合来表示： $\left\{\begin{array}{c}x(r,\mathrm{\θ})=(1-r)*x_i\left(\mathrm{\θ}\right)+r*x_0\left(\mathrm{\θ}\right)\\ y(r,\mathrm{\θ})=(1-r)*y_i\left(\mathrm{\θ}\right)+r*y_0\left(\mathrm{\θ}\right)\end{array}\right..$ 这样就可以将虹膜图像归一化为以θ为横轴，r为纵轴的固定尺寸的灰度矩阵P(x，y)，虹膜归一化处理结果如图3所示。

步骤2)：直方图均衡化

对步骤1)得到的灰度矩阵P(x，y)做直方图均衡化处理得到归一化灰度矩阵PI(x，y)。

步骤6：归一化后灰度矩阵第一次图像分割

在PI(x，y)中，选取灰度矩阵上方16×1024的部分作为虹膜纹理A区域。选取行宽为17～48，列的范围分别是1～128、384～640、896～1024共3小块区域作为虹膜纹理B区域。选取行宽为49～64，列的范围分别是1～64、448～576、980～1024共3小块区域作为虹膜纹理C区域，虹膜图像第一次分割如图4所示。

步骤7：归一化后灰度矩阵第二次图像分割。

把步骤6中的纹理A区域分成8个小块的灰度矩阵：行宽为1～16，列宽分别为1～128、129～256、257～384、385～512、513～640、641～768、769～896、897～1024。把虹膜纹理B区域也分成8个小块的灰度矩阵：行宽为17～32，列宽分别为1～128、385～512、513～640、897～1024；行宽为33～48，列宽分别为1～128、385～512、513～640、897～1024。把虹膜纹理C区域分成2个小块的灰度矩阵：行宽为49～64，列宽为1～64、961～1024的两个灰度矩阵区域结合在一起；行宽为49～64，列宽为449～578的灰度矩阵区域。虹膜图像第二次分割如图5所示。

步骤8：对每个分割区域进行二维小波变换

对步骤7中的每个小分割区域以Haar小波为小波基进行3阶二维小波变换，二维小波变换后共得到10个小波通道。把这些通道分别记为LL³、LH³、HL³、HH³、LH²、HL²、HH²、LH¹、HL¹、HH¹。其中HH¹、HH²、HH³三通道代表了虹膜图像在水平高频和垂直高频下的信息，它们含有大量的噪声，不利于提取虹膜特征。舍弃这三个通道只保留剩余的7个通道，如图6所示。

步骤9：提取小波系数均值和方差作为特征值

对某个小分割区域二维小波变换后的每个小波通道提取其均值和样本方差：

$E_n=\frac{1}{M\×N}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\left|x(i,j)\right|,$ $D_n=\frac{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\[\left|x(i,j)\right|-E_n\]}^2}{M\×N-1}$ 作为特征值。这样该小分割区域每个小波通道能提取出两个特征值。7个小波通道就能提取出14个特征值。对每个小分割区域都重复此过程。每个小分割区域都能提取出14个特征值，共有18个小分割区域，故能提取出252个特征值。

步骤10：用方差倒数加权均值差进行匹配识别。具体包括以下步骤：

步骤1)：提取完整特征做为虹膜特征样本

在系统的学习阶段，对要学习的每幅虹膜图像都按照步骤2到步骤9的处理过程进行处理。每幅虹膜图像能提取出252个特征值，把这些特征值作为该虹膜图像的样本存储在虹膜样本数据库中，用于后面识别判定。

步骤2)：提取均值特征用于识别

在系统的识别阶段，对一幅未知的虹膜图像，对该图像按照步骤2到步骤9的处理过程进行处理。并在步骤9中只提取均值E_n作为特征值。这样用于识别的虹膜图像只需提取出128个特征值。

步骤3)：分3个部分分别进行匹配识别

在匹配识别过程中，把样本中的252个特征值按照步骤6的三个纹理区域分成3个部分：E_A、E_B、E_C。E_A为虹膜纹理A区域经变换提取的特征值。E_B为虹膜纹理B区域经变换提取的特征值。E_c为虹膜纹理C区域经变换提取的特征值。同理也把待识别虹膜图像的128个特征值分成3个部分：e_A、e_B、e_C。按照方差倒数加权求和的匹配识别算法对3个部分分别进行匹配识别运算。即 $P_j=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{{(e_{\mathrm{ji}}-E_{\mathrm{ji}})}^2}{D_{\mathrm{ji}}},$ j＝A，B，C。

步骤4)：用不同的系数加权3个识别结果

由步骤3)可以得到P_A、P_B、P_C。对每部分的识别结果乘以加权系数即得到了最终的识别结果P。即P＝a*P_A+b*P_B+c*P_C。通常加权系数根据虹膜纹理分布的情况则加权系数a∶b∶c＝7∶2∶1。在本发明的具体实现过程中，取a＝0.7，b＝0.2，c＝0.1。

步骤5)：设定阈值T进行识别判决

设定一个阈值T，当P当P<T时，即判定待识别的虹膜图像与数据库中的当前虹膜图像样本来自同一只眼睛。当P>T时，即判定待识别的虹膜图像与数据库中的当前虹膜图像样本来自不同的眼睛。

通过以上10个步骤，就可以通过采集虹膜图像来进行实时性的身份识别。

需要说明的是：

1、步骤3的步骤1)中引入安全系数的目的是为了避免在对虹膜图像二值化过程中引入眼睫毛等噪声干扰。安全系数的值不能设的过大或过小，通常取值为5较适宜。

2、步骤3的步骤2)中找到粗略圆心的目的是为了在步骤3的步骤4)中能以此圆心为坐标原点建立坐标系，进而确定四个可以搜索圆的象限。

3、步骤3的步骤3)在众多的边缘检测算子中，只选用Roberts算子进行边缘检测。是由于二值化的阈值选择较好，使得内变缘已经较明显。故使用简单的Roberts算子已经能精确地检测出内边缘，且由于其计算简单，所以也具有很高的运算速度。

4、在步骤3的步骤2)中按照灰度的累积计算已可以粗略的估计出虹膜内边缘的圆心和半径。但还要进行步骤3)、步骤4)、步骤5)的原因是由于虹膜的纹理大部分都在靠近内边缘的地方，所以粗略的估计会影响识别率。虹膜的内边缘需要精确定位，故通过步骤3)、步骤4)、步骤5)来精确的定位内边缘。

5、步骤4中定位虹膜外边缘的过程之所以与步骤3定位内边缘的方法不同是由于虹膜外边缘和巩膜灰度差值不是很大，不能用二值化的方法来快速的定位。并且由于靠近外边缘的虹膜区域纹理分布较少，还容易受到眼睫毛的干扰，所以外边缘的定位不需要很高的精确度。采用圆形边缘检测器进行来定位外边缘已经达到了定位的精度要求。

6、之所以要通过步骤5做归一化处理是由于每次采集虹膜图像时虹膜的位置都会有所不同。而且由于采集系统光照的影响，会导致瞳孔的放大或收缩，这会使虹膜的大小也跟着变化。所以定位出圆环形的虹膜区域不能直接用来特征提取，必须把圆环形虹膜转换成一个大小固定的灰度矩阵图像。虽然每次采集虹膜图像时，都会造成虹膜的绝对位置的改变，但是虹膜纹理的相对位置一般是不会变化的。故我们采用极坐标变换的方法对虹膜进行归一化。由于虹膜的内外边缘通常不是同心的，所以这种极坐标变换也不是同心的。

7、步骤5中的步骤3)是对归一化后的灰度图像通过直方图均衡化来达到图像增强的目的。该步骤能解决图像采集过程中由于非均匀光照引起的虹膜图像灰度分布不均的问题。

8、步骤6为对归一化后虹膜图像灰度矩阵做第一次图像分割处理。它按照虹膜纹理分布的稀疏把归一化灰度矩阵分成A、B、C3个部分。

9、步骤7为第二次图像分割，它在第一次图像分割的基础上进而把虹膜图像分成18个受噪声干扰可能性较小的区域。

10、步骤8中，每个分割区域经过3阶二维小波变换后都能得到10个小波通道，之所以要舍弃3个包含对角细节的小波通道是由于这3个小波通道受噪声干扰的可能性最大，会影响识别率，所以舍弃这3个小波通道。

11、步骤9中，均值用来衡量小波通道的能量，方差用来衡量小波通道中小波系数偏离均值的程度。

12、步骤10的步骤2)中，在系统的识别阶段只需要提取一半特征值是由于在匹配算法中采用的是方差倒数加权求和的匹配识别算法。并且所用的方差是数据库中虹膜样本特征中的方差，而不需要知道待识别虹膜图像中的方差特征。所以在识别阶段，提取待识别虹膜图像的虹膜特征时，只需提取待识别虹膜图像经过二维小波变换后小波系数的均值特征，而不需要提取小波通道的方差特征。这在一定程度上提高了虹膜识别的运算速度，提高了识别系统的实时性。

13、步骤10的步骤5)中，可以根据安全需要的不同来设置不同的阈值T。在对安全性要求较高的场合，可以相对设置较低的阈值T；在对安全性要求较低的场合，可以相对设置较高的阈值T。

本发明把虹膜定位分成内边缘定位和外边缘定位两个部分，并把重点放在内边缘定位上。接着采用直角坐标与极坐标映射理论把定位后的图像归一化到一个固定的灰度矩阵中。并对归一化后的图像进行两次图像分割，最终分割成18个小区域。再通过对归一化处理后虹膜图像的各个小分割区域进行二维小波变换，提取出主要小波通道的小波系数均值和方差作为特征值。在匹配识别算法中，对归一化处理中第一次图像分割后的3部分区域分别采用方差倒数加权求和的匹配算法进行识别判决，得到3个识别结果。最后用不同的可信度系数来加权这3个识别结果，得到最终的识别结果。采用本发明的虹膜识别算法进行身份识别可以获得较高的识别准确率和较好的抗噪性能，并有较快运行速度，满足系统实时性的要求。

本发明的创新之处在于：

1、在对虹膜图像定位过程中，按照先定位内边缘再定位外边缘的顺序，先对虹膜图像二值化，接着用Roberts算子进行边缘检测，然后以Hough变换的思想来精确定位虹膜内边缘；在定位外边缘时，采用Daugman提出的圆边缘检测模板来定位外边缘，并根据己定位的内边缘的参数，缩小了圆边缘检测模板的搜索范围，提高了运算速度。

2、对归一化矩阵进行两次图像分割。第一次分割把归一化矩阵分割成3个部分。第二次分割把归一化矩阵分割成18个小区域。经过两次分割较大限度的避免了眼皮和眼睫毛的影响，并为后面的特征提取做了充分的准备。接着对归一化后的灰度矩阵图像做了直方图均衡化的图像增强处理，避免了由于光照不均而引起图像强暗差别较大的情况。

3、在特征提取阶段，以Haar小波为小波基对18个小分割区域分别进行3阶二维小波变换得到10个小波通道。在这10个小波通道中，舍弃掉含有噪声干扰的包含水平高频和垂直高频信息的3个小波通道，保留剩余的7个通道。并对这7个通道的二维小波系数求均值和方差作为特征值。

4、采用方差倒数加权求和的匹配算法。并且所用的方差为数据库中已有的方差，使得在识别阶段只需提取待识别虹膜的小波系数的均值作为特征值，即识别系统识别阶段提取的特征值数量仅为学习阶段提取特征值数量的一半。这在一定程度上较好的加快了系统的运算速度，提高了系统的实时性能。

附图说明

图1是虹膜定位结果示意图

该图虹膜原图像使用的是CASIA虹膜数据库(version 1.0)中的一幅虹膜图像。它是一幅高为280个象素，宽为320个象素的灰度图像矩阵。图中两个白色的圆形包含的区域即为环行虹膜区域。

图2为坐标系变换模型示意图

图3为虹膜归一化处理结果

把图1中的环行虹膜区域归一化成一个固定的灰度矩阵，图中灰度矩阵的高为64个象素，宽为1024个象素。

图4为虹膜图像第一次分割

按照虹膜的纹理分布情况把归一化灰度矩阵分为3个区域。其中A区域纹理分布最密，整个虹膜70％的纹理都分布在该区域；B区域包含有20％的虹膜纹理；C区域包含有10％的虹膜纹理。

图5为虹膜图像第二次分割

按照受噪声干扰的可能性，在图4的基础上，继续把图像分割成18个小分割区域(图中(1)至(18))。这18个小区域受眼皮和眼睫毛的干扰较小。.

图6为二维变换小波通道选择示意图

技术方案中的步骤8对步骤7中的每个小分割区域以Haar小波为小波基进行3阶二维小波变换，二维小波变换后共得到10个小波通道，分别记为LL³、LH³、HL³、HH³、LH²、HL²、HH²、LH¹、HL¹、HH¹；舍弃HH¹、HH²、HH³这三个通道只保留剩余的7个通道。其中，LL³为3阶二维小波变换后虹膜图像在水平低频和垂直低频下的信息。LH¹、LH²、LH³分别为1、2、3阶小波变换后虹膜图像在水平低频和垂直高频下的信息。HL¹、HL²、HL³分别为1、2、3阶小波变换后虹膜图像在水平高频和垂直低频下的信息。HH¹、HH²、HH³分别为1、2、3阶小波变换后虹膜图像在水平高频和垂直高频下的信息。

图7为本发明流程示意图。

具体实施方式

采用本发明的算法，在CASIA虹膜数据库(version 1.0)中进行了实验。我们随机选取CASIA数据库中100组虹膜图像。每组虹膜图像中取四幅图像，共400幅虹膜图像进行实验。在学习阶段，对每组虹膜图像中的4幅虹膜图像分别提取252个均值和方差特征值，再对这4幅图像的特征值取均值作为最后的对应于该组虹膜图像的样本特征值存入样本数据库中。同理，我们把这100组虹膜图像的样本特征值都提取出来，并存入样本数据库中。在识别阶段，我们对样本数据库中每组样本特征值都用400幅虹膜图像来进行模式匹配识别运算。这样总共进行40000(400×100)次模式匹配识别运算。共得到40000个匹配识别计算结果p。根据p的经验值来设置阈值T。当T＝2.0时能取得较好的正确识别率，此时正确识别率为98.6％。

Claims (5)

1.基于图像分割和二维小波变换的虹膜识别方法，其特征是，它包含下列步骤：

步骤1：采集虹膜图像

通过虹膜图像的采集装置采集虹膜图像，得到可用于进一步处理的虹膜图像灰度矩阵H(x，y)；

步骤2：虹膜图像中值滤波

对步骤1所得的虹膜图像灰度矩阵H(x，y)做平滑处理得到灰度矩阵I(x，y)；

步骤3：虹膜内边缘定位，具体包括以下步骤：

步骤1)：图像二值化

求出灰度矩阵I(x，y)的灰度直方图，找出灰度直方图内中灰度值在(20～125)范围内的峰值所对应的灰度值M，把灰度值M与一安全系数D相加得到可用于灰度图像二值化的阈值Y；以Y为阈值对虹膜灰度图像I(x，y)进行二值化，得到二值化图像B(x，y)；

步骤2)：找内边缘粗略圆心

计算出二值化矩阵B在x方向的灰度投影量P(x)，在y方向的灰度投影量P(y)；在一维数组P(x)中，找到其最小值，并找到该最小值对应的x₁；同理在一维数组P(y)中，找到其最小值，并找到该最小值对应的y₁；(x₁，y₁)即为虹膜内边缘的粗略圆心；

步骤3)：边缘检测

对二值化后的图像B(x，y)用Roberts算子进行图像的边缘检测，得到一个含有边缘的二值图像BW(x，y)；

步骤4)：将边缘点分成4个象限

在BW(x，y)中以步骤3)中得到的内边缘粗略圆心(x₁，y₁)为坐标原点建立坐标系，把离散边缘点分为四个象限；在每个象限内，在距离原点30个像素和50个像素的扇形范围内，随机选取3个像素相差在10个单位以上的边缘点；同理，在其他象限内也是如此选择3个像素点，这样4个象限共有12个像素点；

步骤5)：4个象限联合精确定位

在步骤4)的12个点中选出3个不在同一直线上的点，由这个3个点构成一个圆，求出其余9个点到这个圆的距离d；12个点最多能构成220个圆，也即有220个距离d，在这220个距离中找到最小的d对应的那个圆即为虹膜的内边缘，设得到的虹膜内边缘的圆心为(x_a，y_a)，半径为r_a；

步骤4：虹膜外边缘定位，具体包括以下步骤：

步骤1)：限定圆边缘检测模板的迭代范围

在步骤2得到I(x，y)中，用圆形边缘检测器进行迭代求灰度积分值；在迭代过程中把内边缘的圆心(x_a，y_a)作为迭代(x₀，y₀)的初始值，并把(x₀，y₀)的搜索范围限定在以(x_a-10，y_a-10)、(x_a-10，y_a+10)、(x_a+10，y_a-10)、(x_a+10，y_a+10)为顶点的矩形内，并把r的搜索范围限定在70～110个单位像素内；在搜索的过程中，并不是在整个圆上对灰度做圆积分，而是在以(x₀，y₀)建立的直角坐标系中，角度为-45°～45°、135°～225°的圆形弧线上求灰度的积分；

步骤2)：迭代中找出外边缘

在步骤1)的范围内迭代参数空间的(r，x₀，y₀)的值，求出灰度积分变化最大的那个圆即为虹膜的外边缘，对应的参数值(r，x₀，y₀)即为虹膜外边缘的半径和圆心值(r_b，x_b，y_b)；

步骤5：虹膜图像归一化，具体包括以下步骤：

步骤1)：建立坐标转换模型、并将虹膜图像归一化

由步骤3和步骤4得到虹膜内、外边缘的圆周参数(r_a，x_a，y_a)、(r_b，x_b，y_b)，以内圆圆心(x_a，y_a)作为坐标系的原点建立将直角坐标转换为极坐标的数学模型；在该模型中从原点开始做与水平线成θ角的射线，它与内、外边界各有一个交点，分别记作B(x_i，y_i)，A(x_o，y_o)；射线上两个交点A，B之间的任何一点的坐标(x，y)都可以用A(x_o(θ)，y_o(θ))，B(x_i(θ)，y_i(θ))的线性组合来表示： $\left\{\begin{array}{c}x(r,\mathrm{\&theta;})=(1-r)*x_i\left(\mathrm{\&theta;}\right)+r*x_o\left(\mathrm{\&theta;}\right)\\ y(r,\mathrm{\&theta;})=(1-r)*y_i\left(\mathrm{\&theta;}\right)+r*y_o\left(\mathrm{\&theta;}\right)\end{array}\right.,$ 这样就可以将虹膜图像归一化为以θ为横轴，r为纵轴的固定尺寸的灰度矩阵P(x，y)；

步骤2)：直方图均衡化

对步骤1)得到的灰度矩阵P(x，y)做直方图均衡化处理得到归一化灰度矩阵PI(x，y)；

步骤6：归一化后灰度矩阵第一次图像分割

在PI(x，y)中，选取灰度矩阵上方16×1024的部分作为虹膜纹理区域(A)；选取行宽为17～48，列的范围分别是1～128、384～640、896～1024共3小块区域作为虹膜纹理区域(B)；选取行宽为49～64，列的范围分别是1～64、448～576、980～1024共3小块区域作为虹膜纹理区域(C)；

步骤7：归一化后灰度矩阵第二次图像分割。

把步骤6中的纹理区域(A)分成8个小块的灰度矩阵：行宽为1～16，列宽分别为1～128、129～256、257～384、385～512、513～640、641～768、769～896、897～1024；把虹膜纹理区域(B)也分成8个小块的灰度矩阵：行宽为17～32，列宽分别为1～128、385～512、513～640、897～1024；行宽为33～48，列宽分别为1～128、385～512、513～640、897～1024；把虹膜纹理区域(C)分成2个小块的灰度矩阵：行宽为49～64，列宽为1～64、961～1024的两个灰度矩阵区域结合在一起；行宽为49～64，列宽为449～578的灰度矩阵区域；

步骤8：对每个分割区域进行二维小波变换

对步骤7中的每个小分割区域以Haar小波为小波基进行3阶二维小波变换，二维小波变换后共得到10个小波通道，把这些通道分别记为LL³、LH³、HL³、HH³、LH²、HL²、HH²、LH¹、HL¹、HH¹；舍弃HH¹、HH²、HH³这三个通道只保留剩余的7个通道；

步骤9：提取小波系数均值和方差作为特征值

对某个小分割区域二维小波变换后的每个小波通道提取其均值和样本方差：

$E_n=\frac{1}{M\&times;N}\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\left|x(i,j)\right|,$ $D_n=\frac{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\left[\right|x(i,j)|-E_n{]}^2}{M\&times;N-1}$ 作为特征值，这样该小分割区域每个小波通道能提取出两个特征值，7个小波通道就能提取出14个特征值；对每个小分割区域都重复此过程，共能提取出252个特征值；

步骤10：用方差倒数加权均值差进行匹配识别，具体包括以下步骤：

步骤1)：提取完整特征做为虹膜特征样本

在系统的学习阶段，对要学习的每幅虹膜图像都按照步骤2到步骤9的处理过程进行处理，每幅虹膜图像能提取出252个特征值，把这些特征值作为该虹膜图像的样本存储在虹膜样本数据库中，用于后面识别判定；

步骤2)：提取均值特征用于识别

在系统的识别阶段，对一幅未知的虹膜图像，对该图像按照步骤2到步骤9的处理过程进行处理，并在步骤9中只提取128个均值E_n用于识别；

步骤3)：分3个部分分别进行匹配识别

在匹配识别过程中，把样本中的252个特征值按照步骤6的三个纹理区域分成3个部分：E_A、E_B、E_C：E_A为虹膜纹理区域(A)经变换提取的特征值，E_B为虹膜纹理区域(B)经变换提取的特征值，E_C为虹膜纹理区域(C)经变换提取的特征值；同理也把待识别虹膜图像的128个特征值分成3个部分：e_A、e_B、e_C；按照方差倒数加权求和的匹配识别算法对3个部分分别进行匹配识别运算，即 $P_j=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{{(e_{\mathrm{ji}}-E_{\mathrm{ji}})}^2}{D_{\mathrm{ji}}},j=A,B,C;$

步骤4)：用不同的系数加权3个识别结果

由步骤3)可以得到P_A、P_B、P_C，对每部分的识别结果乘以加权系数即得到了最终的识别结果P，即P＝a*P_A+b*P_B+c*P_C；

步骤5)：设定阈值T进行识别判决

设定一个阈值T，当P当P＜T时，即判定待识别的虹膜图像与数据库中的当前虹膜图像样本来自同一只眼睛；当P＞T时，即判定待识别的虹膜图像与数据库中的当前虹膜图像样本来自不同的眼睛。

2.根据权利要求1所述的基于图像分割和二维小波变换的虹膜识别方法，其特征是，所述步骤2中对步骤1所得的虹膜图像灰度矩阵H(x，y)做平滑处理过程中，滤波器为采样窗口大小为9个离散象素点的非线性中值滤波器。

3.根据权利要求1所述的基于图像分割和二维小波变换的虹膜识别方法，其特征是，步骤3中的步骤1)图像二值化处理过程中，所述安全系数D的取值范围为：3≤D≤7。

4.根据权利要求3所述的基于图像分割和二维小波变换的虹膜识别方法，其特征是，所述安全系数D为5。

5.根据权利要求1所述的基于图像分割和二维小波变换的虹膜识别方法，其特征是，步骤10的步骤4)中所述加权系数具体比例是：a∶b∶c＝7∶2∶1。

Images