Analyzing the Effect of the Intra-Pixel Position of Small PSFs for Optimizing the PL of Optical Subpixel Localization

Abstract

Subpixel localization techniques for estimating the positions of point-like images captured by pixelated image sensors have been widely used in diverse optical measurement fields. With unavoidable imaging noise, there is a precision limit when estimating the target positions on image sensors, which depends on the detected photon count, noise, point spread function (PSF) radius, and PSF’s intra-pixel position. Previous studies have clearly reported the effects of the first three parameters on the precision limit but have neglected the intra-pixel position information. Here, we develop a localization precision limit analysis framework for revealing the effect of the intra-pixel position of small PSFs. To accurately estimate the precision limit in practical applications, we provide effective PSF (ePSF) modeling approaches and apply the Cramér-Rao lower bound. Based on the characteristics of small PSFs, we first derive simplified equations for finding the best precision limit and the best intra-pixel region for an arbitrary small PSF; we then verify these equations on real PSFs. Next, we use the typical Gaussian PSF to perform a further analysis and find that the final optimum of the precision limit is achieved at the pixel boundaries when the Gaussian radius is as small as possible, indicating that the optimum is ultimately limited by light diffraction. Finally, we apply the maximum likelihood method. Its combination with ePSF modeling allows us to successfully reach the precision limit in experiments, making the above theoretical analysis effective. This work provides a new perspective on combining image sensor position control with PSF engineering to make full use of information theory, thereby paving the way for thoroughly understanding and achieving the final optimum of the precision limit in optical localization.

Content

《1、引言》

1、引言

在光学测量中，类点目标是最基础和典型的光学目标。如图1（a）所示，点光源发出的光子通过光学系统，以斑点的分布形式（称为点扩散函数，PSF）出现在像素化图像探测器上，在像素采样后形成数字图像。根据像素响应值，PSF在探测器上的位置（x*, y*）可以以亚像素精度被估计。此过程称为光学亚像素定位，在不同光学测量领域有广泛应用[1‒9]。在天文观测中，望远镜[10]或星敏感器[11]对恒星进行成像和定位，从而进行黑洞探测[12‒13]、航天器导航[14‒16]等宇宙探测任务。在生物显微中，结合间歇激活的荧光探针，单分子的亚像素定位产生了超越衍射极限（DL）的纳米级重建图像，使生物学家能够以前所未有的分辨率观察细胞结构[17‒20]。在这些光学测量应用中，最基本的问题涉及亚像素定位的精度和精度极限。

《图1》

图1 小PSF像素相位对定位精度极限的影响。（a）点源的光学成像和亚像素定位过程，即根据像素响应估计PSF在像素化图像探测器上的位置(x*, y*)。（b）忽略PSF的像素相位，估计x*时的平均精度极限与高斯半径之间的关系。三个不同PSF半径及其精度极限结果标记在曲线上，三个PSF及其成像显示在右侧。精度极限计算细节如下：假设PSF为二维高斯函数，探测到的光子数为1500，添加与实际应用相仿的典型噪声（见第3.1节）。在这种情况下，半径为0.27像素的PSF达到了0.012像素的精度极限。（c）半径为0.50像素（顶部）、0.27像素（中心）和0.16像素（底部）的三个PSF在不同像素相位上的精度极限。在右侧比较了y = 0.50像素的横截面上三个PSF的精度极限曲线，其中，最小的PSF在像素边界处达到0.005像素的精度极限。PSF：点扩散函数；r：PSF的半径。

由于不可避免的成像噪声，如光子散粒噪声和像素暗噪声等，亚像素定位存在理论精度极限（PL），其主要与探测光子数N、像素暗噪声σ_d²、PSF半径r，以及PSF像素相位（x, y）（表示相对于像素边界的PSF亚像素位置）有关。统计学中的Cramér-Rao下界（CRLB）描述了估计参数时无偏估计量的最小方差，可以用于揭示和优化定位精度极限[21‒23]。现有研究已清楚阐明前三个变量对精度极限函数的影响[23‒25]，但由于假设PSF位置在像素内随机分布，并对所有像素相位的精度极限取均方根值，忽略了PSF像素相位对精度极限的影响。在此种情况下，当探测光子数和噪声水平固定时，中等大小的PSF使精度极限PL (N, σ_d², r)达到最优值[图1（b）]。这些研究适用于尺寸相当于或大于像素尺寸的望远镜或显微镜典型PSF [21,26]，对于此类PSF，在不同的像素相位，精度极限几乎没有变化。然而，对于尺寸很小的PSF（如r＜0.3像素），像素相位效应不能再被忽视，即精度极限随着PSF像素相位不同而发生显著变化，且在像素边界附近可以达到非常高的精度[图1（c）]。涵盖小PSF像素相位效应的定位精度极限PL (N, σ_d², r, x, y)的优化问题尚未得到足够的理论和实验研究。此外，现有精度极限估计通常假设PSF是理想解析函数，如高斯函数[23‒25,27]等，但光学像差、像素采样、像素不均匀响应等非理想成像条件会导致实际PSF发生偏差，降低了实际定位应用中精度极限预测的准确性。

