---
title: "复算点透镜放大因子，再检验一次透镜搜索"
description: "用独立的菲涅耳积分核对放大因子的闭合解，再把有、无透镜的啁啾注入合成噪声。"
canonical_url: "https://sandboxphysics.com/zh/stories/gravitational-wave-lensing/recompute-the-point-lens-amplification/"
language: "zh-CN"
page_type: "story"
alternate_languages:
  en: "https://sandboxphysics.com/stories/gravitational-wave-lensing/recompute-the-point-lens-amplification/"
last_modified: "2026-11-12"
---

[引力波透镜与频率条纹](https://sandboxphysics.com/zh/stories/gravitational-wave-lensing/)

# 复算点透镜放大因子，再检验一次透镜搜索

Sandbox Physics · 2026-11-12

这份指南用两条互相独立的路线计算点透镜放大因子，得到的复数相差不到 $4.2\times10^{-4}$。接着检验低频和高频两个极限，再在合成噪声上跑一次小型透镜搜索。结果是：$y=0.5$ 的透镜一眼可辨，$y=2$ 的透镜真假难分，纯噪声大约每九次就有一次冒充出透镜。下面每个数字都来自同一个脚本，跑一遍大约一分钟。

先读[点质量透镜放大因子的精确解](https://sandboxphysics.com/zh/stories/gravitational-wave-lensing/point-lens-amplification-factor/)和[匹配滤波那一章](https://sandboxphysics.com/zh/stories/gravitational-wave-lensing/matched-filtering-a-lensed-chirp/)。这里不用任何探测器数据，信号和噪声全是合成的。

## 下载文件

把下面八个文件放进同一个空文件夹：

- [README.md](https://sandboxphysics.com/zh/stories/gravitational-wave-lensing/materials/README.md)：运行说明和失败案例
- [point\_lens\_reference.py](https://sandboxphysics.com/zh/stories/gravitational-wave-lensing/materials/point_lens_reference.py)：闭式解和两像极限
- [independent\_fresnel\_integral.py](https://sandboxphysics.com/zh/stories/gravitational-wave-lensing/materials/independent_fresnel_integral.py)：计算 $F$ 的第二条路线
- [signal\_injection.py](https://sandboxphysics.com/zh/stories/gravitational-wave-lensing/materials/signal_injection.py)：噪声、啁啾、注入和搜索
- [reproduce.py](https://sandboxphysics.com/zh/stories/gravitational-wave-lensing/materials/reproduce.py)：全部重跑并比对
- [injection-inputs.json](https://sandboxphysics.com/zh/stories/gravitational-wave-lensing/materials/injection-inputs.json)：每次注入的设置和随机种子
- [expected-outputs.json](https://sandboxphysics.com/zh/stories/gravitational-wave-lensing/materials/expected-outputs.json)：本页引用的全部数字
- [provenance.json](https://sandboxphysics.com/zh/stories/gravitational-wave-lensing/materials/provenance.json)：来源、常数和 SHA-256 校验和

然后运行：

```
python3 -m pip install numpy scipy mpmath
python3 reproduce.py
```

脚本会打印与保存结果的最大差异。所有放大因子相差不超过 $2\times10^{-6}$、所有注入中位数相差不超过 $10^{-6}$ 时，退出状态为 0。我们的环境是 Python 3.9.6、numpy 2.0.2、scipy 1.13.1 和 mpmath 1.3.0。换了版本，随机噪声的末几位可能变化，严格比对就会失败；这时用 python3 reproduce.py --write 保存你自己的结果，中位数应当和我们的接近。

## 单位与约定

**表 5.1** 所有脚本共用的单位与约定。

| 量 | 约定 |
| --- | --- |
| 透镜质量 | 红移后的质量 $M_{Lz}$，以太阳质量为单位，$GM_\odot/c^3=4.925490947\times10^{-6}$ s |
| 频率 | 无量纲的 $w=8\pi GM_{Lz}f/c^3$ |
| 偏移 | $y$，以爱因斯坦半径为单位 |
| 到达时间 | $T=\tfrac12\lvert x-y\rvert^2-\ln\lvert x\rvert-\phi_m(y)$，第一个像在 $T=0$ |
| 傅里叶变换 | $\tilde h(f)=\int h(t)\,e^{+2\pi i f t}\,\mathrm dt$，延迟对应乘上 $e^{+2\pi i f\Delta t}$ |
| 内积 | 从 20 Hz 到 219.9 Hz 求和 $4\,\mathrm{Re}\sum\tilde a^*\tilde b/S_n\,\Delta f$ |

太阳质量参数取 IAU 2015 名义值 $1.3271244\times10^{20}\,\mathrm{m^3\,s^{-2}}$，再除以 $c^3$。

## 第一步：闭式解

point\_lens\_reference.py 用 mpmath 以 60 位有效数字计算 Takahashi 和 Nakamura（2003）的式 (17)，以及他们的两像公式式 (18)。式 (17) 完全按印出的形式使用。[Takahashi 和 Nakamura 2003，第 2 节](https://arxiv.org/abs/astro-ph/0305055v2)

## 第二步：完全不用闭式解的第二条路线

independent\_fresnel\_integral.py 从透镜平面积分重新出发，把透镜平面上的点按到达时间分组。设 $A(\tau)$ 是 $T<\tau$ 的区域面积，则

$$
F(w)=\frac{w}{2\pi i}\int_0^\infty\frac{\mathrm dA}{\mathrm d\tau}\,e^{iw\tau}\,\mathrm d\tau=1-iw\int_0^\infty\bigl[I(\tau)-1\bigr]e^{iw\tau}\,\mathrm d\tau,\qquad I(\tau)=\frac1{2\pi}\frac{\mathrm dA}{\mathrm d\tau}.
$$

第二个等号用到 $\int_0^\infty e^{iw\tau}\mathrm d\tau=i/w$。$I(\tau)$ 是透镜对一个尖脉冲的响应，没有透镜时等于 1，也就是 2003 年论文式 (16) 的时域图像。

脚本分四步算 $I(\tau)$：

1. 对从透镜出发的每个方向，在 $\ln r$ 上用二分法找到 $T=\tau$ 的位置，射线上极小值两侧各找一个。
2. 用 240 点高斯–勒让德求积对方向积分，先做变量替换，消掉等值线闭合处的平方根奇点。
3. 到达时间超过鞍点后，把求积节点集中到透镜另一侧，那里等值线收窄。
4. 对列表中的 $I(\tau)$ 按分段线性函数做精确傅里叶积分，$\tau=400$ 以外的尾部用 $c_1/\tau+c_2/\tau^2$ 拟合补上。

这条路线完全不碰伽马函数和超几何函数，所以实部和虚部都对上，检验的就不只是式 (17) 的模，还有它的相位约定。

**结果。** 在 $y\in\{0.1,0.3,1,3\}$、$w\in\{0.1,1,3,10,30\}$ 上，最大相对差是 $4.2\times10^{-4}$，出现在 $y=0.3$、$w=30$。每个 $y$ 拟合出的尾部系数都是 $c_1=0.499$。

**表 5.2** 20 个测试点中的 4 个：闭式解与独立菲涅耳路线的对比。

| y | w | 菲涅耳路线 | 闭式解 |
| --- | --- | --- | --- |
| 0.1 | 1 | 1.58197 − 0.87469i | 1.58191 − 0.87445i |
| 0.3 | 10 | 1.15235 − 0.99952i | 1.15249 − 0.99890i |
| 1 | 3 | 1.05026 − 0.38996i | 1.05027 − 0.38982i |
| 3 | 1 | 1.09237 + 0.01542i | 1.09236 + 0.01542i |

全部 20 行都在 expected-outputs.json 里。

到达时间里的常数在这里很要紧。如果照 2003 年论文印出的式 (12) 在 $T$ 里写 $+\phi_m(y)$，菲涅耳路线得到的 $F$ 会多出相位 $e^{2iw\phi_m}$，就和式 (17) 对不上了。论文正文说最早到达时间取为零，这要求写 $-\phi_m$。我们按这个做，除此之外不对论文下任何结论。

## 第三步：两个极限

这些检验用的是闭式解，它已经在第二步得到验证。

- **低频。** $\lvert F-1\rvert$ 应当几乎与 $w$ 成正比地减小，因为有 $w\ln w$ 项，减得稍慢一点。$y=1$ 时，$w=10^{-2}$、$10^{-3}$、$10^{-4}$ 对应 $2.2\times10^{-2}$、$3.3\times10^{-3}$、$4.4\times10^{-4}$。$y=0.1$、$w=10^{-4}$ 时是 $5.13\times10^{-4}$，一阶公式 $\pi w/4+(iw/2)[\ln(w/2)+\gamma-2\phi_m]$ 给出同样的值。
- **高频。** 与两像公式的差除以 $\lvert F\rvert$，应当随 $w$ 增大而减小。$y=1$ 时从 $w=30$ 的 $2.6\times10^{-3}$ 降到 $w=1000$ 的 $2.6\times10^{-5}$。$y=0.3$ 时，$w=30$、100、300、1000 依次是 $3.1\times10^{-2}$、$6.5\times10^{-3}$、$8.7\times10^{-3}$、$1.1\times10^{-3}$。100 到 300 之间的回升是真实的：剩余的差本身随 $w$ 振荡。

## 第四步：注入与搜索

signal\_injection.py 生成一段 8 s、每秒 1024 个采样的数据：

1. **噪声。** 按 Ajith（2011）式 (4.7) 着色的高斯噪声，每个频率格的实部和虚部标准差都是 $\sqrt{S_n/(4\Delta f)}$。[Ajith 2011，第 IV 节](https://arxiv.org/abs/1107.1267v4)
2. **信号。** $10+10\,M_\odot$ 的最低阶啁啾，从 20 Hz 到 219.9 Hz，缩放到无透镜信噪比 20；透镜注入再乘上精确的 $F$。[Cutler 和 Flanagan 1994，式 41–42](https://arxiv.org/abs/gr-qc/9402014v1)
3. **无透镜搜索。** 黑洞质量从 9.6 到 $10.4\,M_\odot$ 取 17 个值，每个都对幅度、到达时间和相位取最大。
4. **透镜搜索。** 透镜质量在 10 到 $1500\,M_\odot$ 之间取 36 个值，$y$ 在 0.1 到 3 之间取 9 个值；然后对每个 $y$，在最佳质量的两个相邻格点之间再细扫 48 个质量；最后对最好的两组做有界的一维优化。透镜模板保持真实的黑洞质量。
5. **统计量。** $\Delta\ln L$，即最佳透镜模板与最佳无透镜模板的对数似然比之差，报告门槛取 4。

加噪声之前，脚本先检查搜索能否找回注入的透镜。两种透镜情形都找回了：透镜质量误差在 0.001% 以内，$\Delta\ln L$ 和精确模板的值相差不到 $10^{-8}$。

**结果。**

**表 5.3** 无透镜信噪比为 20 时的注入结果，单位为 $\Delta\ln L$。

| 情形 | 试验次数 | 无噪声 | 中位数 | 第 90 百分位 | 最大 | 超过 4 |
| --- | --- | --- | --- | --- | --- | --- |
| 透镜，$M_{Lz}=1000$，$y=0.5$ | 20 | 115.2 | 110.2 | 153.3 | 165.9 | 100% |
| 透镜，$M_{Lz}=1000$，$y=2$ | 40 | 6.03 | 4.27 | 8.27 | 12.49 | 52.5% |
| 没有透镜 | 120 | 0 | 2.49 | 4.35 | 7.97 | 11.7% |

清楚的透镜不会认错。弱透镜像抛硬币，而纯噪声的尾巴能伸到弱透镜也会给出的数值。

## 有意保留的失败案例

- **格点太粗。** 第一版只在 36 个质量格点上搜索，相邻格点相差约 15%。$y=2$ 时延迟是 90 ms，相邻质量之间 200 Hz 处的条纹相位会差好几个弧度，所以即使没有噪声，最好的格点模板也只有 $\Delta\ln L=0.73$，而不是 6.03。那一版弱透镜的中位数 2.79 基本是噪声。第四步的细化修正了这个问题，无噪声检查也能再次发现它。
- **放大因子表太粗。** 第一版 $F$ 表在 $w$ 上只有 1500 个点。$y=3$、$w$ 约为 41 时，每个条纹只采到约四个点，插值误差达到 3.8%。改成 6000 个点后降到 0.24%。
- **假透镜。** 120 次纯噪声试验中有 14 次越过门槛。最大的一次是种子 3058，数值 7.97，偏好 $M_{Lz}\approx1150$、$y=1.5$。
- **没有包含的东西。** 没有自旋、不等质量、并合、铃宕、探测器天线方向图、突发干扰，也没有透镜周围的星系；$y$ 只在格点上搜索。这些数字用来说明方法，不是任何真实搜索的灵敏度。

## 相信一次重跑之前，逐项核对

- python3 reproduce.py 退出状态为 0，或者你自己用 --write 得到的结果只在与噪声有关的末位上和我们不同。
- provenance.json 里的 SHA-256 校验和与你的文件一致。
- 菲涅耳路线和式 (17) 在实部和虚部上都一致，而不只是模一致。
- 低频时 $\lvert F-1\rvert$ 大致按 $w\ln w$ 减小，高频时两像误差总体在减小。
- 读任何带噪声的结果之前，无噪声搜索先找回了注入的透镜。
- 每个透镜结果旁边都报告了纯噪声的分布，并且用同一个门槛。

## 到实验里试试

[引力波透镜的频率条纹](https://sandboxphysics.com/zh/labs/gravitational-wave-lens-fringes/)：把实验调到注入的透镜参数，对照匹配度。

## 来源

1. [Takahashi 和 Nakamura 2003，第 2 节](https://arxiv.org/abs/astro-ph/0305055v2)
2. [Ajith 2011，第 IV 节](https://arxiv.org/abs/1107.1267v4)
3. [Cutler 和 Flanagan 1994，式 41–42](https://arxiv.org/abs/gr-qc/9402014v1)

## 继续阅读

引力波透镜与频率条纹 · 阅读路径

### 从两条路径到透镜搜索

1. [文章 **同一段啁啾走了两条路，探测器听到什么**](https://sandboxphysics.com/zh/stories/gravitational-wave-lensing/one-chirp-two-paths/)
2. [文章 **透镜为什么让引力波晚到**](https://sandboxphysics.com/zh/stories/gravitational-wave-lensing/diffraction-and-the-lens-time-delay/)
3. [文章 **点质量透镜放大因子的精确解**](https://sandboxphysics.com/zh/stories/gravitational-wave-lensing/point-lens-amplification-factor/)
4. [文章 **探测器分得清透镜化的啁啾吗**](https://sandboxphysics.com/zh/stories/gravitational-wave-lensing/matched-filtering-a-lensed-chirp/)
5. [当前页面 **复算点透镜放大因子，再检验一次透镜搜索**](https://sandboxphysics.com/zh/stories/gravitational-wave-lensing/recompute-the-point-lens-amplification/)
6. [文章 **怎样读第三轮观测的引力波透镜搜索**](https://sandboxphysics.com/zh/stories/gravitational-wave-lensing/reading-the-o3-microlensing-search/)