用gwpy处理引力波数据
文章目录
- gwpy初步
- 简单滤波
gwpy初步
gwpy是一款用于引力波数据处理的Python模块,提供了多种方案,包括conda, pip等,下面用conda安装
conda install -c conda-forge gwpy安装完成后,可以加载引力波数据,并进行可视化,代码如下
fromgwpy.timeseriesimportTimeSeries hdata=TimeSeries.get("H1",1126259446,1126259478)# 获取GW150914hdata.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.120.12fm,作为参考,氢原子半径是53 5353pm,即53000 5300053000fm,也就是说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\to25050→250Hz之间的频段,这个频段是黑洞和中子星合并信号的主频段。
【notch】为陷波滤波,用于剔除某些专门的频段,上述代码中,通过两次陷波滤波,剔除了60Hz和120Hz这两个频段的噪声。其中,60Hz是交流电频率,120Hz为其二次谐波。
在诸多噪声中,< 50 <50<50Hz的部分最需剔除,因为这些噪声主要源于地震噪声、悬挂系统共振、地面微震等,量级较大,在剔除这些低频噪声后,数量级会降至10 − 22 10^{-22}10−22附近。然后再剔除250 250250Hz以上的高频噪声,就可以看到突兀的引力波信号了。后续对电网噪声的剔除,则让引力波信号更加清晰。
仅筛选出50 250 50~25050250Hz数据的代码如下
hdata.bandpass(50,250).plot().show()滤波结果为
