Vincenty法のソースを表示

'''Vincenty法'''（ヴィンセンティほう、{{lang-en-short|Vincenty's formulae}}）とは[[回転楕円体]]上の[[測地学|測地法]]の[[反復法 (数値計算)|反復計算アルゴリズム]]である。{{仮リンク|サディアス・ヴィンセンティ|en|Thaddeus Vincenty}} (1975a)によって考案された。

* 地球形状は[[回転楕円体]]として近似可能であり、Vincenty法の逆解法は、その面上で二点間の最小[[距離]]（[[測地線]]の長さ）を計算する。ただしVincenty法の計算誤差は<ref>地球を[[球体]]として近似すると[[大円距離]]の計算となるが、それよりは精度が良い。</ref>、与える入力点によっては 0.5&nbsp;mm(0.020″)に達する場合もある。
* 近年は Karney (2013) による計算方法がVincenty法よりも計算精度・計算速度ともに高く<ref>Karney (2013) の方法では[[ニュートン法]]を用い、あらゆる入力点に対して収束が速く、[[倍精度浮動小数点数]]限界精度に達する。</ref>、多く用いられる。
<!--これらは[[地球楕円体]]において0.5&nbsp;mm(0.020″)以下という誤差の小ささから広く用いられる。-->

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

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

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

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

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

== 表記法 ==
以下のように定義する。
; <math>a</math>:長軸半径 (赤道半径) (6378137.06&nbsp;m ([[WGS84]]) )
; <math>f</math>: [[扁平率]] (1/298.257223563 (WGS84) )
; <math>b=(1-f)a</math>: 短軸半径 (極半径) ( 6356752.314245&nbsp;m (WGS84) )
; <math>\phi_1,\phi_2</math>: 各点の[[緯度]]
; <math>U_1=\arctan [(1-f) \tan \phi_1]</math> <br/><math>U_2=\arctan [(1-f) \tan \phi_2]</math>: [[緯度#更成緯度 (reduced latitude)|更成緯度]] (補助球上の緯度)
; <math>L=L_2-L_1</math>: 2点間の[[経度]]差
; <math>\lambda_1, \lambda_2</math>: 補助球上の経度
; <math>\alpha_1, \alpha_2</math>: 各点における方位角
; <math>\alpha</math>: 赤道上での方位角
; <math>s</math> : 2点間の楕円体上の距離
; <math>\sigma</math> : 補助球上の弧の長さ

== 逆解法 ==
逆解法は回転楕円体上の2点の座標<math>(\phi_1, L_1), (\phi_2, L_2)</math>が与えられた時に、方位角<math>\alpha_1, \alpha_2</math>と距離<math>s</math>を導く。

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

::<math>\sin \sigma = \sqrt{ (\cos U_2 \sin \lambda)^2 + (\cos U_1 \sin U_2 - \sin U_1 \cos U_2 \cos \lambda)^2}</math>
::<math>\cos \sigma = \sin U_1 \sin U_2 + \cos U_1 \cos U_2 \cos \lambda \,</math>
::<math>\sigma = \arctan\frac{\sin\sigma}{\cos\sigma}\,</math><ref><math>\sigma</math>は赤道及び極付近での精度を保つために<math>\sin \sigma, \cos \sigma</math>から直接計算は行わない。</ref><ref name="atan">arctanの値は2変数関数である<code>[[atan2]]</code>等によって実装されるべきである。</ref>
::<math>\sin \alpha = \frac{\cos U_1 \cos U_2 \sin \lambda}{\sin \sigma} \,</math><ref name="opposite"><math>\sin \sigma=0</math>の時、<math>\sin \alpha</math>が不定となる。これは2点が同じ位置もしくは[[対蹠点]]であることを意味する。</ref>
::<math>\cos^2 \alpha = 1 - \sin^2 \alpha \,</math>
::<math>\cos (2 \sigma_m) = \cos \sigma - \frac{2 \sin U_1\sin U_2}{\cos^2 \alpha} \,</math><ref name="equator">2点がともに赤道上にあるとき、<math>U=0</math>となるので、<math>\cos (2 \sigma_m) </math>は使用されない。The limiting value is <math>\cos (2 \sigma_m) = -1</math>. <!--よくわからない --></ref>
::<math>C = \frac{f}{16} \cos^2 \alpha \big[4 + f(4-3 \cos^2 \alpha) \big] \,</math>
::<math>\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\} \, </math>

<math>\lambda</math>が所望の精度まで収束したら（10<sup>-12</sup>の偏差ならば0.06&nbsp;mmの精度）以下の計算を行う。

:<math>u^2 = \cos^2 \alpha \frac{a^2 - b^2}{b^2} \,</math>
:<math>A = 1 + \frac{u^2}{16384} \left\{ 4096 + u^2 \left[ -768 +u^2 (320 - 175u^2) \right] \right\}</math>
:<math>B = \frac{u^2}{1024} \left\{ 256 + u^2 \left[ -128 + u^2 (74-47 u^2) \right] \right\} </math>
:<math> \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\} </math>

:<math> s = b A(\sigma - \Delta \sigma) \,</math>
:<math> \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) </math><ref name="atan"/>
:<math> \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) </math><ref name="atan"/>

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