仅当实际应用中的定位精度能达到理论精度极限时，精度极限分析才有意义。最大似然估计（MLE）已被应用于实现一些光学测量系统的精度极限[28‒31]。然而，其并未成为小PSF系统如星敏感器定位研究中的共识。与望远镜相比，星敏感器的F数小，曝光时间短，这使星敏成为应用亚像素定位技术的典型小PSF光学系统。传统的星敏感器恒星定位方法主要分为两类：质心法和拟合法。质心法（CG）计算图像像素矩阵的强度矩[32]，是快速但有偏的估计器，会导致S形周期性系统误差（通常为0.03~0.10像素）[33]，已有研究主要围绕补偿系统误差和调整提取阈值进行[33‒37]。拟合法通过最小二乘拟合用PSF模型拟合所测像素数据，如高斯拟合法（GF）[32]，研究主要围绕降低计算量进行，如高斯网格法和高斯分析法的提出[38‒40]。这些GF类方法在理想情况下能实现高精度，但在实际非理想应用中无法实现精度极限。主要原因是真实PSF模型（尤其是真实小PSF）偏离高斯函数，而拟合不准确的模型会引入系统误差（附录A中的S1节）。其次，与MLE不同，最小二乘拟合没有充分利用噪声信息（即噪声概率密度函数）[27]。总结来说，在该领域，实现信息论提供的定位精度极限的技术还未被建立。

在此，本文提出一种定位精度极限分析框架，用于揭示小PSF的像素相位效应。该框架首先提供了有效PSF（ePSF）建模方法来准确重建实际应用中的真实PSF，通过将CRLB应用于ePSF模型，获得定位精度极限的准确估计。其次，基于小PSF的特性，推导了求解任意小PSF最优精度极限和最佳像素相位的简化公式，并以典型高斯PSF为例，进一步推导了精度极限PL (N, σ_d², r, x, y)在PSF半径最小和像素边缘处取得的最优值，从而揭示了最优定位精度极限最终受限于光的衍射极限。最后，在实际应用中利用MLE拟合ePSF，实现目标定位的理论精度极限。这项工作为光学亚像素定位提供了深入的物理见解，将高精度定位研究扩展到小PSF领域，并提供了将图像探测器位置控制与PSF工程相结合以优化定位性能的新视角，从而促进广泛光学测量领域的亚像素定位向最终精度极限发展。

《2、方法》

2、方法

《2.1 ePSF建模和噪声分析》

2.1 ePSF建模和噪声分析

该框架首先提供了系统ePSF建模方法和随机噪声分析，以准确描述实际定位应用中像素响应的概率密度函数，这是后续精度极限估计和实现的重要需求。为避免解析PSF和实际PSF之间常存在的系统偏差，对实际应用中测量的像素强度和目标位置之间的精确关系进行建模，即形成实验ePSF模型[41‒42]：

$P_{i j} = φ_{i j} (Δ x, Δ y) + s_{i j}$ (1)

式中，_ij是像素(i, j)的响应值；_ij是此像素的背景值；_ij(Δx, Δy)是ePSF模型，其中，Δx = x* - (j + 0.5)，Δy = y* - (i + 0.5)，因此，(Δx, Δy)表示目标相对于像素中心的位置。由于光学像差和像素响应的不均匀性，不同像素对应的ePSF略有不同，故_ij(Δx, Δy)用像素下标ij标记。

在此根据应用需求和条件提供两种ePSF建模方法。第一种方法为形成非常精确的ePSF，每个_ij(Δx, Δy)对应于一个目标和一个特定像素，适用于具有充足标定条件的情况。在建模时，目标和图像探测器之间需要精确的亚像素级相对位移。为形成相对位移，可以利用高精度转台旋转光学系统，或者可以使用纳米位移台仅仅移动图像探测器。图2（a）显示了1D定位情况下的ePSF建模过程，1D定位指2D PSF在2D图像探测器上移动，但仅沿着1D移动，并且需要估计该维度中的位置。如图2（a）中的左图所示，假设目标以亚像素步长Δh在x方向上移动，其y*是固定的。在每个位置采集多帧图像，对其进行平均以减少噪声，然后可以获得像素(i, j)的平均像素值和目标位置的离散函数，其间隔为Δh [图2（a）中图]。ePSF随着目标的缓慢移动而平滑变化，因此使用具有非结点边界条件的三次样条插值来重建连续函数。最后，根据式（1）的简单变换，可以获得ePSF模型_ij(Δx, Δy) [图2（a）中的右图]。需要注意的是，x方向上相邻像素对应的ePSF [如φ₂₃(Δx)、φ₂₄(Δx)、φ₂₅(Δx)]略有不同，但在图中并未明显反映出来。标定间隔Δh与PSF大小有关，一般可以设置为0.05~0.25像素（对Δh影响的仿真分析如附录A中的S2节所示）。对于2D定位情况，目标在两个方向上移动，ePSF建模过程不变。这种方法的优点是所建ePSF是涵盖光学与像素不均匀性信息的精确实验模型，但是需要大量的标定工作。在本文中，这种ePSF建模方法被应用于实验室实验，以准确估计和实现定位精度极限，并验证所推导的小PSF像素相位效应公式。

