再 突入 カ プセル の運動 デー タか ら 空 力係数 を推定す る方法

平 成29年 度(2017年)学 位 論 文(修 士)

再 突入 カ プセル の運動 デー タか ら 空 力係数 を推定す る方法

ProcessofAssumingAerodynamicCoefficient fromRe‑EntryCapsule'sActualData

首 都 大 学 東 京 大 学 院

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

学修 番 号16891541

氏名 渡 辺 礼 奈

指導教員 佐原 宏 典 教授

平 成30年(2018年)2月22日

摘要

小 惑 星 探 査 機 は や ぶ さの よ うな 大 気 圏 再 突 入 型 カ プ セ ル は宇 宙 空 間 か ら地 上 に帰 還 す る 間 に,非 常 に 広 い 領 域 を 飛 行 す る.こ の 飛 行 領 域 に お い て 姿 勢 安 定 を実 現 す る た め に は 実 際 の 飛 行 デ ー タ を用 い て 再 突 入 カ プ

セ ル の 空 力 設 計 の 妥 当性 を 検 証 す る必 要 が あ る.こ れ ま で,再 突 入 カ プ セ ル の 空 力 係 数 は 風 洞 試 験,CFDに よ り取 得,推 定 され た が,は や ぶ さ2の 回 収 カ プ セ ル に は 飛 行 計 測 セ ン サ ー モ ジ ュ ー ル で あ るREMMが 搭 載 され て い る.そ の た め,地 球 帰 還 時 に は 回 収 カ プ セ ル の 母 船 か らの 分 離 か らパ ラ シ ュ ー ト解 散 ま で の 全 飛 行 領 域 にお い て3軸 加 速 度 お よび 角 速 度 を 取 得 す る こ とが で き る.そ こで 本 論 文 で は,実 飛 行 デ ー タ か ら再 突 入 カ プ セ ル の 空 力 設 計 の 妥 当性 を検 証 す る た め に,カ ル マ ン フ ィル タ を 用 い て 空 力 係 数 を推 定 す る 方 法 の 確 立 を 目的 とす る,

は じめ に,再 突 入 カ プ セ ル の 支 配 方 程 式 で用 い る 各 パ ラ メ ー タ を定 義 し,観 測 値 か ら推 定 す る値 の 理 論 式 を 示 す.次 に カ ル マ ン フ ィル タ の フ ィル タ リン グ 設 定 お よび 推 定 パ ラ メー タ の 同 定 問 題 に つ い て ノイ ズ の 最 尤 推 定 値 や 共 分 散 行 列 に お け る相 関 関 数 を 算 出 し,検 討 した,そ の 後,カ ル マ ン フ ィル タ に よ る シ ミ ュ レー シ ョン の 運 動 デ ー タ か ら姿 勢 角 ・角 加 速 度,速 度 ・位 置 お よび 各 空 力係 数 の 推 定 結 果 を 示 した.こ の 推 定 方 法 を 検 証 す る に あ た っ て,は や ぶ さカ プ セ ル の 設 計 の た め に 空 力係 数 を 入 力 す る こ とで ノ ミナ ル 計 算 に よ る 再 突 入 カ プ セ ル の 運 動 デ ー タが 算 出 され るプ ロ グ ラ ム を用 い た,算 出 され た 運 動 デ ー タ に ノイ ズ を のせ,そ れ ら を観 測 値 と して カル マ ン フ ィル タの 空 力 係 数 推 定 プ ロ グ ラ ム に 入 力 す る こ と で,観 測 値 の ノイ ズ が 誤 差 補 正 され た 上 で ノ ミナ ル 計 算 時 に 入 力 して い た 空 力 係 数 を 求 め た.こ れ に よ り,本 研 究 に お け る,シ ミュ レー シ ョン に よ る 観 測 デ ー タか ら空 力 係 数 を推 定 す る 方 法 の 妥 当性 を 検 証 した.

図 目次

1231234567123             l112222222333

図図図図図図図図図図図図図

は や ぶ さ2カ プ セ ル の ノ ミナ ル 計 算2り

再 突 入 カ プ セ ル の 母 船 分 離 か ら慣 性 座 標 系 に 固 定 され る ユ2)

は や ぶ さ2カ プ セ ル に 搭 載 され た飛 行 計 測 セ ン サ ー モ ジ ュ ー ルREMM23}

機体重心座標系B系 回転地球座標 系E系 地心座標系C系 ZCま わ りに甲回 転

YB1ま わ り に θ回 転 XB2ま わ り に ψ回 転

再突入 カプセルの空力中心位置

233566788

フ ィル タ リン グ手 法 の 概 要 図 3軸 角 速 度 の フ ィル タ リン グ結 果 κ,且軸 姿 勢 角 の フ ィル タ リ ン グ 結 果 図3.4ア 軸 姿 勢 角 の フ ィル タ リ ン グ結 果 図353軸 加 速 度 の フ ィル タ リン グ結 果 図3.63軸 速 度 の ブ イ ル タ リ ン グ 結 果 図3.73軸 座 標 位 置 の フ ィ ル タ リ ン グ 結 果

図3.8 図3.9 図3.10 図3.ll 図3.12 図3.13

大 気 圏 再 突 入 後200秒 間 の 軸 力係 数 の推 定 結 果 大 気 圏 再 突 入 後200秒 間 の モ ー メ ン ト係 数 の推 定 結 果

大 気 圏 再 突 入 後200秒 間 の ダ ン ピン グ係 数 の 推 定 結 果 はや ぶ さカ プ セ ル の ノ ミナ ル 計 算 に よ る マ ッハ 数 は や ぶ さ カ プ セ ル の ノ ミナ ル 計 算 に よ るα亡o亡

は や ぶ さ カ プ セ ル の ノ ミナ ル 計 算 に よ る空 力係 数(M・1)

3523344566778811222222222222

目次

第1章 序 論

1,1研 究 背 景 L2課 題 1.3研 究 目 的

‑ニー2弓﹂

第2章 空 力 係 数 の 推 定 課題

2.1使 用 す る 無 次 元 量 ・用 語 の 定 義 2.2座 標 変 換 行 列

2.3再 突 入 カ プ セ ル の 支 配 方 程 式 と空 力 係 数 2.4観 測 値 か ら推 定 す る値

44くソ0ノ

2,4.1姿 勢 角 ・角 加 速 度 2.4.2速 度 ・位 置

2,4.3空 力 中 心 を 原 点 と し た 重 心 位 置

12 12 12 13

第3章 カ ル マ ン フ ィル タ を用 い た 推 定

3.1フ ィ ル タ リ ン グ 設 定

3.1.1 フ ィル タ リ ン グ手 法 の概 要 3、1.2カ ル マ ン フ ィ ル タ に 基 づ く漸 化 式

3,1,3姿 勢 角 の 推 定 に お け る 各 パ ラ メ ー タ の 定 義 3,L4速 度 ・位 置 の 推 定 に お け る 各 パ ラ メ ー タ の 定 義 3,2推 定 パ ラ メ ー タ の 同 定 問 題

3,2,1観 測 方 程 式 ・状 態 方 程 式 に お け る ノイ ズ の最 尤 推 定値 3.2.2観 測 値 の 平 滑 化

3.2.3カ ル マ ン ゲ イ ン の 決 定 3.2.4共 分 散 行 列 の 更 新 3.3シ ミ ュ レ ー シ ョ ン 結 果