== 順解法 ==
順解法は始点<math>(\phi_1, L_1)</math>と始点における方位角<math>\alpha_1</math>、距離<math>s</math>が与えられた時の終点の座標<math>(\phi_2, L_2)</math>と方位角<math>\alpha_2</math>を導く。

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

:<math> U_1 = \arctan \left ((1 - f)\tan \phi_1 \right )\, </math>
:<math> \sigma_1 = \arctan \left ( \frac{ \tan U_1}{ \cos \alpha_1} \right ) \, </math><ref name="atan"/>
:<math> \sin \alpha = \cos U_1 \sin \alpha_1 \, </math>
:<math> \cos^2 \alpha = 1 - \sin^2 \alpha \, </math>
:<math> u^2 = \cos^2 \alpha \frac{a^2 - b^2}{b^2} \, </math>
:<math> A = 1 + \frac{u^2}{16384} \left\{ 4096 + u^2 \left[ -768 +u^2 (320 - 175u^2) \right] \right\} </math>
:<math> B = \frac{u^2}{1024} \left\{ 256 + u^2 \left[ -128 + u^2 (74-47 u^2) \right] \right\} </math>

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

::<math> 2 \sigma_m = 2 \sigma_1 + \sigma \, </math>
::<math> \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\} </math>
::<math> \sigma = \frac{s}{bA} + \Delta \sigma \, </math>

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

:<math> \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) \, </math><ref name="atan"/>
:<math> \lambda = \arctan \left( \frac{\sin \sigma \sin \alpha_1}{\cos U_1 \cos \sigma - \sin U_1 \sin \sigma \cos \alpha_1} \right) \, </math><ref name="atan"/>
:<math> C = \frac{f}{16} \cos^2 \alpha \big[4 + f(4-3 \cos^2 \alpha) \big] \, </math>
:<math> 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\} \, </math>
:<math> L_2 = L + L_1 \, </math>
:<math> \alpha_2 = \arctan \left( \frac{\sin \alpha}{-\sin U_1 \sin \sigma + \cos U_1 \cos \sigma \cos \alpha_1} \right) \, </math><ref name="atan"/>

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

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

:<math>A = \frac {1 + \frac {1}{4} (k_1)^2}{1 - k_1}</math>
:<math>B = k_1(1 - \tfrac {3}{8}(k_1)^2)</math>

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

