Vincenty法

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

表記法

以下のように定義する。

a {\displaystyle a}
長軸半径 (赤道半径) (6378137.06 m (WGS84) )
f {\displaystyle f}
扁平率 (1/298.257223563 (WGS84) )
b = ( 1 − f ) a {\displaystyle b=(1-f)a}
短軸半径 (極半径) ( 6356752.314245 m (WGS84) )
ϕ 1 , ϕ 2 {\displaystyle \phi _{1},\phi _{2}}
各点の緯度
U 1 = arctan ⁡ [ ( 1 − f ) tan ⁡ ϕ 1 ] {\displaystyle U_{1}=\arctan[(1-f)\tan \phi _{1}]}
U 2 = arctan ⁡ [ ( 1 − f ) tan ⁡ ϕ 2 ] {\displaystyle U_{2}=\arctan[(1-f)\tan \phi _{2}]}
更成緯度 (補助球上の緯度)
L = L 2 − L 1 {\displaystyle L=L_{2}-L_{1}}
2点間の経度差
λ 1 , λ 2 {\displaystyle \lambda _{1},\lambda _{2}}
補助球上の経度
α 1 , α 2 {\displaystyle \alpha _{1},\alpha _{2}}
各点における方位角
α {\displaystyle \alpha }
赤道上での方位角
s {\displaystyle s}
2点間の楕円体上の距離
σ {\displaystyle \sigma }
補助球上の弧の長さ

逆解法

逆解法は回転楕円体上の2点の座標 ( ϕ 1 , L 1 ) , ( ϕ 2 , L 2 ) {\displaystyle (\phi _{1},L_{1}),(\phi _{2},L_{2})} が与えられた時に、方位角 α 1 , α 2 {\displaystyle \alpha _{1},\alpha _{2}} と距離 s {\displaystyle s} を導く。

U 1 , U 2 {\displaystyle U_{1},U_{2}} 及び L {\displaystyle L} を計算し、 λ {\displaystyle \lambda } を λ = L {\displaystyle \lambda =L} で初期化し、以下の計算を λ {\displaystyle \lambda } が収束するまで反復する。

