Esta página se ha traducido mediante traducción automática. Haga clic aquí para ver la última versión en inglés.

Receptor de navegación GPS heredado de extremo a extremo mediante código C/A

Este ejemplo utiliza:

Abrir script en vivo

Este ejemplo muestra cómo estimar la posición del receptor del sistema de posicionamiento global (GPS) utilizando una forma de onda de banda base GPS multisatélite. Se utiliza el formato de intercambio independiente del receptor (RINEX) y un archivo de almanaque para modelar la constelación GPS y generar una forma de onda de banda base multisatélite. Simule el escenario del satélite para obtener posiciones relativas de los satélites con respecto al receptor modelado. Para este escenario de satélite, modele el desplazamiento Doppler, el retraso y la potencia de la señal recibida. Según estos cálculos, deteriore la señal de banda base generada con desplazamiento Doppler, retraso y ruido. Este ejemplo muestra cómo estimar la posición del receptor simulada a partir de esta señal de banda base de GPS deteriorada.

Introducción

En este ejemplo, comienza con un archivo RINEX y usa rinexread (Navigation Toolbox) para leer el archivo y proporcionar entrada a satelliteScenario (Satellite Communications Toolbox) para simular la constelación GPS. Con base en esta constelación de GPS y una posición determinada del receptor, se calcula el desplazamiento Doppler, el retraso y la pérdida de ruta de la señal desde cada satélite visible hasta el receptor. Con base en las efemérides y los datos de reloj y almanaque de los archivos RINEX y de almanaque, genere bits de datos según el estándar IS-GPS-200 [1]. Para los satélites que son visibles para el receptor, genere un código de adquisición aproximado (código C/A) y un código de precisión (código P). El código P se coloca en la rama en fase (I) y el código C/A se coloca en la rama de fase en cuadratura (Q) de la forma de onda de banda base. La señal de la rama I se atenúa en 3 dB como se indica en IS-GPS-200 [1]. Genere la forma de onda de banda base para todos los satélites visibles y pase esta forma de onda de banda base a través del canal de propagación como se muestra en la siguiente figura. Las características del canal de propagación para cada señal de satélite son únicas porque la posición y velocidad de cada satélite con respecto al receptor es única. En este ejemplo, modelará las características del canal de propagación (es decir, el desplazamiento Doppler, el retardo y el escalamiento de la potencia de la señal en función de la pérdida de la ruta de propagación) y agregará ruido térmico a la señal compuesta. Proporcione esta señal ruidosa como entrada al receptor GPS. Este ejemplo admite el almacenamiento de esta señal ruidosa en un archivo para que pueda probar su receptor.

La compleja forma de onda de banda base generada a partir del canal de propagación se procesa a través de un receptor GPS. El siguiente diagrama muestra los detalles de alto nivel de dicho receptor GPS. Primero, el módulo de sincronización inicial detecta satélites visibles. Estima valores aproximados del desplazamiento y retraso Doppler para los satélites visibles. Para cada uno de los satélites detectados, cree canales de recepción separados (que no debe confundirse con el canal de propagación). Cada canal receptor realiza seguimiento, sincronización de bits, sincronización de tramas, decodificación de datos, estimación del tiempo de transmisión, estimación de la posición del satélite y cálculo de pseudodistancias. Para estimar la posición del receptor, estime la posición de los satélites en el cielo y las distancias de los satélites al receptor. Estime la posición de los satélites a partir de los datos de efemérides decodificados según lo proporcionado en la Tabla 20-IV en IS-GPS-200 [1]. Calcule las distancias de los satélites al receptor estimando el tiempo de propagación y multiplicándolo por la velocidad de la luz. Para calcular este tiempo de propagación en el receptor, calcule la diferencia entre el tiempo del receptor y el tiempo de transmisión. Aunque lo ideal es calcular con precisión la hora del receptor, el receptor GPS no puede calcularla. Sólo puede estimar el tiempo de transmisión de la señal con alta precisión [2]. No se necesita una hora exacta del receptor para estimar la posición del receptor GPS. En un momento dado del receptor, calcule el tiempo de transmisión de cada señal de satélite y estime el pseudorango para cada satélite a partir de este valor. Los pseudorrangos no son rangos verdaderos y tienen grandes errores, pero como los errores son los mismos para todos los satélites, puedes tener en cuenta el error al resolver la posición del receptor y obtener la posición precisa del receptor [2]. Para obtener más detalles sobre la estimación de la posición del receptor, consulte la sección Estimación de la posición del receptor GPS.

