CN102176033B - 一种基于通用图形处理器的卫星导航系统比特压缩跟踪方法 - Google Patents

Abstract

一种基于通用图形处理器的卫星导航系统比特压缩跟踪方法：根据卫星导航系统跟踪通道数量，将每个通道的跟踪任务划分到通用图形处理器的区块进行处理，同时根据系统采样率N和通用处理器处理字长L，将每L个采样点数进行比特压缩并划分到一个通用处理器的线程进行处理；对于码相位跳变点位置的L个和最后不足L个的采样点分别采用一个线程单独处理；所有线程均负责完成跟踪过程中的相关运算，过程包括数据位对齐、屏蔽掩码的计算、分段相关和矢量缩减求和四个步骤；每个区块存在一个线程根据求和结果进行跟踪环路的计算，负责调整扩频码相位、载波频率和载波相位。本发明利用通用图形处理器的并行计算能力，使卫星导航系统的跟踪运算变得高效。'

Description

一种基于通用图形处理器的卫星导航系统比特压缩跟踪方法

技术领域

本发明涉及软件卫星导航接收机跟踪方法，尤其是涉及一种在通用图形处理器上实现的卫星导航系统信号跟踪的方法。

背景技术

卫星导航技术的飞速发展已逐渐取代了无线电导航、天文导航等传统导航技术，而成为一种普遍采用的导航定位技术，并在精度、实时性、全天候等方面取得了长足进步。目前，世界上多个国家和地区已经建立或者正在建设全球卫星定位导航系统，具体包括美国建立的GPS系统、俄罗斯建立的GLONASS(格鲁纳斯)系统、欧盟正在建设的“伽利略”系统和中国正在建设的北斗二代系统。因而，卫星导航系统的接收设备可能需要完成对不同卫星导航系统信号的接收和处理。

在卫星导航接收机中，信号跟踪方法是卫星导航接收机最核心的算法之一。利用通用处理器，采用软件方法实现的卫星导航接收机被称为软件接收机。传统的卫星导航软件接收机通常使用CPU、DSP和ARM作为其核心处理器。其跟踪方法如图1所示，(1)获取上次记录的载波频率和码相位信息；(2)读取对应载波频率和载波相位的载波；(3)读取码相位对应的扩频码；(4)采用相关积分方法，计算不同延迟扩频码的I、Q支路的相关累加数值；(5)使用上述数值进行环路滤波，计算新的载波频率、载波相位和码相位，并记录相关运算结果用于导航电文解算和伪距测量；(6)运行其它任务调度(如卫星捕获和定位解算)；(7)等待下一个数据块的到达。

由于卫星导航系统信号跟踪过程中的相关累加运算的计算量较大，因此很多研究提出了一些方法来对相关运算进行加速。针对基于中央处理器的软件接收机通常使用比特压缩算法(参见B.M.Ledvina，S.P.Powell，P.M.Kintner and M.L.Psiaki，“A12-Channel Real-Time GPS L1 Software Receiver”，ION-GPS 2003和Jiangping Deng，Ruizhi Chen and Jianyu Wang，“An enhanced bit-wise parallel algorithm forreal-time GPS software receiver”，GPS Solution 2010 vol.14)、MMX或者SSE指令集(参见Charkhandeh et al.“Performance Testing of a Real-Time Software-BasedGPS Receiver for x86 Processors”，ION-GNSS 2006和Gregory W.Heckler and JamesL.Garrison，“SIMD correlator library for GNSS software receivers”，GPSSolution 2004 vol.10)、或者通用图形处理器(参见Thomas Hobiger et al，“A GPUbased real-time GPS software receiver”，GPS Solution 2010，vol.14和Cailun Wuet al，“The Optimized Method and Algorithms in the PC-Based GNSS SoftwareReceiver Design”，ION-GNSS 2009)来进行跟踪过程中信号相关运行的加速。

根据文献Riedl，Bernhard，“Efficient GNSS Signal Acquisition with MassiveParallel Algorithms using OpenCL on GPUs”，ION-GNSS 2010，目前上述算法能够在目前高端的CPU处理器(Intel Core i7处理器)上实现20个通道的20MHz采样信号的跟踪运算，但仍然不足采用软件接收机的方法实现多卫星导航系统的多通道同时跟踪。

发明内容

本发明的技术解决问题：克服现有技术的不足，提出一种基于图形处理器的卫星导航系统比特压缩跟踪方法，该方法在图形处理器上使用比特压缩技术，能够将比特压缩算法的优点和图形处理器的并行计算性能进行结合，提高卫星导航系统跟踪算法的性能。

本发明的技术解决方案：一种基于图形处理器的卫星导航系统比特压缩跟踪方法，包括以下步骤：

(1)由接收机处理器负责生成或者读取载波相位查找表、扩频码查找表，并将传递到通用图形处理器中；

