LIGO-Virgo引力波短时信号提取指引

简介

引力波探测是探索宇宙奥秘的重要途径,在这个领域LIGO和Virgo(VGO)两个合作组织(LIGO Scientific Collaboration and the Virgo Collaboration, LVC) ,从代号为GW150914的黑洞合并事件开始,完成了一系列事件的探测。全球的引力波探测器包含位于美国的两个LIGO设备,一个位于意大利的Virgo设备,以及位于德国、日本、印度等国家的设备。

现今为止,LIGO探测器有O1,O2,O3三轮观测,其中Virgo探测器在O2时候加入探测,并与LIGO探测器一同进行O3的探测。多个探测器协同观测,可以有效地过滤掉设备噪音,并进行协同分析。LIGO-Virgo探测的数据通过引力波开放科学中心(Gravitational-Wave Open Science Center, GWOSC)延时数月公布,供全球的研究者共同研究。

LIGO-Virgo数据的特性

LIGO-Virgo探测器本质上都是大型的迈克尔逊干涉仪(Michelson interferometers),当引力波通过探测器,导致探测器臂长发生改变,而两个臂长的长度差异进而体现在探测器的输出数据上。探测器的实际输出(interferometer photodiode output)经过一系列标准化的程序,最终生成时序的引力波数据(gravitational-wave strain data)。对于LIGO来说,最终数据采样率16384Hz,而对于Virgo探测器,最终采样率20kHz。而处理程序也都存在一个有效频率(validity range)的东西,这里就不赘述了。

在这些探测器里面,还有一些辅助通道(auxiliary channels),这些通道中的时序数据,记录了探测器、周围环境中的一些辅助信息,这些数据通过GWOSC的处理,刻画了包括噪声源的若干数据,并且其中由于仪器故障、调试导致的数据错误亦都有所标注。

检测器噪音的特征

LIGO-Virgo仪器的生成数据事实上包含了很多种噪音,包括——

  • 量子传感噪音(quantum sensing noise)
  • 地质运动噪音(seismic noise)
  • 热力学波动噪音(suspension thermal noise)
  • 反射镜热力学噪音(mirror coating thermal noise)
  • 重力梯度噪音(gravity gradient noise)

在此之外,还包含部分人为噪音、天气噪音、仪器检修故障以及未知的噪声源。

下面,我们从数学得到角度,形式化的定义一下检测器的噪音输出——

我们在一个时间窗口中探究探测器的噪音,令探测器输出为n(t)为一个向量,n_i = n(t_i)表示时刻t_i的输出值。可以发现,探测器的输出实际上是被噪音主导的,而我们可以将噪音理解为若干个高斯分布的分量的组合。

对于探测器输出n来说,我们可以将其看做一个随机过程的样本,而这个随机过程本身可以用其概率密度函数p(n)描述。而p(n)本身又可以分解为均值矩阵\mu = E[n]以及协方差矩阵C_{ij} = E[(n_i - \mu)(n_j - \mu)]

当然,对于参数估计来说,我们有如下公式——

    \[\hat{\mu} = \frac{1}{N} \sum_{i=1}^{N} n_i\]

    \[\hat{C_{ij}} = \frac{1}{M-1} (n_i - \hat{\mu})(n_j - \hat{\mu})\]

有了上面的定义,我们可以通过高维正态分布的公式获得概率密度函数p(n)

    \[p(n) = \frac{1}{det(2\pi C)^{1/2}} exp\left[-\frac{1}{2}\sum_{ij} (n_i-\mu)(n_j-\mu)C^{-1}_{ij}\right]\]

我们发现,当我们将所有数据看做一个样本的话,那么N足够大,而M=1,以致于我们甚至无法定义有效的\hat{C},这个问题我们可以加入一个“噪音为稳定(stationary)的”这个限制条件来解决。

