Pulsed Plasma Thruster の超小型衛星における姿勢制御用アクチュエータとしての活用

(1)

平成

30

年度（

2018

年度）学位論文（修士）

Pulsed Plasma Thruster の超小型衛星における姿勢制御用アクチュエータとしての活用

Utilizing Pulsed Plasma Thruster as a attitude control actuator in microsatellite

首都大学東京大学院

システムデザイン研究科システムデザイン専攻航空宇宙システム工学域博士前期課程

学修番号 17891529 氏名森本明浩指導教員佐原宏典教授

平成 31年（2019年）1月25 日

(2)

(3)

摘要

本論文では，Pulsed Plasma Thruster（PPT）を超小型衛星の姿勢制御用アクチュエータとして用い，高精度姿勢制御手法の検討を姿勢制御シミュレータにより行う．PPTによる超小型衛星の制御の具体的な検証を行うために，本研究室で開発を行っている超小型衛星“ORBIS”のモデルを用いて姿勢制御手法の確立を目指す．

まず初めに人工衛星における一般的な制御系の制御フローを示し，本研究で用いたセンサのパラメータや人工衛星のモデルを示す．次に姿勢制御シミュレータの軌道上での衛星の姿勢，位置，速度の算出方法を説明し，軌道上での環境モデルを数学モデルから算出し，軌道上での姿勢の乱れの要因となる環境外乱の妥当性の検討を行った．その後，PPTでの制御アルゴリズムを提案し，そのアルゴリズムを姿勢制御シミュレータに実装，PPTにより高精度姿勢制御を達成し得ることを検証し，その制御精度の評価を行う．また，超小型衛星の電力の制約やミッションの運用期間を考慮し，高精度姿勢制御を達成し得るPPTのパラメータ要求を抽出する．さらにPPTの特性を踏まえ，超小型衛星とPPTの制約双方を満足し高精度姿勢制御を達成する運用方法を提案し，その評価を行った．

(4)

(5)

図目次

図 1.1 姿勢制御系システム構成図 ... 2

図 2.1人工衛星システム構成図⁷⁾ ... 3

図 2.2 姿勢決定センサ⁷⁾ ... 4

図 2.3 AxelStar-3外観⁹⁾ ... 5

図 2.4 CHR02-025外観¹⁰⁾ ... 5

図 2.5 カルマンフィルタ⁸⁾ ... 6

図 2.6姿勢制御アクチュエータ⁸⁾ ... 7

図 2.7 Pulsed Plasma Thrusterの推力発生メカニズム¹¹⁾ ... 8

図 2.8姿勢制御系ブロック線図 ... 9

図 3.1姿勢シミュレータ概略図 ... 10

図 3.2 ORBIS外観...11

図 3.3 座標系の定義 ... 12

図 3.4 環境外乱モデルフローチャート... 17

図 3.5波長10.7cm当たりの電波流速 ... 17

図 3.6 大気モデル計算結果 ... 18

図 3.7 高度500kmにおける大気密度⁸⁾ ... 18

図 3.8 地磁場計算結果（左：NOAA，右：本シミュレータ） ... 19

図 3.9 蝕模式図 ... 20

図 3.10 高度550kmにおける大気抵抗トルク ... 21

図 3.11 高度550kmにおける残留磁気トルク ... 22

図 3.12 高度550kmにおける重力傾斜トルク ... 23

図 3.13 高度550㎞における太陽輻射圧トルク ... 25

図 3.14 高度550kmにおける環境外乱トルク ... 25

図 3.15 500kg級衛星に作用する環境外乱トルク⁸⁾ ... 26

図 4.1 主要座標系間の関係¹⁸⁾ ... 28

図 4.2 STTによる姿勢決定精度（ピッチ/ヨー） ... 30

図 4.3 STTによる姿勢決定精度（ロール） ... 30

図 4.4 PPT出力例(𝑛𝑚𝑎𝑥 = 10) ... 34

図 4.5 PPT制御モデル ... 35

図 4.6 TMU-PPT-5P外観¹¹⁾ ... 35

図 4.7 環境外乱による姿勢変動（左：姿勢角，右：姿勢角速度） ... 36

図 4.8 PPT搭載位置 ... 36

図 4.9 PPTの1放電の出力 ... 36

図 4.10姿勢指向精度の時間分布と確率密度関数（𝑛𝑚𝑎𝑥 = 9，𝐼𝑏𝑖𝑡 = 48.2[μNs]） ... 38

図 4.11姿勢安定度の時間分布と確率密度関数𝑛𝑚𝑎𝑥 = 9，𝐼𝑏𝑖𝑡 = 48.2[μNs]） ... 38

図 5.1 軌道1周分の各サブシステムの消費電力 ... 40

図 5.2 最大放電回数𝑛𝑚𝑎𝑥と指向精度の関係 ... 41

図 5.3 最大放電回数𝑛𝑚𝑎𝑥と姿勢安定度の関係 ... 41

図 5.4姿勢指向精度の時間分布と確率密度関数（𝑛𝑚𝑎𝑥 = 4，𝐼𝑏𝑖𝑡 = 24.1[μNs]） ... 42

