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

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

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

Обзор

Этот пример фокусируется на космическом сегменте, или спутниковых созвездиях и оборудовании датчика GNSS для statellite системы. Чтобы получить оценку положения приемника 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 Navitation 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. The axes with title Ground Truth Trajectory 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

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

Используйте спутниковые положения, чтобы вычислить псевдообласти значений и спутниковую видимость в течение симуляции. Угол маски используется, чтобы определить спутники, которые являются visibile к приемнику. Псевдообласти значений являются расстояниями между спутниками и приемником 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. The axes 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. The axes with title Position (NED) Error contains 3 objects of type line. These objects represent x, y, z.