剛体の姿勢表現 — 座標変換と回転行列の基礎

人工衛星が宇宙空間で正しくアンテナを地上に向けたり、太陽電池パネルを太陽に向けたりするためには、衛星が「今どの方向を向いているか」を正確に記述できなければなりません。この「向き」を数学的に扱うのが 姿勢表現 の問題です。

しかし、「向き」とは一体何でしょうか。日常生活では「右を向いている」「上を向いている」と直感的に言いますが、3次元空間で物体の向きを曖昧さなく記述するには、厳密な数学的枠組みが必要です。

姿勢表現を理解すると、以下のような幅広い応用先に直結します。

  • 人工衛星の姿勢制御: 軌道上で太陽電池パネルやアンテナの指向方向を管理する
  • ロボットアームの関節制御: 各関節の回転を正確に記述し、先端の位置・姿勢を計算する
  • コンピュータグラフィックス: 3Dオブジェクトの回転・描画に回転行列が不可欠
  • 航空機のフライトコントロール: 機体の姿勢をリアルタイムで把握し制御する

本記事の内容

  • 座標系の定義と姿勢の概念
  • 方向余弦行列(DCM)の定義と性質
  • 基本回転行列 $R_x$, $R_y$, $R_z$ の導出
  • 合成回転の計算
  • Pythonでの回転行列による座標変換の可視化

前提知識

この記事を読む前に、以下の知識があると理解がスムーズです。

  • 線形代数の基礎(行列の積、転置、逆行列、行列式)
  • 三角関数の基本的な性質

座標系とは何か

3次元空間で物体の位置や向きを記述するには、まず「基準」が必要です。地図でいえば「北」や「東」にあたるもの、つまり 座標系 です。

座標系とは、空間内の1点(原点)と、互いに直交する3本の単位ベクトル(基底ベクトル)の組のことです。右手系の直交座標系では、3つの基底ベクトル $\bm{e}_1, \bm{e}_2, \bm{e}_3$ が以下を満たします。

$$ \bm{e}_i \cdot \bm{e}_j = \delta_{ij}, \quad \bm{e}_1 \times \bm{e}_2 = \bm{e}_3 $$

ここで $\delta_{ij}$ はクロネッカーのデルタです。

慣性座標系と機体座標系

衛星の姿勢を考える際に、少なくとも2つの座標系が登場します。

慣性座標系(基準座標系) は、空間に固定された座標系です。宇宙工学ではECI(地球中心慣性座標系)がよく使われます。地球の中心を原点とし、春分点方向が $x$ 軸、地球の自転軸方向が $z$ 軸です。

機体座標系 は、衛星や航空機の本体に固定された座標系です。衛星とともに回転するため、衛星の姿勢を直接反映します。通常、衛星の主要な軸(例えばアンテナの指向方向、太陽電池パネルの法線方向)に合わせて定義します。

「姿勢」とは、慣性座標系から見た機体座標系の向き、すなわち 2つの座標系の間の回転関係 を意味します。この回転関係を数学的に表現するのが、次に説明する方向余弦行列です。

ここまでで座標系と姿勢の概念を確認しました。では、2つの座標系の関係をどのように数式で表すのでしょうか。その答えが 方向余弦行列(Direction Cosine Matrix: DCM) です。

方向余弦行列(DCM)の定義

教室で先生が「北を向きなさい」と言ったとき、あなたの「前・右・上」の3方向は東西南北と一定の角度関係を持ちます。方向余弦行列は、この角度関係を完全に記録する3×3の行列です。

慣性座標系 $\mathcal{A}$ の基底ベクトルを $\{\bm{a}_1, \bm{a}_2, \bm{a}_3\}$、機体座標系 $\mathcal{B}$ の基底ベクトルを $\{\bm{b}_1, \bm{b}_2, \bm{b}_3\}$ とします。方向余弦行列 $\bm{C}$ の各成分は、2つの座標系の基底ベクトル同士の内積(余弦)で定義されます。

$$ C_{ij} = \bm{b}_i \cdot \bm{a}_j = \cos\theta_{ij} $$

ここで $\theta_{ij}$ は $\bm{b}_i$ と $\bm{a}_j$ のなす角です。行列として書くと次のようになります。

