軌道上観測による微小デブリ環境推定手法の構築

Kyushu University Institutional Repository

軌道上観測による微小デブリ環境推定手法の構築

古本, 政博

http://hdl.handle.net/2324/1959134

出版情報：Kyushu University, 2018, 博士（工学）, 課程博士 バージョン：

権利関係：

博⼠学位論⽂

軌道上観測による微⼩デブリ環境推定⼿法の構築

古本 政博

九州⼤学⼤学院 ⼯学府 航空宇宙⼯学専攻

2018

論⽂要旨

本論⽂は，地上から観測できない微⼩なスペースデブリの環境を軌道上観測データを⽤いて推定する

⼿法である，動的環境推定モデルの提案および検証を⾏うものである．

軌道上には，地上からの観測が不可能な微⼩なスペースデブリが存在している．そのような微⼩デブ リは宇宙機に対して致命的な損傷を与え得るにも関わらず，各国の提供する環境モデルの推定結果は⼀

致していない．モデルごとの⼿法の違いだけでなく，観測の困難さによるデータの不⾜がこのことに影 響していると考えられる．また，破砕事象によって⼤量に⽣成される微⼩デブリの継続的な環境把握は，

⽇々変動する宇宙環境の変化を捉える上で重要である．そのため微⼩デブリを軌道上で観測する計画が 各国で進められており，九州⼤学で2011年に発⾜されたIDEA Projectもその⼀例である．そこで本論⽂

では，リアルタイムにデータ取得が可能である衛星による軌道上観測の利点を活かし，常に正確かつ迅 速に微⼩デブリ環境を把握するため，観測データを得られる度にモデルを逐次更新可能である動的な環 境推定モデルの構築を⾏った．

本論⽂ではまず，動的環境推定モデルを構築するため，観測衛星により検出され得る微⼩デブリの軌 道の特徴に関する２つの新たな⼿法を提案した．１つ⽬の拘束⽅程式は観測衛星に衝突し得るデブリの 軌道⾯について単純な⽅程式で拘束条件を与えるものであり，衝突データの解析で対象とすべき軌道を 絞り込むことを可能とする．もう１つのトーラスモデルは２つの軌道間の衝突フラックスを近似的に計 算するものである．いずれの⼿法も衝突位置と軌道要素のみから計算が可能であり，軌道を分割して逐 次的に計算を⾏う必要のある既存の⼿法と⽐較して計算の効率化とアルゴリズムの単純化が可能とな った．さらに，既存の数値シミュレーションとの⽐較を⾏うことで提案した⼿法の妥当性を⽰した．

次に，動的環境推定モデルの理論的基礎として，データ同化を導⼊した．データ同化は観測データを

⽤いて推定結果を逐次更新する⼀般的な⼿法を与えるものであり，本論⽂で⽬指す動的環境推定モデル と良く適合する概念である．データ同化の具体的⼿法には様々なものがあるが，本論⽂では逐次モンテ カルロフィルタ（SMCフィルタ）を採⽤した．これは，⾮線形性が強く，観測データがシミュレーショ ン上で確率的に得られるという微⼩デブリの衝突および環境推定に対して SMC フィルタの特徴が適し ているためである．また，SMCフィルタの適⽤のために必要なシステムモデルと観測モデルの定義を，

拘束⽅程式およびトーラスモデルを利⽤することで実現した．

続いて，構築された動的環境推定モデルの検証を⾏った．本論⽂の⼿法が軌道環境を⼗分に推定でき ることを確かめるため，既存の環境モデルであるMASTER-2009を真値とした観測と推定のシミュレー ションを⾏った．その結果，妥当な初期分布を与えることで観測データから微⼩デブリの分布を推定す ることが可能であることが⽰された．特に，設計上重要となる衝突⽅位⾓と衝突フラックスの関係の推 定が正しく⾏えていることで，本論⽂で構築した動的環境推定モデルが実際の宇宙機の安全に貢献でき ることが確認できた．また，複数回のシミュレーションを⾏うことで，本論⽂で提案する⼿法による軌 道環境の推定には最低半年間の観測が必要であることが明らかとなった．さらに，衝突個数の揺らぎに よる影響の評価も⾏い，衝突フラックスの過⼩評価のおそれが⼩さいことも⽰した．

以上のように，軌道上観測による微⼩デブリの環境推定という課題に対し，本論⽂では動的環境推定

モデルを新たに構築し，シミュレーションにより検証することでモデルの有⽤性と特性を明らかにした．

実際の軌道上観測と，本論⽂で構築された動的環境推定モデルを⽤いて微⼩デブリ環境をより正確かつ 迅速に把握することで，宇宙機の微⼩デブリに対する適切な防護設計や，危険が予測された際の回避運

⽤が可能となり，宇宙活動の安全および持続可能性への貢献が期待される．

⽬次

論⽂要旨 3

⽬次 5

図⽬次 7

表⽬次 9

第1章 序論 10

1.1. 背景 10

1. 1. 1 スペースデブリ問題 10

1. 1. 2 微⼩デブリ 13

1. 1. 3 先⾏研究 15

1.2. 論⽂の着眼点と構成 17

第2章 観測され得る微⼩デブリの軌道特性 18

2.1. 緒⾔ 18

2.2. 球形有限要素モデル 18

2.3. 軌道⾯の拘束⽅程式 26

2.4. 衝突フラックスの近似計算法 32

2.5. 結⾔ 37

第3章 動的環境推定モデル 38

3.1. 緒⾔ 38

3.2. データ同化 38

3.3. 逐次モンテカルロフィルタ 41

3.4. 微⼩デブリ環境のモデリング 48

3. 4. 1 システムモデル 49

3. 4. 2 観測モデル 51

3.5. 結⾔ 53

第4章 シミュレーションによる検証 54

4.1. 緒⾔ 54

4.2. MASTER-2009を⽤いた衝突シミュレーション 55

4.3. モデルの初期分布 57

4.4. 推定結果 60

4.5. 個数を基にした⽐較 66

4.7. 衝突個数の揺れによる影響の評価 72

4.8. 結⾔ 75

第5章 結論 76

謝辞 77

参考⽂献 78

付録A データ同化の基礎となる数学 80

付録B 軌道⾯の拘束⽅程式の結果 81

付録C 軌道⾯分布推定結果 87

付録D 複数機による観測 105

図⽬次

図 1.1 軌道上に存在し追跡されている物体の個数推移[1] 11

図 1.2 将来のデブリ個数予測[3] 11

図 1.3 デブリのサイズ・⾼度別観測状況 [11] 12 図 1.4 微⼩デブリの衝突実験により損傷したケーブル[12] 13 図 1.5 微⼩デブリの衝突によるISS太陽電池パネルの破壊[13] 13 図 1.6 推定衝突フラックスの⽐較[22]（左：ISS軌道，右：太陽同期軌道） 15 図 1.7 破砕由来のデブリ軌道⾯の拘束⽅程式[26] 16

図 2.1 球形有限要素モデル 19

図 2.2 検査要素内パラメータ定義[24] 19

図 2.3 衝突⾓の定義 22

図 2.4 観測衛星の⾶⾏姿勢（⻩⾊がセンサ⾯を⽰す） 23

図 2.5 遠地点近地点フィルタ[27] 23

図 2.6 衝突フラックスとして寄与した軌道：𝑢 = 65° 25 図 2.7 衝突フラックスとして寄与した軌道：𝑢 = 94° 25

図 2.10 座標系の定義 26

図 2.11 観測衛星とデブリの軌道⾯の関係 27

