引　言

夏克–哈特曼波前传感器（Shack-Hartmann wavefront sensor，以下简称SH波前传感器）是一种用于测量光学系统中波前畸变的装置。其结构简单，精度高，实用性强，已被广泛应用于天文学、光学和医学等领域，用来测量光学系统中的像差^[1-3]。SH波前传感器是由微透镜阵列和相机两部分组成的光学波前探测器件。其中，微透镜阵列对随机波前进行空间采样，在相机的靶面上形成一个光斑阵列，通过计算光斑质心的相对偏移量，再利用合适的波前重构方法就能得到光学波前的相位分布。光斑质心探测精度对SH波前传感器的波前重建精度具有决定性的影响^[4-5]。波前检测系统通常存在各种噪声，包括相机的电子噪声、光学系统的杂光噪声，以及外部环境扰动引起的噪声等。眼波前像差检测系统是一种用于检测人眼光学像差的设备，已广泛用于眼视光学的临床诊断，其中基于SH波前传感器的眼波前像差检测技术被大多数厂家所采用。但是，传统的眼波前像差系统仅测量视轴方向的眼波前像差，不能反应视网膜周边的像差分布。近年来，由于视网膜周边离焦理论逐渐成为近视防控的主流学说，对视网膜周边的像差检测成为新的临床需求。采用光学扫描和SH波前传感器，人们可快速获得视网膜周边高分辨率的像差分布^[6-7]。但是对于扫描式眼波前像差检测系统，通过光学扫描法测量不同视场方向的波前像差，不同视场方向的杂光往往难以彻底消除。对于传统的波前像差测量系统，人们经常用孔径光阑来消除杂光的影响，但是对于扫描状态下的波前像差检测，透镜表面反射的杂光的大小与位置会随着扫描角度的改变而变化，很难通过在系统中放置光阑进行消除。此类噪声如果不滤除干净，对SH波前传感器的光斑质心计算会造成严重干扰，影响系统的稳定性、可靠性和准确性。在提高夏克–哈特曼波前传感器探测精度方面，夏爱利等^[8]提出了窗口优化的自适应阈值与线性插值结合的方法去除噪声；师亚萍等^[9]提出了Canny边缘检测算法与二阶质心算法等多种方法结合的自适应计算质心方法，有效地提高了质心定位的精度。大噪声会严重影响质心定位的精准度，甚至使光斑无法识别。但这些算法都无法对一些大的噪声进行滤除，从而限制了SH波前传感器的使用。因此，本文提出了一种大噪声情况下准确识别光斑和精准计算质心的新方法，即对窗口选取的噪声面积和信号光斑面积的大小进行比对，使得选取的噪声不纳入波前重建计算，以达到去除噪声影响的目的。由于信号光斑阵列数目较多，去除少数信号光斑（部分信号光可能与大噪声重叠）不会影响波前重建精度，这样就可在扫描式波前像差系统存在大噪声的情况下实现SH波前传感器光斑质心的准确计算，从而实现不同视场方向精准的波前重建，提高扫描式眼波前像差测量的稳定性和精确性。

1 夏克–哈特曼波前传感器的工作原理

夏克–哈特曼波前传感器的结构比较简单，主要由电荷耦合器件（charge coupled device，CCD）相机和微透镜阵列两部分组成。微透镜阵列由多个小透镜组成，当光束通过微透镜阵列时，每个微透镜都会将其分割的光束聚焦到探测器的焦点上，这些聚焦后的小光束会在探测器（如CCD相机）上形成一个焦点阵列^[10]。在理想情况下，若入射光波是平面波，微透镜阵列将形成一个规则的焦点阵列，常称为参考光斑阵列；如果光波前产生畸变，那么某些焦点的位置会偏离其参考光斑位置。通过测量每个焦点的实际位置与参考位置之间的偏差，可以计算出光波前的局部倾斜度。这些局部倾斜信息可以用来重建整个光波前的形状^[11-12]，表达式为

