Vincenty法

Vincenty法（ヴィンセンティほう、テンプレート:Lang-en-short）とは回転楕円体上の測地法の反復計算アルゴリズムである。テンプレート:仮リンク (1975a)によって考案された。

• 地球形状は回転楕円体として近似可能であり、Vincenty法の逆解法は、その面上で二点間の最小距離（測地線の長さ）を計算する。ただしVincenty法の計算誤差は、与える入力点によっては 0.5 mm(0.020″)に達する場合もある。
• 近年は Karney (2013) による計算方法がVincenty法よりも計算精度・計算速度ともに高く、多く用いられる。

背景

測地法（測地線の長さに関する計算）において、順解法(direct method)はある点からの距離と方位角により、他の点を求めるものである。逆解法(inverse method)は2点間の回転楕円体面上距離と方位角を求めるものである。

Legendreは回転楕円体上の測地線は緯度を更成緯度とし、方位角を同じものとすることで補助球の大円上に正確に対応させられることを示した。 楕円体上の経度と距離は単純な積分によって、補助球上の経度とその大円上の弧長から導かれる。BesselとHelmertはそれらの積分について収束の早い数列を与えた。それにより測地計算が任意の精度で行えるようになった。

Vincenty法の工夫

Vincentyの目的は当時すでに存在していた回転楕円体上の測地法アルゴリズムのプログラムの長さを最小化することであった。彼の未発行のレポート(1975b)によれば数キロバイトのメモリしか搭載していないWang 720 電卓を使用した。

長距離において良い精度を得るために、伝統的な補助球を用いた手法（Legendre (1806), Bessel (1825), and Helmert (1880)）を使用している。また、Vincentyの手法はRainsford (1955)による解法にも依存している。

プログラムサイズを最小化するためにVincentyはこれらの数列を初項を小さな値とする拡張し$f^3$のオーダーで打ち切った。これにより、経度と距離の積分式は簡潔なものとなった。式は入れ子型（ホーナー法）にし、一時レジスタを1つだけ用いて多項式計算を行えるようにした。そして、単純な収束計算を陰関数を解くのに用いた。これは収束が遅くinverse法では収束しないこともありうるが、コードサイズは小さいものとなった。

表記法

以下のように定義する。

$a$

長軸半径 (赤道半径) (6378137.06 m (WGS84) )

$f$

扁平率 (1/298.257223563 (WGS84) )

$b=(1-f)a$

短軸半径 (極半径) ( 6356752.314245 m (WGS84) )

${\varphi }_1,{\varphi }_2$

各点の緯度

$U_1=\arctan [(1-f)\tan {\varphi }_1]$
$U_2=\arctan [(1-f)\tan {\varphi }_2]$

更成緯度 (補助球上の緯度)

$L=L_2-L_1$

2点間の経度差

${\lambda }_1,{\lambda }_2$

補助球上の経度

${\alpha }_1,{\alpha }_2$

各点における方位角

$\alpha$

赤道上での方位角

$s$

2点間の楕円体上の距離

$\sigma$

補助球上の弧の長さ

逆解法

逆解法は回転楕円体上の2点の座標$({\varphi }_1,L_1),({\varphi }_2,L_2)$が与えられた時に、方位角${\alpha }_1,{\alpha }_2$と距離$s$を導く。

$U_1,U_2$及び$L$を計算し、$\lambda$を$\lambda =L$で初期化し、以下の計算を$\lambda$が収束するまで反復する。

$\sin \sigma =\sqrt{(\cos U_2\sin \lambda {)}^2+(\cos U_1\sin U_2-\sin U_1\cos U_2\cos \lambda {)}^2}$

$\cos \sigma =\sin U_1\sin U_2+\cos U_1\cos U_2\cos \lambda$

$\sigma =\arctan \frac{\sin \sigma }{\cos \sigma }$

$\sin \alpha =\frac{\cos U_1\cos U_2\sin \lambda }{\sin \sigma }$

