简介：这是一套面向导航定位、测绘与自动驾驶方向学习者的MATLAB GPS定位算法仿真程序，围绕伪距测量、载波相位与最小二乘定位解算等核心原理展开，适合具备一定MATLAB基础、希望从理论走向工程实现的本科生与研究人员。压缩包共129个文件，约2.43MB，以93个.m脚本为主体，配合观测数据文件（.01o、.01n、.09o等）、导航电文（.nav）、.dat数据、.eps与.png图表及.pdf说明文档，覆盖信号模拟、接收机建模、信道延迟、数据解码到定位解算的完整链路。资源中附带的观测与导航文件可直接用于跑通解算流程，便于读者对照代码理解伪距与相位测量、误差源分析及坐标解算步骤，并在此基础上调整参数、验证新算法或模拟不同环境下的定位效果。目前已有1107人学习下载，适合作为课程设计、毕业设计或算法预研的参考素材。
1. 从一串伪距到一条轨迹：matlab_gps 定位算法仿真程序到底在算什么
打开接收机日志，你会看到每个历元有一堆以米为单位的伪距、卫星坐标和钟差，但真正想要的只是「我在哪」。matlab_gps 定位算法仿真程序要干的事，就是把这份原始观测变成可复现的定位结果：从读取 RINEX 观测文件、计算卫星位置、构造观测方程，到最小二乘或卡尔曼滤波解出接收机坐标与钟差，全程用 MATLAB 脚本跑通。它解决的不是「造一个 GPS 模块」，而是让你在没有硬件、没有天空视野的情况下，把导航定位解算原理吃透，能改参数、能注入误差、能看残差。适合两类人：一类是刚接触卫星导航、想搞懂伪距单点定位怎么落地的新手；另一类是要验证新算法（比如抗差估计、组合导航）但不想每次都上实测数据的熟手。热词里反复出现的 matlab、gps、定位算法、仿真程序、导航定位解算，本质都指向同一件事——用可调试的代码替代黑匣子接收机，把定位链路拆开看。
2. 伪距单点定位的数学骨架：观测方程怎么列、未知数怎么定
2.1 从伪距观测到定位方程
GPS 伪距观测量的基本模型是：接收机在某一历元测得的伪距 ρ，等于接收机到卫星的几何距离加上接收机钟差、卫星钟差、电离层与对流层延迟，再叠加测量噪声。写成标量形式：
ρ_i = || r_sat,i − r_rcv || + c·δt_rcv − c·δt_sat,i + I_i + T_i + ε_i
其中 r_sat,i 是第 i 颗卫星在地心地固系（ECEF）下的位置，r_rcv 是接收机位置，δt_rcv 是接收机钟差，δt_sat,i 是卫星钟差，I 和 T 分别是电离层、对流层延迟。仿真程序要做的第一件事，就是把卫星钟差、电离层、对流层这些「已知量」从伪距里扣掉，得到校正后伪距，再对几何距离和接收机钟差做估计。
这里有个容易翻车的点：卫星位置必须和伪距在同一时刻、同一坐标系下。卫星在信号发射时刻的位置，要经过地球自转改正才能和接收机在接收时刻的 ECEF 坐标对齐。很多初学者直接拿广播星历算出的卫星位置去减接收机坐标，结果定位偏差几十米，还以为是算法错了，其实是时间系统没对齐。
2.2 未知数与方程个数
单点定位的未知数有四个：接收机 ECEF 坐标 X、Y、Z，以及接收机钟差对应的距离项 c·δt_rcv。每颗卫星提供一个方程，所以至少需要 4 颗卫星才能解出唯一解。实际仿真中通常有 6 到 12 颗可见星，方程数大于未知数，用最小二乘求解。
线性化是绕不开的一步。几何距离对接收机坐标是非线性的，需要在近似位置 (X0, Y0, Z0) 处做一阶泰勒展开，得到设计矩阵（方向余弦矩阵）H。H 的每一行是接收机到卫星的单位视线向量和 1（对应钟差项）。这一步在 MATLAB 里就是几行矩阵运算，但近似位置的选取会影响收敛速度，一般用上一历元解或地心坐标作为初值。
2.3 最小二乘解算的核心代码
下面这段是单历元最小二乘定位的核心，输入是校正后伪距、卫星 ECEF 坐标和接收机近似位置，输出是位置增量和钟差。
function [pos, dtr, H, res] = ls_spp(pr_corr, sat_ecef, pos0, max_iter) % pr_corr : Nx1 校正后伪距 (m) % sat_ecef : Nx3 卫星 ECEF 坐标 (m) % pos0 : 1x3 接收机近似位置 (m) % max_iter : 最大迭代次数 % pos : 1x3 解算位置 (m) % dtr : 接收机钟差对应的距离 (m) % H : 设计矩阵 % res : 伪距残差 pos = pos0(:)'; dtr = 0; for iter = 1:max_iter N = size(sat_ecef, 1); H = zeros(N, 4); y = zeros(N, 1); for i = 1:N dx = sat_ecef(i,1) - pos(1); dy = sat_ecef(i,2) - pos(2); dz = sat_ecef(i,3) - pos(3); rho0 = sqrt(dx^2 + dy^2 + dz^2); % 方向余弦 H(i,1) = -dx / rho0; H(i,2) = -dy / rho0; H(i,3) = -dz / rho0; H(i,4) = 1; % 观测残差：校正伪距减去几何距离和当前钟差 y(i) = pr_corr(i) - rho0 - dtr; end % 最小二乘 dx_est = (H' * H) \ (H' * y); pos = pos + dx_est(1:3)'; dtr = dtr + dx_est(4); if norm(dx_est(1:3)) < 1e-3 break; end end % 最终残差 res = y - H * dx_est; end逻辑说明：每次迭代用当前估计位置计算几何距离和方向余弦，构造残差向量 y，解出位置增量和钟差增量，累加到当前估计上。收敛判据是位置增量小于 1 毫米。参数方面，max_iter 一般设 10 就够，因为最小二乘收敛很快；pos0 可以用地心坐标 (0,0,0) 或上一历元结果，用上一历元能减少一次迭代。注意 H 的第四列全为 1，对应钟差项，这是伪距定位的标准形式。
2.4 卫星位置计算与地球自转改正
卫星位置由广播星历参数计算，MATLAB 里可以写一个独立的函数，输入星历参数和信号发射时刻，输出 ECEF 坐标。关键步骤包括：计算半长轴、平均角速度、平近点角、偏近点角（开普勒方程迭代）、真近点角、升交点角距、轨道倾角、升交点经度等。开普勒方程 E = M + e·sin(E) 用不动点迭代或牛顿法解，一般 5 次迭代内收敛。
地球自转改正是另一个必做项。信号从卫星传到接收机大约需要 70 毫秒，这期间地球自转了约 0.3 毫角秒，对应地面距离约 30 米。改正方法是将卫星坐标绕 Z 轴旋转 ω·τ，其中 ω 是地球自转角速度，τ 是信号传播时间。传播时间用伪距除以光速近似即可，迭代一次就够。
function sat_ecef_corr = earth_rotation_corr(sat_ecef, pr) % sat_ecef : 卫星 ECEF 坐标 (m) % pr : 伪距 (m) c = 299792458; omega = 7.2921151467e-5; tau = pr / c; theta = omega * tau; R = [cos(theta), sin(theta), 0; -sin(theta), cos(theta), 0; 0, 0, 1]; sat_ecef_corr = (R * sat_ecef')'; end这段代码把卫星坐标绕 Z 轴旋转了 ω·τ 角度。参数 omega 是 WGS-84 地球自转角速度，tau 用伪距近似传播时间。注意旋转方向：信号传播期间地球在转，接收机坐标系在转，所以要把卫星坐标转到接收时刻的坐标系，旋转矩阵的符号不能搞反，否则误差会翻倍。
3. 用 MATLAB 把仿真程序跑起来：数据、流程与参数配置
3.1 仿真数据从哪来
没有实测数据也能跑仿真。常见做法有三种：一是用 RINEX 观测文件和广播星历文件，这些可以从公开的 GNSS 数据中心获取，MATLAB 里用文本解析函数读取；二是自己生成仿真数据，给定接收机真实位置和卫星星座，正向计算伪距，再叠加噪声和误差；三是用 MATLAB 的卫星导航工具箱（如果有授权）直接调用。我一般会先用第二种方式，因为可控性最强，能精确知道真实位置，方便验证算法正确性。
自己生成伪距的流程：设定接收机真实位置（比如某点的 ECEF 坐标），用广播星历或简化星座模型算出每颗卫星在多个历元的位置，计算几何距离，加上接收机钟差、卫星钟差、电离层延迟、对流层延迟和噪声，得到仿真伪距。这样每个历元的真实位置已知，解算结果可以直接对比。
3.2 主流程脚本的骨架
一个完整的仿真程序主流程包括：读取或生成数据、逐历元循环、卫星位置计算、误差改正、最小二乘解算、结果存储与可视化。下面是一个简化的主脚本框架。
% 主流程：伪距单点定位仿真 clear; clc; % 1. 生成仿真数据 [pr_all, sat_all, true_pos] = gen_sim_data(); % 2. 初始化 n_epoch = size(pr_all, 1); pos_est = zeros(n_epoch, 3); dtr_est = zeros(n_epoch, 1); pos0 = true_pos(1,:) + [10, 10, 10]; % 近似位置加偏差 % 3. 逐历元解算 for k = 1:n_epoch pr = pr_all(k, :)'; sat = squeeze(sat_all(k, :, :)); % 剔除无效卫星 valid = ~isnan(pr) & all(~isnan(sat), 2); pr = pr(valid); sat = sat(valid, :); if length(pr) < 4 pos_est(k, :) = NaN; continue; end % 地球自转改正 for i = 1:length(pr) sat(i, :) = earth_rotation_corr(sat(i, :), pr(i)); end % 最小二乘 [pos, dtr] = ls_spp(pr, sat, pos0, 10); pos_est(k, :) = pos; dtr_est(k) = dtr; pos0 = pos; % 下一历元用当前解作为初值 end % 4. 可视化 figure; plot3(true_pos(:,1), true_pos(:,2), true_pos(:,3), 'b-', 'LineWidth', 1.5); hold on; plot3(pos_est(:,1), pos_est(:,2), pos_est(:,3), 'r--', 'LineWidth', 1.5); xlabel('X (m)'); ylabel('Y (m)'); zlabel('Z (m)'); legend('真实轨迹', '解算轨迹'); grid on;逻辑说明：主脚本先生成仿真数据，然后逐历元解算。每个历元先剔除无效卫星，再做地球自转改正，最后调用最小二乘函数。pos0 在历元间传递，用上一历元解作为下一历元初值，能加快收敛。可视化部分用三维轨迹对比，直观看出解算精度。
参数方面，gen_sim_data 里可以控制噪声大小、电离层延迟模型、对流层延迟模型。噪声一般设 0.5 到 2 米，电离层延迟在 L1 频段白天可达 5 到 15 米，对流层延迟在天顶方向约 2.3 米，随高度角变化。这些参数直接影响定位误差，调参时要有依据。
3.3 误差注入与精度评估
仿真程序的价值在于能单独控制每个误差源。比如想验证电离层延迟对定位的影响，就在生成伪距时只加电离层延迟，不加噪声，看解算结果偏差多少。想验证抗差算法，就在某几颗卫星的伪距上加粗差，看最小二乘和抗差估计的区别。
精度评估常用指标：水平误差、垂直误差、三维误差、均方根误差。水平误差是解算位置与真实位置在水平面上的距离，垂直误差是高程方向偏差。在 MATLAB 里就是几行减法加范数。
% 精度评估 err = pos_est - true_pos; hor_err = sqrt(err(:,1).^2 + err(:,2).^2); ver_err = abs(err(:,3)); rms_3d = sqrt(mean(sum(err.^2, 2))); fprintf('水平误差均值: %.2f m\n', mean(hor_err)); fprintf('垂直误差均值: %.2f m\n', mean(ver_err)); fprintf('三维 RMS: %.2f m\n', rms_3d);这段代码计算水平误差、垂直误差和三维 RMS。注意 ECEF 坐标下的 X、Y 不能直接当水平方向，严格来说要转到 ENU 局部坐标系再算水平误差。简化处理时可以用 X、Y 的范数近似，但精度要求高时要做坐标转换。
4. 避坑与排查：仿真程序跑不通时先看这几处
4.1 定位结果发散或迭代不收敛
现象：最小二乘迭代多次后位置增量不减小，或者解算位置飞到几万公里外。原因通常是设计矩阵 H 构造错误，比如方向余弦符号搞反，或者卫星坐标和接收机坐标不在同一坐标系。解决：检查 H 的每一行，方向余弦应该是 (sat - pos) / rho 取负号，因为残差是观测减计算。再检查卫星坐标是否做了地球自转改正，接收机坐标是否在 ECEF 系。
4.2 伪距残差普遍偏大
现象：解算完成后残差在几十米甚至上百米。原因可能是卫星钟差没扣、电离层对流层延迟没改正，或者伪距单位搞错（比如把米当成千米）。解决：逐项检查改正量，卫星钟差从广播星历第一数据块读取，电离层用 Klobuchar 模型，对流层用 Saastamoinen 模型。单位统一用米，光速用 299792458 m/s。
4.3 卫星数足够但解算失败
现象：可见星有 6 颗以上，但最小二乘报矩阵奇异或结果异常。原因可能是某几颗卫星几何分布太差，方向余弦矩阵接近奇异，或者有卫星坐标重复。解决：计算几何精度因子（GDOP），GDOP 大于 10 时定位精度很差，可以剔除低高度角卫星或等待几何构型改善。检查卫星坐标是否有重复行，重复行会导致矩阵秩亏。
4.4 时间系统不一致导致米级偏差
现象：定位结果整体偏移几米到几十米，但残差不大。原因可能是卫星位置对应的时间与伪距观测时间不一致，比如用了接收时刻的卫星位置而不是发射时刻。解决：明确时间标签，卫星位置在信号发射时刻计算，伪距对应接收时刻，传播时间用伪距除以光速迭代一次。GPS 时和 UTC 的跳秒也要注意，仿真中一般统一用 GPS 时。
4.5 MATLAB 版本与函数兼容性
现象：换一台机器跑，提示函数未定义或结果不同。原因可能是用了新版本才有的函数，或者矩阵运算在不同版本下有细微差异。解决：尽量用基础函数，避免依赖工具箱。如果用了 readtable、datetime 等较新函数，在旧版本上要替换。中文注释乱码问题在 matlab 2023 及更早版本常见，保存时选 UTF-8 编码，或者用英文注释。
5. 从单点定位到滤波与抗差：让仿真程序更接近真实场景
单历元最小二乘解算出来的轨迹通常毛刺很大，因为伪距噪声直接反映到位置结果上。真实接收机之所以输出平滑轨迹，是因为用了卡尔曼滤波。在仿真程序里加一个卡尔曼滤波器，状态量取位置、速度、钟差、钟漂，观测量用伪距，能明显改善轨迹平滑度。我一般会先用最小二乘跑一遍，看残差是否正常，再加滤波，对比滤波前后的误差曲线。
卡尔曼滤波的关键是过程噪声和观测噪声的设定。过程噪声反映接收机动态，静态场景可以设很小，动态场景要调大。观测噪声反映伪距精度，一般设 1 到 3 米。这两个参数调不好，滤波要么跟不上动态，要么平滑过度导致滞后。仿真程序的好处就是可以反复调，看不同参数下的轨迹和误差。
抗差估计是另一个值得加的功能。当某颗卫星伪距有粗差时，最小二乘会被带偏，抗差估计通过降低异常观测的权重来抑制影响。常用方法有 Huber 权函数、IGG3 方案。在 MATLAB 里就是在最小二乘迭代中加一个权矩阵，根据残差大小调整权重。
% 抗差最小二乘（Huber 权函数） function [pos, dtr] = robust_spp(pr, sat, pos0, max_iter, k0) pos = pos0(:)'; dtr = 0; for iter = 1:max_iter N = length(pr); H = zeros(N, 4); y = zeros(N, 1); for i = 1:N dx = sat(i,1) - pos(1); dy = sat(i,2) - pos(2); dz = sat(i,3) - pos(3); rho0 = sqrt(dx^2 + dy^2 + dz^2); H(i,1) = -dx / rho0; H(i,2) = -dy / rho0; H(i,3) = -dz / rho0; H(i,4) = 1; y(i) = pr(i) - rho0 - dtr; end dx_est = (H' * H) \ (H' * y); res = y - H * dx_est; % Huber 权 sigma = 1.4826 * median(abs(res - median(res))); w = ones(N, 1); for i = 1:N if abs(res(i)) > k0 * sigma w(i) = k0 * sigma / abs(res(i)); end end W = diag(w); dx_est = (H' * W * H) \ (H' * W * y); pos = pos + dx_est(1:3)'; dtr = dtr + dx_est(4); if norm(dx_est(1:3)) < 1e-3 break; end end end这段代码在标准最小二乘基础上加了 Huber 权。先算残差，用中位数绝对偏差估计尺度 sigma，残差超过 k0 倍 sigma 的观测降权。k0 一般取 1.345，这是 Huber 权函数的经典参数。注意权矩阵是对角阵，每次迭代重新计算。抗差估计对粗差很有效，但计算量比标准最小二乘大，仿真时可以根据需要选择。
验证滤波和抗差效果，最直接的方法是对比误差曲线。我习惯把最小二乘、卡尔曼滤波、抗差最小二乘的结果画在同一张图上，看水平误差和垂直误差的时序。如果滤波后误差反而变大，多半是过程噪声设得太小，滤波器过于信任动力学模型。如果抗差后误差没改善，检查粗差是否加在了高度角很低的卫星上，低高度角卫星本身观测质量就差，抗差权重可能不够。
最后说一个我踩过的坑：仿真数据生成时，卫星钟差和接收机钟差的符号。广播星历里的钟差是卫星钟相对于 GPS 时的偏差，伪距观测方程里是减去卫星钟差、加上接收机钟差。符号搞反，定位结果会整体偏移光速乘以钟差量级的距离，看起来像坐标系统错误，其实是钟差符号问题。每次改代码后，先用无噪声数据跑一遍，确认解算位置和真实位置一致，再加噪声和误差。这个习惯帮我省了很多排查时间。希望帮到你。
本文还有配套的精品资源，点击获取
