基于矩阵数据的视网膜神经纤维层厚度与眼轴长度及年龄的相关性

李青蒨 常征 郭亚文 杨丽丹 黄丽娜 陶文艳 林晨

作者单位：1中南大学爱尔眼科学院 长沙爱尔眼科医院 410000；2深圳爱尔眼科医院 518000；3香港中文大学 999077；4深圳市人民医院 518000

目前，采用光学相干断层扫描（OCT）定量检测视网膜神经纤维层（Retinal nerve fiber layer，RNFL）厚度是区分正常眼或判断青光眼进展的重要手段。国内外已有研究证实，除了青光眼等疾病因素外，RNFL厚度还与年龄、眼轴长度（Axial length，AL）等因素存在相关性[1-2]。这些因素不可避免地对RNFL测量结果产生影响，因此临床医师很难对高度近视患者是否出现青光眼性视神经损伤下定论。目前，已经有相应的研究对环视盘RNFL（Circumpapillary retinal never fiber layer，cpRNFL）厚度与AL、年龄等因素的相关性做了研究[1-6]。这些研究的数据均是基于环视乳头多位置平均值的一维数据。但由于RNFL的厚度在空间上分布是不同的[7-8]，而且年龄、AL等因素对视网膜每个平面空间位置的影响也不明确。故以多位置平均值为计算结果的一维数据研究很难运用到目前以二维厚度图为主要诊断方式的测量结果中。本研究旨在从二维空间分布角度方面研究正常人群年龄、AL等因素对RNFL厚度的影响，通过自行设计的矩阵模型算法计算其相关系数的二维分布图来探究年龄、AL因素与RNFL厚度的相关性，以期进一步校正RNFL厚度图结果，为提高青光眼的诊断效率提供数据基础。

1 对象与方法

1.1 对象

收集2015年10月至2020年7月于长沙爱尔眼科医院门诊进行体检的正常人群（180 例，360 眼）（包括近视眼）的检查结果资料，所有研究对象均曾接受全面的眼科检查，包括：病史回顾，最佳矫正视力（BCVA）的测量，使用非接触眼压计（日本Canon公司）进行眼内压的测量，使用IOLMaster（美国Carl Zeiss公司）进行AL测量。纳入标准：①年龄18～80岁；②BCVA≥0.9；③眼压：10～21 mmHg（1 mmHg=0.133 kPa）；④前房正常，房角开放，视网膜正常；⑤OCT图像信号强度≥6，图像视盘居中清晰完整，眼底视网膜血管及RNFL连续。排除标准：①视盘外观异常（视盘出血，杯盘比（Cup disc ratio，C/D）>0.5，双眼杯盘比C/D比值相差0.2以上）或有RNFL损害；②有黄斑疾病或任何引起视神经改变的眼部或系统性疾病；③有眼内手术史或视网膜激光治疗史；④无晶状体眼或人工晶状体眼；⑤糖尿病等可能影响测量结果的其他全身性疾病；⑥青光眼家族史；⑦OCT图像信号强度<6，图像视盘偏位，算法分割失败，运动伪影。本研究遵守赫尔辛基宣言，获得爱尔眼科集团伦理委员会批准[伦理批号：长爱伦审字（2019）KYPJ013号]。

1.2 RNFL厚度矩阵数据库的建立与数据分析

1.2.1 RNFL厚度图的获取 采用Cirrus HD-OCT（美国Carl Zeiss公司）对研究对象双眼进行检查，所有研究对象均使用200×200 点阵模式进行视盘扫描，测量范围为6 mm×6 mm的区域。经过计算机分析分层后生成RNFL厚度图，所有OCT检查均在同一台设备上进行，避免产生误差。排除信号强度低于6，运动伪影，丢失数据，明显的偏心错位或算法分割失败的图像。

1.2.2 RNFL厚度矩阵数据库的建立 导出OCT系统中研究对象的RNFL厚度图。利用基于Python（Version 3.7）语言自行开发的计算程序，将厚度图转换为厚度矩阵数据。利用自行开发的分析程序对视盘重新进行中心定位，根据视盘中心进行图像对称翻转，左眼图像经对称翻转后全部变为右眼位，将受检者所有RNFL厚度图进行堆叠，形成200×200像素RNFL厚度矩阵受检者数据库。同时读取以视盘中心点为圆心，半径为1.73 mm的圆周上的平均cpRNFL数据。

