matlab-compute-gnss-position

多星座 GNSS 定位

基于 RINEX v3 数据进行伪距单点定位（SPP）。支持 GPS、GLONASS、Galileo、BeiDou、QZSS、NavIC/IRNSS 和 SBAS。

**要求：**Navigation Toolbox R2026a 或更高版本。

适用场景

• 将 RINEX 观测文件／导航文件处理为位置解
• 计算多星座接收机位置
• 分析定位精度、DOP 和卫星几何分布
• 比较单星座与多星座的性能
• 评估处理选项（高度角掩码、C/N0、改正）的影响
• 教授或演示 GNSS 定位概念

不适用场景

• 载波相位定位（RTK、PPP、模糊度解算）
• 与 IMU 进行传感器融合——请使用 matlab-fuse-inertial-sensors 或 INS 滤波器
• 卫星轨道传播或场景仿真——请使用 Aerospace Toolbox
• 实时 NMEA 数据流处理——请直接使用 nmeaParser
• RINEX v4 文件——嵌套结构体格式需要采用不同的处理方式

必须遵循的规则

• 始终先使用 rinexinfo，确认文件为 RINEX v3，并检查可用的星座／观测码，再读取文件
• 每个星座都要使用正确的观测码——错误的代码会导致测量值为空（参见默认代码表）
• 调用 gnssmeasurements 前，先筛除缺少导航数据的卫星——任何卫星找不到匹配的导航数据都会导致报错
• 始终通过 gnssoptions 应用大气改正——未经改正的 SPP 垂直误差超过 20 m
• 多面板图形使用 tiledlayout/nexttile（切勿使用 subplot）
• 为所有坐标轴标注单位，并为每幅图添加标题
• 报告 DOP 质量——HDOP > 6 表示性能下降，HDOP > 20 表示不可用
• 天空图需要使用全尺寸图窗——切勿将 skyplot 放入紧凑型 tiledlayout 的单个图块中。请为天空图单独使用一个 figure；如果必须与其他图共享布局，则应调整图窗大小，确保天空图所在图块至少为 560×560 像素。在较小的图块中，天空图的卫星标记会难以辨认且彼此重叠。

工作流

1. 检查 RINEX 文件——对每个文件调用 rinexinfo，确认其为 v3、识别星座并列出可用观测码。使用 SatelliteSystem 和 Descriptors 遍历 ObservationTypes：

   info = rinexinfo('rover.obs');
   for i = 1:numel(info.ObservationTypes)
       fprintf('%s: %s\n', info.ObservationTypes(i).SatelliteSystem, ...
           strjoin(info.ObservationTypes(i).Descriptors, ', '));
   end
   

2. 确定观测码——确认默认代码存在于 RINEX 数据中；如果不存在，请从文件头中建议替代代码
3. 收集处理选项——星座选择、高度角掩码（度，默认值为 10）、C/N0 阈值（dB-Hz，默认值为 0）、观测码覆盖项
4. 读取数据——调用 rinexread 读取观测文件和导航文件
5. 提取测量值——使用确认的观测码，针对每个星座调用 gnssmeasurements；预先筛选观测数据，移除找不到匹配导航数据的卫星
6. 按 C/N0 筛选——调用 gnssmeasurements 前，从观测时间表中移除信号强度列低于阈值的行
7. 验证星座——先分别解算各个星座；如果某个星座产生的位置明显错误（误差 >1 km、高度异常），则将其从组合中排除
8. 组合星座——垂直拼接已验证的测量时间表：combinedMeas = [gpsMeas; galMeas; ...]
9. 应用高度角掩码——使用一个历元计算初始位置，再用 lookangles 获取卫星高度角，并移除低于掩码的卫星
10. 解算位置——使用带有大气改正选项的 gnssoptions 调用 receiverposition
11. 报告指标——逐历元报告（HDOP、VDOP、卫星数量），并报告汇总指标（散布 RMS、平均 DOP）
12. 可视化——绘制天空图、位置时间序列、DOP 和卫星数量面板
13. （可选）真值对比——如果用户提供了 LLA，则通过 lla2enu 计算 ENU 误差

默认观测码（优先选择频段 1）

