ヴィンセンティの公式

ヴィンセンティの公式は、測地学で回転楕円体の表面上の2点間の距離を計算するために使用される 2 つの関連する反復法であり、タデウス ヴィンセンティ(1975a)によって開発されました。これらは、地球の形状が扁平回転楕円体であるという仮定に基づいているため、大圏距離など、地球が球形であると仮定する方法よりも正確です。

最初の（直接）方法は、別の点から指定された距離と方位（方向）にある点の位置を計算します。 2 番目の（逆）方法は、指定された 2 つの点間の地理的な距離と方位を計算します。 これらの方法は、地球の楕円体 上で 0.5 mm（0.020 インチ）以内の精度であるため、測地学で広く使用されています。

ヴィンセンティの目標は、楕円体上の測地線に関する既存のアルゴリズムを、プログラムの長さが最小限になる形式で表現することだった (Vincenty 1975a)。彼の未発表のレポート (1975b) には、メモリが数キロバイトしかないWang 720 卓上計算機の使用について触れられている。長い直線に対して高い精度を得るために、このソリューションでは、補助球面に基づく、ルジャンドル (1806)、ベッセル (1825)、ヘルメルト (1880) の古典的なソリューションを使用している。ヴィンセンティは、レインズフォード (1955) が示したこの方法の定式化に依存していた。ルジャンドルは、地理緯度を縮約緯度にマッピングし、大円の方位角を測地線の方位角に等しく設定することで、楕円体測地線を補助球面上の大円に正確にマッピングできることを示した。楕円体上の経度と測地線に沿った距離は、球面上の経度と大円に沿った弧の長さによって、単純な積分で与えられます。ベッセルとヘルマートは、これらの積分に対して急速に収束する級数を与え、これにより測地線を任意の精度で計算できるようになりました。

プログラムのサイズを最小化するために、Vincenty はこれらの級数を取り、各級数の最初の項を小さなパラメータとして使用してそれらを再展開し、^{[明確化が必要]}に切り捨てました。これにより、経度と距離の積分のコンパクトな式が得られました。式は、単一の一時レジスタのみを使用して多項式を評価できるため、ホーナー(またはネストされた) 形式になりました。最後に、直接法と逆法で暗黙の方程式を解くために、単純な反復手法が使用されました。これらは低速ですが (逆法の場合は収束しないこともあります)、コード サイズの増加は最小限に抑えられます。 ${\displaystyle O(f^{3})}$

次の表記を定義します。

2点の座標( Φ1、  L1 )と( Φ2、  L2 )_が与えられた場合、逆問題_では方位角_α1、_α2と楕円体距離_sを求めます。

U ₁、U ₂、Lを計算し、λ = Lの初期値を設定します。次に、 λ が収束するまで次の式を繰り返し評価します。

${\displaystyle \sin \sigma ={\sqrt {\left(\cos U_{2}\sin \lambda \right)^{2}+\left(\cos U_{1}\sin U_{2}-\sin U_{1}\cos U_{2}\cos \lambda \right)^{2}}}}$

${\displaystyle \cos \sigma =\sin U_{1}\sin U_{2}+\cos U_{1}\cos U_{2}\cos \lambda \,}$

${\displaystyle \sigma =\operatorname {arctan2} \left(\sin \sigma ,\cos \sigma \right)}$^[1]

${\displaystyle \sin \alpha ={\frac {\cos U_{1}\cos U_{2}\sin \lambda }{\sin \sigma }}}$^[2]

${\displaystyle \cos^{2}\alpha =1-\sin^{2}\alpha }$

${\displaystyle \cos \left(2\sigma _{\text{m}}\right)=\cos \sigma -{\frac {2\sin U_{1}\sin U_{2}}{\cos ^{2}\alpha }}=\cos \sigma -{\frac {2\sin U_{1}\sin U_{2}}{1-\sin ^{2}\alpha }}}$^[3]

${\displaystyle C={\frac {f}{16}}\cos ^{2}\alpha \left[4+f\left(4-3\cos ^{2}\alpha \right)\right]}$

${\displaystyle \lambda =L+(1-C)f\sin \alpha \left\{\sigma +C\sin \sigma \left[\cos \left(2\sigma _{\text{m}}\right)+C\cos \sigma \left(-1+2\cos ^{2}\left(2\sigma _{\text{m}}\right)\right)\right]\right\}}$

