视网膜损伤建模用于阈值预测与概率风险评估

摘要

本文提出了一种针对激光诱导人眼热损伤的建模方法。该方法旨在通过提供一种计算特定辐照条件下损伤程度的手段，改进对激光相关危害的评估。

摘要

随着激光在医学、国防和技术领域应用的不断增加，人眼对激光源的有意或无意暴露已成为一个主要问题。基于模型预测视网膜损伤阈值，可实现更具体场景下的激光安全评估，特别是针对实验数据未涵盖的激光参数。理想情况下，此类模型可根据波长、脉冲持续时间和光斑形状等因素，计算出ED50（损伤概率为50%时的有效剂量）值。这需要对所有损伤机制进行深入理解与建模，以反映关键参数与主导损伤机制之间的依赖关系。

本文讨论了该方法的现状（本文中已针对热效应条件，或简称为"热损伤条件"进行了验证）；重点指出了可能阻碍其成功应用的关键问题，并概述了其潜在优势。这些优势包括提高眼安全标准中激光暴露限值的准确性、优化视网膜激光手术中的剂量学，以及对户外环境中激光使用进行概率性风险评估。

引言

本研究描述了一种人体眼球生理学细节详尽的热损伤模型的建立与验证过程。在热损伤机制范围内，该模型采用阿伦尼乌斯（Arrhenius）损伤公式预测视网膜温度的演变过程及损伤阈值（参见光热损伤部分的解释说明）。典型应用涵盖所有需要预测眼部温度及损伤情况的场景，例如：评估扫描式视网膜辐照的损伤阈值，理解脉冲序列累加效应对损伤阈值的影响，以及将计算所得的阈值与激光安全标准中的安全限值进行比较。在热机制范围之外，本文还介绍了当前正在研究的建模方法，并提出了将该框架扩展至其他损伤机制的发展路线图。

本文所述工作涉及激光辐照引起视网膜损伤的建模与预测。理论上，虽然可以通过使用与人类视网膜相似的动物视网膜进行实验来确定临界剂量，但迫切需要在不开展实验的情况下实现损伤预测。激光参数的变化空间（包括波长、脉冲持续时间和重复频率）极为广泛，这意味着每一套新参数都需要进行数量大到难以承受的动物实验。此外，对于长时间辐照的情况，还必须考虑视网膜血流的影响，而这将需要进行in vivo实验。因此，对激光与眼睛相互作用进行建模似乎是唯一可行的前进方向。

深入理解损伤机制以及因此产生的损伤阈值（可作为ED50的替代指标）的需求，也与眼安全标准（IEC 60825或ANSI Z136.1）的现状密切相关。由于该标准必须涵盖全范围的波长、脉冲持续时间、重复模式和光斑尺寸，因此不可避免地引入了简化的假设、插值方法以及保守的安全系数，以应对不确定性。由于目前可用的ED50值数量有限——主要来源于非人灵长类动物研究——必须通过插值来建立全面的最大允许照射（MPE）限值。尽管这种方法具有广泛的实用性和适用性，但基于对损伤过程的详细机制理解与建模所构建的框架，在物理透明性、特定场景下的准确性以及无需详尽了解激光安全标准即可应用等方面，具有明显优势。

例如，尽管视网膜扫描会引入额外的时间和空间效应，脉冲激光器和扫描激光器仍被评估为脉冲光源。在安全限值推导过程中如何恰当地处理扫描因素，已成为过去十年中学术界持续讨论的议题。即使标准不断更新以反映技术进步，若不对新型激光系统的复杂构型及其新参数集进行简化并引入保守的安全系数，标准也难以涵盖所有情况。因此，仍存在解释上的空间，可能导致安全评估中出现不一致或错误。

基于物理的建模方法可显著减少对插值法和保守边距的依赖，并将该标准的适用性扩展到复杂场景。由于此类模型的开发与验证直接关联于对潜在损伤机制的深入理解，因此所得出的见解还可支持从现有的ED50数据中更透明且具有物理依据地推导出MPE值。

从长远来看，通过集成化的即插即用建模框架，有望简化眼安全性评估流程。该工具可输入相关系统参数，或直接导入光学设计文件（例如 Zemax 模型），而此类文件通常在产品开发过程中即可获得，从而实现一致且针对特定场景的安全性评估¹。

