利用MATLAB实现空间两点距离的计算以及优化

首先在这里声明，本次使用的地球模型采用的是球形。参考的文章是美团的距离计算，说的很清晰，感兴趣的朋友可以参考。我自己使用了MATLAB实现了其中的算法。仅供参考，互相学习。

大致地说，这里实现的方法就是通过将经纬坐标转换为三维的以球心为原点的立体坐标计算地球表面两点之间的距离长度。

地理空间距离计算优化—美团技术团队

function [x y z] = LL2era(lon,lat)
%函数实现将经纬度坐标转换为以地心为原点的三维立体坐标
%lon：输入精度
%lat：输入纬度
%模型采用地球为球形
lat = lat*pi/180;
lon = lon*pi/180;
r = 6367000; %地球的半径
x = r*cos(lat)*cos(lon);
y = r*cos(lat)*sin(lon);
z = r*sin(lat);

function [distance1,distance2] = distance_compute(lon1,lat1,lon2,lat2)
%根据给出的两点之间的经纬度计算地球上两点之间的距离
%计算两点之间的距离
%对于地球球体来说，计算地球表面两点之间的距离就是求这两点相对地心形成的夹角对应的弧度
%那么为了求得这个弧度的长度，我们需要得到夹角的大小，再使用余弦公式得到弧度对应的边
%最后根据弧度计算公式获得距离
r = 6367000; %地球的半径
[x1 y1 z1] = LL2era(lon1,lat1);
[x2 y2 z2] = LL2era(lon2,lat2);
distance_direcr2 = (x1 - x2)^2 + (y1 - y2)^2 + (z1 - z2)^2; %弧长对应的弦长
costheta = (distance_direcr2 - 2*r^2)/(-2*r^2); %夹角的余弦值
distance1 = r*acos(costheta); %半径乘以角度等于弧长

%改进1，为了加快计算的速度，将arccostheta用sintheta代替，再用sqrt(1-costheta^2)变换
distance2 = r*sqrt((1  - costheta^2));

%同时，下面将给出更加一般化的简化函数
function distance = abstract_compute(lon1,lat1,lon2,lat2)
%以一种简化的方式计算距离
%由于计算的距离比较小，所以可将经纬线认为是互相垂直的
%于是就可以根据垂直、直线之间的关系来求解
r = 6367000; %地球的半径
dx = (lon1 - lon2)*pi/180; %经度差并转换为弧度
dy = (lat1 - lat2)*pi/180; %纬度差并转换为弧度
average = ((lat1 + lat2)/2)*pi/180; %平均纬度
SN = r*dy; %南北方向上的距离
SE = r*dx*cos(average); %东西方向上的距离
distance = sqrt(SN^2 + SE^2);

输入输出示例：


最后一步的简化就是将余弦函数用泰勒公式展开表示，只要展开到三次幂就可以得到误差比较小的结果了！