AlanBacker/audio-flyover-geolocation

聲學飛越地理定位（都卜勒 × ADS-B）

此技能複現並產品化了colsto影片 「我發現了一種洩漏你確切位置的（聲音，而且到處都是）」（2026-11-02，https://www.youtube.com/watch?v=UytYrhYLOKI）中的方法：

螺旋槳或旋翼的葉片通過頻率（BPF） 是一種穩定的窄帶音調。即使沒有人能聽到它，它也會在聲譜圖中留下一條平滑下降的都卜勒曲線。曲線的形狀給出了飛機的速度 v、其最近接近斜距 d 和CPA 時間 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. 工作流程（CLI）

慣例：AFG="python <技能目錄>/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 = 槳葉數 × RPM / 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 接收器。
• 只提取目標日期、目標時間窗口 ± 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. 事件 1 走廊：速度窗口，k = v²/d 一致性，斜距窗口，並且 CPA 不得落在軌跡端點。
2. 事件間偏移：tCPA2 − tCPA1，包括每個事件的 d/c 傳播延遲。
3. 靜默過濾器：在剩餘的安靜時間段內，沒有其他螺旋槳或旋翼飛機通過 --d-quiet 範圍內。

選項：

• 模式 B：添加 --rec-start <UTC> --rec-start-tol 2（或 30）。
• --quiet "a-b,c-d" 指定您確定未聽到任何聲音的錄音時間段（預設：事件之外的所有時間）。當錄音中有大量風聲或語音時，減少指定的時間段。靜默即是信號，但錯誤的靜默會扼殺正確答案。
• 預設網格為 200 米。如果候選項太多，請縮小 bbox，而不是增大網格。

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）

1. 疊加圖：模型曲線必須在兩條漸近線和通過 CPA 時都符合觀測結果，而不僅僅是一個區段。
2. 熱力圖：0.9 區域的半徑是要報告的誤差。當以下任一情況成立時，報告一條條帶：
   • 兩個事件：geometry.warning 被設定（交叉角小）；
   • 一個事件：coarse_plateau.length_m > 600 米。
3. 鏡像候選項和替代項：使用視覺線索（地平線、建築物、植被、道路噪音、回聲）、DEM 視域和太陽方位角來排除它們。
4. 否證：答案附近的其他飛越是否也應該被聽到？
   • 聽到但未檢測到 = 檢測閾值問題。
   • 在不應該聽到時出現 = 位置錯誤。

4. 精度規則（及其原因）

1. 發射時間：f_obs(t) = f0 / (1 + Ṙ(τ)/c)，其中 t = τ + R(τ)/c。使用接收時間的 Ṙ 會給出一條形狀錯誤的曲線，而不僅僅是偏移的。
2. 3D 斜距 + DEM：觀測者高程加上麥克風高度（預設 2 米）。將 500 米的地形起伏視為 0 會使 d 和 tCPA 都出現錯誤。
3. 高度基準：未經校正的壓力高度，加上 30 hPa 的 QNH 偏移 ≈ 250 米的高度誤差，這會直接影響 d。
4. 延遲的 CPA：觀察到 f = f0 的時刻是 tCPA_emit + d/c（1 公里 ≈ 2.9 秒）。事件偏移和時鐘都必須據此進行轉換。
5. k = v²/d 比單獨使用 v 或 d 更可靠。
   • 曲線中段的斜率 −f0·v²/(c·d) 被緊密鎖定。漸近線通常會被噪聲、閘控或飛機加速/爬升所污染。
   • 因此，搜尋使用一個寬的速度窗口和一個窄的 k 窗口，並且每個候選飛機自己的地速會預測其 d。
6. 鎖定評分：兩個事件共享一個 δ。否則每個事件會對齊自己的時間，導致錯誤的位置也能獲得高分。
7. 細化至 ≤100 米：200 米網格僅是粗篩；最終點來自子網格爬山法和熱力圖。
8. 靜默即是信號，但要保持閾值寬鬆（預設：螺旋槳 3 公里，旋翼 2.5 公里）。遠距離或低功率的飛機可能真的聽不到。
9. 寬鬆窗口，嚴格評分：在漏斗中，寧可放寬範圍，因為真正的答案如果被過濾掉就永遠無法復原。將區分任務留給精確模型評分。
10. 按飛行對分組並保留鏡像：直線等速通過相對於其軌跡是對稱的。單事件結果應至少報告兩側。
11. ADS-B 間隙：沒有匹配 ≠ 不在該區域。低空通用航空、軍用飛機和一些直升機可能缺乏 ADS-B 或未被接收到（地形遮蔽、低海拔）。
12. 聲速和風：v 和 d 與 c 成比例縮放。有風時，c_eff = c + w·n；在強風中（>8 米/秒），加寬 d 窗口。

5. 報告

使用 assets/report_template.md 輸出，語言為使用者的語言。它必須包含：

• 答案座標

...（為翻譯而截斷）