Оценка положения приемника GNSS с моделируемыми спутниковыми созвездиями
Отслеживайте положение наземного транспортного средства с помощью приемника Глобальной навигационной спутниковой системы (GNSS). Спутники моделируются с помощью satelliteScenario (Satellite Communications Toolbox) объект, обработка спутниковых сигналов приемника моделируются с помощью lookangles и pseudoranges функций, и положение приемника оценивается с помощью receiverposition функция.
Этот пример требует Navigation Toolbox™.
Обзор
Этот пример фокусируется на космическом сегменте, или созвездиях, и датчике GNSS для системы statellite. Чтобы получить оценку положения приемника GNSS, навигационный процессор требует положения спутника от космического сегмента и псевдообласти значений от процессора диапазона в получателе.
Задайте параметры симуляции
Загрузите MAT-файл, который содержит положение «земля-правда» и скорость наземного транспортного средства, перемещающегося к кампусу Natick, 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")
Симулируйте положения спутников с течением времени
The 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. Термин pseudorange используется, потому что это значение расстояния вычисляется путем умножения различия во времени между временем текущего приемника и временными метками спутникового сигнала на скорость света.
% 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")