図 2.10 衝突フラックスとして寄与した軌道：𝑢 = 65° 30 図 2.13 衝突フラックスとして寄与した軌道：𝑢 = 94° 30 図 2.18 フラックスとして寄与した軌道の拘束⽅程式の値 31 図 2.19 検査要素の概念図（左：球形有限要素モデル，右：トーラスモデル） 32 図 2.20 トーラスモデルにおける検査要素の定義 33

図 2.21 Δ𝑓と𝜙，𝜃の関係 34

図 2.22 デブリの速度，観測衛星の速度，相対速度の関係 34 図 2.23 衝突フラックスの近似計算結果：𝑢 = 60° 35 図 2.24 衝突フラックスの近似計算結果：𝑢 = 90° 36 図 2.25 衝突フラックスの近似計算結果：𝑢 = 120° 36

図 3.1 動的環境推定モデルの概念図 38

図 3.2 確率分布の近似の概念図 42

図 3.3 リサンプリングの概念図 46

図 3.4 SMCフィルタのアルゴリズム 47

図 3.5 MASTER-2009によるデブリの離⼼率と衝突フラックスの関係 48

図 3.6 デブリ分布の離散集合近似の模式図 50

図 4.1 動的環境推定モデルのMASTER-2009を⽤いた検証の流れ 54

図 4.2 シミュレートされた観測データ 56

図 4.3 シミュレートされた観測データ（別例） 56

図 4.4 被追跡物体の軌道⾯の分布 57

図 4.5 初期分布 Case A 58

図 4.6 初期分布 Case B 59

図 4.7 推定結果（Case A, Day 5） 61 図 4.8 推定結果（Case B, Day 5） 61 図 4.9 推定結果（Case A, Day 180） 62 図 4.10 推定結果（Case A, Day 365） 62 図 4.11 推定結果（Case A, Day 540） 63 図 4.12 推定結果（Case A, Day 735） 63 図 4.13 推定結果（Case B, Day 180） 64 図 4.14 推定結果（Case B, Day 365） 64 図 4.15 推定結果（Case B, Day 540） 65 図 4.16 推定結果（Case B, Day 735） 65 図 4.17 軌道傾斜⾓分布の⽐較（Case A） 66 図 4.18 軌道傾斜⾓分布の⽐較（Case B） 67

図 4.19 軌道傾斜⾓の累積分布割合 68

図 4.20 推定総個数の時間変化 68

図 4.21 衝突⽅位⾓に対するフラックス分布 Case A 70 図 4.22 衝突⽅位⾓に対するフラックス分布 Case B 70 図 4.23 軌道傾斜⾓に対するフラックス分布 Case A 71 図 4.24 軌道傾斜⾓に対するフラックス分布 Case B 71 図 4.25 50例の衝突データを⽤いた総個数の推定結果 72 図 4.26 50例の推定結果の平均および標準偏差 73 図 4.27 観測されたデブリの個数と推定衝突フラックスの関係 74

図B.1-36 衝突フラックスとして寄与した軌道 81-86

図C.1-50デブリの軌道⾯分布推定結果 Case A 87-95

図C.51-100デブリの軌道⾯分布推定結果 Case B 96-104

図D.1 観測データ Constellation A 106

図D.2 観測データ Constellation B 107

図D.3 衝突⽅位⾓に対するフラックス分布 Constellation A 107

図D.4 衝突⽅位⾓に対するフラックス分布 Constellation B 108

図D.5 軌道傾斜⾓に対するフラックス分布 Constellation A 108

図D.6 軌道傾斜⾓に対するフラックス分布 Constellation B 109

図D.7 推定総個数の時間変化 109

表⽬次

表 2.1 セルの分割条件 19

表 4.1 計算機および乱数発⽣アルゴリズム 54

表 4.2 シミュレーション条件 55

表 4.3 ノイズの標準偏差 60

表 D.1 ２機⽬の観測衛星の軌道⾯ 105

第1章 序論

1. 1.

背景

1. 1. 1 スペースデブリ問題

1957 年に史上初の⼈⼯衛星スプートニク 1 号が打ち上げられて以来，⼈類は継続して宇宙開発・宇 宙利⽤を⾏ってきた．そのような⼈類の宇宙活動に伴って発⽣し軌道上に残留している不要な⼈⼯物は スペースデブリ，または宇宙ゴミと呼称されている．地球周回軌道上の物体は低軌道（⾼度2000km以 下）では秒速7km以上，静⽌軌道（⾼度約36000km）でも秒速3kmで周回しているため，宇宙機に対 してスペースデブリが衝突した場合には致命的な損害を与え得る．⼈⼯衛星は地球観測や通信等，現代 社会に不可⽋なインフラとなっており，それらを脅かしかねないスペースデブリは，宇宙開発を持続可 能なものとするために必ず解決されなければならない課題である．

図 1.1 は⽶国航空宇宙局により発表された地球周囲の軌道が追跡されている物体の個数推移[1]を⽰

したものであり，⻘線で⽰されているSpacecraftの個数が軌道上の物体数に占める割合は３分の1にも 満たないことがわかる．さらに図中のSpacecraftはすでに運⽤を終了しデブリとなった宇宙機を多数含 んでいるため，地球周辺に存在する⼈⼯物の殆どはデブリと⾔える．また，図 1.1はデブリを発⽣原因 ごとに分類しており，個数の多い順にFragmentation, Mission-related, Rocket Bodyとなっている．それぞ れ，破砕によって⽣じた破⽚，投棄された観測機器の蓋などミッションに伴い⽣じたデブリ，打ち上げ 後軌道上に留まっているロケットの上段機体を意味する．

図 1.1中，スペースデブリの個数が急増している部分が2箇所ある．1つは2007年に⾏われた中国 衛星破壊兵器実験，もう1つは2009年に発⽣した⽶露通信衛星衝突事故である．

2007年1⽉11⽇のFengyun1C衛星破壊兵器実験は，中国によって⾏われた弾道ミサイルによる同国

の気象衛星Fengyun1Cを標的とした衛星破壊兵器の実験である．この実験によって発⽣したデブリは，

観測されているもので3,000個以上あり，ESA の解析ツールであるPOEM (Program for Orbital Debris Environment Modeling)によれば，1 cm以上のデブリがおよそ60,000個，1 mm以上では，3,000,000個で あろうと予測されている．

2009年2⽉10⽇には⽶国通信衛星Iridium33とロシア通信衛星Cosmos2251が衝突する事故が発⽣し

た．Cosmos2251は1995年に運⽤を終了していたが，Iridium33はまだ運⽤中であった．Fengyun1Cのよ

うな軍事⽬的の意図的な衛星破壊兵器実験や衛星対破⽚の衝突は過去に数回発⽣した事例があるが，本 例は史上初の偶発的な衛星対衛星の⼤規模な衝突である．この事故によって発⽣したデブリは，観測さ れているもので1,000個以上あり， POEMによれば，1 cm以上のデブリがおよそ80,000個発⽣したと 予測されている．

これらは，意図的か偶発的かに関わらず，破砕事象が軌道環境に⼤きな影響をもたらすことを⽰して いる．このような破砕事象の可能性はデブリ数の増加に伴い増⼤しており，その帰結として現在，ケス

11

図 1.1 軌道上に存在し追跡されている物体の個数推移[1]

図 1.2 将来のデブリ個数予測[3]

www.orbitaldebris.jsc.nasa.gov

removal rate of more than five objects per year must be implemented.

The predicted collision activities in LEO from the three scenarios are shown in Figure 4. The top three curves indicate the total numbers of “all” collisions – catastrophic and non-catastrophic collisions. The bottom three curves show just the catastrophic collisions. A catastrophic collision occurs when the ratio of impact energy to target mass exceeds 40 J/g.

The outcome of a catastrophic collision is the total fragmentation of the target, whereas a non-catastrophic collision only results in minor damage to the target and generates a small amount of debris that has minimal contribution to population growth. Table 1 provides additional details of the catastrophic collisions predicted by the three scenarios.