$$ \bm{C} = \begin{pmatrix} \bm{b}_1 \cdot \bm{a}_1 & \bm{b}_1 \cdot \bm{a}_2 & \bm{b}_1 \cdot \bm{a}_3 \\ \bm{b}_2 \cdot \bm{a}_1 & \bm{b}_2 \cdot \bm{a}_2 & \bm{b}_2 \cdot \bm{a}_3 \\ \bm{b}_3 \cdot \bm{a}_1 & \bm{b}_3 \cdot \bm{a}_2 & \bm{b}_3 \cdot \bm{a}_3 \end{pmatrix} $$

この行列 $\bm{C}$ を使うと、座標系 $\mathcal{A}$ で表されたベクトル $\bm{r}_{\mathcal{A}}$ を座標系 $\mathcal{B}$ の成分に変換できます。

$$ \bm{r}_{\mathcal{B}} = \bm{C} \, \bm{r}_{\mathcal{A}} $$

逆に、$\mathcal{B}$ から $\mathcal{A}$ への変換は $\bm{C}$ の逆行列を使います。この逆行列が簡単に求まることは、次のセクションで確認する直交性から導かれます。

DCMが「方向余弦」行列と呼ばれるのは、各成分が2つの軸のなす角の余弦値だからです。9つの成分がありますが、独立なパラメータは3つだけです(直交条件6つで拘束されるため)。これは、3次元の回転が3つの自由度を持つことと整合しています。

方向余弦行列を使えばベクトルの座標変換ができることがわかりました。では、この行列にはどのような数学的性質があるのでしょうか。回転行列の本質に迫りましょう。

回転行列の性質

回転行列と呼ばれるためには、単に3×3の行列であるだけでは不十分です。物体の回転は「形を変えない変換」であり、ベクトルの長さや角度を保存しなければなりません。この物理的要請が、回転行列の数学的性質を決定します。

直交性

回転行列 $\bm{C}$ は 直交行列 です。すなわち、

$$ \bm{C}^T \bm{C} = \bm{C} \bm{C}^T = \bm{I} $$

が成り立ちます。ここで $\bm{I}$ は3×3の単位行列です。

この性質を確認してみましょう。$\bm{C}^T \bm{C}$ の $(i,j)$ 成分を計算すると、

$$ (\bm{C}^T \bm{C})_{ij} = \sum_{k=1}^{3} C_{ki} C_{kj} = \sum_{k=1}^{3} (\bm{b}_k \cdot \bm{a}_i)(\bm{b}_k \cdot \bm{a}_j) $$

$\{\bm{b}_k\}$ が正規直交基底であることから、この和は $\bm{a}_i \cdot \bm{a}_j = \delta_{ij}$ に等しくなります。よって $\bm{C}^T \bm{C} = \bm{I}$ が示されました。

直交性の重要な帰結として、逆行列が転置行列に等しい ことが挙げられます。

$$ \bm{C}^{-1} = \bm{C}^T $$

逆行列の計算は一般に計算コストが高いですが、直交行列なら単に転置するだけで済みます。これは数値計算において非常に便利な性質です。

行列式が1

回転行列の行列式は必ず +1 です。

$$ \det(\bm{C}) = +1 $$

直交行列の行列式は $\pm 1$ のいずれかですが、$-1$ の場合は鏡映(反転)を含む変換になります。座標系の右手性を保つ純粋な回転では、行列式は常に +1 です。

なぜ行列式が +1 なのかを直感的に理解するには、回転の連続性を考えます。回転角0(恒等変換)のとき $\bm{C} = \bm{I}$ で $\det(\bm{C}) = 1$ です。回転を連続的に変化させると行列式も連続的に変化しますが、$\det$ は $\pm 1$ しか取れないため、連続性により常に +1 のままです。

ベクトルの大きさの保存

直交行列によるベクトル変換は、ベクトルのノルム(長さ)を保存します。

$$ \|\bm{C}\bm{r}\|^2 = (\bm{C}\bm{r})^T(\bm{C}\bm{r}) = \bm{r}^T \bm{C}^T \bm{C} \bm{r} = \bm{r}^T \bm{r} = \|\bm{r}\|^2 $$

