大円距離

大円距離(だいえんきょり、英: great-circular distance）は、球面上の大円に沿う距離を指す。大円の弧の性質により、球面上の経路として2点間の最短距離である。

特に地球上においては大圏距離（たいけんきょり）ともいう。大円の弧（大圏コース）は等角航路やisoazimuthal線と同様に地球上の任意の2点間を結ぶことができる3つの手法の一つであり、航空機や船舶の経路計算で使われる。

なおこの記事では回転楕円体面上の最短距離（測地距離）は扱わない。

概要

ユークリッド空間では、球内部を通り2点間を直線で結ぶユークリッド距離が最小となるが（したがって大円弦長）、球面上には直線が存在しないためこれとは異なる。非ユークリッド空間では、直線を一般化した測地線を使用する。球面においては測地線は球の中心を中心とする円である大円となるため、大円距離は大円上の2点間の弧の長さとなる。

球面上の対蹠点以外の2点を通る大円は一意に定まる。 2点は大円を2つの弧に分割する。そのうち短い方の弧の長さが大円距離となる。

対蹠点に関しては、その2点を通る任意の円が大円となるが、全ての円において2点間の弧の長さは一定である。すなわち半円の円周であり、半径 $r$ の球においては $\pi r$ である。

地球はほぼ球状であるため、2点間の距離を球として計算しても誤差は0.5%以内となる（後述）。^[1]

数式

２点間の大円距離 $d$ （球面上の経路として距離）は、両者間の中心角である弧度 $\Delta \sigma$ および球の半径 $r$ から求められる。

d=r\,\Delta \sigma

また距離以外に、出発点および中間の各点における方位角の計算も行う。海里の距離の計算の際は度数法における分がそのまま海里として用いられる（度で表したものの60倍）。

球面余弦定理

弧度 $\Delta \sigma$ は球面余弦定理を用いて次のように計算できる。ここで $\phi _{1},\lambda _{1}$ と $\phi _{2},\lambda _{2}$ をそれぞれ点1と点2の緯度と経度とし、 $\Delta \phi ,\Delta \lambda$ はそれらの差を表す。

\Delta \sigma =\arccos {\bigl (}\sin \phi _{1}\sin \phi _{2}+\cos \phi _{1}\cos \phi _{2}\cos \Delta \lambda {\bigr )}

ただし歴史上、この $\arccos$ 関数の引数が1近傍（すなわち近距離条件）もしくは-1近傍（すなわち対蹠点近傍条件）の条件では計算が扱いにくいとされた^{[注釈 1]}。

計算機による計算精度

ただし現在用いられる64ビットの浮動小数点数であれば、これらの条件でも、最大丸め誤差は約1mmであり実用上問題を起こさないと考えられている^[2]。

ただし精度（ビット数）が低い浮動小数点数を扱う計算機においては、大きな丸め誤差が発生し精度条件が悪化する。

Haversine公式

歴史的に上記の球面余弦定理の利用の欠点を避けるために、Haversine 関数^{[注釈 2]}を用いる計算式（半正矢公式（英語版））が、haversineの関数表を使って多用された^[3]。

\Delta \sigma =\operatorname {archav} {\bigl (}\operatorname {hav} \Delta \phi +\cos \phi _{1}\cos \phi _{2}\operatorname {hav} \Delta \lambda {\bigr )}

これは下記の $\sin ,\ \arcsin$ を用いる弦長計算を経由する式と等価である。

ただし浮動小数点数計算では $\operatorname {archav}$ の引数が1を超えないよう注意が必要である^{[注釈 3]}（対蹠点間で発生するおそれがある）。

単位球面における大円の弦長と弧度

単位球面における大円の弦長 $\Delta \sigma _{\text{c}}$ と弧度 $\Delta \sigma$ との関係は^{[注釈 4]}、

\Delta \sigma _{\text{c}}=2\sin {\frac {\Delta \sigma }{2}}

\Delta \sigma =2\arcsin {\frac {\Delta \sigma _{\text{c}}}{2}}

ただし上記と同様に浮動小数点数計算では $\arcsin$ の引数が1を超えないよう注意が必要である^{[注釈 5]}。

弦長の計算

弦長 $\Delta \sigma _{\text{c}}$ の計算式は（ $\phi _{\textrm {m}}\triangleq {\frac {\phi _{1}+\phi _{2}}{2}}$ )：

{\begin{aligned}\Delta \sigma _{\text{c}}&=2{\sqrt {\sin ^{2}\left({\frac {\Delta \phi }{2}}\right)+\cos {\phi _{1}}\cdot \cos {\phi _{2}}\cdot \sin ^{2}\left({\frac {\Delta \lambda }{2}}\right)}}\ ,\\&=2{\sqrt {\left(\sin {\frac {\Delta \lambda }{2}}\cos \phi _{\textrm {m}}\right)^{2}+\left(\cos {\frac {\Delta \lambda }{2}}\sin {\frac {\Delta \phi }{2}}\right)^{2}}}\ .\end{aligned}}

逆正接関数を用いる計算式

