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 を呼び出す前に、航法データがない衛星を必ず事前に除外してください — 対応する航法データがない衛星が1つでもあるとエラーになります
• gnssoptions を使って必ず大気補正を適用してください — 補正なしのSPPでは、垂直方向に20 mを超える誤差が生じます
• 複数パネルの図には tiledlayout／nexttile を使用してください（subplot は使用しないでください）
• すべての軸に単位を付けてラベルを付け、各プロットにタイトルを付けてください
• DOPの品質を報告してください — HDOP > 6は精度が低下しており、HDOP > 20は使用できません
• Skyplotには十分なサイズのFigureが必要です — コンパクトな tiledlayout のタイル内に skyplot を配置しないでください。Skyplot専用の figure を使用するか、レイアウト内に配置する必要がある場合は、Skyplotのタイルが少なくとも560×560ピクセルになるようにFigureのサイズを設定してください。小さいタイルでは、衛星マーカーが重なって判読できなくなります。

ワークフロー

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. 仰角マスクを適用する — 1つのエポックから初期位置を計算し、lookangles を使って衛星の仰角を求め、マスクを下回る衛星を削除します
10. 位置を求める — 大気補正用の gnssoptions とともに receiverposition を呼び出します
11. 指標を報告する — エポックごとの指標（HDOP、VDOP、衛星数）と全体の指標（散布のRMS、DOP平均値）
12. 可視化する — skyplot、位置の時系列、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);

仰角マスクによるフィルタリング

仰角フィルタリングには位置推定値が必要です。1 エポックから初期位置を求めた後、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 フィルタリングは、測定値の時刻表ではなく観測時刻表に対して行います。 gnssmeasurements に渡す前に、S1C 列を使って obsData.GPS の行をフィルタリングしてください。
• 仰角マスクには位置の推定値が必要です。 位置を計算する前に仰角でフィルタリングすることはできません。まずフィルタリングなしで測位し、その後 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 Data Centers、CORS Network。

エラー処理

症状	原因	対処方法
gnssmeasurements エラー：「衛星 X の航法メッセージがありません」	観測された衛星に一致する航法データがありません	SatelliteID に対して intersect を使い、航法データに含まれる衛星だけに観測タイムテーブルを事前に絞り込みます
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.