(8)

図 5.7姿勢安定度の時間分布と確率密度関数（𝑛𝑚𝑎𝑥 = 4，𝐼𝑏𝑖𝑡 = 32.13[μNs]，閾値0.00001） .... 46

(9)

表目次

表 1.1 人工衛星の姿勢制御精度の例1), 2), 3), 4),5) ... 1

表 2.1 AxelStar-3諸元⁹⁾ ... 5

表 2.2 CHR02-025諸元¹⁰⁾ ... 5

表 3.1 ORBSI諸元...11

表 3.2 超小型衛星ORBIS想定軌道 ... 15

表 3.3 大気モデル解析パラメータ ... 18

表 3.4 高度550kmにおける環境外乱トルクのオーダー... 26

表 4.1 最適制御パラメータ ... 33

表 4.2 TMU-PPT-5P諸元¹¹⁾ ... 35

表 4.3 PPTの発生するトルク... 37

表 4.4 スケーリングパラメータ，形状パラメータ（𝑛𝑚𝑎𝑥 = 9，𝐼𝑏𝑖𝑡 = 48.2[μNs]） ... 38

表 5.1 発電力計算パラメータ²¹⁾ ... 39

表 5.2 各サブシステムの使用機器と消費電力 ... 40

表 5.3 𝐼𝑏𝑖𝑡を変化させた場合の総放電回数，指向精度及び安定度（𝑛𝑚𝑎𝑥 = 4） ... 42

表 5.4 スケーリングパラメータ，形状パラメータ（𝑛𝑚𝑎𝑥 = 4，𝐼𝑏𝑖𝑡 = 24.1[μNs]） ... 43

表 5.5 PPT放電回数（𝑛𝑚𝑎𝑥 = 4，𝐼𝑏𝑖𝑡=24.1[𝜇𝑁𝑠]）... 43

表 5.6 高精度姿勢制御達成のためのPPTへのシステム要求 ... 44

表 5.7 閾値と総放電回数(𝑛𝑚𝑎𝑥 = 4) ... 45

表 5.8 各PPTの放電回数（𝑛𝑚𝑎𝑥 = 4，𝐼𝑏𝑖𝑡 = 32.13[𝜇𝑁𝑠]，閾値0.00001） ... 46

表 5.9 RWと併用し高精度姿勢制御を達成するためのシステム要求... 47

(10)

(11)

第 1 章序論

1. 1 背景

1950年代に世界初の人工衛星が打ち上げられて以降，地球観測衛星を始めとして様々な用途の人工衛星が打ち上げられ，生命の起源の解明，地上における我々の生活の向上などいろいろな場面で貢献している．そのような人工衛星であるが，限られた打ち上げ機会の有効活用のために機能充実化が行われ，

人工衛星の大型化が進む一方で2000年代に入り超小型衛星と呼ばれる質量 100kg 以下の人工衛星が大学をはじめとする各研究機関における研究が盛んにおこなわれるようになった．超小型衛星は，大型人工衛星と比べて低コストで短期開発が可能というメリットがある他，1 つのミッションへ専念できる専有性を持っている．その性質から新たな技術実証の場として用いられることに加えて，近年では天体観測や地球観測など理学観測の場としての利用も期待されている．理学観測衛星では多くの場合，観測対象を捉えるために人工衛星の姿勢を高精度に制御することが求められる．一方で，超小型衛星では電力の確保等の制約があるため，少ない電力で高精度な姿勢制御を達成することが要求される．表1.1にこれまでの大型衛星，小型衛星で達成，もしくは要求されている姿勢指向精度及び安定度を示す．

表 1.1 人工衛星の姿勢制御精度の例1), 2), 3), 4),5)

衛星名質量[kg] 開発機関打上年指向精度[deg] 姿勢安定度[arcsec/sec]

れいめい 60 JAXA 2005 0.05 14.4

ひので 900 JAXA 2006 - 0.002

きく8号 2，800 JAXA 2006 0.05 -

ほどよし3号 60 東京大学 2014 0.5 -

ひさき 35 JAXA，他 2013 0.025

Nano-JASMINE 35 東京大学，他未定 0.05 0.72

CIBER-2/EXZIT 複数 JAXA,他未定 0.1 0.06

宇宙赤外線背景放射 50 JAXA,他未定 0.1 0.1

ORBIS 50 首都大学東京，他未定 0.1 -

GEO-X 300 首都大学東京，他未定 0.01 1

HiZ-GUNDAM 400 JAXA 未定 0.02 10

通常，人工衛星の姿勢制御系は内界センサ，外界センサを用いて現在の姿勢を決定する姿勢決定系とアクチュエータを用いて所望する姿勢に制御する姿勢制御系に分かれている．各系の詳細な説明は次節で行う．現在の超小型衛星の姿勢制御系のシステムブロックの例を図1.1に示す．

