采用FPGA技术的北斗B3I同步算法

沈周锋

(漳州职业技术学院电子工程学院, 福建 漳州 363000)

0 引言

当前，联合国确定的全球四大导航系统包括我国的北斗卫星导航系统(Beidou navigation satellite system，BDS)、 美国的全球定位系统(global positioning system，GPS)、 俄罗斯的格洛纳斯系统(Global navigation satellite system, GLONASS)和欧盟伽利略(Galileo). 从一代、 二代至三代的升级过程中，BDS系统已经具备全球服务的能力. 2018年BDS向民用设备开放了B3I频段，利用码率更高的测距码将导航电文频谱扩展到更宽的频段，使其具有更强的抗干扰、 更精准的导航定位性能[1]. 信号解调单元作为终端设备核心之一，决定了整个系统的抗干扰能力，大大影响导航电文分析单元能否正常工作. 文献[2]提出了利用串行相干解调算法检测测距码相位，该算法需要加长累加时长才能充分发挥B3I的强抗干扰性能，计算量巨大. 文献[3-6]提出了利用FFT求取互相关序列的方法，其实质是以提高电路复杂度为代价换取计算量的减少，直接应用于B3I频段时，测距码码长较大将导致电路开销急剧上升，难以实现. 文献[7-8]提出非相干累加和差分累加的方法，虽然减少了计算量，但是抗噪声性能仍然无法和同等累加时长的相干累加相比拟. 文献[9]提出了GPS半比特交替的弱信号捕获算法，由于北斗系统导航电文利用NH(neumann-hoffman)码进行二次编码，二次编码跳变的问题使该算法无法正常工作.

累加时长20 ms的相干累加算法具有较高的抗噪声性能，在弱信号环境下能够胜任信号捕获. 该算法在捕获、 跟踪单元被多次使用. 由于计算量巨大，直接实现将导致导航终端信号捕获跟踪耗时很长，且占用大量的芯片资源. 本研究对算法进行改进，使其在保留原始性能的同时减少计算量，且满足FPGA并行化和流水线化结构的要求，具有较大的实用价值.

1 信号模型和算法分析

1.1 北斗B3I信号同步原理

图1 编码示意图Fig.1 Coding schematic diagram

北斗B3I频段D1导航电文码率为50 bit·s-1，调制到NH码上，形成1 kbit·s-1的码流. 然后利用码长10 230 bit的测距码对其进行直接序列扩频，生成10.23 Mbit·s-1的序列. 该序列与1 268.520 MHz载波进行BPSK调制，形成带宽为20.46 MHz的射频信号. NH码码长为20 bit，由0和1两种码组成，码片宽度为1 ms.C(n)、D(n)和NH(n)采用模二加的方式进行调制, 调制方式如图1所示. 在算法设计中，一般采用数字-1表示码片0，采用数字1表示码片1. 不难得出，两次模二加结果与两次乘法结果相同，基带数字信号可以表示为C(n)D(n)NH(n).

射频信号经过射频前端电路下变频生成中频信号. 采用相位差为90°的两路本地载波分别与中频信号混频，得到同相、 正交两路基带信号[10]. 频偏捕获单元合理控制本地载波频率，将频偏补偿掉, 基带信号经过采样并合成复数形式. 频偏完全被补偿的理想情况下，接收序列R(n)可表示为：

R(n)=AC(n)D(n)NH(n)

(1)

其中：A表示接收信号幅度. 令PN(n)=C(n)NH(n)，则PN(n)相当于码率10.23 Mbit·s-1，码长(10 230 ×20) bit，周期20 ms的伪随机序列，其自相关曲线和互相关曲线如图2所示. 横轴表示相关运算滑动位置，纵轴表示相关运算值. PN(n)继承了NH(n)和C(n)良好的自相关特性，同时继承了C(n)的伪正交特性. 图2(a)为自相关曲线，具有尖锐的相关峰. 以1号卫星和2号卫星PN序列为例，互相关曲线如图2(b)所示，曲线上无明显相关峰.

图2 PN(n)自相关和互相关特性曲线Fig.2 PN(n) autocorrelation and cross correlation characteristic curve

接收信号可以看成PN(n)对导航电文D(n)的直接序列扩频，R(n)重写为：

R(n)=AD(n)PN(n)

(2)

本地伪码发生器产生伪随机序列，表示为PNL(n). 与R(n)进行互相关，得到以下序列:

