直接比较两个光谱之间的亚像素偏移——并获得可信的误差

信息处理 采样 互相关 平滑
2021-12-26 16:58:43

我有同一个天文物体的两个光谱。基本问题是:如何计算这些光谱之间的相对偏移并获得该偏移的准确误差?

如果您还和我在一起,请提供更多详细信息。每个光谱将是一个具有 x 值(波长)、y 值(通量)和误差的数组。波长偏移将是亚像素。假设像素是规则间隔的,并且只有一个波长偏移应用于整个光谱。所以最终的答案是这样的:0.35 +/- 0.25 像素。

这两个光谱将是许多无特征的连续谱,被一些相当复杂的吸收特征(倾角)打断,这些特征不容易建模(并且不是周期性的)。我想找到一种直接比较两个光谱的方法。

每个人的第一直觉是进行互相关,但是对于亚像素偏移,您将不得不在光谱之间进行插值(通过先平滑?)——而且,错误似乎很难纠正。

我目前的方法是通过与高斯核卷积来平滑数据,然后对平滑结果进行样条曲线化,并比较两个样条曲线光谱——但我不相信它(尤其是错误)。

有谁知道正确执行此操作的方法?

这是一个简短的 python 程序,它将生成两个移动了 0.4 像素的玩具光谱(写在 toy1.ascii 和 toy2.ascii 中),您可以使用它们。尽管这个玩具模型使用了简单的高斯特征,但假设实际数据无法与简单模型拟合。

import numpy as np
import random as ra
import scipy.signal as ss
arraysize = 1000
fluxlevel = 100.0
noise = 2.0
signal_std = 15.0
signal_depth = 40.0
gaussian = lambda x: np.exp(-(mu-x)**2/ (2 * signal_std))
mu = 500.1
np.savetxt('toy1.ascii', zip(np.arange(arraysize), np.array([ra.normalvariate(fluxlevel, noise) for x in range(arraysize)] - gaussian(np.arange(arraysize)) * signal_depth), np.ones(arraysize) * noise))
mu = 500.5
np.savetxt('toy2.ascii', zip(np.arange(arraysize), np.array([ra.normalvariate(fluxlevel, noise) for x in range(arraysize)] - gaussian(np.arange(arraysize)) * signal_depth), np.ones(arraysize) * noise))
1个回答

我认为使用互相关和插值峰值会很好。互相关之前的上采样是否无用中所述?,在互相关之前进行插值或上采样实际上不会为您提供更多信息。关于子样本峰的信息包含在它周围的样本中。您只需要以最小的错误提取它。我在这里收集了一些笔记

最简单的方法是二次/抛物线插值,我在这里有一个 Python 示例。如果您的光谱基于高斯窗口,或者峰值恰好落在样本之间的中点上,那么这应该是准确的,否则会有一些错误所以在你的情况下,你可能想要使用更好的东西。

这是更复杂但更准确的估计器的列表。“在上述方法中,Quinn 的第二个估计器的 RMS 误差最小。”

我不知道数学,但是这篇论文说他们的抛物线插值的理论精度是 FFT bin 宽度的 5%。

在互相关输出上使用 FFT 插值没有任何偏差误差,因此如果您想要非常好的精度,这是最好的。如果您需要平衡精度和计算速度,建议进行一些 FFT 插值,然后使用其他估计器之一来获得“足够好”的结果。

这只是使用抛物线拟合,但如果噪声低,它会输出正确的偏移值:

def parabolic_polyfit(f, x, n):
    a, b, c = polyfit(arange(x-n//2, x+n//2+1), f[x-n//2:x+n//2+1], 2)
    xv = -0.5 * b/a
    yv = a * xv**2 + b * xv + c

    return (xv, yv)

arraysize = 1001
fluxlevel = 100.0
noise = 0.3 # 2.0 is too noisy for sub-sample accuracy
signal_std = 15.0
signal_depth = 40.0
gaussian = lambda x: np.exp(-(mu-x)**2/ (2 * signal_std))
mu = 500.1
a_flux = np.array([ra.normalvariate(fluxlevel, noise) for x in range(arraysize)] - gaussian(np.arange(arraysize)) * signal_depth)
mu = 500.8
b_flux = np.array([ra.normalvariate(fluxlevel, noise) for x in range(arraysize)] - gaussian(np.arange(arraysize)) * signal_depth)

a_flux -= np.mean(a_flux)
b_flux -= np.mean(b_flux)

corr = ss.fftconvolve(b_flux, a_flux[::-1])

peak = np.argmax(corr)
px, py = parabolic_polyfit(corr, peak, 13)

px = px - (len(a_flux) - 1)
print px

在此处输入图像描述 在此处输入图像描述

样本中的噪声产生的结果变化超过整个样本,所以我减少了它。使用更多的峰值点来拟合曲线有助于在一定程度上收紧估计,但我不确定这在统计上是否有效,而且它实际上使低噪声情况下的估计更糟。

在噪声 = 0.2 和 3 点拟合的情况下,对于偏移 = 0.4,它给出的值类似于 0.398 和 0.402。

在噪声 = 2.0 和 13 点拟合的情况下,对于偏移 = 0.4,它给出的值类似于 0.156 和 0.595。