一般的に超小型衛星において高精度な姿勢制御が求められる場合，姿勢制御アクチュエータとしてリアクションホイール（以下 RW）が使用されている．RW はきめ細やかな制御が可能であり，非常に高精度な制御が可能である．しかし，消費電力が大きく，比較的高価であること，稼働による角運動量の蓄積により飽和した際には他のアクチュエータによるアンローディングが必要となる等運用上の制約が多い．そこで本研究では姿勢制御アクチュエータに RW を用いるのではなく電気推進の一つである

Pulsed Plasma Thruster（以下PPT）を用いて超小型衛星の姿勢制御を行うことを提案する．PPTは低コス

トかつ低電力での運用が可能であり，任意の時間間隔で微小推力を発生させることが出来る．このPPT

(12)

型衛星でしか成しえなかったミッション，大型衛星では成しえなかったミッションが可能となり更なる宇宙利用の拡大と発展につながると期待される．

関連研究では2016年に興津によりPulsed Plasma Thrusterを用いた超小型衛星の高精度姿勢制御について研究が行われており，その際の結果ではPPTを用いて超小型衛星の姿勢制御を行うことにより指向精度についてはRWと同等程度，安定度についてはRWより劣るという結果が得られている．また，PPT へのシステム要求も抽出しており，その中の作動保証回数は1,141,050,973回の噴射が求められている⁶⁾．しかしながらPPTは動作原理の都合上約1000万回の噴射が上限であり，約11.5億回の噴射は不可能であるという結果が得られている．以上を踏まえPPTによる姿勢制御性能の向上と超小型衛星の制約だけでなくPPTの制約も踏まえた運用方法が求められている．

1. 2 目的

本論文では，Pulsed Plasma Thrusterの超小型衛星における姿勢制御用アクチュエータとしての活用として，以下を目的とする．

 PPTを姿勢制御用アクチュエータとして使用し，前研究で達成しえなかった姿勢安定度の向上と指向精度の更なる向上をすること．

 超小型衛星に搭載することを想定し，電力等の制約を踏まえたうえで高精度姿勢制御を達成するためのPPTのシステム要求を提案すること．

 超小型衛星の制約だけでなくPPTの制約も踏まえたうえで実現可能な運用方法を提案し，それを達成するためにPPTのシステム要求を提案すること．

なお，超小型衛星に搭載することを想定し衛星固有のパラメータを含む具体的な検証を行うために本論文における PPT の搭載対象は本研究室で開発を行っている超小型衛星“ORBIS”を前提としている．

ORBISに関する詳細は第2章及び第3章に記載する．

図 1.1 姿勢制御系システム構成図

(13)

第 2 章人工衛星の姿勢制御

本章では，人工衛星のシステム及び姿勢の制御方法についてを示す．

2. 1 人工衛星のシステム構成

人工衛星のシステム構成としては電力，通信，姿勢制御等，人工衛星の基本的な機能として必要なバス機器．そして通信，地球観測や科学データの収集などの人工衛星ごとのミッションのために必要なミッション機器の2つから構成されている．人工衛星のシステム構成の1例を図2.1に示す⁷⁾．

バス部には，様々な役割を担うサブシステムが存在し，例えばバッテリや太陽電池等を用いて電力収支を管理するサブシステムを電源系，通信機器を用いて地上とのコマンドの送受信を行うサブシステムを通信系と呼ぶ．姿勢制御を行うサブシステムは姿勢系と呼ばれ，姿勢系は姿勢決定系と姿勢制御系の 2

図 2.1人工衛星システム構成図⁷⁾

(14)

節に姿勢決定系で用いられるセンサ及びアルゴリズム，姿勢制御系で用いられるアクチュエータやそのアルゴルズムを紹介するとともに，本研究において用いたセンサやアクチュエータ，姿勢決定・制御アルゴリズムについて述べる．

2. 2 姿勢決定系

2. 2. 1 姿勢決定センサ

姿勢決定に用いられるセンサは外界センサと内界センサの2つの種類に分類される．外界センサは地球や恒星といった人工衛星外部の環境を測定し，基準方向に対する衛星の姿勢を検出するものであるが，

外界センサは蝕の間計測できない期間があり，観測データの取得には数秒から数十秒の時間間隔を要する．その間の姿勢決定を補完し，精度を高めるために内界センサを用いて衛星の角速度等のパラメータを計測し，時間経過による姿勢の変化を検出する．図2.2に人工衛星において用いられる姿勢決定センサの例を示す⁷⁾．

姿勢決定センサには数多くの種類が存在するが，本研究において姿勢角の検出は，高精度な姿勢決定が必要な場合に多く用いられている恒星センサ（STT）を用いる．角速度の検出には，理学観測などの高精度な姿勢検出が必要な衛星には，光ファイバジャイロ（FOG：Fiber Optical Gyro）が用いられることが多いが，1軸あたり100万円オーダーと機械式ジャイロに比べて高価であり，比較的消費電力も大きいという超小型衛星には適さない欠点があることに加えて，更に温度変化により計測誤差が生じるといった欠点もある．そこで本研究では，MEMS（Micro Electro Mechanical Systems）ジャイロを用いることを想定している．ここ数十年間において半導体微細加工技術の発達により，電子回路やセンサなどをシリコンチップ上に制作することが可能となってきた．この技術を MEMS と呼び，それを利用して製作されたジャイロをMEMSジャイロと呼ぶ．近年MEMS技術の発達や検出器の性能向上により，FOG に匹敵するMEMSジャイロが開発されるようになってきた．加えて3軸で10万円オーダーとFOGに