sin ⁡ σ = ( cos ⁡ U 2 sin ⁡ λ ) 2 + ( cos ⁡ U 1 sin ⁡ U 2 − sin ⁡ U 1 cos ⁡ U 2 cos ⁡ λ ) 2 {\displaystyle \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 ⁡ σ = sin ⁡ U 1 sin ⁡ U 2 + cos ⁡ U 1 cos ⁡ U 2 cos ⁡ λ {\displaystyle \cos \sigma =\sin U_{1}\sin U_{2}+\cos U_{1}\cos U_{2}\cos \lambda \,}
σ = arctan ⁡ sin ⁡ σ cos ⁡ σ {\displaystyle \sigma =\arctan {\frac {\sin \sigma }{\cos \sigma }}\,} [3][4]
sin ⁡ α = cos ⁡ U 1 cos ⁡ U 2 sin ⁡ λ sin ⁡ σ {\displaystyle \sin \alpha ={\frac {\cos U_{1}\cos U_{2}\sin \lambda }{\sin \sigma }}\,} [5]
cos 2 ⁡ α = 1 − sin 2 ⁡ α {\displaystyle \cos ^{2}\alpha =1-\sin ^{2}\alpha \,}
cos ⁡ ( 2 σ m ) = cos ⁡ σ − 2 sin ⁡ U 1 sin ⁡ U 2 cos 2 ⁡ α {\displaystyle \cos(2\sigma _{m})=\cos \sigma -{\frac {2\sin U_{1}\sin U_{2}}{\cos ^{2}\alpha }}\,} [6]
C = f 16 cos 2 ⁡ α [ 4 + f ( 4 − 3 cos 2 ⁡ α ) ] {\displaystyle C={\frac {f}{16}}\cos ^{2}\alpha {\big [}4+f(4-3\cos ^{2}\alpha ){\big ]}\,}
λ = L + ( 1 − C ) f sin ⁡ α { σ + C sin ⁡ σ [ cos ⁡ ( 2 σ m ) + C cos ⁡ σ ( − 1 + 2 cos 2 ⁡ ( 2 σ m ) ) ] } {\displaystyle \lambda =L+(1-C)f\sin \alpha \left\{\sigma +C\sin \sigma \left[\cos(2\sigma _{m})+C\cos \sigma (-1+2\cos ^{2}(2\sigma _{m}))\right]\right\}\,}

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

u 2 = cos 2 ⁡ α a 2 − b 2 b 2 {\displaystyle u^{2}=\cos ^{2}\alpha {\frac {a^{2}-b^{2}}{b^{2}}}\,}
A = 1 + u 2 16384 { 4096 + u 2 [ − 768 + u 2 ( 320 − 175 u 2 ) ] } {\displaystyle A=1+{\frac {u^{2}}{16384}}\left\{4096+u^{2}\left[-768+u^{2}(320-175u^{2})\right]\right\}}
B = u 2 1024 { 256 + u 2 [ − 128 + u 2 ( 74 − 47 u 2 ) ] } {\displaystyle B={\frac {u^{2}}{1024}}\left\{256+u^{2}\left[-128+u^{2}(74-47u^{2})\right]\right\}}
Δ σ = B sin ⁡ σ { cos ⁡ ( 2 σ m ) + 1 4 B [ cos ⁡ σ ( − 1 + 2 cos 2 ⁡ ( 2 σ m ) ) − 1 6 B cos ⁡ ( 2 σ m ) ( − 3 + 4 sin 2 ⁡ σ ) ( − 3 + 4 cos 2 ⁡ ( 2 σ m ) ) ] } {\displaystyle \Delta \sigma =B\sin \sigma {\Big \{}\cos(2\sigma _{m})+{\tfrac {1}{4}}B{\big [}\cos \sigma {\big (}-1+2\cos ^{2}(2\sigma _{m}){\big )}-{\tfrac {1}{6}}B\cos(2\sigma _{m})(-3+4\sin ^{2}\sigma ){\big (}-3+4\cos ^{2}(2\sigma _{m}){\big )}{\big ]}{\Big \}}}
s = b A ( σ − Δ σ ) {\displaystyle s=bA(\sigma -\Delta \sigma )\,}
α 1 = arctan ⁡ ( cos ⁡ U 2 sin ⁡ λ cos ⁡ U 1 sin ⁡ U 2 − sin ⁡ U 1 cos ⁡ U 2 cos ⁡ λ ) {\displaystyle \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)} [4]
α 2 = arctan ⁡ ( cos ⁡ U 1 sin ⁡ λ − sin ⁡ U 1 cos ⁡ U 2 + cos ⁡ U 1 sin ⁡ U 2 cos ⁡ λ ) {\displaystyle \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)} [4]

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

順解法

順解法は始点 ( ϕ 1 , L 1 ) {\displaystyle (\phi _{1},L_{1})} と始点における方位角 α 1 {\displaystyle \alpha _{1}} 、距離 s {\displaystyle s} が与えられた時の終点の座標 ( ϕ 2 , L 2 ) {\displaystyle (\phi _{2},L_{2})} と方位角 α 2 {\displaystyle \alpha _{2}} を導く。

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