(2)接收机处理器读取需要进行跟踪的卫星导航系统中频信号，并将其传递到通用图形处理器的全局内存中；

(3)接收机处理器读取需要进行跟踪的各个通道的扩频码相位、载波相位、载波频率信息、卫星编号信息，并将其传递到通用图形处理器的全局内存中；

(4)根据跟踪任务划分区块，根据采样率和通用处理器处理字长划分线程，启动通用处理器处理程序，并将线程划分为数据交换线程(组)、相关计算线程和环路计算线程三类功能；

(5)各个通道的数据交换线程(组)至少负责获取本通道的扩频码相位、载波相位、载波频率信息、卫星编号，并存储在本通道的存储器中；

(6)各个通道的相关运算线程根据本通道的扩频码相位信息、载波相位信息和载波频率信息，负责进行跟踪相关运算；

(7)环路计算线程负责计算DLL、PLL/FLL、载噪比等信息，并将上述结果保存在本通道的存储器中；然后回到步骤(5)直到步骤(2)所读取的中频信号处理完毕；

(8)各个通道的数据交换线程(组)至少负责将本通道的扩频码相位、载波相位、载波频率信息、跟踪处理结果存储在图形处理器的全局存储器中；

(9)接收机处理器负责将各个通道的扩频码相位、载波相位、载波频率信息、跟踪处理结果从通用图形处理器存储器中取回后，继续完成卫星导航其他捕获、解码和定位等处理任务；

所述的步骤(4)中根据跟踪任务划分区块、根据系统采样率和通用处理器处理字长划分线程，并将线程划分为数据交换线程、相关计算线程和环路计算线程三类的步骤如下：

(a)划分区块数量为跟踪通道数量的整数倍；

(b)划分线程数量的基数为扩频码周期的采样点数和通用处理器字长除法取整(N/L)，其中N为每扩频码周期采样点个数，L为通用处理器字长。同时在码相位跳变点位置多划分一个线程、对于采样率不是通用处理器处理字长整数倍的情况，再多划分一个线程。当每扩频码周期采样点为N个，通用处理器处理字长为L个。当N为L的整数倍时，划分线程数量为N/L+1；当每扩频码周期采样点为N个，通用处理器处理字长为L个。当N不为L的整数倍时，划分线程数量为N/L+2；

(c)从划分的线程中选择一个或者多个线程作为数据交换线程(组)，该线程的功能是实现接收机处理器和通用图形处理器之间的数据交换功能，每个线程负责传送本通道的一个或者几个交换数据；

(d)所有线程均属于相关计算线程，这些线程都将参与跟踪过程中的相关运算；

(e)从划分的线程中选择一个线程作为环路计算线程，该线程的功能是负责导航接收机跟踪环路及其参数的计算，包括但不限于PLL跟踪环路、DLL跟踪环路、FLL跟踪环路和EKF跟踪环路、跟踪通道的信号载噪比等；

所述的步骤(6)中跟踪相关运算步骤包括：

(a)根据本通道的扩频码相位信息，确定本线程处理的采样点。对于跳变点之前(含跳变点)的线程t，处理的采样点范围为0～t*L，对于跳变点之后(含跳变点)的线程t’，处理的采样点范围为t*L～N；

(b)根据处理的采样点范围，获取本通道该采样点的对应中频信号数据，并进行数据位对齐；

(c)根据处理的采样点范围和本通道扩频码相位信息、载波相位信息和载波频率信息，获取载波数值、扩频码数值，并计算屏蔽掩码；

(d)据屏蔽掩码和上述的中频信号、扩频码数值、载波相位数值进行比特压缩的相关运算。每个线程会得中频信号、扩频码和载波I、Q支路相乘的结果，该结果可能包括但不具有限于：

(d.1)对于PSK调制的卫星导航系统，如果采用超前、即时和滞后码，会得到中频信号与超前码、即时码、滞后码和载波I、Q两个支路相乘的六个结果；

(d.2)对于PSK调制的卫星导航系统，如果采用即时和超前减滞后码，会得到中频信号与即时码、超前减滞后码和载波I、Q两个支路相乘的四个结果；

(d.3)对于BOC调制的卫星导航系统，如果采用超超前、超前、即时、滞后和超滞后码，会得到中频信号与超超前码、超前码、即时码、滞后码和超滞后码和载波I、Q两个支路相乘的十个结果；

(d.4)对于BOC调制的卫星导航系统，如果采用超超前减超滞后码、超前减滞后码和即时码，会得到中频信号与超超前码减超滞后码、超前减滞后码、即时码和载波I、Q两个支路相乘的六个结果；

(e)对于两部分计算线程分别进行矢量缩减求和运算，分别得到两部分的I、Q支路累加结果。

本发明与现有技术相比的优点在于：