Inicializar parámetros

Inicialice los parámetros necesarios para configurar y ejecutar la simulación del receptor GPS de un extremo a otro.

Inicialice la duración de los datos durante la cual se debe ejecutar este ejemplo. Normalmente, un receptor GPS necesita al menos 50 segundos de datos para estimar la posición del receptor. Simular este ejemplo con 50 segundos de datos requeriría mucho tiempo según los recursos informáticos disponibles. Para los fines de este ejemplo, establezca la duración de los datos en 3 segundos. Además, inicialice la frecuencia de muestreo para la generación de formas de onda de banda base.

simulatedDataDuration = 3; % In seconds
samplingRate = 10.23e6;    % In Hz

Si es necesario, habilite la propiedad WriteWaveformToFile para que pueda procesar esta forma de onda a través del receptor de su elección.

WriteWaveformToFile = false;

Este ejemplo procesa datos en fragmentos de 1 milisegundo. Es decir, el generador de formas de onda genera 1 milisegundo de datos y el receptor procesa ese milisegundo de datos. Este ejemplo funciona solo para un tiempo de paso de un milisegundo.

stepTime = 1e-3;                                 % In seconds
numSteps = (simulatedDataDuration/stepTime) + 1;

Inicialice los parámetros necesarios para configurar el canal. Este ejemplo genera una forma de onda GPS a partir de una constelación GPS simulada. Como punto de partida para simular la constelación GPS, lea el archivo RINEX.

rinexFileName = "GODS00USA_R_20211750000_01D_GN.rnx";
almanacFileName = "gpsAlmanac.txt";
rinexdata = rinexread(rinexFileName);

Inicialice la posición del receptor para modelar el canal de propagación en función de su ubicación. Para estimar la posición del receptor, modele con precisión el retardo de propagación. El retraso depende de la posición del receptor. Este ejemplo modela un receptor estacionario y no admite un receptor en movimiento.

rxlat = 39.021;  % Latitude in degrees north
rxlon = -76.827; % Longitude in degrees east; negative longitude is degrees west
rxalt  = 19;      % Altitude in meters

Debido a que la velocidad de bits de datos del GPS es de 50 bits por segundo y cada bloque de código C/A es de 1 milisegundo, cada bit consta de 20 bloques de código C/A. Inicialice este parámetro.

numCACodeBlocksPerBit = 20;

Inicialice las propiedades necesarias para el receptor. El receptor debe esperar algún tiempo antes de poder comenzar a recibir algunos datos significativos debido a los grandes retrasos modelados en el ejemplo. Una señal GPS típica tiene un retraso de 60 a 90 milisegundos. Entonces, si el receptor espera 100 milisegundos, comienza a procesar señales en lugar de ruido puro.

rxWaitTime = 100;       % Milliseconds
performInitSync = 1;    % Initially this must be set to true

% Initialize maximum number of tracking channels. Minimum of 4 tracking
% channels are needed for proper functioning of the GPS receiver.
maxNumTrackingChannels = 8;

% Noise bandwidth of each of the tracking loops
PLLNoiseBW = 90; % In Hz
FLLNoiseBW = 4;  % In Hz
DLLNoiseBW = 3;  % In Hz

% Bit synchronization parameters
isBitSyncComplete = zeros(maxNumTrackingChannels,1);
numBitsForBitSync = 100;
numWaitingStepsForBitSync = numCACodeBlocksPerBit*numBitsForBitSync;
rxcntr = 1;

La estimación de la posición requiere, como mínimo, que las subtramas 2 y 3 estén decodificadas. La duración de una subtrama es de 6 segundos. Para este ejemplo, decodifique 48,5 segundos de datos para garantizar que se incluyan las subtramas 2 y 3.

minTimeForPosEst = 48.5;                         % In seconds
minStepsForPosEst = minTimeForPosEst/stepTime;
subframeDuration = 6;                            % In seconds
numStepsPerSubframe = subframeDuration/stepTime;