U 1 = arctan ⁡ ( ( 1 − f ) tan ⁡ ϕ 1 ) {\displaystyle U_{1}=\arctan \left((1-f)\tan \phi _{1}\right)\,}
σ 1 = arctan ⁡ ( tan ⁡ U 1 cos ⁡ α 1 ) {\displaystyle \sigma _{1}=\arctan \left({\frac {\tan U_{1}}{\cos \alpha _{1}}}\right)\,} [4]
sin ⁡ α = cos ⁡ U 1 sin ⁡ α 1 {\displaystyle \sin \alpha =\cos U_{1}\sin \alpha _{1}\,}
cos 2 ⁡ α = 1 − sin 2 ⁡ α {\displaystyle \cos ^{2}\alpha =1-\sin ^{2}\alpha \,}
u 2 = cos 2 ⁡ α a 2 − b 2 b 2 {\displaystyle u^{2}=\cos ^{2}\alpha {\frac {a^{2}-b^{2}}{b^{2}}}\,}
A = 1 + u 2 16384 { 4096 + u 2 [ − 768 + u 2 ( 320 − 175 u 2 ) ] } {\displaystyle A=1+{\frac {u^{2}}{16384}}\left\{4096+u^{2}\left[-768+u^{2}(320-175u^{2})\right]\right\}}
B = u 2 1024 { 256 + u 2 [ − 128 + u 2 ( 74 − 47 u 2 ) ] } {\displaystyle B={\frac {u^{2}}{1024}}\left\{256+u^{2}\left[-128+u^{2}(74-47u^{2})\right]\right\}}

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

2 σ m = 2 σ 1 + σ {\displaystyle 2\sigma _{m}=2\sigma _{1}+\sigma \,}
Δ σ = B sin ⁡ σ { cos ⁡ ( 2 σ m ) + 1 4 B [ cos ⁡ σ ( − 1 + 2 cos 2 ⁡ ( 2 σ m ) ) − 1 6 B cos ⁡ ( 2 σ m ) ( − 3 + 4 sin 2 ⁡ σ ) ( − 3 + 4 cos 2 ⁡ ( 2 σ m ) ) ] } {\displaystyle \Delta \sigma =B\sin \sigma {\Big \{}\cos(2\sigma _{m})+{\tfrac {1}{4}}B{\big [}\cos \sigma {\big (}-1+2\cos ^{2}(2\sigma _{m}){\big )}-{\tfrac {1}{6}}B\cos(2\sigma _{m})(-3+4\sin ^{2}\sigma ){\big (}-3+4\cos ^{2}(2\sigma _{m}){\big )}{\big ]}{\Big \}}}
σ = s b A + Δ σ {\displaystyle \sigma ={\frac {s}{bA}}+\Delta \sigma \,}

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

ϕ 2 = arctan ⁡ ( sin ⁡ U 1 cos ⁡ σ + cos ⁡ U 1 sin ⁡ σ cos ⁡ α 1 ( 1 − f ) sin 2 ⁡ α + ( sin ⁡ U 1 sin ⁡ σ − cos ⁡ U 1 cos ⁡ σ cos ⁡ α 1 ) 2 ) {\displaystyle \phi _{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)\,} [4]
λ = arctan ⁡ ( sin ⁡ σ sin ⁡ α 1 cos ⁡ U 1 cos ⁡ σ − sin ⁡ U 1 sin ⁡ σ cos ⁡ α 1 ) {\displaystyle \lambda =\arctan \left({\frac {\sin \sigma \sin \alpha _{1}}{\cos U_{1}\cos \sigma -\sin U_{1}\sin \sigma \cos \alpha _{1}}}\right)\,} [4]
C = f 16 cos 2 ⁡ α [ 4 + f ( 4 − 3 cos 2 ⁡ α ) ] {\displaystyle C={\frac {f}{16}}\cos ^{2}\alpha {\big [}4+f(4-3\cos ^{2}\alpha ){\big ]}\,}
L = λ − ( 1 − C ) f sin ⁡ α { σ + C sin ⁡ σ [ cos ⁡ ( 2 σ m ) + C cos ⁡ σ ( − 1 + 2 cos 2 ⁡ ( 2 σ m ) ) ] } {\displaystyle L=\lambda -(1-C)f\sin \alpha \left\{\sigma +C\sin \sigma \left[\cos(2\sigma _{m})+C\cos \sigma (-1+2\cos ^{2}(2\sigma _{m}))\right]\right\}\,}
L 2 = L + L 1 {\displaystyle L_{2}=L+L_{1}\,}
α 2 = arctan ⁡ ( sin ⁡ α − sin ⁡ U 1 sin ⁡ σ + cos ⁡ U 1 cos ⁡ σ cos ⁡ α 1 ) {\displaystyle \alpha _{2}=\arctan \left({\frac {\sin \alpha }{-\sin U_{1}\sin \sigma +\cos U_{1}\cos \sigma \cos \alpha _{1}}}\right)\,} [4]

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