λが所望の精度（10 ^{−12 は} 約 0.06 mm に相当）に収束したら、次のことを評価します。

${\displaystyle {\begin{aligned}u^{2}&=\cos ^{2}\alpha \left({\frac {a^{2}-b^{2}}{b^{2}}}\right)\\A&=1+{\frac {u^{2}}{16384}}\left(4096+u^{2}\left[-768+u^{2}\left(320-175u^{2}\right)\right]\right)\\B&={\frac {u^{2}}{1024}}\left(256+u^{2}\left[-128+u^{2}\left(74-47u^{2}\right)\right]\right)\\\Delta \sigma &=B\sin \sigma \left\{\cos(2\sigma _{\text{m}})+{\frac {1}{4}}B\left(\cos \sigma \left[-1+2\cos ^{2}\left(2\sigma _{\text{m}}\right)\right]-{\frac {1}{6}}B\cos \left[2\sigma _{\text{m}}\right]\left[-3+4\sin ^{2}\sigma \right]\left[-3+4\cos ^{2}\left(2\sigma _{\text{m}}\right)\right]\right)\right\}\\s&=bA(\sigma -\Delta \sigma )\,\\\alpha _{1}&=\operatorname {arctan2} \left(\cos U_{2}\sin \lambda ,\cos U_{1}\sin U_{2}-\sin U_{1}\cos U_{2}\cos \lambda \right)\\\alpha _{2}&=\operatorname {arctan2} \left(\cos U_{1}\sin \lambda ,-\sin U_{1}\cos U_{2}+\cos U_{1}\sin U_{2}\cos \lambda \right)\end{aligned}}}$

ほぼ対蹠的な 2 つの点の間では、反復式が収束しない場合があります。これは、上記の式で計算されたλの最初の推定値の絶対値がπより大きい場合に発生します。

初期点 ( Φ ₁、L ₁ ) と初期方位角α ₁、および測地線に沿った距離sが与えられている場合、問題は終点 ( Φ ₂、L ₂ ) と方位角α ₂を見つけることです。

まず、次の計算を行ってください。

${\displaystyle {\begin{aligned}U_{1}&=\arctan \left[(1-f)\tan \phi _{1}\right]\\\sigma _{1}&=\operatorname {arctan2} \left(\tan U_{1},\cos \alpha _{1}\right)\\\sin \alpha &=\cos U_{1}\sin \alpha _{1}\\u^{2}&=\cos ^{2}\alpha \left({\frac {a^{2}-b^{2}}{b^{2}}}\right)=\left(1-\sin ^{2}\alpha \right)\left({\frac {a^{2}-b^{2}}{b^{2}}}\right)\\A&=1+{\frac {u^{2}}{16384}}\left(4096+u^{2}\left[-768+u^{2}(320-175u^{2})\right]\right)\\B&={\frac {u^{2}}{1024}}\left(256+u^{2}\left[-128+u^{2}\left(74-47u^{2}\right)\right]\right)\end{aligned}}}$

次に、初期値を使用して、 σに大きな変化がなくなるまで次の式を繰り返します。 ${\displaystyle \sigma ={\tfrac {s}{bA}}}$

${\displaystyle {\begin{aligned}2\sigma _{\text{m}}&=2\sigma _{1}+\sigma \\\Delta \sigma &=B\sin \sigma \left\{\cos \left(2\sigma _{\text{m}}\right)+{\frac {1}{4}}B\left(\cos \sigma \left[-1+2\cos ^{2}\left(2\sigma _{\text{m}}\right)\right]-{\frac {1}{6}}B\cos \left[2\sigma _{\text{m}}\right]\left[-3+4\sin ^{2}\sigma \right]\left[-3+4\cos ^{2}\left(2\sigma _{\text{m}}\right)\right]\right)\right\}\\\sigma &={\frac {s}{bA}}+\Delta \sigma \end{aligned}}}$

十分な精度で σが得られたなら、次のように評価します。

