双波长视网膜血氧测量系统

先永利，戴 云，高椿明，杜 睿



双波长视网膜血氧测量系统

先永利1，2，3，4，戴 云1，3，高椿明2，杜 睿1，3

( 1. 中国科学院自适应光学重点实验室，成都 610209； 2. 电子科技大学 光电信息学院，成都 610054；3. 中国科学院光电技术研究所，成都 610209；4. 中国科学院大学，北京 100049 )

为了测量视网膜血氧饱和度，本文设计一种以商业眼底相机为平台的双波长视网膜血氧测量系统。该系统将商业眼底相机输出的视网膜像进行二次成像，同时采集含氧血红蛋白和还原血红蛋白的等吸收波段和非等吸收波段图像。对获取的双波长图像进行图像配准、血管分割、光密度比计算等一系列运算，最终完成视网膜血管血氧饱和度的测量。对健康人眼进行重复性测量，动脉和静脉的血氧饱和度平均值分别是92.85%和56.75%，标准差的平均值分别为2.41%和3.63%，表明系统的准确性和重复性较好。本系统在获取视网膜结构图像的同时获得视网膜血管功能信息，可以为生命科学相关研究、眼底相关疾病的诊断提供有力工具。

眼底相机；视网膜；血氧；光密度比

0 引 言

视网膜血管状态是眼科疾病临床诊断和治疗的重要依据之一，许多视网膜病变，如糖尿病视网膜病[1-2]、青光眼[3-5]、血管阻塞[6-7]等，均会引起血液中氧气的大量消耗，则相应血氧饱和度会发生变化。此外，视网膜组织和血管也是唯一能直接观测到的人体深部结构，是了解某些全身疾病的重要窗口，比如视网膜血管的变化和高血压以及其他心血管疾病之间存在必然的联系[8]。因此，通过对视网膜图像进行处理和分析、获取其结构和血氧饱和度功能信息对疾病的早期诊断和检测有着非常重要的意义。

眼底视网膜血氧无损测量技术最早可追溯至20世纪60年代，Hickam和Frayser首次通过照相技术获取人眼的双波长图像，测量视网膜血氧值[9]。近年来，研究人员一直致力于研究更为可靠的血氧测量技术，大部分方法都是基于含氧血红蛋白(HbO2)和还原血红蛋白(Hb)对不同波长的光吸收差异，采集多个波长的眼底图像，通过图像处理计算得到血氧值[10]。国外已有学者对其技术基础、光路设计、应用价值等进行了一定的研究。1999年，J. Beach等人以商用眼底相机为平台，采用自制的分光装置，把两个波长的眼底图像同时成像在CCD的不同位置[11]。2006年，S. Hardarson等人采用商用的分光装置，使四个波长的眼底图像同时成像在同一个CCD的不同位置[12]。2010年，D.Nakamura等人利用分光镜和特制的滤光片将两个波长的眼底图像成像在两个CCD上[13]。尽管利用上述多波长视网膜成像系统可以对视网膜血管的血氧饱和度进行功能评价，但是诸如人眼运动、成像光路的复杂性、视网膜色素沉积、CCD相机的响应差异等问题，限制了血氧测量技术的发展。近年来，一些非常规的多光谱成像技术也被应用于眼底血氧检测的研究。超光谱成像技术因其可以提供更丰富的视网膜信息，受到研究人员的重视。2004年，Khoobehi利用超光谱成像技术研究视网膜血氧，该技术使用计算机生成的二维全息光栅，产生一系列一阶和二阶衍射图像，但计算这些衍射图像很费力，而且可能依赖于迭代过程的初始条件的选择，同时，由于其成像时间长，需要眼睛固定不动，目前还只能用于动物实验[14]。

国内眼底相机技术研究及产业化起步较晚，目前商业化的眼底相机也仅用于常规眼底视网膜成像检查及基本结构分析，对于视网膜血氧测定等功能性分析暂无相关研究报道。本文提出一种基于商业眼底相机的双波长视网膜血氧测量系统，该系统由视网膜双波长图像采集和血氧计算两部分组成。视网膜图像采集子系统将眼底相机输出的视网膜像一分为二，通过特定波长的滤光片获得视网膜双波长图像，并分别成像在两个CCD上。双波长图像经数据传输至计算机，由血氧计算子系统对其处理。血氧计算子系统对双波长图像进行配准、血管分割、光密度比计算，最终完成血氧值的测量。实验表明，系统可同时实现对眼底视网膜结构和血管血氧饱和度功能成像，为人眼血氧研究和眼底相关疾病的诊断提供有力工具。

1 基于眼底相机的双波长视网膜血氧测量系统