Collisions are separated into three categories – those involving two intact objects (i-i), those involving one intact and one fragment (i-f), and those involving two fragments (f-f). Each number is the average of 100 MC runs. In general, an intact-intact collision will generate more debris than an intact-fragment collision.

The ADR target selection criterion used in the LEGEND simulations was the [mass × collision probability] value of each object. This criterion can be applied to objects in the current environment to identify potential targets for removal in the near future. The altitude-versus- inclination distribution of the top 500 objects identified via this selection criterion is shown in Figure 5. The prograde group is dominated by several well-known classes of vehicles: SL-3 R/Bs (Vostok second stages; 2.6 m diameter by 3.8 m length; 1440 kg dry mass), SL-8 R/Bs (Kosmos 3M second stages; 2.4 m diameter by 6 m length; 1400 kg dry mass), SL-16 R/Bs (Zenit second stages, 4 m diameter by 12 m length; 8900 kg dry mass), and various Meteor- series and Cosmos S/Cs (masses ranging from 1300 to 2800 kg). The fact that many of the SL R/Bs are high on the list is not surprising, based on the mass distribution in the environment (see also Figures 1 and 2). Below 1100 km altitude, the total mass of all SL-3, SL-8, and SL-16 R/Bs is about 500 tons, which accounts for about 20% of the total mass in LEO. Objects in the retrograde region are more diverse. They include, for example, Ariane R/Bs (1700 kg dry mass), CZ-series R/Bs (1700–3400 kg dry mass), H-2 R/Bs (3000 kg dry mass), SL-16 R/Bs, and S/Cs such as Envisat (8000 kg) and

continued from page 4

ADR

Figure 3. The benefits of using ADR to better limit the growth of the future LEO population. Each future projection is the average of 100 LEGEND Monte Carlo (MC) runs. An ADR of five objects per year can stabilize the future environment in LEO.

20 30 40 50

rooisionsobtsLEO

LEO En iron nt ro tion

Reg Launches + 90% PMD (all collisions)

Reg Launches + 90% PMD + ADR2020/02 (all collisions) Reg Launches + 90% PMD + ADR2020/05 (all collisions) Reg Launches + 90% PMD (cat. collisions)

Reg Launches + 90% PMD + ADR2020/02 (cat. collisions) Reg Launches + 90% PMD + ADR2020/05 (cat. collisions)

0 10 0

2010 2030 2050 2070 2090 2110 2130 2150 2170 2190 2210

uuatiubr

ar

Figure 4. Projected LEO collision activities in the next 200 years. Even with an ADR of 5 objects per year, a total of 14 catastrophic collisions are expected in the next 200 years.

10000 12000 14000 16000 18000 20000 22000 24000

ubroObts

LEO En iron nt ro tion a ra s o LE E D M runs

Reg Launches + 90% PMD

Reg Launches + 90% PMD + ADR2020/02 Reg Launches + 90% PMD + ADR2020/05

0 2000 4000 6000 8000

1950 1970 1990 2010 2030 2050 2070 2090 2110 2130 2150 2170 2190 2210

Etiu

ar

Table 1 – Predicted LEO catastrophic collisions in the next 200 years.

i-i collisions i-f collisions f-f collisions Total

90% PMD 10.2 10.9 3.0 24.1

90% PMD+ADR2020/2 8.2 7.0 1.9 17.1

90% PMD+ADR2020/5 6.5 5.5 1.8 13.8

ケスラーシンドロームとは，デブリ同⼠やデブリと宇宙機の衝突によって多数の新たなデブリが⽣ま れることによってデブリの数が⾃⼰増殖を続け，新たな宇宙開発・利⽤が不可能となるという想定であ り，D. Kessler（1978，[2]）によって提唱されたためこのように呼ばれている．

図 1.2は，運⽤を終了した衛星の90％が軌道離脱（Post Mission Disposal, PMD）を⾏い，さらに能動 的除去（Active Debris Removal, ADR）を⾏わなかった場合，2020年以降年間2個除去した場合，年間5 個除去した場合を想定し，将来の軌道上物体数をシミュレーションにより予測したものである．この結 果は，ケスラーシンドロームが既に現実的な想定となっており，これを防ぐためには衛星の正常な運⽤

終了と軌道離脱を確実なものとしながら⼤きなデブリの除去を⾏っていかなくてはならないことを⽰

している．

スペースデブリ対策を⾏うために重要なのは，まず現在の環境を正確に把握することである．スペー スデブリ環境の将来予測を⾏う環境モデルとして⽶国航空宇宙局（National Aeronautics and Space Administration, NASA）によるLEGEND[4]，欧州宇宙機関（European Space Agency, ESA）によるDELTA[5]，

⽇本の九州⼤学と宇宙航空研究開発機構によるNEODEEM[6]等が存在するが，それらの予測の元となる 現在の環境は観測によって得られた情報に基づいている．

現在，およそ10 cm以上の物体は光学観測により軌道の追跡がなされており，⽶国戦略軍統合宇宙運

⽤センター（Joint Space Operations Center, JSpOC）でカタログ化されている[7]．また，⾼度2000 km以 下の低軌道（Low Earth Orbit, LEO）に存在する2 mm以上のデブリはヘイスタックレーダ等により検出 が可能である[8]．⼀⽅，それよりも更に⼩さいデブリについては，スペースシャトルへの衝突痕[9][10]

などによってわずかな情報が得られているのみである．図 1.3は，このような状況を図⽰したものであ る．

図 1.3 デブリのサイズ・⾼度別観測状況 [11]

1. 1. 2 微⼩デブリ

前⼩節に⽰したように，2 mm 以下のデブリを地上から観測する⼿段は存在しない．このようなデブ リを，微⼩デブリと呼ぶ．これらは微⼩であっても⾼速であり⼤きな運動エネルギーをもつため，宇宙 機に対して衝突した場合には致命的な損傷を与え得る．図 1.4は，ケーブルハーネスに対し微⼩デブリ

を模した0.3 mmのステンレス粒⼦を4.01 km/sで衝突させた実験の結果であり，ケーブルの損傷が発⽣

したことが⽰されている[12]．また，実際の宇宙空間においても微⼩デブリの衝突は発⽣しており，図 1.5に⽰すように国際宇宙ステーション（International Space Station, ISS）の太陽電池アレイが損傷した事 例が報告されている[13]．これらの事例が⽰すように，微⼩デブリの衝突は宇宙機の構成要素を破壊す るには⼗分なエネルギーを持っている．そのため，宇宙機の機能喪失を防ぐには電源系統など重要な要 素に対する微⼩デブリの衝突を予測し，防護壁や配置の⼯夫によりその確率を最⼩化するような設計が 求められる．

図 1.4 微⼩デブリの衝突実験により損傷したケーブル[12]

図 1.5 微⼩デブリの衝突によるISS太陽電池パネルの破壊[13]

宇宙機に対する微⼩デブリの衝突を予測することのできる環境モデルとしては，NASA による ORDEM（Orbital Debris Engineering Model）[14]，ESA による MASTER（Meteoroid and Space Debris Terrestrial Environment Reference）[15]が存在する．それぞれの最新バージョンは2018年現在ORDEM 3.0

とMASTER-2009となっている．これらのモデルの詳細は次の⼩節に述べるが，地上から観測できない

微⼩デブリに関しては，表⾯検査によるデータが使⽤されている．表⾯検査とは，軌道上から地上に帰 還した物体の表⾯を検査することにより軌道上で衝突した微⼩デブリを検出する⼿法である．これまで に，LDEF（Long Duration Exposure Facility），SFU（Space Flyer Unit），スペースシャトルの窓やラジエー タなどからデータが得られている．

