引力波天文V：数据分析实践

Posted on 2023-09-01 In Learning , Project Views: Waline:

引力波天文系列文章将介绍引力波能观测到什么，如何观测，以及基本的数据分析方法。侧重于引力波的空间探测和分析，但核心的测量原理、仪器响应以及数据分析方法基本都是通用的。

研究目标：波源概述(基本波源及频谱划分)
探测技术：响应与灵敏度、空间探测技术
数据分析：数据分析方法、数据分析实践

软件

GWOSC：GWOSC数据接口
GWpy：GW数据分析处理
pyCBC：CBC信号探测与参数估计
cWB：基于Wavelets的信号探测
Bilby：贝叶斯参数估计
BayesWave：基于Wavelets的glitch分析
Omicron：基于Q变换的数据分析
ligo.skymap：天图可视化
GWDetChar 引力波探测器

上手

在线性近似下，信号的探测本质上是一个线性变换，而探测器可以视为一个线性时不变系统，而线性时不变系统可表示对信号的卷积操作。

探测器对信号的线性响应对应于卷积操作，同时不可避免地会引入噪声，信号为 $h(t)$ 经探测器响应后变为 $\bar{s}(t) = h(t) * r(t) + \bar{n}(t)$ ，对应到频域域 $\bar{S}(f) = S_h(f)R(f) + S_n(f)$ ，其中 $S_h(f)$ 为信号频谱， $S_n(f)$ 为仪器噪声频谱(噪声曲线)，而 $R(f)$ 为探测器响应函数， $\bar{S}(f)$ 为观测数据频谱。考虑 $R(f)$ 及 $S_n(f)$ 都与仪器相关，除以 $R(f)$ 定义信号的观测频谱 $S(f) = S_h(f) + S_n(f)/R(f)$ ，其中 $S_n(f)/R(f)$ 被称为探测器的灵敏度曲线，对应到时域 $s(t) = h(t) + n(t)$ 。
示意图

h(t) \rightarrow h(t) * r(t) + n(t)

Signal Detection https://www.elenacuoco.com/2020/04/13/lecture-i-signal-detection/

重采样
https://ww2.mathworks.cn/help/signal/ug/resampling-uniformly-sampled-signals.html

探测：匹配滤波

匹配滤波

信号白化后内积

⟨u(t)|v(t)⟩ = \int u(t)v(t) dt = \int \overline{U(f)}V(f) df

⟨u(t)|v(t)⟩ = \int \frac{\overline{U_x(f)}V_x(f)}{S_n(f)} df

尤其是在信号弱噪声强时，模板白化后近似白噪声
$|\hat{X}_{\frak t}| = \left| \frac{\overline{X_{\frak t}(f)}}{\sqrt{S_n(f)}} \right| \sim T$ ，从而 $E = \int T df = T f_s = N$ ！！！

从相关角度：
数据和模板分别白化后相关

y(t) = \mathcal{F}^{-1}\left\{\frac{X(f)}{\sqrt{S_n(f)}} \frac{\overline{X}_s(f)}{\sqrt{S_n(f)}}\right\}

具体数值受模板能量影响

\sqrt{\int \left|\frac{\overline{X}_s(f)}{\sqrt{S_n(f)}} \right|^2 df}

对(白化后)模板能量归一化

\hat{y}(t) = \mathcal{F}^{-1}\left\{\frac{X(f)\overline{X}_s(f)}{S_n(f)}\right\} / \sqrt{\int \frac{|X_s(f)|^2}{S_n(f)} df}

从滤波角度：
数据滤波

y(t) = \mathcal{F}^{-1}\left\{X(f)H(f)\right\}, H(f) = \frac{\overline{X}_s(f)}{S_n(f)}

SNR = \frac{E\{y(t)\}}{\sqrt{E\{y^2_n(t)\}}} = \frac{\mathcal{F}^{-1}\left\{X(f)H(f)\right\} }{ \sqrt{ \int S_n(f)|H(f)|^2 df}} = \mathcal{F}^{-1}\left\{\frac{X(f)\overline{X}_s(f)}{S_n(f)}\right\} / \sqrt{\int \frac{|X_s(f)|^2}{S_n(f)} df}

可以看到白化后模板的能量与滤波后噪声功率是对应的，进而归一化后的相关函数与滤波后信噪比对应。简单思考，滤波与白化后相关对应，不考虑模板的影响，数据在白化后，噪声功率谱密度为1，因此滤波后噪声功率谱密度就对应白化后的模板，噪声功率就对应(白化)模板能量。