Inicialice las constantes físicas necesarias para la simulación.

c = physconst("LightSpeed"); % Speed of light in m/sec
fe = 1575.42e6;              % GPS L1 frequency in Hz
Dt = 12;                     % Directivity of the transmit antenna in dBi
DtLin = db2pow(Dt);
Dr = 4;                      % Directivity of the receive antenna in dBi
DrLin = db2pow(Dr);
Pt = 44.8;                   % Typical transmission power of a GPS satellite in watts
k = physconst("boltzmann");  % Boltzmann constant in Joules/Kelvin
T = 300;                     % Room temperature in Kelvin
rxBW = samplingRate;         % Two-sided Bandwidth
Nr = k*T*rxBW;               % Thermal noise power in watts
seed = 73;

Configuración de simulación

En esta sección, configure el ejemplo en función de los parámetros que se inicializan en la sección Inicializar parámetros.

Configure el escenario satelital basado en el archivo RINEX.

% Initialize satellite scenario
sc = satelliteScenario;

% Set up the satellites based on the RINEX data
sat = satellite(sc,rinexdata,"OrbitPropagator","gps");
rx = groundStation(sc,rxlat,rxlon); % Set up the receiver
ac = access(sat,rx);                % Calculate access between satellites and the receiver

% Get the list of satellites that are considered in satellite scenario from
% the RINEX data
indices = ones(length(sat),1);
for isat = 1:length(sat)
    ele = orbitalElements(sat(isat));

    % Check for match of time of applicability and pseudo-random noise
    % (PRN) IDs so that data from RINEX file is considered for waveform
    % generation
    indices(isat) = find(rinexdata.GPS(:,:).Toe == ele.GPSTimeOfApplicability & ...
        rinexdata.GPS(:,:).SatelliteID == ele.PRN);
end

Generar los datos de navegación a transmitir. Primero genere un objeto de configuración de navegación utilizando la función auxiliar HelperGPSRINEX2Config. Esta función auxiliar asigna los parámetros leídos del archivo RINEX a los parámetros de configuración necesarios para la generación de datos de navegación.

% Generate navigation configuration object
navcfg = HelperGPSRINEX2Config(almanacFileName,rinexdata.GPS(indices,:));

[mintow,locmintow] = min([navcfg(:).HOWTOW]);

% The time of week (TOW) value is such that it must be 1 plus a multiple of
% 5 so that data of subframe 1 is always generated first.
mintow = ceil((mintow-1)/5)*5 + 1;

% HOWTOW is a counter that contains integer values. Incrementing by a value
% of 1 represents 6 seconds of data. The counter increases by a value of 1
% for each new subframe.
[navcfg(:).HOWTOW] = deal(mintow);
% Set the starting of frames based on mintow
firstsubframeID = mod(mintow-1,125) + 1;
frameID = ceil(firstsubframeID/5);
allFrameIDs = [frameID:25,1:(frameID-1)];
[navcfg(:).FrameIndices] = deal(allFrameIDs);

% Generate GPS navigation data
numNavBits = 37500;                         % Full GPS data length is 37500 bits
navdata = zeros(numNavBits,length(navcfg)); 
for isat = 1:length(navcfg)
    navdata(:,isat) = HelperGPSNAVDataEncode(navcfg(isat));
end

Obtenga las posiciones y velocidades de los satélites a lo largo del tiempo para modelar con precisión el Doppler, el retraso y la potencia. Para obtener más información, consulte Calculate Latency and Doppler in a Satellite Scenario (Satellite Communications Toolbox). 

sc.StartTime = HelperGPSConvertTime(navcfg(locmintow).WeekNumber, ...
    mintow*subframeDuration);
sc.StopTime = sc.StartTime + seconds(simulatedDataDuration);
sc.SampleTime = stepTime;

acstats = accessStatus(ac);
% This example runs for a maximum of 2 minutes of data and in that
% duration, access does not change. Hence, consider only the access status
% of the first sample time.
satindices = find(acstats(:,1));
numsat = length(satindices);

Esta figura muestra los satélites de la constelación GPS. Las líneas de puntos verdes indican el acceso entre el receptor y los satélites GPS. Utilice play(sc) para obtener esta cifra.

