【Matlab】根據經緯度計算兩點間的球面距離

來自專欄數學建模與數學實驗14 人贊了文章

做建模或者研究空間數據，可能會遇到「根據經緯度計算兩點間的球面距離」的問題，網上的資料很多，都是各種公式推導，但是一旦按公式編程計算，很可能得不到正確的距離。根本原因是在「角度-弧度的轉化」與「軟體語法」結合上出錯了。

為此，我梳理了兩種常用的計算方法，並用 Matlab 編程實現，保證能計算出正確的距離。

方法一：Great-Circle距離（基於球面餘弦公式）

計算公式為：

其中， 為地球半徑， 分別表示 點的經度、緯度； 分別表示 點的經度、緯度。

公式的具體推導過程略（可參考文獻[1][2]）

Matlab實現：

先定義一個角度轉弧度的函數：

function rad = D2R(theta) rad = theta*pi/180;

再定義計算Great-Circle距離的函數：

function d = SphereDist(x,y,R)%根據兩點的經緯度計算大圓距離(基於球面餘弦公式)%x為A點[經度, 緯度], y為B點[經度, 緯度]if nargin < 3 R = 6378.137;endx = D2R(x);y = D2R(y);DeltaS = acos(cos(x(2))*cos(y(2))*cos(x(1)-y(1))+sin(x(2))*sin(y(2)));d = R*DeltaS;

註：R 設為了默認參數，這裡要吐槽一下 Matlab 在這方面的啰嗦，其他語言直接 SphereDist(x,y,R=6378.137) 就行了。

方法二：基於Haversine公式

前面的球面餘弦公式中有 項，當兩點間距離很短時（比如相距幾百米的兩點），餘弦函數會得出 的結果，會導致較大的舍入誤差。Haversine 方法進行了某種變換消除了 項，能夠避免該問題。

實際上，以現在電腦軟體的一般計算精度，這已經不是問題。用計算器算可能有這個問題。

Haversine公式：

其中，

這裡， 為地球半徑； 表示兩點的緯度； 表示兩點經度的差值。

公式的具體推導過程略（可參考文獻[2]）

Matlab實現：

function d = SphereDist2(x,y,R)%根據兩點的經緯度計算大圓距離(基於Haversine公式)%x為A點[經度, 緯度], y為B點[經度, 緯度]if nargin < 3 R = 6378.137;endx = D2R(x);y = D2R(y);h = HaverSin(abs(x(2)-y(2)))+cos(x(2))*cos(y(2))*HaverSin(abs(x(1)-y(1)));d = 2 * R * asin(sqrt(h));function h = HaverSin(theta) h=sin(theta/2)^2;endend

註：嵌套了一個 HaverSin 函數。

測試：

%計算北京到天津的球面距離SphereDist([116.4,39.9], [117.2, 39.1]) SphereDist2([116.4,39.9], [117.2, 39.1])

運行結果：

參考文獻：

[1] 根據兩點的經緯度求方位角和距離，等

https://www.cnblogs.com/lihaibo-Leao/p/5142154.html

[2] 地理空間距離計算及優化（根據兩個點經緯度計算距離），https://blog.csdn.net/u011001084/article/details/52980834

原創作品，轉載請註明。