另一个应用领域是正在发展的高能激光（HEL）领域，例如用于反制无人机。在此领域中，主要难点在于来自目标（尤其是金属目标）的激光反射，这些反射可能快速且随机地变化，对军事人员和平民构成安全隐患^2,3。这种随机且非确定性的情况需要一种适当的评估机制——通常采用概率方法。该方法可对特定辐照情况（如辐照强度、暴露时间、波长）的发生概率进行描述，但仍需将其转化为造成损伤的可能性。本文所讨论的损伤建模方法能够填补这一空白，将概率性场景转化为损伤发生概率。

详细了解损伤机制的原理，并在软件模型中模拟这些机制，是确定损伤阈值的一种直接方法，无需进行实验。根据脉冲持续时间的不同，视网膜损伤通过不同的激光-组织相互作用机制发生（图1）^4,5,6,7：

图1：损伤机制概述。损伤机制的类型取决于暴露时间和辐照度。从短脉冲持续时间内高辐照度引起的光机械性损伤，到长时间低辐照度导致的光化学损伤，均包含在内。请点击此处查看此图的放大版本。

光机械断裂（10^-12–10^-6 s，10¹⁰–10¹⁶ W/cm²）：
光致破裂：在超短脉冲至纳秒级脉冲持续时间及极高辐照度下，非线性吸收引发光学击穿。由此形成由自由电子和离子组成的致密等离子体，迅速爆炸性膨胀并产生强烈的冲击波。该过程称为光致破裂，可在几乎不引起整体加热的情况下机械性撕裂组织，即使在吸收较弱的区域也能形成边界清晰的病灶^4,8,9,10。

等离子体诱导消融：在纳秒至微秒脉宽范围内且处于极高辐照度时，等离子体效应再次占主导地位。此时，组织去除不仅由冲击波驱动，还受到等离子体膨胀和爆炸性消融的共同作用。该机制被称为等离子体诱导消融，可导致显著的物质喷射^4,11。

光消融（≈ 10^⁻9–10⁻⁶ 秒；≈ 10⁷–10¹⁰ W/cm²）：
当脉冲持续时间在纳秒以上、微秒以下时，会发生光消融。在此过程中，分子键被辐射直接打断。该方法例如可用于通过重塑角膜来矫正眼睛的屈光不正（即所谓的 LASIK 方法）。通常，光消融操作在尚未发生等离子体形成的功率范围内进行^4,11。

热机械损伤（约 10^-9 – 10^-6 秒；约 10⁶ – 10⁸ W/cm²）：
在色素性眼组织中，黑色素体在纳秒至微秒范围内的强吸收可导致迅速过热。当黑色素体表面温度达到约 150 °C 时，微泡成核^12,13。微泡的膨胀与塌陷产生机械应力波，从而损伤视网膜色素上皮（RPE）。这种热机械机制介于光消融与热损伤之间：它并非由等离子体驱动，而是涉及与局部加热耦合的机械瞬态过程。损伤阈值取决于黑色素体的大小、形状、取向以及局部光照条件^4,11。

光热损伤（≈ 1 µs–60 s；≈ 10–10⁶ W/cm²）：
从微秒到秒量级，组织加热占主导作用。能量沉积导致温度升高，首先引起蛋白质变性，更高能量暴露下则导致凝固性坏死和碳化。文献中的时间范围存在差异：Niemz⁴ 提出为 1 µs–60 s，Zuclich¹⁴ 则为 10 µs–60 s。在辐射暴露剂量（~1–1000 J/cm²）下大致符合倒易律，但极短脉冲（热扩散受限）和极长脉冲（灌注冷却增强）情况下会出现偏离^4,6,15。

本研究中热损伤的建模基于阿伦尼乌斯积分，其定义如下：

   (1)

其中 τ 表示暴露时间，E_a 为活化能，R 为通用气体常数，T 为暴露期间的温度，A 为指前指数因子。本研究中采用的参数为 A = 1.3 × 10⁹⁹ s⁻¹ 和 E = 628 kJ/mol¹⁵。当 Ω = 1 时，视为组织损伤的起始条件。有关该模型的详细讨论，请参见先前的出版物⁴。

光化学损伤（约 10 秒 – 10⁴ 秒；约 10⁻³ – 10² W/cm²）：
在低辐照度下长时间暴露时，累积性光化学反应占主导地位，例如视色素漂白或活性氧（ROS）介导的通路。此类损伤发生在暴露时间超过 10 秒⁶（或超过 1 秒⁴）的情况下。因此，存在一个中间范围，在此范围内热损伤和光化学损伤可能同时发生¹⁶。慢性低水平蓝光危害即为典型示例。