% Calculate Doppler shift over time for all the visible satellites
fShift = dopplershift(sat(satindices),rx,Frequency=fe);

% Get the states of satellites in the sky over time. Also, shuffle the
% dimensions such that each row corresponds to one satellite data
initsatpos = permute(states(sat,"CoordinateFrame","ecef"),[3 1 2]);

% Compute the distance between satellites and receiver over time
satdist = zeros(numsat,numSteps);
for istep = 1:numSteps
    satdist(:,istep) = pseudoranges([rxlat,rxlon,rxalt], ...
        initsatpos(satindices,:,istep),"RangeAccuracy",0);
end
delays = satdist/c;

% Power at receiver from free space pathloss equation
Pr = (Pt*DtLin*DrLin)*(1./(4*pi*(fe+fShift).*delays));
SNRs = 10*log10(Pr/Nr); % in dB
timeinweek = mintow*6;
PRNIDs = [navcfg(:).PRNID];
disp("Available satellites - " + num2str(PRNIDs(satindices)))
Available satellites - 5  10  13  15  18  20  23  24  27  29  30

Inicializar el objeto de generación de forma de onda según la configuración inicial.

gpswavegen = gpsWaveformGenerator(SignalType="legacy",PRNID=PRNIDs(satindices),SampleRate=samplingRate);

Configurar el objeto de canal GPS.

gnssChannelObj = HelperGNSSChannel(SampleRate=samplingRate, ...
    RandomStream="mt19937ar with seed",Seed=seed);
caCodeBlockDuration = 1e-3;                     % Constant value
SamplesPerCodeBlock = caCodeBlockDuration*samplingRate; % Number of samples per C/A code block
numcacodeblocks = stepTime/caCodeBlockDuration; % Each C/A-code block is of 1 millisecond duration

Inicialice el objeto del diagrama de constelación.

rxconstellation = comm.ConstellationDiagram(1,ShowReferenceConstellation=false, ...
    Title="Constellation diagram of signal at the output of tracking");

Inicialice el objeto que realiza la sincronización inicial en el receptor GPS.

initialsync = gnssSignalAcquirer;
initialsync.SampleRate = samplingRate;
initialsync.IntermediateFrequency = 0        % Baseband signal
initialsync = 
  gnssSignalAcquirer with properties:

              GNSSSignalType: "GPS C/A"
                  SampleRate: 10230000
       IntermediateFrequency: 0
              FrequencyRange: [-10000 10000]
         FrequencyResolution: 500
    DetectionThresholdFactor: 1.9000

Inicialice el escritor de archivos de banda base.

if WriteWaveformToFile == 1
    bbWriter = comm.BasebandFileWriter("gpsBBWaveform.bb",samplingRate,0);
end

Inicialice las propiedades necesarias para la cadena de simulación de un extremo a otro.

% Properties that store outputs of initial synchronization
[doppleroffsets,codephoffsets] = deal(zeros(1,maxNumTrackingChannels));

% Properties required for storing outputs from tracking module
[accuPh,accuFqy,accuFqyErr,accuPhErr,accuIntegWave,accuDelay,accuDelayErr] = ...
                deal(zeros(numSteps,maxNumTrackingChannels));

% Properties to store outputs of bit synchronization and frame
% synchronization
[maxTransitionLocation,sampleCounter] = deal(zeros(maxNumTrackingChannels,1));
syncidx = zeros(maxNumTrackingChannels,1);

% Property to store output of data decoder
deccfg = cell(maxNumTrackingChannels,1);

% Initialize maximum number of steps for which the simulation chain runs
maxSimSteps = 50/stepTime;

Cadena de simulación de extremo a extremo

Esta sección contiene tres pasos principales:

  1. Generación de forma de onda

  2. Canal de propagación

  3. Algoritmos de procesamiento de señales del receptor.

Como parte del canal de propagación, modele el desplazamiento, el retraso y el ruido Doppler. Como parte del receptor, los procesos involucrados incluyen sincronización inicial, seguimiento, sincronización de bits, sincronización de tramas y decodificación de datos.

