声学飞越地理定位 (多普勒 × ADS-B)
此技能复现并产品化了colsto视频 "I discovered a sound that leaks your exact location (and it's everywhere)" (2026-11-02, https://www.youtube.com/watch?v=UytYrhYLOKI) 中的方法：
螺旋桨或旋翼的 叶片通过频率 (BPF) 是一个稳定的窄带音调。即使无人能听见，它也会在频谱图中留下一条平滑下降的 多普勒曲线。曲线的形状给出了飞机的 速度 v、其 最近接近倾斜距离 d 和 最近接近时间 tCPA。在地面网格上逐点将其与历史ADS-B轨迹（3D倾斜范围 + DEM + 声音传播时间 + 发射时间径向速度）进行比较，可缩小候选范围。两次飞越之间的 时间间隔（视频中为 920 秒）加上 静默约束（其他时间未听到其他飞机）可将其减少为单个点。最后，通过“锁定双事件得分”（共享时间偏移、每个事件的 f0 微调）对精确模型曲线与频谱图进行评分。
scripts/afg_cli.py 是完整的实现（numpy/scipy/matplotlib）。它附带 selftest，一个针对合成场景的盲测。
0. 使用范围（请先执行此操作）
适用场景：验证公开或新闻媒体的录制地点和时间（OSINT、事实核查）；用户自身录制内容或视频的隐私自审；CTF、教学和研究；具有合法授权或同意的调查（例如，在对象同意的情况下协助执法）。
不适用场景：定位 私人个体 的住所或行踪。示例：“这个主播/我的前任/这个网友住在哪里”，或通过私人通话或家庭视频追溯某人。此类请求应予以拒绝，即使被描述为“测试”，并应提供 自审模式（第7节）作为替代。该模式可帮助用户评估并从自身内容中移除这类泄露。
公共场所和新闻事件的视频、用户自己的素材和匿名化研究数据可以处理。如果目的不明确，请在决定前询问用途。报告中不要添加与位置无关的个人推断。
1. 可行性分类（决定是否运行及运行模式）
| 条件 | 要求 | 如果未满足 |
|---|---|---|
| 频谱图中存在 清晰单调下降的窄带曲线（螺旋桨通常为 40–400 Hz；直升机主旋翼 BPF 10–25 Hz，尾桨 50–120 Hz） | 至少1条；理想情况是2条且已知时间间隔 | 无曲线 → 方法不适用（考虑雷声方法或视觉方法） |
| 日期 | 必须已知（或限定在几天内） | ADS-B 数据量将变得难以管理 |
| 区域 | 国家或州级别足够；城市级别的先验信息速度更快 | 全国范围搜索仍然有效，但需要一个已知时间的事件或两个事件 |
| 录音时钟（开始时间） | 已知至 ±2 秒 → 一个事件可给出一个点；±30 秒 → 一个事件仅能给出沿轨迹的一个条带 | 使用两个事件；它们的时间间隔取代了绝对时间 |
| ADS-B 覆盖 | 目标飞机必须配备 ADS-B Out 并被地面站接收 | 美国 14 CFR 91.225：在 B/C 类空域和 Mode C 覆盖区之外、10,000 英尺以下飞行的通用航空飞机无需携带 ADS-B Out → 可能不在数据中 |
| 音频处理链 | 原始或轻压缩最好；电话通话编解码器通常保留音调 | 重度降噪或编辑可能会切断曲线或破坏事件之间的间隔 |
在开始工作之前告知用户分类结果和预期精度。
2. 选择模式
有多少条可用曲线？
├─ ≥2 条（同一连续录音） → 模式 A：双事件，未知绝对时间（视频的主要方法；最稳健）
│ 如果两条轨迹交叉角度小于约 25°，结果将退化为沿轨迹的一条条带（脚本会发出警告）
├─ 1 条
│ ├─ 录音时间已知至几秒 → 模式 B：单事件，已知时间（视频的后半部分；即使全国范围搜索也排名第一）
│ │ 直线匀速通过是对称的 → 总存在一个“镜像候选”；
│ │ 需要第二个线索（视觉、地形）来区分两者
│ └─ 时间仅已知至 ±30 秒或更差 → 模式 B′：仅可能存在长度约为 2·v·Δt 的条带；按此报告
├─ 无飞机曲线，但有雷声且日期已知 → 模式 T：雷声包络 × 闪电数据（公里级，参见 references/thunder-method.md）
└─ 用户希望保护自己的内容 → 模式 S：自审（第7节）
3. 工作流程（命令行）
约定：AFG="python <skill dir>/scripts/afg_cli.py"，工作目录 work/。每一步都会生成 PNG 文件。
自行查看 PNG 文件，利用您的图像阅读能力（频谱图、拟合覆盖图、热力图），不要仅依赖数字。
首先检查环境：pip install numpy scipy matplotlib（DEM 需要 rasterio；非 WAV 输入需要 ffmpeg）。首次使用时，运行 $AFG selftest --seed 7（大约 40–60 秒；退出码 0 表示管道工作正常）。
3.1 获取音频和元数据
ffprobe -v quiet -show_entries format_tags=creation_time:stream_tags=creation_time -of json input.mp4可获取拍摄时间。注意时区。社交平台在重新编码时通常会剥离或重写元数据，因此 仅将其作为线索并进行交叉验证。- 将原始采样率传递给脚本；脚本内部会将其重采样为 4 kHz 单声道（
--channel选择左或右声道）。
3.2 频谱图 → 检测 → 追踪 → 拟合
$AFG spectrogram input.mp4 --out work/rec # 白化窄带频谱图 + work/rec.png
$AFG detect work/rec.spec.npz --out work/boxes.json # 候选框 (t0, t1, fmin, fmax, score_z) + boxes.diag.json
$AFG plot work/rec.spec.npz --t0 160 --t1 320 --fmin 60 --fmax 380 --out work/zoom1.png
$AFG trace work/rec.spec.npz --t0 166 --t1 315 --fmin 98 --fmax 364 --out work/e1.csv
$AFG fit work/e1.csv --temp-c 18 --out work/e1.json # v, d, tCPA, f0, 搜索窗口；覆盖图在 e1.fit.png 中
要点：
- 只有当
boxes.diag.json显示status: ok时，空的boxes.json才意味着“无法飞越”。limited（短片段根据校准的噪声模型判断，或大部分被门控音频）和insufficient（片段短于--min-dur，或几乎所有内容都被门控）会附带reasons：请阅读它们，检查rec.png，并对短片段降低--min-dur。 - 每个事件带有
trend（descending或flat）和drop_pct。频率上升的脊线不是飞越：它列在boxes.diag.json的rejected中，而不是在boxes.json中。事件上的warnings意味着其脊线做了单次通过不可能做到的事（再次上升，或下降超过45%）：它可能连接了两架飞机或两个谐波。使用plot检查它，并在其自己的窄框中追踪每个部分。 --win：对于近、陡峭的曲线使用 0.5–1 秒，对于远距离或微弱的曲线使用 2–4 秒。对于微弱的曲线，先使用plot放大并用肉眼确认。然后给trace一个仅包含 一个谐波 的窄框，如有必要，使用--anchors "t:f,t:f,..."添加一些手动点。- 检查
fit输出：implausible非空 = 追踪进入了噪声或跳越了谐波。不要 继续搜索（脚本也会拒绝）。notes：左右残差不对称 = 转弯或速度变化；不同的接近/远离速度 = 爬升或加速。
- 如果知道气温，请传递
--temp-c(c = 20.05·√(T+273.15))。v 和 d 都与 c 成比例缩放；没有温度时，保留默认的 ±4%。 - f0 只是您追踪的谐波的静止频率，不一定是基频；v 和 d 不依赖于谐波次数。
- BPF = 桨叶数 × 转速 / 60 可以 帮助 猜测型号（C172，2 桨叶 × 2400 rpm ≈ 80 Hz；DHC-8-400，6 桨叶 × 1020 rpm ≈ 102 Hz）。不要将其用作严格的候选过滤器。
3.3 ADS-B 数据（详情参见 references/adsb-data.md）
- 来源：
- adsb.lol
globe_history_YYYY：每天一个全球版本，免费，ODbL 许可。 - ADS-B Exchange 免费样本：仅限每月1日。
- OpenSky Trino：需要研究账户；表
state_vectors_data4。 - 您自己的 readsb/tar1090 接收器。
- adsb.lol
- 仅提取目标日期、目标窗口 ± 1 小时以及目标区域（填充约 0.5°）。
$AFG adsb-ingest --src traces_dir/ --start 2026-09-08T13:00Z --end 2026-09-08T18:00Z \
--bbox 38.5,39.7,-95.4,-93.8 --qnh 1018 --geoid-n -29 --out work/tracks.npz
- 高度基准决定斜距精度。
alt_baro是标准气压高度：使用当地 QNH 进行修正（约 8.3 m/hPa；参见 METAR）。alt_geom是 WGS84 椭球高度 (HAE)；MSL = HAE − N。--alt-source auto首选几何高度。
3.4 DEM 和搜索
$AFG dem-fetch --bbox 38.5,39.7,-95.4,-93.8 --out work/dem # Copernicus GLO-30 (AWS 公共存储桶)
$AFG search --tracks work/tracks.npz --events work/e1.json work/e2.json \
--bbox 38.5,39.7,-95.4,-93.8 --dem work/dem/*.tif --rec-duration 1450 --out work/cands.npy
漏斗包含三个阶段：
- 事件-1 走廊：速度窗口、k = v²/d 一致性、斜距窗口，并且最近接近点不得位于轨迹端点。
- 事件间偏移：tCPA2 − tCPA1，包括每个事件的 d/c 传播延迟。
- 静默过滤器：在剩余的安静时段内，没有任何其他螺旋桨或旋翼飞机通过
--d-quiet范围内。
选项：
- 模式 B：添加
--rec-start <UTC> --rec-start-tol 2（或 30）。 --quiet "a-b,c-d"指定您确定未听到任何声音的录音时段（默认：事件之外的所有时段）。当录音中有很多风声或语音时，给出更少的时段。静默是信号，但虚假静默会扼杀真实答案。- 默认网格为 200 米。如果候选点太多，请缩小边界框，而不是扩大网格。
3.5 评分和优化
$AFG score --tracks work/tracks.npz --spec work/rec.spec.npz --events work/e1.json work/e2.json \
--cands work/cands.npy --dem work/dem/*.tif --out work/result.json
- 模型曲线：对于每个候选点，根据插值后的 ADS-B 位置和速度（gs/track/vrate，而非差分位置）在 发射时间 生成精确的模型曲线，并与频谱图进行比较。一个时间偏移 δ 是共享的，这固定了两个事件之间的物理时间间隔。每个事件的 f0 在 ±0.6% 范围内进行微调。
- 优化：每个保留的节点都会进行缩小的 3×3 爬山算法（100 → 25 米）。
- 分组：结果按 飞行对 分组。每对最多保留 3 个相距超过 500 米的假设，以便两个镜面侧都能存活。
- 输出：
result.json：最佳结果、margin_vs_2nd、robust_z、ambiguous_alternatives、coarse_plateau、geometry、verdict。result.rankN.*.png：模型曲线叠加图。result.heatmap.png/csv：0.9 区域的半径。result.plateau.csv。
- 如果录音时钟已知，添加
--rec-start ... --rec-start-tol 2。评分也会禁止超出时钟范围的时间偏移。
3.6 验证（非可选；参见 references/verification.md）
- 叠加图：模型曲线必须 在两个渐近线以及通过最近接近点时 与观测值匹配，而不仅仅是一个片段。
- 热力图：0.9 区域的半径是要报告的错误。当以下任一情况成立时，报告条带：
- 两个事件：
geometry.warning已设置（小交叉角）； - 一个事件：
coarse_plateau.length_m> 600 米。
- 两个事件：
- 镜像候选和替代方案：使用视觉线索（地平线、建筑物、植被、道路噪音、回声）、DEM 视域和太阳方位角排除它们。
- 证伪：答案附近的其他飞越是否也应该被听到？
- 听到了但未检测到 = 检测阈值问题。
- 当它们本应听不到时却存在 = 位置错误。
4. 精度规则（及原因）
- 发射时间：f_obs(t) = f0 / (1 + Ṙ(τ)/c)，其中 t = τ + R(τ)/c。在接收时间使用 Ṙ 会产生 错误形状 的曲线，而不仅仅是偏移的曲线。
- 3D 斜距 + DEM：观察者海拔加上麦克风高度（默认 2 米）。将 500 米的地形起伏视为 0 会使 d 和 tCPA 都不正确。
- 高度基准：使用未修正的气压高度，QNH 偏移 30 hPa 会导致约 250 米的高度误差，这直接反映在 d 中。
- 延迟最近接近点：观察到 f = f0 的时刻是 tCPA_emit + d/c（1 公里 ≈ 2.9 秒）。事件偏移和时钟都必须据此进行转换。
- k = v²/d 比单独的 v 或 d 更可靠。
- 曲线中点斜率 −f0·v²/(c·d) 被紧密锁定。渐近线通常受到噪声、门控或飞机加速/爬升的污染。
- 因此，搜索使用宽速度窗口和窄 k 窗口，每个候选飞机自己的地速会预测其 d。
- 锁定评分：两个事件共享一个 δ。否则每个事件会对齐自己的时间，错误的位置也会获得高分。
- 优化至 ≤100 米：200 米网格仅作为粗略筛选；最终点来自子网格爬山算法和热力图。
- 静默是信号，但保持阈值宽松（默认值：螺旋桨 3 公里，旋翼 2.5 公里）。远处或低功率的飞机可能确实听不见。
- 宽松窗口，严格评分：在漏斗中，宁可范围宽，因为真实的答案一旦被过滤掉就无法恢复。将判别留给精确模型评分。
- 按飞行对分组并保留镜像：直线匀速通过相对于其轨迹是对称的。单事件结果应至少报告两侧。
- ADS-B 间隙：无匹配 ≠ 不在该区域。低空通用航空飞机、军用飞机和一些直升机可能缺少 ADS-B 或未被接收（地形遮蔽、低空）。
- 声速和风：v 和 d 与 c 成比例缩放。有风时，c_eff = c + w·n；在强风（>8 m/s）中，加宽 d 窗口。
5. 报告
使用 assets/report_template.md 以用户语言输出报告。报告必须包含：
- 答案坐标
…
[为方便翻译而截断]