优先在所有星座中使用频段 1 观测码。 一致使用频段 1 观测码可以避免频间偏差，从而获得更稳健的多星座解。

星座	RINEX 字段	默认码	C/N0 列	备注
GPS	.GPS	"C1C"	S1C	L1 C/A（频段 1）
GLONASS	.GLONASS	"C1C"	S1C	L1 C/A（频段 1）
Galileo	.Galileo	"C1C"	S1C	E1；若 C1C 不可用，则使用 "C1X"
BeiDou	.BeiDou	"C1P" 或 "C1X"	S1P 或 S1X	优先选择频段 1（C1P、C1X 或任何 C1*）
QZSS	.QZSS	"C1C"	S1C	L1 C/A（频段 1）
NavIC/IRNSS	.NavIC	"C5A"	S5A	L5 SPS（唯一可用频段）
SBAS	.SBAS	"C1C"	S1C	L1 C/A（频段 1）

C/N0 列模式： 将观测码中的 C 前缀替换为 S（例如，C1C → S1C、C1P → S1P）。

频段 1 观测码的 NaN 比率可能高于更高频段的替代码，但在多星座解中可提供更稳健的位置。请始终使用 rinexinfo 检查观测码是否可用，并优先选择频段 1。

处理默认值

参数	默认值	单位	有效范围
星座选择	所有可用星座	—	支持系统的任意子集
高度角掩码	10	度	0–90
C/N0 阈值	0（不筛选）	dB-Hz	0–60
观测码	各星座的默认值	—	必须存在于 RINEX 观测字段中

在城市环境中，高度角掩码和 C/N0 筛选并非相互独立。 C/N0 筛选通常会移除相同的低高度角、多路径影响卫星。两者同时使用，效果很少优于单独使用 C/N0 筛选。接收机不报告信号强度时，高度角掩码最有效。

主要函数

函数	用途
rinexinfo	检查 RINEX 文件元数据（版本、系统、观测类型），无需读取文件
rinexread	读取 RINEX v3 观测文件和导航文件
gnssmeasurements	从观测数据和导航数据中提取伪距、卫星位置和时钟偏差
gnssoptions	配置大气修正和偏差精度
gnssIonosphere	Klobuchar 电离层延迟模型
gnssTroposphere	Saastamoinen 对流层延迟模型
receiverposition	加权最小二乘位置解；返回 [pos, vel, hdop, vdop, info]
lookangles	根据接收机位置计算卫星方位角、高度角和可见性
skyplot	按星座分组的卫星位置极坐标图
lla2enu	大地坐标到本地 ENU 坐标的转换（用于计算地面真值误差）

用法模式

仅使用 GPS 定位

dataDir = fullfile(matlabroot, 'toolbox', 'nav', 'positioning', ...
    'core', 'positioningdata');

obsData = rinexread(fullfile(dataDir, ...
    'GODS00USA_R_20211750000_01H_30S_MO.rnx'));
gpsNav = rinexread(fullfile(dataDir, ...
    'GODS00USA_R_20211750000_01D_GN.rnx'));

gpsMeas = gnssmeasurements(obsData.GPS, gpsNav.GPS);
opts = gnssoptions( ...
    Ionosphere=gnssIonosphere("klobuchar"), ...
    Troposphere=gnssTroposphere("saastamoinen"));

[recPos, recVel, hdop, vdop, info] = receiverposition(gpsMeas, opts);

多星座（GPS + Galileo）

galNav = rinexread(fullfile(dataDir, ...
    'GODS00USA_R_20211750000_01D_EN.rnx'));

gpsMeas = gnssmeasurements(obsData.GPS, gpsNav.GPS);
galMeas = gnssmeasurements(obsData.Galileo, galNav.Galileo, "C1X");

combinedMeas = [gpsMeas; galMeas];
[recPos, recVel, hdop, vdop] = receiverposition(combinedMeas, opts);

预先筛除缺少导航数据的卫星

BeiDou 和其他星座可能观测到没有匹配导航消息的卫星。如果任何卫星缺少导航数据，gnssmeasurements 就会报错。

bdsNav = rinexread(fullfile(dataDir, ...
    'GODS00USA_R_20211750000_01D_CN.rnx'));

