打开航旅App看北京飞洛杉矶的航班航线不是一条穿过太平洋的直线而是向北绕一圈经过俄罗斯远东、白令海最后再沿北美西海岸南下。第一次看到的人多半以为飞机在绕远其实这才是真正的近路。地球是圆的地图是平的地图上的直线在球面上往往是一段弯曲的劣弧真正的最短路径叫大圆航线学术上叫测地线。这篇文章想聊的就是把这条最短路径算出来这件事。无论你是写地图服务、做外卖或打车系统的LBS模块还是对GPS轨迹做里程统计都会碰到同一个核心问题已知A点和B点的经纬度怎么拿到地球表面两点间最短距离。适合的人群也很直接刚入门GIS或地图开发的工程师、搞不清Haversine和Vincenty区别的初学者以及想把自己的距离计算从“百度一个公式改改”升级成“知道为什么这么算、怎么避坑”的人。1. 为什么地图上的直线不是最短路径1.1 大圆与小圆球面“直线”的几何本质先做一次思想实验。拿一个地球仪找一根细线把线的一端按在北京另一端按在洛杉矶然后拉紧。线不会平铺在这两个城市所在的同一纬线上而是会微微“拱”向高纬度方向。因为球面上两点之间根本没有欧氏几何里的直线它的“直线”角色由大圆劣弧来扮演。所谓大圆是指过球心的平面与球面相交得到圆。赤道是典型的大圆每一条经线也是大圆的一部分。而普通纬线比如北纬40度那条线它的圆心并不在地球中心所以是小圆。小圆上的弧如果当作最短路径通常都会绕远。可以用切西瓜来帮助记忆一个平面从西瓜正中间劈下去西瓜皮上切出的那道圆环就是大圆如果你斜着切一面切出来的圆环偏在一边那就是小圆。球面上任意两点只要不是对跖点正好在球的相对两侧这两点和球心就能确定唯一一个平面这个平面与球面相交得到的圆就是过这两点的大圆。两点在大圆上把圆周分成两段弧较短的那段就是最短路径。实际生活里最直观的例子就是航班。从北京到旧金山机票上显示的航线基本会先朝东北方向靠近北极圈再向南。很多人以为这是为了避开政治空域或军事禁飞区其实核心原因是数学这条“弯曲”的路线才是真正的大圆劣弧距离更短。走常规纬线反而要多飞几百甚至上千公里。1.2 投影带来的误导墨卡托与航线错觉问题来了既然球面上最短路径是大圆劣弧为什么在普通地图上画出来它反而是一条弧线这就要说到地图投影了。地图是把三维球面强行压到二维平面的结果这个过程必然产生变形。航海常用的墨卡托投影有一个特点把纬线拉成等间距水平线同时在高纬度区域把东西方向也成比例放大于是格陵兰岛看起来和非洲一样大。在这个投影里球面上的大圆航线会被画成一条弧线而地图上看起来笔直的“直线”反而不是球面上的最短路径。对做开发的人来说这里藏着一个容易忽略的问题。如果在Web地图上直接拿两个经纬度点对应的平面像素坐标计算距离结果会偏离真实值纬度越高偏差越离谱。这不是公式不够好而是从一开始就在错误的坐标系里算距离。理解大圆距离之前先要把“投影画面里的直线长度”和“球面上的最短距离”彻底分开。2. 三种主流计算方法原理、公式与精度2.1 球面余弦定理直观但不建议短距离使用最早接触球面距离时最常见的就是球面余弦定理。给定两点经纬度φ1, φ2 是纬度λ1, λ2 是经度Δλ λ2 − λ1球面余弦定理把两点的球面夹角算出来cos c sinφ1·sinφ2 cosφ1·cosφ2·cosΔλ然后用地球半径 R 乘以这个夹角得到距离d R · arccos(c)公式本身非常直观就是把平面余弦定理推广到球面。但有一个重大缺陷当两点很近时cos c 的值非常接近1arccos(1) 的导数趋于无穷数值上极其敏感。浮点数稍微有点舍入误差结果就可能在几米到几十米之间乱跳极端情况下arccos的参数可能因为舍入变成1.0000000001直接返回NaN。所以我的建议是球面余弦定理适合用来理解不适合用来写生产代码尤其不要用在门店签到、短距离匹配这类经常出现“两点只隔几米”的场景。2.2 Haversine公式稳定、经典开发首选Haversine公式是导航界的老牌方案也是目前写距离计算最常用的公式。公式分两步a sin²(Δφ/2) cosφ1 · cosφ2 · sin²(Δλ/2)c 2 · atan2(√a, √(1−a))最终距离d R · c其中 Δφ φ2 − φ1Δλ λ2 − λ1。这里的 φ、λ 都需要先转成弧度。这个公式和球面余弦定理在数学上等价区别在于它绕开了 arccos。计算过程里用的是正弦函数当角度很小时正弦的变化率平缓不会出现“接近1时突然放大误差”的情况。所以无论两点距离是5米还是5000公里Haversine都能保持相对稳定的精度。还有一个细节很多人不注意完整公式里最后用 atan2而不是直接 asin(√a)。如果只用 asin当 a 因为浮点舍入变成1.0000000001你又会得到NaN。而 atan2(y, x) 天然对 x 为0或负值不敏感配合 clamp 就能把边界情况全部兜住。2.3 Vincenty公式椭球面上的高精度解法Haversine把地球当成标准球体但实际上地球是一个赤道半径略大于极半径的椭球体长半轴约6378.137公里短半轴约6356.752公里。对于跨越大洋的航线、精密测量、测绘级应用这21公里的半径差会让球面模型产生约0.3%到0.5%的误差换算成长距离可能是几十公里级别。Vincenty公式就是用来解决这个问题的。它把地球当成旋转椭球体通过迭代逼近求解两点间的测地线长度。核心思想是先把两点按椭球参数投影到辅助球面上用迭代法反复修正经度差直到误差小于阈值再用修正后的参数计算最终距离这里的迭代步骤我在项目里写过几乎所有语言都有现成实现直接贴完整代码反而容易让读者迷失在变量名里。关键是了解它的适用边界精度可以做到毫米级但对近似对跖点两点几乎在球体两端时不收敛可能出现迭代几十次还跳不出来。常规地图业务完全用不到这个精度只有在天文大地测量、海底光缆路径规划这些领域才有必要上。2.4 算法对比精度、复杂度与选型建议用一个表格总结三种方法方便日常选型方法地球模型典型误差适用距离代码复杂度球面余弦定理球体0.3% ~ 0.5%中长距离不适用近距离低Haversine公式球体0.3% ~ 0.5%任意距离近距离更稳低Vincenty公式椭球体毫米级高精度测量对跖点不收敛中高我在实际项目里的选型标准很简单如果只是做LBS、订单距离、运动轨迹里程统计Haversine加上一个合理的地球半径就够了如果做专业测绘或需要跟高精度设备数据比对再上Vincenty。盲目追求毫米级精度大多数业务场景不但用不上还会引入更多迭代异常分支需要处理。3. 从公式到实战完整落地流程与代码3.1 动手前先统一参数经纬度顺序、弧度与地球半径很多人写距离代码第一个坑不在公式而在参数。真实项目里经纬度来源五花八门有手机GPS、有地图SDK反查、有第三方接口返回。我见过好几次线上问题最后发现不是公式错是某个上游把 lat/lon 传反了。先说约定经纬度顺序。GeoJSON的规范是 [经度, 纬度]也就是 lon 在前但不少地图SDK和数据库习惯用 lat, lon。项目里必须统一成一种并且在接口文档里写明。我个人的习惯是内部统一用(lon, lat)因为和GeoJSON、PostGIS的习惯一致。角度与弧度。三角函数的 sin、cos、atan2 全部要求弧度但输入坐标通常是十进制度数。这个换算漏掉结果会离谱到完全不能用。地球半径。Haversine用的R常见取值有6371公里平均半径、6378.137公里赤道半径。两者对结果的影响在千分之几以内用哪个都行但要保证全项目统一最好做成配置项。还有一个容易被忽视的细节经纬度本身可以是小数也可能是度分秒(DMS)格式。比如“116°2345.6”是一分很多地图接口返回的是十进制度数进入公式前必须转统一。这个我在接第三方数据时踩过一次返回来的字段在某个边界值上突然变成度分秒字符串直接导致一整批距离计算失败。3.2 Haversine函数的直接实现Python/JS直接上一份可以用的Python实现import math def haversine_m(lon1, lat1, lon2, lat2): # 地球平均半径单位米 R 6371000.0 phi1 math.radians(lat1) phi2 math.radians(lat2) dphi math.radians(lat2 - lat1) dlambda math.radians(lon2 - lon1) a math.sin(dphi / 2) ** 2 math.cos(phi1) * math.cos(phi2) * math.sin(dlambda / 2) ** 2 # 防止浮点误差导致 a 超出 [0, 1] a min(1.0, max(0.0, a)) c 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a)) return R * cJavaScript版本几乎一样function haversineKm(lon1, lat1, lon2, lat2) { const R 6371.0; const rad x x * Math.PI / 180; const phi1 rad(lat1); const phi2 rad(lat2); const dphi rad(lat2 - lat1); const dlambda rad(lon2 - lon1); let a Math.sin(dphi / 2) ** 2 Math.cos(phi1) * Math.cos(phi2) * Math.sin(dlambda / 2) ** 2; a Math.min(1, Math.max(0, a)); const c 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a)); return R * c; }如果你看网上的老代码有人会用 a Math.asin(Math.sqrt(a)) 然后 d R * a其实也成立只是 asin 在输入接近1时没有 atan2 稳。我写代码时习惯把 clamp 那行单独注释出来因为后来维护的人如果看不懂可能会觉得多余而删掉删掉之后线上就会不定期报NaN。3.3 批量计算从千万次循环里挤出性能业务里经常要算的不是“两个点”而是“一个点对 N 个点”比如找出某个定位点附近500米内所有门店或者计算一条轨迹的总里程。如果每次都用 Python 的 for 循环调 Haversine几十万条数据就慢到用户能感知了。这时候建议用 numpy 做向量化计算。核心思路是把所有点的经纬度放进数组一次性做三角函数运算import numpy as np def haversine_np(lon1, lat1, lon2, lat2): lon1, lat1, lon2, lat2 map(np.radians, [lon1, lat1, lon2, lat2]) dlon lon2 - lon1 dlat lat2 - lat1 a np.sin(dlat / 2.0) ** 2 np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2.0) ** 2 a np.minimum(1.0, np.maximum(0.0, a)) c 2.0 * np.arctan2(np.sqrt(a), np.sqrt(1.0 - a)) return 6371.0 * c调用时 lon1、lat1 是标量或数组lon2、lat2 是数组一次就能返回所有距离。实测下来50万次距离计算在普通笔记本上大概几十毫秒比纯Python循环快两个数量级。如果数据量继续放大到亿级就不该在应用层硬算。通常先按 geohash 或者四叉树索引做粗筛把候选集合缩小到几百个点再做精确的Haversine重排。这一步不是优化是必须采取的设计。3.4 “附近的人”场景先粗筛再精算的经典套路举一个非常常见的业务用户打开App我需要返回他身边5公里内所有门店。如果直接对全量门店算Haversine效率太低。常规做法是先算一个经纬度范围把数据库查询限制在这个框内再用Haversine或者直接用SQL里的距离函数精排。粗筛的近似公式很简单纬度方向1个纬度大约对应110.574公里经度方向1个经度大约对应111.320 × cos(latitude) 公里假设中心点纬度是 lat要查 distance_km 范围内的点查询范围可以这样估算import math def rough_bounds(lat, lon, distance_km): lat_min lat - distance_km / 110.574 lat_max lat distance_km / 110.574 lon_min lon - distance_km / (111.320 * math.cos(math.radians(lat))) lon_max lon distance_km / (111.320 * math.cos(math.radians(lat))) return lat_min, lat_max, lon_min, lon_max要注意用 cos(latitude) 修正经度在低纬地区比较准确到高纬度地区误差会变大。因为经线每度长度是固定的约111公里纬线的实际长度则是随纬度变化的而“111.320 × cos”本身是球面近似90度附近会趋近0。所以这个粗筛方法只适合小范围几十公里内的初步过滤不能直接拿结果当最终距离用。当粗筛已经能排除90%以上无关数据之后再对剩下的点调用Haversine精确计算并排序性能就很从容了。4. 常见问题与排坑实录4.1 距离结果和地图App对不上先查坐标系这是我在社区里遇到求助最多的一类问题代码算出来的距离比地图App显示的远了一百多米或者定位点明明在店门口App却显示在马路对面。大多数情况下不是算法问题而是坐标系问题。GPS芯片直接输出的是WGS84坐标而国内主流地图App出于地理信息保密要求会把坐标加密成GCJ-02火星坐标。你在一个坐标系里拿到的经纬度放到另一个坐标系的地图上已经偏移了几十到几百米。在这个基础上算距离当然对不上。正规的解法是不要在业务层自己写坐标转换算法而是使用地图SDK提供的坐标转换接口或者服务端调用地图服务商的坐标转换API把WGS84转成GCJ-02之后再入库和计算。这样至少把坐标系拉齐了距离才能一致。另外要注意同一款Android手机在室外开阔地和城市高架桥下定位精度能差将近一个数量级。距离计算做得再准输入坐标本身有几十米抖动结果也必然跟着抖。这时候问题已经不是数学而是定位源质量。可以对定位点做简单的速度滤波或者轨迹平滑但那是另一个话题。4.2 短距离计算结果变成NaN或负数浮点越界短距离下最容易出现的两个bug一是用球面余弦定理arccos的参数因为浮点舍入略大于1直接报错二是Haversine里的 a 值因为负零或者其他舍入问题跑到 [0,1] 区间之外。出现这个问题的典型场景是用户原地不动两次定位坐标完全一样或者只差0.000001度。这时候 Δφ、Δλ 都接近0sin 算出来的值非常小再经过加法、乘法浮点误差就会放大。对策就是我在代码里写的 clampa min(1.0, max(0.0, a))不管 a 因为什么原因越界先把它拉回有效区间再用 atan2 计算。这个处理成本几乎为零但能避免线上出现“偶尔一个订单距离异常”的脏数据。另一个建议是函数内部统一用双精度浮点。JavaScript的 Number 就是双精度Python也是默认双精度没问题。但如果某些环境里用了单精度的 Float32短距离下精度损失会明显变大最好不要用。4.3 高纬度区域的距离偏差Web墨卡托投影陷阱做Web GIS的人容易犯一个错在地图上拿到两个marker的像素坐标然后按比例尺换算成实际距离。在小范围、低纬度地区误差还能接受但一旦跑到高纬度误差会快速放大。Web墨卡托投影为了让地图瓦片无缝衔接把椭圆体世界按球体公式投影。这种投影在纬度60度附近的长度拉伸大约是2倍也就是说在莫斯科或者北欧地区直接从平面图上量距离得到的结果可能比真实距离大上一倍或更多。如果你在哈尔滨北纬45度左右做精确测距平面换算的误差也已经有肉眼可见的偏差。正确做法是所有距离计算都回到经纬度上用球面公式或椭球公式投影平面只负责“显示”不负责“测量”。在Leaflet或Mapbox上画一个半径为500米的圆如果想要它代表真实距离也不应该直接换算成固定像素半径而是按照目标纬度重新计算每一帧投影后的点再连线成圆。否则在不同纬度看到的圆大小比例会失真。4.4 经纬度精度等级与测量可靠性“距离算不准”还有一个常被忽视的来源经纬度本身的分辨率。经纬度小数位数和实际分辨率的关系大致如下小数位数纬度方向分辨率经度方向分辨率低纬3位约111米约111米4位约11米约11米5位约1.1米约1.1米6位约0.11米约0.11米手机定位在开阔地带通常只能做到5到10米的精度城市峡谷和多径环境下误差甚至会到50米。所以当你发现Haversine计算结果和地图App相差几十米时不必过度怀疑算法先检查两端坐标分别是从哪里取的、精度等级是多少。我个人的做法是在业务侧把“计算距离”和“定位精度”分开对待。只要是GPS来源距离结果统一保留到整数米界面展示用“米/公里”不做花哨的毫米位输出。这既避免给用户虚假的精确感也让排查问题时更容易接受误差范围。5. 个人实践中的几个细节补充5.1 坐标系先于算法在项目里调距离计算我踩过最深的坑就是坐标系不统一。那次是接第三方门店数据对方给的是WGS84经纬度App端拿到的却是GCJ-02地图坐标两边混在一起算“附近门店”结果用户明明站在店门口匹配到的却是隔壁两条街的门店。排查了两天才定位到问题最后是统一走地图SDK的坐标转换接口才解决。所以我的经验是开始写距离函数之前先把所有数据源的坐标系在入口处统一。所有纬度、经度字段都要确认是WGS84、GCJ-02还是其他平台坐标并且落到数据表里单独建一列存坐标系标识。坐标系没统一再好的公式也白搭。5.2 用一个已知距离做基准测试距离计算代码写完不要只看样例数据顺眼就上线。我一般会准备一组“基准用例”比如上海到北京的实际大圆距离约1067公里在单元测试里跑Haversine断言结果落在正负1%的范围内。另外再加一个“同点距离为0”的用例专门用来验证NaN和越界问题。有这几个用例兜底后续换地球半径、换公式甚至换坐标系转换逻辑回归测试都能直接暴露问题。别看距离计算是个很基础的函数一旦线上出了问题排查成本往往比写一万行业务逻辑还高。最后分享一个从实战里得到的意见距离公式怎么选优先级永远排在“坐标统一”和“误差定位”之后。默认用Haversine封装成公共函数把所有参数、单位、半径、坐标系提前订死比在业务里到处散落Math.sin和Math.cos要省心得多。