技 術 解 説

MATLAB によるクォータニオン数値計算
宍戸 幹夫＊
Mikio Shishido

クォータニオンは、航空機などでよく用いられるオイラー角のように角速度ベクトルとの関係が非線形で特異点をもつようなことはないため、飛行数値シミュレーションでは有用な回転表現である。
一般に、オイラー角の微分方程式を解く数値シミュレーションにおいては、特異点に落ち込んでしまい、数値解が発散してしまう。
クォータニオンについては昨年度技報で「クォータニオン計算便利ノート」として鎌倉事業部矢田部より報告がなされている。本書では、この資料をもとにクォータニオン表現による剛体の運動、およびオイラー表現による剛体の運動をそれぞれ実際にMATLABで実装、数値計算を行い、両回転表現の比較について報告する。
The quaternion is a useful expression for rotation used in flight simulation. The relation between thequaternion and its angular rate is linear and the quaternion doesn t have a singular point like Euler
angles. Euler angles are often used to describe aircraft dynamics.
Mr. Yatabe s paper

Convenient note for calculating Quaternion

appeared in last year s MSS

Technical Report. In his paper, Mr. Yatabe summarized useful equations for attitude calculation using
the quaternion. My paper compares quaternion and Euler angle expressions for rigid body dynamics.
MATLAB is used to formulate and simulate rigid body dynamics as follow :1．まえがきそこで、本報告では、クォータニオン表現による回転表現を実際に数値計算が行えるようMATLABで実装お力学では、位置や姿勢角を時間微分することにより速よび計算結果の確認を行い、そのプログラムソースコー

度、角速度を得られるとは限らない。これは、位置や姿ドを記載する。本報告による実装は、飛行機、ロケット、勢角を定義する座標系の移動、回転を考慮して時間微分人工衛星等の運動シミュレーションや慣性装置、ジンバを行わなければならないためである。このようなことをル機構を有するハードウエア等のモデリングとシミュレ避けるため、位置ベクトルについては、座標系の移動、ーションに適用できると考えられる。

回転を考慮しなくて良い慣性空間上で定義されることが多い。

本報告は、まず取り扱う座標系と方向余弦行列の定義を行い、次に、オイラー角、クォータニオンの定式化にしかし、姿勢角についてはこのようにはならず、航空ついて述べる。この定式化をもとに数値解析ツールであ機などで用いられるオイラー角では角速度ベクトルとのるMATLABで実装を行った後に数値計算結果を示す。

関係が非線形で特異点をもつ。そのため、オイラー角表
現による飛行運動シミュレーションにおいては、上記特異点近傍で姿勢角が発散してしまう。

2．座標系と方向余弦行列の定義
本資料では、以下の2 つの座標系を取り扱う。

オイラー角に変わる回転表現としてクォータニオンが・慣性座標系ある。クォータニオンは、回転軸と回転角よりなる４成・機体座標系分で３次元空間の回転を表現する。クォータニオンでは、慣性座標系は、慣性空間に固定された座標系である。

オイラー角による回転表現で生じるような特異点が存在

機体座標系は、航空機等の機体の重心位置に固定され、しない。そのため、航空機等の姿勢角表現に有効な方法である。

＊鎌倉事業部 第三技術部

慣性空間を移動、回転する座標系である。
両座標系とも右手系の直交座標系とする。機体座標系

MSS技報・Vol.19

44

技 術 解 説
3.2

では、 軸を機体前方にとり、 軸を下方にとることとする。機体座標系を図１に示す。航空力学では一般に機方向余弦行列からオイラー角表現

式（2）より以下のように方向余弦行列を用いてオイラー角を定義することができる。

体座標系の 軸をロール軸、 軸をピッチ軸、 軸をヨー軸と呼ぶ。

（3a）方向余弦行列とは、任意の２つの座標系の関係を変換する行列である。本資料では、慣性座標系から機体座標系への変換行列を方向余弦行列
性座標系の基底ベクトルを底ベクトルを（3b）とする。ここで、慣、機体座標系の基（3c）とすると、以下のような関係
3.3で表される。

オイラー角の微分方程式

機体座標系での各軸の角速度をそれぞれ ， ， とすると、式（2）よりオイラー角の運動方程式は以下のよ（1）

うに表わされる。
（4a）3．オイラー角の定式化
オイラー角は各軸まわりの回転角の組として定義され、（4b）方向余弦行列と同様に座標系の関係を表現する要素である。

（4c）

方向余弦行列が９つのパラメータからなるため剛体の回転を直感的に理解することは難しいが、オイラー角は上式のとおり、オイラー角の微分方程式はピッチ角３軸の回転角の組であるため回転を理解することは容易の±90degに特異点が存在する。これが、オイラー角表である。

現による姿勢角計算の欠点である。オイラー角による姿
勢角計算を行う際は、特異点をさけるため、回転順序を3.1変えてオイラー角を定義する等の工夫が必要である。

