简介面向GIS与航空航天领域的MATLAB坐标转换小工具用于解决地心地固坐标系ECEF与经纬度高度坐标系LLH间的相互转换问题同时提供ENU局部坐标到ECEF的变换能力适合测绘、遥感及导航相关专业的开发者参考。压缩包共2个文件均为.m脚本整体仅2KB包含xyz2llh.m和enu2xyz.m两个核心函数分别实现ECEF转LLH及ENU转ECEF代码精简且易于二次封装调用。已有480人学习下载适合需要快速完成坐标转换或理解椭球模型解算原理的MATLAB使用者。资源基于WGS84椭球模型展开脚本中涉及地球半长轴、半短轴及自转角速度等参数可帮助读者掌握从空间直角坐标到经纬度高程的完整推算流程并在此基础上扩展局部东-北-天坐标的应用场景。整体内容虽小但覆盖了坐标转换中两个常用方向对从事卫星定位、路径规划或遥感图像处理的人员具有直接参考价值。1. 从地固坐标系到经纬高两个 m 文件解决测绘数据的最常见坐标系翻转从事卫星导航、航空测量或 GIS 开发的工程师一定遇到过这种场景天线相位中心给的是 ECEF 直角坐标x、y、z 都是上千万米的大数可甲方图纸上要的却是经度、纬度和海拔反过来无人机机载基站把差分结果写成东向、北向、天向的局部坐标却要和全球卫星轨道比对。按教材公式手推时第一个坑就是地球不是球第二个坑是纬度解算迭代不收敛。这个压缩包里的xyz2llh.m和enu2xyz.m恰好覆盖两条路径前者把地球固连直角坐标转为经纬高后者把局部东-北-天坐标还原到地固坐标系。对刚入门大地测量或者在 MATLAB 里做 GNSS 数据预处理的人来说这就是第一版可用的基础工具。2. WGS84 椭球与迭代求解ECEF 到 LLH 的数学底座要能看懂xyz2llh.m先要理解它背后的椭球模型。2.1 为什么必须用 WGS84 椭球而不是球体地球自转导致赤道隆起赤道半径比极半径长 21 公里左右。用标准球体做正反算纬度误差很容易超过 0.1°地面距离误差能到 11 公里。WGS84 用长半轴 (a6378137) 米、扁率 (f1/298.257223563) 定义旋转椭球短半轴 (ba(1-f))。在 MATLAB 里我会直接把这些常量放到脚本头部避免每个函数重复赋值% wgs84 椭球参数放在公共脚本里 global a b f e2 a 6378137.0; % 长半轴单位 m f 1 / 298.257223563; % 扁率 b a * (1 - f); % 短半轴单位 m e2 (a^2 - b^2) / a^2; % 第一偏心率平方这段代码里的e2在后面法线长度计算中会反复出现。实际使用中我一般只记三个量a、f、e2b需要时从a和f推出来。下面是常见的 WGS84 查表值参数表达式WGS84 数值用途长半轴 a定义值6378137.0赤道半径扁率 f定义值1/298.257223563椭球形状短半轴 ba(1-f)6356752.31424518极半径第一偏心率平方 e²(a²-b²)/a²0.00669437999014法线长度计算到这里只解决了“用哪个椭球”的问题。真正麻烦的是一个空间点在 WGS84 椭球上的法线方向不一定过地心所以“纬度”并不是简单的atan2(z, p)。2.2 经度很简单纬度才是迭代的主战场经度可以直接写lon atan2(y, x); % 返回弧度atan2会自动处理四个象限得到的是 (-\pi) 到 (\pi) 的弧度值。若需要角度制乘180/pi即可。纬度分大地纬度和地心纬度。地心纬度是点位向量和赤道面的夹角用atan2(z, sqrt(x^2y^2))就能算但导航和测绘里要的是大地纬度即点沿椭球法线在椭球上的投影点对应的纬度。由于椭球扁平的影响二者最大可差 0.19°换算成地面距离接近 20 公里所以不能把地心纬度直接当结果。迭代的初值我一般用地心纬度p sqrt(x^2 y^2); % 点到旋转轴的距离 theta atan2(z * a, p * b); % 地心纬度作为初值 lat atan2(z e2 * b * sin(theta)^3, ... p - e2 * a * cos(theta)^3); % Bowring 一次逼近这里第一行p是空间点在赤道平面上的投影半径第三行用的是 Bowring 给出的初值公式它在大多数地面点上的误差小于 0.001 角秒后续迭代一两次就能稳定到毫米级。至于为什么带sin(theta)^3这是从椭球法线方程做泰勒展开得到的工程上先把公式当模板用没有问题想深挖可以去看大地测量教材。2.3 从初值到收敛法线长的双重代入有了lat的初值后要算椭球高h再反过来修正lat。法线长度 (N) 由纬度和第一偏心率共同决定for k 1:5 N a / sqrt(1 - e2 * sin(lat)^2); % 卯酉圈曲率半径 h p / cos(lat) - N; % 第一次求高度 lat atan2(z, p * (1 - e2 * N / (N h))); % 用高度修正纬度 end这个循环的逻辑是先按当前纬度算出法线长度N然后用水平面关系求出高度再更新纬度新的纬度会改变N所以迭代到h和lat稳定。把循环次数写死为 5 是离线批处理里最省心的做法MATLAB 跑几千个点也就是一眨眼的事。若改用while收敛阈值设为1e-9米时实际迭代次数一般不超过 4 次。h在这里是沿法线方向到椭球面的距离也就是椭球高不是海拔。海拔与椭球高之间的差叫大地水准面差距在 WGS84 基准下通常是几米到几十米要不要补偿取决于项目精度目标。这一层弄通了xyz2llh.m的每一行就都能对上号。3. xyz2llh.m 拆解从 x/y/z 到经纬高的完整实现3.1 一个可以直接运行的函数把上面的思路整理成独立函数就可以放到任意工程里使用function [lat, lon, h] xyz2llh(x, y, z, a, e2) % xyz2llh ECEF 直角坐标转经纬高 % 输入: x,y,z 为 WGS84 地固系坐标单位 m % a, e2 为椭球长半轴和第一偏心率平方 % 输出: lat, lon 单位为弧度h 为椭球高单位 m p sqrt(x^2 y^2); lon atan2(y, x); % 经度弧度 b a * sqrt(1 - e2); theta atan2(z * a, p * b); % 地心纬度初值 lat atan2(z e2 * b * sin(theta)^3, ... p - e2 * a * cos(theta)^3); for k 1:5 N a / sqrt(1 - e2 * sin(lat)^2); h p / cos(lat) - N; lat atan2(z, p * (1 - e2 * N / (N h))); end end这段代码把b用a * sqrt(1 - e2)就地算出来避免调用方再传一个参数。输入输出都写在注释里实际复用的时候a和e2用上一章的常量即可。如果你希望接口更简洁也可以把a、e2设为全局变量但那样不利于后续做多椭球对比。注意循环里h p / cos(lat) - N在纬度接近 ±90° 时会出现cos(lat)很小的情况。若点恰好在极轴附近p也接近 0两者相除依然稳定但是在浮点环境下当p小于 1e-8 时我一般直接返回lat sign(z) * pi/2并用h sign(z) * z - b单独算高度避免除零。3.2 单点测试与整文件批处理拿到xyz2llh.m后先做单点测试。以北京附近一个 WGS84 坐标为例经纬度为 116.397°E、39.908°N椭球高取 50 m。先用正变换公式生成 ECEF 坐标再放进函数看能否还原。场景输入坐标 (m)期望经度(°)期望纬度(°)期望高度(m)赤道本初子午线6378137, 0, 0000北京附近约 4391413, 523247, 4186720116.39739.90850南半球同经度上例 z 取反116.397-39.90850第二个测试点生成方式是这样的x_test 4391412.8; y_test 523246.9; z_test 4186719.5; [lat, lon, h] xyz2llh(x_test, y_test, z_test, a, e2); fprintf(lat%.9f deg, lon%.9f deg, h%.4f m\n, ... rad2deg(lat), rad2deg(lon), h);上面代码里fprintf的格式写 9 位小数是因为 1e-9 弧度的纬度对应大约 0.1 mm足以观察收敛噪声。如果你的数据是从 RINEX 文件里批量读出的上万条历元把循环套在外面逐点调用即可。xyz2llh本身没有状态量适合用arrayfun或者直接写for循环。3.3 精度控制在 1e-9 还是 1e-12我在不同项目里把收敛阈值从 1e-6 试到 1e-12 米结论是地面接收机后处理用 1e-6 就够如果做低轨卫星精密定轨需要 1e-9 以下再往下受双精度浮点尾数限制结果差异已经小到没有工程意义。更值得留意的是atan2的象限一致性。当 x 为负、y 接近 0 时经度会从 179.999° 跳到 -179.999°如果下游要把经度连续画图最好在函数出口统一规范到 02π 或 -ππ不要混用。若输入 z 为负南半球或高度为负矿井迭代初值公式依然成立但h p / cos(lat) - N的结果可能为负这时要检查N h是否接近 0否则纬度修正会出现异常。遇到地下场景我优先把循环条件写成abs(delta_h) 1e-6并对N h做下限保护。4. enu2xyz.m 实战局部东-北-天坐标如何回到地球质心4.1 ENU 基准点先有原点才有东和北ENU 坐标系是站心坐标系原点通常选在接收机天线相位中心、基站或者是城市控制点上。它和 ECEF 之间是“先平移、再旋转”的关系。假设基准点 O 的 ECEF 坐标为 (x0, y0, z0)移动站在 ENU 下偏移了 (e, n, u)要求移动站的 ECEF 坐标。这个场景在车载组合导航里非常常见惯导输出的是相对起点的东北天增量组合滤波器要和卫星定位帧合并必须先回到同一个地球固连框架里。常见的错误是只做平移、不转坐标轴这样误差会以公里级随距离放大。4.2 旋转矩阵要从经纬度反推ENU 的三个轴方向并不是固定单位向量而是依赖基准点的经纬度。东轴沿经度方向北轴沿纬度方向天轴沿椭球法线。把基准点的经纬度算出来后可以构造旋转矩阵function ecef enu2xyz(e, n, u, lat0, lon0, h0, a, e2) % enu2xyz 将局部东-北-天坐标转换为 ECEF 坐标 % 输入 e,n,u 为相对基准点的增量单位 m % lat0, lon0 为基准点大地经纬度弧度h0 为椭球高单位 m % 输出 ecef [x; y; z] [x0, y0, z0] llh2ecef(lat0, lon0, h0, a, e2); slat sin(lat0); clat cos(lat0); slon sin(lon0); clon cos(lon0); % ENU 到 ECEF 的旋转矩阵 R [ -slon, -slat*clon, clat*clon; clon, -slat*slon, clat*slon; 0, clat, slat ]; ecef R * [e; n; u] [x0; y0; z0]; end这段代码里的llh2ecef是xyz2llh的反函数按标准公式可以立即补上。旋转矩阵三列分别是东、北、天方向的单位向量在 ECEF 系下的投影。R列的含义对照如下R 列对应 ENU 轴数学表达式作用第 1 列东向 E[-sinlon; coslon; 0]把东向分量投到 ECEF第 2 列北向 N[-sinlatcoslon; -sinlatsinlon; coslat]把北向分量投到 ECEF第 3 列天向 U[coslatcoslon; coslatsinlon; sinlat]把天向分量投到 ECEF这里最关键的一点是R是正交矩阵所以反过来从 ECEF 到 ENU 时直接用R即可不必重新推方向余弦。这个性质在接收机 RTK 解算里经常被利用。4.3 与 xyz2llh 放在一起时的链式处理实际项目里常常是先拿基准点(x0,y0,z0)调xyz2llh得到(lat0,lon0,h0)再调enu2xyz把移动站 NEU 坐标转到 ECEF。整个过程我习惯封装成一个入口函数neu [12.345, -45.678, 0.123]; % 相对基准点偏移单位 m x0 4410000.0; y0 520000.0; z0 4190000.0; [lat0, lon0, h0] xyz2llh(x0, y0, z0, a, e2); P_ecef enu2xyz(neu(1), neu(2), neu(3), lat0, lon0, h0, a, e2); disp(P_ecef.); % 再反算一次验证 ENU 是否还原 neu_check ecef2enu(P_ecef, lat0, lon0, h0);这里ecef2enu是enu2xyz的反变换只要把R换成R即可。验证时neu_check与原本的neu差应在 1e-9 米量级。链路测过一遍后后面换任何基准点都敢直接用。5. 验证、精度陷阱与 MATLAB 里的进阶用法5.1 用 MATLAB 内置函数做对拍验证新版 MATLAB 的 Mapping Toolbox 提供了ecef2lla和lla2ecef。拿到xyz2llh.m以后第一步先和内置函数对拍。做法是生成一组全球均匀格网点lats -80:20:80; lons -180:30:180; for i 1:numel(lats) for j 1:numel(lons) [x, y, z] lla2ecef(lats(i), lons(j), 100); [lat_m, lon_m, h_m] xyz2llh(x, y, z, a, e2); err [abs(rad2deg(lat_m) - lats(i)), ... abs(rad2deg(lon_m) - lons(j)), ... abs(h_m - 100)]; if max(err) 1e-6 fprintf(large diff at (%g, %g): %e\n, ... lats(i), lons(j), max(err)); end end end这段代码把误差阈值放宽到 1e-6 度/米通常我们的迭代函数精度能到 1e-9 量级结果差异主要来自内置函数是否使用相同的椭球参数。若你是用 2023b 之后的版本授权和工具箱配置没问题时ecef2lla还能直接处理数组输入适合大规模点云一次转换。5.2 CSV/TXT 批处理给 coord 类工具很多测绘作业现场并不在 MATLAB 里看结果而是要把经纬高导成 CSV 或 TXT再交给 coord、PANDA 或自研软件。用writematrix能把转好的列一次性落盘data [x, y, z, rad2deg(lat), rad2deg(lon), h]; writematrix(data, converted.csv, Delimiter, ,);如果 coord 主界面只认带表头的 CSV 或 TXT则先用table包一层再writetableT table(x, y, z, rad2deg(lat), rad2deg(lon), h, ... VariableNames, {X_m,Y_m,Z_m,Lat_deg,Lon_deg,H_m}); writetable(T, coord_input.txt, Delimiter, \t);writetable默认把每行数据写到文本文件coord 类工具导入时通常能自动识别表头如果识别失败手动指定第一行为变量名即可。用writetable替代fprintf的好处是在中文 Windows 环境下字段类型保留更好少写不少循环。5.3 让外部脚本或 AI 工具也能调用转换函数如果你想用 Codex 或命令行工具直接操作 MATLAB 任务最省事的方式是把转换函数编成独立.m文件然后在命令行里执行matlab -batch load(data.mat); [lat,lon,h]xyz2llh(x,y,z,a,e2); save(out.mat,lat,lon,h)这样可以不打开桌面 GUICI 流水线、Python 脚本通过subprocess都能调用。核心原则是保持函数无全局副作用、输入输出全部通过参数传递这样无论是人还是工具都能稳定复用。本文还有配套的精品资源点击获取