基于导航观测的GPS接收机位置计算（Matlab代码实现）

 👨‍🎓个人主页：研学社的博客    

💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现

---

💥1 概述

文献来源：

全球定位系统(GPS)是第一个也是唯一一个提供全球覆盖的自主地理空间定位的全功能卫星导航系统。

为了计算用户在地球上的位置，必须确定精确的GPS卫星轨道。

一种适用的方法是利用广播星历计算卫星轨道，其优点是可以在任何给定时间确定轨道信息。然而，这种广播的星历表只能提供大约1.6米的有限精度。另一种方法是直接从国际GNSS服务(IGS)获取卫星轨道信息，也称为精确轨道或精确星历，其精度可达到毫米级。然而，精确的轨道每隔15分钟报告一次，因此需要在定位时间内插入。

本文介绍了GPS卫星轨道确定和插值的两种方法。

分析了广播星历和精密星历的优缺点。最后，利用各方法获得的卫星轨道信息，对GPS定位的性能进行了比较。

GPS接收机位置的初步猜测是从观测文件中获得的，并使用广播轨道，伪距和载波相位测量进行更新。这里从导航文件中计算广播轨道，然后计算用户和卫星之间的粗略距离。通过减去计算的距离和观察到的距离，为最小二乘技术构建残差矩阵。

📚2 运行结果

部分代码：

clc
clear
format long g

% read observation file to get orbit elements
[obs,rec_xyz] = read_rinex_obs('madr2000.06o');

% read navigation file to get orbit elements
ephemeris = read_rinex_nav('brdc2000.06n');

epochs = unique(obs.data(:, obs.col.TOW));
TimeSpan=epochs(1:2880);

% Broadcast Orbit
satOrbits.XS=zeros(1,length(TimeSpan));
satOrbits.YS=zeros(1,length(TimeSpan));
satOrbits.ZS=zeros(1,length(TimeSpan));
satOrbits.VXS=zeros(1,length(TimeSpan));
satOrbits.VYS=zeros(1,length(TimeSpan));
satOrbits.VZS=zeros(1,length(TimeSpan));
satOrbits.clk=zeros(1,length(TimeSpan));
satOrbits.Rel=zeros(1,length(TimeSpan));

% GPS Satellite Measurements
c = 2.99792458e8 ; % speed of light (m/s)
fL1 = 1575.42e6;   % L1 frequency (Hz)
fL2 = 1227.6e6;    % L2 frequency (Hz)
B=fL2^2/(fL2^2-fL1^2);
A=-B+1;
satOrbits.C1=zeros(1,length(TimeSpan));
satOrbits.L1=zeros(1,length(TimeSpan));
satOrbits.P2=zeros(1,length(TimeSpan));
satOrbits.L2=zeros(1,length(TimeSpan));
satOrbits.P3=zeros(1,length(TimeSpan)); % Iono free pseudorange
satOrbits.CorrP1=zeros(1,length(TimeSpan)); % Corrected Pseudorange from broadcast orbit
satOrbits.CorrP2=zeros(1,length(TimeSpan)); % Corrected Pseudorange from precise orbit
satOrbits.TOW=TimeSpan';
satOrbits.PRN=0;

satOrbits = repmat(satOrbits,1,32);
for ii=1:32
    satOrbits(ii).PRN=ii;
end

% Initialize User Position
userPos=zeros(length(TimeSpan),4);

for ii=1:length(TimeSpan)    
    this_TOW = TimeSpan(ii);
    index = find(obs.data(:,obs.col.TOW) == this_TOW);
    curr_obs.data = obs.data(index, :);
    curr_obs.col = obs.col;
    
    for jj=1:size(curr_obs.data,1)        
        PRN_obs.data = curr_obs.data(jj,:);
        PRN_obs.col = curr_obs.col;
        
        % Record Measurements
        satOrbits(PRN_obs.data(PRN_obs.col.PRN)).C1(ii)=PRN_obs.data(PRN_obs.col.C1);
        satOrbits(PRN_obs.data(PRN_obs.col.PRN)).L1(ii)=PRN_obs.data(PRN_obs.col.L1);
        satOrbits(PRN_obs.data(PRN_obs.col.PRN)).P2(ii)=PRN_obs.data(PRN_obs.col.P2);
        satOrbits(PRN_obs.data(PRN_obs.col.PRN)).L2(ii)=PRN_obs.data(PRN_obs.col.L2);
        
        % Calculate Iono Free Measurement
        P1 = satOrbits(PRN_obs.data(PRN_obs.col.PRN)).C1(ii);
        P2 = satOrbits(PRN_obs.data(PRN_obs.col.PRN)).P2(ii);
        P3=A*P1+B*P2;
        satOrbits(PRN_obs.data(PRN_obs.col.PRN)).P3(ii)=P3;
    end
    
    stop = 10;
    while stop ~= 1
        for jj=1:size(curr_obs.data,1)
            PRN_obs.data = curr_obs.data(jj,:);
            PRN_obs.col = curr_obs.col;
            
            % Obtain the broadcast orbits
            PRN_obs = get_broadcast_orbits(PRN_obs,ephemeris,rec_xyz');
            satOrbits(PRN_obs.data(PRN_obs.col.PRN)).XS(ii)=PRN_obs.data(PRN_obs.col.XS);
            satOrbits(PRN_obs.data(PRN_obs.col.PRN)).YS(ii)=PRN_obs.data(PRN_obs.col.YS);
            satOrbits(PRN_obs.data(PRN_obs.col.PRN)).ZS(ii)=PRN_obs.data(PRN_obs.col.ZS);
            satOrbits(PRN_obs.data(PRN_obs.col.PRN)).VXS(ii)=PRN_obs.data(PRN_obs.col.VXS);
            satOrbits(PRN_obs.data(PRN_obs.col.PRN)).VYS(ii)=PRN_obs.data(PRN_obs.col.VYS);
            satOrbits(PRN_obs.data(PRN_obs.col.PRN)).VZS(ii)=PRN_obs.data(PRN_obs.col.VZS);
            satOrbits(PRN_obs.data(PRN_obs.col.PRN)).clk(ii)=PRN_obs.data(PRN_obs.col.satClkCorr);
            satOrbits(PRN_obs.data(PRN_obs.col.PRN)).Rel(ii)=PRN_obs.data(PRN_obs.col.Rel);            
            
            % Calculate corrected pseudorange based on broadcast orbit
            satOrbits(PRN_obs.data(PRN_obs.col.PRN)).CorrP1(ii)=...
                satOrbits(PRN_obs.data(PRN_obs.col.PRN)).P3(ii)+...
                satOrbits(PRN_obs.data(PRN_obs.col.PRN)).clk(ii)+satOrbits(PRN_obs.data(PRN_obs.col.PRN)).Rel(ii);            
        end
        
        % Calculate User Position
        [broadcast_obs,~]=createObs(this_TOW,satOrbits);
        delta_xyz = comp_pos(broadcast_obs,rec_xyz');
        rec_xyz = rec_xyz + delta_xyz(1:3);
        
        stop=stop-1;
    end
    userPos(ii,1:4) = [rec_xyz; delta_xyz(4)]';
    [lon1(ii),lat1(ii),alt1(ii)] = Geodetic(rec_xyz);
end

fprintf('mean of user position:');
fprintf('%16f%16f%16f\n',mean(userPos(:,1)),mean(userPos(:,2)),mean(userPos(:,3)));
fprintf('delta_xyz:');
fprintf('%16f\n',mean(userPos(:,4)));

🎉3 参考文献

部分理论来源于网络，如有侵权请联系删除。

🌈4 Matlab代码实现