オイラー角表現による方向余弦行列

慣性座標を ， ， 軸の順でそれぞれ ， ， 回転させたとき機体座標となる角をオイラー角とする。航空力学では ， ， 軸まわりの回転角をロール角
ピッチ角 、ヨー角 と呼ぶ。
ここで、方向余弦行列はオイラー角を用い以下の、4．クォータニオンの定式化
クォータニオンは、回転軸と回転角よりなる４成分で３次元空間の回転を表現する。クォータニオンでは、オイラー角による回転表現で生じるような特異点が存在しない。そのため、飛行機等の姿勢角表現に有効な方法でように表わされる。

ある。
クォータニオンについて、「クォータニオン計算便利ノート」3）を参考に述べる。
（2）

4.1

クォータニオンの定義

クォータニオン は、回転方向の単位ベクトル と回転角 の４成分で定義される。クォータニオンを図２にy軸示す。
q（5）
pnx軸
θrz軸

図１45

機体座標系の定義

図２クォータニオンここで、クォータニオンの定義から明らかなとおりノルムは１である。

要な統合開発環境を提供するアプリケーションソフトウエアである。さらに、制御系設計、信号処理、画像処理など各種分野に対して利便性の高い数値演算ライブラリ（6）

がToolboxと呼ばれるMATLABを拡張するパッケージの形で提供されている。また、プログラミング言語として、効率的な数値解析プログラム開発を行う事も可能で

4.2

ある。

クォータニオン表現による方向余弦行列

クォータニオン表現による慣性座標系から機体座標系への方向余弦行列

MATLABの特徴を以下に簡単に述べる。
・データアクセスの容易性

は、以下のとおりである。

データの基本構造は行列であり、変数の型宣言や配列の大きさを宣言する必要がない。また行列を直接演算する事ができ、行列の一要素ごと演算する必要がな（7）

い。
・対話型操作環境とプログラミング機能

4.3

統合開発環境により、コマンド実行、ファイル操作

方向余弦行列からクォータニオン表現

方向余弦行列からクォータニオンを求めるには、式（7）およびヘルプの参照などのオペレーションを対話形式を用いる。まず、方向余弦行列の対角項からクォータニで行える。また、MATLABはインタプリタ型式の言オンの１変数を求め、次に、求めた変数を用いてその他語であるため、コンパイル、リンクといったプログラの変数を求める。対角項からムの作成および変更から実行するまでの手間が少ない。

を求めたときは、以下・豊富な数値計算関数とグラフィックス関数のとおりとなる。
（8a）
（8b）
（8c）
（8d）
ただし、方向余弦行列の対角項からクォータニオンの１変数を選ぶ際は、変数がゼロにならないようなものに

MATLABのみで800以上の関数、コマンドが提供されている。拡張パッケージであるToolboxにより、さらに関数、コマンドが提供されている。また、豊富なグラフィックス関数を用いることにより、容易に計算結果を可視化することができる。
5.2

MATLABで作成した関数一覧

数値計算を行うにあたって、MATLABで実際に作成した関数の一覧を表１に示す。
また、クォータニオン、オイラー角、方向余弦行列の関係および作成した関数の関係を図３に示す。図３で示しなければならない。

されるように、方向余弦行列を介して、クォータニオン4.4とオイラー角の変換が可能である。

クォータニオンの微分方程式

クォータニオンの時間微分は、角速度ベクトル

次節以降で作成した各関数についての概要を説明する。

を用いて以下のように表される。
5.3

方向余弦行列計算

オイラー角から方向余弦行列を計算する関数（9）euler2dcm のソースコードを図４に示す。この関数の表１

上式のとおり、クォータニオンの微分方程式には、オイラー角では式（4）のように存在する特異点を持たない。
5．MATLABによる実装
5.1

MATLABとは

MATLABは、データ解析、モデリング、シミュレーション、可視化などの数値解析、アルゴリズム開発に必

関数

関数一覧
処理内容

dcm＝euler2dcm（euler） オイラー角から方向余弦行列を計算
dcm＝quat2dcm（q）

クォータニオンから方向余弦行列を計算 図５

euler＝dcm2euler（dcm） 方向余弦行列からオイラー角を計算
q＝dcm2quat（dcm）

図番号
図４図６

方向余弦行列からクォータニオンを計算 図７dx＝euler_eq（t，x）

オイラー角の微分方程式を計算

図８dx＝quat_eq（t，x）

クォータニオンの微分方程式を計算

図９simdemo

オイラー角とクォータニオンの比較計算 図10

MSS技報・Vol.1946

技 術 解 説
引数はオイラー角３要素をもつベクトルeulerであり、各クォータニオン変数を算出し、その中で最大となる変戻り値は方向余弦行列dcmである。計算式は、式（2）で数を選んで、その他の変数を再計算している。