1.1 血氧饱和度测量原理

视网膜血氧饱和度的测量基于光谱法，它是一种结合组织吸收光谱特性和分光光度计技术的无损测量方法[10]。血氧测量系统需要一个对血氧饱和度不敏感的光和一个对血氧饱和度敏感的光，采集这两种波长光的眼底图像，利用图像处理计算得到血氧值。波长的选择依据见图1，图中给出了含氧血红蛋白(HbO2)与还原血红蛋白(Hb)在不同波长下的消光系数。结合眼底相机照明光源出射光谱，本文选择570 nm和600 nm光进行双波长成像，其中：570 nm光的HbO2和Hb消光系数近似相等，是对血氧饱和度不敏感的光；600 nm光的HbO2和Hb消光系数差异较大，是对血氧饱和度敏感的光。

图1 含氧血红蛋白(HbO2)和还原血红蛋白(Hb)在不同波长下的消光系数

血氧饱和度与血管的光密度值(Optical Density，OD，用表示)有关，OD定义：

血氧饱和度与不同波长下血管光密度比(Optical Density Ratio，ODR，用表示)近似成线性关系[12]：

需要说明的是，要测量血氧值，还需对式(2)中的常数、进行标定。本文参考了Schweitzer D等人的血氧仪标定值[15]，即动脉和静脉的血氧饱和度平均值分别为92.2%和57.9%。对10位健康人右眼进行5次重复测量，共50次，测量动脉和静脉的平均值分别为0.21和0.49，则代入式(2)可得：

由式(3)和式(4)联合可得= -1.225，=1.179。故血氧饱和度()计算公式：

1.2 双波长图像采集子系统

双波长图像采集子系统基于项目合作单位重庆康华瑞明科技股份有限公司生产的眼底相机平台，为了尽可能的减小对眼底相机产品的影响，提高系统兼容性，双波长图像采集子系统仅对其成像端进行改造，将其输出的视网膜图像进行二次成像，利用窄带干涉滤光片获取双波长视网膜图像，其光路原理图如图2所示。眼底相机获取的视网膜像经成像物镜L1成像于无穷远，构成像方远心光路。通过分光棱镜BS将光束分成两路，分别透过(570±5) nm和(600±5) nm的干涉滤光片滤出所需波长的成像光，经成像透镜(L2和L3)后同时成像在CCD1和CCD2靶面上，实现视网膜双波长同时成像。将采集到的视网膜双波长图像输入计算机由血氧计算子系统处理，计算血氧饱和度。二次成像及波长选取带来的光能损失通过扩瞳、适当加大曝光量和提高CCD灵敏度等措施弥补。图3为双波长图像采集子系统实物图及系统集成照片，图4为实际获得的一组视网膜双波长图像。

图2 双波长图像采集子系统光路图

图3 双波长图像采集子系统

图4 双波长眼底图像

从图4(a)、(b)可以看出，570 nm的眼底视网膜图像中，动静脉的光密度基本一致，表现为动脉血管和静脉血管的颜色相当。而600 nm的眼底图像中，静脉血管的光密度小于动脉血管，表现为静脉血管相对于动脉血管颜色更深。

1.3 血氧计算子系统

根据1.1原理部分的阐述，血氧饱和度的测量需要获取血管内像素段的灰度最小值，以及血管外像素段的平均灰度值，并保证双波长图像中用于计算O的血管像素点位于同一眼底位置。由于双波长成像非共光路引入了图像位置偏差，为了满足血氧测量的需求，需要将双波长图像配准到同一坐标系。因此，血氧计算子系统主要包括图像配准、血管分割、血氧计算等三个部分。

1.3.1 图像配准

将570 nm图像血管边界(红线)和600 nm图像血管边界(绿线)直接叠加到570 nm图像上，发现两者的血管边界没有重合(由于600 nm眼底图像中动脉血管与背景对比度低，血管边界提取有很大的难度，因此这里的重合是针对两者的静脉而言，下同)，见图5。这是由于双波长成像非共光路引入的图像位置偏差，因此在进行血氧计算之前必须先对两幅眼底图像进行配准。本文采用基于互信息的方法对双波长图像进行配准。

图5 配准前的图像

互信息描述两幅待配准图像的相似程度，两幅图像相互包含信息量越大，互信息值越高[16]。由于视网膜各种物质在不同波长光下，其消光系数不同，因此视网膜同一位置的不同波长光的视网膜图像之间存在差异。但是不同波长的光都能对眼底相同的结构成像，说明它们在统计学上是相关的。因此，不同波长光的视网膜图像可以使用基于互信息的图像配准算法。

