GPS接收机设计（5）—

1.流程图

（待补充，没人看就不补充了。。。。。。）

2.代码解析

function [navSolutions, eph] = postNavigation(trackResults, settings)
%Function calculates navigation solutions for the receiver (pseudoranges,
%positions). At the end it converts coordinates from the WGS84 system to
%the UTM, geocentric or any additional coordinate system.
%
%[navSolutions, eph] = postNavigation(trackResults, settings)
%
%   Inputs:
%       trackResults    - results from the tracking function (structure
%                       array).
%       settings        - receiver settings.
%   Outputs:
%       navSolutions    - contains measured pseudoranges, receiver
%                       clock error, receiver coordinates in several
%                       coordinate systems (at least ECEF and UTM).
%       eph             - received ephemerides of all SV (structure array).%--------------------------------------------------------------------------
%                           SoftGNSS v3.0settings = initSettings();
if (settings.msToProcess < 36000) || (sum([trackResults.status] ~= '-') < 4)% Show the error message and exitdisp('Record is to short or too few satellites tracked. Exiting!');navSolutions = [];eph          = [];return
end%% 寻找帧起始位置==========================================[subFrameStart, activeChnList] = findPreambles(trackResults, settings);%%解码星历 =====================================================for channelNr = activeChnList%===  将跟踪输出的结果转换为导航电文=======================%--- 从跟踪输出结果复制5个子帧数据---------------navBitsSamples = trackResults(channelNr).I_P(subFrameStart(channelNr) - 20 : ...subFrameStart(channelNr) + (1500 * 20) -1)';navBitsSamples = reshape(navBitsSamples, ...20, (size(navBitsSamples, 1) / 20));%每20个值合成1列%--- Sum all samples in the bits to get the best estimate -------------navBits = sum(navBitsSamples);%对列求和navBits = (navBits > 0);%转换为逻辑值0/1navBitsBin = dec2bin(navBits);%以字符矢量形式返回navBits的二进制表示形式%=== 解调星历表和截短周内时计数（TOW） ================[eph(trackResults(channelNr).PRN), TOW] = ...ephemeris(navBitsBin(2:1501)', navBitsBin(1));%--- 剔除缺少必要导航数据信息的卫星 -----if (isempty(eph(trackResults(channelNr).PRN).IODC) || ...isempty(eph(trackResults(channelNr).PRN).IODE_sf2) || ...isempty(eph(trackResults(channelNr).PRN).IODE_sf3))activeChnList = setdiff(activeChnList, channelNr);end
end%% 检查卫星数量是否超过3=====================
if (isempty(activeChnList) || (size(activeChnList, 2) < 4))% Show error message and exitdisp('Too few satellites with ephemeris data for postion calculations. Exiting!');navSolutions = [];eph          = [];return
end%% Initialization初始化 =========================================================% Set the satellite elevations array to INF to include all satellites for
% the first calculation of receiver position. There is no reference point
% to find the elevation angle as there is no receiver position estimate at
% this point.将卫星高度角设置为无穷大以包含所有卫星用于计算接收机的位置；在没有参考点的情况下计算高度角
satElev  = inf(1, settings.numberOfChannels);readyChnList = activeChnList;%保存已跟踪和具有必要导航信息的卫星transmitTime = TOW;%截短周内计时%##########################################################################
%#   卫星与接收机位置解算                  #
%##########################################################################%% 初始化当前测量 ==================================
for currMeasNr = 1:fix((settings.msToProcess - max(subFrameStart)) / ...settings.navSolPeriod)%每50ms解算一次位置% Exclude satellites, that are belove elevation mask activeChnList = intersect(find(satElev >= settings.elevationMask), ...readyChnList);%选出高度角大于settings.elevationMask的卫星%保存用于位置解算的卫星列表navSolutions.channel.PRN(activeChnList, currMeasNr) = ...[trackResults(activeChnList).PRN]; % These two lines help the skyPlot function. The satellites excluded% do to elevation mask will not "jump" to possition (0,0) in the sky% plot.navSolutions.channel.el(:, currMeasNr) = ...NaN(settings.numberOfChannels, 1);navSolutions.channel.az(:, currMeasNr) = ...NaN(settings.numberOfChannels, 1);%% 伪距计算 ======================================================navSolutions.channel.rawP(:, currMeasNr) = calculatePseudoranges(...trackResults, ...subFrameStart + settings.navSolPeriod * (currMeasNr-1), ...activeChnList, settings);%% 卫星位置和钟差信息 =======================[satPositions, satClkCorr] = satpos(transmitTime, ...[trackResults(activeChnList).PRN], ...eph, settings);%% 确定接收机位置=================================================% 3D receiver position can be found only if signals from more than 3% satellites are available  if size(activeChnList, 2) > 3 %确定卫星数是否大于3颗？%=== 解算接收机位置==================================[xyzdt, ...navSolutions.channel.el(activeChnList, currMeasNr), ...navSolutions.channel.az(activeChnList, currMeasNr), ...navSolutions.DOP(:, currMeasNr)] = ...leastSquarePos(satPositions, ...navSolutions.channel.rawP(activeChnList, currMeasNr)' + satClkCorr * settings.c, ...settings);%最小二乘法%--- Save results -------------------------------------------------navSolutions.X(currMeasNr)  = xyzdt(1);navSolutions.Y(currMeasNr)  = xyzdt(2);navSolutions.Z(currMeasNr)  = xyzdt(3);navSolutions.dt(currMeasNr) = xyzdt(4);%接收机钟差% 更新卫星高度角向量satElev = navSolutions.channel.el(:, currMeasNr);satElev = satElev';%chen 20200215***************************************%=== 用星钟误差校正伪距测量Correct pseudorange measurements for clocks errors ===========navSolutions.channel.correctedP(activeChnList, currMeasNr) = ...navSolutions.channel.rawP(activeChnList, currMeasNr) + ...satClkCorr' * settings.c + navSolutions.dt(currMeasNr);%% 坐标转换 ==================================================%=== Convert to geodetic coordinates转换到大地坐标系 ==============================[navSolutions.latitude(currMeasNr), ...navSolutions.longitude(currMeasNr), ...navSolutions.height(currMeasNr)] = cart2geo(...navSolutions.X(currMeasNr), ...navSolutions.Y(currMeasNr), ...navSolutions.Z(currMeasNr), ...5);%=== Convert to UTM coordinate system =============================navSolutions.utmZone = findUtmZone(navSolutions.latitude(currMeasNr), ...navSolutions.longitude(currMeasNr));[navSolutions.E(currMeasNr), ...navSolutions.N(currMeasNr), ...navSolutions.U(currMeasNr)] = cart2utm(xyzdt(1), xyzdt(2), ...xyzdt(3), ...navSolutions.utmZone);else % if size(activeChnList, 2) < 3 %--- There are not enough satellites to find 3D position ----------disp(['   Measurement No. ', num2str(currMeasNr), ...': Not enough information for position solution.']);%--- Set the missing solutions to NaN. These results will be%excluded automatically in all plots. For DOP it is easier to use%zeros. NaN values might need to be excluded from results in some%of further processing to obtain correct results.navSolutions.X(currMeasNr)           = NaN;navSolutions.Y(currMeasNr)           = NaN;navSolutions.Z(currMeasNr)           = NaN;navSolutions.dt(currMeasNr)          = NaN;navSolutions.DOP(:, currMeasNr)      = zeros(5, 1);navSolutions.latitude(currMeasNr)    = NaN;navSolutions.longitude(currMeasNr)   = NaN;navSolutions.height(currMeasNr)      = NaN;navSolutions.E(currMeasNr)           = NaN;navSolutions.N(currMeasNr)           = NaN;navSolutions.U(currMeasNr)           = NaN;navSolutions.channel.az(activeChnList, currMeasNr) = ...NaN(1, length(activeChnList));navSolutions.channel.el(activeChnList, currMeasNr) = ...NaN(1, length(activeChnList));% TODO: Know issue. Satellite positions are not updated if the% satellites are excluded do to elevation mask. Therefore rasing% satellites will be not included even if they will be above% elevation mask at some point. This would be a good place to% update positions of the excluded satellites.end % if size(activeChnList, 2) > 3%=== Update the transmit time ("measurement time") ====================transmitTime = transmitTime + settings.navSolPeriod / 1000;%下一次解算时间end %for currMeasNr...

3.子程序

3.1 伪距计算

function [pseudoranges] = calculatePseudoranges(trackResults, ...msOfTheSignal, ...channelList, settings)%   Inputs:
%       trackResults    - output from the tracking function
%       msOfTheSignal   - pseudorange measurement point (millisecond) in
%                       the trackResults structure伪距测量起点对应的
%       channelList     - list of channels to be processed
%       settings        - receiver settings
%
%   Outputs:
%       pseudoranges    - relative pseudoranges to the satellites. travelTime = inf(1, settings.numberOfChannels);% 每个扩频码所对应的采样点个数
samplesPerCode = round(settings.samplingFreq / ...(settings.codeFreqBasis / settings.codeLength));clear channelNr;
for channelNr = channelList%--- Compute the travel times -----------------------------------------    travelTime(channelNr) = ...trackResults(channelNr).absoluteSample(msOfTheSignal(channelNr)) / samplesPerCode;%传输的毫秒数
end%--- 截断伪距测量（与实际传输时间差一个常数）---------------------
minimum         = floor(min(travelTime));
travelTime      = travelTime - minimum + settings.startOffset;
%所有通道的传播时间差一个常数不影响距离计算，最终将体现在接收机钟差中
pseudoranges    = travelTime * (settings.c / 1000);

3.2 最小二乘法解算位置

function [pos, el, az, dop] = leastSquarePos(satpos, obs, settings)
%Function calculates the Least Square Solution.
%
%[pos, el, az, dop] = leastSquarePos(satpos, obs, settings);
%
%   Inputs:
%       satpos      - Satellites positions (in ECEF system: [X; Y; Z;] -
%                   one column per satellite)
%       obs         - Observations - the pseudorange measurements to each
%                   satellite:
%                   (e.g. [20000000 21000000 .... .... .... .... ....])
%       settings    - receiver settings
%
%   Outputs:
%       pos         - receiver position and receiver clock error
%                   (in ECEF system: [X, Y, Z, dt])
%       el          - Satellites elevation angles (degrees)高度角
%       az          - Satellites azimuth angles (degrees)方位角
%       dop         - Dilutions Of Precision ([GDOP PDOP HDOP VDOP TDOP])%--------------------------------------------------------------------------
%                           SoftGNSS v3.0
%--------------------------------------------------------------------------
%Based on Kai Borre
%Copyright (c) by Kai Borre
%Updated by Darius Plausinaitis, Peter Rinder and Nicolaj Bertelsen
%
% CVS record:
% $Id: leastSquarePos.m,v 1.1.2.12 2006/08/22 13:45:59 dpl Exp $
%==========================================================================%=== Initialization =======================================================
nmbOfIterations = 7;%迭代次数dtr     = pi/180;
pos     = zeros(4, 1);%初始化接收机位置和钟差
X       = satpos;%卫星位置
nmbOfSatellites = size(satpos, 2);%卫星数量A       = zeros(nmbOfSatellites, 4);%初始化系数矩阵
omc     = zeros(nmbOfSatellites, 1);
az      = zeros(1, nmbOfSatellites);%初始化方位角
el      = az;%初始化高度角%=== Iteratively find receiver position ===================================
for iter = 1:nmbOfIterationsfor i = 1:nmbOfSatellitesif iter == 1%--- Initialize variables at the first iteration --------------Rot_X = X(:, i);trop = 2;else%--- Update equations -----------------------------------------rho2 = (X(1, i) - pos(1))^2 + (X(2, i) - pos(2))^2 + ...(X(3, i) - pos(3))^2;%平方距离traveltime = sqrt(rho2) / settings.c ;%传输时间%--- Correct satellite position (do to earth rotation) --------Rot_X = e_r_corr(traveltime, X(:, i));%--- Find the elevation angel of the satellite ----------------[az(i), el(i), dist] = topocent(pos(1:3, :), Rot_X - pos(1:3, :));if (settings.useTropCorr == 1)%settings.useTropCorr == 1进行电离层校正%--- Calculate tropospheric correction --------------------trop = tropo(sin(el(i) * dtr), ...0.0, 1013.0, 293.0, 50.0, 0.0, 0.0, 0.0);else% Do not calculate or apply the tropospheric correctionstrop = 0;endend % if iter == 1 ... ... else %--- Apply the corrections ----------------------------------------omc(i) = (obs(i) - norm(Rot_X - pos(1:3), 'fro') - pos(4) - trop);%--- Construct the A matrix ---------------------------------------A(i, :) =  [ (-(Rot_X(1) - pos(1))) / obs(i) ...(-(Rot_X(2) - pos(2))) / obs(i) ...(-(Rot_X(3) - pos(3))) / obs(i) ...1 ];end % for i = 1:nmbOfSatellites% These lines allow the code to exit gracefully in case of any errorsif rank(A) ~= 4pos     = zeros(1, 4);returnend%--- Find position update ---------------------------------------------x   = A \ omc;%--- Apply position update --------------------------------------------pos = pos + x;end % for iter = 1:nmbOfIterationspos = pos';%=== Calculate Dilution Of Precision ======================================
if nargout  == 4%--- Initialize output ------------------------------------------------dop     = zeros(1, 5);%--- Calculate DOP ----------------------------------------------------Q       = inv(A'*A);dop(1)  = sqrt(trace(Q));                       % GDOP    dop(2)  = sqrt(Q(1,1) + Q(2,2) + Q(3,3));       % PDOPdop(3)  = sqrt(Q(1,1) + Q(2,2));                % HDOPdop(4)  = sqrt(Q(3,3));                         % VDOPdop(5)  = sqrt(Q(4,4));                         % TDOP
end

GPS接收机设计（5）——定位解算相关推荐

4.6 定位解算和1PPS时标支持
\qquad为了解算用户位置,需要使用统一的时间参考量,合理的选择是接收时间.粗略地说,就是应用多颗卫星同时到达用户接收机的信号,求出各颗卫星发射信号时的位置和距离(伪距),用测量学上的后方交会法,解 ...
GPS接收机-从射频信号到定位解算
GNSS接收机-从特高频信号到定位解算 GPS信号历险记天线极化方式抗干扰处理射频前端基带处理信号搜索和捕获跟踪载波环锁相环锁频环码环(延迟锁定环路DLL) 位同步帧同步定位 ...
GPS接收机总体设计——数据写入、捕获、跟踪、定位解算
1.总体流程图 2. 代码实现 disp ('Starting processing...'); settings=initSettings();%系统初始化 [fid, message] = fop ...
北斗信号服务器解算,北斗导航系统接收机定位解算设计与实现
摘要: 随着北斗导航系统的建设不断推进,其应用范围越来越广,因此北斗接收机需求也越来越大.不同的应用场景的接收机结构和侧重点有所不同,但是其中的定位解算模块都是其关键部分.本文主要对北斗接收机的整体结 ...
北斗信号服务器解算,GPS/北斗定位解算算法的研究
摘要: 卫星导航是一种通过全球卫星导航系统(Global Navigation Satellite System,GNSS)精确的测定地球上任何一点的位置和时间的方法.目前,卫星导航接收机可提供个人定 ...
GNSS说第(七)讲---自适应动态导航定位(三)---序贯导航定位解算原理
GNSS说第(七)讲-自适应动态导航定位(三)-序贯导航定位解算原理序贯导航定位解算原理序贯平差法属于逐步平差法,即递推平差.其基本思想是:在对线性模型的统计性质作某些合适的假设后,基于不同的平差 ...
UWB TDOA一维定位解算
在某些定位场景,比如在隧道.走廊等区域,需要用到一维解算,下面介绍TDOA的长直线解算定位标签位置(当然也可以用TWR实现一维解算).定位模型与已知量如下: 解算不考虑z坐标,基站A和基站B的坐标分别 ...
本科生如何入门GNSS算法（二）- rtklib定位解算过程中的GNSS数据格式以及基本概念
目录 rtklib单点定位命令分析 rtklib日志定位结果pos文件说明以及定位精度评估 rtklib界面API rtkplot使用坐标转换 XYZ->BLH 其他的GNSS数据下载公 ...
GPS接收机设计（4）——帧同步
帧同步的过程比较简单,第一步将已知的同步码与跟踪截断获得的导航电文进行互相关,相关值最大时即子帧的起始位置:第二步是奇偶校验过程,获取当前子帧的前两个字与上一个子帧的最后两个数据比特进行奇偶校验,校验 ...

GPS接收机设计（5）——定位解算

1.流程图

2.代码解析

3.子程序

3.1 伪距计算

3.2 最小二乘法解算位置

GPS接收机设计（5）——定位解算相关推荐

最新文章

热门文章