Оцените положение приемника GNSS с симулированными спутниковыми группировками

Отследите положение наземного транспортного средства с помощью симулированного приемника Глобальной навигационной спутниковой системы (GNSS). Спутники симулированы с помощью satelliteScenario Объект (Satellite Communications Toolbox), спутниковая обработка сигналов приемника симулирована с помощью lookangles и pseudoranges функции и положение приемника оцениваются с receiverposition функция.

Этот пример требует Navigation Toolbox™.

Обзор

Этот пример фокусируется на космическом сегменте, или спутниковых группировках и оборудовании датчика GNSS для спутниковой системы. Чтобы получить оценку положения приемника GNSS, процессор навигации требует спутниковых положений от космического сегмента и псевдообластей значений от располагающегося процессора в приемнике.

Задайте параметры симуляции

Загрузите MAT-файл, который содержит положение основной истины и скорость наземного транспортного средства, перемещающегося к Натику, кампусу MA MathWorks, Inc.

Задайте запуск, остановитесь, и шаг расчета симуляции. Кроме того, задайте угол маски или минимальный угол возвышения, приемника GNSS.

% Load ground truth trajectory.
load("routeNatickMA.mat","lat","lon","pos","vel","lla0");
recPos = pos;
recVel = vel;

% Specify simulation times.
startTime = datetime(2020,10,28,8,0,0,"TimeZone","America/New_York");
simulationSteps = size(pos,1);
dt = 1;
stopTime = startTime + seconds((simulationSteps-1)*dt);

% Specify mask angle.
maskAngle = 5; % degrees

% Convert receiver position from North-East-Down (NED) to geodetic
% coordinates. Requires Navigation Toolbox(TM).
receiverLLA = ned2lla(recPos,lla0,"ellipsoid");

% Set RNG seed to allow for repeatable results. 
rng("default");

Визуализируйте geoplot для траектории основной истины.

figure
geoplot(lat,lon)
geobasemap("topographic")
title("Ground Truth Trajectory")

Figure contains an axes object. The axes object contains an object of type line.

Симулируйте спутниковые положения в зависимости от времени

satelliteScenario объект позволяет вам задать начальные орбитальные параметры и визуализировать их использующий satelliteScenarioViewer Объект (Satellite Communications Toolbox). Этот пример использует satelliteScenario и MAT-файл начальными орбитальными параметрами, чтобы симулировать созвездия GPS в зависимости от времени. В качестве альтернативы вы могли использовать gnssconstellation объект, который симулирует спутниковые положения с помощью номинальных орбитальных параметров, и только текущего времени симуляции, необходим, чтобы вычислить спутниковые положения.

% Create scenario.
sc = satelliteScenario(startTime, stopTime, dt);

load("initialOrbitalParameters.mat","semiMajorAxis","eccentricity",...
    "inclination","rightAscensionOfAscendingNode",...
    "argumentOfPeriapsis","trueAnomaly","prnStr");

% Initialize satellites.    
satellite(sc,semiMajorAxis,eccentricity,inclination,...
    rightAscensionOfAscendingNode,argumentOfPeriapsis,trueAnomaly,...
    "Name",prnStr);

% Preallocate results.
numSats = numel(sc.Satellites);
allSatPos = zeros(numSats,3,simulationSteps);
allSatVel = zeros(numSats,3,simulationSteps);

% Save satellite states over entire simulation.
for i = 1:numel(sc.Satellites)
    [oneSatPos, oneSatVel] = states(sc.Satellites(i),"CoordinateFrame","ecef");
    allSatPos(i,:,:) = permute(oneSatPos,[3 1 2]);
    allSatVel(i,:,:) = permute(oneSatVel,[3 1 2]);
end

Вычислите псевдообласти значений

Используйте спутниковые положения, чтобы вычислить псевдообласти значений и спутниковую видимость в течение симуляции. Угол маски используется, чтобы определить спутники, которые отображаются к приемнику. Псевдообласти значений являются расстояниями между спутниками и приемником GNSS. Термин псевдообласть значений используется, потому что это значение расстояния вычисляется путем умножения разницы во времени между текущим приемником, показывают время и спутниковый сигнал с меткой времени скоростью света.

% Preallocate results.
allP = zeros(numSats,simulationSteps);
allPDot = zeros(numSats,simulationSteps);
allIsSatVisible = false(numSats,simulationSteps);

% Use the skyplot to visualize satellites in view.
sp = skyplot([], []);

for idx = 1:simulationSteps
    satPos = allSatPos(:,:,idx);
    satVel = allSatVel(:,:,idx);
    
    % Calculate satellite visibilities from receiver position.
    [satAz,satEl,allIsSatVisible(:,idx)] = lookangles(receiverLLA(idx,:),satPos,maskAngle);
    
    % Calculate pseudoranges and pseudorange rates using satellite and
    % receiver positions and velocities.
    [allP(:,idx),allPDot(:,idx)] = pseudoranges(receiverLLA(idx,:),satPos,recVel(idx,:),satVel);
    
    set(sp,"AzimuthData",satAz(allIsSatVisible(:,idx)), ...
        "ElevationData",satEl(allIsSatVisible(:,idx)), ...
        "LabelData",prnStr(allIsSatVisible(:,idx)))
    drawnow limitrate
end

Figure contains an object of type skyplot.

Оцените положение приемника из псевдообластей значений и спутниковые положения

Наконец, используйте спутниковые положения и псевдообласти значений, чтобы оценить положение приемника с receiverposition функция.

% Preallocate results.
lla = zeros(simulationSteps,3);
gnssVel = zeros(simulationSteps,3);
hdop = zeros(simulationSteps,1);
vdop = zeros(simulationSteps,1);

for idx = 1:simulationSteps
    p = allP(:,idx);
    pdot = allPDot(:,idx);
    isSatVisible = allIsSatVisible(:,idx);
    satPos = allSatPos(:,:,idx);
    satVel = allSatVel(:,:,idx);
    
    % Estimate receiver position and velocity using pseudoranges,
    % pseudorange rates, and satellite positions and velocities.
    [lla(idx,:),gnssVel(idx,:),hdop(idx,:),vdop(idx,:)] = receiverposition(p(isSatVisible),...
        satPos(isSatVisible,:),pdot(isSatVisible),satVel(isSatVisible,:));
end

Визуализация результатов

Постройте положение основной истины и предполагаемую позицию приемника по geoplot.

figure
geoplot(lat,lon,lla(:,1),lla(:,2))
geobasemap("topographic")
legend("Ground Truth","Estimate")

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent Ground Truth, Estimate.

Постройте абсолютную погрешность в оценке положения. Ошибка сглаживается движущейся медианой, чтобы сделать график более читаемым. Ошибка в x-и оси Y меньше, потому что существуют спутники по обе стороны от приемника. Ошибка в оси z больше, потому что существуют только спутники выше приемника, не ниже его. Ошибка изменяется в зависимости от времени, когда приемник перемещается, и некоторые спутники входят и из представления.

estPos = lla2ned(lla,lla0,"ellipsoid");
winSize = floor(size(estPos,1)/10);
figure
plot(smoothdata(abs(estPos-pos),"movmedian",winSize))
legend("x","y","z")
xlabel("Time (s)")
ylabel("Error (m)")
title("Position (NED) Error")

Figure contains an axes object. The axes object with title Position (NED) Error contains 3 objects of type line. These objects represent x, y, z.