前言

在我们的实验中, 我们是先利用位移的估计子绘出时域信号, 再利用算法来估计波形中的参数, 其中我们主要关心的是频率. 在实验上, 我们是先利用位移的估计子来得到波形, 再从波形得到频率. 我们可以采用最小二乘估计来实现这个实验. 取函数最小值又可以利用梯度下降算法来迭代, 直到达到我们要的目标精度.

另外, 在文献中提到, 简单的求平均的估计子是一个不好的估计子, 容易受到噪声影响. 由于我们需要用位移的估计子来获得时域信号, 同时又由于课题本质上是一个性能对比实验, 这要求 SPADE 和 DI 的估计子都是最好的. 这样 DI 就不能选取简单的求平均的估计子来对比. 所以我使用梯度下降算法实现了 DI 的最大似然估计, 结果发现最大似然估计子能够在有噪声的情况下也给出非常准确的估计值, 但缺点是算法复杂度高, 处理比较慢.

估计波形中的参数

在 Fundamantals of Statistical Signal Processing 书中提到, 一个采集到的信号可以分解为以下形式:

其中 是发送者发送的真实信号, 是采集过程中被带入的噪声. 为了简单起见, 在这里只讨论 是一个零均值的高斯白噪声 (white Gaussian noise, WGN) 的情况, 同时 .

波形中的参数的 CRB

假设信号是依赖于参数 的, 信号就可以写为如下形式:

因为 是 WGN, 可以把 看作是一个均值随着时间变化的随机变量 . 这样似然函数就可以写为:

取对数后求导:

求二阶导:

取数学期望后就得到了经典 Fisher 信息:

例: 计算正弦波参数估计的 CRB

假设采集到的信号为:

书中提到了一组近似并给出了参考文献:

计算要利用前面推导得到的 CFI 计算式:

其中导数容易计算得到:

这里取了 . 综上 CFI 就等于:

我理解中, 如果要应用这个结论, 首先得把频率估计问题拆成两步. 把频率估计的过程考虑为估计了两次. 第一次估计是用位移的估计子得到波形 (时域信号), 第二次估计是用波形估计波形中的频率参数. 这两次估计的过程是分离的, 也就是我可以把第一次估计的过程看作是一个黑匣子模型. 此时黑匣子的特性就决定了波形的特性. 我们知道 DI 和 SPADE 关于位移的估计子都是 MVU 估计子, 每次估计位移都是独立的. 同时因为估计子的渐进正态性, 如果我们暂时不考虑 SPADE 在位移不为 0 的时候 , 如果把波形拆成信号和噪声两项, DI 和 SPADE 噪声会是一样的, 而且都是高斯白噪声.

在实验上, 我们可以令点光源的位移为零, 利用位移的估计子估计出时域信号. 因为点光源的位移为 0, 即信号为 0, 采集到的时域信号可以认为完全是噪声过程. 对这个信号进行 FT 分析并做直方统计, 可以发现这个噪声过程是恰好是一个零均值高斯白噪声, 即便是引入了额外噪声也是如此.

最小二乘估计

对于估计波形中的参数而言, 同样假设信号为:

最小二乘估计 (least squares, LS) 的定义为:

这个估计的方法要求我们要对真实的信号模型有一个比较完整的认知. 另外, 当 是 WGN 的时候, 应用最大似然估计 (maximum likelihood, ML), 考虑负的对数似然:

除去一些不影响极大值选取的常数, MLE (最大似然估计子) 就等价于选取下式最小的 值:

形式上就等同于 LSE (最小二乘估计子). 所以当 是 WGN 时候, LSE 是最小方差无偏 (minimun variance unbiased, MVU) 估计子.

例: 正弦信号频率估计问题

假设信号模型为:

其中频率 是待估参数, 可以写出:

只需要让这个函数最小, 解出 . 注意到在这个例子里面, LS 误差是 的高度非线性函数, 往往无法解出解析解. 所以我们一般需要用到梯度下降法等算法来进行求解.

梯度下降算法