% Find satellite IDs present in both obs and nav
obsIDs = unique(obsData.BeiDou.SatelliteID);
navIDs = unique(bdsNav.BeiDou.SatelliteID);
validIDs = intersect(obsIDs, navIDs);

% Filter observation data to valid satellites only
bdsObs = obsData.BeiDou(ismember(obsData.BeiDou.SatelliteID, validIDs), :);
bdsMeas = gnssmeasurements(bdsObs, bdsNav.BeiDou, "C1X");  % Use C1P or C1X or any C1* (band 1 preferred)

C/N0 信号强度筛选

在提取测量值前筛除弱信号。C/N0 列名遵循以下模式：将观测码中的 C 前缀替换为 S。

% Filter GPS observations with C/N0 < 30 dB-Hz
cn0Threshold = 30;
gpsObs = obsData.GPS;
gpsObs = gpsObs(gpsObs.S1C >= cn0Threshold, :);  % S1C corresponds to C1C

gpsMeas = gnssmeasurements(gpsObs, gpsNav.GPS);

高度角掩码筛选

高度角筛选需要位置估计。先使用一个历元获取初始定位结果，再应用 lookangles 筛除低高度角卫星。

% Remove rows with non-finite satellite positions (can occur in multi-GNSS data)
validRows = all(isfinite(combinedMeas.SatellitePosition), 2);
combinedMeas = combinedMeas(validRows, :);

% Get initial position from first epoch (unfiltered)
[initPos] = receiverposition(combinedMeas);
recLLA = initPos(1,:);  % [lat lon alt]

% Get unique epochs
epochs = unique(combinedMeas.Time);
filteredMeas = timetable();

for i = 1:numel(epochs)
    epochMeas = combinedMeas(combinedMeas.Time == epochs(i), :);
    satPos = epochMeas.SatellitePosition;
    [~, el, vis] = lookangles(recLLA, satPos, elevationMask);
    filteredMeas = [filteredMeas; epochMeas(vis, :)]; %#ok<AGROW>
end

大气修正

opts = gnssoptions( ...
    Ionosphere=gnssIonosphere("klobuchar"), ...
    Troposphere=gnssTroposphere("saastamoinen"));

[recPos, recVel, hdop, vdop, info] = receiverposition(combinedMeas, opts);
% info.ClockBias, info.ClockDrift, info.TDOP also available

质量指标

% Per-epoch metrics from receiverposition outputs
nSats = arrayfun(@(t) sum(combinedMeas.Time == t), unique(combinedMeas.Time));

% Aggregate scatter RMS (no ground truth needed)
meanPos = mean(recPos, 1, 'omitnan');
enu = lla2enu(recPos, meanPos, 'ellipsoid');
hRMS = rms(vecnorm(enu(:,1:2), 2, 2));  % Horizontal scatter
vRMS = rms(enu(:,3));                     % Vertical scatter

fprintf('Horizontal scatter RMS: %.2f m\n', hRMS);
fprintf('Vertical scatter RMS:   %.2f m\n', vRMS);
fprintf('Mean HDOP: %.2f\n', mean(hdop, 'omitnan'));
fprintf('Mean VDOP: %.2f\n', mean(vdop, 'omitnan'));

选择 C/N0 阈值。 测试一组较小的阈值（例如 25、30、35、40 dB-Hz），并使用相应的质量指标选出最佳值：

场景	质量指标	停止条件
有地面真值	相对于真值的 H RMS 和 V RMS	改善趋于平稳，或历元损失不可接受
无地面真值，静态接收机	位置散布 RMS	散布趋于平稳，或历元损失显著
无地面真值，动态接收机	HDOP 和有效历元数	HDOP 变差，或历元损失显著

C/N0 滤波通常对垂直精度的改善大于对水平精度的改善。滤波后持续存在的水平误差属于结构性多路径效应，SPP 无法解决。注意：散布 RMS 仅适用于静态接收机；对于移动接收机，它反映的是车辆轨迹，而非定位误差。

地面真值比较（可选）

% User provides known ground truth LLA
truthLLA = [39.020734 -76.826657 -8.95];  % Example: GODS station