図 2.2 姿勢決定センサ⁷⁾

(15)

比べて非常に低価格であり，小型で温度制御が不要であるなど精度以外の面において MEMS ジャイロの優位性は明らかなものとなっている．

表2.1及び図2.3に本研究で用いたAxelSpace社製の恒星センサであるAxelStar-3の諸元と外観を，表 2.2及び図2.4にシリコンセンシング社製のMEMSジャイロであるCRH02-025の諸元と外観を示す．

表 2.1 AxelStar-3諸元⁹⁾ 精度

ｎ

ピッチ/ヨー 7[arcsec]（3σ）

ロール 77[arcsec]（3σ）

視野角 8[deg]×8[deg]

更新周期 1，0.5，0.2[Hz]

太陽除外角（ロングフード） 35[deg]

質量（ロングフード） 516[g]

寸法（ロングフード） 150×80×75[mm]

消費電力 3.0[W]（typ.）

動作温度 -20～50[℃]

表 2.2 CHR02-025諸元¹⁰⁾

ランダムノイズ 2.26e-3[rad/√sec]

ランダムウォークノイズ 5e-7[rad/(√sec²)]

質量 45[g]

寸法 33×33×25.4[mm]

消費電力 0.3[W]

s計測範囲 ±25[deg/sec]

2. 2. 2 姿勢決定アルゴリズム

高精度姿勢決定が要求される理学衛星では恒星センサとジャイロセンサを用いた姿勢決定が一般的であるが，STT及びジャイロの出力にはランダムノイズやバイアス，ドリフトといったノイズが含まれる．高精度姿勢決定を行うためにはセンサの出力からこれらのノイズを取り除き真の値を推定する必要がある．そこでよく用いられる姿勢決定アルゴリズムとして，カルマンフィルタがある．

カルマンフィルタは1960年代にR. E. Kalmanによって提唱された状態推定アルゴリズムである．時系列的な変化をするデータ履歴から次に取りうる値を予測するフィルタであり，姿勢伝搬と観測更新という2つの処理から構成される．姿勢伝搬は現在の姿勢とジャイロの角速度によって次の姿勢を予測する処理のことであるが，ジャイロが示す角速度にはランダムなノイズが含まれており，これが推定値の誤差として蓄積される．その誤差をSTTの観測値によって修正する処理が観測更新である．これらの姿勢伝搬と観測更新を繰返すことで姿勢推定制度を高めるものである．図2.5にカルマンフィルタによる処理のアルゴリズムを示す．

図 2.3 AxelStar-3外観⁹⁾

図 2.4 CHR02-025外観¹⁰⁾

(16)

2. 3 姿勢制御系

2. 3. 1 姿勢制御用アクチュエータ

姿勢制御に用いられるアクチュエータは，衛星本体と角運動量を交換することを通して，角運動量の時間変化率としてトルクを発生させる内力アクチュエータと，外力としてトルクを発生させる外力アクチュエータの2つに大別される．内力アクチュエータの例としては，衛星構体に固定されたホイールを有するモーメンタムホイールやリアクションホイール，ジンバル機構により角運動量ベクトルの発生方向を変化させることが出来るコントロールモーメンタムジャイロ（CMG）などが挙げられる．また，外力アクチュエータには燃料の放出によりトルクを発生させるスラスタや，衛星の周囲の場との相互作用を利用する磁気トルカやロッドなどがある．図2.6に代表的なアクチュエータの例を示す．

図 2.5 カルマンフィルタ⁸⁾

(17)

従来，高精度な姿勢制御が要求される衛星には，内力アクチュエータの一つであるRWが広く用いられてきた．しかし1.1 節で述べたようにRWは消費電力が大きく，比較的高価であること，角運動量が飽和した際に他のアクチュエータによるアンローディングが必要といった運用上の制約もある．そこで本研究ではRWではなく外力アクチュエータであり，電気推進の一つであるPulsed Plasma Thruster（PPT）

を用いて人工衛星の姿勢制御を行い，その際の制御性能の向上を目的の一つとしている．図2.7にPPT の推力発生メカニズムを示すとともに以下に説明する．

PPTは主に，2枚の電極（アノード及びカソード），固体プロペラント，プロペラント用のフィード機構，点火装置としてのイグナイタ，主放電エネルギー充電用のコンデンサ，イグナイタ用パルス電源，

主放電電源のみで構成されており，構造が非常に簡易であるという特徴がある．

1) イグナイタにより印加された高電圧の火花放電が少量の固体，または液体のプロペラントを昇華させ，その一部をプラズマ化させる（図2.7（a））．

2) 一部のプロペラントの昇華により生じたプラズマは，アノード・カソードの両電極間に広がり，