(1)与传统的基于通用处理器的卫星导航系统跟踪算法相比，本发明采用了比特压缩算法，因此在GPU上启动的线程数量只是传统算法的1/32或者更少。更少的处理线程会加快图形处理器任务调度，并且利用比特压缩算法可以实现多个采样点相关运算的并行性，从而加快跟踪运算的速度。

(2)与传统的基于中央处理器的卫星导航系统跟踪算法相比，在通用图形处理器上能够实现对信号处理的并行性。由于通用图形处理器具有多个处理核心，支持MIMD指令系统，因此可以同时完成对多个比特压缩后的信号进行并行处理，从而加快跟踪运算的速度。

(3)本发明采用独立的线程分别处理跳变点数据和非整数位置数据，不会引起图形处理器内存访问不对齐，内存读取延迟小；基于此方法处理跳变点位置的中频信号处理和非整数位置的中频信号数据基本不会增加算法处理的复杂度。

(4)本发明采用处理掩码使跳变点和非整数位置数据的处理和其他数据处理方法基本一致。因此在图形处理器的指令执行系统中，对这两个特殊的采样点的处理指令和其他采样点的处理指令可以保持一致，从而不会影响图形处理器流水线的执行，保证算法执行的效率。

(5)将传统的相关累加运算分解为两组矢量缩减求和运算，使其能够并行地计算相关累加结果。该算法既保证了相关累加过程不会跨越跳变点，同时算法复杂度与一组矢量缩减求和运算基本相同，因此能够提高跟踪算法的执行效率。

(6)将跟踪环路放在GPU的环路计算线程上进行，可以使跟踪过程中产生的鉴别器数值、IQ支路相关结果以及跟踪过程中的码相位、载波相位和载波频率的变化信息都不需要传送回中央处理器进行处理，从而减少中央处理器和图形处理器的数据交换。由于全局内存拷贝在图形处理器中需要使用较长时间，因此本发明的方法能够极大地缩短跟踪算法执行的时间。

附图说明

图1为传统的卫星导航系统中频信号跟踪方法处理流程图；

图2为本发明的跟踪方法处理流程图；

图3为本发明的基于比特压缩线程任务分配示意图；

图4为本发明的线程功能划分示意图；

图5为本发明的相关运算线程流程图；

图6为本发明的相关运算过程中屏蔽掩码计算示意图；

图7为本发明的两组矢量缩减求和运算示意图。

具体实施方式

本发明中处理的数据是卫星导航系统的中频数字信号，该信号是由卫星导航系统射频信号通过硬件进行下变频后采样得到的。数字中频信号的采样频率必须符合过采样或带通采样定律，通常情况大于5MHz。对数字中信信号的处理是基于数据块进行的，即每次从采集系统或者文件系统中获取一定数量的采样数据(譬如40毫秒的采样数据)。假设扩频码周期为1毫秒，则对于40毫秒的采样数据，跟踪处理算法需要完成40次相关运算处理。接收机会对于每一组采样数据的所有通道进行并行处理，如图2所示。本发明的具体步骤如下：

(1)由接收机处理器负责生成或者读取载波相位查找表、扩频码查找表，并将传递到通用图形处理器中；

扩频码查找表的生成算法与使用的扩频码有关：GPS L1频段的扩频码生成方法参见ICD200c和IS800、GPS L2频段的扩频码生成方法参见ICD200d、GPS L5频段的扩频码生成方法参见IS705、GLONASS系统的扩频码生成方法参见GLONASS ICD v5.1、伽利略系统的扩频码生成方法参加Galileo OS SIS ICD。生成的扩频码查找表以卫星编号和码相位偏移作为索引，通过这两个索引可以查找到指定编号和码相位偏移的一个周期扩频码。以GPS L1频段的C/A码为例，当处理字长为32比特时，其索引为32个导航卫星和32个码相位偏移，其查找表内容为首尾相接的C/A扩频码。载波相位查找表的生成算法与卫星导航系统中频信号的频率有关。生成的载波相位查找表以载波频率和载波相位作为索引，通过这两个索引可以查找指定频率和相位的一个周期的本地复现载波信号。载波相位查找表的载波频率间隔可以设置为50Hz，载波相位间隔可以设置为π/32。在卫星导航接收机运行后，这两个查找表可以保存在通用图形处理器的全局内存中。

(2)接收机处理器读取需要进行跟踪的卫星导航系统中频信号，并将其传递到通用图形处理器的全局内存中。

接收机处理器从采样缓冲器或者文件中读取需要进行跟踪的卫星导航系统中频信号，其读取长度为一次需要处理的数据块长度(譬如40毫秒的中频数据块)，并将其传递到通用图形处理器的全局内存中。在本数据块处理完成后，接收机处理器需要传递下一个需要处理的数据块。当采样频率为12MHz，量化位数为2比特，长度为40毫秒的中频数据需要占用120000字节。