enuError = lla2enu(recPos, truthLLA, 'ellipsoid');
error2D = vecnorm(enuError(:,1:2), 2, 2);
error3D = vecnorm(enuError, 2, 2);

fprintf('2D RMS error: %.2f m\n', rms(error2D));
fprintf('3D RMS error: %.2f m\n', rms(error3D));
fprintf('Mean East error:  %.2f m\n', mean(enuError(:,1)));
fprintf('Mean North error: %.2f m\n', mean(enuError(:,2)));
fprintf('Mean Up error:    %.2f m\n', mean(enuError(:,3)));

可视化

% Satellite skyplot — always in its own dedicated figure (never inside tiledlayout).
figure;
epochs = unique(combinedMeas.Time);
epoch1 = combinedMeas(combinedMeas.Time == epochs(1), :);
satPos1 = epoch1.SatellitePosition;
[az, el] = lookangles(recPos(1,:), satPos1);
skyplot(az, el, string(epoch1.SatelliteID))
title('Satellite Skyplot — First Epoch')

% Multi-panel metrics (separate figure)
figure;
t = tiledlayout(2, 2);
title(t, 'GNSS Positioning Metrics')

nexttile
plot(epochs, hdop, 'b-', epochs, vdop, 'r-', 'LineWidth', 1.5)
legend('HDOP', 'VDOP')
ylabel('DOP')
title('Dilution of Precision')
grid on

nexttile
plot(epochs, nSats, 'k-', 'LineWidth', 1.5)
ylabel('Count')
title('Satellites Used per Epoch')
grid on

nexttile
plot(epochs, recPos(:,1), 'b-', 'LineWidth', 1.5)
ylabel('Latitude (deg)')
title('Position — Latitude')
grid on

nexttile
plot(epochs, recPos(:,3), 'b-', 'LineWidth', 1.5)
ylabel('Altitude (m)')
title('Position — Altitude')
grid on

迭代优化

保存多次运行的结果，以便比较不同的处理配置。

% Initialize or append to run history
if ~exist('runHistory', 'var')
    runHistory = table('Size', [0 6], ...
        'VariableTypes', ["string","double","double","double","double","double"], ...
        'VariableNames', ["Config","MeanHDOP","MeanVDOP","HorizRMS","VertRMS","MeanSats"]);
end

newRow = table(configLabel, mean(hdop,'omitnan'), mean(vdop,'omitnan'), ...
    hRMS, vRMS, mean(nSats), ...
    'VariableNames', runHistory.Properties.VariableNames);
runHistory = [runHistory; newRow];
disp(runHistory)

注意事项