1.2.3 基于矩阵数据的混合线性模型的设计及运算 将经过转换过的OCT矩阵数据整合成50×50的超级像素矩阵并导入运算库，建立与每个矩阵数据相对应的年龄、AL数据。以随机截矩-斜率的混合线性模型yloc（t）=XlocT（t）β1+ZlocT（t）bloc+eloc（t），eloc（t）～N（0，σ2）为基础（loc为位置坐标），将年龄、AL作为协变量，RNFL厚度作为因变量，将每一个矩阵同一超级像素点（RNFL厚度图有200×200个像素点，将其整合成50×50个超级像素点，1个超级像素点等于4×4个像素点）对应的厚度数据作为厚度数组，代入混合线性模型。共计运算2 500 次（50×50 次）。生成年龄协变量、AL协变量的有效斜率图。将因图像对位可能产生明显结果误差的数据矩阵四边2超级像素区域的数值设置为0。为方便大规模的数据计算，本研究建立了基于混合模型矩阵的函数包（包括输入参数：厚度数据、位置坐标、年龄、AL等，输出参数：P值、斜率等），利用循环语句进行几千次的计算。再按照对应的位置信息进行重新排列组合，将有统计学意义的区域进行颜色标记，如红色代表正相关、蓝色代表负相关，颜色深浅代表相关性的大小，无统计学意义的区域设为空白，生成与年龄、AL相关的RNFL分布变化斜率图（见图1）。在正相关区域，RNFL厚度随AL增长而增厚；在负相关区域，RNFL厚度随AL增长而变薄，采用RNFL厚度/眼轴变化（μm/mm）表示变化率。

图像峰值计算系通过索引在斜率图上分别找出正向、负向的最高及最低值，并找出其所在坐标位置。cpRNFL单参数混合线性模型采用与矩阵模型同一方法，计算1次取得结果。

1.2.4 堆叠图及角度测量 将全部受检者的RNFL图像、AL<24 mm及AL≥28 mm受检者的RNFL厚度图分别进行堆叠后生成总厚度图。以视盘中心为轴心，计算在360°辐轴上RNFL厚度堆叠图上方、下方的视网膜神经纤维（Retinal nerve fiber，RNF）束中总厚度最大值所在的角度线，测量上、下方束连线所形成的夹角。

1.3 统计学方法

系列病例研究。本研究使用STATA软件（version 14.0，stata Corp.，College StationTX）及基于Python（version 3.7）语言的程序进行数据处理。利用混合线性模型，对所有正常人群的全RNFL厚度图资料进行矩阵运算。其中，每1 个位点的RNFL厚度作为因变量，年龄、AL作为协变量，计算得到该位点的变化率及P值。将每次运算得到的变化率、P值整合成为对应50×50 的斜率图、P值图，以P<0.05为阈值，设置斜率图中P≥0.05的对应位置斜率为0，从而形成年龄协变量、AL协变量的有效斜率图。使用单混合线性模型分析年龄、AL与半径1.73 mm处的cpRNFL厚度的关系，以截距、变化率、P值描述其线性关系。

图1.RNFL厚度变化率算法流程图利用混合线性模型对所有受检者的全RNFL厚度图资料进行矩阵运算，生成与RNFL厚度图相对应的RNFL变薄速率、P值图Figure 1.Diagram of steps of the algorithm for the rate of change in RNFL thickness.A mixed linear model is used to perform matrix operations on the data from the full RNFL thickness map for all subjects.Then the rate of RNFL thinning and P value maps corresponding to the RNFL thickness map were generated.RNFL,retinal nerve fiber layer.

2 结果

2.1 一般资料

由于图像视盘偏位、算法分割失败等原因，有147眼被排除，最终123例（213眼）纳入研究，其中男31例（55眼），女92例（158眼）；RNFL厚度图213张。123例研究对象年龄为18～69（30.9±13.0）岁，SE为-14～+5（-6.70±3.12）D，AL为22.07～30.00（26.10±1.42）mm。

2.2 cpRNFL的厚度与年龄、AL的相关性

整体RNFL的厚度值为（105.1±9.5）μm；年龄与cpRNFL厚度呈负相关（r=-0.146，P=0.023），在1.73 mm半径处，全周cpRNFL厚度随年龄的变化率为-0.146 μm/年。AL与cpRNFL厚度呈负相关（r=-1.012，P=0.026），全周cpRNFL厚度随AL的变化率为-1.012 μm/mm。

2.3 RNFL厚度与AL、屈光度的相关性及空间分布

