Estimar la posición del receptor GNSS con constelaciones de satélites simuladas
Realice un seguimiento de la posición de un vehículo terrestre utilizando un receptor simulado del Sistema Global de Navegación por Satélite (GNSS). Los satélites se simulan usando el objeto satelliteScenario
(Satellite Communications Toolbox) , el procesamiento de la señal satelital del receptor se simula usando las funciones lookangles
y pseudoranges
, y la posición del receptor se estima con la función receiverposition
.
Visión general
Este ejemplo se centra en el segmento espacial o constelaciones de satélites y el equipo de sensores GNSS para un sistema de satélites. Para obtener una estimación de la posición del receptor GNSS, el procesador de navegación requiere las posiciones de los satélites del segmento espacial y las pseudodistancias del procesador de determinación de distancias del receptor.
Especificar parámetros de simulación
Cargue el archivo MAT que contiene la posición ground-truth y la velocidad de un vehículo terrestre que viaja hacia el campus de The MathWorks, Inc en Natick, MA.
Especifique el inicio, la parada y el tiempo de muestreo de la simulación. Además, especifique el ángulo de máscara, o ángulo de elevación mínimo, del receptor GNSS. Finalmente, especifique el archivo de mensajes de navegación RINEX para los parámetros orbitales iniciales del satélite.
% Load ground truth trajectory. load("routeNatickMA.mat","lat","lon","pos","vel","lla0"); recPos = pos; recVel = vel; % Specify simulation times. startTime = datetime(2021,6,24,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. receiverLLA = ned2lla(recPos,lla0,"ellipsoid"); % Specify the RINEX file. rinexFile = "GODS00USA_R_20211750000_01D_GN.rnx"; % Set RNG seed to allow for repeatable results. rng("default");
Visualice el geoplot
para la trayectoria ground-truth.
figure geoplot(lat,lon) geobasemap("topographic") title("Ground Truth Trajectory")
Simular posiciones de satélites a lo largo del tiempo
El objeto satelliteScenario
le permite especificar parámetros orbitales iniciales y visualizarlos utilizando el objeto satelliteScenarioViewer
. Este ejemplo utiliza el satelliteScenario
y un archivo RINEX con parámetros orbitales iniciales para simular las constelaciones GPS a lo largo del tiempo. Alternativamente, puede usar el objeto gnssconstellation
que simula las posiciones de los satélites usando parámetros orbitales nominales o un archivo RINEX, y solo se necesita el tiempo de simulación actual para calcular las posiciones de los satélites.
% Create scenario. sc = satelliteScenario(startTime, stopTime, dt); % Initialize satellites. navmsg = rinexread(rinexFile); satellite(sc,navmsg); satID = sc.Satellites.Name; % 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
Calcular pseudorangos
Utilice las posiciones de los satélites para calcular los pseudorangos y la visibilidad de los satélites a lo largo de la simulación. El ángulo de máscara se utiliza para determinar los satélites que son visibles para el receptor. Los pseudodistancias son las distancias entre los satélites y el receptor GNSS. El término pseudodistancia se utiliza porque este valor de distancia se calcula multiplicando la diferencia de tiempo entre la hora actual del reloj del receptor y la señal del satélite con marca de tiempo por la velocidad de la luz.
% Preallocate results. allP = zeros(numSats,simulationSteps); allPDot = zeros(numSats,simulationSteps); allIsSatVisible = false(numSats,simulationSteps); % Use the skyplot to visualize satellites in view. sp = skyplot([],[],MaskElevation=maskAngle); 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",satID(allIsSatVisible(:,idx))) drawnow limitrate end
Estimar la posición del receptor a partir de pseudodistancias y posiciones de satélites
Finalmente, utilice las posiciones de los satélites y los pseudodistancias para estimar la posición del receptor con la función 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
Visualizar resultados
Trace la posición ground-truth y la posición estimada del receptor en un geoplot
.
figure geoplot(lat,lon,lla(:,1),lla(:,2)) geobasemap("topographic") legend("Ground Truth","Estimate")
Trace el error absoluto en la estimación de la posición. El error se suaviza mediante una mediana móvil para que el gráfico sea más legible. El error en los ejes xey es menor porque hay satélites a ambos lados del receptor. El error en el eje z es mayor porque sólo hay satélites encima del receptor, no debajo de él. El error cambia con el tiempo a medida que el receptor se mueve y algunos satélites aparecen y desaparecen de la vista.
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")