${\displaystyle {\begin{aligned}\phi _{2}&=\operatorname {arctan2} \left(\sin U_{1}\cos \sigma +\cos U_{1}\sin \sigma \cos \alpha _{1},(1-f){\sqrt {\sin ^{2}\alpha +\left(\sin U_{1}\sin \sigma -\cos U_{1}\cos \sigma \cos \alpha _{1}\right)^{2}}}\right)\\\lambda &=\operatorname {arctan2} \left(\sin \sigma \sin \alpha _{1},\cos U_{1}\cos \sigma -\sin U_{1}\sin \sigma \cos \alpha _{1}\right)\\C&={\frac {f}{16}}\cos ^{2}\alpha \left[4+f\left(4-3\cos ^{2}\alpha \right)\right]\\L&=\lambda -(1-C)f\sin \alpha \left\{\sigma +C\sin \sigma \left(\cos \left[2\sigma _{\text{m}}\right]+C\cos \sigma \left[-1+2\cos ^{2}\left(2\sigma _{\text{m}}\right)\right]\right)\right\}\\L_{2}&=L+L_{1}\\\alpha _{2}&=\operatorname {arctan2} \left(\sin \alpha ,-\sin U_{1}\sin \sigma +\cos U_{1}\cos \sigma \cos \alpha _{1}\right)\end{aligned}}}$

初期点が北極または南極にある場合、最初の方程式は不確定です。初期方位角が真東または真西にある場合、2 番目の方程式は不確定です。標準の 2 引数アークタンジェントatan2関数を使用すると、通常、これらの値は正しく処理されます。^{[説明が必要]}

1976年にSurvey Reviewに宛てた手紙の中で、VincentyはAとBの級数表現を、ヘルマートの展開パラメータk _1を使ったより単純な式に置き換えることを提案した。

${\displaystyle A={\frac {1+{\frac {1}{4}}{k_{1}}^{2}}{1-k_{1}}}}$

${\displaystyle B={k_{1}}\left(1-{\frac {3}{8}}{k_{1}}^{2}\right)}$

どこ

${\displaystyle k_{1}={\frac {{\sqrt {1+u^{2}}}-1}{{\sqrt {1+u^{2}}}+1}}}$

上で述べたように、逆問題の反復解法は、ほぼ対蹠点に対しては収束しないか、収束が遅い。収束が遅い例としては、WGS84 楕円体の場合の ( Φ ₁、  L ₁ ) = (0°、 0°) および ( Φ ₂、  L ₂ ) = (0.5°、 179.5°) がある。この場合、1 mm の精度の結果を得るには約 130 回の反復が必要である。逆法の実装方法に応じて、アルゴリズムは正しい結果 (19936288.579 m)、誤った結果、またはエラー インジケーターを返す可能性がある。誤った結果の例としては、NGS オンライン ユーティリティが挙げられ、約 5 km 長い距離が返される。Vincenty は、このような場合に収束を加速する方法を提案している (Rapp、1993)。

逆法が収束しない例としては、WGS84 楕円体の場合、( Φ ₁ ,  L ₁ ) = (0°, 0°) および ( Φ ₂ ,  L ₂ ) = (0.5°, 179.7°) があります。未発表のレポートで、Vincenty (1975b) は、このようなケースを処理するための代替反復スキームを示しました。これは、約 60 回の反復後に正しい結果 19944127.421 m に収束しますが、他のケースでは何千回もの反復が必要になります。

Karney (2013) は、逆問題を 1 次元の根探索問題として再定式化しました。これは、入力点のすべてのペアに対してニュートン法で迅速に解くことができます。

• 地理的距離
• 大圏距離
• 子午線弧
• 楕円体上の測地線
• タデウス・ヴィンセンティ
• 測地学

• Geoscience Australiaのオンライン計算機:
  • ヴィンセンティ ダイレクト（目的地）
  • ヴィンセンティ逆数（点間の距離）
• 米国国立測地測量局の計算機:
  • 2 次元と 3 次元の両方で順方向 (直接) 問題と逆問題を含む、オンラインおよびダウンロード可能な PC 実行可能計算ユーティリティ (2011 年 8 月 1 日にアクセス)。
• Chris Veness による JavaScript ソース コード付きオンライン計算機 (Creative Commons Attribution ライセンス):
  • ヴィンセンティ ダイレクト（目的地）
  • ヴィンセンティ逆数（点間の距離）
• GeographicLib は、直接測地線問題および逆測地線問題を解くためのユーティリティ GeodSolve (MIT/X11 ライセンスのソース コード付き) を提供します。Vincenty と比較すると、これは約 1000 倍正確 (誤差 = 15 nm) で、逆解も完全です。こちらに GeodSolve のオンライン バージョンがあります。
• ソース コード付きの Vincenty の直接および逆数式の完全な実装、Tomasz Jastrzębski による Excel VBA 実装