方案

第1节描述了模型构建所需的一般步骤，这些步骤可借助多种软件完成。第2节以Altair Hypermesh（有限元建模与分析[FEA]软件）和Ansys Fluent（流体模拟软件）的特定示例情况，给出了具体操作说明。此外，还提供了一个补充文件（补充文件 1^{17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45}) 提供了建模方法的描述，并以非逐步的方式给出了理论背景（图1 - F图5，表1 - 表3).

1. 建模方法的实施——通用操作步骤

1. 根据上述几何结构创建三维模型。
   1. 在本研究中，根据补充文件1中描述的参数，在三维机械计算机辅助设计软件（版本16）中创建几何结构。
2. 根据眼球各部分的属性和延伸范围创建体网格，具体参数见补充文件1。
   1. 使用有限元分析软件（版本11）和计算机辅助工程（CAE）分析软件（版本11）进行网格划分和后处理。该模型包含761,766个四面体单元。网格单元的尺寸介于37 µm至491 µm之间。
      注意：若采用更粗糙的网格，将导致计算误差或偏差。增加单元数量虽可行，但会延长计算时间。本研究证实该网格尺寸已达到独立性，进一步细化单元尺寸对计算结果的影响可忽略不计。
   2. 为保持合理的计算时间，建议采用此网格，仅在特定区域需要更高分辨率时进行调整。相应地，将视网膜重新划分网格，最小特征尺寸设为5 µm，以实现后续协议中所述的视网膜色素上皮（RPE）层内的能量吸收。
      注意：图2展示了隐藏玻璃体后以显露视网膜的网格示例图像。
3. 集成边界条件。
   1. 在巩膜处，使用仅考虑传导的边界系数 h_scl = 20 W/m²K，角膜值 h'_corn,1 = 14.21 W/m²K，该值如上文详细推导所得。使用这些数值分别定义代表巩膜和角膜表面的边界条件。
4. 实施血流模拟。
   注意：该场量为每个网格单元提供流动方向和大小，并在能量方程中加入对流项。
   1. 按以下构建步骤进行：
      1. 设定生理上的流入/流出边界——前部血流来自虹膜大动脉环，后部血流来自短睫状动脉，静脉回流通过涡静脉排出。
      2. 计算所得的流场分布。
      3. 将所有局部流速按Peyman⁴⁵所述统一重缩放至5 mm/s，同时保持方向不变，通过流体模拟软件的用户自定义函数（UDF）实现。结果是在脉络膜区域内获得一个方向正确且幅值恒定的速度场。
5. 计算温度分布。
   1. 通过定义激光束随时间变化的位置和形状，开始温度计算。为获得视网膜上（或更广义地说，眼内）这些随时间变化的位置，需并行使用眼的光线追踪软件模型与热力学模型。
      注意：基于对所评估激光系统的光线追踪建模，该方法可提供激光光斑位置和形状的时间演化过程。出于可读性考虑，本文未详细说明此过程，但此前已有描述^1,18。
   2. 随时间将计算所得的光斑位置和形状映射到相应网格单元的区域，并为这些单元分配相应的激光功率。
   3. 通过UDF将这些数据导入流体模拟软件，软件以循环方式处理这些数据。UDF识别相关网格单元，并根据指定激光功率和对应照射时间引入能量源项。需同时考虑连续波（cw）和脉冲照射，以及静态或随时间变化的光束运动。
   4. 对每个位置，利用给定激光波长下的吸收系数计算被吸收的能量，并将其分配给受影响的网格单元，这些单元将被同步加热。
      注意：例如，在532 nm波长照射下，厚度为5 µm的RPE层中发生51%的吸收⁷，如图3所示。布鲁赫膜（Bruch's membrane）内不发生吸收^7,46。剩余激光能量在脉络膜中按朗伯-比尔定律衰减，假设脉络膜厚度为100 µm，吸收系数⁷为 α_Ch = 270 cm^-1。UDF计算每个网格单元厚度内的吸收能量，并将该值分配给该单元。从角膜到视网膜的透射率取20%⁷。若改变激光波长，仅需相应调整吸收系数和角膜至视网膜的透射率即可。
   5. 此外，需定义单个时间步长的持续时间，因为光斑形状数据（来自光线追踪模型，或若无需光线追踪模型时直接在UDF中定义）仅描述照射序列，而无时间尺度。
   6. 为加速计算，将所有无照射的时间段合并，并分配更大的时间步长，因为在这些时段无需精细的时间分辨率。
      注意：应用上述设置后，求解器可计算眼组织的瞬态热响应。由于当前模型仅需考虑热传导和血流，求解器仅限于求解相应的微分方程。示例性的用户自定义函数参见补充文件2。
