多波段数据的Lomb-Scarger周期图#
Lomb-Scarger周期图(在Lomb之后 [1], 和斯卡格尔 [2]) 是一种常用的统计工具，用于检测不均匀间隔观测中的周期信号。基地 LombScargle 为Lomb-Scarger周期图的多个实现提供接口。然而， LombScargle 仅处理单波段数据。这个 LombScargleMultiband 类调整此接口以处理多波段数据(其中存在多个波段/滤镜)。
这里的代码改编自 astroml 包装 ([3], [4]) 以及 gatspy 包装 ([5], [6]) 中建立的设计范例，但与 LombScargle 。有关指导本课程开发的多波段Lomb-Scarger周期图的详细实践讨论，请参见 Periodograms for Multiband Astronomical Time Series [6].
基本用法#
备注
如中所示 LombScargle 、频率中的 LombScargleMultiband 是 not 角度频率，而不是振荡频率(即每单位时间的周期数)。
Lomb-Scarger多频段周期图旨在检测存在多个频段数据的不均匀间隔观测中的周期信号。
例子#
要在不均匀分布的观测中检测周期信号，请考虑以下多频段数据，其中5个频段(u、g、r、i和z)各有60个数据点。
>>> import numpy as np
>>> t = []
>>> y = []
>>> bands = []
>>> dy = []
>>> N=60
>>> for i, band in enumerate(['u','g','r','i','z']):
... rng = np.random.default_rng(i)
... t_band = 300 * rng.random(N)
... y_band = 3 + 2 * np.sin(2 * np.pi * t_band)
... dy_band = 0.01 * (0.5 + rng.random(N)) # uncertainties
... y_band += dy_band * rng.standard_normal(N)
... t += list(t_band)
... y += list(y_band)
... dy += list(dy_band)
... bands += [band] * N
在根据输入数据自动选择的频率上计算的Lomb-Scarger周期图可以使用以下方法计算 LombScargleMultiband 类，并使用 bands 参数是与 LombScargle 接口：
>>> from astropy.timeseries import LombScargleMultiband
>>> frequency,power = LombScargleMultiband(t, y, bands, dy).autopower()
使用Matplotlib绘制结果可以得到：
周期图显示了一个明显的峰值，频率为每单位时间1个周期，正如我们从构建的数据中所预期的那样。产生的功率是单个阵列，每个频段的组合输入取决于在 method 关键字。
周期图来自 TimeSeries 对象#
LombScargleMultiband 能够进行手术 TimeSeries 对象，但前提是 TimeSeries 对象满足格式要求。要求在单独的列中提供每个波段的通量(或幅度)和误差。如果相反，您的 TimeSeries 对象具有具有关联的波段标签列的单个通量列，则这些列可以直接传递到 LombScargleMultiband 作为一维阵列。
例子#
考虑下面的生成器代码 TimeSeries 填充三个光度学波段(g、r、i)的时间序列数据的对象。
>>> from astropy.timeseries import LombScargleMultiband, TimeSeries
>>> from astropy.table import MaskedColumn
>>> import numpy as np
>>> import astropy.units as u
>>> rng = np.random.default_rng(1)
>>> deltas = 240 * rng.random(180)
>>> ts1 = TimeSeries(time_start="2011-01-01T00:00:00",
... time_delta=deltas*u.minute)
>>> # g band fluxes
>>> g_flux = [0] * 180 * u.mJy
>>> g_err = [0] * 180 * u.mJy
>>> y_g = np.round(3 + 2 * np.sin(10 * np.pi * ts1['time'].mjd[0:60]),3)
>>> dy_g = np.round(0.01 * (0.5 + rng.random(60)), 3) # uncertainties
>>> g_flux.value[0:60] = y_g
>>> g_err.value[0:60] = dy_g
>>> ts1["g_flux"] = MaskedColumn(g_flux, mask=[False]*60+[True]*120)
>>> ts1["g_err"] = MaskedColumn(g_err, mask=[False]*60+[True]*120)
>>> # r band fluxes
>>> r_flux = [0] * 180 * u.mJy
>>> r_err = [0] * 180 * u.mJy
>>> y_r = np.round(3 + 2 * np.sin(10 * np.pi * ts1['time'].mjd[60:120]),3)
>>> dy_r = np.round(0.01 * (0.5 + rng.random(60)), 3) # uncertainties
>>> r_flux.value[60:120] = y_r
>>> r_err.value[60:120] = dy_r
>>> ts1['r_flux'] = MaskedColumn(r_flux, mask=[True]*60+[False]*60+[True]*60)
>>> ts1['r_err'] = MaskedColumn(r_err, mask=[True]*60+[False]*60+[True]*60)
>>> # i band fluxes
>>> i_flux = [0] * 180 * u.mJy
>>> i_err = [0] * 180 * u.mJy
>>> y_i = np.round(3 + 2 * np.sin(10 * np.pi * ts1['time'].mjd[120:]),3)
>>> dy_i = np.round(0.01 * (0.5 + rng.random(60)), 3) # uncertainties
>>> i_flux.value[120:] = y_i
>>> i_err.value[120:] = dy_i
>>> ts1["i_flux"] = MaskedColumn(i_flux, mask=[True]*120+[False]*60)
>>> ts1["i_err"] = MaskedColumn(i_err, mask=[True]*120+[False]*60)
>>> ts1
<TimeSeries length=180>
time g_flux g_err r_flux r_err i_flux i_err
mJy mJy mJy mJy mJy mJy
Time float64 float64 float64 float64 float64 float64
----------------------- ------- ------- ------- ------- ------- -------
2011-01-01T00:00:00.000 3.0 0.012 ——— ——— ——— ———
2011-01-01T02:02:50.231 3.891 0.009 ——— ——— ——— ———
2011-01-01T05:50:56.909 4.961 0.007 ——— ——— ——— ———
2011-01-01T06:25:32.807 4.697 0.014 ——— ——— ——— ———
2011-01-01T10:13:13.359 4.451 0.005 ——— ——— ——— ———
2011-01-01T11:28:03.732 4.283 0.008 ——— ——— ——— ———
2011-01-01T13:09:39.633 1.003 0.015 ——— ——— ——— ———
2011-01-01T16:28:18.550 3.833 0.008 ——— ——— ——— ———
2011-01-01T18:06:31.018 1.02 0.013 ——— ——— ——— ———
... ... ... ... ... ... ...
2011-01-15T16:03:17.207 ——— ——— ——— ——— 4.656 0.014
2011-01-15T17:29:38.139 ——— ——— ——— ——— 1.423 0.01
2011-01-15T20:03:35.935 ——— ——— ——— ——— 4.805 0.008
2011-01-15T21:35:02.069 ——— ——— ——— ——— 3.042 0.007
2011-01-15T23:06:35.567 ——— ——— ——— ——— 1.162 0.01
2011-01-16T01:07:30.330 ——— ——— ——— ——— 4.99 0.009
2011-01-16T01:11:31.138 ——— ——— ——— ——— 5.0 0.011
2011-01-16T03:09:58.569 ——— ——— ——— ——— 1.314 0.01
2011-01-16T07:03:09.586 ——— ——— ——— ——— 3.383 0.005
我们的时间序列数据被设置为异步的，其中给定的时间戳对应于单个频带中的测量。但是，如果您的数据对于多个波段测量或混合测量具有一个时间戳， LombScargleMultiband 仍然可以对其进行手术。
对该示例进行操作 TimeSeries ， LombScargleMultiband 具有加载器功能，如下所示：
>>> ls = LombScargleMultiband.from_timeseries(ts1, signal_column=['g_flux', 'r_flux', 'i_flux'],
... uncertainty_column=['g_err', 'r_err', 'i_err'],
... band_labels=['g', 'r', 'i'])
signal_column 需要与每个波段中的通量或震级测量对应的列的列表。 uncertainty_column 和 band_labels 是可选的，但如果指定，则必须是大小与 signal_column 。 uncertainty_column 指定包含每个波段的关联错误的列，而 band_labels 提供要用于每个光度学标注栏的标签。从这里开始， LombScargleMultiband 可以照常使用。例如：
>>> frequency,power = ls.autopower()
与以下各项保持一致 LombScargle#
LombScargleMultiband 是继承的类 LombScargle ，并被开发为提供类似于的接口 LombScargle 尽可能的。从这个角度来看，有几个核心方面 LombScargle 这一点对 LombScargleMultiband 。
测量不确定度#
这个 LombScargleMultiband 接口还可以处理具有测量不确定性的数据。如上例所示。
周期图和单位#
这个 LombScargleMultiband 接口正确处理 Quantity 附加了单位的对象，并将验证输入以确保单位正确。
指定频率栅格#
如上所示， autopower() 方法自动确定频率网格，使用 autofrequency() 。可调参数与显示的参数相同 LombScargle 。同样地，自定义频率网格可以被直接提供给 power() 功能。
例子#
>>> frequency = np.linspace(0, 2, 1000)
>>> power = LombScargleMultiband(t, y, bands, dy).power(frequency)
周期图实现#
多波段Lomb-Scarger周期图的两个实现可在 LombScargleMultiband ， flexible 和 fast ，它们可通过 power() 方法的 method 参数。 flexible 是Gatspy中使用的LombScargleMultiband算法的直接端口 gatspy 包裹。它构建了一个公共模型，以及每个单独波段的偏移量模型。然后，它将正则化应用于结果模型以限制复杂性，从而为任何给定的多波段时间序列数据集产生灵活的模型。顾名思义， fast 是潜在的更快的替代方案，独立适合每个波段，并按权重组合。独立的逐个频段适合杠杆作用 LombScargle 。因此， sb_method 参数在以下位置提供 power() 选择中使用的单波段方法的步骤 power() 每个乐队都有。请记住， fast 取决于选择的潜在速度 sb_method 。
例子#
flexible ：
>>> frequency, power = LombScargleMultiband(t,y,bands,dy).autopower(method='flexible')
fast 同 fast 也被选为 power() 方法：
>>> frequency, power = LombScargleMultiband(t,y,bands,dy).autopower(method='fast', sb_method='fast')
多频带Lomb-Scarger模型#
这个 model() 方法在选定的频率下将正弦模型与数据进行拟合。正弦模型的复杂性可以通过 nterms_base 和 nterms_band 参数。这些控制可用于基本模型的正弦项的数量(对所有频带通用)和可用于每个频带偏移模型的正弦项的数量。
备注
以下任一项 nterms_base 和 nterms_band 可以设置为0，但不能同时设置为0。在以下情况下 nterms_base =0和 nterms_band =1是一种特殊情况，称为 multi-phase model ，其中基本模型被简化为简单的偏移量，因此波段被独立求解(单波段拟合)。中进一步讨论了 Periodograms for Multiband Astronomical Time Series [6]
例子#
下面的示例使用与上面相同的数据。 autopower() 用于返回周期图，我们可以为我们的模型选择功率最大的频率：
>>> model = LombScargleMultiband(t, y, bands, dy, nterms_base=1, nterms_band=1)
>>> frequency, power = model.autopower(method='flexible')
>>> freq_maxpower = frequency[np.argmax(power)]
然后，我们可以基于找到的频率和时间(按频率相控)进行建模：
>>> t_phase = np.linspace(0, 1/freq_maxpower, 100)
>>> y_fit = model.model(t_phase, freq_maxpower)
然后，所得到的拟合是形状(频带数、时间步数)，或者在这种特定情况下为(5,100)。通过绘制结果图，我们看到模型已经恢复了以正确频率恢复的预期正弦：
虚警概率#
不像 LombScargle ， LombScargleMultiband 没有实施错误警报概率。可用的算法可用于 LombScargle 仅对单项周期图有效，这对于多波段情况下的模型很少有效。