${\cos }^2\alpha =1-{\sin }^2\alpha$

$\cos (2{\sigma }_m)=\cos \sigma -\frac{2\sin U_1\sin U_2}{{\cos }^2\alpha }$

$C=\frac{f}{16}{\cos }^2\alpha [4+f(4-3{\cos }^2\alpha )]$

$\lambda =L+(1-C)f\sin \alpha \{\sigma +C\sin \sigma [\cos (2{\sigma }_m)+C\cos \sigma (-1+2{\cos }^2(2{\sigma }_m))]\}$

$\lambda$が所望の精度まで収束したら（10^-12の偏差ならば0.06 mmの精度）以下の計算を行う。

$u^2={\cos }^2\alpha \frac{a^2-b^2}{b^2}$

$A=1+\frac{u^2}{16384}\{4096+u^2[-768+u^2(320-175u^2)]\}$

$B=\frac{u^2}{1024}\{256+u^2[-128+u^2(74-47u^2)]\}$

$\mathrm{\Delta }\sigma =B\sin \sigma \{\cos (2{\sigma }_m)+\frac{1}{4}B[\cos \sigma (-1+2{\cos }^2(2{\sigma }_m))-\frac{1}{6}B\cos (2{\sigma }_m)(-3+4{\sin }^2\sigma )(-3+4{\cos }^2(2{\sigma }_m))\left]\right\}$

$s=bA(\sigma -\mathrm{\Delta }\sigma )$

${\alpha }_1=\arctan \left(\frac{\cos U_2\sin \lambda }{\cos U_1\sin U_2-\sin U_1\cos U_2\cos \lambda }\right)$

${\alpha }_2=\arctan \left(\frac{\cos U_1\sin \lambda }{-\sin U_1\cos U_2+\cos U_1\sin U_2\cos \lambda }\right)$

両極付近の2点間の計算は最初の$\lambda$を設定した際に$\pi$以上となると、収束しない。

順解法

順解法は始点$({\varphi }_1,L_1)$と始点における方位角${\alpha }_1$、距離$s$が与えられた時の終点の座標$({\varphi }_2,L_2)$と方位角${\alpha }_2$を導く。

まず、以下の計算を行う。

$U_1=\arctan ((1-f)\tan {\varphi }_1)$

${\sigma }_1=\arctan \left(\frac{\tan U_1}{\cos {\alpha }_1}\right)$

$\sin \alpha =\cos U_1\sin {\alpha }_1$

${\cos }^2\alpha =1-{\sin }^2\alpha$

$u^2={\cos }^2\alpha \frac{a^2-b^2}{b^2}$

$A=1+\frac{u^2}{16384}\{4096+u^2[-768+u^2(320-175u^2)]\}$

$B=\frac{u^2}{1024}\{256+u^2[-128+u^2(74-47u^2)]\}$

$\sigma$を$\sigma =\frac{s}{bA}$で初期化し、以下の収束するまで反復計算を行う。

$2{\sigma }_m=2{\sigma }_1+\sigma$

$\mathrm{\Delta }\sigma =B\sin \sigma \{\cos (2{\sigma }_m)+\frac{1}{4}B[\cos \sigma (-1+2{\cos }^2(2{\sigma }_m))-\frac{1}{6}B\cos (2{\sigma }_m)(-3+4{\sin }^2\sigma )(-3+4{\cos }^2(2{\sigma }_m))\left]\right\}$

$\sigma =\frac{s}{bA}+\mathrm{\Delta }\sigma$

$\sigma$が所望の精度まで収束したら以下の計算を行う。

${\varphi }_2=\arctan \left(\frac{\sin U_1\cos \sigma +\cos U_1\sin \sigma \cos {\alpha }_1}{(1-f)\sqrt{{\sin }^2\alpha +(\sin U_1\sin \sigma -\cos U_1\cos \sigma \cos {\alpha }_1{)}^2}}\right)$

