GNSS接收机设计：(5)导航解算

姜鱼yuu

已于 2024-01-10 20:15:07 修改

阅读量2.4k

点赞数 41

分类专栏： GNSS接收机文章标签：笔记 matlab

于 2023-12-29 21:57:54 首次发布

本文链接：https://blog.csdn.net/Kingkingkingdom/article/details/135137900

版权

GNSS接收机专栏收录该内容

6 篇文章 30 订阅

订阅专栏

postNavigation.m

本文目录

整体流程图
（1）首先检查数据是否足够
（2）找到导航数据帧的起始位置
（3）提取导航数据
（4）解码星历
（5）验证星历数据完整性
（6）检查剩余卫星数量
（7）初始化
（8）开始解算

整体流程图

简单画了个流程图：
在这里插入图片描述

（1）首先检查数据是否足够

根据GNSS信号结构，我们知道，1个帧是由5个子帧构成，示意图如下：

或者看下面这个图更直观。

总之，前三个子帧包含的信息是最重要的，而sub4和5存放的主要是全部卫星的历书（也就是不精确的星历），原则上不需要4、5也可以。譬如接收机冷启动一般需要十几秒，也就是获取前三个子帧所需要的时间（1个子帧是6s）。
但是在我们这套软件接收机中，星历解码函数ephemeris需要输入全部的5个子帧，而5个子帧=30s=30000ms。并且如果跟踪在子帧中间开始，还要额外再增加几秒钟（为什么？），这样的话，最终我们是需要36000ms的数据。
因此第一部分代码如下：

function [navSolutions, eph] = postNavigation(trackResults, settings)

%% 检查是否有足够的数据 ===========

settings = initSettings();
if (settings.msToProcess < 36000) || (sum([trackResults.status] ~= '-') < 4)
    % 如果数据不足3600ms或者卫星不足4颗，都无法继续进行
    disp('Record is to short or too few satellites tracked. Exiting!');
    navSolutions = [];
    eph          = [];
    return
end

（2）找到导航数据帧的起始位置

主要是调用了findPreambles函数，可以看我之前的文章帧同步函数

[subFrameStart, activeChnList] = findPreambles(trackResults, settings);

（3）提取导航数据

根据ephemeris函数的要求，输入需要5个子帧，并且是二进制形式，所以据此对导航数据进行处理

for channelNr = activeChnList

    % 在导航数据I_P中，取5子帧长度的数据；5个子帧组成一个帧，所以相当于取了1帧
    % 1子帧300比特，5子帧就是1500比特，1比特=20个数据单元
    % subFrameStart(channelNr)相当于数组索引，来指示从 导航数据I_P 中取信息的开始位置
    navBitsSamples = trackResults(channelNr).I_P(subFrameStart(channelNr) - 20 : ...
                               subFrameStart(channelNr) + (1500 * 20) -1)';
    
    % 把取出的数据重新形状为20行矩阵（1比特=20个数据单元）再求和
    % 所以navBits存放的是以比特为单位的数据
    navBitsSamples = reshape(navBitsSamples, ...
                             20, (size(navBitsSamples, 1) / 20));
    navBits = sum(navBitsSamples);

    %阈值为0，判决，得到01序列
    navBits = (navBits > 0);
    
    % dec2bin：该函数作用是十进制转换成二进制**字符串**，得到navBitsBin.
    % 这里navBits已经是二进制了，所以这个转换其实只为了把数值变成字符串
    navBitsBin = dec2bin(navBits);

（4）解码星历

