人工衛星の回帰日数計算ツールを作成
背景
人工衛星の軌道にはいくつかの種類がありますが、その中で、1日で元の地点の上空に戻ってくる回帰軌道と、ある程度の日数をかけて元の地点の上空に戻る準回帰軌道があり、戻るまでの日数を回帰日数(Repeat Cycle)と呼びます*1*2。
回帰日数は、Web掲載されていることも多いのですが、未掲載の場合もあり、故あって、回帰日数を調べる必要があり、回帰日数計算ツールを作成しました*3。
自前で軌道計算 => 失敗
当初、人工衛星の軌道情報である2行軌道要素形式(TLE: Two Line Element*4 )から衛星の軌道を計算しようと、Webで見つかった計算方法*5をExcelで試してみたのですが、何故か、何回確認しても計算が合わず、困っていました*6。
ライブラリで軌道計算
そこで、無料で利用できる軌道計算ソフトウェアやライブラリがないか調べていたところ、丁度良いライブラリpyehemが見つかりました。
pyehemは、衛星のTLEを与えると、指定日時での衛星の位置を計算でき、衛星直下の緯度(longitude)と経度(latitude)なども分かります。
回帰日数の計算方法
衛星が地球を一周する時間はTLEから分かるので、当初の日時から何回か回った後で緯度・経度を求め、当初の緯度・経度との差がある程度小さければ、元の地点の上空に戻ってきたと分かるはずです。なお、緯度は、何回周回しても当初と同じとの想定です。
上記想定で、回帰日数計算ツールを作ってみたのですが、衛星が何回か回ると、緯度が少しずつズレてくることが分かりました。1周当たりの緯度のズレは約0.2度以下ですが、1日で約3.2度もズレる為、補正が必要と分かりました。そこで、衛星が1周回る度に、当初の緯度と同じになる日時を計算して補正しました*7。
テスト結果：COSMO-SkyMed-4(回帰日数16日)
テストとして、回帰日数が既知のCOSMO-SkyMed-4(回帰日数16日)*8で試した結果が下記です。16日後の緯度・経度が、当初の緯度・経度の約250m以内(-0.2431Km)に収まっており、ツールの動作が確認できました。COSMO-SkyMed-4は、非常に高い精度で回帰しているようです。
COSMO-SKYMED 4
1 37216U 10060A 18258.17806865 -.00000011 00000-0 51372-5 0 9997
2 37216 97.8902 80.4717 0001356 76.0322 284.1025 14.82152713425139
Lat(+N) diff[deg] diff[Km] | Long(+E) diff[deg] diff[Km]
[2018-09-17 00:39:00 = 0.00(days)] 39.4522 +0.0000 +0.0000 | -96.8107 +0.0000 +0.0000
[2018-10-03 00:39:00 = 16.00(days)] 39.4522 -0.0000 -0.0008 | -96.8136 -0.0028 -0.2431
[2018-10-19 00:39:00 = 32.00(days)] 39.4522 -0.0000 -0.0004 | -96.8119 -0.0012 -0.1051
[2018-11-04 00:38:59 = 48.00(days)] 39.4522 -0.0000 -0.0004 | -96.8053 +0.0054 +0.4662
[2018-11-20 00:38:57 = 64.00(days)] 39.4521 -0.0000 -0.0011 | -96.7930 +0.0178 +1.5276
[2018-12-06 00:38:53 = 80.00(days)] 39.4522 -0.0000 -0.0008 | -96.7747 +0.0360 +3.0975
テスト結果：WorldView-1(回帰日数不明)
次に、回帰日数がWeb未掲載のWorldView-1で試した結果が下記です。回帰日数は119日と非常に長く、また、元の場所から500m程度ズレていることが分かます(-0.5007Km)。
WORLDVIEW-1 (WV-1)
1 32060U 07041A 18258.02295190 .00000790 00000-0 35109-4 0 9993
2 32060 97.3879 16.4069 0002228 57.6683 54.8987 15.24397549611548
Lat(+N) diff[deg] diff[Km] | Long(+E) diff[deg] diff[Km]
[2018-09-18 10:46:01 = 0.00(days)] -0.1840 +0.0000 +0.0000 | 40.9467 +0.0000 +0.0000
[2018-10-05 10:47:27 = 17.00(days)] -0.1840 +0.0000 +0.0001 | 40.5894 -0.3573 -39.7724
[2018-10-18 10:42:41 = 30.00(days)] -0.1840 +0.0000 +0.0000 | 41.7862 +0.8395 +93.4509
[2018-10-22 10:48:26 = 34.00(days)] -0.1840 +0.0000 +0.0000 | 40.3508 -0.5959 -66.3345
[2018-11-04 10:43:17 = 47.00(days)] -0.1840 +0.0000 +0.0000 | 41.6392 +0.6925 +77.0936
[2018-11-08 10:48:55 = 51.00(days)] -0.1840 +0.0000 +0.0000 | 40.2321 -0.7145 -79.5426
[2018-11-21 10:43:25 = 64.00(days)] -0.1840 +0.0000 +0.0000 | 41.6122 +0.6655 +74.0819
[2018-11-25 10:48:57 = 68.00(days)] -0.1840 +0.0000 +0.0000 | 40.2332 -0.7135 -79.4297
[2018-12-08 10:43:05 = 81.00(days)] -0.1840 +0.0000 +0.0000 | 41.7037 +0.7570 +84.2650
[2018-12-12 10:48:30 = 85.00(days)] -0.1840 +0.0000 +0.0000 | 40.3523 -0.5944 -66.1650
[2018-12-25 10:42:17 = 98.00(days)] -0.1840 +0.0000 +0.0000 | 41.9123 +0.9656 +107.4907
[2018-12-29 10:47:35 = 102.00(days)] -0.1840 +0.0000 +0.0000 | 40.5885 -0.3582 -39.8736
[2019-01-15 10:46:12 = 119.00(days)] -0.1840 +0.0000 +0.0000 | 40.9422 -0.0045 -0.5007
[2019-02-01 10:44:20 = 136.00(days)] -0.1840 +0.0000 +0.0000 | 41.4151 +0.4684 +52.1456
[2019-02-05 10:49:24 = 140.00(days)] -0.1840 +0.0000 +0.0000 | 40.1541 -0.7926 -88.2335
[2019-02-22 10:46:57 = 157.00(days)] -0.1840 +0.0000 +0.0000 | 40.7760 -0.1707 -19.0025
[2019-03-11 10:44:01 = 174.00(days)] -0.1840 +0.0000 +0.0000 | 41.5179 +0.5712 +63.5844
[2019-03-15 10:48:49 = 178.00(days)] -0.1840 +0.0000 +0.0000 | 40.3199 -0.6268 -69.7701
テスト結果：WorldView-1(回帰日数不明)、別TLE
更に、上記より少しだけ新しいTLE*9で試した結果が下記です。回帰日数が140日と更に長く、また、元の場所からのズレも9.4Kmと大きいことが分かります(+9.4144Km)。当初の場所と回帰場所のズレとして、どの程度のまで許容できるか気になり、Webで調べましたが、残念ながら基準は見つかりませんでした。もう少しズレが大きくても許容できるならば(+13.8680Km)、回帰日数は102日と考えられるので、何らかの基準が欲しいですね....
WORLDVIEW-1 (WV-1)
1 32060U 07041A 18258.81766689 .00000958 00000-0 41975-4 0 9990
2 32060 97.3884 17.1911 0002170 64.8045 86.2416 15.24400435611664
Lat(+N) diff[deg] diff[Km] | Long(+E) diff[deg] diff[Km]
[2018-09-18 10:46:01 = 0.00(days)] -0.2185 +0.0000 +0.0000 | 40.9431 +0.0000 +0.0000
[2018-10-05 10:47:22 = 17.00(days)] -0.2185 +0.0000 +0.0001 | 40.6084 -0.3347 -37.2598
[2018-10-18 10:42:28 = 30.00(days)] -0.2185 +0.0000 +0.0000 | 41.8378 +0.8947 +99.5998
[2018-10-22 10:48:10 = 34.00(days)] -0.2185 +0.0000 +0.0000 | 40.4152 -0.5279 -58.7642
[2018-11-04 10:42:50 = 47.00(days)] -0.2185 +0.0000 +0.0000 | 41.7539 +0.8108 +90.2596
[2018-11-08 10:48:24 = 51.00(days)] -0.2185 +0.0000 +0.0000 | 40.3651 -0.5780 -64.3473
[2018-11-21 10:42:37 = 64.00(days)] -0.2185 +0.0000 +0.0000 | 41.8132 +0.8701 +96.8600
[2018-11-25 10:48:03 = 68.00(days)] -0.2185 +0.0000 +0.0000 | 40.4580 -0.4851 -54.0061
[2018-12-12 10:47:09 = 85.00(days)] -0.2185 +0.0000 +0.0000 | 40.6926 -0.2505 -27.8884
[2018-12-29 10:45:41 = 102.00(days)] -0.2185 +0.0000 +0.0000 | 41.0677 +0.1246 +13.8680
[2019-01-15 10:43:40 = 119.00(days)] -0.2185 +0.0000 +0.0000 | 41.5835 +0.6404 +71.2906
[2019-01-19 10:48:40 = 123.00(days)] -0.2185 +0.0000 +0.0000 | 40.3357 -0.6074 -67.6159
[2019-02-05 10:45:57 = 140.00(days)] -0.2185 +0.0000 +0.0000 | 41.0277 +0.0846 +9.4144
[2019-02-22 10:42:39 = 157.00(days)] -0.2185 +0.0000 +0.0000 | 41.8637 +0.9206 +102.4796
[2019-02-26 10:47:21 = 161.00(days)] -0.2185 +0.0000 +0.0000 | 40.6916 -0.2515 -27.9979
[2019-03-15 10:43:21 = 178.00(days)] -0.2185 +0.0000 +0.0000 | 41.7050 +0.7618 +84.8073
テスト結果：COSMO-SkyMed-4(日々の軌道のズレ確認)
本ツールでは、経度方向の軌道のズレの許容値を指定できるようになっており、「1°」程度が適切だが、「10°」程度まで許容して計算すると、最初の軌道を基準とした東西方向への軌道のズレを、日々確認することもできる。下記がその結果(プラス符号が東方向)。
COSMO-SKYMED 4
1 37216U 10060A 18258.17806865 -.00000011 00000-0 51372-5 0 9997
2 37216 97.8902 80.4717 0001356 76.0322 284.1025 14.82152713425139
Lat(+N) diff[deg] diff[Km] | Long(+E) diff[deg] diff[Km]
[2018-09-17 00:39:00 = 0.00(days)] 39.4522 +0.0000 +0.0000 | -96.8107 +0.0000 +0.0000
[2018-09-18 00:57:13 = 1.01(days)] 39.4521 -0.0000 -0.0011 | -101.3680 -4.5573 -391.7246
[2018-09-19 01:15:27 = 2.03(days)] 39.4522 -0.0000 -0.0008 | -105.9252 -9.1145 -783.4468
[2018-09-21 00:14:42 = 3.98(days)] 39.4522 -0.0000 -0.0008 | -90.7359 +6.0748 +522.1689
[2018-09-22 00:32:55 = 5.00(days)] 39.4522 -0.0000 -0.0008 | -95.2931 +1.5176 +130.4502
[2018-09-23 00:51:09 = 6.01(days)] 39.4522 +0.0000 +0.0000 | -99.8503 -3.0395 -261.2661
[2018-09-24 01:09:23 = 7.02(days)] 39.4522 -0.0000 -0.0008 | -104.4074 -7.5967 -652.9820
[2018-09-26 00:08:37 = 8.98(days)] 39.4522 -0.0000 -0.0008 | -89.2179 +7.5928 +652.6479
[2018-09-27 00:26:51 = 9.99(days)] 39.4521 -0.0000 -0.0011 | -93.7750 +3.0357 +260.9368
[2018-09-28 00:45:05 = 11.00(days)] 39.4522 -0.0000 -0.0008 | -98.3321 -1.5214 -130.7731
[2018-09-29 01:03:19 = 12.02(days)] 39.4521 -0.0000 -0.0011 | -102.8892 -6.0785 -522.4813
*2:Repeat Cycle以外に、Revisit Timeという指標もあります。Revisit Timeは、同じ場所を再度撮影できるまでの期間であり、通常1日程度です。なお、撮影場所から見える衛星の方向は真上とは限らず、毎回変わるので、衛星方向を傾けて撮影します。そのため、撮影場所の見え方は異なります。参考：What is the difference between a repeat cycle and revisit time for a satellite?
*4:CelesTrak: Current NORAD Two-Line Element Setsから最新のTELを入手可能
*5:http://www.infra.kochi-tech.ac.jp/takagi/Geomatics/5Estimation2.pdf
*6:離心近点角 E をExcelのソルバーで求めた値E=200.9819218が、記載数値E=200.9785378と一致しない。つまり、記載のE、e、Mの値で計算すると方程式 「E − e sin E = M」が成立しない。具体的には、「E − e sin E - M = 200.9785378 - 0.0001679 * sin( RADIANS(200.9785378) ) - 200.9819819 = -0.0033840」となる。なお、E=200.9819218の場合は「E − e sin E - M = 0」となる。
*9:軌道情報の元期(起点)が、上記よりも約0.8日新しい。