(3)接收机处理器读取需要进行跟踪的各个通道的扩频码相位、载波相位、载波频率信息、卫星编号信息，并将其传递到通用图形处理器的全局内存中。

接收机处理器将各个跟踪通道的扩频码相位τ、载波相位φ、载波频率信息f、卫星编号信息n保存到一个数组中，并将其传递到通用图形处理器的全局内存中。跟踪算法利用步骤(2)传递的中频信号数据块和步骤(3)传递的跟踪初始信息对中频信号进行处理。待处理完成后，再使用步骤(2)和步骤(3)传递中频信号数据和下一组跟踪通道信息。

(4)根据跟踪任务划分区块，根据采样率和通用处理器处理字长划分线程，启动通用处理器处理程序，并将线程划分为数据交换线程(组)、相关计算线程和环路计算线程三类功能。

图形处理器计算任务的区块数量等于所需要处理的通道数，而划分线程数量的基数为扩频码周期的采样点数和通用处理器字长除法取整(N/L)，其中N为每扩频码周期采样点个数，L为通用处理器字长。同时在码相位跳变点位置多划分一个线程、对于采样率不是通用处理器处理字长整数倍的情况，再多划分一个线程。即当N为L的整数倍时，划分线程数量为N/L+1；当N不为L的整数倍时，划分线程数量为N/L+2。譬如，当采样率为12MHz、图形处理器处理字长为32比特、扩频码周期为1毫秒时，划分线程数量为376个(12000/32+1)；当采样率为5.172MHz、图形处理器处理字长为32比特、扩频码周期为1毫秒时，划分线程数量为163个(5172/32+2)。其划分方法如图3所示。

线程按照功能被分为数据交换线程(组)、相关计算线程和环路计算线程三类，其分类方法如图4所示。当一个区块存在376个线程的时候，可以将0～3号线程划分为数据交换线程组，每个线程负责通道内一个跟踪参数的交换，也可以只将0号线程划分为数据交换线程，负责通道内所有跟踪参数的交换；将0号线程划分为环路计算线程，负责计算跟踪环路参数；所有线程(0～375号线程)均属于相关计算线程，均参与信号相关积分的计算过程。

(5)各个通道的数据交换线程(组)至少负责获取本通道的扩频码相位、载波相位、载波频率信息、卫星编号，并存储在本通道的存储器中；

数据交换线程负责读取本通道的扩频码相位τ、载波相位φ、载波频率信息f、卫星编号n，并存储在本通道的存储器中。同时该线程还负责将所获取到的参数值换算为索引编号后，保存在寄存器中，以便于相关运算线程进行查询。譬如，如果处理器字长为32位、载波相位被划分为64个单元、载波频率的间隔为50Hz，此时扩频码相位索引为τ/32、载波相位索引为

载波频率信息索引为f/50、卫星编号索引为n。

(6)各个通道的相关运算线程根据本通道的扩频码相位τ、载波相位φ、载波频率信息f、卫星编号n，进行跟踪相关运算。具体跟踪相关运算的过程分为划分线程、数据位对齐、掩码计算、相乘运算、矢量缩减求和五个步骤，其流程如图5所示：

(a)每个通道的相关运算线程需要处理L个采样点数据(L为通用图形处理器字长)，其采样点范围为t_a～t_b。对于编号为t的线程而言，采样点范围的两侧位置分别用a和b表示，其取值分别为：

