文章目录
-
- Rayleigh统计量
- 时频域Rayleigh分析
- 小结
用gwpy处理引力波系列
Rayleigh统计量
在LIGO引力波数据分析中，瑞利统计量定义为
R
(
f
)
=
σ
(
f
)
μ
(
f
)
R(f)=\\frac{\\sigma(f)}{\\mu(f)}
R(f)=μ(f)σ(f)
式中，
μ
,
σ
\\mu, \\sigma
μ,σ分别是功率谱密度(PSD)的均值和标准差。当
R
≈
1
R\\approx1
R≈1时，说明数据为高斯噪声，符合探测器设计假设。
在gwpy中，【rayleigh_spectrum】是时间序列用于频域 Rayleigh 分析的方法，其具体的计算过程为
R
(
f
)
=
1
K
−
1
∑
k
=
1
K
(
P
k
(
f
)
−
P
~
k
(
f
)
)
2
P
~
(
f
)
,
P
~
(
f
)
=
1
K
∑
k
P
k
(
f
)
R(f)=\\frac{\\sqrt{\\frac{1}{K-1}\\sum^K_{k=1}\\left(P_k(f)-\\tilde P_k(f)\\right)^2}}{\\tilde P(f)}, \\quad \\tilde P(f)=\\frac{1}{K}\\sum_kP_k(f)
R(f)=P~(f)K−11∑k=1K(Pk(f)−P~k(f))2
,P~(f)=K1k∑Pk(f)
计算结果如下
其中，
R
≪
1
R\\ll1
R≪1的部分，一般对应强相干谱线，例如交流电的
60
60
60Hz及其倍频等区域；
R
>
1
R>1
R>1的区域，对应非平稳噪声，例如散射光、机械共振等。而在
R
≈
1
R\\approx1
R≈1的区域，为良好的高斯噪声区域，适合进行引力波搜索。
代码为
from gwpy.timeseries import TimeSeries
# hdata = TimeSeries.get("H1", 1126259446, 1126259478) # 获取GW150914
hdata = TimeSeries.read("gw150914.hdf5")
rayleigh = hdata.rayleigh_spectrum(2, 1)
# 与 ASD 对比绘图
from gwpy.plot import Plot
plot = Plot(hdata.asd(2, 1), rayleigh, geometry=(2, 1), sharex=True)
plot.axes[0].set_yscale('log')
plot.axes[0].set_ylabel('PSD')
plot.axes[1].set_ylim(0, 2) # Rayleigh 值通常在 0-2 范围
plot.axes[1].set_ylabel('Rayleigh statistic')
plot.show()
时频域Rayleigh分析
【rayleigh_spectrogram】有个gram后缀，暗示用于绘制时频图，示例如下，其横坐标为时间，纵坐标为频率，伪彩对应
R
R
R值。
示例代码如下
rayleigh = hdata.rayleigh_spectrogram(5, fftlength=2, overlap=1)
# 可视化
plot = rayleigh.plot(norm='log', vmin=0.25, vmax=4, cmap='coolwarm')
ax = plot.gca()
ax.set_yscale('log')
ax.set_ylim(30, 1500)
ax.set_title('Rayleigh Statistic vs Time and Frequency')
ax.colorbar(label='Rayleigh statistic')
plot.show()
小结
本文介绍了Rayleigh统计量在LIGO引力波数据分析中的应用。该统计量定义为功率谱密度(PSD)的标准差与均值之比，当R≈1时表明数据符合高斯噪声假设。文章详细说明了Rayleigh谱的计算公式，并通过gwpy库的rayleigh_spectrum方法实现分析，指出R≪1和R>1区域分别对应强相干谱线和非平稳噪声。同时介绍了时频域分析方法rayleigh_spectrogram，可生成R值分布的时频图。文末提供了完整的Python代码示例，展示如何计算Rayleigh统计量并进行可视化分析。