《图2》

图2 ePSF建模和噪声分析。（a）第一种建模方法获得精确ePSF的过程。目标位于间隔为Δh的不同位置的图像如左图所示，10个像素（第2~3行，第2~6列）用红色标记，它们的像素值与目标位置的关系如中图所示。通过式（1）进行简单变换，计算了6个像素（第2~3行，第3~5列）的ePSF，如右图所示。（b）第二种建模方法获得近似ePSF的过程。中图表示了从三张图像中获得的ePSF模型φ(Δx, Δy)的样本，评估图中黑色的网格点(0, 0)时，对红色窗口内的样本进行平均，得到φ(0, 0)。（c）基于标准像素响应模型的噪声分析。本文所作假设和近似用红色文字标记。DN：数字值；FES：波动电子信号；η：量子效率；l_ij：像素(i, j)接收的光信号；d_ij：暗信号；d_0,_ij，d_1,_ij：暗信号的常数部分和泊松部分；σ_q²：量化噪声；K：像素增益。

第二种方法是形成对应于一个目标的近似ePSF，即φ(Δx, Δy)，而忽略光学和像素响应的不均匀性。实际应用中很多情况不能满足目标以亚像素步长在图像探测器上运动的条件，目标在图像中的位移可能是几个像素，并且通常是未知的[图2（b）中的上图]。经实验测试，在此采用了哈勃望远镜星图处理方法[41]的简化方法。对于每帧图像，提取感兴趣区域（ROI）内目标的像素值，通过传统CG方法估计目标的初始位置，忽略像素不均匀性，通过公式（1）的简单变换，将像素值变成φ(Δx, Δy)的采样样本。例如，在图2（b）中，使用CG方法估计图像1中目标的位置为（1.4, 1.4）。通过P₁₁ ‒ s₁₁ = φ (1.4-1.5, 1.4-1.5)，可以将像素（1, 1）的值映射到ePSF的φ (-0.1, -0.1)位置，通过P₁₂ ‒ s₁₂ = φ (1.4-2.5, 1.4-1.5)可以将像素值P₁₂映射到φ (-1.1, -0.1)。如果目标的ROI区域包含3 × 3像素，则每帧图像可以获得φ (Δx, Δy)的9个采样点。

在获得了来自多帧图像的所有样本后，在q × q网格点以Δh的间隔评估ePSF。如图2（b）中图所示，Δh设置为0.25像素，绘制间隔为Δh的蓝色虚线，蓝色虚线的交点表示评估ePSF的网格点。对于每个网格点，对Δx和Δy方向Δh距离内的样本进行平均，并去除偏离平均值2.5σ的样本（其中，σ是这些样本的标准差），最后的平均值用作该网格点处的ePSF值。最后，进行三次样条插值以获得最终的ePSF [图2（b）下图]。本文将这种方法应用于真实恒星观测。为方便，这里的ePSF模型包括目标能量，使得每个ePSF对应于一个目标，可以对ePSF进行能量归一化，使其适应具有相同强度分布但不同能量的多个目标。

接下来，基于欧洲机器视觉协会（EMVA）1288标准[43]中的像素响应模型对成像噪声进行分析。如图2（c）所示，光子到达像素产生_ij数量的电子（_ij平均值为 $μ_{l_{i j}}$ ，方差为 $σ_{l_{i j}}^{2}$ ）。光生电子的数量具有随机波动的特性，称为“散粒噪声”，并遵循泊松分布，因此 $μ_{l_{i j}} = σ_{l_{i j}}^{2}$ 。除了光生信号之外，像素值还包含暗信号_ij（其平均值为 $μ_{d_{i j}}$ ，方差为 $σ_{d_{i j}}^{2}$ ）。两个信号的和通过放大增益K，并经过模数转换器（ADC）进行数字转换，成为像素值_ij。为方便分析，本文做了一些假设和近似：ADC量化噪声（方差为σ_q²）与光噪声和暗噪声相比非常小，在这里被忽略；由于暗信号的噪声主要由热诱导的泊松分布电子引起，假设暗信号包括两个部分：常数部分d_0,_ij（平均值为 $μ_{d_{0, i j}}$ ，方差为0）和泊松部分d_1,_ij。在像素值侧，暗信号产生的背景_ij也有两个相应部分，即常数部分诱导的s_0,_ij和泊松部分诱导的s_1,_ij。在该模型中，由_ij + d_1,_ij或(_ij - s_0,_ij)/K表示的光信号_ij和暗泊松信号d_1,_ij的和被称为波动电子信号（FES），FES遵循泊松分布，其平均值（μ_FES）和方差（σ $_{F E S}^{2}$ ）与ePSF模型_ij(Δx, Δy)之间的关系为：

$μ_{F E S} = σ_{F E S}^{2} = μ_{{l_{i}}_{j}} + μ_{{d_{1}}_{, i j}} = φ_{i j} (Δ x, Δ y) / K + σ_{d_{i j}}^{2}$ (2)

因此，对于像素值_ij，有以下概率密度函数：

