Vincenty法(ヴィンセンティほう、: Vincenty's formulae)とは回転楕円体上の測地法反復計算アルゴリズムである。サディアス・ヴィンセンティ英語版 (1975a)によって考案された。

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

背景

編集

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

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

Vincenty法の工夫

編集

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

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

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

表記法

編集

以下のように定義する。

 
長軸半径 (赤道半径) (6378137.06 m (WGS84) )
 
扁平率 (1/298.257223563 (WGS84) )
 
短軸半径 (極半径) ( 6356752.314245 m (WGS84) )
 
各点の緯度
 
 
更成緯度 (補助球上の緯度)
 
2点間の経度
 
補助球上の経度
 
各点における方位角
 
赤道上での方位角
 
2点間の楕円体上の距離
 
補助球上の弧の長さ

逆解法

編集

逆解法は回転楕円体上の2点の座標 が与えられた時に、方位角 と距離 を導く。

 及び を計算し、  で初期化し、以下の計算を が収束するまで反復する。

 
 
 [3][4]
 [5]
 
 [6]
 
 

 が所望の精度まで収束したら(10-12の偏差ならば0.06 mmの精度)以下の計算を行う。

 
 
 
 
 
 [4]
 [4]

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

順解法

編集

順解法は始点 と始点における方位角 、距離 が与えられた時の終点の座標 と方位角 を導く。

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

 
 [4]
 
 
 
 
 

  で初期化し、以下の収束するまで反復計算を行う。

 
 
 

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

 [4]
 [4]
 
 
 
 [4]

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

Vincentyによる修正

編集

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

 
 

ここで、 である。

対蹠付近の点

編集

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

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

関連項目

編集

注釈

編集
  1. ^ 地球を球体として近似すると大円距離の計算となるが、それよりは精度が良い。
  2. ^ Karney (2013) の方法ではニュートン法を用い、あらゆる入力点に対して収束が速く、倍精度浮動小数点数限界精度に達する。
  3. ^  は赤道及び極付近での精度を保つために から直接計算は行わない。
  4. ^ a b c d e f g arctanの値は2変数関数であるatan2等によって実装されるべきである。
  5. ^  の時、 が不定となる。これは2点が同じ位置もしくは対蹠点であることを意味する。
  6. ^ 2点がともに赤道上にあるとき、 となるので、 は使用されない。The limiting value is  .

参考文献

編集

外部リンク

編集
  • Online calculators from Geoscience Australia:
  • Calculators from the U.S. National Geodetic Survey:
  • Online calculators with JavaScript source code by Chris Veness (Creative Commons Attribution license):
  • 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.