• 在紧凑的 tiledlayout 中，skyplot 难以阅读。 当坐标区小于默认图窗尺寸时，卫星标记和标签会重叠。始终为 skyplot 单独创建一个图窗。如果必须将其放入 tiledlayout，请将图窗尺寸设为确保 skyplot 所在的平铺区域至少为 560×560 像素。
• 缺少导航数据时，gnssmeasurements 会报错。 如果任何观测卫星 ID 不在导航时间表中，该函数就会报错。始终使用卫星 ID 上的 intersect，预先筛选观测数据，仅保留具有匹配导航消息的卫星。
• 多 GNSS 数据中的卫星位置可能为非有限值。 gnssmeasurements 返回的某些行在 SatellitePosition 或 Pseudorange 中可能包含 NaN/Inf。在调用 lookangles 或 receiverposition 前筛除这些行：validRows = all(isfinite(meas.SatellitePosition), 2);
• 合并前，分别验证各个星座。 伪距质量较差的星座（例如受到严重城市多路径效应影响）可能会破坏整个多星座解算结果。先单独测试每个系统——如果单独解算的位置误差很大（> 1 km 或高度异常），请将其排除或尝试使用其他观测码。
• 星座单独解算失败并非总能恢复。 即使其他星座在同一地点给出一致的位置，任何星座也可能产生极其错误的高度（例如 ±数十 km）。根本原因可能是接收机特有的伪距偏差、卫星几何分布不佳或时钟改正不完整——仅凭 RINEX 文件不一定能诊断出来。如果某星座未通过单独验证，且 C/N0 筛选也无法解决问题，请将其排除，而不要继续调整参数。
• 所有星座都优先使用频段 1 的观测码。 使用一致的频段 1 观测码（GPS/GLONASS/Galileo/QZSS 使用 C1C，BeiDou 使用 C1P 或 C1X）可获得更稳健的多星座解算结果。较高频段的观测码可能有更多数据，但会降低精度。
• 观测码因接收机而异。 根据接收机的不同，Galileo 可能使用 "C1X" 或 "C1C"。务必先使用 rinexinfo 检查。
• C/N0 筛选针对观测时间表，而非测量时间表。 将 obsData.GPS 行中的 S1C 列作为筛选依据，再传递给 gnssmeasurements。
• 仰角掩码需要位置估计值。 在计算位置之前，无法按仰角筛选。先进行一次不筛选的定位，然后使用 lookangles 找出低仰角卫星。
• 对于卫星少于 4 颗的历元，receiverposition 会返回 NaN。 检查输出中的 NaN 行，并向用户报告跳过的历元。
• RINEX v4 混合导航文件包含嵌套结构体（例如 navData.GPS.LNAV.EPH），与 v3 的扁平时间表不同。本技能仅针对 v3。
• GLONASS 轨道有效期为 ±15 分钟，而 GPS/Galileo/BeiDou 为 ±2 小时。请使用与观测时间接近的导航数据。
• 不同星座间的时间系统对齐（GPST、GLONASS UTC、GST、BDT）由 MATLAB 在时间表中通过 datetime 索引处理——无需手动转换。
• HDOP 质量等级： < 1 = 理想，1–2 = 极佳，2–5 = 良好，5–10 = 中等，10–20 = 一般，> 20 = 较差。超过 6 时标记为质量下降。

示例数据

用于测试的 RINEX v3 文件位于 fullfile(matlabroot, 'toolbox', 'nav', 'positioning', 'core', 'positioningdata')：

文件	类型	系统
GODS00USA_R_20211750000_01H_30S_MO.rnx	观测	GPS、GLONASS、Galileo、BeiDou
GODS00USA_R_20211750000_01D_GN.rnx	GPS 导航	GPS
GODS00USA_R_20211750000_01D_EN.rnx	Galileo 导航	Galileo
GODS00USA_R_20211750000_01D_RN.rnx	GLONASS 导航	GLONASS
GODS00USA_R_20211750000_01D_CN.rnx	BeiDou 导航	BeiDou

GODS 站的大致地面真值：[39.020734, -76.826657, -8.95]（来自 RINEX 文件头的 LLA）。

公开 RINEX 数据来源：IGS 数据中心、CORS 网络。

错误处理

症状	原因	修复方法
gnssmeasurements 错误：“Missing navigation message for satellite X”	观测到的卫星没有匹配的导航数据	使用 intersect 按 SatelliteID 预先筛选观测时间表，仅保留导航数据中存在的卫星
gnssmeasurements 返回空时间表	观测码错误或没有匹配的时间跨度	使用 rinexinfo 检查可用的观测码；确认观测数据与导航数据的时间范围有重叠
receiverposition 返回的结果全为 NaN	每个历元的卫星数量少于 4 颗	放宽高度角掩码、降低 C/N0 阈值，或增加星座数量
HDOP/VDOP 很高（> 20）	卫星几何分布不佳	增加星座数量、降低高度角掩码，并检查天空中受遮挡的区域
位置散布 > 50 m	未应用大气修正	添加带有 Klobuchar 电离层模型和 Saastamoinen 对流层模型的 gnssoptions
rinexread 读取文件时出错	RINEX v2、v4 或文件已损坏	使用 rinexinfo 检查版本；此技能要求 v3

技能反馈

如果此技能的指导有误、不完整或无法奏效，请在回复中输出纯文本 反馈标记，然后继续执行任务。

格式：NAV-SKILL-FEEDBACK [matlab-compute-gnss-position] [category]: what failed; what you did instead 类别：wrong-command（命令不正确或已过时）、missing-guidance（遗漏步骤）、 discovery-failure（建议的方法无效）、platform-issue（特定于环境的问题）

版权 2026 The MathWorks, Inc.