ここで、$\bm{C}^T \bm{C} = \bm{I}$(直交性)を使いました。回転は向きだけを変え、大きさを変えない変換であることが数学的に裏付けられています。

回転群 SO(3)

以上の性質を満たす行列の集合 — 直交かつ行列式が1の3×3実行列 — を 特殊直交群 SO(3) と呼びます。

$$ SO(3) = \{ \bm{C} \in \mathbb{R}^{3 \times 3} \mid \bm{C}^T \bm{C} = \bm{I}, \; \det(\bm{C}) = +1 \} $$

SO(3) は「群」としての構造を持ちます。つまり、

  • 閉性: 2つの回転行列の積も回転行列(2回の回転を合成しても回転)
  • 結合則: $(\bm{C}_1 \bm{C}_2) \bm{C}_3 = \bm{C}_1 (\bm{C}_2 \bm{C}_3)$
  • 単位元: $\bm{I} \in SO(3)$(回転しない変換)
  • 逆元: $\bm{C} \in SO(3) \Rightarrow \bm{C}^{-1} = \bm{C}^T \in SO(3)$(逆回転も回転)

この群構造は、回転の合成や逆回転が自然に定義できることを保証しています。

回転行列の一般的な性質を理解しました。次に、3次元回転の最も基本的な構成要素である「単一の座標軸まわりの回転」を具体的に導出しましょう。

基本回転行列の導出

3次元の任意の回転は、座標軸まわりの回転の組み合わせで表現できます。ここでは、$x$ 軸、$y$ 軸、$z$ 軸それぞれの周りの回転行列を導出します。

$z$ 軸まわりの回転 $R_z(\psi)$

まず最も直感的な $z$ 軸まわりの回転から始めましょう。紙の上で点を原点中心に回す操作がこれに対応します。

$z$ 軸まわりに角度 $\psi$ だけ回転すると、$x$-$y$ 平面内で回転が起こり、$z$ 成分は変化しません。元の座標 $(x, y, z)$ が新しい座標 $(x’, y’, z’)$ に変換されます。

$x$-$y$ 平面での2次元回転を考えます。点 $(x, y)$ を角度 $\psi$ だけ反時計回りに回転すると、三角関数の加法定理から次が得られます。

$$ \begin{align} x’ &= x \cos\psi – y \sin\psi \\ y’ &= x \sin\psi + y \cos\psi \end{align} $$

$z$ 成分は変わらないので $z’ = z$ です。これを行列形式で書くと、

$$ R_z(\psi) = \begin{pmatrix} \cos\psi & -\sin\psi & 0 \\ \sin\psi & \cos\psi & 0 \\ 0 & 0 & 1 \end{pmatrix} $$

第3行・第3列が $(0, 0, 1)$ であることは、$z$ 軸方向が不変であることを反映しています。

$x$ 軸まわりの回転 $R_x(\phi)$

同じ考え方を $x$ 軸まわりに適用します。今度は $y$-$z$ 平面内で回転が起こり、$x$ 成分は不変です。

$y$-$z$ 平面で角度 $\phi$ だけ回転すると、

$$ \begin{align} y’ &= y \cos\phi – z \sin\phi \\ z’ &= y \sin\phi + z \cos\phi \end{align} $$

$x$ 成分は $x’ = x$ のまま変化しないので、行列形式は次のようになります。

$$ R_x(\phi) = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \cos\phi & -\sin\phi \\ 0 & \sin\phi & \cos\phi \end{pmatrix} $$

$y$ 軸まわりの回転 $R_y(\theta)$

$y$ 軸まわりの回転では、$z$-$x$ 平面内で回転が起こります。ここで注意が必要なのは、右手系の規約に従うと $z$ から $x$ への回転が正方向となることです。

$$ \begin{align} z’ &= z \cos\theta – x \sin\theta \\ x’ &= z \sin\theta + x \cos\theta \end{align} $$

$(x’, y’, z’)$ の順に整理して行列形式で書くと、

$$ R_y(\theta) = \begin{pmatrix} \cos\theta & 0 & \sin\theta \\ 0 & 1 & 0 \\ -\sin\theta & 0 & \cos\theta \end{pmatrix} $$