高導電性の領域を作る．この領域により両電極間が短絡され，両電極間につながれたキャパシタ内の電荷が一斉に流れ，主放電が形成される（図2.7（b））．

3) この手法電による電流成分はジュール加熱によってプロペラントにエネルギーを与え昇華させる．消化したプロペラントは，高温気体の膨張による空気力学的加速を受け，一部は電離してプラズマとなり，主放電電流とその自己誘起磁場が作るローレンツ力による電磁的加速を受ける

（図2.7（c））

4) 電磁気的・空気力学的加速を受けたプラズマは，下流方向に加速され，その放電領域を広げつつ，

スラスタ外へ排出される．

5) 消費量と同量のプロペラントがフィード機構によって自動的に供給され，新たな放電チャンネルを形成する．

以上のサイクルを1ショットとし，任意の時間間隔で駆動して推力を発生することができる．1サイクルとなる主放電は数10μsecの減衰振動波形を示す．

本研究では，このPPTを用いることにより，超小型衛星の高精度姿勢制御の実現を目指す．

図 2.6姿勢制御アクチュエータ⁸⁾

(18)

2. 3. 2 姿勢制御アルゴリズム

ここでは，人工衛星の姿勢制御アルゴリズムについて説明する．モーメンタムホイールのバイアス角運動量の変化により1軸を制御し，他の軸はジャイロ効果によって安定化を図る制御方式をバイアスモーメンタム方式と呼ぶのに対して，複数のアクチュエータを用いてロール軸（x軸），ピッチ軸（y軸），

ヨー軸（z 軸）を独立に制御する方式をゼロモーメンタム方式と呼ぶ．本研究の場合スラスタを制御用アクチュエータとして用いており，3 軸独立して制御する必要があるため，制御方式としてゼロモーメンタム方式を用いる．以下にゼロモーメンタム方式の制御アルゴリズムについて述べる．なお，座標系の定義やダイナミクスに関しては第3章で説明を行う．

衛星の姿勢角を[𝛷 𝜃 𝜓]^𝑇，軌道レートを𝜔₀，主慣性モーメントを[𝐼𝑥 𝐼_𝑦 𝐼_𝑧]，スラスタにより生じるトルクを[𝑇𝑐𝑥 𝑇𝑐𝑦 𝑇𝑐𝑧]^𝑇，環境外乱トルクを[𝜏𝑥 𝜏_𝑦 𝜏_𝑧]^𝑇とする．衛星が定常状態にあり姿勢角の変動が微小であると仮定し，運動方程式を線形化すると次式のように示される⁸⁾．

[

𝐼_𝑥𝛷^̈ 𝐼_𝑦𝜃^̈ 𝐼_𝑧𝜓^̈

] + [

𝜔₀²(𝐼_𝑦−𝐼_𝑧)𝛷

0

𝜔₀²(𝐼_𝑦−𝐼_𝑥)𝜓

] + [

𝜔₀(𝐼_𝑦−𝐼_𝑧−𝐼_𝑥)𝜓^̇

0

𝜔₀(𝐼_𝑦−𝐼_𝑥− 𝐼_𝑧)𝛷^̇

] + [ 0 0 0 ]

= [

−3𝜔₀²(𝐼_𝑦−𝐼_𝑧)𝛷

−3𝜔₀²(𝐼_𝑥−𝐼_𝑧)𝛷

0

] + [

𝜏_𝑥 𝜏_𝑦 𝜏_𝑧^{] + [}

𝑇_𝑐𝑥 𝑇_𝑐𝑦 𝑇_𝑐𝑧

]

(2.1)

軌道レートを微小として無視すると式(2.1)は次式のように表される．

[

𝐼_𝑥𝛷^̈ 𝐼_𝑦𝜃^̈ 𝐼_𝑧𝜓^̈

] = [

𝜏_𝑥 𝜏_𝑦 𝜏_𝑧^{] + [}

𝑇_𝑐𝑥 𝑇_𝑐𝑦 𝑇_𝑐𝑧

] (2.2)

バイアスモーメンタム方式では，ロール軸とヨー軸のカップリングがあり，互いに干渉しあうが，ゼロモーメンタム方式では式(2.2)に示すようにカップリング項及び，重力傾斜トルクが無視できるとすると，各軸を独立して設計，制御することが可能である．更にゼロモーメンタム方式は3軸が同一の式の形で表されているため，1軸周りに注目して設計を進めることができる．例としてロール軸に注目し，

その運動方程式を抽出すると次式で表される．

𝐼_𝑥𝛷̈ =𝜏_𝑥+𝑇_𝑐𝑥 (2.3)

式(2.3)が衛星のダイナミクスであり，姿勢制御系は姿勢決定系が推定する姿勢角をフィードバックし，

所望の姿勢に制御するために必要なトルクの算出を行う．図2.8に姿勢制御系のブロック線図を示す．

図 2.7 Pulsed Plasma Thrusterの推力発生メカニズム¹¹⁾

(19)

図 2.8姿勢制御系ブロック線図

(20)

第 3 章姿勢制御シミュレータ