しかしながら，表⾯検査にはいくつかの課題がある．まず，宇宙機が軌道上にあった期間のうちどの 時刻で衝突が発⽣したのかが不明であり，データの時間分解能が低い．また，機体を地上に帰還させな ければならないためデータを得られる軌道が限られており，⾼度650 km以上については観測できてい ない．更に，スペースシャトルが既に退役しているため，2018年現在では宇宙空間から⼤きな物体を地 上へ帰還させる⼿段が存在しない．

⼀⽅で，デブリの環境は常に変動している．微⼩デブリの環境は宇宙空間での破砕事象を反映したも のであり，ケスラーシンドロームの懸念に対する状況認識や⼤規模破砕による劇的環境変動の迅速な検 知のため，微⼩デブリの常時観測が重要であると考えられる．

以上に述べた表⾯検査の時空間分解能の低さと常時観測の必要性から，微⼩デブリの軌道上観測が近 年進められており，世界での例としては ESA の GORID[16]，英国の DEBIE[17][18]，フランスの

SODAD[19]，NASAのSDS[20]などが挙げられる．九州⼤学においても，IDEA Projectとして2011年よ

り超⼩型衛星による微⼩デブリ観測計画を進めてきた．IDEA Projectは「デブリのその場環境認識」を 意味する“In-situ Debris Environmental Awareness”の頭字語から命名され，以下のような⽬標を持つ．

1. 微⼩デブリの軌道上観測

IDEA衛星はSpace Debris Monitor (SDM)[21]を搭載する．これにより，IDEA衛星は軌道上での

微⼩デブリの衝突時刻・位置を記録する．

2. 動的環境モデルの構築

衛星による軌道上観測の利点はデータを観測後直ちに得られることである．そこで，IDEA

Projectでは常に最新の環境に更新される動的な環境モデルを構築する．

3. 環境変動の検知

破砕事象により発⽣した微⼩デブリを検出し，環境変動を迅速に検知する．

IDEA Project の衛星としては，九州⼤学でフィージビリティスタディとして研究・開発を進めた

“IDEA-1”を元に，計画の意義・理念に共感しオーエスジー株式会社の⽀援を得た株式会社アストロスケ ールが“IDEA OSG 1”を開発・製造した．IDEA OSG 1は2017年11⽉28⽇にロシアより打ち上げられた が，ロケットのトラブルにより⼤⻄洋上に落下し失われた．なお，IDEA-1の軌道は⾼度798 kmの太陽

1. 1. 3 先⾏研究

本⼩節では，既存の環境モデルであるMASTERおよびORDEMによる微⼩デブリ環境の推定および

IDEA Projectに関しての九州⼤学における先⾏研究について述べる．

既存の環境モデル

ESAのMASTERとNASAのORDEMはいずれも，これまでの観測・実験から得られたデータを元に

構築されたスペースデブリの環境モデルであり，微⼩デブリを含むデブリの宇宙機に対する衝突を衝突 フラックスとして予測する機能を持つ．２つのモデルの間には⼿法および結果に違いがあり，Kriskoら

（2015，[22]）による⽐較が⾏われている．

⼿法の違いとしては，ORDEMはデブリを密度により分類しているのに対し，MASTERでは発⽣源に より分類しているという点が挙げられる．ORDEMが密度に着⽬するのは，開発当初の⽬的が有⼈宇宙 機の安全解析であり，衝突時の危険性が密度に⼤きく影響されるためである．⼀⽅，MASTERはこれま でに発⽣した爆発，衝突等のイベント履歴に基づいたモデリングを⾏っているため，デブリはそれに合 わせて発⽣源ごとに分類されている．なお，どちらのモデルも，モデリングしたデブリ環境を元にして 観測データと⽐較し，擦り合わせを⾏っているという点では共通している．

図 1.6はMASTERとORDEMにより推定される衝突フラックスのサイズ別分布を，国際宇宙ステー ションの軌道および⾼度約850 kmの太陽同期軌道で⽐較したものである．モデル間の差はISS 軌道よ りも太陽同期軌道で⼤きくなっており，特に2 mm以下の微⼩デブリでは最⼤で10倍以上の顕著な差が 現れている．この差は主に両モデルの解析⼿法の差に起因していると考えられるが，⼤型の物体や ISS 軌道のように実観測データが得られている場合には⽐較的⼀致しており，実観測データが正確な環境推 定のために重要であることは明らかである．特に，⾼度800 km付近では⼤規模な破砕事象である中国 衛星破壊兵器実験と⽶露通信衛星衝突事故が発⽣しており，地球観測衛星の需要が⾼い軌道であるため，

衛星による微⼩デブリの軌道上観測とそれによる正確な環境把握が不可⽋であると⾔える．

図 1.6 推定衝突フラックスの⽐較[22]（左：ISS軌道，右：太陽同期軌道）

九州⼤学における先⾏研究

1.1.2節で挙げたIDEA Projectの⽬的である「動的環境モデルの構築」および「環境変動の検知」に関

して，これまでに⾏われた研究について述べる．

環境変動の検知については，⼟井（2013，[23]），⽥﨑（2014，[24]），藤⽥（2016，[25]），兒⽟（2017，

[26]）らにより研究が進められた．まず，⼟井の中国衛星破壊兵器実験をテストケースとするシミュレ ーションにより，破砕を起源とする微⼩デブリは主として観測衛星と破砕親物体の軌道⾯の交線上で検 出されることが⽰された．次に，⽥﨑および藤⽥によって，2機の観測衛星と最⼩⼆乗法を⽤いること で破砕軌道⾯を推定する⼿法が提案された．また藤⽥らは，本論⽂で導⼊する軌道⾯の拘束⽅程式を⽤

いることで，破砕軌道⾯の推定において正しい昇交点移動率を⾒出すことの重要性を明らかにした．こ れは，正しい昇交点移動率を与えて衝突データそれぞれに対応する拘束⽅程式（本論⽂第２章に詳述）

を描いたとき図 1.7左のように⼀本の曲線となることから⾒出された．更に兒⽟らは，図 1.7左のよう になるのは観測衛星と破砕軌道⾯の昇交点移動率が近い場合のみであり，昇交点移動率の差が⼤きい場 合には拘束⽅程式が図 1.7右のようにただ２つの交点を持つことを明らかにした．そして，昇交点移動 率に関する⽅程式を導出することで，昇交点移動率の差が⼤きい場合には衛星１機のみによる観測で破 砕軌道⾯の推定を可能とし，IDEA衛星による環境変動の検知が現実的であることを⽰した．

図 1.7 破砕由来のデブリ軌道⾯の拘束⽅程式[26]

（左：昇交点移動率が近い場合，右：昇交点移動率の差が⼤きい場合）

⼀⽅で，観測データを⽤いた環境全体の推定については先⾏研究では完了していない．この⽬的に関し て阿江ら（2013，[27]，2014，[28]）は，遠地点近地点フィルタと球形有限要素モデルを⽤いて，被追跡 物体の軌道を対象に解析を⾏った．その結果，ミッション期間全体を通して衝突フラックスとして寄与 する軌道の特徴として，円軌道が⽀配的であることが⽰されている．しかし，阿江らの解析はミッショ ン期間全体を対象としたものであり，観測衛星によりリアルタイムに取得されるデータを有効に活⽤す るためには，衝突の瞬間に着⽬し，その瞬間に衝突フラックスとして寄与する軌道の解析を⾏う必要が ある．

1. 2.

論⽂の着眼点と構成