针对眼底图像进行分析，由于眼底图像中存在大量相似血管信息，直接对两幅图像进行互信息计算不利于提高配准精度，而且由于大量无关背景信息参与计算，增大了运算量。由于视盘区域存在明显特征且包含大量互信息计算可用数据，因此通过定位视盘实现配准小区域的获取，有利于提高配准速度[17]。

基于视盘定位的关键块图像配准算法步骤如下：

1) 处理对象选择信息丰富的基准图像即570 nm图像，首先对图像进行自适应二值化处理，即类间方差最大的阈值计算；

2) 对二值化图像进行八连通区域标记，并分别求取面积；

3) 针对所有区域面积进行比较，得最大面积连通区包含视盘；

4) 针对选择连通区进行质心求取，并以其为中心对两幅待处理图像进行开窗处理，得待处理图像块，图像块长宽取图像对应长宽的1/8；

5) 分析系统采集的大量双波长图像数据，平移范围在-20 pixels~20 pixels，旋转范围在-10°~10°，算法实现时，步长均取1。根据经典互信息法对两个图像块进行计算，得配准参数，即旋转量和平移量，流程图如图6所示。

图6 基于互信息的图像配准

使用上述算法对图4进行配准，获得的配准参数为：偏移，向左平移14.413 9 pixels；偏移，向上平移8.403 0 pixels；角度旋转，顺时针旋转2.005 1°。

配准后，选取600 nm眼底图像的血管边界，叠加在570 nm眼底图像的血管边界上。可以看出600 nm血管边界(绿线)与570 nm血管边界(红线)基本重合，表明上述算法能够很好地配准双波长眼底图像，如图7所示。

图7 配准结果

1.3.2 血管分割

对配准后的图像进行视网膜血管分割，确定需要计算血氧的区域。血管分割采用基于Hessian矩阵的多尺度视网膜图像分割算法[18]。对于二维图像中的一点，Hessian矩阵定义为

由于图像二阶导数对噪声非常敏感，且眼底血管粗细不一，不适合使用单一尺度的分割方法，故在应用Hessian矩阵时引入低通滤波和多尺度来解决上述问题[19]。低通滤波一般采用高斯模糊，可以减少孤立噪声点被误判为血管的情况，实现为，即每个像素点的二阶导数依次与高斯核卷积。本文使用的高斯核函数为。尺度大小由高斯核的标准差大小决定，大则运算像素区域半径大，反之亦然。尺度小易检测小血管，但可能会引入更多噪声、导致大血管被分割；尺度大则仅能检测大血管。因此，合适的尺度因子范围非常重要，根据文献[18]可知，若一幅血管图像中被检测血管的直径范围为，则空间尺度因子的取值范围应为。

1) 输入待处理原始视网膜图像，尺度因子范围为，迭代步长。

2) 依次遍历图像中的每个像素点；

5) 将特征值带入血管相似性函数[19]，得当前像素点在该尺度下的响应；

8) 像素遍历完成，得新的血管图像。

由于1.3.1已对双波长眼底图像完成配准，即双波长图像中用于计算OD的血管像素点已位于同一眼底位置，故只需对570 nm眼底图像进行血管分割即可，利用上述算法对其进行血管分割，见图8。

图8 570 nm血管分割

1.3.3 血氧计算

针对分割后的血管，首先实现平均光密度比(OD)计算，具体为：血氧计算子系统自动搜索血管内像素段的最暗点，提取其灰度值，即式(1)中值，然后自动选择靠近血管边缘的像素段，计算其平均灰度值，即式(1)中的0。为了提高OD的计算精度，选取多个最暗点计算OD然后取平均。最后根据式(5)计算得到血管的血氧值，并用伪彩色在血管对应位置表示出来，结果如图9所示。图中红色表示血氧浓度高，向蓝色变化表示血氧浓度越来越低，可清晰地观察到视网膜血管血氧饱和度的动静脉差异。

图9 血氧计算结果

2 实验与分析

为了验证系统的可靠性和可重复性，我们对10位受试者进行眼底双波长图像采集(6名男性，4名女性，平均年龄：23.8岁，标准差：1.6岁)，每位重复测量5次，共50次。对每位受试者的5次眼底图像，以视盘中心为圆心画两个不同半径的圆，分别计算环形区域内的动脉和静脉的血氧饱和度的平均值和标准差，这样可以减少视盘对血氧计算的影响。手动确定血氧分析区域，血氧计算子系统会自动给出各段动脉和静脉的血氧值。表1给出了所有受试者进行5次重复测量的动脉和静脉血氧饱和度的平均值和标准差，以及各自的动态范围和均值。为了便于观察，图10只标注了其中一名受试者在环形区域内的一段动脉和一段静脉的血氧值。