6. 计算组织损伤。
   1. 为基于热响应推导损伤阈值，使用基于C++的工具从Fluent结果文件中提取计算所得的温度值（详见^1,18中的细节和代码列表）。
   2. 采用阿伦尼乌斯积分法，由这些温度值计算各个网格单元的阿伦尼乌斯值，并将其存储为与后处理软件兼容的格式。这使得在分析过程中可方便地在温度值与损伤概率之间切换。
   3. 最终，为每个网格单元获得二元损伤结果（是/否）。通过阿伦尼乌斯积分的阈值1⁴确定“是/否”分类。通过固定照射时间并测试不同输入功率，迭代确定该损伤阈值。
7. 设置求解器时间步长。
   1. 使用流体模拟软件（版本14.5）作为求解器，采用有限体积法。时间步长应满足库朗-弗里德里希斯-勒维（Courant-Friedrichs-Lewy）条件，该条件描述了时间步长与网格单元尺寸之间的比例关系，以确保模拟的稳定性。
   2. 逐步减小时间步长，直至结果不再随步长变化。例如，在总时长不超过10 s的计算中，使用1 ms和100 ms时间步长时，温度变化均被发现小于<1%。

图 2：穿过网格的截面图（玻璃体液被遮蔽）。该图显示了巩膜（白色）、脉络膜（红色）、视网膜（黄色）、晶状体（白色）、虹膜（绿色）以及房水（蓝色）¹⁸。 请点击此处查看此图的放大版本。

图3：眼底的吸收特性。在视网膜色素上皮（RPE）中，532 nm激光辐射的51%被吸收；Bruch膜中无吸收发生，脉络膜则根据Lambert-Beer定律行为建模¹⁸。 请点击此处查看该图的放大版本。

2. 示例说明——具体步骤

1. 根据上述几何结构创建一个三维模型。
   注意：三维模型是通过为补充表1中列出的所有元素创建对象而构建的。以下步骤以巩膜为例进行说明。
   1. 创建一个新的零件文件（.ipt）。
   2. 在XY平面上开始一个二维草图。
   3. 绘制一条水平构造线（这将作为旋转轴）。
   4. 为外表面绘制一个半径为12 mm的圆弧。
   5. 在中心（后极部），从外圆弧向内0.99 mm处创建一个点。
   6. 从该点开始绘制内表面曲线。
   7. 在角膜缘处，将内表面与外表面之间的距离设为0.75 mm。
   8. 在边缘处连接内外曲线，形成一个封闭轮廓。
   9. 选择3D模型 > 旋转。
   10. 选择曲面并确认。
2. 创建体网格。
   注意：网格基于三维对象创建。以下是网格划分的一个示例：
   1. 导入在三维机械计算机辅助设计软件中创建的三维模型。
   2. 点击需要划分网格的对象表面，例如巩膜。
   3. 首先对表面进行网格划分——点击该界面中的“Mesh”按钮（在选中表面后出现）。在界面中选择surfs，并选择自动自适应网格划分。
   4. 三维网格（用于定义对象）基于表面网格生成。要创建三维网格，请选择表面，然后从左上角菜单中点击：Mesh > Create > Solid Map Mesh。
   5. 此时将出现一个界面，点击mesh。
3. 集成边界条件。
   注意：必须为巩膜和角膜定义边界条件，因为这些表面与外部环境接触。以下以巩膜为例：
   1. 从顶部菜单进入Analysis > BCs > Create > Constraints (或 Loads)。
   2. 选择实体类型为表面，并从菜单中选择巩膜。
   3. 在同一面板中找到数值输入字段，直接在magnitude字段中输入数值。
   4. 点击create以分配边界条件。