以上のような背景から，軌道上で観測された微⼩デブリのデータに基づき環境を推定する⼿法が求め られている．衛星を⽤いた微⼩デブリの軌道上観測の利点はデータをリアルタイムに得られることであ るため，本論⽂ではこれに着⽬し，データを得られる度に更新できる動的な環境推定⼿法の構築を⽬的 とする．環境推定を逐次更新することで常に最新かつ正確な環境を予測できれば，宇宙機のデブリに対 する防護設計をより適切に⾏うことが可能となる．さらに，すでに軌道上にある宇宙機に関しても，環 境の変動により特定の⽅向から⾶来するデブリが増加すると予測されればその⽅向に対する投影⾯積 を⼩さくする姿勢を取るなど，動的な環境推定によって微⼩デブリに対する対策を⾏うことができる．

このように，本論⽂で⽬的とする動的な環境推定⼿法の構築により，宇宙開発・宇宙利⽤の安全および 持続可能性への貢献が期待される．

以下に本論⽂の構成を⽰す．

第１章 序論

第１章（本章）では，本論⽂の背景，⽬的，構成について述べる．

第２章 観測され得るデブリの軌道特性

第２章では，軌道上観測により検出され得るデブリの軌道を解析する．まず，軌道⾯の拘束⽅程式に より，観測衛星とデブリの衝突可能性の有無を容易に判断できることを述べる．次に，トーラスモデル を導⼊し，衝突確率を少ない計算量で近似的に計算する．さらに，拘束⽅程式とトーラスモデルによる 計算結果を既存の⼿法によるシミュレーション結果と⽐較することで検証する．

第３章 動的環境推定モデル

第３章では，本論⽂で⽬的とする微⼩デブリの軌道上観測による環境推定⼿法を動的環境推定モデル として定義し，そのための数学的理論としてデータ同化と逐次モンテカルロフィルタについて述べる．

そして，第２章で導⼊したデブリの軌道特性に関する⼿法をフィルタに適⽤し，動的環境推定モデルを 理論的に定義する．

第４章 シミュレーションによる検証

第４章では，既存のモデルを真値としたシミュレーションを⾏うことで，第３章で定義した動的環境 推定モデルを検証する．シミュレートされた観測データから動的環境推定モデルを⽤いて推定されたデ ブリ分布を，デブリの個数および衝突フラックスについて真値とした既存のモデルと⽐較することによ り本論⽂で構築した環境推定の⼿法の妥当性を検証する．

第５章 結論

第５章では，本論⽂で得られた成果についてまとめる．

第2章 観測され得る微⼩デブリの軌道特性

2. 1.

緒⾔

本論⽂では，前章で述べたように，IDEA Projectにおける軌道上観測を想定している．具体的には，微

⼩デブリが衝突した時刻および位置を記録し，本研究はそのデータをもとに軌道上のデブリ分布を推定 することを⽬標とする．そこで，従来の表⾯検査ではわからなかった正確な衝突時刻・位置を得られる

というIDEA Projectの特⻑を最⼤限活⽤するために，本章では軌道上の観測衛星の位置とその瞬間に衝

突フラックスとして寄与する軌道の関係を解析する．

まず2.2節では，先⾏研究と同様に球形有限要素モデルを⽤いた解析を⾏う．球形有限要素モデルを

⽤いて抽出された，観測衛星に衝突フラックスとして寄与する軌道を観測衛星の位置ごとに解析する，

続いて2.3節では，観測衛星に衝突フラックスとして寄与する軌道⾯の分布を説明する拘束⽅程式を導 出する．この拘束⽅程式により，デブリが観測衛星に衝突し得るかどうかを容易に決定することが可能 となる．また，2.4 節では，数値シミュレーションによらず微⼩デブリと観測衛星の衝突フラックスを 求める近似的な⼿法であるトーラスモデルを提案する．このトーラスモデルにより，拘束⽅程式で求め た軌道⾯にあるデブリの観測衛星に対する衝突フラックスを近似的に求めることが可能となる．

新たに提案する⼿法はいずれも従来の⼿法である球形有限要素モデルと⽐較して効率的なものであ り，本研究で構築する数理モデルの基礎となるものである．

2. 2.

球形有限要素モデル

観測衛星とデブリそれぞれの軌道から衝突の可能性をシミュレーションにより求めるため，本研究で は球形有限要素モデルによる衝突フラックス法を⽤いた．

球形有限要素モデルとは ORDEM，MASTER および先⾏研究[24][27]でも⽤いられている⼿法で，図

2.1，図 2.2に⽰すように宇宙空間を⾚経，⾚緯，地球中⼼からの距離によって検査体積（セル）ごとに

区切ったモデルのことである．この⼿法では，観測衛星と物体の軌道が共通して通過するセルがあると き，そのセル内での衝突確率を表す衝突フラックス（collision flux）を計算する．本研究では，観測衛星 に対しある物体が衝突フラックスとして寄与するとき，その物体の軌道を検出され得る軌道として取り 扱った．

19

図 2.1 球形有限要素モデル

図 2.2 検査要素内パラメータ定義[24]

表 2.1 セルの分割条件

Δa[deg] Δd[deg] Δr[km]

1 1 20

Cell

Center of the Earth

Earth equatorial plane Geocentric right ascension

Declination

29 29

29

t t+dt N

ODRbin

p

lrel

v

( )

(

)

(

)

^{( )}

N ^{N ODRbin rel}

l lrel ODRbin

!!"

#

$$%

& ⋅ ⋅ ⋅

!! =

"

#

$$%

& ⋅

= +

−

∑

(1)

( ) ^t ( ^N ( ) ^{t N} ( ) ^{t dt} ) ^N ( ) ^{t dt}

N − −

= −

Right&ascension Declination

Radius

α_in α_out

δ_in

δ_out r_out

r_in

Δr Δα

Δδ

M_out v_out

M_in v_in

衝突フラックスの算出⽅法を以下に⽰す[24]．まず，図 2.2 の は⾚経， は⾚緯， は軌道半径，

添字の は検査要素に進⼊するとき， は出るときを表している．また，要素への進⼊時の速度と平 均近点離⾓をそれぞれ ， ，退出時を ， とする．本研究におけるセルの分割は表 2.1に⽰

す通りである．また，プライムを付したものを対象物体（デブリ）の，そうでないものを観測衛星のパ ラメータとする．

1つの要素内における衝突フラックス は観測衛星の存在確率 と物体の空間密度 と物 体の相対速度 を⽤いて，

(2.1)

と表せる．

ここで，観測衛星及び物体の存在確率 は，平均近点離⾓ を⽤いて，

(2.2)

で求められる．また，平均近点離⾓ は離⼼近点⾓ と離⼼率 から，

(2.3)

で求められる．さらに は と真近点離⾓ を⽤いて，式（2.4）で表される．

(2.4)

検査要素内への進⼊時，および退出時の を求める⽅法は，それぞれでどの⾯から進⼊・退出するか によって異なる．図2. 2中の軌道の進⼊時のように⾚経⽅向に垂直な⾯から進⼊・退出する場合は，進 退出の⾚経 ，昇交点⾚経 ，軌道傾斜⾓ ，近地点引数 を⽤いて，式（2.5）で求めることができ る．

(2.5)

α δ r

in

v

M

v

M

Φ

p

ρ

Δv

Φ_bin=p_SCbin

∑ ρ

p

p_bin =dM 2

π

M_in−M_out 2

π

M E

e

M = E − esin E

E

e

E=2 tan⁻¹ 1−e 1+etan f

2

"

#$ %

&

'

"

#$$ %

&

''

f

α Ω i ω

f =tan⁻¹ tan

( α

)

cosi

!

"

# $

%& −

ω

⽤いて球⾯三⾓法から，

(2.6)

で を求め，式（2.5）を⽤いて を求める．更に，軌道半径⽅向に垂直な⾯から進退出する場合は，

進退出時の ，軌道⻑半径 ， を⽤いて，式（2.7）で求められる．