3.3.1角 速 度 ・姿 勢 角 3.3.2加 速 度 ・速 度 ・位 置 3,3,3空 力 係 数

3.4考 察

55567899011224591111111122222222

第4章 結 論

4.1結 論

4.2今 後 の課 題 と展 望

0ハU(U3内﹂﹃﹂

参考文献 ³¹

表 目次

1

表表表表表表表 ーへ∠つ﹂﹂4戸)!0¶且

       2222う臼23

座 標 系 の 定 義24)̲̲̲̲.,,̲

座 標 変 換 行 列TCEで 使 用 す る 記

並 進 運 動 の 各 パ ラ メ ー タ 分 別̲̲̲....,̲̲.̲̲̲.̲.̲.̲

回 転 運 動 の 各 パ ラ メ ー タ 分 別...,̲̲̲̲̲̲....̲

用 い る 既 知 の 値̲.̲̲̲̲̲̲̲̲̲̲.̲̲

空 力 中 心 を 原 点 と し た 重 心 位'置 の 計 算 で 使 用 す る 記 号̲̲̲̲.̲

各 空 力 係 数 の 推 定 値 と ノ ミ ナ ル 計 算 の 比 較(M=D.̲

●響,● ●,■ ■●●■●,■ ■●,● ●■●,■ ■●.■ ・層,・ ・,噛・,・,.・,唱 ・..・.

万'̲̲̲̲̲.,...師̲̲̲̲̲.̲̲・ 岬 ^{●● ● ● ●,■●,■}●■ ● ●● ■ ● ●

●昏 ■ ●● ● ●● ● ■● 噛 ●,■● 唱,● ●● 幽.■.■.■ ・.・.

●■ ●● ●,●● ■■,●●,●,■ ■ ●■ ・● ■ ・■

● ●● ●噛 ● ●■ ●,■●,■,■ ●,■●,■,■ 」,

,●● ■ ●●,■● ■● 唱,●,■ ●,.■ 」 ■ ■・

● ■● ■● ● ■,■,,●,幽 ■,■,● ■,,」 ●

..5 7 9 10 ..̲̲̲IO

13 29

第1章 序 論

本 章 で は,再 突 入 カ プ セ ル の 空 力 設 計 の 変 遷 と,母 船 分 離 か ら地 球 帰 還 ま で の 飛 行 中 に お け る 主 な 現 状 の 課 題 を 二 点 挙 げ る,そ の 後,そ れ らの 課 題 の解 決 策 と して 本 研 究 の 目的 を 示 す,

1.1研 究 背 景

近 年,宇 宙 活 動 に お け る 物 資 回 収 や 小 惑 星 な どの 天 体 か らサ ン プル を 地 球 へ 持 ち帰 る大 気 圏 再 突 入 カ プ セ ル 形 状 の検 討 が 複 数 行 わ れ て い る,国 際 宇 宙 ス テ ー シ ョン(ISS)か らの物 資 回 収 を 目的 と して,宇 宙 ス テ ー シ ョン 補 給 機 こ うの と り(HTV)に 回 収 機 能 を 付加 したHTV‑Rの 検 討 は 将 来 の 有 人 宇 宙 活 動 を行 う 上 で非 常 に 重 要 な 技 術 で あ る1},

一 方,天 体 か らサ ン プ ル を 持 ち 帰 るサ ン プ ル リタ ー ン に つ い て,惑 星 の 中 で も特 に 火 星 は 生命 が 存 在 す る 可能 性 が 最 も 高 く,欧 米 露 で は 精 力 的 に 探 査 が 行 われ て き た.日 本 で も現 在,探 査 計 画 が 進 め られ て お り,そ の 中 で も火 星 か らサ ン プ ル を 持 ち 帰 る火 星 サ ン プル リタ ー ン 計 画 等 に用 い られ る火 星 大 気 圏 突 入 カ プ セ ル 形 状 の 検 討 も複 数 行 わ れ て い る2).再 突 入 カ プ セ ル の 形 状 の 中 で 最 も頻 繁 に 用 い られ るsphere‑cone 形 状 は,搭 載 シ ス テ ム を 単 純 化 す るた め に,受 動 的 な 空 力 安 定性 に よ る 安 定 な 飛 行 の 達 成 を 要 求 され る こ

と が 多 い,こ の 動 的 な 安 定 性 を評 価 す る た め に は2種 類 の 実 験 的 方 法 が 考 え られ る.

第 一 は 風 洞 設 備 を 用 い て 模 型 を 振 動 させ る振 動 法 で あ り,さ らに 大 別 す る と 自 由振 動 法,自 由 回 転 法 及 び 強 制 振 動 法 の 三 っ に な る,Wehrend3Pt}s}は 惑 星 探 査 プ ロー ブ の 空 力 特 性 を調 査 す る 目的 で,半 頂 角10。

の 鈍 頭 円 錐 体 の 静 的 な 空 力 係 数 及 び ピ ッチ 軸 回 りの動 的 な ダ ン ピ ン グ空 力 係 数 を強 制 振 動 法 に よ り測 定 し た.O.65〈M〈2.2の 領 域 で迎 角 αニ18⑪ ま で の 空 力 係 数 を取 得 して い るが,マ ッハ 数Mニ1付 近 で 迎 角 が0。

〈α<5。 で は 動 的 に 不 安 定 に な る こ と,特 に 円 錐 状 の 背 面 形 状 を 有 す る もの に つ い て は フ ラ ッ トな 背 面 形 状 の もの と比 較 して そ の 傾 向 が 顕 著 に 見 られ,測 定 した ほ とん どの マ ッハ 数 領 域 で動 的 に 不 安 定 で あ っ た.

ま た,前 面 形 状 を鈍 頭 に す る こ と及 び 背 面 形 状 を球 状 に す る こ とは ど ち ら も動 不 安 定 で あ る こ とが わ か っ た,さ らに,面 に 円 状 の 切 り込 み を設 け る こ とで 亜 音 速 領 域 で 動 安 定 で あ る が,音 速 領 域 で は 動 不 安 定 と な っ た.Wrightら6)はMercurifGeminiに つ い て 強 制 振 動 法 を用 い て05<M<4.63に お け る ピ ッチ,ヨL・一一F軸回 りの 動 的 な ダ ン ピ ン グ空 力 係 数 を 調 べ た.そ の 結 果,ダ ン ピ ン グ 空 力 係 数 は迎 角 に 対 して非 線 形 な 変 化 を 示 した.Apollo計 画 の 後 期 に は 司 令 船(CommandModule)に つ い て 自由 回 転 法 を 用 い て 試 験 した と こ ろ, 0.3〈M<0.8に お い て 再 突 入 時 の 姿 勢 で あ る 迎 角180。 で は 動 的 に 不 安 定 で あ っ た7},動 的 な 安 定 性 を 評価 す