(3)

式(3)中，当PN(n)和PNL(n)对齐时，Z(n)出现相关峰, 峰值位置即为NH码周期的起始位置.

在实际系统中，下采样得到采样率为10.23 MHz的码流，从任意时刻开始截取409 199个采样点与PNL(n)做滑动互相关运算. 最多经过204 600次运算，必能检测到相关峰. 直接计算Z(n)计算量相当大，且难以将运算电路进行流水线化处理，需要对算法进行改进.

1.2 分部相关算法原理

图3 互相关运算示意图Fig.3 Cross-correlation operation schematic diagram

式(3)互相关运算示意图如图3所示. 本地PN序列由20个周期的测距码连接而成，每个测距码周期由一个格子表示，填充阴影的格子表示该测距码周期对应的NH码片等于-1，其他格子表示对应的NH码片等于1. 本地PN码与接收序列R(n)进行相乘并累加的过程，可以看成将R(n)分成20段，分别与20个周期的测距码相乘并累加，再将20个累加结果加权求和. PN码与接收序列R(n)滑动互相关过程中，存在很多相同的运算，浪费系统资源. 因此提出分部相关算法，适当记录每个段的运算结果，达到减少重复运算的目的.

Z(n)的值可由累加时长1 ms的互相关结果进行加减运算求得.Z(n)可整理为：

(4)

Z(n)由20个累加项与NH码相乘并求和得到. 每个累加项对应一个累加时长1 ms的互相关运算，表示为：

(5)

代入式(4)，则：

(6)

只需计算出Y(n)序列，抽取特定位置的数据与NH码进行乘法和加法运算即可得到Z(n). 计算Z(n)过程中，Y(n)序列的每个数值被多次调用. 计算过程的中间结果被充分利用，大大减小了计算量. 分部相关算法有如下3个具体步骤：

步骤一按照式(5)计算累加时长1 ms的互相关序列Y(n).n的取值从0至(10 230×39-1).

步骤二按照式(6)抽取特定位置的Y值与NH码相乘并累加，计算Z(n).n的取值从0至(10 230×20-1).

步骤三求Z(n)虚部和实部的平方和，与特定阈值对比, 超过阈值即认定为相关峰，该点位置即为NH码周期起始位置，或称为导航电文比特起始位置.

加长累加时长至20 ms，使Z(n)序列信噪比相比接收序列提升(10 230×20)倍，同时相关峰峰值对频偏更加敏感[11]. 频偏存在时，相关峰峰值表示为：

(7)

2 分部相关算法的FPGA实现

图4 相干累加电路框图Fig.4 Block diagram of coherent accumulation circuit

接收序列与本地伪码相乘并累加，求取序列互相关值. 在FPGA中，乘法器电路复杂度远大于加法器，且片上集成的复数乘法单元数量有限. 本地伪码的值为-1和+1，调用乘法器单元将会极大地浪费芯片资源，也没必要用电路复杂度更高的FFT IP核. 只需在累加时根据本地伪码的符号，切换加减运算即可，如图4所示. 当本地伪码为正时，模拟选通开关向上接通，复数序列进入加法器，与累加结果相加后存入累加结果寄存器. 反之，当本地伪码为负时，选通开关向下接通，复数序列进入减法器，与累加结果相减后存入累加结果寄存器. 电路在时钟信号驱动下，执行相干累加运算. 达到既定累加次数后，输出使能引脚置一，累加结果从数据端口输出.

NH码相位检测实现框图如图5所示. 射频前端将接收信号下变频并采样为数字基带信号后，送入本检测电路. 复振荡信号发生器产生特定频率的信号与接收序列共轭相乘，进一步补偿剩余频偏，然后存入信号缓存. 10 230次相干累加器从信号缓存中读取所需接收序列和本地测距码，生成Y(n)信号，存入Y(n)序列缓存中. 20次相干累加器读取Y(n)的值，与NH码进行相干累加，生成Z(n). 相关峰检测单元求取Z(n)的实部和虚部的平方和，与阈值对比，检测当前值是否为相关峰. 最终根据相关峰位置换算出NH码周期起始位置.

图5 NH码相位检测实现框图Fig.5 Realization block diagram of NH code phase detection

