人工衛星や国際宇宙ステーション（ISS）が今どこにいて、今の自分の位置から見えるか（可視）かどうかを計算する機会は多くあると思います。
今回は、地球を周回する人工衛星が地上局から見えている時間（可視時間）を幾何学的に計算し、Pythonでシミュレーションする方法を解説します。
本記事の内容
- 衛星と地上局の可視条件の幾何学
- 仰角と可視弧の計算式の導出
- 可視時間に影響するパラメータ
- Pythonでの可視時間シミュレーション
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
可視条件の幾何学
地上局から衛星が見えるためには、衛星が地上局の「地平線」よりも高い位置にある必要があります。これを仰角（Elevation angle）$\varepsilon$ で定量化します。
仰角の定義
地上局の水平面と衛星方向のなす角を仰角と呼びます。仰角が正のとき衛星は地平線上にあり、一般に通信を行うには最低仰角 $\varepsilon_{\min}$（通常5度から10度）以上が必要です。
幾何学的関係
地球の中心を原点とし、地上局の位置ベクトルを $\bm{r}_g$、衛星の位置ベクトルを $\bm{r}_s$ とします。
地球の半径を $R_E$、衛星の軌道高度を $h$ とすると、円軌道の場合に衛星が可視となる最大の天頂角 $\eta_{\max}$ は、
$$ \sin \eta_{\max} = \frac{R_E}{R_E + h} \cos \varepsilon_{\min} $$
で与えられます。ここで天頂角 $\eta$ は、衛星から見た地球中心方向と地上局方向のなす角です。
地心角
衛星のサブサテライトポイント（直下点）と地上局の間の地心角 $\rho$ は、
$$ \cos \rho = \sin \varphi_g \sin \varphi_s + \cos \varphi_g \cos \varphi_s \cos(\lambda_s – \lambda_g) $$
ここで $(\varphi_g, \lambda_g)$ は地上局の緯度・経度、$(\varphi_s, \lambda_s)$ は衛星の直下点の緯度・経度です。
可視条件は $\rho \leq \rho_{\max}$ であり、最大地心角は、
$$ \rho_{\max} = \arccos\left(\frac{R_E}{R_E + h} \cos \varepsilon_{\min}\right) – \varepsilon_{\min} $$
で計算できます。
仰角の計算
衛星と地上局の間の仰角 $\varepsilon$ は地心角 $\rho$ から次のように求まります。
$$ \tan \varepsilon = \frac{\cos \rho – \frac{R_E}{R_E + h}}{\sin \rho} $$
可視時間の理論
単一パスの可視時間
円軌道の衛星が地上局の上空を通過する際の可視時間は、衛星の軌道速度と可視弧の大きさで決まります。最大パスの場合（衛星が地上局の真上を通過する場合）、可視時間 $T_{\text{vis}}$ は次のようになります。
$$ T_{\text{vis}} = \frac{2 \rho_{\max}}{n} $$
ここで $n = 2\pi / T$ は衛星の平均運動（角速度）、$T$ は軌道周期です。
高度と可視時間の関係
軌道高度が高いほど可視範囲（フットプリント）が広がるため、1回のパスでの可視時間は長くなります。ただし、軌道周期も長くなるため、1日あたりのパス数は減少します。
Pythonでの可視時間シミュレーション
import numpy as np
import matplotlib.pyplot as plt
# 定数
Re = 6371.0 # 地球の平均半径 [km]
def max_geocentric_angle(h, eps_min_deg):
"""可視範囲の最大地心角を計算"""
eps_min = np.radians(eps_min_deg)
rho_max = np.arccos(Re / (Re + h) * np.cos(eps_min)) - eps_min
return rho_max
def max_visibility_time(h, eps_min_deg):
"""最大可視時間を計算（真上通過の場合）"""
mu = 398600.4418 # [km^3/s^2]
a = Re + h
T = 2 * np.pi * np.sqrt(a**3 / mu) # 軌道周期 [s]
n = 2 * np.pi / T # 平均運動 [rad/s]
rho_max = max_geocentric_angle(h, eps_min_deg)
t_vis = 2 * rho_max / n
return t_vis
def elevation_profile(h, pass_distance_km, n_points=200):
"""パスの仰角プロファイルを計算
pass_distance_km: 最接近時の地上距離 [km]
"""
a = Re + h
mu = 398600.4418
T = 2 * np.pi * np.sqrt(a**3 / mu)
n = 2 * np.pi / T
rho_max = max_geocentric_angle(h, 0)
t_half = rho_max / n
times = np.linspace(-t_half, t_half, n_points)
elevations = []
for t in times:
# 衛星の地心角（簡易モデル）
rho = np.sqrt((n * t)**2 + (pass_distance_km / Re)**2)
if rho > rho_max:
elevations.append(0)
else:
eps = np.arctan2(np.cos(rho) - Re / (Re + h), np.sin(rho))
elevations.append(np.degrees(eps))
return times, elevations
# 高度と最大可視時間の関係
altitudes = np.linspace(200, 2000, 100)
eps_mins = [0, 5, 10, 15]
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
for eps_min in eps_mins:
t_vis = [max_visibility_time(h, eps_min) / 60 for h in altitudes]
axes[0].plot(altitudes, t_vis, label=f'eps_min = {eps_min} deg')
axes[0].set_xlabel('Altitude [km]')
axes[0].set_ylabel('Max Visibility Time [min]')
axes[0].set_title('Maximum Visibility Time vs Altitude')
axes[0].legend()
axes[0].grid(True)
# 仰角プロファイル
h = 500 # km
for d_km in [0, 200, 500, 1000]:
t, el = elevation_profile(h, d_km)
axes[1].plot(t, el, label=f'd = {d_km} km')
axes[1].set_xlabel('Time [s]')
axes[1].set_ylabel('Elevation [deg]')
axes[1].set_title(f'Elevation Profile (h = {h} km)')
axes[1].legend()
axes[1].grid(True)
axes[1].set_ylim(0, 90)
plt.tight_layout()
plt.show()
# 具体的な数値を表示
print("=== 高度ごとの最大可視時間 ===")
for h in [300, 500, 700, 1000, 1500]:
for eps_min in [5, 10]:
t = max_visibility_time(h, eps_min)
rho = np.degrees(max_geocentric_angle(h, eps_min))
print(f" 高度 {h:5d} km, 最低仰角 {eps_min:2d} deg: "
f"可視時間 {t/60:.1f} min, 最大地心角 {rho:.1f} deg")
まとめ
本記事では、人工衛星と地上局の可視時間の計算方法について解説しました。
- 可視条件は仰角 $\varepsilon \geq \varepsilon_{\min}$ で定量化される
- 最大地心角 $\rho_{\max}$ から可視範囲を幾何学的に求められる
- 軌道高度が高いほど1パスの可視時間は長くなるが、パス数は減る
- 最低仰角の設定は通信品質と可視時間のトレードオフとなる
次のステップとして、以下の記事も参考にしてください。