== 対蹠付近の点 ==
先に触れたように、逆解法の反復計算は対蹠点付近において、収束しなかったり、収束が遅かったりする。例えば、WGS84測地系において<math>(\phi_1, L_1)=(0^\circ, 0^\circ),\ (\phi_2, L_2)=(0.5^\circ, 179.5^\circ)</math>であるとする。この計算は1&nbsp;mmの精度に達するまでに130回の計算が必要になる。逆解法がどのように実装されるかによって、計算結果は正しく19936288.579&nbsp;mを返したり、間違った結果を返したり、エラーとなったりする。例えば[http://www.ngs.noaa.gov/TOOLS/Inv_Fwd/Inv_Fwd.html NGS online utility]による計算結果は5&nbsp;kmも長い。Vincentyが提案した方法はこのような場合に収束を速める(Rapp, 1973)。

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

== 関連項目 ==
*[[経緯度]]
*[[大円距離]]
*[[子午線弧]]
*[[測地法]]
* [[測地線#回転楕円体面上の測地線]]

== 注釈 ==
{{Reflist}}

==参考文献==
* {{cite journal
 |first=F. W. |last=Bessel |authorlink=Friedrich Bessel
 |title=The calculation of longitude and latitude from geodesic measurements (1825)
 |journal=Astron. Nachr.
 |year=2010
 |volume=331 |issue=8 |pages=852&ndash;861
 |doi=10.1002/asna.201011352
 |arxiv=0908.1824 |postscript=. English translation of Astron. Nachr. '''4''', 241&ndash;254 (1825).
}}
* {{cite book
 |first=F. R. |last=Helmert |authorlink=Helmert
 |title=Mathematical and Physical Theories of Higher Geodesy, Part 1 (1880)
 |publisher=Aeronautical Chart and Information Center
 |year=1964
 |location=St. Louis
 |url=http://geographiclib.sf.net/geodesic-papers/helmert80-en.html
 |accessdate=2011-07-30
 |postscript=. English translation of ''Die Mathematischen und Physikalischen Theorieen der Höheren Geodäsie'', Vol. 1 (Teubner, Leipzig, 1880).
}}
* {{Cite journal | ref = harv| last1 = Karney | first1 = C. F. F. | doi = 10.1007/s00190-012-0578-z | title = Algorithms for geodesics | journal = Journal of Geodesy | volume = 87 | pages = 43–42 | year = 2013| issue = 1| pmid =  | pmc = | postscript = &nbsp;(open access). [http://geographiclib.sf.net/geod-addenda.html Addenda].|arxiv = 1109.4448 |bibcode = 2013JGeod..87...43K }}
* {{cite journal
 |first=A. M. |last=Legendre |authorlink=Adrien-Marie Legendre
 |title=Analyse des triangles tracės sur la surface d'un sphėroïde
 |journal=Mém. de l'Inst. Nat. de France
 |year=1806
 |issue=1st sem. |pages=130&ndash;161
 |url=https://books.google.co.jp/books?id=-d0EAAAAQAAJ&pg=PA130-IA4&redir_esc=y&hl=ja
 |accessdate=2011-07-30
}}
* {{Cite journal | ref = harv| last1 = Rainsford | first1 = H. F. | doi = 10.1007/BF02527187 | title = Long geodesics on the ellipsoid | journal = Bulletin géodésique | volume = 37 | pages = 12–22 | year = 1955 | pmid =  | pmc = | bibcode = 1955BGeod..29...12R}}
* {{cite techreport
 |first=R. H. |last=Rapp
 |title=Geometric Geodesy, Part II
 |institution=Ohio State University
 |date=March 1993
 |url=https://hdl.handle.net/1811/24409
 |accessdate=2011-08-01
}}
* {{cite journal
 |first=T. |last=Vincenty |authorlink=Thaddeus Vincenty
 |title=Direct and Inverse Solutions of Geodesics on the Ellipsoid with application of nested equations
 |journal=Survey Review
 |volume=XXIII (misprinted as XXII) |issue=176 |date=April 1975a |pages=88&ndash;93
 |doi = 10.1179/sre.1975.23.176.88
 |url=http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf |accessdate=2009-07-11
}}
* {{cite journal
 |first=T. |last=Vincenty |authorlink=Thaddeus Vincenty
 |title=Correspondence
 |journal=Survey Review
 |volume=XXIII |issue=180 |date=April 1976 |pages=294
}}
* {{cite techreport
 |first=T. |last=Vincenty |authorlink=Thaddeus Vincenty
 |title=Geodetic inverse solution between antipodal points
 |institution=DMAAC Geodetic Survey Squadron
 |date=August 1975b
 |doi = 10.5281/zenodo.32999
}}
* {{cite book
 |publisher=Intergovernmental committee on survey and mapping (ICSM)
 |date=February 2006
 |isbn=0-9579951-0-5
 |title=Geocentric Datum of Australia (GDA) Reference Manual
 |url=http://www.icsm.gov.au/gda/gdatm/index.html
 |format=PDF
 |accessdate=2009-07-11
}}{{リンク切れ|date=January 2018}}

==外部リンク==
* Online calculators from Geoscience Australia:
** [http://www.ga.gov.au/geodesy/datums/vincenty_direct.jsp Vincenty Direct] (destination point)
** [http://www.ga.gov.au/geodesy/datums/vincenty_inverse.jsp Vincenty Inverse] (distance between points)
* Calculators from the U.S. National Geodetic Survey:
** [http://www.ngs.noaa.gov/TOOLS/Inv_Fwd/Inv_Fwd.html 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):
** [http://www.movable-type.co.uk/scripts/latlong-vincenty-direct.html Vincenty Direct] (destination point)
** [http://www.movable-type.co.uk/scripts/latlong-vincenty.html Vincenty Inverse] (distance between points)
* [http://geographiclib.sourceforge.net/ 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&nbsp;nm) and the inverse solution is complete. Here is an [http://geographiclib.sourceforge.net/cgi-bin/GeodSolve online version of GeodSolve].

{{DEFAULTSORT:ういんせんていほう}}
[[Category:測量]]
[[Category:アルゴリズム]]
[[Category:エポニム]]