AlanBacker/audio-flyover-geolocation

音響フライオーバー位置特定（ドップラー×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 + 音速伝播時間 + 放射時の動径速度）、候補を絞り込みます。2回のフライオーバー間の時間間隔（動画では920秒）と無音制約（他の時間帯に他の航空機が聞こえないこと）を加えることで、単一の地点にまで絞り込まれます。最後に、「ロックされた2イベントスコア」（共有時間オフセット、イベントごとの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データ量が膨大になり管理不能
地域	国または州レベルで十分；都市レベルが事前情報としてあれば高速	全国検索も可能だが、時刻が既知の1イベントまたは2イベントが必要
録音時計（開始時刻）	±2秒以内で既知→1イベントで地点特定可能；±30秒以内→1イベントでは航跡に沿った帯域のみ	2イベントを使用；その時間間隔が絶対時刻の代わりになる
ADS-Bカバレッジ	対象航空機がADS-B Outを搭載し、地上局で受信されていること	米国14 CFR 91.225: クラスB/C空域外かつモードCベール外の10,000フィート以下の一般航空機はADS-B Out非搭載可→データに含まれない可能性あり
音声処理チェーン	オリジナルまたは軽度圧縮が最適；電話コーデックでも通常はトーン維持	強いノイズ除去や編集により、曲線が切断されたりイベント間のギャップが破壊される可能性あり

作業開始前に、トリアージ結果と予想精度をユーザーに伝えてください。

2. モード選択

使用可能な曲線の数？
├─ ≥2（同一連続録音内）→ モードA: 2イベント、絶対時刻不明（動画の主要手法；最もロバスト）
│     2つの航跡が約25°未満で交差する場合、結果は航跡に沿った帯域に劣化する（スクリプトが警告）
├─ 1
│   ├─ 録音時刻が数秒以内で既知 → モードB: 1イベント、時刻既知（動画の後半；米国全域でも最優先）
│   │     等速直線飛行は鏡像対称→常に「鏡像候補」が存在する；
│   │     2番目の手がかり（視覚、地形）が必要
│   └─ 時刻が±30秒程度の精度のみ → モードB′: 長さ≈ 2·v·Δtの帯域のみ；そのように報告
├─ 航空機曲線なし、ただし雷＋日付既知 → モードT: 雷鳴エンベロープ×落雷データ（kmレベル、参照/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.json が空の場合、「フライオーバーなし」を意味するのは boxes.diag.json の status が ok の場合のみ。limited（キャリブレーション済みノイズモデルに対して判断された短いクリップ、またはほとんどがゲートされた音声）および insufficient（--min-dur より短いクリップ、またはほぼすべてがゲートされた）には reasons が付属：これを読み、rec.png を確認し、短いクリップでは --min-dur を下げてください。
• 各イベントには trend（descending または flat）と drop_pct が含まれます。周波数が上昇するリッジはフライオーバーではありません：それは boxes.diag.json の rejected にリストされ、boxes.json には含まれません。イベントに warnings がある場合、そのリッジは1回の通過ではありえない動き（再上昇、または45%超の下降）をしていることを意味します。おそらく2機の航空機または2つの高調波を結合しています。plot で確認し、各部分を独自の狭いボックスでトレースしてください。
• --win: 近くて急峻な曲線には0.5～1秒、遠くてかすかな曲線には2～4秒を使用。かすかな曲線の場合は、plot でズームインして目視で確認してください。その後、1つの高調波のみを含む狭いボックスを 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: 1日1回のグローバルリリース、無料、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

漏斗は3段階で構成:

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 m。候補が多すぎる場合は、グリッドを拡大するのではなく、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、差分位置ではない）から放射時刻における正確なモデル曲線を生成し、スペクトログラムと比較します。1つの時間オフセットδが共有され、これにより2つのイベント間の物理的な時間間隔が固定されます。各イベントのf0は±0.6%以内で微調整されます。
• 精緻化: 保持された各ノードは、縮小する3×3の山登り法（100→25 m）を受けます。
• グループ化: 結果はフライトペアごとにグループ化されます。各ペアは500 m以上離れた最大3つの仮説を保持するため、両方の鏡像側が生き残ります。
• 出力:
  • 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を通して観測値に適合しなければなりません（1つのセグメントだけでは不十分）。
2. ヒートマップ: 0.9領域の半径が報告すべき誤差です。以下のいずれかに該当する場合は帯域として報告:
   • 2イベント: geometry.warning が設定されている（交差角度が小さい）；
   • 1イベント: coarse_plateau.length_m > 600 m。
3. 鏡像候補と代替案: 視覚的手がかり（地平線、建物、植生、道路騒音、エコー）、DEMビューズヘッド、太陽方位を用いて除外。
4. 反証: 答えの近くで他のフライオーバーも聞こえるはずか？
   • 聞こえるはずだが検出されなかった = 検出閾値の問題。
   • 聞こえるはずがないのに存在する = 位置が間違っている。

4. 精度ルール（およびその理由）

1. 放射時刻: f_obs(t) = f0 / (1 + Ṙ(τ)/c)、ここで t = τ + R(τ)/c。受信時刻のṘを使用すると、シフトされただけでなく形状が異なる曲線になります。
2. 3D斜距離 + DEM: 観測者標高＋マイク高（デフォルト2 m）。500 mの地形起伏を0として扱うと、dとtCPAの両方が誤ります。
3. 高度データム: 30 hPaのQNHオフセットを補正しない気圧高度は約250 mの高度誤差となり、これがdに直接影響します。
4. 遅延CPA: f = f0が観測される瞬間は tCPA_emit + d/c（1 km ≈ 2.9 s）。イベントオフセットと時計は両方ともこれを用いて変換する必要があります。
5. k = v²/d はvやd単独よりも信頼性が高い。
   • 曲線中間の傾き −f0·v²/(c·d) は正確に決定されます。漸近線はノイズ、ゲーティング、または航空機の加速・上昇によって汚染されることがよくあります。
   • そのため、探索では広い速度ウィンドウと狭いkウィンドウを使用し、各候補航空機の対地速度がそのdを予測します。
6. ロックスコアリング: 両方のイベントが1つのδを共有します。そうしないと、各イベントが独自の時刻に合わせられ、誤った位置も高いスコアを獲得します。
7. ≤100 mへの精緻化: 200 mグリッドは粗いスクリーニングに過ぎません。最終地点はサブグリッドの山登り法とヒートマップから得られます。
8. 無音は信号ですが、閾値は寛大に保ってください（デフォルト: プロペラ機3 km、ローター機2.5 km）。遠方または低出力の航空機は実際に聞こえない場合があります。
9. 緩いウィンドウ、厳格なスコアリング: 漏斗では広めに取ってください。真の答えがフィルタリングされると回復できません。識別は正確なモデルスコアリングに委ねてください。
10. フライトペアでグループ化し、鏡像を保持: 等速直線飛行はその航跡に関して対称です。単一イベントの結果は少なくとも両側を報告する必要があります。
11. ADS-Bギャップ: マッチがない ≠ そのエリアにいない。低高度の一般航空機、軍用機、一部のヘリコプターはADS-Bを搭載していないか、受信されない場合があります（地形による遮蔽、低高度）。
12. 音速と風: vとdはcに比例します。風がある場合、c_eff = c + w·n；強風（>8 m/s）ではdウィンドウを広げてください。

5. レポート

assets/report_template.md をユーザーの言語で出力します。以下を含める必要があります:

• 答えの座標

…

[翻訳のため一部省略]