る 第 二 の 試 験 法 は 実 際 に 飛 翔 体 を飛 翔 させ る 自 由飛 行 法 で あ る.Krumins8)9)は 飛 翔 体 を 高 速 で 射 出 す る 弾 道 飛 行 装 置(バ リス テ ィ ッ ク レン ジ)を 用 い て 半 頂 角60。 の 鈍 頭 円錐 体 カ プ セ ル に つ い て 動 安 定 性 を 調 査 し,M=1付 近 で 動 不 安 定 で あ り,M=4以 上 で 動 安 定 で あ る とい う傾 向 が 得 られ た.ま た,Sammondsi。 川 は バ リステ ィ ッ ク レ ン ジ を 用 い て 半頂 角55。 と60。 の2つ の 鈍 頭 円 錐 体 カ プ セ ル に っ い て 動 安 定 性 を 調 査 し,1.0〈M<1.4で 動 不 安 定 で あ っ た.河 本12)は超 音 速 風 洞 内 にMercury型 の カ プ セ ル を 初 期 迎 角30。 で 打 ち 出 し,1.5〈M<3.0で 振 幅 は 発 散 して動 不安 定 で あ った.ま た,こ れ らの 動 的 な 安 定性 を 評 価 す る2種 類 の 実 験 法 が 確 立 され て く る と,両 者 の結 果 を 比 較 す る調 査 が 行 わ れ た,吉 永 ら13)14)IS)16}17)IS}は,軌道 再 突 入 実 験 機OPEX(OrbitalReentryExperiment)に つ い て,遷 音 速 か ら極 超 音 速 の 広 い マ ッハ 数 範 囲 に わ た っ て 自 由 回 転 法 に よ る動 的 風 洞 試 験 を 行 っ た.OREXは 極 超 音 速 領 域 のM=7.1で は動 的 に ほ ぼ 安 定 とみ な す こ とが で きた.し か しな が ら,0.7<M<2、Oの 遷 音 速 か ら超 音 速 領 域 で は振 動 は 発 散 し,あ る振 幅 で リ ミ

ッ トサ イ クル に な っ た.振 動 の 最 大 振 幅 は1.OくM<1.1で 観 察 され,20。 程 度 ま で 達 した,ま た,OREXの 実 機 飛 翔 結 果 の うち推 定 飛 行M=1付 近 で は 風 洞 試 験 と同 程 度 の 振 動 が 発 生 した,こ の よ うに,こ れ ま で 行 わ れ た 再 突 入 カ プ セ ル の 安 定性 の 研 究 か ら,静 的 に は安 定 で あ っ て も亜 音 速 か ら遷 音 速 に お い て 動 的 に 不 安 定 とな る傾 向 が 報 告 され て い る,こ の 動 不 安 定 性 は,カ プセ ル 背 面 の 非 定 常 的 な圧 力 変 動 に依 存 して い

る と考 え られ て い る が,未 だ 十 分 な 知 見 は得 られ て い な い,そ の た め,再 突 入 カ プ セ ル の 運 用 の た め に は 亜 音 速 か ら遷 音 速 に か け て の 動 安 定 特性 を把 握 す る こ とが 重 要 で あ る,さ らに,Useltenig)に よ る火 星 探 査 プ ロー ブ(Viking)の 風 洞 試 験 及 び 飛 翔 試 験 で は動 不 安 定 性 の 存 在 を 明確 に す る と同 時 に,風 洞 試 験 で は 模 型 を支 持 す る ス テ ィ ン グ の 太 さ,長 さ をパ ラ メー タ と してVikingの 動 的 不 安 定 性 へ 与 え る影 響 に つ い て の 調 査 も行 っ た.大 き な剥 離 を伴 う亜音 速 流 れ に お い て,非 定 常 現 象 に お け る 支 持 干 渉 の影 響 は 顕 著 で あ る.特 に後 部 が 丸 い カ プ セ ル 形 状 は 実 飛 行 で 剥 離位 置 が 固 定 され ず,物 体 の 運 動 と相 互 に影 響 し合 うた め, この よ うに 従 来 の 風 洞 試 験 で は 支 持 装 置 の 影 響 が 大 き く正 確 な測 定 が 困難 で あ っ た.そ こで 近 年 で は,動 安 定 特 性 調 査 の第 一 歩 と して,宇 宙 航 空 研 究 開発 機構(JAXA)の 支 持 干 渉 が 無 い磁 力 支 持 天 秤 装 置(MSBS) 用 い て,再 突 入 型 カ プ セ ル の 軸 力係 数 お よび モ ー メ ン ト係 数 を算 出 す る試 験 も行 わ れ て い る が,磁 力 支 持 の 可 能 範 囲 の 制 約 に よ り,M=0.05未 満 に お け る非 常 に 低 い 速 度 領 域 にお け る測 定 が 実 現 され て い る こ とが 現 状 で あ る20),

1.2課 題

小 惑 星 探 査 機 は や ぶ さの よ うな 大 気 圏 再 突 入 型 カ プ セ ル は 宇 宙 空 間 か ら地 上 に 帰 還 す る問 に,図1,1に 示 す 極 超 音 速 か ら亜 音 速 領 域,自 由分 子 流 か ら連 続 流 領 域 ま で の 非 常 に広 い領 域 を 飛 行 す る.

1⑪ 10 OOOOO

﹂￡彗Z岩省5藷

105 106

腿 湘 購'匹 蜘1"̀u''"1̀二

A

一

Tran5i髄on臼1nOW 轡

鞠

一

一

i

一

一

一

‑

一

A

旧

■ Con血uou5Flow

辱

一

'圃,h嘲1,1,,、,、 曜』,†1,,

㎞ 晦 」Win曲 螂1

、 』HIhIl‑、

o ¹⁰ 2030

MachNumber

40

図1,1は や ぶ さ2カ プ セ ル の ノ ミ ナ ル 計 算2D

こ こ で,再 突 入 カ プ セ ル の 飛 行 領 城 に お い て 姿 勢 安 定 を実 現 す る た め に検 証 す べ き 主 な課題 が 二 点 あ る.

一一点 目は,自 由分 子 流 領 域 に お い て 母 船 か らカ プ セ ル が 分 離 され る 際 に慣 性 座 標 空 間 に 固 定 され,分 離 時 の 姿 勢 の ま ま 大 気 圏 に 再 突 入 す る こ とで あ る.(図L2)こ こ で は,空 力 加 熱 の 影 響 を 小 さ くす る た め に 迎 角 0。 の 姿勢 に 固 定 され る こ とが 理 想 で あ るが,0.3Hzの ロー ル ス ピ ン状 態 で カ プ セ ル は母 船 か ら分 離 され る, そ こで,ロ ー ル ダ ン ピ ン グ係 数 を は じめ とす る空 力 係 数 を求 め る こ とで,分 離 擾 乱 の 影 響 を検 証 す る 必 要 が あ る.二 点 目は,従 来 の 風 洞 試 験 な どか ら挙 げ られ る遷 音 速 領 域 で の カ プ セ ル の動 的 不 安 定 性 に っ い て,亜 音 速 領 域 で の パ ラ シ ュー ト開 傘 へ の影 響 を検 証 す る 必 要 が あ る.

＼

/・'

/ノL

̲!

一"Hl

図1.2再 突 入 カ プ セ ル の 母 船 分 離 か ら慣 陸座 標 系 に 固 定 され る22)

1.3研 究 目 的