上記のHaversine型計算式は球面上のほとんどの点の間において正確だが、対蹠点間に対しては $\arcsin$ 関数の計算精度が低下する（ただし最大丸め誤差は約1mmであり実用上問題を起こさないと考えられている）。対策として代わりに逆正接関数 $\arctan$ ^{[注釈 6]}を用いて下記のように書き換えることができる。

\Delta \sigma =2\arcsin y=2\arctan {\frac {y}{\sqrt {1-y^{2}}}},\ y\triangleq {\frac {\Delta \sigma _{\text{c}}}{2}}.

また同様の目的だが、弦長 $\Delta \sigma _{\text{c}}$ を経由しない形の下記の式がある^{[注釈 7]}。

\Delta \sigma =\arctan {\frac {\sin \Delta \sigma }{\cos \Delta \sigma }},

\cos \Delta \sigma =\sin \phi _{1}\cdot \sin \phi _{2}+\cos \phi _{1}\cdot \cos \phi _{2}\cdot \cos \Delta \lambda ,

\sin \Delta \sigma ={\sqrt {\left(\cos \phi _{2}\cdot \sin \Delta \lambda \right)^{2}+\left(\cos \phi _{1}\cdot \sin \phi _{2}-\sin \phi _{1}\cdot \cos \phi _{2}\cdot \cos \Delta \lambda \right)^{2}}}.

直交座標系による弦長の計算

単位球面における大円弦長 $\Delta \sigma _{\text{c}}$ は直交座標系を経由しても計算できる。

{\begin{aligned}\Delta \sigma _{\text{c}}&={\sqrt {(\Delta {X})^{2}+(\Delta {Y})^{2}+(\Delta {Z})^{2}}}\\\Delta {X}&\triangleq \cos \phi _{2}\cdot \cos \lambda _{2}-\cos \phi _{1}\cdot \cos \lambda _{1}\\\Delta {Y}&\triangleq \cos \phi _{2}\cdot \sin \lambda _{2}-\cos \phi _{1}\cdot \sin \lambda _{1}\\\Delta {Z}&\triangleq \sin \phi _{2}-\sin \phi _{1}\end{aligned}}

ベクトル表現

弧度 $\Delta \sigma$ は、ユークリッド空間上の単位長3次元ベクトルの内積と外積により以下のように表すことができる^[5]。 $\mathbf {n} _{1}$ と $\mathbf {n} _{2}$ は球面上の2点の単位法線ベクトルを表す。計算精度条件の比較は上記と同様となる。

{\begin{aligned}\Delta \sigma &=\arccos(\mathbf {n} _{1}\cdot \mathbf {n} _{2})\\\Delta \sigma &=\arcsin \left|\mathbf {n} _{1}\times \mathbf {n} _{2}\right|\\\Delta \sigma &=\arctan {\frac {\left|\mathbf {n} _{1}\times \mathbf {n} _{2}\right|}{\mathbf {n} _{1}\cdot \mathbf {n} _{2}}}\\\end{aligned}}\,\!

回転楕円体面上の距離

→詳細は「経緯度 § 回転楕円体面に沿う最短距離の式」を参照

地球半径

地球の形状は回転軸方向に潰れた球（扁球）とみなすことができる。このとき赤道半径 $a$ は6378.137 km、極半径 $b$ は6356.752 kmとなる。赤道付近の短い南北方向の線においては半径 $b^{2}/a$ (6335.439 km)とした際が最も良い近似となり、極においては半径 $a^{2}/b$ (6399.594 km) が最も良い。この差は1%である。つまり、地球を球体と仮定した計算においては、地球上の任意の2点間の距離に対する1つの計算による誤差は 0.5% 以内とすることができる（ただし、限られた地域に関してはより誤差の少ない値を使用することもできる。）。

この値として平均地球半径を用いるとよく^[6]、その値は $R_{1}={\frac {1}{3}}(2a+b)\approx 6371\,\mathrm {km}$ (値はGRS 80における回転楕円体近似に対して)である。扁平率が小さい場合はこの値が平均自乗誤差を最小化する。

脚注

[脚注の使い方]

注釈

^ 例えば地球上においては1 kmの距離に対して $\arccos$ 関数の引数は0.99999999となる。
^ $\operatorname {hav} (\theta )\triangleq \sin ^{2}(\theta /2),\ \operatorname {archav} x=2\arcsin {\sqrt {x}}$
^ $\operatorname {archav} 1=\pi$
^ $\Delta \sigma _{\text{c}}=0$ における級数展開形は、 $\Delta \sigma =\Delta \sigma _{\text{c}}\left(1+{\frac {\Delta \sigma _{\text{c}}^{2}}{24}}+{\frac {3\Delta \sigma _{\text{c}}^{4}}{640}}+\cdots \right)$ 。
^ $\arcsin 1={\frac {\pi }{2}}$
^ プログラミングの際は通常の逆正接関数atan(y/x)よりも atan2(y,x) 関数を用いたほうが全象限で出力されるため優れている。
^ 導出法の一つとしてVincenty formula^[4]内で全ての軸の経が等しい楕円体からも導かれる。