🙌

Juliaで人工衛星の見える時間を調べる

2021/02/20に公開

Julia

はじめに

最近Juliaの勉強を始めました。
なにか面白いことできないかなと思って、人工衛星周りの処理をやってみたくて調べるとSatellite Toolboxと言うライブラリがあったので簡単な使い方をまとめてみます。
Pythonとかだとpyorbital等便利なライブラリの記事いっぱいありますがJuliaだと日本語はおろか英語ですらほとんど情報なかったので勉強がてら少しずつ書いていこうと思います。
Juliaも衛星軌道も詳しくないので間違い等あればご指摘ください。

実行環境について

PCにJulia環境入れてIDEでも開発できるようにしていますが、軽く勉強するのだったらJupyterLab環境が楽そうだったのでDocker上に構築してブラウザでアクセスできるようにしました。
導入は下記記事を参考にしています。

Satellite Toolboxとは？

Juliaで軌道計算など人工衛星や天体に関する処理が行えるライブラリです。

下記記事によると開発はブラジル国立宇宙研究所が行っているため専門的な領域でも使えそうです。

ISSの可視時間を計算してみる

衛星が地平線から現れてから消えるまでの事をパスと言います。
地上のアンテナのある位置から通信できる時間とか変わってくるためとても基本的な情報なので、基本的にどのライブラリでも簡単に計算できる関数などが用意されています。
それではまず基本として国際宇宙ステーション(ISS)が自分のいる場所から見える時間帯を計算してみたいと思います。

TLEの取得

TLE(2行軌道要素形式)とはNASAを始め各国宇宙開発関連で標準的に利用されている、人工衛星の軌道要素を表すフォーマットです。

一部軍事衛星以外のTLEデータは北アメリカ航空宇宙防衛司令部(NORAD)により公開されており下記のサイトで検索することができます。
https://celestrak.com/

例えばISSの場合だと下記ページで閲覧でき、名前の通り2行(+名前)で記述されています。(2月20日記事執筆時点)

ISS (ZARYA)             
1 25544U 98067A   21050.88691362  .00002417  00000-0  52052-4 0  9993
2 25544  51.6444 202.9022 0003123  33.5384 127.6558 15.48985045270415

URLで指定している25544の部分がカタログ番号と呼ばれるもので各衛星に固有で割り振られています。
このフォーマットで計算される予想軌道自体は数週間もすると実際の軌道とずれてくるので、データは定期的に更新されます。
実際に利用する際には最新のデータを取得しましょう。

記述されている値の意味は理解しなくてもこれをライブラリに入力してあげればよしなに計算してくれるのでまずは上記ページから自動でデータを持ってこれるようにします。
HTTPで取得したいのでまずはHTTP通信用のライブラリを導入します。

using Pkg
Pkg.add("HTTP")
using HTTP

次にGETリクエストでTLEの文字列を取得します。

req = HTTP.request("GET", "https://celestrak.com/satcat/tle.php?CATNR=25544")
tle_string =String(req.body)

Satellite Toolbox導入

次にSatellite Toolboxのライブラリを導入します。

Pkg.add("SatelliteToolbox")
using SatelliteToolbox

ライブラリにTLEデータを読み込ませます。

tle = read_tle_from_string(tle_string)

tleは次のようになっています。

1-element Array{TLE,1}:
 TLE: ISS (ZARYA) (Epoch = 2021-02-19T08:33:35.794)

TLEは単一の衛星データだけでなく複数のものが入っていることもあるため配列として格納されています。
今回はISSのカタログ番号を指定しているため一件です。

地上の位置情報を調べる

次に地球上の特定の地点から見える位置を調べたいので、その場所の情報が必要です。
緯度・経度・高度の情報が必要なのですが、幸い現代ではGoogle Earthを使用することで簡単に調べることができます。
今回は自分が福岡に住んでいるので博多駅屋上の広場から見ると想定しましょう。
Google Earthを開いて左側のメニューから観測地点を検索します。
住所や地点を確定するとその場所が表示されます。デフォルトではその地点を中心に画像が回転するのでクリックするなどして回転を止めます。
観測地点に見えている状態でその場所の上にマウスカーソルを持っていくと下記画像のように右下に緯度、経度、標高が表示されます。

今回の場合

緯度 33.3523°
経度 130.2512°
標高 60m
ということがわかりましので変数として格納しましょう。
経緯はradianで使用するので先に変換しておきます。

lat = deg2rad(33.3523) #緯度
lon = deg2rad(130.2512) #経度
alt = 60.0 #標高

これで観測地の情報が揃いました。

パスの取得

それでは実際にパスを計算していきましょう。
Satellite Toolboxでパスを取得するにはground_station_accessesという関数を使用すれば良いようです。
GroundStation(地上局)というのは衛星と通信するためのアンテナの置いてある施設のことです。
地上局と衛星が通信できる時間帯を求める関数ということですね。

リンクに関数の説明が載っているのですが、知識が足りないためよくわからないパラメータもあります。

ground_station_accesses(orbp, [(WGS84)], Δt, ECI, ECEF, vargs...; kwargs...)

引数のうちorbpはOrbit Propagatorと言うもので軌道要素のことでしょうか。
今回はTLEから求めるのでinit_orbit_propagator関数を使用して下記のように書けば良いようです。

orbp = init_orbit_propagator(Val{:sgp4}, tle[1])

SGP4

SGP4というのは軌道計算アルゴリズムの一つでTLEで標準的に使用されているものです。
Juliaでは配列は1始まりなので先程求めた配列の1つ目をデータとして渡してpropagatorを取得しています。

