飛行機が「30°右に旋回してから20°上昇した」と聞けば、誰でもおおよその姿勢を想像できます。このように3つの角度で物体の向きを表すのが オイラー角 です。回転行列の9要素に比べてたった3パラメータで姿勢を記述でき、物理的な直感とも一致するため、航空宇宙工学で広く使われています。
しかし、オイラー角には致命的な弱点があります。ある特定の姿勢に達すると、3つの回転軸のうち2つが重なってしまい、1つの自由度が失われる ジンバルロック という現象が発生するのです。この問題は、1969年のアポロ11号のミッションでも実際に懸念されたことで知られています。
オイラー角を理解すると、以下の分野で役立ちます。
- 航空機の姿勢記述: ロール・ピッチ・ヨーという3つの角度で機体姿勢を直感的に表現できる
- ロボット工学: 関節角度の表現やエンドエフェクタの姿勢指定に使われる
- 衛星の姿勢制御: 地球指向衛星の姿勢角の微小変動を記述する際に利用される
- ジンバルロックの回避: 問題の本質を理解することで、クォータニオンなどの代替手法が必要な場面を判断できる
本記事の内容
- オイラー角の定義と複数の回転順序
- 回転行列とオイラー角の対応関係
- ジンバルロックの幾何学的説明と数学的原因
- 航空機・衛星での使い分け
- Pythonによるオイラー角とジンバルロックの可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
オイラー角とは
回転行列は3×3の行列（9要素）で姿勢を表しますが、直交条件 $C^TC = I$ が6つの拘束を与えるため、独立なパラメータは3つしかありません。「ならば、最初から3つのパラメータだけで姿勢を表現できないか」という発想は自然です。
オイラー角は、3次元の任意の回転を、座標軸まわりの3回の連続回転の組み合わせで表す という考え方です。18世紀の数学者レオンハルト・オイラーにちなんで名付けられました。
直感的に言えば、首を動かして任意の方向を向くとき、「左右に回す」「上下に回す」「首を傾ける」という3つの動作を組み合わせるのと同じです。この3つの角度がオイラー角に対応します。
回転順序の分類
3回の連続回転をどの軸順で行うかによって、複数のオイラー角の定義が存在します。大きく2つのカテゴリに分類されます。
固有オイラー角（Proper Euler angles）: 最初と最後の回転軸が同じもの。例えば $z$-$x$-$z$（3-1-3）、$z$-$y$-$z$（3-2-3）など6通りあります。ジャイロスコープの記述や天体力学でよく使われます。
テイト-ブライアン角（Tait-Bryan angles）: 3つの回転軸がすべて異なるもの。例えば $z$-$y$-$x$（3-2-1）、$x$-$y$-$z$（1-2-3）など6通りあります。航空工学ではロール・ピッチ・ヨーとしておなじみです。
合計12通りの回転順序が存在しますが、どれを使っても任意の回転を表現できます。以下では、航空宇宙工学で最もよく使われる 3-2-1回転（$z$-$y$-$x$ 順）を中心に解説します。
この多様な回転順序の中から、航空宇宙で最も一般的な3-2-1回転を詳しく見ていきましょう。
3-2-1 オイラー角（ヨー・ピッチ・ロール）
定義
3-2-1オイラー角は、以下の3回の連続回転で姿勢を記述します。
- ヨー（$\psi$）: $z$ 軸まわりに角度 $\psi$ 回転
- ピッチ（$\theta$）: 新しい $y$ 軸（$y’$）まわりに角度 $\theta$ 回転
- ロール（$\phi$）: 新しい $x$ 軸（$x”$）まわりに角度 $\phi$ 回転
航空機で考えると、
- ヨー $\psi$: 機首を左右に振る動き（方位角の変化）
- ピッチ $\theta$: 機首を上下に振る動き（仰角の変化）
- ロール $\phi$: 翼を傾ける動き（バンク角の変化）
に対応します。パイロットが「右30°旋回、機首20°上げ、左バンク10°」と報告すれば、$\psi=30°$、$\theta=20°$、$\phi=-10°$ と表されます。
回転行列との対応
前記事で導出した基本回転行列を使うと、3-2-1回転に対応する全体の回転行列は次のようになります。
$$ \bm{C}_{321}(\phi, \theta, \psi) = R_x(\phi) \, R_y(\theta) \, R_z(\psi) $$
各基本回転行列を代入します。
$$ R_z(\psi) = \begin{pmatrix} c_\psi & -s_\psi & 0 \\ s_\psi & c_\psi & 0 \\ 0 & 0 & 1 \end{pmatrix}, \quad R_y(\theta) = \begin{pmatrix} c_\theta & 0 & s_\theta \\ 0 & 1 & 0 \\ -s_\theta & 0 & c_\theta \end{pmatrix}, \quad R_x(\phi) = \begin{pmatrix} 1 & 0 & 0 \\ 0 & c_\phi & -s_\phi \\ 0 & s_\phi & c_\phi \end{pmatrix} $$
ここで $c_\alpha = \cos\alpha$, $s_\alpha = \sin\alpha$ と略記しています。
まず $R_y(\theta) R_z(\psi)$ を計算します。
$$ R_y R_z = \begin{pmatrix} c_\theta c_\psi & -c_\theta s_\psi & s_\theta \\ s_\psi & c_\psi & 0 \\ -s_\theta c_\psi & s_\theta s_\psi & c_\theta \end{pmatrix} $$
次に、この結果に左から $R_x(\phi)$ を掛けます。$(2,1)$ 成分を例にとると、
$$ (R_x)_{21} \cdot (R_y R_z)_{11} + (R_x)_{22} \cdot (R_y R_z)_{21} + (R_x)_{23} \cdot (R_y R_z)_{31} $$
$= 0 \cdot c_\theta c_\psi + c_\phi \cdot s_\psi + (-s_\phi)(-s_\theta c_\psi)$ $= c_\phi s_\psi + s_\phi s_\theta c_\psi$
同様にすべての成分を計算すると、最終的な回転行列が得られます。
$$ \bm{C}_{321} = \begin{pmatrix} c_\theta c_\psi & -c_\theta s_\psi & s_\theta \\ s_\phi s_\theta c_\psi + c_\phi s_\psi & -s_\phi s_\theta s_\psi + c_\phi c_\psi & -s_\phi c_\theta \\ -c_\phi s_\theta c_\psi + s_\phi s_\psi & c_\phi s_\theta s_\psi + s_\phi c_\psi & c_\phi c_\theta \end{pmatrix} $$
この式は重要なので、しっかり記憶しておく価値があります。
回転行列からオイラー角への逆変換
回転行列 $\bm{C}$ が与えられたとき、オイラー角を逆算する方法を導出しましょう。$\bm{C}_{321}$ の各成分を $C_{ij}$ と書きます。
$(1,3)$ 成分に着目すると、
$$ C_{13} = \sin\theta $$
よって、ピッチ角は直接得られます。
$$ \theta = \arcsin(C_{13}) $$
$\cos\theta \neq 0$（すなわち $\theta \neq \pm 90°$）のとき、$(2,3)$ と $(3,3)$ 成分の比からロール角が求まります。
$$ \frac{C_{23}}{C_{33}} = \frac{-\sin\phi \cos\theta}{\cos\phi \cos\theta} = -\tan\phi $$
したがって、
$$ \phi = \arctan2(-C_{23}, C_{33}) $$
同様に、$(1,2)$ と $(1,1)$ 成分の比からヨー角が求まります。
$$ \frac{C_{12}}{C_{11}} = \frac{-\cos\theta \sin\psi}{\cos\theta \cos\psi} = -\tan\psi $$
$$ \psi = \arctan2(-C_{12}, C_{11}) $$
ここで $\arctan2(y, x)$ は4象限逆正接関数で、$x$ と $y$ の符号から正しい象限を判定します。通常の $\arctan(y/x)$ では象限が不定になりますが、$\arctan2$ ならば $[-\pi, \pi]$ の範囲で一意に角度を求められます。
ただし $\cos\theta = 0$（$\theta = \pm 90°$）の場合、分母がゼロになるためこの逆変換は使えません。これがジンバルロックの数学的な兆候です。
3-2-1回転以外のオイラー角も同様の構造を持ちますが、特異点の位置が異なります。次に、他の回転順序も見ておきましょう。
その他の回転順序
3-1-3 オイラー角（古典オイラー角）
天体力学や回転する剛体の運動で伝統的に使われる回転順序です。
- $z$ 軸まわりに $\alpha$ 回転
- 新しい $x’$ 軸まわりに $\beta$ 回転
- 新しい $z”$ 軸まわりに $\gamma$ 回転
対応する回転行列は、
$$ \bm{C}_{313} = R_z(\gamma) \, R_x(\beta) \, R_z(\alpha) $$
各成分を展開すると、
$$ \bm{C}_{313} = \begin{pmatrix} c_\alpha c_\gamma – s_\alpha c_\beta s_\gamma & -s_\alpha c_\gamma – c_\alpha c_\beta s_\gamma & s_\beta s_\gamma \\ c_\alpha s_\gamma + s_\alpha c_\beta c_\gamma & -s_\alpha s_\gamma + c_\alpha c_\beta c_\gamma & -s_\beta c_\gamma \\ s_\alpha s_\beta & c_\alpha s_\beta & c_\beta \end{pmatrix} $$
この場合、特異点は $\beta = 0$ または $\beta = \pi$（$\sin\beta = 0$）のときに発生します。$\beta = 0$ は回転軸が重なる姿勢に対応します。
回転順序の選び方
12通りの回転順序のうち、どれを使うかは応用分野によって慣習が異なります。
| 回転順序 | 表記 | 主な応用分野 | 特異点 |
|---|---|---|---|
| 3-2-1 | $z$-$y$-$x$ | 航空工学（ヨー・ピッチ・ロール） | $\theta = \pm 90°$ |
| 3-1-3 | $z$-$x$-$z$ | 天体力学、剛体のオイラー運動 | $\beta = 0, \pi$ |
| 1-2-3 | $x$-$y$-$z$ | ロボット工学の一部 | $\theta = \pm 90°$ |
| 3-2-3 | $z$-$y$-$z$ | 量子力学（ウィグナーD行列） | $\beta = 0, \pi$ |
重要なのは、どの回転順序を選んでも、必ずどこかに特異点が存在する ということです。これは3パラメータで SO(3) の全体を覆おうとする限り避けられない数学的な制限です（「毛の生えた球は梳かせない」というトポロジーの定理に関連します）。
各回転順序にはそれぞれ特異点がありますが、その物理的意味を理解するためにジンバルロックを詳しく見ていきましょう。
ジンバルロックの幾何学的理解
ジンバルとは
ジンバルロックの「ジンバル」は、コンパスやジャイロスコープに使われる回転支持機構です。3つのリング（外側ジンバル、中間ジンバル、内側ジンバル）が、それぞれ異なる軸で回転できるように入れ子になっています。
3つのジンバルが異なる方向を向いているとき、内側の物体は3次元空間のどの方向にも向けることができます。3つの回転自由度が独立に機能しているからです。
軸が重なるとき
ところが、中間ジンバルを $90°$ 回転させると、外側ジンバルの回転軸と内側ジンバルの回転軸が平行になってしまいます。2つのジンバルが同じ軸を共有するため、見かけ上1つの自由度が失われ、3軸の独立な回転が2軸分の回転にしかなりません。
具体的に、3-2-1オイラー角でピッチ角 $\theta = 90°$ の場合を考えましょう。
$\theta = 90°$ を $\bm{C}_{321}$ に代入すると、$\cos\theta = 0$、$\sin\theta = 1$ なので、
$$ \bm{C}_{321}\big|_{\theta=90°} = \begin{pmatrix} 0 & 0 & 1 \\ s_\phi c_\psi + c_\phi s_\psi & -s_\phi s_\psi + c_\phi c_\psi & 0 \\ -c_\phi c_\psi + s_\phi s_\psi & c_\phi s_\psi + s_\phi c_\psi & 0 \end{pmatrix} $$
ここで三角関数の加法定理 $\sin(\phi+\psi) = s_\phi c_\psi + c_\phi s_\psi$、$\cos(\phi+\psi) = c_\phi c_\psi – s_\phi s_\psi$ を使うと、
$$ \bm{C}_{321}\big|_{\theta=90°} = \begin{pmatrix} 0 & 0 & 1 \\ \sin(\phi+\psi) & \cos(\phi+\psi) & 0 \\ -\cos(\phi+\psi) & \sin(\phi+\psi) & 0 \end{pmatrix} $$
驚くべきことに、$\phi$ と $\psi$ が個別には現れず、和 $\phi + \psi$ のみで回転行列が決まります。$\phi = 30°, \psi = 20°$ と $\phi = 40°, \psi = 10°$ は同じ回転行列を与えます。つまり、$\phi$ と $\psi$ を個別に決定することが不可能になり、1つの自由度が実質的に失われるのです。
ジンバルロックの直感的なまとめ
ジンバルロックを一言で表すと、「3つの回転パラメータのうち2つが同じ効果を持ってしまい、独立に制御できなくなる状態」 です。
車のステアリングに例えるなら、前輪を真横に向けてしまうと（$90°$ステアリング）、アクセルを踏んでもステアリングを切っても車は同じ方向にしか動かない — そんな状況に似ています。
この現象は数学的にどう記述できるのでしょうか。次に、ヤコビアン行列を通じてジンバルロックを解析的に理解しましょう。
ジンバルロックの数学的原因 — ヤコビアンの特異性
角速度とオイラー角の関係
機体の角速度ベクトル $\bm{\omega} = (\omega_x, \omega_y, \omega_z)^T$ は、オイラー角の時間微分 $\dot{\phi}, \dot{\theta}, \dot{\psi}$ と次の関係で結ばれます。
$$ \bm{\omega} = \bm{S}(\phi, \theta) \begin{pmatrix} \dot{\phi} \\ \dot{\theta} \\ \dot{\psi} \end{pmatrix} $$
ここで $\bm{S}$ は角速度変換行列で、3-2-1オイラー角の場合は以下のようになります。
なぜこの行列が必要かを直感的に理解しましょう。3つの回転の軸はそれぞれ別の座標系に属しています。$\dot{\psi}$ は元の $z$ 軸、$\dot{\theta}$ は1回回転後の $y’$ 軸、$\dot{\phi}$ は2回回転後の $x”$ 軸まわりの角速度です。これらを同一の座標系（機体座標系）で表すために変換行列が必要なのです。
機体座標系での角速度を計算するため、各軸まわりの回転レートを機体座標系に変換します。
$\dot{\phi}$ は $x”$ 軸（= 機体の $x$ 軸）まわりなので、機体座標系ではそのまま $(1, 0, 0)^T \dot{\phi}$ です。
$\dot{\theta}$ は $y’$ 軸まわりです。機体座標系に変換するには $R_x(\phi)$ で変換する必要があります。
$$ R_x(\phi) \begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix} \dot{\theta} = \begin{pmatrix} 0 \\ \cos\phi \\ \sin\phi \end{pmatrix} \dot{\theta} $$
$\dot{\psi}$ は元の $z$ 軸まわりです。機体座標系に変換するには $R_x(\phi) R_y(\theta)$ で変換します。
$$ R_x(\phi) R_y(\theta) \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix} \dot{\psi} = \begin{pmatrix} \sin\theta \\ -\sin\phi \cos\theta \\ \cos\phi \cos\theta \end{pmatrix} \dot{\psi} $$
これらを合わせると、角速度変換行列は次のようになります。
$$ \bm{S} = \begin{pmatrix} 1 & 0 & \sin\theta \\ 0 & \cos\phi & -\sin\phi\cos\theta \\ 0 & \sin\phi & \cos\phi\cos\theta \end{pmatrix} $$
ヤコビアンの行列式と特異性
この変換行列の行列式を計算しましょう。第1行で余因子展開すると、
$$ \det(\bm{S}) = 1 \cdot (\cos\phi \cdot \cos\phi\cos\theta + \sin\phi \cdot \sin\phi\cos\theta) + \sin\theta \cdot (0) $$
$\cos^2\phi + \sin^2\phi = 1$ を使うと、
$$ \det(\bm{S}) = \cos\theta $$
したがって、$\theta = \pm 90°$ のとき $\det(\bm{S}) = 0$ となり、行列 $\bm{S}$ は特異（逆行列が存在しない）になります。
これは何を意味するでしょうか。$\bm{S}$ が特異ということは、角速度 $\bm{\omega}$ が与えられてもオイラー角レート $(\dot{\phi}, \dot{\theta}, \dot{\psi})^T$ を一意に求められないということです。制御の観点からは、望みの角速度を実現するためのオイラー角の変化率が計算できないため、姿勢制御が破綻します。
もう少し具体的に言えば、$\theta = 90°$ のとき、
$$ \bm{S}\big|_{\theta=90°} = \begin{pmatrix} 1 & 0 & 1 \\ 0 & \cos\phi & 0 \\ 0 & \sin\phi & 0 \end{pmatrix} $$
第1列と第3列が同一になっています。$\dot{\phi}$ を増やす効果と $\dot{\psi}$ を増やす効果が同じになり、独立に制御できないことがわかります。
特異点は避けられない
オイラー角のジンバルロックは、3パラメータで SO(3) を覆おうとする限り原理的に避けられません。位相幾何学の観点から、SO(3) は3次元球面と局所的に同相な多様体ですが、3次元球面を1つの座標系で覆うことはできません（ちょうど地球表面の全ての点を単一の経度・緯度で表そうとすると北極と南極で問題が生じるのと同じです）。
この本質的な限界を克服するために、4パラメータを使うクォータニオン表現が導入されます。しかし、その話は次の記事に譲り、まずはPythonでオイラー角とジンバルロックを可視化してみましょう。
Pythonでの実装
オイラー角から回転行列への変換
まず、3-2-1オイラー角と回転行列の変換を実装します。
import numpy as np
def euler_to_dcm_321(phi, theta, psi):
"""3-2-1オイラー角(ロール, ピッチ, ヨー)から回転行列(DCM)を計算"""
# 基本回転行列
Rx = np.array([
[1, 0, 0],
[0, np.cos(phi), -np.sin(phi)],
[0, np.sin(phi), np.cos(phi)]
])
Ry = np.array([
[ np.cos(theta), 0, np.sin(theta)],
[ 0, 1, 0 ],
[-np.sin(theta), 0, np.cos(theta)]
])
Rz = np.array([
[np.cos(psi), -np.sin(psi), 0],
[np.sin(psi), np.cos(psi), 0],
[0, 0, 1]
])
return Rx @ Ry @ Rz
def dcm_to_euler_321(C):
"""回転行列(DCM)から3-2-1オイラー角を逆算"""
theta = np.arcsin(np.clip(C[0, 2], -1.0, 1.0))
if np.abs(np.cos(theta)) > 1e-10:
phi = np.arctan2(-C[1, 2], C[2, 2])
psi = np.arctan2(-C[0, 1], C[0, 0])
else:
# ジンバルロック: phi+psiのみ決定可能
phi = 0.0 # phiを0に固定（慣習）
if C[0, 2] > 0:
psi = np.arctan2(C[1, 0], C[1, 1])
else:
psi = -np.arctan2(C[1, 0], C[1, 1])
return phi, theta, psi
# 動作確認
phi, theta, psi = np.radians(30), np.radians(45), np.radians(60)
C = euler_to_dcm_321(phi, theta, psi)
phi_r, theta_r, psi_r = dcm_to_euler_321(C)
print("元のオイラー角 [deg]:", np.degrees([phi, theta, psi]))
print("逆変換結果 [deg]: ", np.degrees([phi_r, theta_r, psi_r]))
print("差のノルム [deg]: ", np.linalg.norm(
np.degrees(np.array([phi, theta, psi]) - np.array([phi_r, theta_r, psi_r]))))
元のオイラー角と逆変換結果が一致していることから、順変換と逆変換が正しく実装できていることが確認できます。差のノルムはマシンイプシロン程度に収まります。ジンバルロック付近では dcm_to_euler_321 が $\phi = 0$ に固定する分岐に入るため、逆変換の結果が元の角度と異なる場合があります。
オイラー角による姿勢変化の可視化
3つのオイラー角を順に適用したときの姿勢変化を可視化します。
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
def euler_to_dcm_321(phi, theta, psi):
Rx = np.array([[1,0,0],[0,np.cos(phi),-np.sin(phi)],[0,np.sin(phi),np.cos(phi)]])
Ry = np.array([[np.cos(theta),0,np.sin(theta)],[0,1,0],[-np.sin(theta),0,np.cos(theta)]])
Rz = np.array([[np.cos(psi),-np.sin(psi),0],[np.sin(psi),np.cos(psi),0],[0,0,1]])
return Rx @ Ry @ Rz
def plot_frame(ax, R, origin=np.zeros(3), alpha=1.0, length=1.0):
colors = ['r', 'g', 'b']
labels = ['x', 'y', 'z']
for i in range(3):
d = R[:, i] * length
ax.quiver(*origin, *d, color=colors[i],
arrow_length_ratio=0.1, alpha=alpha, linewidth=2.5)
ax.text(*(origin + d*1.15), labels[i], color=colors[i],
fontsize=11, alpha=alpha)
phi, theta, psi = np.radians(30), np.radians(45), np.radians(60)
# 段階的に回転を適用
Rz_only = euler_to_dcm_321(0, 0, psi)
RyRz = euler_to_dcm_321(0, theta, psi)
RxRyRz = euler_to_dcm_321(phi, theta, psi)
fig = plt.figure(figsize=(18, 5))
titles = [
f'Step 1: Yaw ψ={np.degrees(psi):.0f}°',
f'Step 2: + Pitch θ={np.degrees(theta):.0f}°',
f'Step 3: + Roll φ={np.degrees(phi):.0f}°'
]
rotations = [Rz_only, RyRz, RxRyRz]
for idx, (R, title) in enumerate(zip(rotations, titles)):
ax = fig.add_subplot(1, 3, idx+1, projection='3d')
plot_frame(ax, np.eye(3), alpha=0.2)
plot_frame(ax, R)
ax.set_title(title, fontsize=12)
ax.set_xlim([-1.3,1.3]); ax.set_ylim([-1.3,1.3]); ax.set_zlim([-1.3,1.3])
ax.set_xlabel('X'); ax.set_ylabel('Y'); ax.set_zlabel('Z')
plt.tight_layout()
plt.savefig('euler_angles_steps.png', dpi=150, bbox_inches='tight')
plt.show()
上のグラフでは、3-2-1回転の各段階での座標系の変化を表示しています。
- Step 1（ヨー60°）: $z$ 軸まわりに60°回転。$x$, $y$ 軸が $z$ 軸を中心に回転し、$z$ 軸は不変です。
- Step 2（+ ピッチ45°）: さらに新しい $y’$ 軸まわりに45°回転。$x$ 軸が上を向き始め、$z$ 軸が下を向きます。
- Step 3（+ ロール30°）: 最後に新しい $x”$ 軸まわりに30°回転。$y$ 軸と $z$ 軸がさらに回転します。
3つの段階を経て、最終的な姿勢に至る過程が視覚的にわかります。
ジンバルロックの可視化
ジンバルロックの核心部分を可視化します。ピッチ角 $\theta$ を $0°$ から $90°$ まで変化させたときの、角速度変換行列のヤコビアン行列式の変化を調べます。
import numpy as np
import matplotlib.pyplot as plt
# ヤコビアンの行列式 = cos(theta)
theta_range = np.linspace(-np.pi/2, np.pi/2, 500)
det_S = np.cos(theta_range)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# (a) det(S)の変化
axes[0].plot(np.degrees(theta_range), det_S, 'b-', linewidth=2)
axes[0].axhline(y=0, color='r', linestyle='--', alpha=0.7, label='det(S) = 0')
axes[0].axvline(x=90, color='r', linestyle=':', alpha=0.5)
axes[0].axvline(x=-90, color='r', linestyle=':', alpha=0.5)
axes[0].fill_between(np.degrees(theta_range), det_S, 0,
where=(np.abs(det_S) < 0.1), color='red', alpha=0.15,
label='Gimbal lock zone')
axes[0].set_xlabel('Pitch angle θ [deg]', fontsize=12)
axes[0].set_ylabel('det(S)', fontsize=12)
axes[0].set_title('Jacobian Determinant vs Pitch Angle', fontsize=13)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
# (b) 条件数（逆行列の安定性指標）
phi = np.radians(30)
cond_numbers = []
for th in theta_range:
S = np.array([
[1, 0, np.sin(th)],
[0, np.cos(phi), -np.sin(phi)*np.cos(th)],
[0, np.sin(phi), np.cos(phi)*np.cos(th)]
])
cond_numbers.append(np.linalg.cond(S))
axes[1].semilogy(np.degrees(theta_range), cond_numbers, 'b-', linewidth=2)
axes[1].axvline(x=90, color='r', linestyle=':', alpha=0.5, label='θ = ±90°')
axes[1].axvline(x=-90, color='r', linestyle=':', alpha=0.5)
axes[1].set_xlabel('Pitch angle θ [deg]', fontsize=12)
axes[1].set_ylabel('Condition number (log scale)', fontsize=12)
axes[1].set_title('Condition Number of S Matrix', fontsize=13)
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('gimbal_lock_analysis.png', dpi=150, bbox_inches='tight')
plt.show()
上の2つのグラフから、ジンバルロックの数学的性質が明確に読み取れます。
- 左図（ヤコビアン行列式）: $\det(\bm{S}) = \cos\theta$ が $\theta = \pm 90°$ でゼロに達します。赤い領域はジンバルロックの影響が大きい「危険ゾーン」を示しています。ゼロに近づくほど角速度からオイラー角レートへの変換が不安定になります。
- 右図（条件数）: 条件数は行列の「逆行列の安定性」を示す指標で、$\theta = \pm 90°$ に近づくと急激に増大（対数スケールで直線的に上昇）します。条件数が大きいほど数値計算が不安定になり、わずかな角速度の誤差がオイラー角レートの巨大な誤差に増幅されます。
ジンバルロック時の自由度損失の可視化
ジンバルロック状態で、ロール角 $\phi$ とヨー角 $\psi$ を個別に変えても同じ効果しか生まないことを可視化します。
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
def euler_to_dcm_321(phi, theta, psi):
Rx = np.array([[1,0,0],[0,np.cos(phi),-np.sin(phi)],[0,np.sin(phi),np.cos(phi)]])
Ry = np.array([[np.cos(theta),0,np.sin(theta)],[0,1,0],[-np.sin(theta),0,np.cos(theta)]])
Rz = np.array([[np.cos(psi),-np.sin(psi),0],[np.sin(psi),np.cos(psi),0],[0,0,1]])
return Rx @ Ry @ Rz
def plot_frame(ax, R, origin=np.zeros(3), alpha=1.0, length=1.0):
colors = ['r', 'g', 'b']
labels = ['x', 'y', 'z']
for i in range(3):
d = R[:, i] * length
ax.quiver(*origin, *d, color=colors[i],
arrow_length_ratio=0.1, alpha=alpha, linewidth=2.5)
theta_lock = np.radians(90)
fig = plt.figure(figsize=(16, 5))
# 異なる(phi, psi)の組で同じ合計角の場合
combos = [
(0, 60, 'φ=0°, ψ=60°'),
(20, 40, 'φ=20°, ψ=40°'),
(40, 20, 'φ=40°, ψ=20°'),
(60, 0, 'φ=60°, ψ=0°')
]
for idx, (phi_d, psi_d, label) in enumerate(combos):
ax = fig.add_subplot(1, 4, idx+1, projection='3d')
phi, psi = np.radians(phi_d), np.radians(psi_d)
C = euler_to_dcm_321(phi, theta_lock, psi)
plot_frame(ax, np.eye(3), alpha=0.2)
plot_frame(ax, C)
ax.set_title(f'{label}\nφ+ψ={phi_d+psi_d}°', fontsize=11)
ax.set_xlim([-1.3,1.3]); ax.set_ylim([-1.3,1.3]); ax.set_zlim([-1.3,1.3])
ax.set_xlabel('X'); ax.set_ylabel('Y'); ax.set_zlabel('Z')
ax.view_init(elev=25, azim=45)
plt.suptitle('Gimbal Lock (θ=90°): Different (φ,ψ) with same φ+ψ',
fontsize=14, y=1.02)
plt.tight_layout()
plt.savefig('gimbal_lock_freedom_loss.png', dpi=150, bbox_inches='tight')
plt.show()
この結果は非常に印象的です。4つのパネルすべてで、最終的な座標系の向きが完全に同じです。$\phi$ と $\psi$ の値はそれぞれ異なる（0°/60°、20°/40°、40°/20°、60°/0°）のに、$\phi + \psi = 60°$ が同じであるため、回転行列は全く同一になります。
これがジンバルロックの本質です。$\theta = 90°$ のとき、$\phi$ と $\psi$ は独立なパラメータではなくなり、回転行列に影響するのは和 $\phi + \psi$ だけです。3つのパラメータで3自由度を持つはずが、実質2自由度に退化してしまうのです。
航空機姿勢の3D可視化
最後に、航空機モデルの姿勢をオイラー角で表現する応用例を示します。
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
def euler_to_dcm_321(phi, theta, psi):
Rx = np.array([[1,0,0],[0,np.cos(phi),-np.sin(phi)],[0,np.sin(phi),np.cos(phi)]])
Ry = np.array([[np.cos(theta),0,np.sin(theta)],[0,1,0],[-np.sin(theta),0,np.cos(theta)]])
Rz = np.array([[np.cos(psi),-np.sin(psi),0],[np.sin(psi),np.cos(psi),0],[0,0,1]])
return Rx @ Ry @ Rz
# 簡易航空機の形状定義
def aircraft_vertices():
"""簡易的な航空機の頂点を返す"""
# 胴体
fuselage = np.array([
[1.0, 0, 0], [-0.5, 0, 0], # 機首〜尾部
[0.5, 0.08, 0.05], [0.5, -0.08, 0.05],
[0.5, 0.08, -0.05], [0.5, -0.08, -0.05],
[-0.3, 0.06, 0.04], [-0.3, -0.06, 0.04],
[-0.3, 0.06, -0.04], [-0.3, -0.06, -0.04],
])
# 主翼
wing = np.array([
[0.2, 0, 0], [0.0, 0, 0],
[0.2, 0.8, 0], [0.0, 0.8, 0],
[0.2, -0.8, 0], [0.0, -0.8, 0],
])
# 尾翼
tail = np.array([
[-0.3, 0, 0], [-0.5, 0, 0],
[-0.3, 0.3, 0], [-0.5, 0.3, 0],
[-0.3, -0.3, 0], [-0.5, -0.3, 0],
[-0.3, 0, 0.25], [-0.5, 0, 0.25],
])
return fuselage, wing, tail
fig = plt.figure(figsize=(16, 5))
cases = [
(0, 0, 30, 'Yaw 30° (heading change)'),
(0, 20, 0, 'Pitch 20° (climb)'),
(25, 0, 30, 'Yaw 30° + Roll 25° (banked turn)')
]
for idx, (phi_d, theta_d, psi_d, title) in enumerate(cases):
ax = fig.add_subplot(1, 3, idx+1, projection='3d')
C = euler_to_dcm_321(np.radians(phi_d), np.radians(theta_d), np.radians(psi_d))
# 主翼を描画（三角形のポリゴンで表現）
wing_pts = np.array([
[0.2, 0, 0], [0.0, 0.8, 0], [0.0, -0.8, 0], # 主翼（三角形）
])
tail_h = np.array([
[-0.3, 0, 0], [-0.5, 0.3, 0], [-0.5, -0.3, 0], # 水平尾翼
])
tail_v = np.array([
[-0.3, 0, 0], [-0.5, 0, 0], [-0.5, 0, 0.25], # 垂直尾翼
])
body = np.array([
[1.0, 0, 0], [0.2, 0.06, 0.04], [0.2, -0.06, 0.04],
[0.2, 0.06, -0.04], [0.2, -0.06, -0.04],
])
# 回転を適用
wing_r = (C @ wing_pts.T).T
tail_h_r = (C @ tail_h.T).T
tail_v_r = (C @ tail_v.T).T
body_r = (C @ body.T).T
# 描画
ax.add_collection3d(Poly3DCollection(
[wing_r], alpha=0.6, facecolor='#2196F3', edgecolor='black'))
ax.add_collection3d(Poly3DCollection(
[tail_h_r], alpha=0.6, facecolor='#FF9800', edgecolor='black'))
ax.add_collection3d(Poly3DCollection(
[tail_v_r], alpha=0.6, facecolor='#4CAF50', edgecolor='black'))
# 機首方向を矢印で表示
nose = C @ np.array([1.0, 0, 0])
ax.quiver(0, 0, 0, *nose, color='red', arrow_length_ratio=0.08, linewidth=2)
ax.set_title(title, fontsize=11)
ax.set_xlim([-1, 1]); ax.set_ylim([-1, 1]); ax.set_zlim([-0.5, 0.5])
ax.set_xlabel('X (forward)'); ax.set_ylabel('Y (right)')
ax.set_zlabel('Z (down)')
ax.view_init(elev=20, azim=-60)
plt.tight_layout()
plt.savefig('aircraft_euler.png', dpi=150, bbox_inches='tight')
plt.show()
上の3つのパネルは、航空機の典型的な姿勢をオイラー角で表現したものです。
- 左図（ヨー30°）: 機首が右に30°旋回しています。水平飛行のまま方位だけが変わっています。
- 中図（ピッチ20°）: 機首が上に20°傾いた上昇姿勢です。翼の傾きはありません。
- 右図（ヨー30° + ロール25°）: 旋回しつつバンク角25°で傾いた状態です。実際の航空機はバンクをつけて旋回するため、この姿勢が現実的です。
オイラー角は、このように航空機の姿勢を人間にとって直感的にわかりやすい形で表現できます。通常の飛行ではピッチ角が $\pm 90°$ に達することはめったにないため、ジンバルロックは実用上の問題になりにくいのです。
航空機と衛星での使い分け
航空機の場合
航空機では、3-2-1オイラー角（ヨー・ピッチ・ロール）が標準的に使われます。理由は明確です。
- 通常の飛行では $\theta$（ピッチ角）が $\pm 30°$ 程度に収まり、ジンバルロックの $\pm 90°$ から十分に離れている
- ヨー・ピッチ・ロールの3パラメータが操縦桿の操作と直感的に対応する
- フライトデータの記録・解析・通信で人間が読みやすい
ただし、戦闘機の急激な機動（ピッチ $90°$ に近づくコブラ機動など）では問題が生じる可能性があります。
衛星の場合
人工衛星では状況が異なります。
地球指向衛星（常に地球を向く衛星）の場合、姿勢の微小変動をオイラー角で記述するのは便利です。正常運用時のロール・ピッチ・ヨーは数度以内に収まるため、ジンバルロックの心配はありません。線形化した姿勢制御則（PD制御など）も書きやすいです。
一方、大角度の姿勢変更（太陽追尾モードへの切り替え、セーフモードからの復帰など）では、オイラー角の使用は危険です。ピッチが $90°$ を通過する可能性がある場合、ジンバルロックを避けるために クォータニオン を使います。
| 場面 | 推奨手法 | 理由 |
|---|---|---|
| 航空機の通常飛行 | 3-2-1 オイラー角 | 直感的、ピッチが小さい |
| 地球指向衛星の微小姿勢制御 | オイラー角（線形化） | 制御則が書きやすい |
| 衛星の大角度姿勢変更 | クォータニオン | 特異点回避 |
| ジャイロスコープ | 3-1-3 オイラー角 | 回転対称性と相性が良い |
| ゲーム・CG | クォータニオン | SLERP補間が滑らか |
まとめ
本記事では、オイラー角による姿勢表現とその限界であるジンバルロックについて解説しました。
- オイラー角 は3つの回転角で姿勢を表す方法で、12通りの回転順序がある。航空宇宙では3-2-1（ヨー・ピッチ・ロール）が代表的
- 回転行列との対応: $\bm{C}_{321} = R_x(\phi) R_y(\theta) R_z(\psi)$ で展開でき、逆変換には $\arctan2$ を使う
- ジンバルロック: 特定の姿勢（3-2-1では $\theta = \pm 90°$）で2つの回転軸が重なり、1自由度が失われる
- 数学的原因: 角速度変換行列のヤコビアン $\det(\bm{S}) = \cos\theta$ がゼロになり、逆変換が不可能になる
- 原理的限界: 3パラメータで SO(3) 全体を特異点なしに覆うことはできない
ジンバルロックは3パラメータ表現の本質的な限界であり、回転順序を変えても特異点の位置が変わるだけです。この限界を根本的に解決するのが、4パラメータを使う クォータニオン（四元数） 表現です。
次のステップとして、以下の記事も参考にしてください。