Vincenty法

Vincenty法（英: Vincenty's formulae）とは回転楕円体上の2点間の距離（測地線の長さ）を計算する測地法の反復計算アルゴリズムである。サディアス・ヴィンセンティ（英語版） (1975a)によって考案された。 地球形状は回転楕円体として近似可能なので、地球面上で二点間距離の計算に用いられる。ただしVincenty法の計算誤差は 0.5 mm(0.020″)近くなる場合もある。なお地球を球体として近似する大円距離より精度が良い。 測地線の長さに関して、順解法(direct method)はある点からの距離と方位角により、他の点を求めるものである。逆解法(inverse method)は2点間の回転楕円体面上距離と方位角を求めるものである。 なお計算方法としては、近年は Karney (2013) による方法が計算精度・計算速度ともに高く、Vincenty法よりも多く用いられる。 Vincentyの目的は当時すでに存在していた回転楕円体上の測地法アルゴリズムのプログラムの長さを最小化することであった。彼の未発行のレポート(1975b)によれば数キロバイトのメモリしか搭載していないWang 720 電卓を使用した。長距離において良い精度を得るために、伝統的な補助球を用いた手法（Legendre (1806), Bessel (1825), and Helmert (1880)）を使用している。また、Vincentyの手法はRainsford (1955)による解法にも依存している。Legendreは回転楕円体上の測地線は緯度を更成緯度とし、方位角を同じものとすることで補助球の大円上に正確に対応させられることを示した。 楕円体上の経度と距離は単純な積分によって、補助球上の経度とその大円上の弧長から導かれる。BesselとHelmertはそれらの積分について収束の早い数列を与えた。それにより測地計算が任意の精度で行えるようになった。 プログラムサイズを最小化するために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 }}\,} sin ⁡ α = cos ⁡ U 1 cos ⁡ U 2 sin ⁡ λ sin ⁡ σ {\displaystyle \sin \alpha ={\frac {\cos U_{1}\cos U_{2}\sin \lambda }{\sin \sigma }}\,} 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 }}\,} 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)} α 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)} 両極付近の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)\,} 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)\,} λ = 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)\,} 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)\,} 始点が極点である時、最初の等式は不定となる。 始点の方位角が真西か真東の場合2番目の等式は不定となる。 しかし、2変数関数のatan2等を用いることでこれらの値は正しく計算できる。 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という正しい計算結果に達した。しかし、この手法は他の多くの場合において何千回もの反復が必要となる。 ニュートン法は全ての入力点に対して収束が速い(Karney, 2013)。 大円距離 子午線弧 測地法 測地線#回転楕円体面上の測地線