tic % Start of simulation
for istep = 1:numSteps
    %% Generate waveform
    bitidx = floor((istep-1)/numCACodeBlocksPerBit)+1;
    % Generate GPS signals for one bit for all satellites indicated by satindices
    CodeBlockNumber = mod(istep-1,numCACodeBlocksPerBit);
    if CodeBlockNumber == 0
        OnebitWaveform = gpswavegen(navdata(bitidx,satindices));
    end
    % Select the signals coresponding to one code block for further processing
    iqsig = OnebitWaveform(CodeBlockNumber*SamplesPerCodeBlock + (1:SamplesPerCodeBlock),:);
    
    %% Pass the signal through a channel to introduce frequency offset, delay and noise
    gnssChannelObj.FrequencyOffset = fShift(:,istep).';
    gnssChannelObj.SignalDelay = delays(:,istep).'; 
    gnssChannelObj.SignalToNoiseRatio = SNRs(:,istep).';
    waveform = gnssChannelObj(iqsig);

    if WriteWaveformToFile == 1
        bbWriter(waveform);
    end

    timeinweek = timeinweek + stepTime;

    %% Receiver
    % Because there are large delays (70 to 80 milliseconds) modeled on the
    % signal, start the receiver after some time only so as to process
    % valid signal instead of pure noise.
    if istep > rxWaitTime
        if performInitSync == 1
            performInitSync = 0;
            [y,corrval] = initialsync(waveform,1:32); % Initial synchronization

            PRNIDsToSearch = y(y(:,4).IsDetected==1,1).PRNID.';
            doppleroffsets = y(y(:,4).IsDetected==1,2).FrequencyOffset;
            codephoffsets = y(y(:,4).IsDetected==1,3).CodePhaseOffset;

            numdetectsat = length(PRNIDsToSearch);
            if numdetectsat > maxNumTrackingChannels
                % Set this value to limit the number of tracking channels
                numdetectsat = maxNumTrackingChannels;
            end

            % Perform the required initialization of tracking modules for each
            % channel and buffers to store the data.
            disp("The detected satellite PRN IDs: " + num2str(PRNIDsToSearch))

            % Plot the correlation plot for the first satellite
            figure;
            mesh(-10e3:500:10e3, 0:size(corrval,1)-1, corrval(:,:,PRNIDsToSearch(1)));
            xlabel("Doppler Offset")
            ylabel("Code Phase Offset")
            zlabel("Correlation")
            msg = ["Correlation Plot for PRN ID: " num2str(PRNIDsToSearch(1))];
            title(msg)
            
            framesyncbuffer = cell(1,numdetectsat);

            % Create a cell array, where each element corresponds to a carrier
            % tracking object.
            carrierCodeTrack = cell(numdetectsat,1);
            framesync = cell(numdetectsat,1);
            prev_cntr = ones(numdetectsat,1);

            % Update properties for each tracking loop
            for isat = 1:numdetectsat
                carrierCodeTrack{isat} = HelperGPSCACodeCarrierTracker;
                carrierCodeTrack{isat}.SampleRate = samplingRate;
                carrierCodeTrack{isat}.CenterFrequency = 0;
                carrierCodeTrack{isat}.PLLNoiseBandwidth = PLLNoiseBW;
                carrierCodeTrack{isat}.FLLNoiseBandwidth = FLLNoiseBW;
                carrierCodeTrack{isat}.DLLNoiseBandwidth = DLLNoiseBW;
                carrierCodeTrack{isat}.PLLIntegrationTime = 1; % In milliseconds
                carrierCodeTrack{isat}.PRNID = PRNIDsToSearch(isat);
                carrierCodeTrack{isat}.InitialDopplerShift = doppleroffsets(isat);
                carrierCodeTrack{isat}.InitialCodePhaseOffset = codephoffsets(isat);

                % Initialize frame synchronization object
                framesync{isat} = HelperGPSLNAVFrameSynchronizer;
            end
        end

        % Because it would be sufficient to get receiver position after
        % running the simulation for 50 seconds of data, stop the loop
        % after executing for 50 seconds of data. In the default run of the
        % example, only 3 seconds of data is processed and this line is not
        % needed. This check is when running the example for large data
        % duration only (for at least 50 seconds of data).
        if istep > maxSimSteps
            break;
        end

        for isat = 1:numdetectsat % Perform tracking for each satellite

            [integwave,fqyerr,fqyoffset,pherr,phoffset,derr,dnco] = ...
                carrierCodeTrack{isat}(waveform);

            % Accumulate the values to see the results at the end
            accuFqyErr(rxcntr,isat) = fqyerr;
            accuFqy(rxcntr,isat) = fqyoffset;
            accuPhErr(rxcntr,isat) = pherr;
            accuPh(rxcntr,isat) = phoffset;
            accuIntegWave(rxcntr,isat) = sum(integwave);
            accuDelayErr(rxcntr,isat) = derr;
            accuDelay(rxcntr,isat) = dnco;
        end

        % Perform bit synchronization, frame synchronization, and data
        % decoding if numWaitingStepsForBitSync of receiver steps are
        % complete.
        if rxcntr > numWaitingStepsForBitSync
            % For each detected satellite, perform bit synchronization,
            % frame synchronization, and data decoding
            for isat = 1:numdetectsat
                if ~isBitSyncComplete(isat)
                    maxTransitionLocation(isat) = ...
                        gnssBitSynchronize( ...
                        imag(accuIntegWave(1:numWaitingStepsForBitSync,isat)), ...
                        numCACodeBlocksPerBit);
                    isBitSyncComplete(isat) = 1;
                    sampleCounter(isat) = rxcntr - maxTransitionLocation(isat) + 1;
                    framesyncbuffer{isat} = accuIntegWave( ...
                        maxTransitionLocation(isat):end,isat);
                else % Perform frame synchronization and data decoding
                    sampleCounter(isat) = sampleCounter(isat) + 1;
                    framesyncbuffer{isat}(sampleCounter(isat)) = accuIntegWave(rxcntr,isat);
                    if mod(sampleCounter(isat),numStepsPerSubframe) == 0
                        samples = framesyncbuffer{isat}(sampleCounter(isat) - ...
                            numStepsPerSubframe+1:sampleCounter(isat));
                        sym = mean(reshape(samples,numCACodeBlocksPerBit,[]));
                        bits = imag(sym)<0;
                        [syncidx(isat),rxsubframes,subframeIDs] = framesync{isat}(bits(:));
                        if ~isempty(rxsubframes) % Then perform data decoding
                            deccfg{isat}.PRNID = PRNIDsToSearch(isat);
                            deccfg{isat} = HelperGPSLNAVDataDecode(rxsubframes,deccfg{isat});
                            prev_cntr(isat) = rxcntr;
                        end
                    end
                end
            end
        end

        if mod(rxcntr,1000) == 0
            disp("Processed " + (rxcntr/1000) + " sec of data at the receiver.")
            rxconstellation(accuIntegWave(rxcntr-999:rxcntr,1)/ ...
                rms(accuIntegWave(rxcntr-999:rxcntr,1)))
        end

        % Update rxcntr
        rxcntr = rxcntr + 1;
    end