$R_y$ では $(1,3)$ 成分が $+\sin\theta$ で $(3,1)$ 成分が $-\sin\theta$ となっています。これは $R_x$ や $R_z$ とは符号の位置が異なるため、注意が必要です。右手系の規約で $z \to x$ が正回転方向であることから来ています。

基本回転行列の性質の検証

導出した3つの基本回転行列が、先ほど確認した回転行列の性質を満たしているか確認しましょう。$R_z(\psi)$ を例にとります。

直交性の確認:

$R_z^T R_z$ を計算すると、

$$ R_z^T R_z = \begin{pmatrix} \cos\psi & \sin\psi & 0 \\ -\sin\psi & \cos\psi & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} \cos\psi & -\sin\psi & 0 \\ \sin\psi & \cos\psi & 0 \\ 0 & 0 & 1 \end{pmatrix} $$

$(1,1)$ 成分は $\cos^2\psi + \sin^2\psi = 1$、$(1,2)$ 成分は $-\cos\psi\sin\psi + \sin\psi\cos\psi = 0$ となり、確かに $\bm{I}$ が得られます。

行列式の確認:

$R_z$ の行列式を第3行で余因子展開すると、

$$ \det(R_z) = 1 \cdot (\cos^2\psi + \sin^2\psi) = 1 $$

基本回転行列はいずれも回転行列の定義を満たすことが確認できました。

3つの基本回転行列を手に入れました。しかし、実際の回転は1つの座標軸まわりだけとは限りません。次に、複数の回転を組み合わせる「合成回転」について見ていきましょう。

合成回転

回転の合成は行列の積

物体がまず $z$ 軸まわりに $\psi$ 回転し、次に $y$ 軸まわりに $\theta$ 回転し、最後に $x$ 軸まわりに $\phi$ 回転する場合、全体の回転行列 $\bm{C}$ は各回転行列の積で表されます。

$$ \bm{C} = R_x(\phi) \, R_y(\theta) \, R_z(\psi) $$

ここで重要なのは、行列の積の順序です。行列の掛け算は一般に交換法則が成り立ちません($\bm{A}\bm{B} \neq \bm{B}\bm{A}$)。回転の物理的意味でも、「先に右を向いてから上を向く」と「先に上を向いてから右を向く」では、最終的な姿勢は異なります。

固定軸回転と回転軸回転

合成回転には2つの解釈があります。

固定軸回転(外因性回転): 回転の基準となる軸が空間に固定されたまま、各回転を順に適用します。ベクトル $\bm{r}$ に対して、まず $R_z$、次に $R_y$、最後に $R_x$ を掛けます。

$$ \bm{r}’ = R_x(\phi) \, R_y(\theta) \, R_z(\psi) \, \bm{r} $$

この場合、行列は 右から左 に適用される順序で書きます。

回転軸回転(内因性回転): 各回転の軸が、直前の回転で動いた後の座標系に固定されています。物体に取り付けた座標軸が回転とともに動いていくイメージです。興味深いことに、内因性回転を $z’$-$y’$-$x’$ の順に適用した場合の結果は、外因性回転を $x$-$y$-$z$ の順に適用した場合と同じ行列になります。

合成回転の具体例

$z$ 軸まわりに $90°$、次に $x$ 軸まわりに $90°$ 回転する場合を考えます。

まず $R_z(90°)$ を計算します。