(2.7)

以上より検査要素内の存在確率 を求めることが可能となる．物体の空間密度 は，検査要素の体 積 と存在確率 から，式（2.8）から算出される．

(2.8)

ここで，検査要素の体積 は，幾何学的に式（2.9）から求めることができる．

(2.9)

観測衛星と対象物体との相対速度は，観測衛星と物体の速度ベクトル ， を⽤いて，

(2.10)

で求めることができ，検査要素内を通過する際の平均速度ベクトルは，進⼊時の位置ベクトル ，退 出時の位置ベクトル ，検査要素内の滞在時間 を⽤いて，

(2.11)

としている．よって，式（2.10）と式（2.11）から観測衛星と対象物体との相対速度を算出することが できる．以上より，式（2.1）の右辺のパラメータをすべて求めることができ，検査要素内の衝突フラッ クスを計算できる．

さらに，この検査要素内の衝突フラックスを軌道⼀周回に渡って合計することで，観測衛星と対象物 体の⼀周回分の衝突フラックスが求まる．

α = sin

tan δ tan i

"

# $ %

&

' + Ω

α

r a e

f = cos

a ( 1− e

) ^{− r}

r ⋅ e

#

$

% %

&

' ( (

p

ρ

V

p

ρ

= p

V

V

V

= 2

3 3r

+ 1 4 ( ) Δr

!

"

# $

% & cos δ sin Δδ

2

!

"

# $

% &ΔαΔr v v '

x

x

Δt

Δv = x

− x

Δt

⾯に対する衝突⾓の考慮

また，以上の計算は対象の衛星を球体として考えた場合のものであり，平⾯に対する衝突フラックス の計算においては以下のような補正を⾏う．

まず，衛星に対するデブリの衝突⽅位⾓（Azimuth, Az）と衝突仰⾓（Elevation, El）を図 2.3のように 定義する．このとき，デブリの相対速度vrelは次のように表される．

𝒗₀₁₂= 𝑣₀₁₂4 − sin 𝐸𝑙

− cos 𝐸𝑙 cos 𝐴𝑧 cos 𝐸𝑙 sin 𝐴𝑧

? (2.12)

図 2.3 衝突⾓の定義

また，平⾯に対するデブリの⼊射⾓をθとすると，衝突フラックスは

Φ_ABC= 𝑝_E/GABC𝜌_ABC𝑣₀₁₂𝑐𝑜𝑠𝜃 (2.13)

と計算される．

ここで cosθは平⾯から機体中⼼に向かう単位法線ベクトルと，相対速度⽅向の単位ベクトルの内積

として計算できる．例として，図 2.4に⽰すような姿勢で観測を⾏うIDEA衛星の左右2枚のセンサ⾯

に対するフラックスを計算する．このときセンサ⾯の単位法線ベクトルは右⾯が(0, − 1 √2⁄ , 1 √2⁄ )，左

⾯が(0, − 1 √2⁄ , −1 √2⁄ )であり，相対速度⽅向の単位ベクトルは(−𝑠𝑖𝑛𝐸𝑙, −𝑐𝑜𝑠𝐸𝑙𝑐𝑜𝑠𝐴𝑧, 𝑐𝑜𝑠𝐸𝑙𝑠𝑖𝑛𝐴𝑧)であ る．したがって，式（2.13）のcosθは右⾯と左⾯でそれぞれ次のように計算される．

𝑐𝑜𝑠𝜃_U= 1

√2𝑐𝑜𝑠𝐸𝑙𝑐𝑜𝑠𝐴𝑧 − 1

√2𝑐𝑜𝑠𝐸𝑙𝑠𝑖𝑛𝐴𝑧

また，式（2.14）の値が負になる場合，センサ⾯の裏側，つまり衛星の観測⾯でない側からデブリが

⾶来することを表す．そのような場合には，衝突の可能性がある軌道であっても検出される可能性の無 い軌道といえる．これを，左右それぞれのセンサ⾯について計算し，判定することでセンサ⾯に対する 衝突⾓を考慮できる．

図 2.4 観測衛星の⾶⾏姿勢（⻩⾊がセンサ⾯を⽰す）

図 2.5 遠地点近地点フィルタ[27]

遠地点近地点フィルタ

さらに，計算の効率化のために本研究では図 2.5 に⽰す遠地点近地点フィルタ[27]を使⽤した．これ は，宇宙機の近地点がデブリの遠地点より外側である場合，または宇宙機の遠地点よりもデブリの近地 点が内側である場合には⾼度が交わることがなく，衝突しないことが⾃明であることを利⽤するもので ある．このフィルタにより衝突の可能性がある軌道のみを詳細に計算した．

衝突フラックスとして寄与するデブリの軌道⾯分布

以上の⼿法により，観測衛星に衝突フラックスとして寄与する物体の軌道の抽出を⾏う．観測衛星の 軌道は⾼度798 km，軌道傾斜⾓98.6°の太陽同期軌道とし，SSN (Space Surveillance Network) により追 跡されている⼤型物体の軌道をデブリの軌道として⽤いた．これは，微⼩デブリの軌道が不明であるた め，現実的に存在し得る軌道として追跡されている物体の軌道を参照したものである．

第１章で述べたように，阿江らの先⾏研究により観測衛星が検出するデブリはほとんどが円軌道にあ ることが報告されている．また，観測衛星とデブリがいずれも円軌道にあるならば，衝突したデブリと 観測衛星の⾼度は等しいと考えられる．そのため，衝突フラックスとして寄与するデブリ軌道の形状お よび⼤きさはほぼ同⼀であるとみなせるため，デブリの軌道は軌道⾯によって特徴付けられる．

そこで，シミュレーションにより抽出された衝突フラックスとして寄与するデブリの軌道⾯を，検出 時の観測衛星の緯度引数ごとに傾斜⾓ベクトル図上にプロットした．これにより，衝突の瞬間に着⽬し た解析が可能である．結果の⼀部を図 2.6, 7に⽰す．また，その他の緯度引数での結果は付録Bに⽰す．

これらの結果から，観測衛星に対して衝突フラックスとして寄与するデブリの軌道⾯は，傾斜⾓ベク トル上で何らかの曲線に沿って分布していることがわかる．また，その曲線は観測衛星の緯度引数によ って異なる形状を取っている．次節では，この軌道⾯の分布の原因について考察する．

図 2.6 衝突フラックスとして寄与した軌道：𝑢 = 65°

図 2.7 衝突フラックスとして寄与した軌道：𝑢 = 94°

2. 3.

軌道⾯の拘束⽅程式

本節では，観測衛星に対し衝突フラックスとして寄与する物体の軌道⾯について軌道⼒学に基づく考 察を⾏い，前節の結果を説明する軌道⾯に関する⽅程式を⾒出す．

座標系の定義

まず，以降の議論のため図 2.8に⽰す２つの直交座標系を定義する．

1つは地⼼⾚道⾯基準座標系であり，原点を地球中⼼，基準⾯を⾚道⾯，第⼀軸を春分点⽅向にとる．

このとき，第三軸は北極⽅向と⼀致する．この座標系を以後IJK座標系と呼ぶ．IJK座標系は慣性系で ある．

もう１つは地⼼軌道⾯基準座標系であり，原点を地球中⼼，基準⾯を衛星の軌道⾯，第⼀軸を地⼼か ら衛星に向かうベクトルと⼀致するようにとる．このとき，第三軸は衛星の⾓運動量ベクトル⽅向と⼀

致する．また，衛星が円軌道の場合には第⼆軸と衛星の速度ベクトルが平⾏となる．この座標系を以後 RSW座標系と呼ぶ．RSW座標系は対象とする衛星の位置・軌道によって変化する．

図 2.8 座標系の定義