$ \frac{\text{∂} \overline{w}}{\text{∂} x}=\frac{\mathrm{\Delta }x}{f} \text{，} \frac{\text{∂} \overline{w}}{\text{∂} y}=\frac{\mathrm{\Delta }y}{f} $

(1)

式中：$ {\text{∂} \overline{w}}/{\text{∂} x} $，$ {\text{∂} \overline{w}}/{\text{∂} y} $分别表示子孔径区域波前在x方向和y方向的平均斜率；f表示透镜阵列的焦距。SH波前传感器工作原理如图1所示。

2 算法流程

本文算法的总流程如图2所示。

2.1 图像预处理

对于波前重建来说，图像的预处理尤为重要。预处理不仅提高了输入图像的质量，而且为后续质心计算和波前重建提供了更加清晰、准确的基础^[13]。波前重建的精度直接依赖于初始图像数据的质量和准确性。本文提出的算法，首先对全局使用中值滤波和先开运算后闭运算的数学形态学处理，然后对局部的光斑区域使用改进的最大类间方差法（简称OTSU），自适应阈值。它的计算方法是根据图像的灰度分布特性，将图像分成背景和目标两部分。分割的依据是两部分之间的间类方差最大，即类别内的差异最小化（使得每个子孔径都能获得一个最佳阈值）。

中值滤波用来减少图像中的噪声，特别是随机噪声，如散斑等。这一步骤在去除噪声的同时，有效地保持了图像边缘和关键结构的清晰度。开运算帮助去除图像中的小尺度噪点和不必要的细节，而闭运算则填补了焦点区域的小空洞和间隙，进一步增强了焦点的连续性和完整性。自适应阈值最大类间方差法通过自动计算，并应用最佳阈值，使得每个子孔径实现自适应阈值处理。

这一系列精细的预处理步骤能够显著提升输入图像的质量，确保后续波前重建算法的精准度和有效性。

2.2 大噪声情况下真实光斑的识别

1）子孔径自适应阈值的改进

传统的子孔径自适应阈值算法通常都是对子孔径内所有的灰度值进行排序，再选取第N位的灰度值作为阈值对图像进行二值化，或者是找出子孔径内最大亮度的光斑，然后根据亮度的百分比选取出灰度值作为阈值对图像进行二值化^[14]。以上算法中，所有的子孔径灰度值都需要减去当前的阈值，来达到噪声去除的目的。当二值化完毕后，选取出该子孔径体积最大的光斑作为真实光斑进行质心计算。但是这些算法都存在一定的不足，如果当背景噪声特别大的时候，背景噪声可能会被当作真实光斑进行选取，而真正的信号可能会被滤除。

根据上述问题，本文提出一种改进的算法。先是对每个子孔径使用OTSU阈值算法，从而得到每个子孔径的最佳阈值，然后再将每个子孔径的灰度值减去其对应的最佳阈值来达到去除每个子孔径噪声的目的。子孔径第i行、j列像素I_ij取阈值的计算式为