$p [P_{i j} | (x^{*}, y^{*})] = p [F E S = (P_{i j} - s_{0, i j}) / K | (x^{*}, y^{*})] = \frac{[φ_{i j} (Δ x, Δ y) / K + σ_{d_{i j}}^{2}]^{(P_{i j} - s_{0, i j}) / K}}{[(P_{i j} - s_{0, i j}) / K]!} e^{- [φ_{i j} (Δ x, Δ y) / K + σ_{d_{i j}}^{2}]}$ (3)

《2.2 精度极限估计》

2.2 精度极限估计

基于描述测量像素值的准确概率密度函数，可以根据CRLB理论准确估计实际定位应用中的精度极限[23]。基于式（3），任意ROI中像素值矩阵的联合概率密度函数为：

$\begin{array}{l} p [P | (x^{*}, y^{*})] = \\ \prod_{i j \in R O I} \frac{[φ_{i j} (Δ x, Δ y) / K + σ_{d_{i j}}^{2}]^{(P_{i j} - s_{0, i j}) / K}}{[(P_{i j} - s_{0, i j}) / K]!} e^{- [φ_{i j} (Δ x, Δ y) / K + σ_{d_{i j}}^{2}]} \end{array}$ (4)

无偏估计方法估计x*的CRLB为：

$C R L B_{x} = \frac{1}{\sum_{i j, i j \in R O I} \frac{[\partial φ_{i j} (Δ x, Δ y) / \partial x^{*}]^{2}}{K φ_{i j} (Δ x, Δ y) + K^{2} σ_{d_{i j}}^{2}} - {\sum_{i j, i j \in R O I} {\frac{[\partial φ_{i j} (Δ x, Δ y) / \partial x^{*}] [\partial φ_{i j} (Δ x, Δ y) / \partial y^{*}]}{K φ_{i j} (Δ x, Δ y) + K^{2} σ_{d_{i j}}^{2}}}}^{2} / \sum_{i j, i j \in R O I} \frac{[\partial φ_{i j} (Δ x, Δ y) / \partial y^{*}]^{2}}{K φ_{i j} (Δ x, Δ y) + K^{2} σ_{d_{i j}}^{2}}}$ (5)

计算细节如附录A中的第S3节所示。估算x*的精度极限为 $P L_{x} = \sqrt[]{C R L B_{x}}$ ，估计y*类似[只需在式（5）中交换x*和y*的位置]。在实际应用中，可以使用光子转移法获得图像传感器的K和 $σ_{d_{i j}}^{2}$ [43]。值得注意的是，_ij(Δx, Δy)及其导数随目标位置的变化而变化，因此，目标的像素相位对定位精度极限会产生影响。

《2.3 小PSF像素相位效应分析》

2.3 小PSF像素相位效应分析

由于公式（5）的复杂形式，很难在解析层面研究像素相位对精度极限的影响。然而，对于非常小的PSF，可以基于小PSF特性对公式进行简化。

在估计小PSF的x*时，y*的影响很小，因此假设PSF能量仅分布在像素阵列的某一行中。当PSF位于像素(i, j)内除像素边缘以外的大多数位置时，PSF的能量I（单位：数字值DN）几乎集中在该像素中，有以下近似：