end
The detected satellite PRN IDs: 24   5  18  20  27  30  13  15  23  29  10

Figure contains an axes object. The axes object with title Correlation Plot for PRN ID: 24, xlabel Doppler Offset, ylabel Code Phase Offset contains an object of type surface.

Processed 1 sec of data at the receiver.
Processed 2 sec of data at the receiver.


if WriteWaveformToFile == 1
    release(bbWriter);
end
e2eTime = toc;
disp("End-to-End chain ran for " + e2eTime + " seconds.")
End-to-End chain ran for 151.7801 seconds.

Estimación de la posición del receptor GPS

Para estimar la posición del receptor, necesita conocer las distancias de al menos 4 satélites al receptor y la ubicación de cada satélite en el cielo en el momento de la transmisión. El cálculo de la posición del receptor a partir de la salida del decodificador de datos implica estos pasos.

  1. Estimar el tiempo de transmisión de los satélites visibles en un tiempo de recepción determinado: Las estimaciones del tiempo de transmisión son medidas naturales de un receptor GPS, no pseudodistancias. El tiempo de transmisión de una señal es la culminación de los cálculos de temporización realizados en la sincronización inicial, el seguimiento, la sincronización de bits, la sincronización de tramas y la decodificación de datos.

  2. Calcular pseudodistancias de todos los satélites: Debido a que el reloj del receptor no es preciso, se pueden esperar grandes desfases en los pseudorangos estimados. Por la misma razón que estas estimaciones no son rangos verdaderos, se denominan pseudorangos.

  3. Estimar la ubicación de los satélites visibles en el cielo en un momento de transmisión determinado: De los datos decodificados, se obtienen los datos de efemérides. Estos datos de efemérides son válidos durante 2 a 6 horas [2] y para calcular la posición precisa de un satélite en el cielo a partir de estos datos de efemérides en un momento determinado, utilice gnssconstellation (Navigation Toolbox).

  4. Estimar la posición del receptor: Estime la posición del receptor utilizando receiverposition (Navigation Toolbox).