一个噪音稳定当且仅当C_{ij}只与时间间隔|i-j|有关。稳定噪音可以使用相关性矩阵C(\tau)描述,其中\tau为时间间隔。当然,既然这个过程和波有关,人们习惯上都会为了降低复杂性,将它变换到傅里叶空间中,自此,我们将i,j修改为频域f_i, f_j,将稳定噪音的相关性矩阵设置为C_{ij} = \delta_{ij}S_n{f_i},而这里面的S_n{f_i}就是功率谱(power spectral density)函数。而振幅谱(Amplitude spectral density)函数则是功率谱函数的平方根,以Hz^{-1/2}为单位。

在上面的定义中,如果C_{ij} = \delta_{ij}\sigma^2,那么噪音为白噪声,毕竟每一个时间分量独立且方差皆为\sigma^2。白噪声是LIGO-Virgo探测器噪声的一个不那么准确的近似。

LIGO-Virgo探测器输出数据,在时域和频域上面,都有非常显著的结构特征。而就引力波事件GW170817分析,噪音几乎主导了数据信噪比。这使得对噪音的分析正确与否直接决定了检测的正确性。

使用傅里叶变换进行噪音分析

由于对于噪声稳定的假设,噪声在每一个频率区间是假设独立的。对于大部分LVC分析来说,都首先使用FFT将数据进行一次变换。由于FFT假设数据是循环的,我们需要使用窗函数对FFT后数据的频域进行修正,这里我们使用土耳其窗(Turkey-Window)函数进行修正。

接下来,我们将FFT后的频谱强度除以噪声的频谱强度,进行标准化,以便让噪声强的区域和噪声弱的区域拥有近似相等的影响度。在可视化的时候,数据还要额外经过一个[35Hz, 350Hz]的滤波器,以过滤区间外的噪音,人为保留重要的区间供研究人员观察。

    \[d(t) \xrightarrow{FFT} \bar{d}(f(t)) \xrightarrow{Whiten} \bar{d}(f) = \frac{\bar{d}(f)}{S_n^{1/2} (f)} \xrightarrow{iFFT} d_w(t)\]

噪声频谱强度S_n(f)我们事先并不知道,需要从数据本身估计出来。一般来说,我们可以在待检测信号周围取一个足够大的时间范围,并使用该范围内的数据流进行计算。不过这样的计算有一个问题,就是对于每一个频率,只会算出一个复数(实部+虚部)强度。因此,对于不同的频率位置,估计值的偏差会比较大。

对于上述问题,有两种可能的解决办法:1)应用某种平均的方法降低误差;2)使用某种物理学中更加准确的算法,如Welch方法

上图是一个信号序列,一次经过一个Tukey窗函数,除以噪声强度谱,并经过[35Hz, 350Hz]的滤波器后,揭露了引力波信号GW150914的存在。

参考文献:A guide to LIGO-Virgo detector noise and extraction of transient gravitational-wave signals

使用Bartlett和Welch法进行频谱估计

在频谱分析中,我们通常希望能够通过给定的波形数据流,拟合出不同的频率的强度。而拟合的过程用的是离散傅里叶变换。

上述的拟合有一个小问题,就是实验发现,上述频谱分析算法,给每一个频率计算了一个强度(一个实部+一个虚部)。实际上收到数据噪音的影响比较大。一个直观的解决方法是使用某种平均来减小估计误差。

Barlett算法步骤为:

  • 将最初的N个数据点均分为不相交的K段,每段长度M
  • 对于每一段,使用离散傅里叶变换(DFT)计算振幅谱,并平方除M
  • 对于上一步得到的K个结果,取平均来减小误差。

Welch算法步骤为:

  • 原信号被分割为L个互相重叠的段,段长M,重叠部分长度为D
    • 如果D=M/2,那么重叠部分就占50%
    • 如果D=0,那么所有段不重叠,退化为Barlett算法。
  • 每一个单独的数据段,先过某一个窗函数。(窗函数通常使得数据段中心的权值更高,故使用互相重叠的段,保证不会丢失边缘数据点的信息)。
  • 使用将每段数据的离散傅里叶变换结果,平方,求均值,最终获得 能力-频率 的谱函数。

参考链接:

Welch’s method:https://en.wikipedia.org/wiki/Welch%27s_method
Barlett’s method:https://en.wikipedia.org/wiki/Bartlett%27s_method