ある。
次に、クォータニオンから方向余弦行列を計算する関5.5

微分方程式計算数quat2dcmのソースコードを図５に示す。この関数のオイラー角の微分方程式を計算する関数euler_eqのソ

引数はクォータニオン４要素をもつベクトルquatであースコードを図８に示す。次に、クォータニオンの微分り、戻り値は方向余弦行列dcmである。計算式は、式（7）方程式を計算する関数quat_eqのソースコードを図９にである。

示す。また、この微分方程式を実際に解く関数が図10で示されるsimdemoである。この関数は、MATLAB組み

5.4方向余弦行列からオイラー角、クォータニオン計算方向余弦行列からオイラー角から計算する関数dcm2込み関数である常微分方程式のソルバode23を使用して両微分方程式を解いている。

eulerのソースコードを図６に、クォータニオンから方向余弦行列を計算する関数dcm2quatのソースコードを図７に示す。これらは、式（3)、式（8）を実装したものである。ただし、dcm2quat 関数については、対角項から

quat2dcm
クォータニオン

47dcm2euler 関数

図７dcm2quat 関数

図８euler_eq 関数

quat2euler
方向余弦行列dcm2quat

図３

図６

オイラー角
euler2dcm

クォータニオン、オイラー角および方向余弦行列の関係

図４eular2dcm 関数

図５quat2dcm 関数

表２シミュレーション条件図９quat_eq 関数

ロール角φ［deg］

条件１

条件２

初期ロール姿勢角−30deg

０deg

初期ピッチ姿勢角−20deg

80deg

初期ヨー姿勢角−10deg

０deg

ロール角速度

５deg/s０deg/s

ピッチ角速度

10deg/s

５deg/s

ヨー角速度

15deg/s

０deg/s100500

クォータニオン表現による計算
オイラー表現による計算

ピッチ角θ［deg］−500

1234567891001234567891001234567891060
40
200−20ヨー角Ψ［deg］

60
40200−20時間【s】

図11

条件１計算結果：オイラー角q110−101234

5678

9100

12345678910012345678910012345678910

図10

simdemo 関数

q210−1q310−16．数値計算結果
数値計算は、表２に示す２条件で行った。条件１として、特異点をとらない場合、オイラー角、クォータニオ

q410−1ンによる両手法が一致することを確認する条件を選んだ。
次に、条件２として、オイラー角が特異点をとる条件を時間【s】

図12

条件１計算結果：クォータニオン選んだ。
まず、条件１での数値計算結果を図11、12に示す。
図11で示すクォータニオン表現による姿勢角は、クォーが90degになる条件である。そのため、図13からわかるタニオンを方向余弦行列に変換し、さらにオイラー角にようにオイラー角による計算では、2sにロール姿勢角と変換したものである。条件１は、オイラー角が特異点にヨー姿勢角が発散している。これに対してクォータニオよる発散がない条件であるため、図11からわかるようにンによる計算では、図14で示されるように発散せず良好オイラー角、クォータニオンによる計算結果が一致する。

な結果である。クォータニオンから変換したピッチ姿勢オイラー角に変換する前のクォータニオンの各パラメー

角 が2s以降減少しているのは、方向余弦行列からオイ

タは図12を参照されたい。

ラー角に変換する際に式（3）のとおりピッチ姿勢角の範条件２での数値計算結果を図13、14に示す。条件２は、時間2sでオイラー角の特異点であるピッチ姿勢角囲が−90deg90degとなるためである。ただし、クォータニオン自体に発散はないため問題ではない。

MSS技報・Vol.1948

200

1100q1

ロール角φ［deg］

技 術 解 説0
クォータニオン表現による計算
オイラー表現による計算−100−200
123456789q2200

123456789100123456789100123456789100123456789100100−101−100−2000123456789102000−110010q4

ヨー角Ψ［deg］

0110q3ピッチ角θ［deg］

00−1−100−20001234567890−110

時間【s】

時間【s】

図13条件２計算結果：オイラー角

図14条件２計算結果：クォータニオン

7．むすび
剛体の回転表現であるオイラー角とクォータニオンに

参考文献
盧James R. Wertz: Spacecraft Attitude Determi

ついて、定式化とMATLABによる実装および数値計算

nation and Control, Kluwer Academic Publishers，

を行った。

1978.

数値計算結果から、特異点をとらない飛行条件ではオ盪

イラー角とクォータニオンの両手法の結果が一致することを確認した。また、特異点をとる飛行条件においては、クォータニオンの優位性を確認した。これにより、クォータニオンが飛行数値シミュレーションで有用な回転表現であることを示せた。

49

加藤寛一朗：航空機力学入門，東京大学出版会，
1982.

蘯矢田部学：クォータニオン計算便利ノート，MSS
技報vol.18，2007.