地震の解析では、震源の断層がどのように動いたかという「発震機構」（focal mechanism）を求める。この解を可視化する際は「ビーチボール」の形で図を描く。大きな地震の後の気象庁会見でもチラッと出てきたりする¹。

実用上は The Generic Mapping Tools (GMT) などのツールを使えば簡単に描けるものの、遊びの一環で自分でコードを組んで描いてみたくなった。とはいえ、初動発震機構解であれば白黒が4象限にきれいに分かれるので何とか描けそうだが、セントロイドモーメントテンソル（CMT）解では白黒の境目が微妙な曲線になるため描き方が見当つかない。

そこで、GMTのソースコードを読んで勉強し、JavaScriptで実装してみた。

デモ

CodePenの埋め込みを表示する

See the Pen Beachball Viewer by hmmnrst (@hmmnrst) on CodePen.

冒頭の図を真似するには、例えば以下のように入力する。

Format: -Sa / Parameters: 130 83 -6
Format: -Sm / Parameters: -7 150 -142 31 61 -65

ソースコード： https://github.com/hmmnrst/beachball-viewer
（TODO: READMEを作成する）

実装

以下の流れで処理する。

発震機構解のパラメーターを読み込む
0. 例えば、テキストボックスのchangeイベントを契機に開始する
0. データ形式に応じてパースする
0. 最終的に主軸/主値3組と節面2つを求める
押し引きの分布を描く
0. 中立軸の押し引きに応じて、下地と上塗りの色を決める
0. 下地の円を circle で描く
0. 上塗り部のパスを求める
0. パスを下半球と上半球に分割し、それぞれ polygon で描く
節面の曲線を polyline で描く
縁の円周を circle で描く

GMTからは主に以下の点を変えてある。

軸や節面などは極座標（方位角 azimuth ・傾斜角 plunge ）でなく直交直線座標（深さ成分・北成分・東成分）²で管理
ダブルカップル（DC）成分のみであっても、モーメントテンソルと同じ方法で描画
その他
- ビーチボールの大きさは固定（地震モーメントあるいはマグニチュードを無視）
- パスの計算間隔を設定可能 → 15°にして、描画が折れ線であることを強調
- ヤコビ法は3x3専用にして、冗長になったfor文は展開

押し引きの分布の描画

今回勉強するにあたって最大の疑問だった部分。

ベクター画像として押し引きの分布を描くには、境界を求めてパスを描くことになる。この境界線は言い換えれば「ビーチボール表面で力の垂直成分がゼロとなっている点の集合」である。3x3のモーメントテンソルを M 、方向ベクトルを r とすれば、境界線は $\boldsymbol{r} \boldsymbol{M} \boldsymbol{r} = \sum M_{ij} r_i r_j = 0$ を満たす r として表せる。ただ、これをそのまま解くのは難しそう。

ここで、 M が対角化済み（対角成分 v_i 、他はゼロ）なら方程式は簡単になる： $\sum v_i r_i^2 = 0$ 。これは主軸に沿った座標系で考えていることになり、解を元の座標系に変換すれば本来の位置がわかる。

入力がモーメントテンソルなら、実対称行列と捉えられるので固有値は全て実数、固有ベクトルの組は直交行列に必ずなる。GMTではヤコビ法を用いて求めている。
入力が断層パラメーターなら、直交する節面が先に求まるので主軸を改めて求める必要がある。地震モーメントを M₀ とすると、主値は +M₀, 0, -M₀ を割り当てればいい³。

方程式は r_i が全て2乗になっているので、主値 v₁ ≥ v₂ ≥ v₃ の符号が重要になる。

符号が全て同じ → 境界なし、つまり全方向に押し（膨張）か引き（収縮）
v₂ がマイナス → 境界は張力軸を1周
v₂ がプラス → 境界は圧力軸を1周
v₂ がゼロ → 境界は中立軸を通る大円2つで⁴、見方を変えれば張力軸も圧力軸も1周

境界の無い最初のパターンを除けば、中立軸の主値 v₂ の符号によって境界が1周する軸が決まる。この軸を中心にして、極座標で方位毎に境界の位置を求めれば、上塗り部のパスが球面上で求まる（1周しない軸を中心にすると、境界の無い方位があり計算が破綻する）。例えば v₂ がマイナスの場合は次式のようになり、全方位 θ に対して傾斜角 α が求まることが確かめられる。

\tan^2{\alpha} = \frac{- v_2 \cos^2{\theta} - v_3 \sin^2{\theta}}{v_1} \ge 0

