하버사인 공식(Haversine formula)

지구를 완전한 공으로 보면, 지구 위 두 점 사이의 가장 짧은 길은 지구 중심을 지나는 평면이 지표를 자른 원, 곧 대원⁠(great circle)⁠을 따라가는 측지선⁠(geodesic)⁠입니다(구면기하⁠, spherical geometry⁠). 그 길이는 지구 반지름 R(평균⁠(mean)⁠ 약 6,371 km)에 두 점이 지구 중심에서 이루는 각 θ를 곱한 $d = R\theta$입니다. 문제는 위도 φ와 경도 λ만 알 때 θ를 구하는 것이고, 하버사인 공식이 그 답입니다. 하버사인은 $\operatorname{hav}\theta = \sin^2(\theta/2) = (1 - \cos\theta)/2$로, 두 번째 등호는 배각 공식에서 나옵니다.

가로가 경도, 세로가 위도인 지도에서 두 점 A, B를 끌어 보세요. 서울–샌프란시스코 런던–시드니 서울에서 1 km 중심각⁠(central angle)⁠은 , 대원 거리는 입니다. 반면 지도 위의 직선을 지도 눈금대로(경도 1°와 위도 1°를 모두 적도에서의 길이로 쳐서) 재면 가 나옵니다.

지도 위의 직선은 지구 위에서 가장 짧은 길이 아니고, 지도 눈금으로 잰 길이도 참 거리가 아닙니다. 이 지도는 높은 위도일수록 가로를 늘려 그립니다. 위도 φ에서 경도 1°의 실제 길이는 적도의 $\cos\varphi$배뿐이기 때문입니다(코사인⁠, cosine⁠). 그래서 대원 항로는 지도에서 극 쪽으로 불룩하게 휘어 보입니다. 가우스가 보였듯 구면의 가우스 곡률⁠(Gaussian curvature)⁠이 0이 아니니, 모든 거리를 그대로 담는 평면 지도는 없습니다.

같은 θ는 구면 코사인 법칙⁠(law of cosines)⁠ $\cos\theta = \sin\varphi_1\sin\varphi_2 + \cos\varphi_1\cos\varphi_2\cos\Delta\lambda$로도 구할 수 있습니다. 식은 맞지만, 가까운 두 점에서는 $\cos\theta \approx 1 - \theta^2/2$가 1에 너무 가까워 앞자리 숫자들이 모두 9로 채워지고, 정작 필요한 차이는 뒷자리에서 잘려 나갑니다. 하버사인은 1에서 뺄 필요 없이 작은 수 $\theta^2/4$ 자체를 다루니 자릿수를 잃지 않습니다. 아래 그래프는 유효숫자가 7자리 남짓인 단정밀도⁠(single precision)⁠로 두 공식을 계산한 상대 오차입니다. 코사인 법칙은 1.5 km쯤 아래에서 완전히 무너집니다. 지금 두 점을 단정밀도로 계산하면 하버사인은 , 코사인 법칙은 입니다.

계산기가 없던 시절 항해사들도 같은 문제를 겪었습니다. 자릿수가 정해진 표를 찾아 계산했으니, 1에 가까운 코사인 값은 쓸모가 없었습니다. 게다가 하버사인은 음수가 되지 않아 부호를 따질 필요가 없고, 로그를 취해 곱셈을 덧셈으로 바꾸기도 좋았습니다. '하버사인'이라는 이름은 1835년 영국의 수학자이자 해군 교육자 제임스 인먼이 항해용 표를 펴내며 붙였다고 전해집니다. 컴퓨터 시대에도 같은 이유로 쓰입니다. 유효숫자가 16자리쯤인 배정밀도⁠(double precision)⁠에서도 코사인 법칙은 1 m 거리에서 방향에 따라 상대 오차가 0.4% 안팎까지 벌어지고, 10 cm 안팎부터는 아예 0을 내놓기 시작합니다.

이어지는 곳.

• $\operatorname{hav}\theta$는 사실 지구 속을 곧게 뚫은 현의 길이 c로 $(c/2R)^2$입니다. 현의 길이는 3차원 유클리드 거리⁠(Euclidean distance)⁠이니, 하버사인 공식은 피타고라스 정리⁠(Pythagorean theorem)⁠로 현을 재고 그것을 호의 길이로 바꾸는 셈입니다.
• 아주 가까운 두 점이라면 좁은 땅은 평평해 보이니(국소 선형성⁠, local linearity⁠), 경도 차에 $\cos\varphi$를 곱해 가로 거리로 바꾼 뒤 평면의 피타고라스 정리를 써도 충분합니다.
• 거꾸로 거의 정반대편인 두 점에서는 하버사인도 정밀도⁠(precision)⁠를 잃습니다. hav θ가 1에 가까워져, 이번에는 1과의 작은 차이가 뒷자리에서 잘려 나가기 때문입니다. 점 (x, y)의 각도를 사분면까지 가려 돌려주는 함수⁠(function)⁠ atan2로 $\theta = 2\operatorname{atan2}(\sqrt h, \sqrt{1-h})$(h는 hav θ)라 쓰는 흔한 형태는, 반올림 때문에 h가 1을 살짝 넘어 계산이 멈추는 일은 막아 주지만 잃은 자릿수를 되살리지는 못합니다. 어느 거리에서나 정밀하려면 두 점의 위치 벡터⁠(vector)⁠로 사잇각의 사인⁠(sine)⁠ 쪽 값(외적의 크기)과 코사인 쪽 값(내적⁠, dot product⁠)을 따로 구해 atan2에 함께 넣습니다.
• 지구는 적도 반지름이 극 반지름보다 약 21 km 긴, 살짝 납작한 타원체입니다. 그래서 R을 평균 반지름으로 둔 구면 공식은 실제 타원체 위의 거리와 최대 0.5% 정도 어긋나고, 더 정밀하게는 1975년 측지학자 타데우스 빈센티가 만든 빈센티 공식⁠(Vincenty's formulae)⁠ 같은 타원체 공식을 씁니다.
• 대원 거리는 세 규칙을 지키는 거리 함수⁠(metric)⁠이니, 가장 가까운 공항을 찾는 최근접 이웃⁠(nearest neighbor)⁠ 검색이나 구면 위의 보로노이 다이어그램⁠(Voronoi diagram)⁠에도 그대로 쓰입니다.