Skyfieldを使った地上駅パススケジュールの作成

ローカルパススケジューラが必要な理由

SatNOGS ステーション、アマチュア無線シャックトラッキング AO-91/ISS/NOAA、または初期段階のキューブサット地上セグメントを運営している場合は、自動パス予測が必要です。商用ツール (STK、NOVA) はやり過ぎです。無料のウェブプレディクター (AMSAT、Heavens-Above) はうまくスクリプト化されません。STK Cloud が2026年3月にサービスを終了する中、最もクリーンで無料の代替手段は、PyEphemの後継である最新のPythonである Skyfieldです。

この記事では、約80行のPythonで完全な地上局スケジューラについて説明します。新しいTLEを取得し、衛星リストの次の24時間のパスを計算し、最小標高でフィルタリングし、ダウンリンク頻度のドップラーを計算し、パスがカレンダーアプリに表示されるようにiCalendarを送信します。

インストール

pip install skyfield requests icalendar

それでおしまい。スカイフィールドは初回実行時に独自のエフェメライド (DE421、IERS) を取得します。約 17 MB がローカルにキャッシュされます。

ステップ 1: 現在の TLE を取得する

TLE ドリフト — エポックから 7 日以内に使用してください。セレストラックから現在のエレメントセットを取り出してください:

from skyfield.api import load

stations_url = 'https://celestrak.org/NORAD/elements/gp.php?GROUP=amateur&FORMAT=tle'
sats = load.tle_file(stations_url)
by_name = {s.name: s for s in sats}
print(f'Loaded {len(sats)} amateur-radio satellites')

ターゲットを名前で選択します（Celestrakファイルでは大文字と小文字が区別されます）。

targets = [
    by_name['AO-91 (FOX-1B)'],
    by_name['ISS (ZARYA)'],
    by_name['NOAA 19'],
    by_name['GreenCube (IO-117)'],
]

ステップ 2: 地上局の定義

from skyfield.api import wgs84, load

ts = load.timescale()
my_station = wgs84.latlon(40.0150, -105.2705, elevation_m=1624)  # Boulder, CO
min_elevation_deg = 10  # horizon mask to skip low passes

elevation_mパラメータを使用してください。これは、低い仰角での傾斜範囲に大きく影響します。

ステップ 3: 上昇/最高潮/セットイベントを検索する

Satellite.find_events()はイベントのフラットリスト (0 = 上昇、1 = クライミング、2 = セット) を返します。これらを 3 つに分けて処理します。

from datetime import timedelta

t0 = ts.now()
t1 = ts.from_datetime(t0.utc_datetime() + timedelta(hours=24))

for sat in targets:
    t, events = sat.find_events(my_station, t0, t1,
                                altitude_degrees=min_elevation_deg)
    # Group into rise/peak/set triples
    triples = zip(t[0::3], t[1::3], t[2::3])
    for rise, culm, setting in triples:
        # Peak elevation
        topo = (sat - my_station).at(culm)
        alt, az, dist = topo.altaz()
        print(f'{sat.name:30s}  {rise.utc_iso():20s}  peak {alt.degrees:4.1f}°  slant {dist.km:4.0f} km')

t0ですでに地平線上にある衛星は無視してください（イベントリストは最高潮とセットのペアから始まります）。

ステップ 4: パス上でのドップラーシフト

アマチュアのUHF (435 MHz) では、パスによってキャリアが±10 kHzずれる可能性があります。レシーバーはトラッキングする必要があります。レンジレートから瞬時ドップラーを計算:

import numpy as np

FREQ_HZ = 435_800_000  # AO-91 downlink
C = 299_792_458.0      # m/s

def doppler_shift(sat, observer, t):
    """Positive when satellite approaching (received frequency is higher)."""
    # Numerical range rate from 1s apart
    dt = 1.0  # seconds
    t_next = ts.from_datetime(t.utc_datetime() + timedelta(seconds=dt))
    r0 = (sat - observer).at(t).distance().m
    r1 = (sat - observer).at(t_next).distance().m
    range_rate = (r1 - r0) / dt  # m/s, positive = moving away
    return -FREQ_HZ * range_rate / C  # Hz, positive = approaching