两种相干累加运算由不同运算单元独立完成. 合理控制Y(n)序列生成顺序，两种运算可同时执行，实现流水线化，加快NH码相位检测速度[12]. 由式(6)可知，每个Z(n)是由20个间隔10 230个点的Y(n)序列与本地NH码相干累加求得. 例如：Y(10 230×0)，Y(10 230×1)，…，Y(10 230×38)等39个点用于计算Z(10 230×0)，Z(10 230×1)，…，Z(10 230×19)等20个点. 39个Y(n)值存入序列缓存后，立即开始下一组Y(n)序列的生成，同时启动20次相干累加器计算Z(n)序列，实现两种相干运算流水线化结构. 由于Z(n)序列生成速度大于Y(n)序列，Y(n)序列缓存被新数据覆盖前Z(n)序列已经计算完成. 计算完Z(n)后，Y(n)序列缓存数据可丢弃，用于存储下一轮的Y(n)序列，大大节省内部存储空间.

10 230次相干累加相比20次相干累加，计算量较大，逐个计算Y(n)数据将大大降低运算速度. 须采用并行化结构解决此“瓶颈”，加快运算速度. 如图5所示，实例化39个10 230次相干累加器，同时计算39个间隔为10 230个采样点的Y(n)的值，存入Y(n)序列缓存中. 39个电路同步运行，同一时刻所需的接收序列点间隔10 230个采样点，接收序列缓存须保证间隔10 230个采样点的接收序列数据并行输出. 将间隔10 230个采样点的接收序列值存入信号缓存的同一个字的不同位. 当接收序列缓存读取到地址信号，则输出对应字节，39个数据包含在同一字节的不同位中，并行输出. 同一时刻所需的本地测距码相同，只需要一个本地测距码发生器即可满足39个相干累加器的运算需求.

具体的计算流程图如图6所示，最下方的箭头为时间轴，各轮运算沿着时间轴从左到右顺序执行. 每一轮运算得到20个Z(n)数据，最多经过10 230轮，可完成NH码相位捕获. 总之，上述结构占用的运算单元主要有1个复数乘法器、 80个复数加法器，大大降低了电路结构的复杂度和存储器开销. 同时实现了并行化和流水线化结构，加快算法执行速度，无损地保留了算法的高抗干扰性能.

图6 算法流程图Fig.6 Flow diagram of algorithm

3 并行搜星结构的FPGA实现

北斗导航系统采用2000中国大地坐标系(CGCS2000). 该坐标系以地球的质心为原点，Z轴指向国际地球自转服务组织(IERS)定义的参考极(IRP)方向，X轴为IERS定义的参考子午面(IRM)与通过原点且正交于Z轴的赤道面的交线.Y轴与Z轴、X轴构成右手直角坐标系. 导航系统中，卫星按照特定的星历运行，并定时进行校正. 因此，在特定时间卫星所处三维坐标值是明确的，作为接收机定位的参考点. 理论上，接收机锁定三个参考点的信号，分别测量三个参考点的距离，建立方程组，即可换算出接收机的三维坐标值，达到导航定位的目的.

图7 并行化搜星电路框图Fig.7 Parallel circuit diagram of satellite search

北斗导航系统采用的无线电信号测距方法称为单向测距法. 接收机通过解调出卫星信号中包含的系统时间信息，与接收终端的本地时间进行对比，得到该电文在空间中传输的时长. 该时长乘以无线电的传播速度，即可换算出卫星与接收机的距离. 所有北斗卫星的星上时间定期与地面控制站进行校对，可以达到极高的一致性. 但是，由于接收机时钟精度和成本原因，接收机无法与导航系统严格对时，与系统时间存在一定的误差，称为“钟差”. 因此所测得的参考点距离并不是真实距离，而是包含了钟差项，称为“伪距”. 接收机测量锁定四路卫星信号，分别测量四个参考点的伪距，建立四个方程，通过解方程组消除钟差影响，求得三维坐标.

导航系统需要锁定4路以上的卫星信号才能完成定位功能. 在可见卫星数量较多时，搜索更多卫星信号有助于提高定位精度. 定位终端的高速位移或者周围环境的干扰，已锁定信号失锁后可能导致定位精度的下降. 此时须快速锁定其他替补信号，避免定位服务中断. 因此，捕获模块需要保持激活状态，不断搜索新的可用的卫星信号. NH码相位检测电路简单，适合实例化多个电路，同时搜索多路卫星信号，实现快速捕获.