$\{\begin{matrix} R O I \approx 1 p i x e l \\ φ_{i j} \approx I \\ \frac{\partial φ_{i j}}{\partial x^{*}} \approx 0 \end{matrix}$ (6)

为表达方便，在此将_ij(Δx, Δy)写为_ij。公式（5）的倒数，表示包含在测量像素数据中的目标位置的信息量，可以化简为：

$\frac{1}{C R L B_{x}} \approx \frac{[\partial φ_{i j} / \partial x^{*}]^{2}}{K φ_{i j} + K^{2} σ_{d_{i j}}^{2}} - \frac{[\partial φ_{i j} / \partial x^{*}]^{2} [\partial φ_{i j} / \partial y^{*}]^{2}}{(K φ_{i j} + K^{2} σ_{d_{i j}}^{2})^{2}} / \frac{[\partial φ_{i j} / \partial y^{*}]^{2}}{K φ_{i j} + K^{2} σ_{d_{i j}}^{2}} = 0$ (7)

因此，不可能精确地定位能量只分布在一个像素中的小PSF，实际小PSF系统常通过轻微离焦来避免这种情况的发生。

然而，当小PSF位于像素边缘附近时，结果完全不同。此时PSF能量分布在两个相邻的像素中，假设较大的像素值小于较小像素值的十倍，否则可以近似为上述的能量集中在一个像素的情况。在此，有以下近似：

$\{\begin{matrix} R O I \approx 2 p i x e l s \\ φ_{i j} + φ_{i, j + 1} \approx I \\ \partial φ_{i j} / \partial x^{*} \approx - \partial φ_{i, j + 1} / \partial x^{*} \\ [\partial φ_{i, j + 1} / \partial y^{*}] / [\partial φ_{i, j} / \partial y^{*}] \approx φ_{i, j + 1} / φ_{i j} = α, α \in [1 / 10,10] \end{matrix}$ (8)

式中，α是两个像素值的比值。另外，假设对于不同的像素，暗噪声是均匀的，其方差由σ_d²表示，则公式（5）的倒数化简为：

$\frac{1}{C R L B_{x}} \approx [\partial φ_{i j} / \partial x^{*}]^{2} {\frac{I / K + 2 σ_{d}^{2}}{(φ_{i j} + K σ_{d}^{2}) (φ_{i, j + 1} + K σ_{d}^{2})} - \frac{K [(1 - α) σ_{d}^{2}]^{2}}{[α^{2} (φ_{i j} + K σ_{d}^{2}) + φ_{i, j + 1} + K σ_{d}^{2}] (φ_{i j} + K σ_{d}^{2}) (φ_{i, j + 1} + K σ_{d}^{2})}}$ (9)

如果光信号比暗噪声大得多，则公式（9）可以简化如下：

$\frac{1}{C R L B_{x}} \approx [\partial φ_{i j} / \partial x^{*}]^{2} \frac{I / K}{φ_{i j} φ_{i, j + 1}}$ (10)

那么，估计x*的精度极限为：

$P L_{x} = \sqrt[]{C R L B_{x}} \approx \{\begin{matrix} \frac{1}{| \partial φ_{i j} / \partial x^{*} | \sqrt[]{\frac{I / K + 2 σ_{d}^{2}}{(φ_{i j} + K σ_{d}^{2}) (φ_{i, j + 1} + K σ_{d}^{2})} - \frac{K [(1 - α) σ_{d}^{2}]^{2}}{[α^{2} (φ_{i j} + K σ_{d}^{2}) + φ_{i, j + 1} + K σ_{d}^{2}] (φ_{i j} + K σ_{d}^{2}) (φ_{i, j + 1} + K σ_{d}^{2})}}}, \begin{matrix} w i t h & σ_{d}^{2} \end{matrix} \\ \frac{1}{| \partial φ_{i j} / \partial x^{*} |} \sqrt[]{\frac{φ_{i j} φ_{i, j + 1}}{N}}, \begin{matrix} w i t h o u t & σ_{d}^{2} \end{matrix} \end{matrix}$ (11)

式中，N表示探测到的光电子数量，N = I/K； $|\frac{\partial φ_{i j}}{\partial x^{*}}|$ 为 $\frac{\partial φ_{i j}}{\partial x^{*}}$ 的绝对值。值得注意的是，小PSF在像素边缘附近的 $|\frac{\partial φ_{i j}}{\partial x^{*}}|$ 会非常大，因此定位精度很高。式（11）可用于计算任意小PSF的最优精度极限和对应的最佳像素相位。此外，由于ePSF为准确的实验模型，式（11）可以在实际小PSF定位应用中提供准确的精度极限估计。

接下来，以典型高斯函数作为PSF，对精度极限的最优值作进一步分析。因上文假设PSF能量分布在像素阵列的一行中，则_ij可以用1D高斯函数的积分表示，即：

$\{\begin{matrix} φ_{i j} = \frac{I}{r \sqrt[]{2 π}} \int_{j}^{j + 1} e^{- \frac{(x - x^{*})^{2}}{2 r^{2}}} d x \\ [\partial φ_{i j} / \partial x^{*}]_{x^{*} \approx j + 1} = \frac{I}{r \sqrt[]{2 π}} [e^{- \frac{(j - x^{*})^{2}}{2 r^{2}}} - e^{- \frac{(j + 1 - x^{*})^{2}}{2 r^{2}}}] \approx - \frac{I}{r \sqrt[]{2 π}} e^{- \frac{(j + 1 - x^{*})^{2}}{2 r^{2}}} \end{matrix}$ (12)

式中，x* ≈ j + 1表示PSF位于像素边界j + 1附近。结合式（11）和式（12），可以计算高斯函数在像素边界附近的精度极限。可以发现，当r → 0和x* = j + 1时，式（9）达到其最大值，也就是说，此时PSF正好位于像素边界，即_ij = _i_,_j₊₁ = I/2，式（9）最大值为：

${[\frac{1}{C R L B_{x}}]}_{m a x} \approx \frac{2 {(I / K)}^{2}}{π r^{2} (\frac{I}{K} + 2 σ_{d}^{2})} = {\frac{2 N^{2}}{π r^{2} (N + 2 σ_{d}^{2})}}_{m a x}$ (13)

因此，精度极限的最优值为：

$[P L_{x}]_{o p t i m u m} \approx \{\begin{matrix} \frac{r}{N} \sqrt[]{\frac{π (N + 2 σ_{d}^{2})}{2}}, \begin{matrix} w i t h & σ_{d}^{2} \end{matrix} \\ r \sqrt[]{\frac{π}{2 N}}, \begin{matrix} w i t h o u t & σ_{d}^{2} \end{matrix} \end{matrix}$ (14)

到此，已经解决了高斯PSF的PL_x(N, σ_d², r, x)优化问题。在探测到一定数量光子的情况下，当PSF半径尽可能小时，精度极限在像素边界处达到最优，式（14）可用于估计该最优值。

但在物理上，由于光的衍射效应，PSF半径r不可能无限小，衍射极限DL（单位：像素）为：

$[D L] = \frac{1.22 λ F}{a}$ (15)

式中，λ表示波长；F是光学系统的F数；a为像素尺寸。DL描述了从艾里斑模型的峰值到第一个零点的距离，而艾里斑常近似为高斯函数，对于该函数，距离峰值三倍高斯半径内的区域几乎覆盖了所有能量。因此，艾里斑的DL和其近似高斯函数的半径之间的关系为：

$[D L] \approx 3 r$ (16)

结合式（14）和式（16）可以发现，定位精度极限的最优值最终受到衍射极限的限制，可由下式表示：

$[P L_{x}]_{o p t i m u m} \approx \{\begin{matrix} \frac{[D L]}{3 N} \sqrt[]{\frac{π (N + 2 σ_{d}^{2})}{2}}, \begin{matrix} w i t h & σ_{d}^{2} \end{matrix} \\ \frac{[D L]}{3} \sqrt[]{\frac{π}{2 N}}, \begin{matrix} w i t h o u t & σ_{d}^{2} \end{matrix} \end{matrix}$ (17)

《2.4 基于MLE的定位精度极限实现方法》

2.4 基于MLE的定位精度极限实现方法

最后，将MLE应用于所建立的ePSF模型，以在每个像素相位达到定位精度极限。使测量像素值概率 $p [P | (x^{_{*}}, y^{_{*}})]$ 最大化的位置即是定位结果，测量像素矩阵概率的负对数用作代价函数：

$χ = - l n p [P | (x^{*}, y^{*})] = \sum_{i j \in R O I} [φ_{i j} (Δ x, Δ y) / K + σ_{d_{i j}}^{2}] + \sum_{i j \in R O I} l n [((P_{i j} - s_{0, i j}) / K)!] - \sum_{i j \in R O I} [(P_{i j} - s_{0, i j}) / K) l n (φ_{i j} (Δ x, Δ y) / K + σ_{d_{i j}}^{2})]$ (18)

迭代牛顿-拉斐逊方法用于最小化代价函数，求解目标位置：

$[\begin{matrix} δ x^{*}_{}^{[n]} \\ δ y^{*}_{}^{[n]} \end{matrix}] = - (H^{[n]})^{- 1} J^{[n]}$ (19)

式中，和分别为代价函数χ的雅可比矩阵和黑森矩阵（计算细节如附录A中的第S3节所示），初始迭代值 $x^{* [0]}$ 、 $y^{* [0]}$ 可以使用CG方法来计算，ePSF模型_ij(Δx, Δy)及其导数可以提前标定并制成表格。通常，定位过程很快，在很少的迭代次数中即可达到收敛（在本文实验中通常为2~4次迭代）。

《3、结果》

3、结果

《3.1 数值仿真》

3.1 数值仿真

首先使用仿真数据来表征精度极限分析框架的性能。对具有不同半径的高斯函数进行像素积分，在不同的像素相位生成模拟目标，引入泊松分布的光噪声和暗噪声，通量I设置为300DN，像素系统增益K设置为0.2，因此探测光电子数N为1500，噪声1的暗噪声方差 $σ_{d_{i j}}^{2}$ 设置为25（根据噪声模型，对应于像素值的1DN标准差），噪声2的暗噪声方差设置为250，所设参数（除了噪声2）与实验中测得的参数接近。

使用第2.1节中的第一种方法进行ePSF建模，利用式（5）在每个像素相位计算x*估计的精度极限。对于每个PSF，计算所有像素相位的平均精度极限（通过取RMS），并计算最优精度极限和最佳像素相位。图3（a）红色曲线显示了先前研究所述的平均精度极限和PSF半径的关系，蓝色曲线表示了最优精度极限与PSF半径的函数。结果表明，尽管非常小的PSF具有较差的平均精度极限，但它可以在像素边缘处实现非常高的定位性能，并且随着PSF尺寸的减小，性能不断提高，直到PSF尺寸受到衍射极限的限制。由式（14）估计的最优精度极限如图3（a）中的黑色曲线所示，在此，式（14）中的第一式用于噪声2情况（高噪声水平），第二式用于噪声1情况（低噪声水平）。对于r < 0.3像素的小PSF，黑色曲线与蓝色曲线显示出很好的一致性，验证了该公式描述小尺寸高斯PSF最优精度极限的有效性。

《图3》

图3 仿真结果。（a）平均PL，在像素边界处获得的最优PL，以及式（14）估计的最优PL与PSF半径的关系。仿真测试了两种暗噪声水平。（b）沿着x = 0~1像素和y = 0.5像素的轨迹，高斯半径为0.3像素的PSF的PL估计结果和定位结果，PL的放大视图显示在右侧，其中也显示了式（11）和式（14）的估计结果。（c）高斯半径为0.2像素的PSF的PL估计结果和定位结果。

在仿真图像上测试所开发的MLE定位方法。在此直接将定位精度与精度极限进行比较，因为精度极限虽为精密度极限，但也反映了无偏估计器的精确度。图3（b）和（c）为典型噪声水平（噪声1）下两个目标（r = 0.2像素，r = 0.3像素）在轨迹x = 0~1像素、y = 0.5像素上x*的定位结果。在轨迹上的每个位置定位多帧图像（重复次数n = 500），计算RMS误差以表示精度，利用式（5）计算精度极限（红色曲线）。结果表明，与CG方法相比，本文所提出的基于MLE和ePSF建模的方法大大提高了定位精度，在每个像素相位达到了精度极限。在图3（b）和（c）定位结果的放大图中，显示了式（14）和式（11）第二式的估计结果。结果表明，式（14）准确预测了在像素边界处获得的最优精度极限，式（11）描述了小PSF在像素边界附近的精度极限情况，在该结果中，式（11）在距离像素边界0.25个像素以内有效。更多目标的定位结果如附录A中的第S4节所示，附录中还验证了所提方法在每个像素相位的估计无偏性。

《3.2 实验结果》

3.2 实验结果

在实验室中进行定位实验测试该框架。实验平台如图4（a）所示，使用小PSF光学导航仪器（清华大学皮型星敏感器）进行图像采集，利用高精度三轴转台生成光学仪器与点源之间的相对运动。来自点源的光子经过平行光管，到达固定在转台上的星敏感器的CMOS图像探测器。实验中转台的旋转方向使光斑在x方向上移动，星敏感器的焦距为25 mm，图像探测器的像素大小为5.3 µm × 5.3 µm，转台旋转步长设置为0.0005°，因此图像中光斑的步长约为25 mm × tan(0.0005°)/5.3 µm = 0.041像素。通过多次旋转转台，使光斑的总位移达到10个像素左右。在每个位置采集30帧图像，将一半位置采集的图像用于ePSF建模，另一半位置的图像用于定位测试。经测量，实验中像素背景值_ij为37DN，系统增益K为0.14，暗信号方差 $σ_{d_{i j}}^{2}$ 为54.40，这里忽略了不同像素的这些参数的差异。目标的计算窗口设置为以最亮像素为中心的5×5像素区域，真实位置基准由高精度转台的运动信息确定（如附录A中的第S5节所示）。

《图4》

图4 定位实验平台及结果。（a）实验装置和十分聚焦的PSF图像；（b）使用CG法、GF法、所提方法得到的定位精度结果，并与精度极限估计值及式（11）进行了比较，曲线的RMS值被标记在图上；（c）三种方法通过30次重复测量得到的平均定位误差，右图为所提方法定位精度曲线的放大视图。

首先对一个能量非常集中的光斑进行测试[图4（a）]，利用2.1节中的第一种方法获得ePSF模型，利用式（5）估计每个像素相位的精度极限，通过传统定位方法及所提出的基于MLE和ePSF的方法来定位目标，计算每个位置30次重复测量的定位RMS误差和平均误差。结果表明，与CG和GF方法相比，MLE方法将精度从约0.100像素提高到0.011像素，实现了一个数量级的提升，而精度极限曲线的RMS值为0.010像素[图4（b）]，MLE定位结果与精度极限估计有很好的一致性，为理论精度极限分析赋予了意义和有效性。式（11）第二式的估计结果如图4（b）中的右图所示，结果验证了此简化公式能够有效估计真实PSF的最佳像素相位和最优精度极限。在此实验中，成功预测和实现了小PSF在像素边界附近优于0.005像素的单帧超高定位精度。

此外，与传统方法相比，利用所提出方法的无偏性可以有效地降低定位误差。本方法通过30次重复测量实现了0.004像素的平均误差[图4（c）]。然而，在平均测量结果中像素相位效应并不明显，并且从n次测量平均中获得的定位性能没有提高 $\sqrt[]{n}$ 倍，这主要是因为转台会不可避免地会引入额外的位置误差（附录A中的第S6节），是实验平台而不是所提出的方法阻碍了定位性能的进一步提高。

通过调整平行光管，产生另一个能量较发散的光斑[图5（a）]进行测试。对于该目标，传统方法定位精度达到约0.05像素，MLE方法单次测量的定位误差为0.01像素[图5（b）]，30次重复测量的平均误差为0.003像素[图5（c）]，且MLE结果显示出与精度极限估计的良好一致性，式（11）的有效性也获得了验证。与图4中的目标相比，MLE单帧定位精度尽管仍然为0.01像素，但没有局部精度可以达到图4（b）中的0.005像素。与预期一致，尺寸稍大的目标像素相位效应不明显。

《图5》

图5 另一个PSF的定位实验结果。（a）较为发散的PSF图像；（b）使用CG法、GF法、所提方法得到的定位精度结果，并与精度极限估计值及式（11）进行了比较；（c）三种方法通过30次重复测量得到的平均定位误差。

《3.3 真实星空观测》

3.3 真实星空观测

基于地面观测的真实星空实验常用于测试星敏感器的导航性能[44]。与在轨不同，地面观测受到大气视宁度的影响，恒星图像存在额外的复杂噪声。由于所提框架中不包括大气噪声模型，大气噪声会对定位结果尤其是单次定位结果产生较大的影响。此外，由于地球自转，恒星在图像探测器上的位置不可固定，不可在某位置进行重复测量，也就无法获得该位置处的定位精度统计值，难以与精度极限估计结果进行准确对比。因此在这种情况下，本文仅验证两个方面：①MLE和ePSF建模的结合可以提高星敏感器的定位和导航性能，多个位置的平均定位精度可以接近平均理论精度极限；②恒星的像素相位对定位精度有影响，可以近似估计出低误差区域。在实验中，星敏感器固定在平台上，指向天顶，对真实夜空进行采样，采样间隔为366 ms，在此分析1000帧。利用第二种ePSF建模方法获得每颗恒星的近似ePSF [图6（a）]，ePSF的ROI设置为7×7像素，在29×29个网格点（即Δh为0.25像素）评估ePSF。

《图6》

图6 真实恒星观测结果。（a）实验装置和真实恒星图像。将每帧图像中的同一颗恒星提取出来，建立近似ePSF模型。（b）通过传统方法（CG）和MLE方法重建了两颗恒星的轨迹。（c）一颗恒星在y方向上定位时的精度极限估计结果。根据精度极限大小，像素被划分为三个区域。（d）八颗恒星在三个子像素区域内的x方向（左图）和y方向（右图）定位精度。RLE：低误差区域；RME：中误差区域；RHE：高误差区域。

利用所提出的方法计算图像中恒星的2D位置，其中两颗恒星的轨迹结果如图6（b）所示。结果表明，MLE和ePSF建模结合的定位方法有效地降低了定位随机误差和系统误差。使用二次多项式来拟合x和y方向上的定位结果，并使用残差的RMS来表示平均定位精度。对于这两颗恒星，传统方法（CG）在x方向上的精度分别为0.065和0.084像素，所提方法的精度分别为0.021像素和0.026像素，与0.024像素和0.026像素的平均精度极限估计值一致（其他恒星的结果见附录A中的第S7节），从而大大提高了星敏感器的姿态确定精度（附录A中的第S7节）。

此外，将精度极限估计结果均匀划分为三个级别，相应地，每个像素被划分为三个区域：低误差区域（RLE）、中误差区域（RME）和高误差区域（RHE）。例如，如图6（c）中，给出了一颗恒星在y方向上的精度极限及其区域划分结果。计算所有恒星在三个不同区域的定位性能，其中八颗恒星的结果如图6（d）所示。尽管受到大气的影响，光斑尺寸没有预期小，但RLE的精度仍然优于RHE [RHE和RLE曲线的RMS值在图6（d）中的左图中分别为0.025像素和0.020像素，在右图分别为0.026和0.020像素]，像素相位效应存在，所有恒星的结果如附录A中的第S7节所示。实验表明MLE和ePSF建模的结合极大地提高了星敏感器的定位和导航性能，在具有复杂噪声的真实应用中，验证了在特定像素相位定位目标以获得更好性能的潜力。

《4、讨论》

4、讨论

在亚像素定位中，已有研究已充分阐述了精度极限PL (N, σ_d², r)的优化问题，即在高探测光子数N、低像素暗噪声σ_d²和中等PSF半径r的情况下，精度极限达到最优。在此问题中，目标的像素相位(x, y)通常被忽略，因为它对望远镜或显微镜中使用的典型PSF的精度极限几乎没有影响。然而，对于小PSF系统，目标定位精度极限随像素相位的变化而显著变化，并且在像素边缘附近可以达到非常高的精度，这是因为像素边缘处的像素强度梯度非常大，有利于识别目标微小的位置变化。

本文为指导实际定位应用提出了一种体现像素相位效应的精度极限分析框架。为了准确估计实际应用中的定位精度极限，提供了实验ePSF准确建模方法，并将CRLB应用于ePSF。基于小PSF的特性，导出了描述任意小PSF的像素边界附近最优精度极限的简化公式（11），并在实验室真实PSF上获得了验证。其次，使用典型的高斯PSF进一步分析PL (N, σ_d², r, x, y)的优化问题，通过推导公式（14），发现当高斯半径尽可能小时，PL的最优值在像素边界处实现，而高斯半径又受到衍射效应的限制。虽然亚像素定位方法用于超越衍射极限DL是一种共识，但本文工作揭示，物理衍射极限最终限制了定位精度极限的终极优化结果[式（17）]。最后，应用MLE方法，它与ePSF模型的结合使真实PSF数据成功达到了定位精度极限，为理论精度极限分析赋予了意义和有效性。本工作为亚像素定位理论提供了深入的物理见解，为实际定位应用提供了准确详细的指导，也适用于普通相机且不限于点源，而且它可以扩展到具有锐利图像特征的一般光学目标。

这项工作的主要局限是只在理论和实验上揭示了小PSF的像素相位效应，而使用固定的光学系统并不总能获得最优的定位性能。因此，开发一种具有精确运动调制的图像探测器的光学测量系统，如同基于像素移位的高分辨率相机[45]，是我们未来的研究方向。这项工作提供和论证了将PSF工程与图像探测器位置调制相结合的新型测量原理，为充分利用信息论蕴含的潜力、实现光学测量的最优定位精度极限奠定了重要基础。

Show More