Vincenty法 （英 : Vincenty's formulae ）とは回転楕円体 上の測地法 の反復計算アルゴリズム である。サディアス・ヴィンセンティ （英語版 ）
地球形状は回転楕円体 として近似可能であり、Vincenty法の逆解法は、その面上で二点間の最小距離 （測地線 の長さ）を計算する。ただしVincenty法の計算誤差は[1]
近年は 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法では収束しないこともありうるが、コードサイズは小さいものとなった。
逆解法は回転楕円体上の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 }
-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による修正
編集
対蹠付近の点
編集
先に触れたように、逆解法の反復計算は対蹠点付近において、収束しなかったり、収束が遅かったりする。例えば、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 })}
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 })}
関連項目
編集
^ 地球を球体 として近似すると大円距離 の計算となるが、それよりは精度が良い。
^ Karney (2013) の方法ではニュートン法 を用い、あらゆる入力点に対して収束が速く、倍精度浮動小数点数 限界精度に達する。
^
σ
{\displaystyle \sigma }
sin
σ
,
cos
σ
{\displaystyle \sin \sigma ,\cos \sigma }
^ a b c d e f g arctanの値は2変数関数であるatan2
^
sin
σ
=
0
{\displaystyle \sin \sigma =0}
sin
α
{\displaystyle \sin \alpha }
対蹠点 であることを意味する。^ 2点がともに赤道上にあるとき、
U
=
0
{\displaystyle U=0}
cos
(
2
σ
m
)
{\displaystyle \cos(2\sigma _{m})}
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) . http://geographiclib.sf.net/geodesic-papers/helmert80-en.html 2011年7月30日 閲覧。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 . 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 ISBN 0-9579951-0-5 . http://www.icsm.gov.au/gda/gdatm/index.html 2009年7月11日 閲覧。 [リンク切れ 外部リンク
編集
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 . 
![{\displaystyle U_{1}=\arctan[(1-f)\tan \phi _{1}]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/672aa17db392c77a316983d9d698bd0496902796)
![{\displaystyle U_{2}=\arctan[(1-f)\tan \phi _{2}]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/048ee561b37307d44ab1f65cf6c17a0921b45c5f)
![{\displaystyle C={\frac {f}{16}}\cos ^{2}\alpha {\big [}4+f(4-3\cos ^{2}\alpha ){\big ]}\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f031265d4bd2496585dc5e8e0e07c785aa9f6c77)
![{\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\}\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/5cd81f22ec10cb045321efd48cfb4257c695dfad)
![{\displaystyle A=1+{\frac {u^{2}}{16384}}\left\{4096+u^{2}\left[-768+u^{2}(320-175u^{2})\right]\right\}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/903dc7be839ff4bd7b9a13a5d892a87b8e521e5f)
![{\displaystyle B={\frac {u^{2}}{1024}}\left\{256+u^{2}\left[-128+u^{2}(74-47u^{2})\right]\right\}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ec4796e5ba3bab3a3829919e1a8d282d75e66eab)
![{\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 \}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c859929310e12cace30dda7ae7b7a43e1fcf98b1)
![{\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\}\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/758d94ec8122294d894b17d802c3cae89fa1e960)