本研究ではPPTをアクチュエータとして用いた超小型衛星の姿勢制御の検証をMATLAB/Simulink^®を用いて開発した姿勢シミュレータを用いて行う．姿勢シミュレータとは，軌道上での姿勢の決定及び制御を模擬するシミュレータであり，任意の軌道上，開始時刻からの軌道上の環境外乱を考慮し，衛星ダイナミクスから算出した姿勢角を姿勢決定センサの値からフィルタリングし，姿勢決定を行い，RW などのアクチュエータを用いて姿勢制御を模擬するものとなっている．これを用いることにより，軌道パラメータを任意の値に設定することで様々な軌道上外乱やその軌道における衛星の姿勢変化，消費電力の時間履歴を求めることが可能となる．

このシミュレータは環境モデル，衛星ダイナミクスモデル，センサモデル，姿勢決定・制御アルゴリズム，アクチュエータモデルから構成され，環境モデルと衛星ダイナミクスのブロックにより，環境外乱から衛星の真の姿勢情報を得ることができる．実際に運用される衛星では真の姿勢情報ではなく，センサによって観測された値から得られる推定値を利用するため真の姿勢情報を取得するセンサを模擬したセンサモデルを組み込み，それらのノイズやバイアス，ドリフトといった雑音情報を真の姿勢情報に付加しセンサの値を模擬している．この値を用いて姿勢決定アルゴリズムにおいて衛星の姿勢情報の推定を行っている．姿勢制御アルゴリズムでは，姿勢決定アルゴリズムにより決定した衛星の姿勢情報の推定値から各種アクチュエータへ所望の姿勢へ制御するための制御パラメータをアクチュエータモデルへ出力し，アクチュエータモデルにおいて，本研究ではPPTを模擬したものであるが要求されたトルクを姿勢制御アルゴリズムから送られてきたパラメータに準ずる駆動を行う．図3.1に姿勢シミュレータの概略図を示す．

本章では，本研究でPPTの搭載を仮定している超小型衛星ORBISの衛星モデルについて述べ，座標系の定義ともに，衛星のダイナミクスについて説明を行う．また，姿勢変化の原因となる環境モデルについても述べる．

図 3.1姿勢シミュレータ概略図

(21)

3. 1 人工衛星のダイナミクス

3. 1. 1 衛星モデル

本研究では，衛星モデルとして本研究室で開発を行っている超小型衛星 ORBIS のモデルを用いている．いかにORBISの諸元及び外観を示す．

表 3.1 ORBSI諸元

寸法 498×498×444[mm]

（パドルサイズ：420×440[mm]）

質量 50[kg]

重心^*1 X：3.4，Y：6.3，Z：238.5[mm]

*1) X軸については体心から観測窓方向，Z軸は分離面基準，

Y軸については体心からX軸，Z軸の右手系正方向

慣性モーメントIは太陽電池パドル展開後の値を3DCADを用いて計算した値を用いる．

𝐼 = [

1.508 −0.0105 0.0126

−0.0105 1.4630 0.0079 0.0126 0.0079 1.3910

] [kg・m²] (3.1)

3. 1. 2 座標系の定義

軌道上での姿勢角や一の表現のために，次の3つの座標系を定義する．

1) J2000.0赤道面座標系

地球を中心とする地心慣性座標系であり，西暦2000年の春分点方向をX軸正方向，北極点方向を Z軸正方向としてX軸，Z軸を基準とした右手系座標となるようY軸を定めたもの．慣性空間における衛星の軌道上における位置はこの座標系で表現される

2) 軌道座標系

軌道上の速度ベクトル方向をX軸，地球中心方向をZ軸とし，これらを基準としてY軸を定めたもの．

3) 機体座標系

機体の重心を中心とし，観測窓方向をX軸正方向，太陽電池パドル方向をZ軸正方向とし，それらを基準として右手系座標となるようY軸を定めた機体に固定された座標系．

本論文では次項より座標系の表現にJ2000.0赤道面座標系を{N}系，軌道座標系を{O}系，機体座標系を{B}系と定義する．

図 3.2 ORBIS外観

(22)

3. 1. 3 姿勢表現

衛星の姿勢表現にはいくつかの方法があるが，ここでは代表的な方向余弦行列，オイラー角，クオータニオン表現について説明する．

まず方向余弦行列であるが，これは任意の2つの座標系の変換を表す行列である．任意の2つの座標系については，ここではJ2000.0赤道面座標系と機体座標系を例にして説明を行う．J2000.0赤道面座標系の基底ベクトルを[𝒙𝒏 𝒚_𝒏 𝒛_𝒏]，機体座標系の基底ベクトルを[𝒙𝒃 𝒚_𝒃 𝒛_𝒃]とすると，これら 2つの基底ベクトルと方向余弦行列Cの関係は次式となる．ここで本論において基準座標系{N}系に対して{B}

系の方向余弦行列を𝑪_𝑵^𝑩と定義する．