$\left\{\begin{array}{cc}a=t\×L,b=\min ((t+1)\×L,N)& t\≤\frac{\mathrm{\τ}}{L}+1\\ a=(t-1)\×L,b=\min (t\×L,N)& t>\frac{\mathrm{\τ}}{L}+1\end{array}\right.$

采样点范围为t_a～t_b的信号在中频信号和查找表中的位置为p，p＝a/32。对于编号为t的线程而言，p的取值为：

$\left\{\begin{array}{cc}p=t& t\≤\frac{\mathrm{\τ}}{L}+1\\ p=t-1& t>\frac{\mathrm{\τ}}{L}+1\end{array}\right.;$

(b)根据处理的采样点范围在查找表中的位置p，获取本通道该采样点的对应中频信号数据的第p个数值，并进行数据位对齐。数据位对齐主要针对非1比特采样的中频信号进行，通过数据位对齐将多比特采样的数据分成几个不同存储数组；

(c)根据处理的采样点范围和本通道扩频码相位信息、载波相位信息和载波频率信息，获取载波数值、扩频码数值，并计算屏蔽掩码；

使用步骤(5)计算得出的索引值获取扩频码查找表、载波相位查找表对应的信号数组，并取出其中第p个信号作为扩频码数值和载波数值，其中p值在步骤(a)中计算得到。计算屏蔽掩码m如图6所示。计算屏蔽掩码仅在跳变点位置的两个线程和非整数部分的一个线程中存在，其计算过程需要先计算非对齐数量r。当图形处理器字长为32比特时候，r的计算方法如下：

$\left\{\begin{array}{cc}r=\mathrm{\τ}\%32& t=\frac{\mathrm{\τ}}{L}+1,t=\frac{N}{L}+1\\ r=32-(\mathrm{\τ}\%32)& t=\frac{\mathrm{\τ}}{L}+2\end{array}\right.$

计算屏蔽掩码m的计算方法如下，

$\left\{\begin{array}{cc}m=0\mathrm{xFFFFFFFF}<<r& t=\frac{\mathrm{\&tau;}}{L}+1,t=\frac{N}{L}+1\\ m=0\mathrm{xFFFFFFFF}>>r& t=\frac{\mathrm{\&tau;}}{L}+2\end{array}\right.;$

(d)据屏蔽掩码和上述的中频信号、扩频码数值、载波相位数值进行比特压缩的相乘运算。每个线程会得中频信号、扩频码和载波I、Q支路相乘的结果。由于采用比特压缩的相关运算，因此相关运算中的乘法运算可以使用位异或等位运算代替，一个处理字长内的加法运算可以使用位计数运算代替。

以GPS系统PSK调制的C/A码为例，传统的32个点的超前码、及时码和滞后码相关积分运算如下所示，

$I_E=\underset{i=1}{\overset{32}{\mathrm{\&Sigma;}}}s\left(p\right)C(p-\mathrm{\&tau;}+{\mathrm{\&Delta;}}_1)\cos (2\mathrm{\&pi;fp}+\mathrm{\&phi;})$

$I_P=\underset{i=1}{\overset{32}{\mathrm{\&Sigma;}}}s\left(p\right)C(p-\mathrm{\&tau;})\cos (2\mathrm{\&pi;fp}+\mathrm{\&phi;})$

$I_L=\underset{i=1}{\overset{32}{\mathrm{\&Sigma;}}}s\left(p\right)C(p-\mathrm{\&tau;}-{\mathrm{\&Delta;}}_1)\cos (2\mathrm{\&pi;fp}+\mathrm{\&phi;})$

$Q_E=\underset{i=1}{\overset{32}{\mathrm{\&Sigma;}}}s\left(p\right)C(p-\mathrm{\&tau;}+{\mathrm{\&Delta;}}_1)\sin (2\mathrm{\&pi;fp}+\mathrm{\&phi;})$

$Q_P=\underset{i=1}{\overset{32}{\mathrm{\&Sigma;}}}s\left(p\right)C(p-\mathrm{\&tau;})\cos (2\mathrm{\&pi;fp}+\mathrm{\&phi;})$

$Q_L=\underset{i=1}{\overset{32}{\mathrm{\&Sigma;}}}s\left(p\right)C(p-\mathrm{\&tau;}-{\mathrm{\&Delta;}}_1)\cos (2\mathrm{\&pi;fp}+\mathrm{\&phi;})$

其中，I_E、I_P、I_L、Q_E、Q_P、Q_L分别表示超前码、及时码和滞后码的I、Q支路相关结果，s(p)表示p时刻的信号采样，C表示卫星导航扩频码，τ表示卫星到接收载体的传输延迟，Δ₁表示PSK扩频方法的相关器间隔，f表示载波频率，φ表示载波相位。而在本发明的运算过程可以被替换为：

I_E＝popcnt(m×s(p)×C(p-τ+Δ₁)×cos(2πfp+φ))

I_P＝popcnt(m×s(p)×C(p-τ)×cos(2πfp+φ))

I_L＝popcnt(m×s(p)×C(p-τ-Δ₁)×cos(2πfp+φ))

Q_E＝popcnt(m×s(p)×C(p-τ+Δ₁)×sin(2πfp+φ))

Q_P＝popcnt(m×s(p)×C(p-τ)×sin(2πfp+φ))

Q_L＝popcnt(m×s(p)×C(p-τ-Δ₁)×sin(2πfp+φ))

其中，popcnt为比特求和运算。

以GPS系统PSK调制的C/A码为例，传统的32个点的超前码减滞后码和及时码相关积分运算如下所示，

$I_{E-L}=\underset{i=1}{\overset{32}{\mathrm{\&Sigma;}}}s\left(p\right)C_{E-L}(p-\mathrm{\&tau;}+{\mathrm{\&Delta;}}_1,p-\mathrm{\&tau;}-{\mathrm{\&Delta;}}_1)\cos (2\mathrm{\&pi;fp}+\mathrm{\&phi;})$

$I_P=\underset{i=1}{\overset{32}{\mathrm{\&Sigma;}}}s\left(p\right)C(p-\mathrm{\&tau;})\cos (2\mathrm{\&pi;fp}+\mathrm{\&phi;})$

$Q_{E-L}=\underset{i=1}{\overset{32}{\mathrm{\&Sigma;}}}s\left(p\right)C_{E-L}(p-\mathrm{\&tau;}+{\mathrm{\&Delta;}}_1,p-\mathrm{\&tau;}-{\mathrm{\&Delta;}}_1)\sin (2\mathrm{\&pi;fp}+\mathrm{\&phi;})$

$Q_P=\underset{i=1}{\overset{32}{\mathrm{\&Sigma;}}}s\left(p\right)C(p-\mathrm{\&tau;})\cos (2\mathrm{\&pi;fp}+\mathrm{\&phi;})$

其中，C_E-L表示卫星导航扩频码超前码减滞后码的结果。而在本发明的运算过程可以被替换为：

I_E-L＝popcnt(m×s(p)×C_E-L(p-τ+Δ₁，p-τ-Δ₁)×cos(2πfp+φ))

I_P＝popcnt(m×s(p)×C(p-τ)×cos(2πfp+φ))

Q_E-L＝popcnt(m×s(p)×C_E-L(p-τ+Δ₁，p-τ-Δ₁)×sin(2πfp+φ))

Q_P＝popcnt(m×s(p)×C(p-τ)×sin(2πfp+φ))

以伽利略系统BOC调制的扩频码为例，传统的32个点的超超前码、超前码、及时码、滞后码和超滞后码相关积分运算如下所示：

$I_{\mathrm{VE}}=\underset{i=1}{\overset{32}{\mathrm{\&Sigma;}}}s\left(p\right)C(p-\mathrm{\&tau;}+{2\mathrm{\&Delta;}}_2)\cos (2\mathrm{\&pi;fp}+\mathrm{\&phi;})$

$I_E=\underset{i=1}{\overset{32}{\mathrm{\&Sigma;}}}s\left(p\right)C(p-\mathrm{\&tau;}+{\mathrm{\&Delta;}}_2)\cos (2\mathrm{\&pi;fp}+\mathrm{\&phi;})$

$I_P=\underset{i=1}{\overset{32}{\mathrm{\&Sigma;}}}s\left(p\right)C(p-\mathrm{\&tau;})\cos (2\mathrm{\&pi;fp}+\mathrm{\&phi;})$

$I_L=\underset{i=1}{\overset{32}{\mathrm{\&Sigma;}}}s\left(p\right)C(p-\mathrm{\&tau;}-{\mathrm{\&Delta;}}_2)\cos (2\mathrm{\&pi;fp}+\mathrm{\&phi;})$

$I_{\mathrm{VL}}=\underset{i=1}{\overset{32}{\mathrm{\&Sigma;}}}s\left(p\right)C(p-\mathrm{\&tau;}-{2\mathrm{\&Delta;}}_2)\cos (2\mathrm{\&pi;fp}+\mathrm{\&phi;})$

$Q_{\mathrm{VE}}=\underset{i=1}{\overset{32}{\mathrm{\&Sigma;}}}s\left(p\right)C(p-\mathrm{\&tau;}+{2\mathrm{\&Delta;}}_2)\sin (2\mathrm{\&pi;fp}+\mathrm{\&phi;})$

$Q_E=\underset{i=1}{\overset{32}{\mathrm{\&Sigma;}}}s\left(p\right)C(p-\mathrm{\&tau;}+{\mathrm{\&Delta;}}_2)\sin (2\mathrm{\&pi;fp}+\mathrm{\&phi;})$

$Q_P=\underset{i=1}{\overset{32}{\mathrm{\&Sigma;}}}s\left(p\right)C(p-\mathrm{\&tau;})\cos (2\mathrm{\&pi;fp}+\mathrm{\&phi;})$

$Q_L=\underset{i=1}{\overset{32}{\mathrm{\&Sigma;}}}s\left(p\right)C(p-\mathrm{\&tau;}-{\mathrm{\&Delta;}}_2)\cos (2\mathrm{\&pi;fp}+\mathrm{\&phi;})$

$Q_{\mathrm{VL}}=\underset{i=1}{\overset{32}{\mathrm{\&Sigma;}}}s\left(p\right)C(p-\mathrm{\&tau;}-{2\mathrm{\&Delta;}}_2)\cos (2\mathrm{\&pi;fp}+\mathrm{\&phi;})$

其中，Δ₂表示BOC扩频方法的相关器间隔。而在本发明的运算过程可以被替换为：

I_VE＝popcnt(m×s(p)×C(p-τ+2Δ₂)×cos(2πfp+φ))

I_E＝popcnt(m×s(p)×C(p-τ+Δ₂)×cos(2πfp+φ))

I_P＝popcnt(m×s(p)×C(p-τ)×cos(2πfp+φ))

I_L＝popcnt(m×s(p)×C(p-τ-Δ₂)×cos(2πfp+φ))

I_VL＝popcnt(m×s(p)×C(p-τ-2Δ₂)×cos(2πfp+φ))

Q_VE＝popcnt(m×s(p)×C(p-τ+2Δ₂)×sin(2πfp+φ))

Q_E＝popcnt(m×s(p)×C(p-τ+Δ₂)×sin(2πfp+φ))

Q_P＝popcnt(m×s(p)×C(p-τ)×sin(2πfp+φ))

Q_L＝popcnt(m×s(p)×C(p-τ-Δ₂)×sin(2πfp+φ))

Q_VL＝popcnt(m×s(p)×C(p-τ-2Δ₂)×sin(2πfp+φ))

以伽利略系统BOC调制的扩频码为例，传统的32个点的超超前减超滞后码、超前减滞后码、及时码相关积分运算如下所示：

$I_{\mathrm{VE}-\mathrm{VL}}=\underset{i=1}{\overset{32}{\mathrm{\&Sigma;}}}s\left(p\right)C_{E-L}(p-\mathrm{\&tau;}-2{\mathrm{\&Delta;}}_2,p-\mathrm{\&tau;}+{2\mathrm{\&Delta;}}_2)\cos (2\mathrm{\&pi;fp}+\mathrm{\&phi;})$

$I_{E-L}=\underset{i=1}{\overset{32}{\mathrm{\&Sigma;}}}s\left(p\right)C_{E-L}(p-\mathrm{\&tau;}-{\mathrm{\&Delta;}}_2,p-\mathrm{\&tau;}+{\mathrm{\&Delta;}}_2)\cos (2\mathrm{\&pi;fp}+\mathrm{\&phi;})$

$I_P=\underset{i=1}{\overset{32}{\mathrm{\&Sigma;}}}s\left(p\right)C(p-\mathrm{\&tau;})\cos (2\mathrm{\&pi;fp}+\mathrm{\&phi;})$

$Q_{\mathrm{VE}-\mathrm{VL}}=\underset{i=1}{\overset{32}{\mathrm{\&Sigma;}}}s\left(p\right)C_{E-L}(p-\mathrm{\&tau;}-2{\mathrm{\&Delta;}}_2,p-\mathrm{\&tau;}+{2\mathrm{\&Delta;}}_2)\sin (2\mathrm{\&pi;fp}+\mathrm{\&phi;})$

$Q_{E-L}=\underset{i=1}{\overset{32}{\mathrm{\&Sigma;}}}s\left(p\right)C_{E-L}(p-\mathrm{\&tau;}-{\mathrm{\&Delta;}}_2,p-\mathrm{\&tau;}+{\mathrm{\&Delta;}}_2)\sin (2\mathrm{\&pi;fp}+\mathrm{\&phi;})$

$Q_P=\underset{i=1}{\overset{32}{\mathrm{\&Sigma;}}}s\left(p\right)C(p-\mathrm{\&tau;})\cos (2\mathrm{\&pi;fp}+\mathrm{\&phi;})$

而在本发明的运算过程可以被替换为：

I_VE-VL＝popcnt(m×s(p)×C_E-L(p-τ+2Δ₂，p-τ-2Δ₂)×cos(2πfp+φ))

I_E-L＝popcnt(m×s(p)×C_E-L(p-τ+Δ₂，p-τ-Δ₂)×cos(2πfp+φ))

I_P＝popcnt(m×s(p)×C(p-τ)×cos(2πfp+φ))

Q_VE-VL＝popcnt(m×s(p)×C_E-L(p-τ+2Δ₁，p-τ-2Δ₁)×sin(2πfp+φ))

Q_E-L＝popcnt(m×s(p)×C_E-L(p-τ+Δ₁，p-τ-Δ₁)×sin(2πfp+φ))

Q_P＝popcnt(m×s(p)×C(p-τ)×sin(2πfp+φ))；

(e)对于两部分计算线程分别进行矢量缩减求和运算，分别得到两部分的I、Q支路累加结果；

为了避免积分跨越跳变点，对步骤(a)所示的两组线程的相乘结果分别进行矢量缩减求和运算。其相乘结果可以从步骤(d)中获得，而矢量缩减求和运算后即可得到两部分的I、Q支路累加结果。两组矢量缩减求和运算的方法如图7所示。

(7)环路计算线程负责计算DLL、PLL/FLL、载噪比等信息，并将上述结果保存在本通道的存储器中；然后回到步骤(5)直到步骤(2)所读取的中频信号处理完毕；

环路计算线程利用步骤(6)中计算得到的I、Q支路累加结果，通过DLL、PLL/FLL跟踪环路的计算得到扩频码相位修正量Δτ、载波频率修正量Δf、载波相位修正量Δφ。有关不同卫星导航系统的跟踪算法参见Kai Borre et al，“A Software-Defined GPS andGalileo Receiver：A Single-Frequency Approach”。将该修正量与输入的参数相加后，得到下一个扩频码周期数据跟踪过程所需要使用的跟踪参数。然后返回步骤(5)重复进行跟踪运算，直到步骤(2)所读取的中频信号处理完成。

(8)各个通道的数据交换线程(组)至少负责将本通道的扩频码相位、载波相位、载波频率信息、跟踪处理结果存储在图形处理器的全局存储器中；

(9)接收机处理器负责将各个通道的扩频码相位、载波相位、载波频率信息、跟踪处理结果从通用图形处理器存储器中取回后，继续完成卫星导航其他捕获、解码和定位等处理任务。

本发明未详细阐述部分属于本领域公知技术。

Claims (1)

1.一种基于通用图形处理器的卫星导航系统比特压缩跟踪方法，其特征在于步骤如下：

(1)由接收机处理器负责生成或者读取载波相位查找表、扩频码查找表，并将载波相位查找表、扩频码查找表传递到通用图形处理器中；

(2)接收机处理器读取需要进行跟踪的卫星导航系统中频信号，并将中频信号传递到通用图形处理器的全局内存中；

(3)接收机处理器读取需要进行跟踪的各个通道的扩频码相位、载波相位、载波频率信息和卫星编号信息，并将扩频码相位、载波相位、载波频率信息和卫星编号信息传递到通用图形处理器的全局内存中；

(4)根据跟踪任务划分区块、采样率和通用处理器处理字长划分线程，并将线程划分为数据交换线程、相关计算线程和环路计算线程三类；

(5)各个通道的数据交换线程至少负责获取本通道的扩频码相位、载波相位、载波频率信息和卫星编号，并存储在本通道的存储器中；

(6)各个通道的相关计算线程根据本通道的扩频码相位信息、载波相位信息和载波频率信息，负责进行跟踪相关运算；

(7)环路计算线程负责计算DLL、PLL/FLL和载噪比信息，并将所述DLL、PLL/FLL和载噪比信息保存在本通道的存储器中；然后回到步骤(5)直到步骤(2)所读取的中频信号处理完毕；

(8)各个通道的数据交换线程至少负责将本通道的扩频码相位、载波相位、载波频率信息、跟踪处理结果存储在通用图形处理器的全局存储器中；

(9)接收机处理器负责将各个通道的扩频码相位、载波相位、载波频率信息和跟踪处理结果从通用图形处理器的全局存储器中取回后，继续完成卫星导航其他捕获、解码和定位处理任务；

所述的步骤(4)中根据跟踪任务划分区块、根据系统采样率和通用处理器处理字长划分线程，并将线程划分为数据交换线程、相关计算线程和环路计算线程三类的步骤如下：

(4a)划分区块数量为跟踪通道数量的整数倍；

(4b)划分线程数量的基数为扩频码周期的采样点数和通用处理器字长按照N/L除法取整，其中N为每扩频码周期采样点个数，L为通用处理器字长，同时在码相位跳变点位置多划分一个线程、对于采样率不是通用处理器处理字长整数倍的情况，再多划分一个线程；当每扩频码周期采样点为N个，通用处理器处理字长为L个；当N为L的整数倍时，划分线程数量为N/L+1；当每扩频码周期采样点为N个，通用处理器处理字长为L个；当N不为L的整数倍时，划分线程数量为N/L+2；

(4c)从划分的线程中选择一个或者多个线程作为数据交换线程，该线程的功能是实现接收机处理器和通用图形处理器之间的数据交换功能，每个线程负责传送本通道的一个或者几个交换数据；

(4d)所有线程均属于相关计算线程，这些线程都将参与跟踪过程中的相关运算；

(4e)从划分的线程中选择一个线程作为环路计算线程，该线程的功能是负责导航接收机跟踪环路及其参数的计算，包括PLL跟踪环路、DLL跟踪环路、FLL跟踪环路和EKF跟踪环路、跟踪通道的信号载噪比；

所述的步骤(6)中跟踪相关运算步骤包括：

(6a)根据本通道的扩频码相位信息，确定本线程处理的采样点，对于跳变点之前的线程t，处理的采样点范围为0-t*L，对于跳变点之后的线程t’，处理的采样点范围为t*L-N；

(6b)根据处理的采样点范围，获取本通道该采样点的对应中频信号数据，并进行数据位对齐；

(6c)根据处理的采样点范围和本通道扩频码相位信息、载波相位信息和载波频率信息，获取载波数值、扩频码数值，并计算屏蔽掩码；

(6d)根据屏蔽掩码和上述的中频信号数据、扩频码数值、载波数值进行比特压缩的相关运算；

(6e)对于两部分计算线程分别进行矢量缩减求和运算，分别得到两部分的I、Q支路累加结果。

Images