上一部分的最后，我们说使用dec2bin函数的目的是类型转换，究其原因，还是因为ephemeris这个函数要求输入是字符串类型，这里就不展开说了。总之是得到星历eph、周内时TOW（应该都是十进制）

 [eph(trackResults(channelNr).PRN), TOW] = ...
                            ephemeris(navBitsBin(2:1501)', navBitsBin(1));

（5）验证星历数据完整性

如果当前这颗卫星的星历数据不完整，那么就把它扔出activeChnList，然后继续处理下一颗星。

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

（6）检查剩余卫星数量

注意，上面需要遍历可用列表里的所有卫星，解码它们的星历，全部处理完后再进行这一步：

%% 处理完不可用的卫星之后，检查剩下的卫星是否>3颗 =====================
if (isempty(activeChnList) || (size(activeChnList, 2) < 4))
    % Show error message and exit
    disp('Too few satellites with ephemeris data for postion calculations. Exiting!');
    navSolutions = [];
    eph          = [];
    return

end

（7）初始化

% 将卫星仰角阵列设置为 INF，以便首次计算接收机位置时包含所有卫星。
% 由于此时没有接收机位置估计值，因此没有参考点来查找仰角。
satElev  = inf(1, settings.numberOfChannels);

% 保存活动频道列表。该列表包含已跟踪并拥有所需星历数据的卫星。
% 下一步，列表将取决于每颗卫星的仰角，而仰角会随时间变化。
readyChnList = activeChnList;

transmitTime = TOW;  %TOW 周内时

（8）开始解算

首先进行一些准备工作

for currMeasNr = 1:fix((settings.msToProcess - max(subFrameStart)) / ...
                                                     settings.navSolPeriod)

    % intersect计算交集，把仰角小于阈值的卫星也扔出去
    activeChnList = intersect(find(satElev >= settings.elevationMask), ...
                              readyChnList);
    
    % Save list of satellites used for position calculation
    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);

调用calculatePseudoranges计算伪距

    navSolutions.channel.rawP(:, currMeasNr) = calculatePseudoranges(...
            trackResults, ...
            subFrameStart + settings.navSolPeriod * (currMeasNr-1), ...
            activeChnList, settings);

这个函数较简单，这里就不展开解释了，只提一下以前不理解的一个地方。
在calculatePseudoranges代码最后，用了下面这两个式子来计算传播时间：

minimum         = floor(min(travelTime));
travelTime      = travelTime - minimum + settings.startOffset;

解释：
因为缺少时间参考点，所以我们只能采用相对的方法来计算传播时间。
首先假设前提条件：所有卫星的子帧1都在同一时刻发射。因此就可以用接收时间的差值来表示接收机到不同卫星的伪距的差值（相对伪距）。而第一子帧每30s出现一次，两颗卫星之间最大的时间差是19ms，这就确保子帧1之间是可比的，也就是说很容易将不同卫星第一次和第二次发射的子帧1区分开来。卫星到用户接收机的时间延迟范围在67-86ms，所以在代码中取：
settings.startOffset = 68.802; % 单位ms
把它作为初始传播时间加在travelTime – minimum后面，就把相对时间转化为绝对时间，从而进行伪距的计算。

调用satpos解算卫星位置

根据星历，先计算各种轨道参数，再经过两次旋转得到卫星在地心地固坐标系下的位置。

 [satPositions, satClkCorr] = satpos(transmitTime, ...
                                        [trackResults(activeChnList).PRN], ...
                                        eph, settings);

最小二乘解算接收机位置

详见位置解算原理及leastSquarePos.m注释

 if size(activeChnList, 2) > 3

        %=== 调用leastSquarePos计算位置 ==================================
        [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);

 %--- 保存结果 -------------------------------------------------
        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';

用钟差修正伪距

  %==============
        navSolutions.channel.correctedP(activeChnList, currMeasNr) = ...
                navSolutions.channel.rawP(activeChnList, currMeasNr) + ...
                satClkCorr' * settings.c + navSolutions.dt(currMeasNr);

坐标转换

        %=== 将接收机位置转换为大地坐标（经纬高） ==============================
        [navSolutions.latitude(currMeasNr), ...
         navSolutions.longitude(currMeasNr), ...
         navSolutions.height(currMeasNr)] = cart2geo(...
                                            navSolutions.X(currMeasNr), ...
                                            navSolutions.Y(currMeasNr), ...
                                            navSolutions.Z(currMeasNr), ...
                                            5);

        %=== 将接收机位置转换为UTM坐标 =============================
        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);

一些赋值

%将缺失解设为 NaN。这些结果将在所有绘图中自动排除。
        % 对于 DOP 而言，使用 0 更为简单。
        % 在某些进一步处理中，可能需要从结果中排除 NaN 值，以获得正确的结果。
        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));

       

    end % if size(activeChnList, 2) > 3

更新传输时间

    %=======================
    transmitTime = transmitTime + settings.navSolPeriod / 1000;

end %for currMeasNr...