解 α(θ) は正負の2つある（＝境界が2本ある）が、片方だけ求めておけばいい。

※ GMTでは α を求めるのはもっと複雑な計算式だったが、結局読み解けなかった…。

境界を元の座標系に変換すると、図で表示しない上半球へ飛び出してしまうことがあるので、クリッピングが必要になる。ビーチボールの対称性からすると、飛び出した部分は捨てるのではなく球面の正反対の位置に描けばいい。

最後に紙面に描く際は、球面上で求めたパスを「下半球等積投影」で平面上の座標に変換し、SVGの polygon で描いて塗りつぶす。折れ線であってもパスを細かくすればほぼ曲線に見える（今回は折れ線を強調するために粗くしている）。

節面の管理

DC成分のみの場合は押し引きが4象限に分かれ、その境界は2つの直交する面になる。実際の断層面に対応するのは片方だけだが、発震機構解からは判別できない。

これらの節面は、断層運動を表すパラメーター（走向 strike ・傾斜角 dip ・すべり角 slip/rake ）だけで決められる。ただ、2面が対称的なことを考えると、すぐに面を切り替えられる方法で管理した方が楽だと感じた。

上盤のすべり方向をベクトルで表す
- これは断層パラメーターの模式図そのまま
- ベクトルのz成分によって正断層・逆断層が決まる
断層面は法線ベクトルで表す
- 走向と傾斜角のみから求まる
- 上盤のある方向を選ぶ（z成分が0以下）

こうすると2つのベクトルはビーチボールの黒い部分を挟む形になる。もうひとつの面を表すには、単に2つのベクトルを入れ替えればいい（上盤を表すよう、必要に応じて両ベクトルを符号反転する）。張力軸・圧力軸と成す角は45°なので、主軸との相互変換も簡単。

節面の曲線を描く際は、押し引きの境界と同様に、球面上でのパスを計算してから投影する。GMTは極座標で計算しているため、方位ごとに傾斜角を求めている。今回の実装では直交直線座標で計算しているため、断層面上ですべり角を 0° ～ -180° と変化させて方向ベクトルを求めている。

扱うデータ形式

GMTでは5種類の入力形式があるので、それら全てを扱えるようにしたい。最終的にモーメントテンソルの主軸/主値を求めることで、描画処理は前述の方法に統一できる。

非DC成分⁵を含む（自由度6、うち1つは大きさ）
- モーメントテンソル：名前通りテンソルの成分を指定する。対称なので6成分 (mrr, mtt, mff, mrt, mrf, mtf) のみ⁶。固有値と固有ベクトルを求めればいい⁷。
- 主軸：テンソルの3つの軸の向き（方位・傾斜角）と主値の大きさを直接指定する。9成分あり、軸が互いに直交している保証は無い。
DC成分のみ（自由度3、大きさは別）
- Aki and Richards convention：普通に断層パラメーターで指定する。ちょうど3成分。
- Global CMT convention：断層パラメーターを両方の節面について指定する。倍の6成分で、整合的な保証は無い。
- 部分的なデータ：節面1の走向と傾斜角、節面2の走向、断層の正/逆を指定する。3成分＋フラグなので自由度ぴったりだが、節面2の傾斜角を決定できないパターンがある。

成分の数が自由度を超えている形式については、GMTと同様の対応はできてなく、値が整合的と仮定して計算・描画している。

参考

気象庁｜地震・津波｜発震機構解（精査後）
- 地震月報（カタログ編）
- 発震機構解とは何か
- 震源球と断層面の関係
The Generic Mapping Tools
- ソースコード
  - src/seis/
  - src/gmt_vector.c 内の gmt_jacobi()
- Documentation » GMT Modules » meca
SVG.js

報道発表資料のPDFを見るのが確実。第2報などでなければ「CMT解（速報）」のページがあるはず。 ↩
右手系になっている。方位角・傾斜角との変換も楽。GMTでは処理の一部でのみ使っている。 ↩
モーメントテンソルとは M₀² = 1/2 Σ M_ij M_ji の関係がある（座標に依らない）。モーメントマグニチュード M_w はこれの対数を元に計算される。 ↩
双曲線の漸近線が2本の直線であることと似ている。 ↩
等方成分とCLVD成分。 ↩
r は鉛直上向き、 t (=θ) は南向き、 f (=φ) は東向き。 ↩
なお、固有値が縮退している場合は固有ベクトルに任意性が出る。押し引きの分布を描くには適当に選べばよく、ヤコビ法はうまく機能してくれる。 ↩