平时做轨迹分析、配送调度或者地图相关的开发最容易被“经纬度直接套欧氏距离”这个操作坑到。拿到两对经纬度坐标想也不想就sqrt((lat1-lat2)^2(lon1-lon2)^2)再乘个 111000 当成米数结果短距离时好像还行一旦跨纬度、跨经度稍微大一点误差能大到让人怀疑数据是不是错了。原因很简单经纬度是球面坐标不是平面坐标度本身也不是线性的长度单位直接当平面坐标去算欧氏距离等于默认地球是一张可以随意展开的纸。这篇博客就把这件事彻底讲明白为什么不能直接算、正确的转化思路是什么、在 MATLAB 里怎么落地成代码以及不同场景下投影误差到底有多大、什么时候该换方法。内容偏工程实践适合做GIS数据处理、写定位算法或者刚接触轨迹分析的同学参考可以直接照抄代码也顺便把背后的原理补上。1. 直接拿经纬度差算“距离”误差到底从哪里来1.1 经度1度在不同纬度的地面长度完全不同先看一个最简单的事实。纬度方向还好说1度纬度对应的地面距离在大部分地区大约是 110 到 111 公里这个数值相对稳定。但经度方向的 1 度就“看人下菜碟”了赤道上 1 度经线大约是 111.3 公里到了北纬 60 度同样的 1 度经度差只有大约 55.7 公里而在南北极附近直接趋近于 0。用 MATLAB 跑一小段就特别直观lats [0 30 60 80]; for i 1:length(lats) dist1degLon 111320 * cosd(lats(i)); fprintf(纬度 %.0f°: 1度经度对应的地面距离 %.1f m\n, ... lats(i), dist1degLon); end输出结果纬度 0°: 1度经度对应的地面距离 111320.0 m 纬度 30°: 1度经度对应的地面距离 96413.8 m 纬度 60°: 1度经度对应的地面距离 55660.0 m 纬度 80°: 1度经度对应的地面距离 19333.4 m同一个 1 度经度差在赤道附近是 111 公里在北纬 80 度却只剩下 19 公里。如果你拿着平面公式统一乘 111 公里那在北纬 60 度以上算出来的东西方向距离会被成倍夸大。这就是“直接把经纬度当成平面坐标”的第一个致命伤度不是米而且度在不同位置对应的米数还不一样。1.2 球面曲率带来的不可忽略误差第二个问题在于地球表面是一个球面或者说是一个椭球面而任何平面展开都会引入变形。即便你已经修正了经度方向上的长度系数只要两点跨度够大球面的弯曲本身也会让“短程线”和“平面上画的直线”产生差异。想象地球是一个半径为 R 的球球面上两点之间的最短路径是大圆航线对应弧长是R * Δσ其中Δσ是两点对球心的夹角。如果你把它们投影到切平面上再算欧氏距离得到的是弦长近似为2R * sin(Δσ/2)。当Δσ很小的时候弧长和弦长非常接近但当Δσ变大比如两点之间隔着 3000 公里时弧长和弦长之间的差异能到 1% 以上。这个误差在路径规划、油耗估算这种对距离敏感的任务里是不能忽略的。1.3 一个量化演示平面公式变形到什么程度我看过很多人用一句话经验公式做快速估算纬度差乘 111 公里经度差乘 111 公里再乘cos(平均纬度)然后勾股定理算出距离。这个公式在小范围内是好用的本质上它就是本地切平面投影的一种近似但它有一个隐藏假设两点所在的纬度范围足够窄经度差也足够小。举个例子假设 A 点在某华东城市附近北纬 31.2 度、东经 121.5 度B 点在华北某市附近北纬 39.9 度、东经 116.4 度。两者纬度差约 8.7 度、经度差约 5.1 度。如果直接按 111 公里/度做平面估算dLat 39.9042 - 31.2304; dLon 116.4074 - 121.4737; dEst sqrt(dLat^2 dLon^2) * 111; % 单位: km算出来大约是 1115 公里。但用后面要讲的 Haversine 大圆距离真实地面距离大约是 1087 公里左右误差接近 2.6%。如果是跨国、跨大洲的数据这种误差会膨胀到几十甚至上百公里。所以结论很明确平面公式只能用于局部小区间它不是“近似正确”而是“在小范围恰好近似正确”。2. 两条路线怎么选要距离标量还是要平面坐标理解完误差来源接下来就要选技术路线。很多人不知道的一点是“算距离”和“得到欧氏距离”其实是两码事取决于你的下游任务到底是什么。2.1 路线A只要距离标量用大圆距离就好如果你的目标就是单纯得到两个地理位置之间的地面距离比如计算骑行里程、判断两点是否在 5 公里范围内那根本不需要把经纬度转化成平面坐标。直接用球面大圆距离公式也就是 Haversine 公式MATLAB 几行代码就能搞定。Haversine 公式基于球面三角学通过两点经纬度直接算出球面弧长。它在中短距离几米到几千公里范围内精度足够并且不依赖任何投影参数。需要注意它假设地球是半径 6371 公里左右的球体而真实地球是一个椭球所以在高精度测量场景下会有 0.1% 到 0.5% 的系统误差。不过对于绝大多数工程任务这个误差远小于其他误差来源。2.2 路线B要平面坐标才需要做投影转换但另一类场景是必须要“真正”的平面坐标的。我举几个典型的做 DBSCAN 或 KMeans 聚类算法假设特征空间是欧氏空间直接把经纬度塞进去聚类出来的“距离”没有任何物理意义。做栅格化统计把轨迹点分配到固定大小的网格里需要明确的以米为单位的横纵坐标。做最近邻检索需要把点集放入空间索引比如 R 树平面坐标更利于建立包围盒和距离下界。做轨迹相似度比对比如 DTW 或 Fréchet 距离中间大量依赖位移向量的加减必须是统一的米制平面坐标。这些场景讲究的是“坐标”本身而不只是一个距离标量。所以你需要先把经纬度投影到一个二维平面上得到以米为单位的 (x, y)然后再计算欧氏距离。这个过程才是标题里说的“从经纬度坐标到欧式距离的转化”。2.3 常见投影方案对比把球面展开成平面世界上的手段五花八门但工程里常用的大概就下面几类投影方式变形特征实现复杂度适用场景等距圆柱投影高纬度面积和距离变形大极简单快速可视化、演示高斯-克吕格投影中央经线附近无变形越远越大需要带号与椭球参数国内测绘、地形图UTM 投影中央经线缩放 0.9996带宽 6 度全球统一需计算带号跨国、大范围轨迹本地切平面ENU参考点附近几乎无变形远离参考点误差增大最简单城市级、区域级工程自实现时最推荐的是本地切平面投影因为它不需要复杂的椭球迭代公式只凭“以参考点为中心把球面近似成切平面”的思路就能得到亚米级到米级的精度特别适合覆盖范围在几百公里内的业务场景。大范围的数据则建议用 UTM 或者直接回退到大圆距离。3. MATLAB 三件套大圆距离、本地切平面、UTM 投影的实现这一部分给完整可复制的代码。为了避免教程打架我的建议是三个函数都写下来放进一个工具包里按场景调用。3.1 Haversine 距离函数一个文件解决“距离基准”function d haversineDist(lat1, lon1, lat2, lon2, unit) % 大圆距离计算默认返回米 % 输入均为十进制度lat/lon 可支持标量或等长向量 % unit 可选 m米或 km千米 if nargin 5 unit m; end R 6371000; % 地球平均半径单位米 lat1 deg2rad(lat1); lon1 deg2rad(lon1); lat2 deg2rad(lat2); lon2 deg2rad(lon2); dLat lat2 - lat1; dLon lon2 - lon1; a sin(dLat/2).^2 cos(lat1) .* cos(lat2) .* sin(dLon/2).^2; c 2 * atan2(sqrt(a), sqrt(1 - a)); d R * c; if strcmpi(unit, km) d d / 1000; end end注意这里用atan2(sqrt(a), sqrt(1-a))而不是直接用asin(sqrt(a))主要是数值稳定性考虑避免a由于浮点舍入出现比 1 大一点点、进而导致asin报 NaN 的情况。这个函数可以当成其他方法的“参照基准”当你想知道投影法的误差有多大就把投影距离和它对比。3.2 本地切平面ENU投影多数工程场景的最佳自选本地切平面投影的思路很朴素在参考点(refLat, refLon)处做一个与地球表面相切的平面然后把目标点投影到这个平面上去。得到的东北坐标东向、北向就是相对参考点的平面位移。function [xEast, yNorth] enuFromLatLon(lat, lon, refLat, refLon) % 将经纬度转为以参考点为中心的本地东北坐标单位米 % 输入均为十进制度lat/lon 可支持向量 R 6371000; % 地球平均半径单位米 refLatRad deg2rad(refLat); refLonRad deg2rad(refLon); latRad deg2rad(lat); lonRad deg2rad(lon); dLat latRad - refLatRad; dLon lonRad - refLonRad; yNorth R * dLat; xEast R * cos(refLatRad) * dLon; end这个近似在局部是几何直观的它认为目标点的纬线平行于参考点纬线且球面曲率小到可以忽略。于是北向位移就是子午线方向上的弧长R * dLat东向位移近似为参考点所在纬线圈上的弧长R * cos(refLat) * dLon。有了东北坐标任意两点之间的欧氏距离直接就是dist sqrt((x1 - x2)^2 (y1 - y2)^2);也可以写成一个更友好的封装函数function d euclideanDistFromLatLon(lat1, lon1, lat2, lon2, refLat, refLon) % 基于本地切平面投影的两点欧氏距离单位米 % 如果未指定参考点默认取两点中点的纬度/经度 if nargin 6 refLat (lat1 lat2) / 2; refLon (lon1 lon2) / 2; end [x1, y1] enuFromLatLon(lat1, lon1, refLat, refLon); [x2, y2] enuFromLatLon(lat2, lon2, refLat, refLon); d sqrt((x1 - x2).^2 (y1 - y2).^2); end注意点同一个数据集里的所有点必须用同一个参考点投影。如果每个点都拿自己做参考点坐标基准不一致减法没有任何意义。参考点的选择原则很简单选数据覆盖范围的几何中心这样整个区域内的投影误差是均匀分布的。3.3 使用地图工具箱进行全局投影大跨度轨迹的备选如果数据的跨度到了几百公里甚至上千公里本地切平面的误差会开始明显累积。这时候有两种选择一是用大圆距离直接算二是做全局投影比如 UTM。MATLAB 从 R2019b 开始提供了projcrs配合projfwd的用法对于 UTM 投影来说非常方便% 假设数据大致在东经 114-120 度、北纬 29-40 度对应 UTM 50N % EPSG:32650 表示 WGS84 坐标系下的 UTM 50N 带 utmProj projcrs(32650, Authority, EPSG); [x, y] projfwd(utmProj, lat, lon);得到x, y已经是以米为单位的平面坐标直接算欧氏距离即可。326XX是北半球 UTM 带号327XX是南半球。国内通常用 49 到 52 带具体带号要根据经度换算带号约为floor((lon 180) / 6) 1。如果没有安装 Mapping Toolbox自实现 UTM 投影会涉及椭球参数、墨卡托级数展开这类比较繁琐的数学不建议手写。更务实的做法是如果跨带或者大跨度直接用 3.1 节的大圆距离或者把数据分块后分别用本地切平面。3.4 批量实战一组 GPS 点到参考点的距离计算假设有一组经纬度点想计算它们到某个中心点的距离并对比投影距离与大圆距离的差异可以直接这样写pts [ 31.2304, 121.4737; 39.9042, 116.4074; 30.5728, 114.3162; 23.1291, 113.2644 ]; refLat 34; refLon 118; % 本地切平面投影 [x, y] enuFromLatLon(pts(:,1), pts(:,2), refLat, refLon); dProj sqrt(x.^2 y.^2); % 大圆距离基准 dTrue haversineDist(refLat, refLon, pts(:,1), pts(:,2), km); % 查看差异 disp(table(dProj/1000, dTrue, (dProj/1000 - dTrue)./dTrue*100, ... VariableNames, {投影距离km, 大圆距离km, 相对误差%}));在几百公里跨度下投影距离与大圆距离的差异通常在 0.1% 到 0.5% 之间表现是相当可用的。如果误差超预期大概率是参考点选得太偏或者数据跨度过大。4. 投影误差的边界多少公里内可以放心用拿到代码之后心里得有一杆秤我的数据跨度多大用本地切平面还能不能顶住这一章就用脚本把边界测出来免得以后拍脑袋。4.1 基于参考点投影的距离误差随跨度变化以参考点北纬 30 度、东经 110 度为例让目标点在东西南北各个方向上偏移 0 到 10 度逐点比较本地切平面投影距离与大圆距离。lat0 30; lon0 110; northSteps -10:0.5:10; eastSteps -10:0.5:10; err zeros(length(northSteps), length(eastSteps)); for i 1:length(northSteps) for j 1:length(eastSteps) lat lat0 northSteps(i); lon lon0 eastSteps(j); [x, y] enuFromLatLon(lat, lon, lat0, lon0); dProj sqrt(x^2 y^2); dTrue haversineDist(lat0, lon0, lat, lon); % 已返回米 err(i,j) (dProj - dTrue) / dTrue * 100; end end % 观察最大绝对误差范围 fprintf(最大相对误差(%%): %.3f\n, max(abs(err(:))));基于球面几何可以推导正东、正西方向偏移时的误差主要由弧长与弦长的差异决定正北、正南方向偏移时在球体模型下的投影距离近似等于弧长误差反而很小。最大的误差出现在同时包含东西和南北方向的大跨度位移因为此时投影平面上的直线相当于球面大圆的一条弦距离会系统性偏小。实测下来在 5 度跨度约 550 公里内最大相对误差通常在 1% 到 2% 以内在 2 度跨度约 220 公里内误差普遍小于 0.2%。4.2 高纬度与参考点选择对误差的影响本地切平面投影的东向系数是R * cos(refLat)如果你的目标点实际纬度与参考点纬度相差很大这个系数就会失真。比如参考点在北纬 30 度目标点在北纬 60 度目标点附近真实的 1 度经线长度只有约 55.7 公里但投影公式按cos(30°)算出来是约 96.6 公里东西方向会偏差接近 74%。这种极端情况很容易遇到全国范围的配送数据、跨省车辆轨迹如果硬要用同一个本地切平面误差会非常难看。所以我的经验判断是50 公里以内本地切平面误差通常小于 0.05%放心用。200 公里以内误差一般小于 0.3%绝大多数业务可以接受。500 公里以内误差可能到 1% 到 2%需要结合精度要求判断。1000 公里以上不要再用本地切平面改用大圆距离或 UTM 投影。4.3 什么时候必须放弃平面近似需要特别警惕的是两类场景。第一类是高纬度区域比如北纬 70 度以上的极区附近投影变形本身就很敏感处理这类数据时优先选择极地方位投影或者干脆用球面距离。第二类是跨国跨洋的长距离轨迹本地切平面已经完全失去意义此时投影距离与真实大圆距离的偏差可达数个百分点而且由于地球曲率投影坐标下的直线方向也会失真。如果既想要平面坐标来跑下游算法又想要大跨度距离精度实际工程里有个折中方案按区域分块每块单独建立本地切平面参考点。比如按 2 度 x 2 度的格网切块每个网格内的点用网格中心做参考点投影既能延续欧氏空间计算又能把误差控制在很小的范围内。5. 坐标转换里最容易翻车的五个细节即使掌握了上面的核心代码实际工程中还是很常见各种奇葩问题。下面这几个细节几乎每个做轨迹数据处理的人都踩过。5.1 坐标顺序lat/lon 还是 lon/latMATLAB 地图相关函数比如geoplot的习惯是纬度在前、经度在后也就是(lat, lon)。但很多数据源恰好相反比如 GeoJSON 的标准是[lon, lat]部分设备输出也是先把经度放在前面。一旦写反所有计算全部错乱而且很难通过肉眼发现。更隐蔽的是经度规范化问题。比如一个点在东经 179.9 度另一个点在西经 -179.9 度它们的实际经度差只有 0.2 度但如果直接做减法会得到 359.8 度算出来距离远得离谱。这种情况在跨太平洋数据里特别常见一定要先做规范化function dlon normalizedLon(dLonDeg) dlon mod(dLonDeg 180, 360) - 180; end5.2 度分秒与十进制度混用很多老式 GPS 设备输出的是度分格式比如31.2304可能代表 31 度 23.04 分而计算时必须用十进制度。直接把度分当成十进制度参与计算误差大概是 60 倍量级距离一下子偏出去十分离谱。转换公式是十进制度 度 分/60 秒/3600。MATLAB 里也有现成函数dms2degrees但关键是入库和清洗阶段就要统一不能在计算时才发现单位不一致。5.3 坐标系不统一数据源之间可能存在系统性偏移GPS 原始数据通常是 WGS84 坐标但国内很多地图服务返回的坐标是经过偏移变换的本地坐标系两者之间有几十米到几百米的系统性偏差。如果在同一个工程里混用了不同坐标系的数据去做距离筛选看起来每个点单独没问题但两两之间的距离会引入不可控误差边界过滤、聚合判断全都不可信。我处理这种问题时都会要求数据接入的第一步做基准校验拿同一个物理位置在不同数据源里的坐标做差确认坐标系是否统一。如果确实不一致必须用官方或公开的转换算法统一到同一个基准后再入库。注意转换过程也有精度损失定量评估后再决定是否接受。5.4 投影带号与参考点选错用 UTM 投影时如果数据跨了两个带强行使用同一个带号距离会异常偏大。比如东经 117 度和东经 123 度分别落在 UTM 50N 和 51N直接用 50N 投影两个点横坐标会被拉开几千公里显然不可用。处理方式要么分带投影要么回到大圆距离。本地切平面的参考点选择同样关键。参考点应该放在数据覆盖的中心而不是数据集边界点或者随便一个点。参考点放偏投影误差只在一侧膨胀距离计算会表现出明显的系统性偏差。用中心点时误差围绕真实值均匀分布而且最大值最小化。5.5 地球半径取值的玄学代码里我统一用 6371000 米这是国际大地测量与地球物理联合会推荐的平均半径。有的实现会用地球赤道半径 6378137 米有的会用 6356752 米极半径三者相差最多约 0.7%。对一般工程来说0.1% 到 0.7% 的差距几乎无感知但在高精度基线计算里就会成为问题。如果需要更高精度建议直接换用基于椭球模型的 Vincenty 公式而不是纠结平均半径到底取哪个数。MATLAB 里如果安装了 Mapping Toolboxdistance函数默认就用椭球模型在长距离上的表现比球面 Haversine 更可靠但短距离内两者差距微乎其微。最后说点我自己的习惯。我处理一批轨迹或者配送点位数据时第一步永远是看覆盖范围和用途如果所有点集中在几十公里范围内直接本地切平面投影参考点选数据范围中心然后放心大胆地用欧氏距离如果跨到省一级就考虑 UTM 或者大圆距离如果既有投影需求又要跨大范围就分块建局部坐标系绝不在全局硬切。另外我每次写完投影代码都会顺手做一次“基准校验”随机抽几对点用haversineDist和投影距离对比偏差超过预期就说明参考点或者带号有问题早发现早解决。这套流程帮我避开了很多坐标相关的坑希望能给同样被经纬度折磨的人一点参考。 
长河君 · 2026/10/15 17:07:43 · 阅读 1.2k 经纬度如何转成欧氏距离？MATLAB 投影与大圆距离实战
平时做轨迹分析、配送调度或者地图相关的开发最容易被“经纬度直接套欧氏距离”这个操作坑到。拿到两对经纬度坐标想也不想就sqrt((lat1-lat2)^2(lon1-lon2)^2)再乘个 111000 当成米数结果短距离时好像还行一旦跨纬度、跨经度稍微大一点误差能大到让人怀疑数据是不是错了。原因很简单经纬度是球面坐标不是平面坐标度本身也不是线性的长度单位直接当平面坐标去算欧氏距离等于默认地球是一张可以随意展开的纸。这篇博客就把这件事彻底讲明白为什么不能直接算、正确的转化思路是什么、在 MATLAB 里怎么落地成代码以及不同场景下投影误差到底有多大、什么时候该换方法。内容偏工程实践适合做GIS数据处理、写定位算法或者刚接触轨迹分析的同学参考可以直接照抄代码也顺便把背后的原理补上。1. 直接拿经纬度差算“距离”误差到底从哪里来1.1 经度1度在不同纬度的地面长度完全不同先看一个最简单的事实。纬度方向还好说1度纬度对应的地面距离在大部分地区大约是 110 到 111 公里这个数值相对稳定。但经度方向的 1 度就“看人下菜碟”了赤道上 1 度经线大约是 111.3 公里到了北纬 60 度同样的 1 度经度差只有大约 55.7 公里而在南北极附近直接趋近于 0。用 MATLAB 跑一小段就特别直观lats [0 30 60 80]; for i 1:length(lats) dist1degLon 111320 * cosd(lats(i)); fprintf(纬度 %.0f°: 1度经度对应的地面距离 %.1f m\n, ... lats(i), dist1degLon); end输出结果纬度 0°: 1度经度对应的地面距离 111320.0 m 纬度 30°: 1度经度对应的地面距离 96413.8 m 纬度 60°: 1度经度对应的地面距离 55660.0 m 纬度 80°: 1度经度对应的地面距离 19333.4 m同一个 1 度经度差在赤道附近是 111 公里在北纬 80 度却只剩下 19 公里。如果你拿着平面公式统一乘 111 公里那在北纬 60 度以上算出来的东西方向距离会被成倍夸大。这就是“直接把经纬度当成平面坐标”的第一个致命伤度不是米而且度在不同位置对应的米数还不一样。1.2 球面曲率带来的不可忽略误差第二个问题在于地球表面是一个球面或者说是一个椭球面而任何平面展开都会引入变形。即便你已经修正了经度方向上的长度系数只要两点跨度够大球面的弯曲本身也会让“短程线”和“平面上画的直线”产生差异。想象地球是一个半径为 R 的球球面上两点之间的最短路径是大圆航线对应弧长是R * Δσ其中Δσ是两点对球心的夹角。如果你把它们投影到切平面上再算欧氏距离得到的是弦长近似为2R * sin(Δσ/2)。当Δσ很小的时候弧长和弦长非常接近但当Δσ变大比如两点之间隔着 3000 公里时弧长和弦长之间的差异能到 1% 以上。这个误差在路径规划、油耗估算这种对距离敏感的任务里是不能忽略的。1.3 一个量化演示平面公式变形到什么程度我看过很多人用一句话经验公式做快速估算纬度差乘 111 公里经度差乘 111 公里再乘cos(平均纬度)然后勾股定理算出距离。这个公式在小范围内是好用的本质上它就是本地切平面投影的一种近似但它有一个隐藏假设两点所在的纬度范围足够窄经度差也足够小。举个例子假设 A 点在某华东城市附近北纬 31.2 度、东经 121.5 度B 点在华北某市附近北纬 39.9 度、东经 116.4 度。两者纬度差约 8.7 度、经度差约 5.1 度。如果直接按 111 公里/度做平面估算dLat 39.9042 - 31.2304; dLon 116.4074 - 121.4737; dEst sqrt(dLat^2 dLon^2) * 111; % 单位: km算出来大约是 1115 公里。但用后面要讲的 Haversine 大圆距离真实地面距离大约是 1087 公里左右误差接近 2.6%。如果是跨国、跨大洲的数据这种误差会膨胀到几十甚至上百公里。所以结论很明确平面公式只能用于局部小区间它不是“近似正确”而是“在小范围恰好近似正确”。2. 两条路线怎么选要距离标量还是要平面坐标理解完误差来源接下来就要选技术路线。很多人不知道的一点是“算距离”和“得到欧氏距离”其实是两码事取决于你的下游任务到底是什么。2.1 路线A只要距离标量用大圆距离就好如果你的目标就是单纯得到两个地理位置之间的地面距离比如计算骑行里程、判断两点是否在 5 公里范围内那根本不需要把经纬度转化成平面坐标。直接用球面大圆距离公式也就是 Haversine 公式MATLAB 几行代码就能搞定。Haversine 公式基于球面三角学通过两点经纬度直接算出球面弧长。它在中短距离几米到几千公里范围内精度足够并且不依赖任何投影参数。需要注意它假设地球是半径 6371 公里左右的球体而真实地球是一个椭球所以在高精度测量场景下会有 0.1% 到 0.5% 的系统误差。不过对于绝大多数工程任务这个误差远小于其他误差来源。2.2 路线B要平面坐标才需要做投影转换但另一类场景是必须要“真正”的平面坐标的。我举几个典型的做 DBSCAN 或 KMeans 聚类算法假设特征空间是欧氏空间直接把经纬度塞进去聚类出来的“距离”没有任何物理意义。做栅格化统计把轨迹点分配到固定大小的网格里需要明确的以米为单位的横纵坐标。做最近邻检索需要把点集放入空间索引比如 R 树平面坐标更利于建立包围盒和距离下界。做轨迹相似度比对比如 DTW 或 Fréchet 距离中间大量依赖位移向量的加减必须是统一的米制平面坐标。这些场景讲究的是“坐标”本身而不只是一个距离标量。所以你需要先把经纬度投影到一个二维平面上得到以米为单位的 (x, y)然后再计算欧氏距离。这个过程才是标题里说的“从经纬度坐标到欧式距离的转化”。2.3 常见投影方案对比把球面展开成平面世界上的手段五花八门但工程里常用的大概就下面几类投影方式变形特征实现复杂度适用场景等距圆柱投影高纬度面积和距离变形大极简单快速可视化、演示高斯-克吕格投影中央经线附近无变形越远越大需要带号与椭球参数国内测绘、地形图UTM 投影中央经线缩放 0.9996带宽 6 度全球统一需计算带号跨国、大范围轨迹本地切平面ENU参考点附近几乎无变形远离参考点误差增大最简单城市级、区域级工程自实现时最推荐的是本地切平面投影因为它不需要复杂的椭球迭代公式只凭“以参考点为中心把球面近似成切平面”的思路就能得到亚米级到米级的精度特别适合覆盖范围在几百公里内的业务场景。大范围的数据则建议用 UTM 或者直接回退到大圆距离。3. MATLAB 三件套大圆距离、本地切平面、UTM 投影的实现这一部分给完整可复制的代码。为了避免教程打架我的建议是三个函数都写下来放进一个工具包里按场景调用。3.1 Haversine 距离函数一个文件解决“距离基准”function d haversineDist(lat1, lon1, lat2, lon2, unit) % 大圆距离计算默认返回米 % 输入均为十进制度lat/lon 可支持标量或等长向量 % unit 可选 m米或 km千米 if nargin 5 unit m; end R 6371000; % 地球平均半径单位米 lat1 deg2rad(lat1); lon1 deg2rad(lon1); lat2 deg2rad(lat2); lon2 deg2rad(lon2); dLat lat2 - lat1; dLon lon2 - lon1; a sin(dLat/2).^2 cos(lat1) .* cos(lat2) .* sin(dLon/2).^2; c 2 * atan2(sqrt(a), sqrt(1 - a)); d R * c; if strcmpi(unit, km) d d / 1000; end end注意这里用atan2(sqrt(a), sqrt(1-a))而不是直接用asin(sqrt(a))主要是数值稳定性考虑避免a由于浮点舍入出现比 1 大一点点、进而导致asin报 NaN 的情况。这个函数可以当成其他方法的“参照基准”当你想知道投影法的误差有多大就把投影距离和它对比。3.2 本地切平面ENU投影多数工程场景的最佳自选本地切平面投影的思路很朴素在参考点(refLat, refLon)处做一个与地球表面相切的平面然后把目标点投影到这个平面上去。得到的东北坐标东向、北向就是相对参考点的平面位移。function [xEast, yNorth] enuFromLatLon(lat, lon, refLat, refLon) % 将经纬度转为以参考点为中心的本地东北坐标单位米 % 输入均为十进制度lat/lon 可支持向量 R 6371000; % 地球平均半径单位米 refLatRad deg2rad(refLat); refLonRad deg2rad(refLon); latRad deg2rad(lat); lonRad deg2rad(lon); dLat latRad - refLatRad; dLon lonRad - refLonRad; yNorth R * dLat; xEast R * cos(refLatRad) * dLon; end这个近似在局部是几何直观的它认为目标点的纬线平行于参考点纬线且球面曲率小到可以忽略。于是北向位移就是子午线方向上的弧长R * dLat东向位移近似为参考点所在纬线圈上的弧长R * cos(refLat) * dLon。有了东北坐标任意两点之间的欧氏距离直接就是dist sqrt((x1 - x2)^2 (y1 - y2)^2);也可以写成一个更友好的封装函数function d euclideanDistFromLatLon(lat1, lon1, lat2, lon2, refLat, refLon) % 基于本地切平面投影的两点欧氏距离单位米 % 如果未指定参考点默认取两点中点的纬度/经度 if nargin 6 refLat (lat1 lat2) / 2; refLon (lon1 lon2) / 2; end [x1, y1] enuFromLatLon(lat1, lon1, refLat, refLon); [x2, y2] enuFromLatLon(lat2, lon2, refLat, refLon); d sqrt((x1 - x2).^2 (y1 - y2).^2); end注意点同一个数据集里的所有点必须用同一个参考点投影。如果每个点都拿自己做参考点坐标基准不一致减法没有任何意义。参考点的选择原则很简单选数据覆盖范围的几何中心这样整个区域内的投影误差是均匀分布的。3.3 使用地图工具箱进行全局投影大跨度轨迹的备选如果数据的跨度到了几百公里甚至上千公里本地切平面的误差会开始明显累积。这时候有两种选择一是用大圆距离直接算二是做全局投影比如 UTM。MATLAB 从 R2019b 开始提供了projcrs配合projfwd的用法对于 UTM 投影来说非常方便% 假设数据大致在东经 114-120 度、北纬 29-40 度对应 UTM 50N % EPSG:32650 表示 WGS84 坐标系下的 UTM 50N 带 utmProj projcrs(32650, Authority, EPSG); [x, y] projfwd(utmProj, lat, lon);得到x, y已经是以米为单位的平面坐标直接算欧氏距离即可。326XX是北半球 UTM 带号327XX是南半球。国内通常用 49 到 52 带具体带号要根据经度换算带号约为floor((lon 180) / 6) 1。如果没有安装 Mapping Toolbox自实现 UTM 投影会涉及椭球参数、墨卡托级数展开这类比较繁琐的数学不建议手写。更务实的做法是如果跨带或者大跨度直接用 3.1 节的大圆距离或者把数据分块后分别用本地切平面。3.4 批量实战一组 GPS 点到参考点的距离计算假设有一组经纬度点想计算它们到某个中心点的距离并对比投影距离与大圆距离的差异可以直接这样写pts [ 31.2304, 121.4737; 39.9042, 116.4074; 30.5728, 114.3162; 23.1291, 113.2644 ]; refLat 34; refLon 118; % 本地切平面投影 [x, y] enuFromLatLon(pts(:,1), pts(:,2), refLat, refLon); dProj sqrt(x.^2 y.^2); % 大圆距离基准 dTrue haversineDist(refLat, refLon, pts(:,1), pts(:,2), km); % 查看差异 disp(table(dProj/1000, dTrue, (dProj/1000 - dTrue)./dTrue*100, ... VariableNames, {投影距离km, 大圆距离km, 相对误差%}));在几百公里跨度下投影距离与大圆距离的差异通常在 0.1% 到 0.5% 之间表现是相当可用的。如果误差超预期大概率是参考点选得太偏或者数据跨度过大。4. 投影误差的边界多少公里内可以放心用拿到代码之后心里得有一杆秤我的数据跨度多大用本地切平面还能不能顶住这一章就用脚本把边界测出来免得以后拍脑袋。4.1 基于参考点投影的距离误差随跨度变化以参考点北纬 30 度、东经 110 度为例让目标点在东西南北各个方向上偏移 0 到 10 度逐点比较本地切平面投影距离与大圆距离。lat0 30; lon0 110; northSteps -10:0.5:10; eastSteps -10:0.5:10; err zeros(length(northSteps), length(eastSteps)); for i 1:length(northSteps) for j 1:length(eastSteps) lat lat0 northSteps(i); lon lon0 eastSteps(j); [x, y] enuFromLatLon(lat, lon, lat0, lon0); dProj sqrt(x^2 y^2); dTrue haversineDist(lat0, lon0, lat, lon); % 已返回米 err(i,j) (dProj - dTrue) / dTrue * 100; end end % 观察最大绝对误差范围 fprintf(最大相对误差(%%): %.3f\n, max(abs(err(:))));基于球面几何可以推导正东、正西方向偏移时的误差主要由弧长与弦长的差异决定正北、正南方向偏移时在球体模型下的投影距离近似等于弧长误差反而很小。最大的误差出现在同时包含东西和南北方向的大跨度位移因为此时投影平面上的直线相当于球面大圆的一条弦距离会系统性偏小。实测下来在 5 度跨度约 550 公里内最大相对误差通常在 1% 到 2% 以内在 2 度跨度约 220 公里内误差普遍小于 0.2%。4.2 高纬度与参考点选择对误差的影响本地切平面投影的东向系数是R * cos(refLat)如果你的目标点实际纬度与参考点纬度相差很大这个系数就会失真。比如参考点在北纬 30 度目标点在北纬 60 度目标点附近真实的 1 度经线长度只有约 55.7 公里但投影公式按cos(30°)算出来是约 96.6 公里东西方向会偏差接近 74%。这种极端情况很容易遇到全国范围的配送数据、跨省车辆轨迹如果硬要用同一个本地切平面误差会非常难看。所以我的经验判断是50 公里以内本地切平面误差通常小于 0.05%放心用。200 公里以内误差一般小于 0.3%绝大多数业务可以接受。500 公里以内误差可能到 1% 到 2%需要结合精度要求判断。1000 公里以上不要再用本地切平面改用大圆距离或 UTM 投影。4.3 什么时候必须放弃平面近似需要特别警惕的是两类场景。第一类是高纬度区域比如北纬 70 度以上的极区附近投影变形本身就很敏感处理这类数据时优先选择极地方位投影或者干脆用球面距离。第二类是跨国跨洋的长距离轨迹本地切平面已经完全失去意义此时投影距离与真实大圆距离的偏差可达数个百分点而且由于地球曲率投影坐标下的直线方向也会失真。如果既想要平面坐标来跑下游算法又想要大跨度距离精度实际工程里有个折中方案按区域分块每块单独建立本地切平面参考点。比如按 2 度 x 2 度的格网切块每个网格内的点用网格中心做参考点投影既能延续欧氏空间计算又能把误差控制在很小的范围内。5. 坐标转换里最容易翻车的五个细节即使掌握了上面的核心代码实际工程中还是很常见各种奇葩问题。下面这几个细节几乎每个做轨迹数据处理的人都踩过。5.1 坐标顺序lat/lon 还是 lon/latMATLAB 地图相关函数比如geoplot的习惯是纬度在前、经度在后也就是(lat, lon)。但很多数据源恰好相反比如 GeoJSON 的标准是[lon, lat]部分设备输出也是先把经度放在前面。一旦写反所有计算全部错乱而且很难通过肉眼发现。更隐蔽的是经度规范化问题。比如一个点在东经 179.9 度另一个点在西经 -179.9 度它们的实际经度差只有 0.2 度但如果直接做减法会得到 359.8 度算出来距离远得离谱。这种情况在跨太平洋数据里特别常见一定要先做规范化function dlon normalizedLon(dLonDeg) dlon mod(dLonDeg 180, 360) - 180; end5.2 度分秒与十进制度混用很多老式 GPS 设备输出的是度分格式比如31.2304可能代表 31 度 23.04 分而计算时必须用十进制度。直接把度分当成十进制度参与计算误差大概是 60 倍量级距离一下子偏出去十分离谱。转换公式是十进制度 度 分/60 秒/3600。MATLAB 里也有现成函数dms2degrees但关键是入库和清洗阶段就要统一不能在计算时才发现单位不一致。5.3 坐标系不统一数据源之间可能存在系统性偏移GPS 原始数据通常是 WGS84 坐标但国内很多地图服务返回的坐标是经过偏移变换的本地坐标系两者之间有几十米到几百米的系统性偏差。如果在同一个工程里混用了不同坐标系的数据去做距离筛选看起来每个点单独没问题但两两之间的距离会引入不可控误差边界过滤、聚合判断全都不可信。我处理这种问题时都会要求数据接入的第一步做基准校验拿同一个物理位置在不同数据源里的坐标做差确认坐标系是否统一。如果确实不一致必须用官方或公开的转换算法统一到同一个基准后再入库。注意转换过程也有精度损失定量评估后再决定是否接受。5.4 投影带号与参考点选错用 UTM 投影时如果数据跨了两个带强行使用同一个带号距离会异常偏大。比如东经 117 度和东经 123 度分别落在 UTM 50N 和 51N直接用 50N 投影两个点横坐标会被拉开几千公里显然不可用。处理方式要么分带投影要么回到大圆距离。本地切平面的参考点选择同样关键。参考点应该放在数据覆盖的中心而不是数据集边界点或者随便一个点。参考点放偏投影误差只在一侧膨胀距离计算会表现出明显的系统性偏差。用中心点时误差围绕真实值均匀分布而且最大值最小化。5.5 地球半径取值的玄学代码里我统一用 6371000 米这是国际大地测量与地球物理联合会推荐的平均半径。有的实现会用地球赤道半径 6378137 米有的会用 6356752 米极半径三者相差最多约 0.7%。对一般工程来说0.1% 到 0.7% 的差距几乎无感知但在高精度基线计算里就会成为问题。如果需要更高精度建议直接换用基于椭球模型的 Vincenty 公式而不是纠结平均半径到底取哪个数。MATLAB 里如果安装了 Mapping Toolboxdistance函数默认就用椭球模型在长距离上的表现比球面 Haversine 更可靠但短距离内两者差距微乎其微。最后说点我自己的习惯。我处理一批轨迹或者配送点位数据时第一步永远是看覆盖范围和用途如果所有点集中在几十公里范围内直接本地切平面投影参考点选数据范围中心然后放心大胆地用欧氏距离如果跨到省一级就考虑 UTM 或者大圆距离如果既有投影需求又要跨大范围就分块建局部坐标系绝不在全局硬切。另外我每次写完投影代码都会顺手做一次“基准校验”随机抽几对点用haversineDist和投影距离对比偏差超过预期就说明参考点或者带号有问题早发现早解决。这套流程帮我避开了很多坐标相关的坑希望能给同样被经纬度折磨的人一点参考。