I J

K

Orbit

Equator

(Vernal Equinox)

Satellite

R W

S

(Angular Momentum Vector)

Earth Center

拘束⽅程式の導出

観測衛星とデブリの軌道⾯の関係を図 2.9に⽰す．図は2物体の衝突が起こり得るのは軌道⾯の交線 上であり，衝突時には観測衛星とデブリ双⽅の地⼼からの位置ベクトルは同⼀となることを⽰している．

ここで，衝突位置の地⼼からのベクトルに平⾏な単位ベクトルを𝒆S，デブリ軌道の⾓速度ベクトルに 平⾏な単位ベクトルを𝒆_W^X とする．𝒆_W^X はその定義から常にデブリの位置ベクトルに直交するため，衝突 時には衝突位置の単位ベクトル𝒆Sとの間に式(2.15)が成⽴する．

𝒆_S∙ 𝒆_W^X = 0 (2.15)

図 2.9 観測衛星とデブリの軌道⾯の関係

また，単位ベクトル𝒆_W^X は IJK 座標系においてデブリの昇交点⾚経Ω^Xおよび軌道傾斜⾓𝑖^Xを⽤いて式

（2.16）で表されるため，式（2.15）から式（2.17）が求まる．

𝒆′_W= \ sin Ω^Xsin 𝑖^X

− cos Ω^Xsin 𝑖^X cos 𝑖^X

] (2.16)

𝑥

𝑟sin Ω^Xsin 𝑖^X−𝑦

𝑟cos Ω^Xsin 𝑖^X+𝑧

𝑟cos 𝑖^X= 0 (2.17)

ここで，(𝑥, 𝑦, 𝑧)は衝突位置のベクトルであり，𝑟はその⼤きさを表す．衝突位置ベクトルは観測衛星 によって記録される観測データであるから，式（2.17）は既知量である衝突位置によって未知量である デブリの昇交点⾚経および軌道傾斜⾓が拘束されることを意味する．そこで以後，式（2.17）を拘束⽅

程式と呼称する．

観測衛星によって記録される衝突位置は，複数の⽅法で表現され得る．まず，観測衛星の軌道要素で 表現した場合には，衝突位置の地⼼からのベクトルに平⾏な単位ベクトル𝒆Sは式（2.18）で表される．

𝒆_S = \cos Ω cos 𝑢 − sin Ω sin 𝑢 cos 𝑖 sin Ω cos 𝑢 + cos Ω sin 𝑢 cos 𝑖

sin 𝑢 sin 𝑖

] (2.18)

ここで，Ωは昇交点⾚経，𝑖は軌道傾斜⾓，𝑢は緯度引数であり，緯度引数𝑢は真近点離⾓𝑓および近点引 数ωの和として計算されるパラメータである．したがって，観測データが観測衛星の軌道要素として与 えられた場合には，拘束⽅程式は式（2.17），（2.18）から式（2.19）のように求まる．

sin(Ω^X− Ω) cos 𝑢 sin 𝑖^X− cos(Ω^X− Ω) cos 𝑖 sin 𝑢 sin 𝑖^X+ sin 𝑢 sin 𝑖 cos 𝑖^X= 0 (2.19)

また，衝突位置ベクトルを⾚経𝛼および⾚緯𝛿で表現することも可能であり，その場合𝒆_Sは式（2.20）

となる．

𝒆_S = \cos 𝛿 cos 𝛼 cos 𝛿 sin 𝛼

sin 𝛿

] (2.20)

このときの拘束⽅程式は式（2.21）である．

cos 𝛿 sin 𝑖^Xsin(𝛺^X− 𝛼) + sin 𝛿 cos 𝑖^X = 0 (2.21)

以上の式をデータフォーマットに応じて⽤いることで，観測された微⼩デブリが取り得る軌道を観測 データと拘束⽅程式によって絞り込むことが可能である．

軌道⾯間の⾓度

デブリの昇交点⾚経と軌道傾斜⾓が拘束⽅程式によって拘束されるということは，これらをただ1つ の未知変数によって表現できることを意味する．ここでは，デブリの昇交点⾚経および軌道傾斜⾓を２ つの軌道⾯がなす⾓（図 2.9中のθ）によって表現する．デブリの軌道⾯基準RSW系は観測衛星のRSW 系をθだけ回転させたものであるから，IJK系における𝒆_e^X は観測衛星のRSW系からIJK系への変換⾏

列Rを⽤いて次のように表すことができる．

[𝑅]i𝐶_k^lm \0 0 1] = 𝒆_e^X [𝑅] ≡ [𝐶_op^q ]i𝐶_oB^lm[𝐶_or^q ]

= \cos Ω cos 𝑢 − sin Ω sin 𝑢 cos 𝑖 − cos Ω sin 𝑢 − sin Ω cos 𝑢 cos 𝑖 sin Ω sin 𝑖 sin Ω cos 𝑢 + cos Ω sin 𝑢 cos 𝑖 − sin Ω sin 𝑢 + cos Ω cos 𝑢 cos 𝑖 − cos Ω sin 𝑖

sin 𝑢 sin 𝑖 cos 𝑢 sin 𝑖 cos 𝑖

]

これを変形することで，デブリの昇交点⾚経および軌道傾斜⾓と軌道⾯間の⾓度との関係を式（2.22），

（2.23）と表すことができる．なお，R12，R13……はそれぞれ⾏列Rの成分を表す．

tan 𝛺^X = 𝑅_lusin 𝜃 + 𝑅_lqcos 𝜃

−(𝑅_uusin 𝜃 + 𝑅_uqcos 𝜃) ^(2.22)

cos 𝑖^X= 𝑅_qusin 𝜃 + 𝑅_qqcos 𝜃 (2.23)

また，軌道⾯間の⾓度は⾓運動量ベクトルのなす⾓と⼀致するため，デブリの軌道⾯がわかっている場 合にはθは式（2.24）のように求まる．

cos 𝜃 = 𝒆_W∙ 𝒆′_W= cos 𝑖 cos 𝑖^′+ sin 𝑖 sin 𝑖^′cos(Ω^X− Ω) ^(2.24)

シミュレーション結果との⽐較

球形有限要素モデルによる結果との⽐較により，拘束⽅程式の妥当性を検証する．

図 2.10, 11は観測衛星がそれぞれに⽰す緯度引数のとき衝突フラックスとして寄与した物体の軌道を 傾斜⾓ベクトルで表現したものである．⻘で⽰された線は拘束⽅程式から導かれる検出し得る軌道⾯の 分布を，⾚で⽰された点は球形有限要素モデルで衝突フラックスとして寄与した物体の軌道⾯を表して いる．すべての図において，衝突フラックスとして寄与した物体の軌道⾯は拘束⽅程式の曲線上に存在 している．この結果から，本節で⾒出された軌道⾯の拘束⽅程式はシミュレーション結果をよく説明で きていることが確認できる．

図 2.10 衝突フラックスとして寄与した軌道：𝑢 = 65°

図 2.11 衝突フラックスとして寄与した軌道：𝑢 = 94°

続いて，拘束⽅程式からのずれについて考える．これは，球形有限要素モデルにおいては同じセルを 通過した場合にフラックスとして寄与すると考えるため，軌道が交わらない場合があるためである．

図 2.12 は，シミュレーション上において物体がフラックスとして寄与した場合の拘束⽅程式左辺の 値を横軸に，寄与した回数を縦軸にとりヒストグラムとして表したものである．

この結果から，衝突フラックスとして寄与した物体の軌道と観測衛星の軌道を⽤いて拘束⽅程式を計 算した場合，左辺の⼤きさはほぼ0.01以下となることがわかる．拘束⽅程式の左辺は衝突位置⽅向の単 位ベクトルと軌道⾯に垂直な単位ベクトルの内積であり，2ベクトルの間の⾓度の余弦である．

