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 utilizando el objeto satelliteScenario (Satellite Communications Toolbox), el procesamiento de la señal del satélite del receptor se simula utilizando 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élite. 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 y la velocidad ground-truth de un vehículo terrestre que viaja hacia el campus de Natick, MA de The MathWorks, Inc.
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 de la ground-truth .
figure geoplot(lat,lon) geobasemap("topographic") title("Ground Truth Trajectory")
Simular posiciones de satélite a lo largo del tiempo
El objeto satelliteScenario le permite especificar parámetros orbitales iniciales y visualizarlos utilizando el objeto satelliteScenarioViewer. Este ejemplo utiliza satelliteScenario y un archivo RINEX con parámetros orbitales iniciales para simular las constelaciones GPS a lo largo del tiempo. Como alternativa, puede utilizar el objeto gnssconstellation que simula las posiciones de los satélites utilizando 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. Se utiliza el término pseudodistancia 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
Por último, utilice las posiciones del satélite 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
Grafique 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 x y y 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")
