用gwpy处理引力波数据

文章目录gwpy初步简单滤波gwpy初步gwpy是一款用于引力波数据处理的Python模块提供了多种方案包括conda, pip等下面用conda安装conda install -c conda-forge gwpy安装完成后可以加载引力波数据并进行可视化代码如下fromgwpy.timeseriesimportTimeSeries hdataTimeSeries.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个数量级的尺度变化。滤波与可视化代码如下。filteredhdata.bandpass(50,250).notch(60).notch(120)plotfiltered.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 5050Hz的部分最需剔除因为这些噪声主要源于地震噪声、悬挂系统共振、地面微震等量级较大在剔除这些低频噪声后数量级会降至10 − 22 10^{-22}10−22附近。然后再剔除250 250250Hz以上的高频噪声就可以看到突兀的引力波信号了。后续对电网噪声的剔除则让引力波信号更加清晰。仅筛选出50 250 50~25050250Hz数据的代码如下hdata.bandpass(50,250).plot().show()滤波结果为