为了解决用穷举法求函数最小值和最小值点算法的复杂度过高的情况, 一般采用梯度下降 (gradient descent, GD) 算法来降低算法的复杂度, 从而将算力用于提高精度. 这是机器学习领域最常用的优化方法, 是一个成熟的算法.

这个算法的做法是人为给出一个步长 , 在机器学习领域常被称作学习率 (learning rate), 计算函数在当前位置的梯度并用该位置的梯度和步长的乘积来更新当前位置, 用公式写如下形式:

计算上式, 再用 替代 , 如此重复, 直到 满足所设置的精度要求或者达到最大的迭代次数. 拓展到多个参数的情况也是简单的, 只需要对每个参数求偏导, 再分别计算 GD 即可.

另外, 在应用最小二乘法的时候, 如果数据点较多 ( 比较大), 需要计算大量数据的求和, 这会消耗太多计算机的内存, 所以就有随机梯度下降 (SGD) 算法和小批量 SGD 算法. 思路是每次只随机选取一部分数据点用于计算 LS 误差并使用 GD 算法进行最优化. 由于目前为止我们应该不涉及大量数据的处理, 所以应该用不上 SGD 算法.

DI 的不同位移估计子对比

在 Single Particle Tracking: From Theory to Biophysical Applications 文章中提到: 简单的求平均的估计子 (后续简称简单估计子) 极易受到噪声的影响, 一般常用的做法是最小二乘估计和最大似然估计来确定点光源的位置. 但是由于后者两个方法一般无法解出解析解, 在计算性能开销方面会远远大于简单估计子. 这使得简单估计子在很多时候仍然有效.

同样在文献中提出了三种不同的估计子, 分别是 LSE、MLE 和一个计算对称性的估计子. 其中第三个估计子是用于降低算法复杂度的做法, 原理应该类似霍夫变换找圆心. 由于我们期望估计子是最好的, 所以我认为直接选取 MLE 即可. 但是文中也提到了 LSE 则对噪声模型的要求更低, 同时对点扩展函数的不完全匹配容忍度更高. 在后续如果有需要也可以尝试使用 LSE 进行估计.

点光源定位问题中的 MLE

无噪声的情况

在不考虑噪声的情况下, LSE、MLE 和简单估计子将会是一样的. 考虑光子落在像平面上位置的 PDF:

由于后续要计算对数似然, 正的常数不会影响最值的选取, 所以这里简记为 A. 继续求取对数连续似然:

去掉不影响最值选取的常数并取负, MLE 就可以写为:

形式上和 LSE 是一样的. 同时如果要求解这个方程, 只需要对 求导后令导数为 0 即可, 计算可以得到:

这和简单估计子是完全一致的, 所以在不考虑噪声情况下三个估计子是完全一样的.

在有噪声的环境

考虑有噪声情况下的 PDF:

其中, 表示噪声, 是归一化常数, 具体的取值和噪声强度以及信号强度有关. 假设噪声在 ROI 范围内是服从均匀分布的, 所以 . 同时假设 ROI 足够大, 远大于 的范围. 这样归一化条件就可以写为:

其中 表示 ROI. 所以 PDF 就可以写为:

其中 . 这其实和常见的处理方法是一样的. 实验上, 由于我们采集到的数据已经是信号和光子叠加在一起的了, 噪声的强度一般要通过先验的信息来确定, 比如仿照《光子与神经元》一书中的做法, 取一个先验知道没有点光源的区域, 计算多个样本在该区域内的平均值 (后续会用这个办法来作为噪声的真值). 我的处理方法是将噪声强度也作为一个待估参数, 利用 SD 算法来同时算出噪声强度和位移参数 . 所以对 做单位上的变换:

其中 表示样本的每个像素点上的平均光子数, 在实验上是好测的, 作为一个已知的参数处理. 而 表示噪声的每个像素点上的平均光子数, 作为一个待估计参数来处理. 下图是使用实验数据估计出来的背景噪声的强度. 而噪声强度的真值是由一组由实验上采集到的, 无信号的纯背景数据计算得到的. 由于我们引入噪声的方法是用一个恒压源给小灯泡供电, 所以小灯泡的亮度 (也就是引入噪声的强度) 是可以准确控制的.