$\lambda =\arctan \left(\frac{\sin \sigma \sin {\alpha }_1}{\cos U_1\cos \sigma -\sin U_1\sin \sigma \cos {\alpha }_1}\right)$

$C=\frac{f}{16}{\cos }^2\alpha [4+f(4-3{\cos }^2\alpha )]$

$L=\lambda -(1-C)f\sin \alpha \{\sigma +C\sin \sigma [\cos (2{\sigma }_m)+C\cos \sigma (-1+2{\cos }^2(2{\sigma }_m))]\}$

$L_2=L+L_1$

${\alpha }_2=\arctan \left(\frac{\sin \alpha }{-\sin U_1\sin \sigma +\cos U_1\cos \sigma \cos {\alpha }_1}\right)$

始点が極点である時、最初の等式は不定となる。 始点の方位角が真西か真東の場合2番目の等式は不定となる。 しかし、2変数関数のatan2等を用いることでこれらの値は正しく計算できる。

Vincentyによる修正

Vincenty自身の論文(1979)によると、$A,B$をHelmertの展開係数$k_1$を用いて以下のより簡潔な形式で表すことを提案した。

$A=\frac{1+\frac{1}{4}(k_1{)}^2}{1-k_1}$

$B=k_1(1-\frac{3}{8}(k_1{)}^2)$

ここで、$k_1=\frac{\sqrt{(1+u^2)}-1}{\sqrt{(1+u^2)}+1}$である。

対蹠付近の点

先に触れたように、逆解法の反復計算は対蹠点付近において、収束しなかったり、収束が遅かったりする。例えば、WGS84測地系において$({\varphi }_1,L_1)=(0^{\circ },0^{\circ }),({\varphi }_2,L_2)=(0.5^{\circ },179.5^{\circ })$であるとする。この計算は1 mmの精度に達するまでに130回の計算が必要になる。逆解法がどのように実装されるかによって、計算結果は正しく19936288.579 mを返したり、間違った結果を返したり、エラーとなったりする。例えばNGS online utilityによる計算結果は5 kmも長い。Vincentyが提案した方法はこのような場合に収束を速める(Rapp, 1973)。

逆解法が収束しない例として$({\varphi }_1,L_1)=(0^{\circ },0^{\circ }),({\varphi }_2,L_2)=(0.5^{\circ },179.7^{\circ })$がある。Vincentyの未出版のレポート(1975b)によると、異なった手法がそのような状況に適しており、60回の反復計算の後に19944127.421 mという正しい計算結果に達した。しかし、この手法は他の多くの場合において何千回もの反復が必要となる。

関連項目

• 経緯度
• 大円距離
• 子午線弧
• 測地法
• 測地線#回転楕円体面上の測地線

参考文献

• テンプレート:Cite journal
• テンプレート:Cite book
• テンプレート:Cite journal
• テンプレート:Cite journal
• テンプレート:Cite journal
• テンプレート:Cite techreport
• テンプレート:Cite journal
• テンプレート:Cite journal
• テンプレート:Cite techreport
• テンプレート:Cite bookテンプレート:リンク切れ

外部リンク

• Online calculators from Geoscience Australia:
  • Vincenty Direct (destination point)
  • Vincenty Inverse (distance between points)
• Calculators from the U.S. National Geodetic Survey:
  • Online and downloadable PC-executable calculation utilities, including forward (direct) and inverse problems, in both two and three dimensions (accessed 2011-08-01).
• Online calculators with JavaScript source code by Chris Veness (Creative Commons Attribution license):
  • Vincenty Direct (destination point)
  • Vincenty Inverse (distance between points)
• GeographicLib provides a utility GeodSolve (with MIT/X11 licensed source code) for solving direct and inverse geodesic problems. Compared to Vincenty, this is about 1000 times more accurate (error = 15 nm) and the inverse solution is complete. Here is an online version of GeodSolve.