[

𝑥_𝑏 𝑦_𝑏 𝑧_𝑏^{] = [}

𝐶₁₁ 𝐶₁₂ 𝐶₁₃ 𝐶₂₁ 𝐶₂₂ 𝐶₂₃ 𝐶₃₁ 𝐶₃₃ 𝐶₃₃

] [

𝑥_𝑛 𝑦_𝑛

𝑧_𝑛^{] =}𝑪_𝑵^𝑩[

𝑥_𝑛 𝑦_𝑛

𝑧_𝑛^] ^(3.2)

衛星の姿勢は基準となる座標系と衛星重心座標系，ここでは機体座標系との方向余弦行列として定義される．しかし，この行列そのものからは姿勢を直感的に把握しにくいため，オイラー角やクオータニオンと呼ばれる姿勢表現がある．J2000.0赤道面座標系を Z軸周り，Y 軸周り，X 軸周りに回転させると，回転後の軌道座標系は機体座標系と一致する．このときのZ軸，Y軸，X軸それぞれについての回転角をオイラー角とし，また，これらの角度をX，Y，Zの順にΦ（ロール角），θ(ピッチ角)，ψ（ヨー角）と呼ぶ．オイラー角は回転軸の選び方と順番によって 12 通りの表現が考えられるが，特に上記の順序で回転させるものを3-2-1系のオイラー角と呼ぶ．

方向余弦行列と3-2-1系のオイラー角との間には以下の関係が成立する．

𝑪 = [

1 0 0

0 cos𝛷 sin𝛷 0 −sin𝛷 cos𝛷

] [

cos𝜃 0 −sin𝜃

0 1 0

0 −sin𝜃 cos𝜃 ] [

cos𝜓 sin𝜓 0

−sin𝜓 cos𝜓 0

0 0 1

]

= [

cos𝜃cos𝜓 cos𝜃sin𝜓 −sin𝜃

sin𝛷sin𝜃cos𝜓 − cos𝛷sin𝜓 sin𝛷sin𝜃sin𝜓 + cos𝛷cos𝜓 sin𝛷cos𝜃 cos𝛷sin𝜃cos𝜓 + sin𝛷sin𝜓 cos𝛷sin𝜃sin𝜓 − sin𝛷cos𝜓 cos𝛷cos𝜃 ]

(3.3)

方向余弦行列と3-2-1系との関係は以下のとおりである．

図 3.3 座標系の定義

(23)

Φ =tan⁻¹𝐶₂₂

𝐶₃₃ (−𝜋 ≦ 𝛷 ≦ 𝜋) θ＝tan⁻¹ 𝐶₁₃

√𝐶₂₃²+𝐶₃₃²

(−𝜋

2≦ 𝜃 ≦𝜋 2)

ψ＝tan⁻¹𝐶₁₂

𝐶₁₁ (−𝜋 ≦ 𝜓 ≦ 𝜋)

(3.4)

cosθ = 0のとき，すなわちθ = ±π/2のとき特異点を持ち，計算することができない．オイラー角を使用する場合には，この特異点に留意する必要がある．

オイラー角の 3 変数で表現したときに生ずる数学的特異点を回避することができる姿勢表現として，

クオータニオン（オイラーパラメータもしくは四次元数と呼ばれることもある）が存在する．クオータニオンは座標系の回転がある直線周りの一つの回転とみなして表現するオイラーの定理を更に発展させて表現される．このときの直線をオイラー軸と呼ぶが，オイラー軸の方向単位ベクトルλを

𝛌 = 𝜆_𝑥𝑒_𝐵𝑥+ 𝜆_𝑦𝑒_𝐵𝑦+ 𝜆_𝑧𝑒_𝐵𝑧 (3.5) とし，この直線周りの回転角度（1直線に向きを与え，右回りを正の角度とする）をθとすると，クオータニオン𝐪 = [𝑞1 𝑞₂ 𝑞₃ 𝑞₄]^𝑇は次式で定義される．

𝑞₁= 𝜆_𝑥𝑠𝑖𝑛𝜃 2 𝑞2= 𝜆𝑦𝑠𝑖𝑛𝜃 2 𝑞₃= 𝜆_𝑧sin𝜃

2 𝑞₄= cos𝜃

2

(3.6)

式(3.6)に示すようにクオータニオンは4個のパラメータを用いて定義されるが，制約条件として

𝑞₁²+ 𝑞₂²+ 𝑞₃²+ 𝑞₄²= 1 (3.7)

があるため，実質的には3個の自由度を有することとなる．また，方向余弦行列とクオータニオンには以下の関係がある．

C = [

𝑞12− 𝑞22− 𝑞32+ 𝑞42 2(𝑞1𝑞2+ 𝑞3𝑞4) 2(𝑞1𝑞3− 𝑞2𝑞4) 2(𝑞₂𝑞₁− 𝑞₃𝑞₄) 𝑞₁²− 𝑞₃²− 𝑞₁²+ 𝑞₄² 2(𝑞₂𝑞₃+ 𝑞₁𝑞₄) 2(𝑞₃𝑞₁+ 𝑞₂𝑞₄) 2(𝑞₃𝑞₁− 𝑞₁𝑞₄) 𝑞₃²− 𝑞₁²− 𝑞₂²+ 𝑞₄²

] (3.8)