.assets/output2.svg)

可以看到噪声的估计值非常准确, 标准差也非常小. 或许后续可以考虑把噪声作为是一个先验的参数放到 PDF 里面, 如果这样做算法的复杂度将得到大幅减小, 额外的算力可以用于提高对位移参数估计的精度.

以一个实验上的样本为例子, 可以看到连续对数似然的极大值点和真值点符合得相当好, 这可以初步看出 MLE 能够给出比简单估计子准确得多的估计值.

噪声强度不变, 不同位移情况下的 MLE 和简单估计子的对比

在这里使用了一组无噪声情况下的实验数据作为参考对比. 在整个实验过程中, 激光功率不变. 在这个激光功率下, DI 信号的平均光子数约为 550 个左右.

.assets/output1.svg)

可以看到简单估计子受到背景噪声的巨大影响, 使其完全无法分辨点光源的位移. 但是注意到简单估计子的估计值方差很小, 这是因为简单估计值完全是由背景噪声决定的, 而背景噪声在实验中是不变的, 所以会使得估计值的方差变小.

位移固定, 不同噪声强度下的 MLE 和简单估计子的对比

这组数据是让位移不变, 通过改变给小灯泡的电压来改变引入噪声的强度. 在不同的噪声强度下, 可以看到 MLE 的估计值虽然方差变大了, 但是不太影响估计值的偏差.

.assets/output3.svg)

总结

由上面的一系列实验数据可以看出, MLE 在有噪声的情况下确实比简单估计子要好得多. 但是由于 MLE 的计算量大, GD 算法调参费时费力, 简单估计子还是有用武之地的. 尤其是在无噪声的条件下, 简单估计子和 MLE 是一样的, 在这种情况下使用简单估计子会方便得多. 但是如果在背景噪声强度比较大的情况下, 简单估计子会因为容易受到噪声影响而无法使用.

附: SPADE 和 DI 关于频率的 CFI

SPADE

由于 , 所以我认为可以直接将 作为待估参数来考虑.

其中:

所以有:

取对数:

记 , 求一阶偏导:

求二阶偏导:

求负期望得到经典 Fisher 信息:

其中:

利用 PMF 的归一化条件, 求和中只有当 有需要单独考虑外, 其他都为 1, 连乘后就消去了.

考虑第一项, 同样利用 PMF 的归一化条件:

类似的, 考虑第二项:

所以:

注意到, 这和前面从波形出发计算的 CFI 形式是一样的, 但是 是不同的. 前面的 是噪声的标准差, 而这里的 是特征宽度. 这就涉及到我们计算 CFI 的时候的假设: 假设光子是一个一个到达的, 每个样本中只有一个光子. 如果考虑这种情况, 从波形中的全部噪声都来源于测量的不确定性, 此时从波形出发的 CFI 中的 就等于 QCRB, 也就是特征宽度. 如此, 这两个 CFI 的表达式就是完全等价的了.

DI

求对数导数:
$$
\pdv{\ln p(\vb x;\omega)}{\omega}=\sum^{N-1}_{n=0}\qty[\frac{2(x-s)\partial_\omega s}{2\sigma^2}]=-\frac A{\sigma^2}\sum_{n=0}^{N-1}n(x-A\cos\alpha)\sin\alpha
$$
求二阶导:
$$
\begin{align}
\pdv[2]{\ln p(\vb x;\omega)}{\omega}&=-\frac A{\sigma^2}\sum_{n=0}^{N-1}n^2(x\cos\alpha+A\sin^2\alpha-A\cos^2\alpha)\\
&=-\frac A{\sigma^2}\sum_{n=0}^{N-1}(n^2x\cos\alpha+n^2A\cos2\alpha)

\end{align}
$$

取负期望得到经典 Fisher 信息:

积分第一项是求 的期望的形式, 第二项是 PDF 在全空间内积分的形式