$$ R_z(90°) = \begin{pmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{pmatrix} $$

次に $R_x(90°)$ を計算します。

$$ R_x(90°) = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 0 & -1 \\ 0 & 1 & 0 \end{pmatrix} $$

合成回転は $\bm{C} = R_x(90°) R_z(90°)$ なので、

$$ \bm{C} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 0 & -1 \\ 0 & 1 & 0 \end{pmatrix} \begin{pmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{pmatrix} = \begin{pmatrix} 0 & -1 & 0 \\ 0 & 0 & -1 \\ 1 & 0 & 0 \end{pmatrix} $$

この結果が直感に合うか確認しましょう。元の $x$ 軸方向のベクトル $(1,0,0)^T$ に $\bm{C}$ を掛けると $(0,0,1)^T$ つまり $z$ 方向になります。まず $z$ 軸回転で $x \to y$ 方向に移り、次に $x$ 軸回転で $y \to z$ 方向に移ったことと整合しています。

順序を入れ替えた場合、$R_z(90°) R_x(90°)$ を計算してみると、

$$ R_z(90°) R_x(90°) = \begin{pmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 & 0 \\ 0 & 0 & -1 \\ 0 & 1 & 0 \end{pmatrix} = \begin{pmatrix} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{pmatrix} $$

明らかに先ほどの結果と異なります。回転の順序が結果に影響することが、数値で確認できました。

合成回転の理論を理解したところで、これをPythonで実装し、回転行列の効果を視覚的に確認しましょう。

Pythonでの実装と可視化

ここでは、回転行列の定義をPythonで実装し、3Dプロットで座標変換の効果を視覚的に確認します。

基本回転行列の実装

まず3つの基本回転行列を関数として定義します。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def Rx(phi):
    """x軸まわりの回転行列"""
    c, s = np.cos(phi), np.sin(phi)
    return np.array([
        [1,  0,  0],
        [0,  c, -s],
        [0,  s,  c]
    ])

def Ry(theta):
    """y軸まわりの回転行列"""
    c, s = np.cos(theta), np.sin(theta)
    return np.array([
        [ c, 0, s],
        [ 0, 1, 0],
        [-s, 0, c]
    ])

def Rz(psi):
    """z軸まわりの回転行列"""
    c, s = np.cos(psi), np.sin(psi)
    return np.array([
        [c, -s, 0],
        [s,  c, 0],
        [0,  0, 1]
    ])

# 直交性の検証
phi = np.radians(30)
C = Rx(phi)
print("C^T @ C =")
print(np.round(C.T @ C, 10))
print(f"det(C) = {np.linalg.det(C):.10f}")

このコードで生成される出力から、$R_x(30°)$ の転置と自身の積が単位行列になること、および行列式が正確に1であることが確認できます。数値誤差は $10^{-16}$ 程度のマシンイプシロンのオーダーに収まっており、理論通りの性質を満たしています。

座標系の回転の可視化

次に、基本回転行列による座標系の回転を3Dで可視化します。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def Rx(phi):
    c, s = np.cos(phi), np.sin(phi)
    return np.array([[1,0,0],[0,c,-s],[0,s,c]])

def Ry(theta):
    c, s = np.cos(theta), np.sin(theta)
    return np.array([[c,0,s],[0,1,0],[-s,0,c]])

def Rz(psi):
    c, s = np.cos(psi), np.sin(psi)
    return np.array([[c,-s,0],[s,c,0],[0,0,1]])

def plot_frame(ax, R, origin=np.zeros(3), label_prefix='',
               colors=('r','g','b'), alpha=1.0, length=1.0):
    """座標系の3軸を矢印で描画"""
    labels = ['x', 'y', 'z']
    for i in range(3):
        direction = R[:, i] * length
        ax.quiver(*origin, *direction, color=colors[i],
                  arrow_length_ratio=0.1, alpha=alpha, linewidth=2)
        tip = origin + direction * 1.15
        ax.text(*tip, f'{label_prefix}{labels[i]}', color=colors[i],
                fontsize=10, alpha=alpha)

fig = plt.figure(figsize=(16, 5))

# (a) z軸まわり45度回転
ax1 = fig.add_subplot(131, projection='3d')
plot_frame(ax1, np.eye(3), label_prefix='', alpha=0.3)
plot_frame(ax1, Rz(np.radians(45)), label_prefix="'")
ax1.set_title(r'$R_z(45°)$', fontsize=13)
ax1.set_xlim([-1.2, 1.2]); ax1.set_ylim([-1.2, 1.2]); ax1.set_zlim([-1.2, 1.2])
ax1.set_xlabel('X'); ax1.set_ylabel('Y'); ax1.set_zlabel('Z')

# (b) x軸まわり45度回転
ax2 = fig.add_subplot(132, projection='3d')
plot_frame(ax2, np.eye(3), label_prefix='', alpha=0.3)
plot_frame(ax2, Rx(np.radians(45)), label_prefix="'")
ax2.set_title(r'$R_x(45°)$', fontsize=13)
ax2.set_xlim([-1.2, 1.2]); ax2.set_ylim([-1.2, 1.2]); ax2.set_zlim([-1.2, 1.2])
ax2.set_xlabel('X'); ax2.set_ylabel('Y'); ax2.set_zlabel('Z')

# (c) y軸まわり45度回転
ax3 = fig.add_subplot(133, projection='3d')
plot_frame(ax3, np.eye(3), label_prefix='', alpha=0.3)
plot_frame(ax3, Ry(np.radians(45)), label_prefix="'")
ax3.set_title(r'$R_y(45°)$', fontsize=13)
ax3.set_xlim([-1.2, 1.2]); ax3.set_ylim([-1.2, 1.2]); ax3.set_zlim([-1.2, 1.2])
ax3.set_xlabel('X'); ax3.set_ylabel('Y'); ax3.set_zlabel('Z')

plt.tight_layout()
plt.savefig('basic_rotations.png', dpi=150, bbox_inches='tight')
plt.show()

上のグラフでは、3つの基本回転のそれぞれの効果を視覚化しています。薄い色の矢印が回転前(基準座標系)、濃い色の矢印が回転後の座標系です。

  1. $R_z(45°)$(左図): $x$ 軸と $y$ 軸が $z$ 軸のまわりに45°回転し、$z$ 軸は不変です。上から見下ろすと反時計回りの回転に見えます。
  2. $R_x(45°)$(中図): $y$ 軸と $z$ 軸が $x$ 軸のまわりに45°回転し、$x$ 軸は不変です。
  3. $R_y(45°)$(右図): $z$ 軸と $x$ 軸が $y$ 軸のまわりに45°回転し、$y$ 軸は不変です。

いずれの場合も、回転軸に対応する基底ベクトルが動かないことが確認できます。

合成回転の順序依存性の可視化

次に、回転の順序を入れ替えると結果が変わることを視覚的に確認します。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def Rx(phi):
    c, s = np.cos(phi), np.sin(phi)
    return np.array([[1,0,0],[0,c,-s],[0,s,c]])

def Ry(theta):
    c, s = np.cos(theta), np.sin(theta)
    return np.array([[c,0,s],[0,1,0],[-s,0,c]])

def Rz(psi):
    c, s = np.cos(psi), np.sin(psi)
    return np.array([[c,-s,0],[s,c,0],[0,0,1]])

def plot_frame(ax, R, origin=np.zeros(3), label_prefix='',
               colors=('r','g','b'), alpha=1.0, length=1.0):
    labels = ['x', 'y', 'z']
    for i in range(3):
        direction = R[:, i] * length
        ax.quiver(*origin, *direction, color=colors[i],
                  arrow_length_ratio=0.1, alpha=alpha, linewidth=2)
        tip = origin + direction * 1.15
        ax.text(*tip, f'{label_prefix}{labels[i]}', color=colors[i],
                fontsize=10, alpha=alpha)

angle = np.radians(90)

fig = plt.figure(figsize=(12, 5))

# (a) Rx(90) @ Rz(90): まずz軸回転、次にx軸回転
ax1 = fig.add_subplot(121, projection='3d')
C1 = Rx(angle) @ Rz(angle)
plot_frame(ax1, np.eye(3), label_prefix='', alpha=0.3)
plot_frame(ax1, C1, label_prefix="'")
ax1.set_title(r'$R_x(90°) \cdot R_z(90°)$' + '\n(z回転→x回転)', fontsize=12)
ax1.set_xlim([-1.5,1.5]); ax1.set_ylim([-1.5,1.5]); ax1.set_zlim([-1.5,1.5])
ax1.set_xlabel('X'); ax1.set_ylabel('Y'); ax1.set_zlabel('Z')

# (b) Rz(90) @ Rx(90): まずx軸回転、次にz軸回転
ax2 = fig.add_subplot(122, projection='3d')
C2 = Rz(angle) @ Rx(angle)
plot_frame(ax2, np.eye(3), label_prefix='', alpha=0.3)
plot_frame(ax2, C2, label_prefix="'")
ax2.set_title(r'$R_z(90°) \cdot R_x(90°)$' + '\n(x回転→z回転)', fontsize=12)
ax2.set_xlim([-1.5,1.5]); ax2.set_ylim([-1.5,1.5]); ax2.set_zlim([-1.5,1.5])
ax2.set_xlabel('X'); ax2.set_ylabel('Y'); ax2.set_zlabel('Z')

plt.tight_layout()
plt.savefig('rotation_order.png', dpi=150, bbox_inches='tight')
plt.show()

上の2つのグラフを比較すると、同じ回転角($90°$ずつ)を使っていても、適用する順序を入れ替えると最終的な座標系の向きが全く異なることがわかります。左図では $x’$ が $z$ 方向を向いていますが、右図では $x’$ が $y$ 方向を向いています。これは回転が非可換($R_x R_z \neq R_z R_x$)であることの視覚的な証拠です。

物体の回転の可視化

最後に、回転行列による物体(直方体)の座標変換を可視化します。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection

def Rx(phi):
    c, s = np.cos(phi), np.sin(phi)
    return np.array([[1,0,0],[0,c,-s],[0,s,c]])

def Ry(theta):
    c, s = np.cos(theta), np.sin(theta)
    return np.array([[c,0,s],[0,1,0],[-s,0,c]])

def Rz(psi):
    c, s = np.cos(psi), np.sin(psi)
    return np.array([[c,-s,0],[s,c,0],[0,0,1]])

# 直方体の8頂点を定義(衛星を簡易的に表現)
L, W, H = 1.0, 0.6, 0.4  # 長さ、幅、高さ
vertices = np.array([
    [-L/2, -W/2, -H/2], [ L/2, -W/2, -H/2],
    [ L/2,  W/2, -H/2], [-L/2,  W/2, -H/2],
    [-L/2, -W/2,  H/2], [ L/2, -W/2,  H/2],
    [ L/2,  W/2,  H/2], [-L/2,  W/2,  H/2]
])

# 6つの面を定義
faces = [
    [0,1,2,3], [4,5,6,7],  # 底面・上面
    [0,1,5,4], [2,3,7,6],  # 前面・背面
    [0,3,7,4], [1,2,6,5]   # 左面・右面
]

face_colors = ['#2196F3', '#2196F3', '#FF9800', '#FF9800', '#4CAF50', '#4CAF50']

def plot_box(ax, verts, title='', alpha=0.6):
    """直方体をポリゴンで描画"""
    for face_idx, face in enumerate(faces):
        poly = [[verts[v] for v in face]]
        ax.add_collection3d(Poly3DCollection(
            poly, alpha=alpha, facecolor=face_colors[face_idx],
            edgecolor='black', linewidth=0.5))

fig = plt.figure(figsize=(16, 5))

# (a) 回転前
ax1 = fig.add_subplot(131, projection='3d')
plot_box(ax1, vertices, alpha=0.4)
ax1.set_title('Original', fontsize=13)
ax1.set_xlim([-1,1]); ax1.set_ylim([-1,1]); ax1.set_zlim([-1,1])
ax1.set_xlabel('X'); ax1.set_ylabel('Y'); ax1.set_zlabel('Z')

# (b) z軸まわり30度回転
C1 = Rz(np.radians(30))
v_rot1 = (C1 @ vertices.T).T
ax2 = fig.add_subplot(132, projection='3d')
plot_box(ax2, vertices, alpha=0.15)
plot_box(ax2, v_rot1, alpha=0.6)
ax2.set_title(r'After $R_z(30°)$', fontsize=13)
ax2.set_xlim([-1,1]); ax2.set_ylim([-1,1]); ax2.set_zlim([-1,1])
ax2.set_xlabel('X'); ax2.set_ylabel('Y'); ax2.set_zlabel('Z')

# (c) さらにx軸まわり45度回転
C2 = Rx(np.radians(45)) @ Rz(np.radians(30))
v_rot2 = (C2 @ vertices.T).T
ax3 = fig.add_subplot(133, projection='3d')
plot_box(ax3, vertices, alpha=0.15)
plot_box(ax3, v_rot2, alpha=0.6)
ax3.set_title(r'After $R_x(45°) \cdot R_z(30°)$', fontsize=13)
ax3.set_xlim([-1,1]); ax3.set_ylim([-1,1]); ax3.set_zlim([-1,1])
ax3.set_xlabel('X'); ax3.set_ylabel('Y'); ax3.set_zlabel('Z')

plt.tight_layout()
plt.savefig('box_rotation.png', dpi=150, bbox_inches='tight')
plt.show()

上の3つのパネルから、回転行列の効果を直感的に理解できます。

  1. 左図(Original): 直方体は座標軸に平行に配置されています。
  2. 中図($R_z(30°)$後): $z$ 軸まわりに30°回転しました。上から見ると反時計回りに回っていますが、高さ方向には変化がありません。薄い直方体が回転前の位置を示しています。
  3. 右図($R_x(45°) \cdot R_z(30°)$後): さらに $x$ 軸まわりに45°回転を加えました。直方体が2軸分の回転で傾いた姿勢を取っています。

このように、回転行列を順に適用することで、任意の姿勢を表現できることが視覚的にも確認できます。

回転行列の性質の数値検証

最後に、合成回転行列の直交性・行列式・ベクトルのノルム保存を数値的に検証しておきます。

import numpy as np

def Rx(phi):
    c, s = np.cos(phi), np.sin(phi)
    return np.array([[1,0,0],[0,c,-s],[0,s,c]])

def Ry(theta):
    c, s = np.cos(theta), np.sin(theta)
    return np.array([[c,0,s],[0,1,0],[-s,0,c]])

def Rz(psi):
    c, s = np.cos(psi), np.sin(psi)
    return np.array([[c,-s,0],[s,c,0],[0,0,1]])

# 任意角度の合成回転
phi, theta, psi = np.radians(30), np.radians(45), np.radians(60)
C = Rx(phi) @ Ry(theta) @ Rz(psi)

print("合成回転行列 C:")
print(np.round(C, 6))
print(f"\nC^T @ C (単位行列に一致すべき):")
print(np.round(C.T @ C, 10))
print(f"\ndet(C) = {np.linalg.det(C):.15f}")

# ベクトルのノルム保存
r = np.array([1.0, 2.0, 3.0])
r_rot = C @ r
print(f"\n元のベクトルのノルム:   {np.linalg.norm(r):.15f}")
print(f"回転後のベクトルのノルム: {np.linalg.norm(r_rot):.15f}")

# 逆行列 = 転置の検証
C_inv = np.linalg.inv(C)
print(f"\nC^(-1) と C^T の差のノルム: {np.linalg.norm(C_inv - C.T):.2e}")

この出力から、合成回転行列についても以下の性質が数値的に確認できます。

  1. 直交性: $C^T C$ が単位行列に一致し、数値誤差は $10^{-16}$ オーダー
  2. 行列式: $\det(C) = 1.000…$ で理論値と一致
  3. ノルム保存: 回転前後でベクトルの長さが完全に保存
  4. 逆行列=転置: $C^{-1}$ と $C^T$ の差が事実上ゼロ

これらの結果は、回転行列の理論的性質が数値計算でも精度よく保たれることを示しています。

まとめ

本記事では、剛体の姿勢表現の基礎として、座標変換と回転行列について解説しました。

  • 座標系と姿勢: 姿勢とは2つの座標系の間の回転関係であり、方向余弦行列(DCM)で記述される
  • 回転行列の性質: 直交性($C^T C = I$)、行列式1、ベクトルのノルム保存を持ち、特殊直交群 SO(3) に属する
  • 基本回転行列: $R_x$, $R_y$, $R_z$ の3種類を三角関数で具体的に導出した
  • 合成回転: 複数の回転を行列の積で合成できるが、順序に注意が必要(非可換性)
  • 逆行列 = 転置: 回転行列の逆変換は転置するだけで得られ、計算コストが低い

回転行列は9個の要素を持ちますが、独立なパラメータはわずか3個です。このパラメータを明示的に3つの角度で表す方法として、次の記事ではオイラー角を扱います。オイラー角は直感的で便利ですが、「ジンバルロック」と呼ばれる特異点の問題を抱えています。

次のステップとして、以下の記事も参考にしてください。