これ ま で,広 い 領 域 を飛 行 す る 再 突 入 カ プ セ ル の 空 力係 数 は 風 洞 試 験,流 体 解 析(ComputationalFluid Dynamics:CFD)に よ り取 得,推 定 され て き た.小 惑 星 探 査 機 はや ぶ さの 回 収 カ プセ ル は 設 計 上 の 制 約 か ら 運 動 計 測 セ ン サ 類 は 非 搭 載 で あ っ た が,は や ぶ さ2の 回 収 カ プセ ル に は飛 行 計 測 セ ンサ ー モ ジ ュ ー ル で あ る REMM(図1.3)が 搭 載 され て い る.そ の た め,地 球 帰 還 時 に は 回 収 カ プ セ ル の 母 船 か らの 分 離 か らパ ラ シ ュ ー ト解 散 ま で の 全 飛 行 領 域 に お い て3軸 加 速 度 お よび 角 速 度 を取 得 す る こ とが で き る.

以 上 よ り本 研 究 の 目的 は,実 飛 行 デ ー タ か ら再 突 入 カ プ セ ル の 空 力 設 計 の 妥 当性 を検 証 す る た め に,第 一 段 階 と して

,シ ミュ レー シ ョ ンの 運 動 デ ー タ か ら空 力 係 数 を推 定 す る 方 法 を確 立 す る こ と とす る.

二軌n血

「「ntrn

図1,3は や ぶ さ2カ プ セ ル に 搭 載 さ れ た 飛 行 計 測 セ ン サ ー モ ジ ュ ー ルREMM2a‑}

第2章 空 力係 数 の 推 定課 題

本 章 で は,再 突 入 カ プ セ ル の 支 配 方 程 式 の 各 パ ラ メー タ お よび そ れ らを 推 定 す る に あ た っ て 必 要 と な る 値 を 示 し,本 研 究 で推 定 す る軸 力 係 数 モ ー メ ン ト係 数 お よび ダ ン ピン グ係 数 を 定 義 す る,ま た,観 測 値 の 座 標 系 と再 突 入 カ プセ ル の 支 配 方 程 式 にお け る既 知 の値 の 座 標 系 が 異 な る こ とか ら,そ れ らの 座 標 変 換 行 列 の 計 算 方 法 に っ い て も示 す.

2.1使 用 す る 無 次 元 量 ・用 語 の 定 義

本 研 究 で は 以 下 の 無 次 元 量 を使 用 して 解 析,考 察 を行 う.

マ ツ ハ 数:班

流 体 の 圧 縮 性 の 尺 度 と な る 無 次 元 量 で あ り,次 式 で 定 義 さ れ る.

こ こ で,u=流 速[m/s],a=音 速[m/s]で あ る,

び

、M=一 α

(2,1)

レ イ ノ ル ズ 数:Rε 数

物 体 回 り の 流 速 を レ イ ノ ル ズ 数 を 川 い て,無 次 元 化 す る.レ イ ノ ル ズ 数 は 以 下 の 式 で 定 義 す る, こ こ で,v:一 様 流 速 度[m/s],L=再 突 入 カ プ セ ル の 最 大 直 径[m],ユ ・:空 気 の 動 粘 性 係 数[m/s2]で あ る,

VL Reニ ー Ψ

(2.2)

抗 力 係 数:CD

軸 力 係 数 に つ い て,一 例 と して 抗 力 係 数CDは 以 下 の 式 で 定 義 す る,同 様 に し て 揚 力 係 数 免 を 定 義 す る, こ こ で,D=抗 力 囲,q二 三ρv2=動 圧[7V/m2】,S=カ プ セ ル の 前 方 投 影 面 積[m2]で あ る,

2

̀・ 「 葬D (2,3)

モ ー メ ン ト係 数=̀拠

モ ー メ ン ト係 数 に つ い て,凡 例 と し てCmは 以 下 の 式 で 定 義 す る.

こ こ で,M;モ ー メ ン ト[N・m],d=カ プ セ ル の 全 長[m],S=カ プ セ ル の 前 方 投 影 面 積[m2]で あ る.

̀肌=顧M (2.4)

2.2座 標 変 換 行 列

本 節 に お い て 座 標 変 換 行 列 を 示 す に あ た っ て,各 座 標 系 の 定 義 を 表2.

に 示 す24).

表2.1座 標 系 の 定 義24}

お よ び 図2」,図2.2,図2.3

座 標 系名(英) 座標 系名(和) 原 点,特 徴 方 向 ・エ ポ ッ ク

Body‑flxed B

(Rotating)Ea!寸11CS E

LocalEarthCS, GeocentricCS C

機 体重 心座標 系

回転地 球座標 系

地 心座標 系

機 体重心

地球 と共 に回転 回転 軸:E2

地球 中心

機 体 と共 に移 動 機 体重 心中心

Bl:mainaxis(nosedirection) B2=normaltosymmetricplane

(right‑direction) B3:BlxB2

El=赤 道 面 内,Greenwich子 午

面 を 貫 く 方 向

E2:地 球 ス ピ ン 軸 北 側 向 き E3:ElxE2

Cl:North

C3=地 球 中 心 方 向 C2=C3xCl

、 Pitch

y'≡z'×Xl'一 、

、ぐ ・」'『 一ゴ ∵ プ

Roll

Yaw

X,Z1

図2.1機 体 重 心 座 標 系B系

Spinaxis

汀.

唱L

㌔ 卜 乳

扇

qノ の。

z。 、

Ψ滑

"姦

P

Equator i藍hmεridian

図2.2回 転 地 球 座 標 系E系

晶Z

で,

」●」

,㌧

盈「ド7

指

=

、 zε

'

'ヘ ド

」 、

「

λ

x。 Eqロator

y.,

Greehd

図2.3地 心 座 標 系C系

観 測 値 で あ る3軸 加 速 度 お よび 角 速 度 は 機 体 座 標 系 のB系 で 取 得 され る た め,支 配 方 程 式 の 各 パ ラ メ ー タ をB系 表 示 し,シ ミ ュ レー シ ョ ン を 行 っ た,既 知 の 重 力 加 速 度 は 回 転 地 球 座 標 系E系 にお け る値 で あ る た め,ま ず 第 一 に,E系 か ら機 体 重 心 中 心 の 地 球 固 定 で あ る座 標 系C系 へ 座 標 変 換 し,そ の 後 オ イ ラ ー 角 の 方 程 式 か ら24),C系 か らB系 へ 座 標 変 換 した,

E系 か らC系 へ の座 標 変 換 行 列TCE

使 用 す る記 号 の 定 義 を表2.2に 示 し,以 下1‑3の 各 軸 回 りの各 回 転 に よ っ て 座 標 変 換 行 列TCEを 式(25) に 示 す.こ こ で,そ れ ぞれ の 回転 角 度 は機 体 の 角 運 動 に よ る も の で,観 測 値 の 角 速 度 か ら推 定 で き る.

LYEま わ り に γR回 転, 2.XEま わ り にdμ 回 転, 3.YEま わ り に 一γ回 転

表2、2座 標 変 換 行 列7"CEで 使 用 す る 記 号

記号 定義

YR

μR

γμμd

再 突 入 点 緯 度=oe 再 突 入 点 経 度 μR=

β亡(βは 地 球 時 点 角 速 度) 再 突 入 飛 行 中 の 現 時 点 緯 度

再 突 入 飛 行 中 の 現 時 点 経 度 再 突 入 点 か らの 経 度 差(μ 一 μR)

TCE・:溜1=腐)[is:細

):農1甥:灘 謂 一[1:鷲

)i:器:i轍):塊 甑:窯:;に:蹴)