所有受检者的RNFL厚度图经过左右方向翻转对位后叠加形成全RNFL厚度堆叠图，在视网膜各区域RNFL厚度分布不均，上、下方RNFL厚度较鼻侧、颞侧厚（见图2）。以视盘为中心，距离视盘越近则该区域的RNFL越厚。图3 为RNFL厚度随AL、屈光度变化的有效斜率图，图中存在2 种相关性相反的区域。在颞上、颞下方，RNFL厚度与AL呈负相关（蓝色区域，P<0.05），即随着AL增长，该区域的RNFL厚度变薄；但在颞下方偏颞侧区域，RNFL厚度与AL呈正相关（红色区域，P<0.05），即随着AL增长，该区域的RNFL厚度变厚（见图3A）。

在有效斜率图中，各位置的斜率（RNFL厚度/AL变化）也不一致。在正相关区域（红色区域），RNFL增厚最大及最小速率分别为6.292 μm/mm（蓝色三角，P=0.0024）和1.389 μm/mm（黑色菱形，P=0.0224）；在负相关区域（蓝色区域），最大及最小速率分别为-8.186 μm/mm（白色星形，P<0.001）和-0.155 μm/mm（黄色圆点，P=0.0496）。其中最大、最小值的位置大致分布在靠近视盘及远离视盘的方向（见图3A）。

图2.全RNFL厚度堆叠图Figure 2.Stacking map of all maps of the RNFL thickness.RNFL,retinal nerve fiber layer.

由于屈光度与AL高度相关，所以本研究也将受检者的屈光度代入模型进行计算，结果显示屈光度与RNFL厚度相关性及分布与上述利用AL计算的结果基本一致，同时存在正负相关区域：颞上、颞下2个正相关区域（红色区域），即随着屈光度增加（远视增加）而RNFL增厚，RNFL增厚最大及最小速率分别为2.302 μm/mm（P=0.009）和1.384 μm/mm（P=0.047）；而以上区域颞侧，存在着2 处负相关区域（蓝色区域），随着屈光度降低（近视增加）而RNFL变薄，在该区域RNFL变薄的最大及最小速率分别为-2.677 μm/mm（P=0.007）和-0.379 μm/mm（P=0.032）（见图3B）。

2.4 堆叠图的分析

图3.RNFL厚度随AL、屈光度变化的速率分布图A：RNFL厚度随AL变化的速率分布图，白色星号表示随AL增长RNFL变薄的速率的最大值；黄色圆圈表示最小值；蓝色三角形表示随AL增长RNFL增厚的速率的最大值；黑色菱形表示最小值。B：RNFL厚度随屈光度变化的速率分布图Figure 3.The rate of distribution of RNFL thickness changes with axial length and spherical equivalent.A:The rate distribution of RNFL thickness changes with axial length.The white asterisk indicates that the maximum value of the RNFL rate of change decreases with axial length and the yellow circle indicates the minimum.The blue triangle indicates that the maximum value of RNFL rate of change increases with axial length and the black diamond indicates the minimum.B:The rate of distribution of RNFL thickness changes with spherical equivalent.RNFL,retinal nerve fiber layer;AL,axial lengths.

将AL<24 mm及AL≥28 mm的受检者的RNFL厚度图分别进行堆叠生成堆叠图，在以视盘为圆心，各轴向RNFL总厚度结果比较中发现，AL≥28 mm的受检者颞上和颞下的RNF束之间的夹角明显比AL<24 mm的受检者小。将这2 种不同夹角的RNFL堆叠图与图3A所示的结果进行比较，结果提示图3A中随着AL变长而出现RNFL增厚的正相关区域，也位于RNF束向颞侧夹角移位后所增厚的区域（见图4，虚线框）。

2.5 RNFL厚度与年龄、AL的相关性

图4.不同AL受检者的RNFL厚度堆叠图A：AL<24 mm的受检者的全RNFL厚度堆叠图，颞上和颞下的RNFL束之间的夹角为133.5°；B：眼轴长度≥28 mm的受检者的全RNFL厚度堆叠图颞上和颞下的RNFL束之间的夹角为111.7°Figure 4.RNFL thickness stack maps of subjects with different ALs.A:The total RNFL thickness map of subjects with an axial length <24 mm,the angle between the superior temporal and inferior temporal retinal nerve fiber bundles is 133.5°.B:The total RNFL thickness map of subjects with an axial length ≥28 mm,the angle between the superior temporal and inferior temporal retinal nerve fiber bundles is 111.7°.RNFL,retinal nerve fiber layer;AL,axial lengths.