4. 实施血流模拟。
   1. 定义前部流入（虹膜大动脉环）、后部流入（短睫状动脉）和流出（涡静脉）的表面。
   2. 为这些表面分配边界条件。
   3. 点击Analysis > BCs > Create > Loads。
   4. 选择inlet surfaces > assign velocity direction vectors（近似生理方向）。
   5. 选择outlet surfaces > assign pressure outlet（或根据求解器设置选择零约束）。
   6. 将模型导出至流体模拟软件——在此处将调整流速矢量的长度。
   7. 在流体模拟软件中使用附录中的ScaleVelocity UDF来调整矢量长度（UDF是一种标准格式，可加载到流体模拟软件中）。
5. 计算温度分布（设置求解器）。
   注意：使用流体模拟软件按以下方式计算温度，并将其保存为CAE分析软件可读取的文件格式。
   1. 点击File > Read > Mesh > Select mesh > Mesh > Check。
   2. 点击 General > Solver: Pressure-Based > Time: Transient。
   3. 点击 Models > Energy > Enable > OK。
   4. 点击 Materials > define ρ, p, kρ > Cell Zone Conditions > 为所有区域分配材料。
   5. 点击 Boundary Conditions > 定义thermal BCs（角膜、巩膜、环境）。
   6. 点击 Solution Initialization > Initialize。
   7. 点击 Run Calculation > 设置Time Step Size + Number of Time Steps + Iterations per step。
   8. 点击 Run Calculation > Calculate。
   9. 点击 File > Export > Solution Data… > File Type: Ensight Gold > 选择Temperature > 选择zones > Write。
6. 计算损伤。
   1. 在包含温度数据的有限元分析软件文件所在文件夹中，执行AddingArrhenius.exe（补充文件3和补充文件4）。该程序基于有限元分析软件文件中存储的温度数据，简单计算阿伦尼乌斯积分。

结果

本节首先展示所提出模型的验证与确认。随后，演示三个示例性应用。

验证与确认
本节首先将所提出的模型与已建立的热生理眼部模型^{19,27,28,30,32}进行比较，以验证其预测的温度分布（稳态情况）的一致性。在第二步中，通过将模拟得到的温度场和阿伦尼乌斯损伤预测结果与激光辐照下获得的实验数据（瞬态温度分布）进行对比，对模型进行验证。

由于长时间照射下的定量人体视网膜温度 体内 数据不可获得，因此血流效应的验证更具挑战性。对于短时间照射（数百毫秒量级），可利用视网膜手术期间获得的温度数据进行验证。这些测量结果与预测的温度升高值表现出良好的一致性。如前所...

讨论

在本方案中，模型构建过程没有关键步骤。关键在于网格尺寸和类型与时间步长的组合选择。就视网膜上采用其他扫描模式而言，建模技术的调整是相关的，但模型本身不应被修改。该技术仅限于相应章节中所述的损伤范围。其意义在于能够在无需动物实验的情况下预测损伤。其中一个最重要的未来应用是在涉及人眼激光照射（HEL）的场景中进行视网膜损伤的概率分析。

一旦该模型针对特定损伤机制得到验证（目前已被验证适用于热损伤），即可用于预测损伤阈值，其作用相当于实验所提供的ED50值。该ED50值随后可用于推导眼安全计算，而无需依赖现行的眼安全标准。为此，必须将ED50值与特定的安全系数相结合，这与IEC 60825标准中已有的假设一致。不同之处在于，由于ED50值现已针对所有参数组合被具体确定，安全系数可在所有情况下保持恒定，而无需反映ED50值本身存在的不确定性；它仅需反映将ED50值按安全系数进行缩放后仍可能导致损伤的概率。该安全系数的具体定义应由激光安全领域相关专家进一步讨论确定。

对于其他情况，由于场景具有不可预测性，标准方法不适用，此时将概率命中模...

材料

本文使用的材料清单

姓名	公司	目录编号	评论
Altair Hypermesh 11.0	Altair Engineering Inc.	版本 11	网格创建与显示；有限元建模与分析 [FEA] 软件
Ansys Fluent 14.5	Ansys Inc.	版本 14.5	热力学求解器；流体模拟软件
Autodesk Inventor	Autodesk	版本 16	三维机械计算机辅助设计软件
Hyperview	Altair Engineering Inc.	版本 11	计算机辅助工程 (CAE) 分析软件
Optic Studio 13 (Zemax)	Zemax Development Corporation，现为 Ansys Inc.	版本 13	光线追迹软件
Optocon FOTEMP2	Optocon	https://comem.com/en/optocon/	光谱仪；眼组织温度测量
Optocon TS2	Optocon	https://comem.com/en/optocon/	测量探针；眼组织温度测量