并行化搜星电路框图如图7所示. 20 ms长的信号序列数据量巨大，各检测单元共用信号缓存可以节约大量的存储空间. 检测单元1至n同步执行，同一时刻所需的接收序列数据相同，实现信号缓存共用. 数字基带信号是多路卫星信号的叠加，每个卫星相对速度不同导致频率偏移各不相同，不同的检测单元需要补偿不同频偏. 各检测单元设置了独立的频偏补偿电路，各自独立运行，补偿不同卫星信号的频偏. 频偏补偿控制电路根据信号缓存数据输出顺序，补偿数据频偏. 各检测单元本地测距码发生器独立运行，根据捕获目标编号产生对应的测距码.

4 算法性能仿真

如表1所示，假设采样率等于测距码码率，对直接计算Z(n)[13]和分部相干算法计算量进行对比. 滑动互相关次数设置为(10 230×20)次，保证Z(n)必定出现相关峰. 对比结果表明，改进算法的计算量降低了一个数量级.

表1 算法计算量对比

不同信噪比(SRN)和不同频偏环境下，对分部相关算法性能进行仿真. 如图8 所示，仿真了信噪比-34 dB环境下，相关峰模值与频偏大小的关系. 横轴表示接收序列频偏(fd)，单位Hz，纵轴表示相关峰的模值，单位为1. 仿真曲线与式(7)变化趋势基本一致. 由于高斯噪声的影响，导致仿真曲线有毛刺存在. 频偏小于±30 Hz时，相关峰仍然明显. 如图9所示，对分部相关算法的捕获概率进行仿真. 采样率设置为10 230 MHz，每个点进行一万次仿真，对比同步结果与真实值换算得到捕获概率. 横轴表示信噪比，单位dB，纵轴表示捕获概率，单位为1. 随着信噪比的下降，捕获概率呈下降趋势. 信噪比大于-34 dB时，基本保持100%的捕获概率；信噪比在-34 dB至-37 dB之间，频偏为30 Hz的曲线率先出现捕获失败，频偏0 Hz捕获概率基本保持100%；信噪比小于-37 dB时，两条曲线均出现捕获失败，频偏越低捕获概率越高，可见剩余频偏一定程度上降低了码流的抗噪声性能. 因此，频偏控制在±30 Hz以内，信噪比在-34 dB以上，分部相关算法能够保证可靠捕获.

图8 分部相关算法抗频偏性能图 图9 分部相关算法捕获概率

Fig.8 Anti frequency offset performance of correlation algorithm Fig.9 Acquisition probability of correlation algorithm

利用Quartus Prime V18.1软件实现了捕获算法，并下载到开发板中进行硬件仿真. 采用MATLAB产生信噪比-30 dB、 频偏0 Hz、 时长40 ms的接收序列，写入FPGA开发板中作为测试信号. 为了验证捕获所需最大时钟周期个数，将NH编码起始位置设置在204 599位置. 硬件仿真时序图如图10所示，Clk表示时钟信号，运算单元跟随时钟信号的步调执行运算. 时钟频率越高，则分部相关算法执行速度越快. 芯片外部采用50 MHz的晶振，通过FPGA芯片内部PLL倍频至500 MHz[14]. Position用于输出相关峰位置，PeakValue用于输出相关峰绝对值. Ena作为输出信号使能端，输出高电平时Position和PeakValue两端口输出数据有效. ClkCount作为时钟周期计数单元，用于验证算法执行所需的时钟周期个数. ClkCount等于104 653 453时，Ena由低电平跳变为高电平，捕获算法结束，得到相关峰峰值为204 694，相关峰位置204 599，与测试数据参数一致，捕获成功.

图10 FPGA仿真时序图Fig.10 Sequence diagram of FPGA simulation

5 结语

根据BDS-B3I 频段D1导航电文的特点，理论演绎了累加时长20 ms互相关运算用于信号捕获的可行性. 针对其计算量大，捕获慢的缺点，设计了分部相关算法，保持原有性能的同时减少90%计算量. 使用FPGA设计了一个并行处理且流水线化的结构，利用加法器这类开销较小的运算器单元实现了该算法. 多个算法单元实例可共用信号缓存并行搜星. MATLAB仿真和FPGA硬件仿真表明，本研究算法能够在频偏小于±30 Hz、 信噪比大于-34 dB环境下可靠完成信号同步，大约10.5×107个时钟周期即可完成捕获. 分部相关算法对跟踪单元的位同步和频偏跟踪单元的设计具有一定借鉴意义，为将来提供重要的研究方向.