将AL和年龄同时作为混合线性模型的协变量，生成AL与年龄的有效斜率图。在图5中，正负相关的区域类似，其中RNFL厚度随年龄丢失速率的最大值及最小值均位于颞下象限，分别为-1.4714 μm/年（图5A红点处，P=0.04385）和-0.1438 μm/年（图5B绿点处，P<0.001）。考虑年龄变化率远小于AL变化率，部分超长AL个例可能作为混淆因素对年龄的结果产生影响，导致误差。研究中将AL>28 mm的16 眼排除，将AL≤28 mm的197 眼重新进行基于二维矩阵的混合线性模型计算，结果见图5C—D。可见AL的有效斜率图保持稳定（图5Dvs.图5B），而年龄的有效斜率图发生明显变化，影响区域明显变小（图5Avs.图5C）。提示在同时考虑年龄和AL这2种协变量时，年龄在空间上对RNFL厚度影响微小。

3 讨论

图5.全体受检者与AL≤28 mm受检者的眼轴、年龄有效斜率对比图A：所有研究对象（213 眼）的RNFL厚度与年龄相关性的有效斜率图；B：所有研究对象（213 眼）的RNFL厚度与AL相关性的有效斜率图；C：剔除掉AL>28 mm的16眼，将AL≤28 mm的197眼纳入矩阵模型后计算所得RNFL厚度与年龄相关性的有效斜率图；D：剔除掉AL>28 mm的16眼，将AL≤28 mm的197眼纳入矩阵模型后计算所得RNFL厚度与AL相关性的有效斜率图Figure 5.Comparison of the effective slope of AL and age between all subjects and subjects with AL ≤28 mm.A:The effective slope map of the correlation between RNFL thickness and age of all subjects (213 eyes).B:The effective slope map of the correlation between RNFL thickness and axial length of all subjects (213 eyes).C:The effective slope map of the correlation between RNFL thickness and age after 16 eyes with AL >28 mm were excluded and 197 eyes with AL ≤28 mm were included in the matrix model.D:The effective slope map of the correlation between RNFL thickness and AL after 16 eyes with AL >28 mm were excluded and 197 eyes with AL ≤28 mm were included in the matrix model.RNFL,retinal nerve fiber layer;AL,axial lengths.

OCT是临床上诊断青光眼的重要工具，它能定量监测视盘及RNFL形态的细微变化。临床上主要利用该设备分析RNFL厚度改变，进而预估青光眼的进展程度。目前现有的OCT数据库是涵盖各年龄段、正常AL的健康人群RNFL厚度的标准数据库，当OCT将采集的RNFL厚度图与标准数据库进行自动对比分析时，并未对AL、年龄等影响因素同时进行自动校正，因此在使用正常人群RNFL厚度数据库对青光眼进行诊断及随访时，如果不考虑RNFL厚度与年龄、AL等因素的相关性，可能会造成青光眼的误诊。本研究纳入AL、年龄因素，分别采用平均cpRNFL厚度一维数据及RNFL厚度二维数据分析正常人群RNFL厚度与AL、年龄的相关性，以期将非疾病相关性的RNFL变薄与疾病相关性的RNFL变薄区分开来，从而提高青光眼的检测效率。

借助自行开发的程序，本研究针对全厚度矩阵数据进行了大规模的运算，研究发现全RNFL厚度图各位点的RNFL厚度数值与AL、年龄的相关性有所不同，最终形成的有效斜率图提示RNFL厚度随年龄、AL的变化斜率在视网膜的空间分布不同。在相关性趋势相同的区域，其变化斜率值的分布也不是一致的。总体来说，距离视盘越近，RNFL厚度随AL、年龄变化的变化率越明显；距离视盘越远，变化率越低平。此外，本研究还发现在颞下方靠近颞侧处的RNFL厚度与AL增加呈正相关。