‑r潔灘 鞭 膿 瓢 二::誰鵜:1鷺 繍 繍 蹴

C系 か らB系 へ の 座 標 変 換 行 列TBC

オ イ ラ ー 角 の 方 程 式 か ら,以 下1‑3を 計 算 し,座 標 変 換 行 列TBCを 式(2.6)に 示 す.

LZcま わ り にep回 転(図2.4) 2.YB1ま わ り に θ回 転(図25) 3.XBZま わ り に ψ回 転(図2.6)

修 ・…

、.■e1￠

図2.4Zcま わ り に ψ回 転

灸 ≒ ら 1

(2.5)

秘2 ^{Pel≒ 」e}

e2

*el

ee2

ち モ ら

、 図2.5YBIま わ り に θ回 転

*e2

劣e3

ee3

Pe3

㌦ ・、・

弟ei

ee5

ち モ2 ei 図2,6XB2ま わ り に ψ 回 転

[7'・]BC=【離 巽1篶 署 糠 器

sinth̀osθ siniPsinθsi冗q+cosψcosep

sinψsinθ̀os(P‑cosiPsin{P 熱1 ^(2.6)

従 っ て,式(2.5),(2.6)か らE系 か らB系 へ の 座 標 変 換 行 列T朋 は,

TBE=TCETBC (2,7)

と 求 め られ る,

2.3再 突 入 カ プセ ル の 支配 方程 式 と空 力 係 数

再 突 入 カ プ セ ル の 飛 行 中 の 支 配 方 程 式 は 並 進 お よび 回 転 運 動 に つ い て そ れ ぞ れ 検 討 し,そ れ ぞ れ の パ ラ メー タ を観 測 値,観 測 値 か ら推 定 す る値 に 分 別 した 後 にそ れ ら の 関 係 か ら,各 空 力 係 数 の 式 を 定 義 す る,

並 進 運 動

並 進 運 動 方 程 式 は,

轍 肌囲 州1+m[9]B

で 与 え られ る.こ こ で,観 測 値,観 測 値 か ら推 定 す る値,既 知 の 値 の 分 別 を 表23に 示 す.

表2.3並 進 運 動 の 各 パ ラ メー タ分 別

(2.8>

観測値 観測値か ら推定する値 既知の値

[α彦]B,[nBE]B ["吾]E,[9]B m,q,s

ま た,[ρBE】Bは 回 転 地 球 座 標 系E系 に 対 す るB系 の 相 対 角 速 度 の 交 代 行 列 をB系 表 示 し た

0‑‑rq [nBE]且 二rO‑P

‑‑qPO

で あ る,

2.2節 の 式(2.5),(2.6),(2,7)よ り,既 知 のE系 に お け る 重 力 加 速 度 は,B系 表 示 し,

9xB

gy=[T]日E[9]E=[T]̀E[T]BC[9]E gz

c。s(γR)COS(γ)+sin(γR〕C。S(dμ)sin(y)‑sin(dμ 〕sin(y)‑sin(γ 月)cos(γ)+COS(γR)c。S(α μ)sin(γ)

=sin(γR)sin伽)c。s(dμ)cos(γR)sin(dμ)

‑c。s(YR)stn(γ)+sin(γR)COS(dμ)COS(y)一 ・sin(dμ)cos(y)sin(γR)sin(γ)+COS(γR)C。s(dμ)sin(y)

cosψcosesiniPcosθ 一sinθ9xE

・cosiPsinθstnep‑sinψcosqsiniPsinθsinOP+COSψcosqCOSθsinOPgy cosiPsinθcosq+sinψsinqsinthsineCOSOP‑cosiPsinOPCOSθcosqgz

さ ら に,E系 に つ い てB系 表 示 し た 角 速 度 を オ イ ラ ー 角 の 微 分 方 程 式 で 示 す と,

「11=

1s̀η ψ̀α れθ̀o∫ ψ 亡αnθ Ocosrp‑sinep

ost冗q/cosθcosep/cosθ

斜F可r

(29)

(2.10)

(2.11)

が 導 出 さ れ る24).

回転運動

回転運動方程式は,

躍繍+一 騨 一欄 一匪F榛}(2・12)

で 与 え られ る,右 辺 第 一 項 は モ ー メ ン ト係 数 で 示 す モ ー メ ン ト項 で あ り,右 辺 第 二 項 は ダ ン ピ ン グ係 数 で 示 す 減 衰 項 で あ る.こ こで,観 測 値,観 測 値 か ら推 定 す る 値,既 知 の 値 の分 別 を 表2.4に 示 す.

表2.4回 転 運 動 の 各 パ ラ メ ー タ分 別

観測値 観測値から推定す る値 既知の値

[ωBE】B,[nBE]B

Bd ωBE

d亡 [18]B,q,S,d

ま た,[1創Bは カ プ セ ル が 軸 対 称 で あ る こ と か ら,

で あ る.

昭F=齎

ま た,以 降 の 計算 で 用 い る 各 既 知 の 値 は 表25を 用 い る.

表25用 い る既 知 の 値

(2.13)

既知量 代入値

質 量 搬[kg]

翼 面 積 ∫[m2]

代 表 長d[川

千'覧衰 巾畠̀[襯]

慣 性 モ ー メ ン ト∬[kg・m2]

165

0,126

0.2

0.4

∬x==023741γ ニ0.!711 1zニ0.1687

以 上 よ り,各 空 力 係 数 の式 を 定 義 す る.

軸 力 係 数

式(2.8),(2.9),(2.10)よ り,

抗 力 係 数 ・揚 力 係 数 式(2.14)お よ び

圏一器[雛1徽]

CDニ ーCxcosα 一Czsina(空 力 中 心 の κ 軸 方 向 空 気 力)

CL・Cxsinα 一Czcosα(空 力 中 心 のz軸 方 向 空 気 力)

よ り,

[豊1一器[准『郷鷲 望 潔瓢1㍉顎協望l!糊

モ ー メ ン ト係 数

式(2,12),(2,13)よ り,

【多1‑1≡llii≡箋1

こ こ で,(XCG,YCG,zCσ)は 空 力 中 心 を 原 点 と し た 重 心 位 置 で あ り,2.4,3で 示 す,

ダ ン ピ ン グ 係 数

式(2.12),(2.13),(2.18)よ り,

ま た,d=機 体 横 基 準 長,

E(一 らy、、+CyZ、 の

1刻 一1噺+c・XCG)

F(‑CyXCG‑CxYca)

̀=機 体 縦 基 準 長 と す る,

1ψ 十qr(iz‑ly) P llyq十rp(1x‑1且) q∫q

∬zfa十pq(1y‑lx) r

(2」4)

(2.15)

(2,16)

(2.17)

(2.18)

(2」9)

2.4観 測 値 か ら 推 定 す る 値

2.3節 に お い て,表2.3や 表i2Aで 示 した,観 測 値 か ら推 定 す る 各 値 に っ い て,以 下 の 方 針 を 検 討 す る, 2.4.1姿 勢 角 ・角 加 速 度

姿 勢 角[γ]Bニ[̀ipθ ψ工B

t=Oの と き,Poニ0ψ 。ニ0と す る と,そ の 後 の 時 間 経 過 で 角 速 度 は 計 測 で き る た め,

亡一sの と きq

、=q(s.、)+P,{s‑(s‑1)}(2・20)

で あ り,亡 二 〇 か らy,z軸 に 関 し て も 同 様 に 各 姿 勢 角 が 連 続 的 に 求 め られ,こ れ に よ り,式(2,10)のB系 に お け る 重 力 加 速 度 が 算 出 可 能 で あ る.

  ほ 角 加 速 度[ωBE]ニ[Pqナ]B

亡=0の と き,p。 ニ0戸 。ニ0と す る と,そ の 後 の 時 閲 経 過 で 角 速 度 は 計 測 で き る た め,

亡一sの ときPs‑

{1≒ ξ笥(2・21) で あ り,亡 二 〇 か らy,z軸 に 関 して も 同 様 に 各 角 加 速 度 が 連 続 的 に 求 め られ る.

2.4.2速 度 ・位 置 速 度 レ 釧Bニ[Vx"y1列B

亡=0の と き,Vxo=Oax。=0で あ り,そ の 後 の 時 間 経 過 で 加 速 度 は 計 測 で き る た め,

亡=sの と きV

xs=Vx(s̲1)+αxs{s‑(s‑1)}(2・22)

で あ り,tニ0か らy.x軸 に 関 し て も 同 様 に 各 軸 速 度 が 連 続 的 に 求 め られ,こ れ に よ り,カ プ セ ル の 飛 行 中 の 迎 角 α就

α・・1=lc・s‑・[y儲B](2・23)

が 求 め られ る.こ こ で,[κB]はB系 に お け る κ軸 方 向 座 標 の 単 位 ベ ク トル[XB]=[100]で あ る.

位 置[5劉Bニ[Sx3アSz]B

亡二 〇 の と き,Sxe二 初 期 値axeニ0で あ り,そ の 後 の 時 間 経 過 で 加 速 度 は 計 測 で き る た め,

亡 一 ・ の と きs

。、=・ 。(、一、)+v。,{s‑(・ 一 ・)}+ia.,{s‑(・ 一 ・)}・(2・24)

2.4.3空 力 中 心 を 原 点 と し た 重 心 位 置

空 力 中 心 を原 点 と した 重 心 位 置 を求 め る に あ た っ て,座 標 系 原 点 は機 体 背 面 中心 と し,使 用 す る 記 号 の 定 義 を 表2.6,そ れ ぞ れ の 中 心 位 置 お よび 迎 角 を 図2.7に 示 す.

表2.6空 力 中 心 を 原 点 と した 重 心 位 置 の 計 算 で使 用 す る記 号

記号 定義 求 め方

Xgc・YgClZgc

Xac、Yac,2「ac

XCG,YCG,ZCG

機体の重心位置 機体の空力申心位置

既 知 の値

姿 勢 角 か ら推 定 す る値

空 力 中 心 を 原 点 と し た 重 心 位 置Xg。‑Xac,Ygc‑Yac,Zgc‑Zac

AZB

y

yB K

gc

、(タ ＼ こ㌔̀L、

v ' ,

ac ὰ0藍

γF 1̀D1

1 ¹

1‑1 6

1工 ご6

XB

図2,7再 突 入 カ プ セル の 空 力 中 心 位 置

ま た,図2.7に お い て,

γ=短鴫) ^(2.25)

で あ る,機 体 の 空 力 中 心 位 置 κα。は,

x。c‑x、c+XcG(2・26)

さ ら に,

(2、27) ZCG=Zgc

以上よ り,空 力中心 を原点 とした重心位置は

x・ ・=

tan(一 髪1、 一 γ)一 一t、n(Zgcαtot十y)(228)

同 様 に して,YCGは

γ・・=tan(Zgc 一 β 一 γ)=‑tan誤 γ)(2・29)

従 っ て,(2.23),(2.25),(2.27),(2,28),(2、29)か ら,モ ー メ ン ト係 数(2.18)お よ び ダ ン ピ ン グ 係 数(2.19)が 求 め ら れ る.

第3章 カ ル マ ン フ ィル タ を用 い た 推 定

本 章 で は,観 測 値 の ノイ ズ が 同 時 に積 分 され る こ と を避 け る た め,カ ル マ ン フ ィル タ に よ る 誤 差 補 正 を 検 討 す る.フ ィル タ リ ン グ 設 定 手 法 を 示 し,観 測 方 程 式 ・状 態 方 程 式 の ノイ ズ設 定 後,共 分 散 行 列 か ら相 関 係 数 を 求 め る こ とで 誤 差 分 散 の 関係 を 検 証 す る.さ らに,誤 差 補 正 を して デ ー タ を 平 滑 化 す る こ とで ノ イ ズ を 減 少 させ,各 パ ラ メ ー タ の シ ミュ レー シ ョ ン結 果 を 考 察 す る.

3.1フ ィ ル タ リ ン グ 設 定

カ ル マ ン フ ィル タ とは,状 態 空 間 モ デ ル にお い て 「状 態 」を効 率 的 に 推 定 す る 計 算 ア ル ゴ リズ ム で あ る.

カ ル マ ン フ ィル タ を用 い て,「 観 測 値 」 と観 測 方 程 式,状 態 方 程 式 か ら 「状 態 」 を 推 定 す る こ とが で き る, ま た,カ ル マ ン フ ィル タ の 計 算 に用 い るパ ラ メ ー タ は事 前 に与 え る 必 要 が あ り,具 体 的 に は ノイ ズ の 大 き

さ(分 散)を 推 定 す る必 要 が あ る,

さ らに,状 態 空 間 モ デ ル を 推 定 す る た め に は, 1.パ ラ メ ー タ を 決 め る 方 法

2.パ ラ メ ー タ を使 っ て カ ル マ ン フ ィル タ を 実 行 し,「 状 態 」 を推 定 す る方 法 を検 討 す る 必 要 が あ る25).

3.1.1フ ィ ル タ リ ン グ 手 法 の 概 要

フ ィル タ リン グ手 法 の概 要 図 を 以 下 の 図3,1に 示 す

センサ

1囮 歴翫 引囮 噺■

1角速度センサ1 悔 卿91 ^, 脚1[石 ^ψ,θ,ψ

' 一一一

ヨ

; ^、

一̲̲一L、 一 一‑『

1一

●●,,● ●

x,y,z

1加速度セン司 ^ヨL轟

1

毛

■ ●,

κ,γ,z

一

x,y,z 図3.1フ ィル タ リ ン グ 手 法 の 概 要 図

上 記 流 れ に お い て 疑 似 観 測 量 を 導 入 し,フ ィル タ リ ン グ を進 め る, 疑 似 観 測 量 を 導 入 した カ ル マ ン フ ィル タ

ノイ ズ や拘 束 条 件 を観 測 して い るか の よ うに 人為 的 に作 り出 した観 測 量 を疑 似観 測 量 とす る.

i7(の ニf[t,x(亡),θ]+{systemnoise}

y(̀)ニh[亡,X(亡)1十{0正}seγ 「巳フ{ztionnò5θ}

上 記 の 推 定 問 題 を 考 え る と き,θ は 未 知 の パ ラ メ ー タ ベ ク トル で あ る,状 態 量 敏 δ と θ を 同 時 に 推 定/

同 定 す る に は 上 式 に 対 して フ ィ ル タ を 構 成 す る だ け で は θ を 精 度 よ く 同 定 で き な い.例 え ば θ が 一 定 値 で あ る な ら ば,そ れ を 一 っ の 状 態 量 と み て θ(∂ と し,

θ(t)ニ0 を シ ス テ ム 方 程 式 に 組 み 入 れ,拡 大 シ ス テ ム 方 程 式

兎(t)=f[t,x(t),θ]+{sys亡emnotse}

θ(t)={noise}

と 観 測 方 程 式 で カ ル マ ン フ ィル タ を 構 成 す る.更 に,拘 束 条 件

̀[x(の]ニ0 に 対 して

y,(t)≡0ニc[x(亡)](+{notse}) を 常 に 観 測 して い る と み な し て,観 測 モ デ ル に 付 加 し て い る,

γ(の ニh[t,X(の]+{observαtionnoise}

yp(t)≡0ニc[x(の](+{noise})

3.1.2カ ル マ ン フ ィ ル タ に 基 づ く 漸 化 式 カ ル マ ン フ ィ ル タ に お け る 状 態 方 程 式 は

Xt=F亡Xt̲̲エ 『一十一G亡 「v亡 (3.1)

観測方程式は

yt=HtXt十Wt (3.2)

で 与 え ら れ,Xt=状 態 ベ ク トル,yt=観 測 ベ ク トル,Ft=遷 移 行 列 で あ る.

カ ル マ ン フ ィル タ に 基 づ い て,予 測 の漸 化 式

x亡1亡一エ ニFtX亡 一1に一1

V、1、一、 ニF、Vt‑、1,.、F,T+6、q、G『

(3,3)

お よび,ブ イル タ の漸 化 式

κ亡 ・=Vtlt‑1H1(HtVtlt‑1H『+Rt)‑1

κ 亡1亡二Xtl亡 一一1+K亡(y亡 一H亡Xtl亡 一工)

レ』1亡二Vtlt‑1十K亡HtIltl亡 一1

(3.4)

が 与 え られ る,こ こ で,Vは 状 態 ベ ク トル の 分 散 共 分 散 行 列 で あ り,添 え 字 の'固 は 時 刻t‑1か らtへ の 予 測 分 布 を 示 す,こ の 結 果 か ら得 られ る状 態 ベ ク トル 銑 が 推 定 結 果 で あ る,

3.1.3姿 勢 角 の 推 定 に お け る 各 パ ラ メ ー タ の 定 義 式(2.11)か ら,

際卜∫/1砒一∫li

∫̀冗ψ 亡απ θ cosep

sinOP/cosθ 馨 鶴1[Ild亡 ^(3.5)

こ の 式 中 の 観 測ltl{であ る 角 速 度 を カ ル マ ン フ ィ ル タ に よ り 誤 ノ曽 甫正 す る.

カ ル マ ン フ ィ ル タ に お け る 状 態 方 程 式(3.1)お よ び 観 測 方 程 式(3.2)に お い て,各 パ ラ メ ー タ を

θψPqrl

二亡

ψ α θα ψ"y亡=

P q r

Ftニ

翫00100

nUnUイ⊥nUハUAU

nU4⊥∩UO∩り∩U 0配0010 00説001

(3.6)

と 設 定 す る.こ こ で,

[菱H藷難ll ^(3.7)

よ り,

θ。 一 ψ。 一 アαれ一些,ψ 。一 ταη一・ 一梅

Zq (3.8)

3.1.4速 度 ・位 置 の 推 定 に お け る 各 パ ラ メ ー タ の 定 義 速 度

観 測 され る加 速 度 につ い て 重 力 加 速 度 を考 慮 し,

團 一一(llyl]B‑r嘲

従 っ て,速 度 は

陰闇 一∫図   ∫C訂+[T]・・[lt]E)dt

で 求 め られ る,こ の 式 中 の 観 測 イ1}㌔:であ る 加 速 度 を カ ル マ ン フ ィル タ に よ り誤 差 補lllす る 、 カ ル マ ン フ ィ ル タ に お け る 状 態 方 程 式(3.Dお よ び 観 測 方 程 式(3.2)に お い て,各 パ ラ メ ー タ を

『・

著・ Xtニ 『9

惣 禦ρ

Za

子・

物 yt=㍗ 馳 )iF

Za

Ft=

00翫001

0翫0010

説00100

nU∩U4⊥ハUAUハ∪

ハU‑⊥AUOOハU

‑←nUnU凸UnUnU

と設 定 す る, 位 置

位 置 は速 度 と同様 に して,

[難1‑∫lvi]dt=∫[≧1砒+∬[矧砒=6t【箋1+叢1

で 求 め ら れ る.

二の 、Wの 観1期llLI:であ るJJI1速 度 を 一/」ル マ ン フ ィル タ にiこ り 、;呉差{{lilll{する 、

カ ル マ ン フ ィ ル タ に お け る 状 態 方 程 式(3,1)お よ び 観 測 方 程 式(32)に お い て,各 パ ラ メ ー タ を

(3.9)

(3.10)

(3,ll)

(3.12)

Xt=

皿皿"皿".",侃.侃κγZ.κ,y.Z.κ,V凶.Z

,二ytニ

α配α侃己配.乱.口・"κyZ.X.Ψ凶,Z.γ⁝,y,Z

,」F亡 二

2/亡♂00δ00100δ

00翫001000

0説0010000

説00100000

001占AUAUOOAUO

O10000000

100000000

00 δ亡2/20

0δt2/2 00 δ亡0

0δ 亡

00 10 01

(3,13)

と設 定 す る.

3.2推 定 パ ラ メ ー タ の 同 定 問 題

実 際 に 計 測 す る観 測 値 に は ノイ ズ が 含 まれ て お り,そ れ らを 除 去 す る た め に は カ ル マ ン フ ィル タ にお い て,3.1.2内 で 述 べ た シ ス テ ム ノイ ズ お よ び 観 測 ノ イ ズ の 設 定 お よび そ れ に伴 うカ ル マ ン ゲ イ ン の 決 定 が 非 常 に 重 要 で あ る,そ こで 本 節 で は,観 測 方 程 式 ・状 態 方 程 式 に 示 す ノイ ズ を最 尤 推 定 法 に よ っ て 設 定 し, 共 分 散 行 列 か ら相 関係 数 を 求 め る こ とで 誤 差 分 散 の 関係 を 示 す.特 に,位 置 の 推 定 に お い て は,観 測 値 で

あ る加 速 度 を2階 積 分 して 求 め るた め,1階 積 分 に よ っ て 求 め られ た結 果 に つ い て,カ ル マ ン ゲ イ ン を1 階 積 分 に お い て行 っ た 誤 差 補 正 よ り閾 値 を 大 幅 に 狭 く設 定 す る こ とで 精 度 を維 持 す る 必 要 が あ る.こ の よ

うに して 誤 差 補 正 を 行 い,デ ー タ を 平 滑 化 す る こ とで推 定 値 の 改 善 を検 討 す る,

3.2.1観 測 方 程 式 ・状態 方 程 式 に お け る ノイ ズ の 最 尤 推 定 値

最 尤 法 は,対 数 尤 度 を未 知 パ ラ メー タ の 関 数 と考 え て 対 数 尤 度 を最 大 にす る未 知 パ ラ メー タ の 値 を パ ラ メ ー タの 推 定 量 と して 採 用 す る推 定 手 法 で あ る,

観 測 方 程 式 お よび 状 態 方 程 式 にお け る ノイ ズ パ ラ メ ー タQ,,R亡 を 最 尤 法 に よ り推 定 す る,ま ず,初 期 状 態 分 布ai〜N(al,Pl)が 既 知 で あ る場 合,γ1,̲..,ytの 同 時 密 度 に つ い て,周 辺 密 度 と条 件 付 き 密 度 で

P(Y1,̲.,γ 亡)ニP(Yt̲エ)P(yt1Yt̲1)

と 分 解 で き る こ と を 用 い て,Yl,.,...,Ynの 同 時 密 度 関 数 は

Pω 一H・(y・1Y・ 一・)

亡=1 (3.14)

た だ し,P(yilYo)=P(Y1)

で 与 え ら れ る.こ こ で,p(ytlYt‑1)は 正 規 分 布N(αt,Ft)の 密 度 関 数 と な る た め,v亡=yt‑atの 密 度 関 数p(Vt) と

P(り 亡)=P(γ 亡1}を一1) (3.15)

の 関 係 に あ る.従 っ て,同 時 密 度p(y)の 対 数 で あ る対 数 尤 度 は

1・gL=1・9P(y)=一 卸 一1重(1・gFt+跨)

亡二1

(3,16)

で 与 え ら れ る,ま た,Ftは カ ル マ ン フ ィル タ か ら 得 ら れ る た め,カ ル マ ン フ ィ ル タ の 結 果 を 用 い て 尤 度 の 評 価 が 可 能 で あ る26).

式(3.16)の 対 数 尤 度 か ら 数 値 的 に ノ イ ズ パ ラ メ ー タq,,R亡 の 最 尤 推 定 量 を 求 め る に あ た っ て,分 散 は 正 の 実 数 に 限 定 され る と い う数 値 最 適 化 上 の 制 約 を 外 す た め に,ψvニlogqt,ψWニ10gRtに よ り全 て の 実 数 を と り う る パ ラ メ ー タ ψ,ψ に 変 換 し た 上 で シ ミ ュ レ ー シ ョ ン を 行 う.

vw

3.2.2観 測 値 の 平 滑 化

各 時 点 で観 測 され た ノ イ ズ の デ ー タ を 平 均 化 し,時 系 列 が 与 え られ た も とで の 状 態 α1,.,̲,αnの条 件 付 き 平 均 お よ び 分 布 を 求 め る,こ こ で,atを 平 滑化 状 態,Vtを 平 滑 化 状 態 分 散 と し,そ れ ら を 計 算 す る こ と で状 態 平 滑 化 が行 わ れ る.

平 滑 化 状 態 名 を フ ィル タ化 推 定 量 α雄 か ら導 出す る と,

δ・一α・1・+P・1・(舞1}+L・+・ 舞ll+…+L・+・ …Ln‑・ 鷺)

」P亡

L・=1一 κ・・ κ・=耳

(3.17)

で 与 え られ る.た だ し,Ptは 各11寺点tに 対 す る 状 態 α亡のVVI先 予 測atの 予 測 誤 差 分 散 で あ る.

式(3.17)の 括 弧 のill身 をrtと お く と,1期 前 のrt‑1に つ い て 後 ろ 向 き の 漸 化 式 が 得 ら れ る こ と か ら,平 滑 化 状 態 傷 は

rt‑、 ニFド1ツ 亡+ム 亡r亡,δtニ α 亡1亡+P亡1亡rt,亡 二 η,…,1・ (3.18)

の 状 態 平 滑 化 漸 化 式 に よっ て逐 次 的 に計 算 す る,た だ し,最 終 地 点 亡二 ηよ り後 に は 観 測 値 が 得 られ な い た め,rnニ0で あ る.

次 に,平 滑化 状 態 分 散Vtは 多 変 量 正 規 分 布 に お け る 条 件 付 き分 散 の 定 理 を用 い て,

UtニP亡lt‑P2tlt(Ft‑.1i+L曇+、1㌫ 亀+…+L曇+、 …L監 一、㌃1) (3」9)

で 与 え ら れ る.

式(3」9)の 括 弧 の 中 身 をNtと お く と,式(3,18)のrtの 漸 化 式 に つ い て 分 散 を と る こ と に よ り,1期 前 の 後 ろ 向 き 漸 化 式 ハ1亡一1二 耳+L蕃+11V亡 が 得 られ る こ と か ら,平  滑 化 状 態 分 散Vtは

N,.、 ニFr1+即 、,VtニP、1、‑P2、1,N、,亡 二n,…,1・ (320)

い た め,Nn=Oで あ る.

以 上 よ り,カ ル マ ン フ ィ ル タ の 成 果 物atlt,P亡lt,v亡 を 用 い て 平 滑 化 状 態 と 平 滑 化 状 態 分 散 を 求 め る 状 態 平 滑 化 の 逐 次 計 算 式 が 示 され た,

さ ら に,観 測 値 が 欠 測 し て い る と 仮 定 して 欠 測 あ り の 状 況 下 に お け る フ ィ ル タ リ ン グ お よ び 状 態 平 滑 化 を 行 う.欠 測 時 点tニ τ,τ+1,,、.tτ'にお け る カ ル マ ン フ ィ ル タ の 逐 次 計 算 式 が

α 亡i亡=α レ 」PtitニP亡

"ε+1=α 亡1亡'P̀+1=PtI亡 十R亡

(3.21)

で 与 え られ る.こ の と き,欠 測 期 間 が 複 数 に 分 か れ て い る 場 合 で も 状 態 平 滑 化 に よ っ て 自動 的 に 欠 測 時 点 t=τ,τ+1,̲tτ ‡の 状 態 が 補 間 さ れ る.

3.2.3カ ル マ ン ゲ イ ン の 決 定

カ ル マ ン ゲ イ ンKtの 決 定 に つ い て,観 測 値ytの 事 後 状 態 推 定 誤 差Xtlt一エと 直 行 す る 原 理 よ り,

E[x,1、.、y、.、‑1工=0 (3.22)

で あ る.従 っ て,事 後 状 態 誤 差 ベ ク トルXtltは,

X亡1亡 二Xti亡 一1+κ 亡(Vt‑HtXtlt‑、) (3.23)

以 上 よ り,式(3.22)お よ び(3.23)か ら,

Kt=Vtl亡 一1H『(HtVqt‑1H『+R亡)‑1 (3,24>

が 与 え られ,カ ル マ ン ゲ イ ン を決 定 す る,こ の た め に は 事 前 共 分 散 行 列Vtltを計 算 す る 必 要 が あ り,次 節 で 方 法 を 示 す.

3.2.4共 分 散 行 列 の 更 新

1期 前t‑1に お け る 事 後 共 分 散 行 列Vt‑"t‑1が 得 られ て い る と す る と,Vt‑11t‑1を 用 い て 時 刻tに お け る 事 前 共 分 散 行 列Vtlt‑1を 計 算 す る 予 測 ス テ ッ プ を 検 討 し,そ の 後Vtlt‑1を 用 い て 時 刻tに お け る 事 後 共 分 散 行 列1知 を 計 算 す る フ ィ ル タ リ ン グ ス テ ッ プ を 検 討 す る.

予 測 ス テ ッ プ は,

V亡it‑1=F亡Vt‑11亡 一1F1+σ 亡G亡G『 (3.25)