cos^-1(0.01)を計算するとおよそ89.4°となり，球形有限要素モデルのセルの幅が1°であることがシミュ

レーションにおいて拘束⽅程式を厳密に満たさない領域が現れた原因であると考えられる．

球形有限要素モデルで衝突し得る軌道に幅を持たせているのは軌道推定等に含まれる誤差を吸収し て考えるためであるから，拘束⽅程式を⽤いて衝突し得るデブリの軌道を絞り込む際にも同程度の幅を 持たせるべきであると考えられる．

図 2.12 フラックスとして寄与した軌道の拘束⽅程式の値

e

e

2. 4.

衝突フラックスの近似計算法

本節では，拘束⽅程式により限定される軌道⾯上にあるデブリの，観測衛星に対する衝突フラックス の計算について述べる．球形有限要素モデルでは，ある軌道にあるデブリの観測衛星に対する衝突フラ ックスを計算できる．しかし，球形有限要素モデルでは⼆物体の軌道をそれぞれセルに分割した上で各 セルにおける衝突フラックスを計算するため，⼤きな計算コストが必要となる．そこで本節では，球形 有限要素モデルよりも⼩さい計算コストで近似的に衝突フラックスを求める⼿法として，新たにトーラ スモデルを提案する．

図 2.13にトーラスモデルの基本的なアイディアを⽰す．球形有限要素モデルでは図 2.13左のように 衛星とデブリが共に通過する全てのセルで衝突フラックスを計算し合計しているため，検査要素は衝突 点近傍のセルの合計となる．そこで，これらのセルを図 2.13 右のように観測衛星とデブリの軌道を含 むような⼀つの検査要素として近似する⼿法がトーラスモデルである．

図 2.13 検査要素の概念図（左：球形有限要素モデル，右：トーラスモデル）

2.2節でも述べたように，検査要素内の衝突フラックス𝛷は式（2.25）で与えられる．

𝛷 =𝑝𝑝^X𝛥𝑣

𝑉_G122 ^(2.25)

ここで，𝑝は検査体積（セル）内の物体の存在確率，Δ𝒗はデブリと衛星の相対速度，𝑉G122は検査要素の 体積を表す．以下に，各要素の計算について詳述する．なお，デブリと衛星の軌道は同⼀⾼度の円軌道

検査要素の定義

図 2.14に⽰すように検査要素を定義する．まず，半径Rの円を𝒆_W+ 𝒆_W^X の周囲に回転させることで，

地球を中⼼とし2物体の軌道⾯の中間をとる円環体（トーラス）を考える．このトーラスの中⼼半径は 2物体の軌道と等しい 𝑟 とする．ここで，2物体の地⼼から⾒た離⾓が𝜙となるときの，衝突位置から の離⾓をΔ𝑓とし，検査要素を衝突位置からΔ𝑓までとして定義する．この検査体積は衝突位置に関して対 称に定義できるため，図 2.14のような検査要素を衝突位置の両側に定義することで図 2.13に⽰すよう な検査要素を定義できる．さらに，2物体が同⼀⾼度の円軌道にあるときは図に⽰す衝突位置から180°

対称な位置でも同じように衝突するため，図 2.14 に⽰す検査要素で求められた衝突フラックスの 4 倍 が軌道⼀周回あたりの衝突フラックスとなる．

図 2.14 トーラスモデルにおける検査要素の定義

各要素の導出

まず，図 2.15にΔ𝑓と𝜙，𝜃の関係を⽰す．左は𝜃 ≤ 𝜋 2⁄ ，右は𝜋 2⁄ < 𝜃の場合である．ここから，Δ𝑓 は球⾯三⾓法を⽤いて式（2.26）のように求まる．

sin Δ𝑓 sin𝜃

2= sin𝜙

2 𝜃 ≤𝜋 2 sin Δ𝑓 cos𝜃

2= sin𝜙 2

𝜋 2< 𝜃 }0 ≤ 𝛥𝑓 ≤𝜋

2 ~

(2.26) debris

φ θ

Δf r

R

Earth center

V

satellite

図 2.15 Δ𝑓と𝜙，𝜃の関係

続いて相対速度Δ𝒗を求める．トーラスモデルにおける2物体の検査要素内での相対速度は位置によっ て変化するが，ここでは衝突点における相対速度で代表させる．衝突点での相対速度Δ𝒗はデブリと観測 衛星がともに円軌道であるという仮定であれば衛星の速度𝒗とデブリの速度𝒗^X，相対速度Δ𝒗の関係は図 2.16のような⼆等辺三⾓形になるため，Δ𝒗は式（2.27）によって求まる．

∆𝑣 = 2𝑣 sin𝜃

2 ^(2.27)

図 2.16 デブリの速度，観測衛星の速度，相対速度の関係

1°

θ

1°

θ

Δf Δf Δf Δf

Δ v

v’

v

θ

次に，衛星が検査要素内に存在する確率𝑝は，円軌道であれば容易に式（2.28）のように求まる．

𝑝 =Δ𝑓

2𝜋 ^(2.28)

ただし，𝜃 < 𝜙 または 𝜃 > 𝜋 − 𝜙の場合には∆𝑓 = 𝜋 2⁄ となる．

最後に，検査要素の体積𝑉G122は式（2.29）により計算される．

𝑉_G122= (𝜋𝑅^u)(𝑓_•‚ƒ𝑟) (2.29)

以上の計算により，Rとφを定義すれば観測衛星およびデブリの軌道要素から衝突フラックスを近似 的に計算することが可能となった．

シミュレーション結果との⽐較

図 2.17-19に，トーラスモデルおよび球形有限要素モデルを⽤いて計算された衝突フラックスの⽐較 を⽰す．図はそれぞれ，観測衛星の緯度引数が 60°，90°，120°近傍のとき検出されたデブリの昇交 点⾚経と衝突フラックスの関係を⽰している．実線の⽰すトーラスモデルにより近似計算された衝突フ ラックスは，プロットの⽰すシミュレーション結果の傾向・オーダーを捉えることができていることが わかる．球形有限要素モデルにおけるセルサイズと近くなるよう，R = 60 km, φ=1°として計算した．

図 2.17 衝突フラックスの近似計算結果：𝑢 = 60°

図 2.18 衝突フラックスの近似計算結果：𝑢 = 90°

図 2.19 衝突フラックスの近似計算結果：𝑢 = 120°

2. 5.

結⾔

本章では，軌道上観測衛星により検出される微⼩デブリの軌道の性質について，拘束⽅程式とトーラ スモデルの２つの新たな解析⼿法を新たに提案した．拘束⽅程式は，観測されたデブリの昇交点⾚経お よび軌道傾斜⾓が取り得る値を衝突位置により拘束するものである．この⽅程式を⽤いることにより，

任意の位置における観測衛星に対し衝突し得るデブリの軌道を絞り込むことが可能となる．また，トー ラスモデルは2物体間の衝突フラックスを軌道要素と衝突位置を⽤いて近似的に計算する⼿法である．

本章ではまた，拘束⽅程式およびトーラスモデルについて球形有限要素モデルを⽤いたシミュレーシ ョンの結果との⽐較を⾏い，それぞれの⼿法がこれまでの解析と⽐較して妥当であることを確認した．

本章で新たに提案された⼿法は衝突可能性および衝突フラックスに関して数学的な記述を可能とす るため，環境推定のための数理モデルを構築する上で重要となる．また，球形有限要素モデルでは2物 体の軌道を多数の検査要素に分割し計算を⾏っていたのに対して，本章で新たに提案した⼿法はいずれ も軌道の分割を必要としないため計算量の削減に⼤きく貢献する．