$ {I}_{ij}{{'}}=\left\{\begin{array}{c}{I}_{ij}-T,{I}_{ij}\geqslant T\\ 0\text{，}{I}_{ij} < T\end{array}\right. $

(2)

式中，T为OTSU阈值算法计算出的最佳阈值。但是当图像中存在大噪声时，该算法不再适用。

如图3所示，当单个子孔径内存在很大的噪声时，不管是传统的方法还是现在普遍使用的OTSU算法都无法正确地计算出最佳阈值。因为此时杂光相对比较强烈，OTSU计算出的最佳阈值可能是120，但是信号光强度明显低于这个值，当该子孔径错误地使用了该最佳阈值则会使得信号光被滤除，从而导致信号光的缺失，严重影响波前重建的准确性与稳定性。

本文的算法针对上述问题进行了改进。首先计算出所有子孔径的阈值，然后对所有子孔径的阈值进行筛选，筛选出异常阈值，并将异常阈值修正成正常的阈值。异常阈值的判定是对所有阈值进行排序，然后做差值计算，当其差值明显高于平均差值时，则判定其为异常阈值。虽然该操作并不能立刻消除大的噪声，但可以配合面积中值法来达到消除噪声的目的。

2）探测窗口的获取

探测窗口的大小和位置选取会直接影响质心定位的精度，过小的探测窗口可能会对光斑进行截断而遗漏有用的光斑信息，而过大的探测窗口可能会引入很多噪声导致质心定位的不准确^[15-16]，所以本文使用了Suzuki的轮廓跟踪算法来选取最合适的窗口来进行质心的定位（见图4）。

该算法先是运用了边缘检测提取出轮廓，然后使用Suzuki算法遍历图像中的所有边缘像素，将它们组织成一个链表。每个链表代表图像中的一个轮廓信息。然后对轮廓信息进行分析，提取出x和y方向最边界的点进行框选，得到最贴合轮廓的窗口。尽管预处理已经去除了大部分杂光，但是在某些情景下，预处理无法去除太大与顽固的杂光。这些杂光会被Suzuki轮廓跟踪算法框选。利用中值面积法进行进一步筛选，以达到去除杂光的目的。

3）中值面积法

在获得每个子孔径区域的光斑轮廓信息后，就可以计算每个轮廓的面积。

首先，把光斑轮廓面积按照大小进行排序，并储存在列表中。随后，取出排在中间的光斑轮廓面积作为参考。然后，将列表中所有的光斑轮廓面积A_s与参考光斑面积A_{s, med}进行对比，则有A_{s, med} n < A_s < A_{s, med} m（n、m为设置的比例系数）。符合这个范围的光斑将被筛选出来作为真实光斑。

4）质心的计算

使用Suzuki轮廓跟踪算法对光斑进行窗口选取后，再对每个窗口进行质心的计算

$ {X}_{\mathrm{h}}=\frac{{\displaystyle\sum} _{i,j}{x}_{i}{I}_{ij}}{{\displaystyle\sum} _{i,j}{I}_{ij}},{Y}_{\mathrm{v}}=\frac{{\displaystyle\sum} _{i,j}{y}_{i}{I}_{ij}}{{\displaystyle\sum} _{i,j}{I}_{ij}} $

(3)

式中：X_h，Y_v为当前窗口内的质心的横纵坐标；x_i，y_i是（i，j）像素点的横纵坐标；I_ij为（i，j）像素点的像素值。

3 实验与结果分析 3.1 实验平台搭建

本实验选用中心波长为850 nm的单模光纤输出超辐射激光二极管（SLD）（Inpenix Inc., USA）作为光源，选用焦距为150 mm的消色差透镜（Edmund Optics, USA）作为准直透镜。把输出光纤端面放在透镜的焦点处，使其经过透镜的光束平行射出，装置如图5所示。SH波前传感器的参数为：微透镜阵列的焦距为14.2 mm；CCD的像元尺寸为3.45 μm；图像的视场像素为1 600 × 1 200，在这个视场下的子孔径数量为18 × 13。

光束最后进入夏克–哈特曼传感器，获取点阵图像。撤去透镜，测得光纤端面到传感器的距离为538.4 mm，将其作为初始球面波半径。再通过光学滑轨，水平移动光纤端面，每次移动间隔为50 mm（读取精度为1 mm），共取了10个球面波半径，并对10个球面波半径生成的10个点阵图像加入相同的噪声，如图6所示。

Zernike多项式在直角坐标系下表示波前像差函数的表达式为

$ W\left(x,y\right)=\sum _{i=0}^{\infty }{c}_{i}{Z}_{i}(x,y) $

(4)

上述表达式在实际应用中通常取其有限项，当选取项数k后，则i→∞变为i→k。Z_i（x, y）为Zernike多项式第i项的表达式，当确定了项数k后就有一组系数c（波前重建计算得到）对应k项多项式的线性组合。

使用峰谷值（V_pv）与均方根值（V_rms）进行算法的稳定性与准确性评估，计算式为

$ V_{{\mathrm{rms}}}\left(c\right)=\sqrt{\sum _{i=0}^{k}{\left({c}_{i}\right)}^{2}} $

(5)

式中：k为Zernike取的项数；c_i为Zernike系数的第i项。

$ V_{{\mathrm{pv}}}={W(x,y)}_{\mathrm{m}\mathrm{a}\mathrm{x}}-{W(x,y)}_{\mathrm{m}\mathrm{i}\mathrm{n}} $

(6)

式中：W（x,y）_max为波前像差函数最大值；W（x,y）_min为波前像差函数最小值。为了统一分析，将Zernike系数的常数项以及x与y的倾斜项置0，再进行峰谷值的计算。

3.2 结果分析

为了验证本文算法的性能，首先在平面波获取的点阵图像（见图7（a））中添加椒盐噪声，然后加入大的背景噪声，如图7（b）～（e）所示。由于原图像分辨率较大，所示图像被压缩，椒盐噪声无法清晰显示。图7（f）为在图7（c）中随机截取部分光斑，放大显示出的椒盐噪声效果。然后对图7（b）～（e）分别应用传统算法与本文的算法进行处理，具有代表性效果的截取部分见图8（a）和（b）。结果显示，传统的算法错误地把大噪声进行了窗口选取，并且对大噪声附近的有用信号光进行了错误的滤除。而在使用了Suzuki轮廓跟踪算法、本文的中值面积算法和改进的子孔径OTSU算法后，在大噪声的情景下仍能正确地避开大噪声，并对真实有效的光斑进行正常的窗口选取。这种方法不仅能选出强度较弱的光斑，而且还能使得在大噪声附近的子孔径阈值不受大噪声的影响，从而精准地选取到每个真实光斑进行质心的计算。

最后，选用4 mm的瞳孔直径对原始点阵图像以及添加了噪声的点阵图像进行波前重建，通过数据分析评估本文算法的性能。结果如表1所示。

通过分析得出，本文算法的稳定性和准确度都高于传统算法，并且误差基本保持在0.015λ之内，V_rms的误差也基本保持在0.000 1λ之内。

为了进一步验证本文算法的可靠性，同样选取4 mm瞳孔直径，再次使用本文算法对加入了大噪声的球面波进行验证，并采用理论峰谷值（V_{pv, th}）进行对照。

$ V_{\mathrm{pv, th}}=R-\sqrt{{R}^{2}-{r}^{2}} $

(7)

式中：R为球面波半径；r为选取的瞳孔半径。

对于球面波而言，计算实际V_pv和平面波一样需要将前3项置0。并且为了更准确地评估噪声对波前重建带来的影响，在计算球面波V_rms时，需要将常数项、x与y的倾斜项，还有离焦项置0。

从表2与表3可以发现，本文算法在球面波的大噪声滤除方面依旧保持着良好的稳定性与准确性，原图的V_pv与V_{pv, th}基本一致（存在测量误差），并且原图经过了大噪声与椒盐噪声之后，计算的V_pv误差在0.015λ之内，V_rms误差也基本保持在0.001λ之内，即本文算法的稳定性和准确度都再次得到了验证。

4 结　论

本文介绍了一种组合算法，其在较大噪声情况下依然能够进行可靠、稳定、准确的质心计算和波前重建，为扫描式眼波前像差检测系统无法同时消除不同视场方向系统杂光的问题提供了一种解决方案。

目前很多算法的基本流程都是先进行窗口的选取，然后使用相应算法对每个窗口内的光斑区域进行判别，去掉窗口内的噪声，找出真实的光斑，再进行质心的计算。而本文算法不需要对每个光斑区域的杂光进行一一判别，只需要对每个光斑区域进行一次改进的OTSU自适应阈值，然后使用Suzuki轮廓跟踪算法，选出最贴合的窗口，最后通过本文的中值面积法进行统一的处理，极大地提升了波前重建的效率，并且能够在大的噪声条件下保证其准确性与稳定性。

尽管本文的算法能够做到大小杂光的滤除，但是仍然存在一些局限性。本文算法实现的依据是，经过预处理后的图像在经过排序后的中间值的光斑面积一定是真实光斑面积。所以在进行图像处理，选取重建区域的时候，一定要让区域内的真实光斑数量尽可能的多于杂光数量，这样才能保证算法的准确性。