Оцените положение приемника 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")
Симулируйте спутниковые положения в зависимости от времени
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
Оцените положение приемника из псевдообластей значений и спутниковые положения
Наконец, используйте спутниковые положения и псевдообласти значений, чтобы оценить положение приемника с 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")
Постройте абсолютную погрешность в оценке положения. Ошибка сглаживается движущейся медианой, чтобы сделать график более читаемым. Ошибка в 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")