(WGS84)

2つ目の(WGS84)と書いているものは地上局の情報です。
WGS84とは測地系の一種でGPS等で用いられている一般的な規格です。
地理座標系なので緯度経度で地点を表します。()になっているのはタプルなので、今回の場合は先ほど求めた観測地点の情報を(lat,lon,alt)のように入れてあげればいいです。

Δt

3つ目のΔt、これはどの程度先の情報まで取得するかの時間です。
秒数なので1日分取得するとしてt=1*24*60*60にしておきましょう。

ECI, ECEF, Args

ECIとECEFは座標系のことでそれぞれ日本語で「地球中心座標系」、「地球中心固定座標系」の事です。

ECI: Earth-Centered Inertial frame in which the state vector of the propagator is represented.
ECEF: Earth-Centered, Earth-fixed frame to be used for the analysis. It must be the same frame used to compute the ground station position vector.
vargs...: list of additional arguments to be passed to the function rECItoECEF when converting the ECI frame to the ECEF.

ECI,
ECEF

座標系については下記ページに日本語でわかり易く説明があります。

ECIは衛星側の座標系でECEFは地上局側の座標系を示します。
WikipediaによるとTLEで使用されている座標系はTEMEなのでTEME()を渡してあげれば良いようです。
ECEFなのですが、今回はWGS84の地理座標系(経緯)で渡していてこの時どう書けば良いのかわからず、この関数を使用している例もウェブ上で一件しか見つからなかったためそれのとおりにPEF()と記述しした。(Pseudo-Earth Fixed reference frame？疑似地球固定座標って訳せば良いのですかね、情報が出てこない……)
【追記】開発者の書いているサイトに下記の説明がありました。

For this example, we will convert TEME into the International Terrestrial Reference Frame (ITRF) for more accurate computation. This kind of conversion requires the Earth Orientation Data (EOP) that is provided by IERS. SatelliteToolbox.jl can easily load and use this data as follows:

The EOP data is measured and published by IERS. This means that the ITRF cannot be used here to predict what will be the ISS position in the future since the EOP data will not be available. In this scenario, if high accuracy is not needed, then the Pseudo-Earth Fixed (PEF) reference frame can be used. The conversion between TEME and PEF using IAU-76/FK5 does not require external data.
https://www.ronanarraes.com/2019/01/the-satellite-toolbox-for-julia/

どうやらEOPというデータが利用できるときはITEF使ったほうがいいようですが、高精度な予測が不要な場合はPEFを指定しておいても問題なさそうです。

Argsは特に何も書かなくても良さそうです。

Keywords

最後にキーワード
細かい設定を名前付き引数で書けるようです。

θ: Minimum elevation angle for communication between the satellite and the ground stations [rad]. (Default = 10ᵒ)
reduction: A function that receives a boolean vector with the visibility between the satellite and each ground station. It must return a boolean value indicating if the access must be computed or not. This is useful to merge access time between two or more stations. (Default = v->|(v...) i.e. compute the access if at least one ground station is visible)

θには角度で何度以上だと見えていることにするのか設定できます。デフォルトは10度以上ですが今回はθ=deg2rad(5)として地平線から5度出るタイミングにしておきます。
引数にギリシャ文字θを使用するとは中々斬新ですね。Juliaだとそういう文化あるのでしょうか？
~~(IMEのない英語圏のユーザーだとどうするんですかね)~~
【追記】JupyterLab、IntelliJ IDEAのJuliaプラグインともに\thetaと入力してTabキー押すと記号入力できるんですね！
その他数学記号も同様に入力することができます。便利だ。

計算

さてやっと引数について確認できたので実際に計算します。
改めて引数設定を合わせて記述します。

lat = deg2rad(33.3523) #緯度
lon = deg2rad(130.2512) #経度
alt = 60.0 #標高
orbp = init_orbit_propagator(Val{:sgp4}, tle[1])
ground_station_accesses(orbp, (lat, lon, alt), 1*24*60*60, TEME(), PEF(), θ=0)

結果として下記のような二次元配列が取得されました。
各行1列目が見え始める時間、2列目が消える時間です。

5×2 Array{DateTime,2}:
 2021-02-20T04:51:32.515  2021-02-20T04:59:57.61
 2021-02-20T19:55:15.013  2021-02-20T20:03:19.451
 2021-02-20T21:32:30.259  2021-02-20T21:39:48.836
 2021-02-21T02:28:21.231  2021-02-21T02:33:44.935
 2021-02-21T04:04:05.737  2021-02-21T04:12:36.921

今回は24時間の間に5回5度以上の軌道で出現することがわかりました。
時間はUTCで出てきます。
日本時間(JST)に変換したい場合はTimeZonesライブラリ等を使って変換しましょう。

おわりに

今回はここらへんまでにしておきます。
初めにも書きましたが情報がなさすぎて辛いですね。
個人的にはJuliaはこのような軌道計算などの分野にとても向いていると思うので、今後より使われていくことを期待して少しでも日本語情報を書いていければと思います。
(現在宇宙分野ではSTKと言う商用ツールやNASAのGMATと言うツールが使用されており、それらではMATLABを利用しているようです。)

Pythonのpyorbitalやephem等の同等の軌道計算ライブラリだとパスの時間と一緒に方角や最大仰角等も一緒に取得することができるのですが、Satellite Toolboxのground_station_accessesだとそれはできないようです。
次回はそのような情報の取得や軌道のプロットなどについてを調査・説明していきたいと思います。
何か質問や指摘がございましたらご気軽にコメントください。