逆に方向余弦行列からクオータニオンを求める際にはまず式(3.7)の制約条件から方向余弦行列の対角項を用いてクオータニオンの1変数を算出し，その求めた1変数を用いて残りの変数を求める．1例として対角項から𝑞4を求めた際の関係を示す¹²⁾

(24)

𝑞₁=1

2√1 + 𝐶11+ 𝐶₂₂+ 𝐶₃₃ 𝑞2= 1

4𝑞₄(𝐶₂₃− 𝐶32) 𝑞₃= 1

4𝑞₄(𝐶₃₁− 𝐶₁₃) 𝑞₄= 1

4𝑞4

(𝐶₁₂− 𝐶₂₁)

(3.6)

3. 1. 4 衛星のキネマティクス方程式

衛星の角速度ベクトルをω = [𝜔𝑥 𝜔𝑦 𝜔𝑧]^𝑇とすると，角速度ベクトルωとクオータニオンqの間には次式の関係が成立する．

𝑑 𝑑𝑡𝒒̂ =1

2𝛀(𝛚)𝒒̂

=1 2[

−𝑞₁𝜔_𝑥− 𝑞₂𝜔_𝑦− 𝑞₃𝜔_𝑧 𝑞₄𝜔_𝑥− 𝑞₃𝜔_𝑦+ 𝑞₂𝜔_𝑧 𝑞₃𝜔_𝑥+ 𝑞₄𝜔_𝑦− 𝑞₁𝜔_𝑧

−𝑞2𝜔𝑥+ 𝑞1𝜔𝑦− 𝑞4𝜔𝑧

]

(3.10)

ここで，𝒒̂及び 𝛀(𝛚)は次式で定義される

𝒒̂ = [𝑞4 𝑞₁ 𝑞₃ 𝑞₃]^𝑇 (3.11)

𝛀(𝛚) = [ 0 −𝝎^𝑻 𝝎 −𝝎̃]

＝[

0 −𝜔_𝑥 𝜔_𝑥 0

−𝜔𝑦 −𝜔𝑧

𝜔_𝑧 −𝜔_𝑦 𝜔𝑦 −𝜔𝑧

𝜔_𝑧 𝜔_𝑦 0 𝜔_𝑥

−𝜔_𝑥 0

] (3.12)

式(3.12)をキネマティクス（運動学）方程式と呼ぶ．ここで𝝎̃は以下の式で定義される¹²⁾．

𝝎̃ = [

0 −𝜔_𝑧 𝜔_𝑦 𝜔_𝑧 0 −𝜔_𝑥

−𝜔_𝑦 𝜔_𝑥 0

] (3.13)

3. 1. 5 衛星のダイナミクス方程式

衛星に作用するトルクを𝝉 = [^𝜏𝑥 𝜏_𝑦 𝜏_𝑧]^𝑇，慣性モーメントをIとすると，トルクτによって衛星に生じる角速度は次式を解くことにより求められる．

𝑰𝝎̇ + 𝝎̃ 𝑰𝝎 = 𝝉 (3.14)

式(3.14)をダイナミクス（動力学）方程式と呼ぶ¹²⁾．

(25)

3. 1. 6 姿勢角の算出

本研究で用いたシミュレータでは式(3.14)のダイナミクス方程式から入力されるトルク τ を用いて衛星の角速度ωを算出しており，その角速度ωをもとに，式(3.10)のキネマティクス方程式からクオータニオンを算出している．算出したクオータニオンを式(3.4)，(3.8)を用いて方向余弦行列やオイラー角に適宜変換し，用いている．

3. 2 軌道計算

人工衛星は地球の周りを周回しており，ORBIS では軌道高度を550km としている．軌道上では地球の重力ポテンシャルの影響や太陽，月，惑星といった他天体の引力を受けるため一定の高度で一定の速度で周回することはできない．この現象のことを軌道の摂動というが，本シミュレーションでは摂動の影響も加味し，軌道上の位置，速度を TLE（Two-Line Element）¹³⁾の要素を基準にシミュレーションの時間経過とSGP4（Simplified General Perturbations Satellite Orbit Model 4）¹⁴⁾を用いて算出している．この説ではTLEとSGP4について説明し，本シミュレーションにおける軌道位置，速度の算出方法について述べる．表3.2に本シミュレーションで用いているORBISの想定軌道を示す．

表 3.2 超小型衛星ORBIS想定軌道

高度 550[km]

傾斜角 31[deg]

昇交点赤経 0[deg]

離心率 0

平均回転運動 15.05534468[周/日]

軌道1周周期 5760[sec]

3. 2. 1 TLE

TLEは2行要素とも呼ばれ，平均軌道要素の一つであり，2行のデータにより軌道情報を表すことからそう名付けられている．以下にORBISのTLEを例として挙げる．なおORBISの打ち上げは未定であり，まだ打ち上げられていないため次項で説明するSGP4に必要な要素だけをシミュレータ開始条件に合わせて設定し，シミュレータに組み込んでいる．