# Sample across a pass
for pct in [0, 25, 50, 75, 100]:
    t_sample = ts.from_datetime(
        rise.utc_datetime() + (setting.utc_datetime() - rise.utc_datetime()) * pct / 100
    )
    d = doppler_shift(sat, my_station, t_sample)
    print(f'  {pct:3d}% pass:  doppler = {d:+8.0f} Hz')

一般的な 6 分間の AO-91 オーバーヘッドパスでは、ドップラーが約 +10 kHz（上昇）から 0（天頂）、-10 kHz（セット）まで流れるのがわかります。これをラジオのメイン VFO でトラッキングすれば、信号が途切れることはありません。

ステップ 5: 各パスのリンクバジェット

傾斜範囲がわかったら、それを rftools 自由空間経路損失計算ツールに入力するか、インラインで計算します。

def fspl_db(distance_m, freq_hz):
    return 20 * np.log10(4 * np.pi * distance_m * freq_hz / C)

# At zenith of an AO-91 pass (~800 km slant at 40° elevation from Boulder)
fspl = fspl_db(800_000, FREQ_HZ)
print(f'FSPL = {fspl:.1f} dB')  # ~143 dB

ITU-R 伝搬モデルによるモンテカルロリンクバジェットの完全版については、衛星リンクバジェットアナライザーを使用し、シナリオの URL をステーションログブックと共有してください。

ステップ 6: iCalendar を発行してカレンダーにパスが表示されるようにする

これがスケジューラが人間にとって役立つところです。

from icalendar import Calendar, Event
from datetime import datetime
import pytz

cal = Calendar()
cal.add('prodid', '-//rftools ground station//')
cal.add('version', '2.0')

for sat in targets:
    t, events = sat.find_events(my_station, t0, t1, altitude_degrees=10)
    triples = zip(t[0::3], t[1::3], t[2::3])
    for rise, culm, setting in triples:
        topo = (sat - my_station).at(culm)
        alt, az, _ = topo.altaz()
        ev = Event()
        ev.add('summary', f'{sat.name} (peak {alt.degrees:.0f}°)')
        ev.add('dtstart', rise.utc_datetime())
        ev.add('dtend', setting.utc_datetime())
        ev.add('description', f'Culminate at {culm.utc_iso()}, AZ {az.degrees:.0f}°')
        cal.add_component(ev)

with open('passes.ics', 'wb') as f:
    f.write(cal.to_ical())
print('Wrote passes.ics — import into Google Calendar / Apple Calendar / Outlook')

passes.icsを Apple カレンダー/Google カレンダー/Outlook にインポートすると、各パスの5分前にプッシュ通知が届きます。

ステップ 7: cron で実行してください

~/bin/station-passes.pyにスクリプトをドロップして、1 日に 1 回 cron してください:

# Refresh pass schedule nightly at 03:00 local
0 3 * * * /usr/bin/python3 ~/bin/station-passes.py > ~/passes.ics 2>> ~/station.log

curl --upload-fileとWebdavカレンダーをペアリングすると、携帯電話にすべてのパスが自動的に表示されます。

STKクラウドのパーツごとの交換

アマチュア/キューブサット/小規模の商業事業では、統合ツールキットがSTK Cloudの日常的な機能に取って代わります。

ミッションデザインの貿易研究や機関プログラムについては、やはりデスクトップSTKまたはNASA GMATが必要です。それ以外は、AMSATの設計レビューを含め、すべてSkyfield + rftoolsスタックが揃っています。

さらに読む

• スカイフィールド・ドキュメンテーション
• セレストラック TLE カタログ
• rftools 衛星パス予測ツール— 最も一般的なアマチュア衛星用のこのスケジューラのホストバージョン
• STK クラウドからの移行:無料の代替案— STK ワークフロー全体の置き換えに関する関連記事
• UHF CubeSat リンクバジェットのサイズ設定— エンドツーエンドの計算例