Vincentyによる修正

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

A = 1 + 1 4 ( k 1 ) 2 1 − k 1 {\displaystyle A={\frac {1+{\frac {1}{4}}(k_{1})^{2}}{1-k_{1}}}}
B = k 1 ( 1 − 3 8 ( k 1 ) 2 ) {\displaystyle B=k_{1}(1-{\tfrac {3}{8}}(k_{1})^{2})}

ここで、 k 1 = ( 1 + u 2 ) − 1 ( 1 + u 2 ) + 1 {\displaystyle k_{1}={\frac {{\sqrt {(1+u^{2})}}-1}{{\sqrt {(1+u^{2})}}+1}}} である。

対蹠付近の点

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

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

関連項目

注釈

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

参考文献

  • Bessel, F. W. (2010). “The calculation of longitude and latitude from geodesic measurements (1825)”. Astron. Nachr. 331 (8): 852–861. arXiv:0908.1824. doi:10.1002/asna.201011352. English translation of Astron. Nachr. 4, 241–254 (1825). 
  • Helmert, F. R. (1964). Mathematical and Physical Theories of Higher Geodesy, Part 1 (1880). St. Louis: Aeronautical Chart and Information Center. http://geographiclib.sf.net/geodesic-papers/helmert80-en.html 2011年7月30日閲覧。. English translation of Die Mathematischen und Physikalischen Theorieen der Höheren Geodäsie, Vol. 1 (Teubner, Leipzig, 1880). 
  • Karney, C. F. F. (2013). “Algorithms for geodesics”. Journal of Geodesy 87 (1): 43–42. arXiv:1109.4448. Bibcode: 2013JGeod..87...43K. doi:10.1007/s00190-012-0578-z (open access). Addenda. 
  • Legendre, A. M. (1806). “Analyse des triangles tracės sur la surface d'un sphėroïde”. Mém. de l'Inst. Nat. de France (1st sem.): 130–161. https://books.google.co.jp/books?id=-d0EAAAAQAAJ&pg=PA130-IA4&redir_esc=y&hl=ja 2011年7月30日閲覧。. 
  • Rainsford, H. F. (1955). “Long geodesics on the ellipsoid”. Bulletin géodésique 37: 12–22. Bibcode: 1955BGeod..29...12R. doi:10.1007/BF02527187. 
  • Rapp, R. H. (March 1993). Geometric Geodesy, Part II (Technical report). Ohio State University. 2011年8月1日閲覧。
  • Vincenty, T. (April 1975a). “Direct and Inverse Solutions of Geodesics on the Ellipsoid with application of nested equations”. Survey Review XXIII (misprinted as XXII) (176): 88–93. doi:10.1179/sre.1975.23.176.88. http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf 2009年7月11日閲覧。. 
  • Vincenty, T. (April 1976). “Correspondence”. Survey Review XXIII (180): 294. 
  • Vincenty, T. (August 1975b). Geodetic inverse solution between antipodal points (Technical report). DMAAC Geodetic Survey Squadron. doi:10.5281/zenodo.32999。
  • (PDF) Geocentric Datum of Australia (GDA) Reference Manual. Intergovernmental committee on survey and mapping (ICSM). (February 2006). ISBN 0-9579951-0-5. http://www.icsm.gov.au/gda/gdatm/index.html 2009年7月11日閲覧。 [リンク切れ]

外部リンク

  • 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.