--- title: "拟合周期，再测一次假峰风险" description: "推导加权浮动均值搜索，校准整个频率网格的最高峰，再把信号探测与选对混叠分开。" canonical_url: "https://sandboxphysics.com/zh/stories/observing-cadence/period-search/" language: "zh-CN" page_type: "story" alternate_languages: en: "https://sandboxphysics.com/stories/observing-cadence/period-search/" last_modified: "2026-10-28" --- [决定什么时候拍下变化的天空](https://sandboxphysics.com/zh/stories/observing-cadence/) AC07 · L3 # 拟合周期，再测一次假峰风险 Sandbox Physics · 2026-10-28 周期搜索在每个频率上拟合一个模型，峰告诉我们哪个模型最配这批样本。要判断假探测风险，还得知道同一套搜索遇到噪声时，最高峰会爬多高。 ## 开始前需要什么 完成“数清两次曝光之间转了几圈”。加权平均让精确测量占更多分量；最小二乘选择残差平方和最小的系数。下面用消元解二元方程，不要求先学矩阵。 ## 每个频率都重新拟合常数 把每个观测写成常数、余弦、正弦与误差之和。固定频率后，三个未知系数都线性出现；余弦和正弦一起允许任意初相位。误差独立且为高斯时，给定标准误差的逆方差给出权重。 如果拍摄漏掉了一部分相位，样本均值就可能偏离源的真实均值。只减一次样本均值，再强迫波通过零，可能把真信号压下去。让常数随每个试探频率重新拟合，就得到广义的、或叫浮动均值的 Lomb–Scargle 拟合。 时间用天，频率用每天周数。浏览器拟合归一化线性强度，真实教程拟合 I 波段星等。两者都能作数值正弦拟合，振幅却表示不同物理量。 $$ y_j=c+a\cos(2\pi ft_j)+b\sin(2\pi ft_j)+\epsilon_j,\qquad \chi^2=\sum_j\frac{[y_j-c-a\cos(2\pi ft_j)-b\sin(2\pi ft_j)]^2}{\sigma_j^2} $$ 每个 f 都选择让无量纲误差量最小的 a、b、c。误差假设决定它的统计解释。 对应阅读：Zechmeister & Kürster, generalized Lomb–Scargle, arXiv v1 · VanderPlas, Understanding the Lomb–Scargle Periodogram, arXiv v1 ## 先消常数，再解两个方程 导数描述微调一个系数时，误差量会怎样变化。光滑最小值处，一阶变化为零。残差平方求导，得到两倍残差，再乘被微调项前面的负系数；让各项加权和为零，就得到拟合方程。 把逆方差权重归一，让总和为一。用 Y、C、S 表示数据、余弦、正弦的加权均值。对常数求导并令其为零，得到常数等于数据均值，减掉两个拟合三角项的均值。 代回原式，就把两条三角函数列都中心化。CC、SS 是加权离差平方，CS 是交叉乘积，YC、YS 是中心化数据与这两列的交叉乘积。对余下两个系数求导，就得到一对二元一次方程。 第一式乘 SS，第二式乘 CS，再相减，就能消去 b、解出 a。交换消元顺序便得到 b。若行列式为零，表示这些采样在该频率上没提供独立的两个分量，程序会跳过这个试探点。 $$ w_j=\frac{\sigma_j^{-2}}{\sum_k\sigma_k^{-2}},\quad Y=\sum_jw_jy_j,\quad C=\sum_jw_j\cos(2\pi ft_j),\quad S=\sum_jw_j\sin(2\pi ft_j),\quad c=Y-aC-bS $$ 加权中心化直接处理不完整的相位分布，不假定三角函数在样本中平均为零。 $$ CC\,a+CS\,b=YC,\qquad CS\,a+SS\,b=YS,\qquad D=CC\,SS-CS^2 $$ 例如 CC 是余弦减去 C 后的加权平方和，YC 是数据减去 Y 与余弦减去 C 的加权乘积和，其他量按同样方式定义。 $$ a=\frac{SS\,YC-CS\,YS}{D},\qquad b=\frac{CC\,YS-CS\,YC}{D} $$ 这是二元消元结果。保存的 Python 重跑另外采用加权正交投影计算。 对应阅读：Zechmeister & Kürster, generalized Lomb–Scargle, arXiv v1 · VanderPlas, Understanding the Lomb–Scargle Periodogram, arXiv v1 ## 把拟合改善画成周期图 先用相同权重拟合一个恒定参考。在每个频率，把最佳正弦、余弦、常数模型的残差平方和，与恒定参考比较，再把改善除以参考值。标准功率就在零到一之间：零表示没改善，一表示样本被完美配出。 功率 0.9，不代表周期有九成概率正确。点很少时，灵活的正弦也能在噪声上取得很大改善；频率搜得更多，碰巧遇到高峰的机会也更多。 有限时间跨度会形成宽度大致为跨度倒数的频率结构。网格应在这段宽度里放几个点，真实教程取五分之一，并保存完整网格。上下界是按问题声明的选择，不能由平均间隔 Nyquist 规则包办。 $$ P(f)=1-\frac{\chi^2(f)}{\chi^2_0},\qquad \Delta f_{\rm grid}=\frac1{5T_{\rm span}} $$ 细网格采出有限基线峰的形状，可用信息仍由基线和混叠图样决定。 对应阅读：Astropy 8.0.2 Lomb–Scargle documentation · VanderPlas, Understanding the Lomb–Scargle Periodogram, arXiv v1 ## 校准整个搜索的最大值 重跑的零假设，是恒定源加上独立高斯误差，大小采用原来给出的误差。每个真实时刻抽一个噪声值，把全部声明频率搜一遍，保存最高功率，再重复很多次。这样得到的最大值分布，同时包含观测窗口和搜索范围的作用。 数一数有多少零假设最大值，至少与实测最大值一样高。B 次试验、K 次超过时，下面的有限蒙特卡洛秩给出非零结果，即使一次也没超过。数值分辨约为试验次数的倒数；短模拟不能保证任意稀有的尾概率。 用另一批独立噪声测试阈值。实测教程先用 4096 次零假设搜索校准百分之五阈值，再用 1024 次新搜索验证。完整与稀疏时刻的实际验证误警报率为 4.39%、4.49%。检查针对声明的独立高斯模型；红噪声或共享标定误差，需要另换零假设。 $$ \mathrm{FAP}(z)=\Pr\!\left[\max_{f\in\mathcal G}P(f)\ge z\mid\text{constant source and stated noise}\right],\qquad p_{\rm MC}=\frac{K+1}{B+1} $$ 条件放在噪声模型上，既不是给定数据后源恒定的概率，也不是所选频率正确的概率。 对应阅读：Astropy 8.0.2 Lomb–Scargle documentation · VanderPlas, Understanding the Lomb–Scargle Periodogram, arXiv v1 ## 例题：更高的峰，反而更弱的证据 固定 OGLE 首季有二十五个测量，最高功率 0.6418，高斯零假设尾秩为 0.004882。每三点只保留一点，剩九点后，功率升至 0.8992，尾秩却升至 0.2060。 稀疏模型配出了更多所选样本，九个噪声样本却也很容易在宽搜索里配出漂亮结果。只比峰高、不按采样重新校准，就会把意思读反。两套真实时刻分别模拟零假设，并使用同一个明确的搜索网格。 下方合成浏览器例子作另一种受控比较。已知源周期，可以直接给每夜一次的混叠贴上正确标签。 合成曝光，种子 707；源周期固定为 0.635 天，由 freeze-synthetic.ts 生成。 | 观测计划 | 曝光数 | 最强周期（天） | 闪光触发数 | | --- | --- | --- | --- | | 每夜一次 | $24$ | $2.352941$ | $0$ | | 每夜三次 | $72$ | $0.634921$ | $0$ | | 抖动与缺口 | $19$ | $0.634921$ | $0$ | | 漏掉闪光 | $96$ | $0.634921$ | $0$ | | 拍到闪光 | $96$ | $0.634921$ | $4$ | 对应阅读：OGLE-IV LMC RR Lyrae official data README · Zechmeister & Kürster, generalized Lomb–Scargle, arXiv v1 · Astropy 8.0.2 Lomb–Scargle documentation ## 选错混叠与没探测到，分开数 误警报率问噪声多常被认作探测，没告诉我们真信号多常漏掉、或周期认错。要回答后两个问题，需要在相同时刻注入指定信号与新噪声，运行不变的搜索及阈值，再与已知输入比较。 真实教程做 128 次合成星等正弦注入，半振幅 0.1 星等，相位随机，频率分布在每天 0.3 到 4.8 周。恢复频率距离输入小于一个基线倒数，算作正确。二十五点计划中，112 次越过阈值，其中 36 次选中正确频率附近，76 次选成其他混叠；十六次没过阈值。 这些计数对应指定信号族、振幅、频率范围与误差模型。换一类源，要另做注入实验。即使还没定好源种群，这项区分仍然有用：探测阈值与周期识别规则回答不同问题。 $$ \text{correct detected}+\text{aliased detected}+\text{missed}=\text{all injections} $$ 先定义分类，再测恢复率。假阳性来自另外的零假设样本。 对应阅读：VanderPlas, Understanding the Lomb–Scargle Periodogram, arXiv v1 · Astropy 8.0.2 Lomb–Scargle documentation ## 练习与解释性答案 1. 同学在小频率范围校准了百分之五阈值，随后扩大搜索，能沿用吗？ 展开解释 更宽搜索给噪声最大值更多机会。需要按扩大的网格、相同观测时刻与噪声模型重新校准。阈值属于整套搜索程序。 2. 一百次注入，九十次过阈值，其中三十五次频率正确、五十五次混叠；十次漏掉。探测率与正确周期恢复率各多少？ 展开解释 探测率九成，正确且探测到的恢复率只有三十五次除以一百。若只问探测样本内部，则三十五次除以九十选中了声明的正确频率邻域。 3. 相邻时刻的残差总倾向同号，信任独立噪声 FAP 之前，应改什么？ 展开解释 先研究时间相关噪声或趋势。独立抽样的高斯模拟没有保留这类相关性，需要有物理理由的协方差或经过验证的相关噪声零假设，再校准整个搜索。 ## 动手检验预测 ### [诊断一个很有说服力的错误峰](https://sandboxphysics.com/zh/labs/observing-cadence-aliases/) 1. 每夜一次时，比较最强周期与已知人工源周期，先找出其他混叠。 2. 再次选择同一预设，导出两份样本；时刻保持不变，探测种子改变，比较峰高变化与混叠族是否持续存在。 3. 用 OGLE 重跑分别比较全频搜索误警报和正确周期注入回收，不从峰高独自推断任何一个。 [上一篇：数清两次曝光之间转了几圈](https://sandboxphysics.com/zh/stories/observing-cadence/sampling-windows/) [下一篇：用真实 OGLE 变星数一次周期混叠](https://sandboxphysics.com/zh/stories/observing-cadence/ogle-lightcurve/) ## 来源与计算范围 1. [VanderPlas, Understanding the Lomb–Scargle Periodogram, arXiv v1](https://arxiv.org/pdf/1703.09824v1) 2017 年 3 月 28 日第一版；3.2、4.2 节与图 7–14 讲窗口及混叠，6.1–6.2 节讲加权浮动均值，7.1、7.4.2 节与图 26 讲网格和误警报。期刊版本为 ApJS 236,16（2018），课程固定阅读此预印本版本。 2. [Zechmeister & Kürster, generalized Lomb–Scargle, arXiv v1](https://arxiv.org/html/0901.2573v1) 2009 年 1 月 16 日第一版，第 2 节式 4–18 与附录 A.1：加权正弦、余弦、常数项的正规方程与归一化卡方改善。重跑使用正交投影这一独立数值写法。 3. [Astropy 8.0.2 Lomb–Scargle documentation](https://docs.astropy.org/en/stable/timeseries/lombscargle.html) 2026 年 10 月 28 日核读，页面更新于 10 月 15 日。阅读附加参数、标准归一化、峰显著性与 RR Lyrae 示例。重点是浮动均值、按相同时刻与网格校准 FAP，以及不规则采样没有通用平均 Nyquist 截止。 4. [OGLE-IV LMC RR Lyrae official data README](https://www.astrouw.edu.pl/ogle/ogle4/OCVS/lmc/rrlyr/README) 定义 HJD minus 2450000、I 波段星等和误差，并说明公开数据引用要求。单对象文件 OGLE-LMC-RRLYR-00001.dat 与目录行均固定保存，附 SHA-256 和 HTTP 元数据。