y(t) = \mathcal{F}^{-1}\left\{\frac{X(f)\overline{X}_s(f)}{S_n(f)}\right\} = \int \frac{X(f)\overline{X}_s(f)}{S_n(f)} e^{i2\pi ft} df

在信噪比最大时， $x(t) = s(t) + n(t), X(f) = kX_s(f)e^{-i2\pi ft}$ ，此时

E\{y(t)\} = \int \frac{k X_s(f)\overline{X}_s(f)}{S_n(f)} df = kE\{y^2_n(t)\}

从而信噪比最大值为 $k\sqrt{E\{y^2_n(t)\}}$ 。

关于白化
用自身功率谱白化，白化后功率谱接近1，这也是白化的本意。数据和模板频率用噪声功率谱加权严格说并非“白化”，得到的频谱仍可能存在较大波动，不过实际中通常不会严格。

注意：这里白化时用的是单边谱，理论上应除以2，否则白化后PSD是1/2，考虑到白化后PSD同样使用单边谱，因此最终还是1。

白化后信号在频域的谱密度为1，而非时域的震荡方差(功率)为1！功率谱密度为1，对应总功率为 $f_s$ ，因此功率为1需要将白化信号除以 $\sqrt{f_s}$ 。直接用单边谱白化时，谱密度为 $1/2$ ，总功率 $f_s/2$ ，需要除 $\sqrt{f_s/2}$ ，下面会看到这个归一化因子。

最后，如果要功率与原信号相当还要再乘以原信号标准差，即乘以 $\sigma/\sqrt{f_s}$ 。

计算等效距离 $D_{\rm eff}$

等效距离是实际距离与探测器响应的综合，考虑了波源相对探测器的位置和朝向，具体可参考。引力波响应幅度与等效距离成反比( $\propto \frac{1}{D_{\rm eff}}$ )：

h(f) = \frac{1 {\rm Mpc}}{D_{\rm eff}} A_{\rm 1 Mpc}(M, \mu) f^{-7/6} e^{-i\Phi(f; M, \mu)}

而在生成模板时，等效距离通常取为1Mpc，即：

h_{\rm template}(f) = A_{1 \rm Mpc}(M, \mu) f^{-7/6} e^{-i\Phi(f; M, \mu)}

y(t) = \mathcal{F}^{-1}\left\{\frac{X(f)\overline{X}_s(f)}{S_n(f)}\right\} = \int \frac{X(f)\overline{X}_s(f)}{S_n(f)} e^{i2\pi ft} df

$y(t) \propto X(f) \propto h(f) \propto \frac{1}{D_{\rm eff}}$ ，在信噪比最大时 $E\{y\}_\max = \frac{1\rm Mpc}{D_{\rm eff}}\sigma^2$ 。
从而

{\rm SNR}_m = \frac{ E\{y\}_\max }{\sigma} \propto \frac{1 \rm Mpc}{D_{\rm eff}} \sigma, ~~~ D_{\rm eff} \sim \frac{\sigma}{\rm SNR} {\rm Mpc}

相同强度的模板信号，距离越远信噪比越低，取信噪比阈值为8，对应特定引力波系统的极限观测距离为 $\frac{\sigma}{8}$

Q变换

引力波数据处理胡一鸣第五章

白化：观测信号数据除以数据的振幅谱密度

DFT v.s. CWT v.s. CQT

经过滤波白化单峰会变成震荡(谱泄露？)
拟合优度goodness of fit
PRL 116, 061102(2016)
pastro
p值完全只考虑噪声
p_astro同时考虑信号事件率模型、仪器选择效应(蒙卡仿真)

基于天体族群演化模型(population model)的事件率以及仪器选择效应
数据标定、噪声的物理机制
源自信号的概率 v.s. 源自噪声的概率

对数据质量的要求，影响考虑的具体数据

“trigger set” gravitational wave false alarm rate time shift
cwb
辐射引力波，向内吸收，向外耗散，阻尼震荡(准正则模式)
lm 空间几何角分布(角模式) (l,m,n) = (2,2,0)
n 时间衰减模式(泛音)
黑洞谱
根据质量自旋只能确定不同模式的频率，无法确定不同模式的幅值，
难点：初始时间的确定；
双黑洞(旋进并和)与单黑洞(铃宕)一致性分析
Savage-Dickey density ratio
bayesian factor
只保留铃阶段，加窗后填补数据，最小化对协方差矩阵的影响
质量、自旋、不同模式的初始相位和幅值为自由参数
检验相对论：不同时间分析的一致性、不同模式分析的一致性
4s