Estimar los retrasos de transmisión de las señales correspondientes a los satélites visibles en un tiempo de receptor determinado.

caChipRate = 1.023e6; % In Hz
codeOffsetTime = codephoffsets(1:numdetectsat)/caChipRate;
rxcntr(rxcntr<=1) = 2; % So that rxcntr-1 is at least 1
trackingOffsetTime = accuDelay(prev_cntr(1),1:numdetectsat)/caChipRate;
bitsyncTime = (maxTransitionLocation(1:numdetectsat)-1)*caCodeBlockDuration;
framesyncTime = (syncidx(1:numdetectsat)-1)*numCACodeBlocksPerBit*caCodeBlockDuration;

% Calculate transmission delay from these parameters
delayEst = codeOffsetTime(:) - trackingOffsetTime(:) + bitsyncTime + framesyncTime;
validDelay = delayEst(syncidx~=0);

Estimar los pseudorangos.

rho = validDelay*c;


% Include time of week decoded from the received navigation message
tow = zeros(maxNumTrackingChannels,1);
for isat = 1:maxNumTrackingChannels
    if isfield(deccfg{isat},"HOWTOW")
        tow(isat) = (deccfg{isat}.HOWTOW + 1)*6; % Each subframe has a duration of 6 seconds
    end
end
tow = tow(syncidx~=0);

Estimar la posición del satélite en el momento de la segunda transmisión.

deccfg1 = deccfg(syncidx~=0);
% Estimate the satellite position only if receiver has processed at least
% minStepsForPosEst of data
isFullDataDecoded = true(length(tow),1);
if rxcntr > minStepsForPosEst
    [timeofweek,SatelliteID,Delta_n,M0,Eccentricity,sqrtA,Toe,Toc,Cis,Cic,Crs,Crc,Cus,Cuc, ...
        OMEGA0,i0,omega,OMEGA_DOT,IDOT,GPSWeek] = deal(zeros(4,1));
    reqCEIFields = ["WeekNumber","ReferenceTimeOfEphemeris","ReferenceTimeOfClock", ...
        "SemiMajorAxisLength","MeanMotionDifference","MeanAnomaly","Eccentricity", ...
        "ArgumentOfPerigee","Inclination","InclinationRate","HarmonicCorrectionTerms", ...
        "RateOfRightAscension","LongitudeOfAscendingNode"];
    for isat = 1:length(tow)
        timeofweek(isat) = tow(isat);
        isFullDataDecoded(isat) = (isfield(deccfg1{isat},"PRNID")) && ...
            (all(isfield(deccfg1{isat},reqCEIFields)));

        if isFullDataDecoded(isat)
            SatelliteID(isat) = deccfg1{isat}.PRNID;
            GPSWeek(isat) = deccfg1{isat}.WeekNumber;
            Toe(isat) = deccfg1{isat}.ReferenceTimeOfEphemeris;
            Toc(isat) = deccfg1{isat}.ReferenceTimeOfClock;
            sqrtA(isat) = sqrt(deccfg1{isat}.SemiMajorAxisLength);
            Delta_n(isat) = deccfg1{isat}.MeanMotionDifference*pi;
            M0(isat) = deccfg1{isat}.MeanAnomaly*pi;
            Eccentricity(isat) = deccfg1{isat}.Eccentricity;
            omega(isat) = deccfg1{isat}.ArgumentOfPerigee*pi;
            i0(isat) = deccfg1{isat}.Inclination*pi;
            IDOT(isat) = deccfg1{isat}.InclinationRate*pi;
            hterms = num2cell(deccfg1{isat}.HarmonicCorrectionTerms);
            [Cis(isat),Cic(isat),Crs(isat),Crc(isat),Cus(isat),Cuc(isat)] = deal(hterms{:});
            OMEGA_DOT(isat) = deccfg1{isat}.RateOfRightAscension*pi;
            OMEGA0(isat) = deccfg1{isat}.LongitudeOfAscendingNode*pi;
        end
    end

    transmissionTime = HelperGPSConvertTime(GPSWeek(isFullDataDecoded),timeofweek(isFullDataDecoded));
    refTime = HelperGPSConvertTime(GPSWeek(:),Toc(:));
    rxtimetable = timetable(refTime,SatelliteID,Delta_n,M0,Eccentricity,sqrtA, ...
        Toe,Cis,Cic,Crs,Crc,Cus,Cuc,OMEGA0,i0,omega, ...
        OMEGA_DOT,IDOT,GPSWeek);
    [satpos, satvel] = gnssconstellation(transmissionTime(1),rxtimetable);
