文章目录
-
- gwpy初步
- 简单滤波
gwpy初步
gwpy是一款用于引力波数据处理的Python模块，提供了多种方案，包括conda, pip等，下面用conda安装
conda install -c conda-forge gwpy
安装完成后，可以加载引力波数据，并进行可视化，代码如下
from gwpy.timeseries import TimeSeries
hdata = TimeSeries.get("H1", 1126259446, 1126259478) # 获取GW150914
hdata.plot().show()
其中，
【TimeSeries.get】用于获取引力波数据，其输入的三个参数分别代表引力波探测器和起止GPS时间，H1代表的是Ligo在汉福德的引力波干涉仪。
【TimeSeries】是gwpy的主要数据类型，其内部封装了大量的数据处理和可视化方法。
在上述代码中，通过【plot】绘制hdata中的数据，并调用【show】弹出图像窗口，结果如下
其中横坐标为时间，单位是秒，这段数据从2015年9月14日的9:50:29开始，总计33秒。其纵坐标为应变，代表的是空间尺度变化的百分比，无量纲。这段原始数据发现不出任何问题，几乎和噪声没有区别。原因在于，引力波的强度仅有
10
−
21
10^{-21}
10−21左右，已经淹没在了噪声中，为了提取处数据，需要进行滤波操作。
简单滤波
【TimeSeries】中封装了许多便捷的数据处理操作，下面对其进行双边滤波，效果如下。
为了便于观察，这里将信号从9:50:44截取，可以看到9:50:45.4出现了一个疑似信号的东西，这是是人类有史以来第一次观测到引力波。LIGO的干涉仪臂长为4km，光在F-P腔内往返约300次，有效臂长约为1200km，则
10
−
22
10^{-22}
10−22表示
1.2
×
10
−
16
1.2\\times10^{-16}
1.2×10−16m，也就是
0.12
0.12
0.12fm，作为参考，氢原子半径是
53
53
53pm，即
53000
53000
53000fm，也就是说LIGO观测到了比氢原子还要小6个数量级的尺度变化。
滤波与可视化代码如下。
filtered = hdata.bandpass(50, 250).notch(60).notch(120)
plot = filtered.plot(
xlim=(1126259461, 1126259463),
ylim=(–1e-21, 1e-21),
)
hdata.bandpass(50, 250).plot().show()
plot.show()
【bandpass】为带通滤波，用于保留
50
→
250
50\\to250
50→250Hz之间的频段，这个频段是黑洞和中子星合并信号的主频段。
【notch】为陷波滤波，用于剔除某些专门的频段，上述代码中，通过两次陷波滤波，剔除了60Hz和120Hz这两个频段的噪声。其中，60Hz是交流电频率，120Hz为其二次谐波。
在诸多噪声中，
<
50
<50
<50Hz的部分最需剔除，因为这些噪声主要源于地震噪声、悬挂系统共振、地面微震等，量级较大，在剔除这些低频噪声后，数量级会降至
10
−
22
10^{-22}
10−22附近。然后再剔除
250
250
250Hz以上的高频噪声，就可以看到突兀的引力波信号了。后续对电网噪声的剔除，则让引力波信号更加清晰。
仅筛选出
50
250
50~250
50 250Hz数据的代码如下
hdata.bandpass(50, 250).plot().show()
滤波结果为 
