「逆変換の式」は存在しない — Newton-Raphson反復法を数式で理解する

⚠ 本記事は、「GeoCoreJPの高精度座標変換を支える設計思想」で紹介した「逆変換は順変換の反復で解く」という設計判断の数式特化版です。前回はSwiftのコードで「どう実装しているか」を示しましたが、今回は「なぜその方法しか選べないのか」「反復計算の中で何が起きているのか」を数式で追います。また、連載㉑（往復変換の誤差は数cmではなくmm未満だった）で報告した実測値の、理論的な裏付けでもあります。

なぜ「逆向きの式」を書けないのか

平面直角座標系への変換や、TKY2JGD・jgd2024bのようなグリッド補正は、緯度経度 $(\phi, \lambda)$ を入力として、変換後の座標 $(\phi’, \lambda’)$ を返す関数 $F$ として書けます。

$$
(\phi’, \lambda’) = F(\phi, \lambda)
$$

一見、これを逆に解いて $(\phi, \lambda) = F^{-1}(\phi’, \lambda’)$ という「逆向きの式」を作れそうに思えます。しかし $F$ は次の2つの性質を持つため、閉形式（有限個の四則演算・三角関数などで書き下せる式）の逆関数 $F^{-1}$ は存在しません。

グリッド補正が区分的・非線形：TKY2JGD・jgd2024bなどの補正量は、格子点ごとの補正値をバイリニア補間して求めます。補正量そのものが「変換前の位置」に応じて連続的に変化する非線形関数であり、代数的に逆を解ける形をしていません。
投影が三角関数・楕円積分を含む：ガウス・クリューゲル投影は緯度経度から平面座標への変換式ですが、その逆（平面座標から緯度経度）は、子午線弧長の逆関数を含み、閉形式では書けません（国土地理院が公開する平面直角座標→緯度経度の算式も、伝統的には反復計算を用いてきました。2026年7月に反復不要の新算式が公開されましたが、これは反復を避けるための別の近似式であり、$F$ の厳密な逆関数を解析的に導いたものではありません）。

つまり「順変換の式」は書けても、「逆変換の式」は原理的に書けない場面が実務では珍しくありません。

素朴な方法とその誤差

閉形式が無いとき、よくある簡易的な回避策は「補正量をそのまま逆向きに引く」方法です。

$$
(\phi, \lambda) \approx (\phi’, \lambda’) – \Delta(\phi’, \lambda’)
$$

ここで $\Delta$ は変換後の座標の位置で評価した補正量です。しかし補正量 $\Delta$ は本来「変換前の位置」で評価されるべきものなので、この方法は評価点がずれた近似にすぎません。補正量の空間的な変化が小さい場所ではほぼ問題になりませんが、地殻変動量が大きい地震の震源近傍などでは、この評価点のずれがそのまま誤差として現れます。

Newton-Raphson反復法 — 順変換だけで逆変換を解く

GeoCoreJPが採用しているのは、逆向きの式を作る代わりに、順変換 $F$ を繰り返し呼び出して目標値に収束させる方法です。

1変数版のNewton-Raphson法（考え方の確認）

一般に $f(x) = 0$ の解を求めるNewton-Raphson法は、次の更新式で反復します。

$$
x_{n+1} = x_n – \frac{f(x_n)}{f'(x_n)}
$$

これを座標変換に当てはめると、「目標座標 $(\phi’, \lambda’)$ に対して、$F(\phi, \lambda) – (\phi’, \lambda’) = 0$ となる $(\phi, \lambda)$ を求める」問題になります。

2変数版 — ヤコビアンによる更新

緯度経度は2変数なので、実際には2×2のヤコビアン行列 $J$ を使います。

$$
J = \begin{pmatrix} \dfrac{\partial \phi’}{\partial \phi} & \dfrac{\partial \phi’}{\partial \lambda} [4pt] \dfrac{\partial \lambda’}{\partial \phi} & \dfrac{\partial \lambda’}{\partial \lambda} \end{pmatrix}
$$

更新式は次のとおりです（$n$ 回目の推定値を $(\phi_n, \lambda_n)$、残差を $e_n = (\phi’, \lambda’) – F(\phi_n, \lambda_n)$ とする）。

$$
\begin{pmatrix} \phi_{n+1} \ \lambda_{n+1} \end{pmatrix} = \begin{pmatrix} \phi_n \ \lambda_n \end{pmatrix} + J^{-1} e_n
$$

GeoCoreJPの実装（前回記事で紹介したコード）では、ヤコビアン $J$ を解析的に微分するのではなく、数値微分（$\phi, \lambda$ をごくわずかに動かして $F$ を2回余分に評価する差分近似）で求めています。解析的な偏微分をすべての補正・投影の組み合わせごとに導出・保守するコストを避け、「順変換関数さえあれば逆変換が動く」という汎用性を優先した設計判断です。

ステップ	内容
① 初期値	目標座標 $(\phi’, \lambda’)$ をそのまま初期推定値とする
② 順変換評価	現在の推定値を $F$ に通し、残差 $e_n$ を計算
③ 収束判定	残差をメートル換算し、閾値未満なら終了
④ ヤコビアン更新	数値微分でヤコビアンを求め、$J^{-1}e_n$ だけ推定値を更新
⑤ ①〜④を繰り返す	最大反復回数に達したら打ち切り

収束判定と発散対策 — なぜ「0.1mm・20回」なのか

GeoCoreJPのJGD補正逆変換は収束閾値 0.1mm・最大反復20回、ガウス・クリューゲル逆投影は収束閾値 1mm・最大反復12回で実装されています。この数値には理由があります。

0.1mm（JGD補正）：測量成果として求められる精度（mm〜cm級）に対して、計算誤差を一桁以上小さく抑えるための余裕を持たせた値です。厳しすぎる閾値は反復回数を無駄に増やし、緩すぎる閾値は精度を損ないます。
最大反復回数：Newton-Raphson法は初期値が真の解に十分近ければ数回で収束しますが、理論上は発散する初期値も存在します。座標変換の場合、初期値（目標座標そのもの）と真の解の差は、補正量やゾーン境界付近でも高々数百m程度に収まるため、実務上はほぼ確実に数回〜十数回で収束します。上限を設けるのは、万一収束しない異常な入力（座標系の取り違えなど）が来たときに、無限ループを防ぐための安全装置です。