本研究中，由于在眼轴的有效斜率图上同时存在正、负2种截然相反的相关区域，这一结果与预计的认知有所不同，所以我们引进了堆叠图来探讨这一现象出现的原因，通过分析RNF束的夹角变化，我们认为正相关区域可能与Leung等[9]学者提出的颞上和颞下的RNF束随AL的增加而逐渐向颞侧区域聚拢有关。但是，颞上区域未出现RNFL厚度随AL增长而增厚的现象，也许是由于颞上区域的数值变化较不均衡，所以差异未在模型结果中出现明显的统计学意义，从而导致AL有效斜率图中未出现相应标记。另外，在同时引入年龄和AL作为协变量的结果中，年龄的有效斜率图分布与AL过于相似，我们考虑可能与纳入少部分超长眼轴导致的模型误差有关，在剔除小量超长眼轴样本（7.5%，16眼）后，发现年龄协变量的影响并不明显，进一步剔除的结果也与此类似。这一结果提示在考虑AL变化的影响时，年龄因素的影响较为微小，或许可以忽略年龄因素以简化模型。

关于RNFL厚度随年龄增加而变薄这一结论已被国内外多个研究所证实[2，4，10-11]。但关于近视对RNFL厚度影响的研究存在着相互矛盾的结论。有学者研究提示RNFL厚度与AL呈负相关，随着近视患者AL的增长，RNFL厚度变薄[12-14]。相反，部分研究并没有发现AL与RNFL厚度之间的相关性[15-16]，也许与这些研究使用分辨率及灵敏度较低的一代OCT及样本量较少有关。目前关于近视患者RNFL变薄的机制存在多种说法，AL增长可导致脉络膜及视网膜的机械拉伸[17-18]，也同样影响视网膜小动静脉的管腔及视网膜血流量[19-20]。这些都可能是导致近视患者RNFL较薄和空间分布改变的原因。

本研究的优势：采用了计算机图像识别处理的方法，使RNFL厚度图各像素点的数据信息可视化，便于充分利用RNFL厚度图各位点的厚度数据进行大规模的矩阵运算，针对近视人群RNFL厚度变化的分析，与cpRNFL一维数据研究对比，具有较高的特异性及准确性。本研究的局限性：系单中心横断面研究，其中针对年龄因素的研究结果较其他长期随访的队列研究而言，可能存在一定偏倚；由于在Python程序中缺乏可以矫正双眼纳入因素的Shared frailty模型，部分源于双眼的数据可能导致轻微结果偏倚；而且本研究受检者年龄跨度较小，中青年居多。屈光状态为近视的受检者数量占比较大，影响因素的分布不均可能导致结果存在一定偏倚，未来研究可通过收集大量不同年龄及AL受检者的RNFL厚度图数据对本研究算法模型做进一步验证。

综上所述，本研究同时纳入年龄和AL因素，使用cpRNFL厚度一维数据和RNFL厚度图二维数据进行基于矩阵的混合模型计算，从空间分布特征角度研究正常人RNFL厚度与AL、年龄的相关性。将这2 种研究方法进行对比，cpRNFL厚度一维数据研究结果提示RNFL厚度与AL、年龄的相关性为单一负相关，而RNFL厚度图二维数据研究发现年龄、AL对视网膜不同区域RNFL厚度的影响并不均一，相比较之下，AL的影响甚大。在同一相关区域，AL、年龄的斜率值在不同位点有所不同。由于AL增长导致视网膜神经纤维束在颞侧移位，可能使AL因素的有效斜率图会同时存在正、负2 种相关区域。

本研究的发现给临床医师带来重要提示：针对近视人群，尤其是对高度近视人群的RNFL厚度变化分析应进行谨慎解读。由于目前现有的OCT数据库是来源于轻度近视或者正视眼人群，OCT设备软件所提供的指导RNFL厚度变化的评估算法是有所偏差的，因为它不涵盖调整眼轴长度和屈光不正的算法，在进行高度近视人群的RNFL厚度变化的检测时出现了较高的假阳性率。本研究在前期研究基础上进行了更深入的研究，纳入AL、年龄干预因素的空间特征，对二者同时进行矫正，生成AL、年龄矫正因素分布图，其呈现了干预因素所导致的RNFL变薄的分布特征。在临床应用中，该空间因素算法分别可以进一步引入偏差图算法中，生成更具有特异性的精细矫正偏差图，指导临床青光眼的判别，可在一定程度上进一步降低高度近视患者中RNFL变化的假阳性率，从而为临床医师提供更为准确的检查报告。

利益冲突申明本研究无任何利益冲突

作者贡献声明李青蒨：收集数据；参与选题及资料的分析与解释；撰写论文。常征：既往病例及图像的提供。郭亚文：提供技术支持及统计分析支持。杨丽丹：参与资料的分析及解释。黄丽娜、陶文艳：对文章的知识性内容作批判性审阅。林晨：参与选题，设计及资料的分析及解释，提供技术支持，修改论文中关键性结果、结论