end

Calcule la posición del receptor, lo que requiere dos entradas: pseudodistancias y la posición exacta del satélite en el momento de la transmisión, que se estima en un momento determinado del receptor. Puede calcular estos valores solo si el ejemplo se ejecuta durante al menos 48 segundos de datos. Debido a que la simulación de dichos datos llevaría mucho tiempo, los pseudorangos y las propiedades de posición del satélite (obtenidas al ejecutar el ejemplo durante 50 segundos) se almacenan en los archivos MAT. Si ejecuta la simulación durante menos de minStepsForPosEst, se cargan estos archivos MAT y se calcula la posición del receptor. De lo contrario, las propiedades requeridas se toman directamente de la cadena de simulación.

if rxcntr <= minStepsForPosEst
    % The parameters that are loaded here are valid for the default
    % configuration of this example
    load receiverPositionProperties;

    % When loading the parameters for default configuration, update the
    % isFullDataDecoded variable with all ones to compute the receiver
    % position.
    isFullDataDecoded = true(length(rho),1);

    defaultRINEXFileName = "GODS00USA_R_20211750000_01D_GN.rnx";
    defaultRxPos = [39.021, -76.827, 19];
    exampleRxPos = [rxlat, rxlon, rxalt];
    
    if ~(strcmp(rinexFileName,defaultRINEXFileName) && isequal(defaultRxPos,exampleRxPos))
        warning("satcom:EndToEndGPSLNAVReceiverExample:InsufficientData", ...
            "Estimated receiver position may be different from what you provided" + ...
            " as the simulation didn't run for entire data." + ...
            " To get accurate receiver position, run the example" + ...
            " for at least 50 seconds of navigation data.");
    end
end

% Compute the receiver position
rxposest = receiverposition(rho(isFullDataDecoded),satpos(isFullDataDecoded,:))
rxposest = 1×3

   39.0209  -76.8271   21.9329


estRxPosNED = lla2ned(rxposest,[rxlat,rxlon,rxalt],'ellipsoid');
distanceError = vecnorm(estRxPosNED) % In meters
distanceError = 
17.7598

Exploración adicional

Este ejemplo muestra cómo realizar el procesamiento del receptor GPS para la constelación de satélites simulada durante solo 3 segundos. Amplíe el ejemplo a 50 segundos y calcule la posición del receptor para ver cómo funciona el receptor GPS.

Utilice su propio archivo RINEX para configurar el ejemplo y estimar la posición del receptor a partir de la forma de onda generada en este ejemplo.

Apéndice

Este ejemplo utiliza estos datos y archivos auxiliares:

Referencias

[1] IS-GPS-200, Rev: M. Interfaces de usuario de navegación/segmento espacial GPS NAVSTAR. 21 de mayo de 2021; Código de identificación: 66RP1.

[2] Elliott D. Kaplan y C. Hegarty, eds., Entendiendo GPS/GNSS: Principios y aplicaciones , Tercera edición, Serie de tecnología y aplicaciones GNSS (Boston; Londres: Casa Artech, 2017).