潘毅，张天骐，张刚，马宝泽

（重庆邮电大学通信与信息工程学院，重庆 400065）

1 引言

全球导航卫星系统（GNSS, global navigation satellite system）能够提供全天候、连续且实时的高精度地理位置信息及导航、授时信息、短分组通信等服务，其中北斗导航定位系统（BDS, Beidou system）是我国独立自主研发建设的[1-2]。目前，北斗导航定位系统除了采用传统的二进制相移键控（BPSK, binary phase shift keying）信号外，提供的L1导航信号等还采用了能够实现频谱共享及跟踪性能更优越的新信号制式二进制偏移载波（BOC,binary offset carrier）信号[3-4]。

实际的通信环境日益复杂，飞行器速度快、加速度大，甚至加加速度大，接收机往往是在高动态环境下工作，使接收到的信号具有大多普勒频偏和一阶、二阶多普勒变化率。这对BOC信号的捕获提出了较大的挑战，而大部分已经提出的捕获算法[5-7]只能在载波频率完全消除的前提下使用，不能对加速度和加加速度分别引起的一阶、二阶多普勒变化率分量进行有效的补偿。为了解决高动态下的捕获问题，文献[8]采用部分匹配滤波器（PMF, partial matched filter）和快速傅里叶变换（FFT, fast Fourier transform）相结合的算法，提高了多普勒频偏估计的范围和分辨率，但是无法对多普勒变化率进行补偿。文献[9]通过离散chirp傅里叶变换对多普勒频偏和一阶多普勒变化率进行估计，但在高速径向运动中效果不尽如人意。文献[10]采用最大似然算法来估计码相位、多普勒频偏和载波相位误差，但同时也引入了巨大的计算量，难以在要求实时性、结构简易的接收机中使用。当伪码对齐时，扩频信号的相关结果可以看作线性调频（LFM, linear frequency modulation）信号，文献[11-12]引入了分数阶傅里叶变换（FRFT, fractional Fourier transform）算法来解决高动态GNSS信号的捕获问题。

目前，几乎所有的高动态捕获算法都是针对大多普勒频偏或存在一阶多普勒变化率的情况，因此本文采用离散多项式相位变换（DPT, discrete polynomial-phase transform）算法来处理存在二阶多普勒变化率的BOC信号，并通过先对BOC信号定阶再根据信号动态阶数分别进行处理，最后通过FRFT算法进行捕获，最终形成了一套完整的高动态BOC信号捕获结构。

2 高动态BOC信号模型

接收到的高动态BOC信号可以表示为

其中，A是信号幅度；d(t)是导航电文；s(t)=c(t)sc(t)是基带信号，c(t)是伪码（PN, pseudo-noise）信号，sc(t)=sgn[sin(2πfsct)]或 sc(t)= sgn[cos(2πfsct)]是副载波信号，fsc是副载波频率，φ0是初始相位，n(t)是均值为 0、方差为σ2的加性高斯白噪声。载波信号φ(t)可以表示为

其中，fi是中频频率，fd是由接收机和卫星之间做径向运动产生的多普勒频偏，kd是由加速度引起的一阶多普勒变化率。

经过相关和积分处理后，载波信号φ(t)的输出可以表示为

其中，T表示相干积分时间，Δφ0表示载波相位误差，R(·)表示BOC信号的自相关函数，Δτ表示码相位时延误差。假设多普勒变化率在估计过程中不变，则可以得到一个平均角速度为

那么，式(3)可以改写为[13]

其中，Δτ可以表示为

其中，fr表示L1导航信号的射频频率。

3 捕获算法原理

3.1 FRFT算法

时频域(t,f)经过旋转角α转换到了变换域(u,v)，信号在这个域的表示则由FRFT给出，则信号x(t)的FRFT定义为

其中，k为整数，为旋转角度，p为FRFT的阶数，Fp[·]为FRFT的算子符号。最初的时频域(t,f)经过旋转角α转换到了变换域(u,v)，能够使原本在(t,f)平面发散的频谱在(u,v)平面得到聚集，如图 1所示，从而能够估计出多普勒频偏fd和多普勒一阶变化率kd。

图1 FRFT：(t,f)平面变换到(u,v)平面

那么S(n)在p阶分数阶傅里叶域上的输出为

其中，S(n)是一个带噪的线性调频信号，ω=UΔt，Δt是变换域的采样间隔。离散信号S(n)关于α的N点离散FRFT可以表示为

分数阶傅里叶域的能量峰值所在位置可以表示为

其中，是关于加速度的FRFT旋转角度的估计值，是关于频率误差的FRFT输出阶数估计。多普勒频偏估计值多普勒一阶变化率估计值和归一化FRFT输出分别为[11]

其中，是旋转角度估计值那么，关于多普勒频偏误差的FRFT搜索范围可以表示为

3.2 二阶多普勒变化率

r(n)的M阶瞬态矩定义为

其中，其中，τ是时延长度。那么，r(n)的M阶DPT可以表示为[14]

其中，N1为采样长度，Ts为采样间隔，ω为角频率。M阶多项式相位信号可以表示为

其中，ak为实相位系数。

将式(18)代入式(17)，可以得到

当载波存在二阶多普勒变化率时，采样后的BOC信号r(n)可以表示为

其中，kf表示加加速度引起的二阶多普勒变化率。通常情况下，接收机和卫星径向加速度为-100g～100g，所引起的最大一阶多普勒变化率为±5.15 kHz/s，而加加速度的范围是-40g/s～40g/s，其中，g为重力加速度。即最大二阶多普勒变化率大约为±2.06 kHz/s2。当kf≠0时，FRFT算法不适用于具有加加速度的信号，所以本文采用瞬时互相关的方法进行降阶。由瞬态矩性质可得，P2[s(n),τ]可认为是原 BOC基带信号s(n)经过瞬时互相关运算变成新BOC基带信号s′(n)，因此可以得到

将接收到的信号首先经过一个平方环处理，消除数据位和伪码调制的影响，只剩下含有高动态项的载波信号，再对其进行DPT处理。通过式(19)估计其阶数，得到的二阶DPT幅值谱波形具有如下特征[15]。

1) 如果二阶DPT幅值谱是具有一定带宽的凌乱谱线，如图2(a)所示，说明载波信号是三阶多项式相位信号并含有二阶多普勒变化率，则需要对信号先进行降阶处理再采用FRFT估计其动态参数。

2) 如果二阶DPT幅值谱是非0频的单线谱或者只有0频的谱线，如图2(b)或图2(c)所示，说明载波信号分别是二阶或一阶多项式相位信号，并不含有二阶多普勒变化率，可以直接采用FRFT算法对接收到的信号参数进行估计（如果信号只含有多普勒频偏时，则FRFT算法的阶数为1，变为传统的傅里叶变换）。

图2 二阶DPT幅值谱波形

4 算法性能和算法流程

4.1 降阶处理和FRFT的信噪比

通过瞬时互相关的方法对含有高动态项的信号进行降阶处理，根据式(1)和式(14)进行一次瞬时互相关后，信号和噪声的输出功率可以分别表示为

其中，L为采样后信号的长度。结合式(24)和式(25)可得输出信噪比为

根据FRFT的Parseval关系，其FRFT域的检测信噪比可以表示为

其中，var表示方差，S(α,U)表示信号在FRFT域的幅值，W(α,U)表示噪声在 FRFT的幅值。FRFT算法能够使降阶处理后的信号在最佳旋转角度的变换域中能量得到聚集，即var[S(α,U)]取得最大值。但对于高斯白噪声在任意阶次的分数阶傅里叶域内能量分布是均匀的，即var[W(α,U)]是不变的，所以在估计动态参数的同时抑制了噪声并提高了信号的检测信噪比，增强了系统的抗噪性能，更易于实现高动态信号的快速捕获。另一方面，对于存在二阶多普勒变化率的BOC信号，所提算法只需要进行一次瞬时互相关就可以达到降阶的目的，而传统的采用DPT和PMF-FFT相结合的捕获算法（以下简称为DPF算法）[17]需要进行2次瞬时互相关才可以进行捕获。可见，在采用相同的延时长度的情况下，所提算法的输出信噪比损失更小，降阶计算量更小。

4.2 检测概率

由于载波信号中的多普勒变化率和多普勒频偏能通过FRFT精确地被估计，信号的高动态项不会影响FRFT输出峰值的检测，因此基于FRFT捕获算法的检测概率和基于 FFT捕获算法的检测概率是相同的。接收信号和本地信号相关后的输出服从 Rice分布，其虚警概率Pfa和检测概率Pd分别表示为

其中，Pn(x)是只含有噪声时的概率密度函数，Ps(x)是含有有用信号和噪声时的概率密度函数，I0是 0阶修正贝塞尔函数，V是门限值，可以表示为

非中心参数λ可以写成

其中，Tcoh是积分间隔。对于一个给定的系统虚警概率Pfa，检测概率为

其中，Q(x,y)是MarcumQ函数。平均捕获时间也是捕获算法一个重要的指标，可以表示为[18]

其中，K是惩罚因子，q是搜索单元数，τd是整个相关所需的积分时间。

4.3 所提算法流程

所提算法步骤如下。

步骤1接收端将接收到的高动态BOC信号下变频，经过模数转换器（ADC, analog-to-digital converter），以伪码速率进行采样，得到离散信号r。

步骤 2采样后的信号通过本地载波数字控制振荡器（NCO, numerically controlled oscillator）进行载波剥离。

步骤3对仍存在残余频偏的信号进行定阶运算，根据阶数确定是否进行降阶处理，得到的信号为s。

步骤4对本地产生的PN序列进行副载波调制，并根据阶数确定是否对本地 BOC序列进行降阶处理，然后与s进行相关处理。

步骤5对相关结果进行FRFT，根据式(10)进行峰值二维搜索。

步骤6将相关结果的最大值与门限V进行比较（V值由式(30)得到）。如果最大值超过预设门限值V，则说明本地PN序列与接收BOC信号的伪码相位对齐，进入跟踪阶段（本文主要解决捕获问题，跟踪过程不予以讨论）；如果没有超过预设门限值V，则通过逻辑控制本地载波 NCO，根据式(11)得到的动态参数对多普勒频偏进行补偿，同时滑动本地PN相位。

步骤7重复步骤3～步骤6，直到伪码相位对齐。如此就形成了一整套捕获方案，如图3所示。

5 仿真实验

为了检测所提捕获算法的性能，本文采用Matlab对 sinBOC(1,1)信号进行仿真实验，其中信号的采样频率为16.368 MHz，码速率为1.023 MHz，码片长度为1 024。FRFT运算时间为51.2 ms。仿真参数设置为多普勒频率fd=50 kHz，一阶多普勒变化率为 5 kHz/s，二阶多普勒变化率为 2 kHz/s2，SNR=-5 dB。实验中采用DPF算法进行对比仿真，其中PMF时间为16 ms，码片长度取3 069，分为3段，每段长度为1 023。图4为Matlab实验操作界面。

图3 本文所提算法捕获流程

图4 Matlab操作界面

5.1 FRFT输出峰值谱

由于高斯函数的Wigner-Ville分布仍然是高斯函数，高斯白噪声在任何阶数FRFT下仍然是高斯形状。而载波信号能够在相应阶数的FRFT下实现能量聚集，比较易于信号的检测与参数估计。FRFT输出峰值谱如图5所示。由图5可以明显看出，输出谱在p=1.051和u=444处出现了一个峰值。根据式(11)和式(23)可以计算出，多普勒频偏fd=51.32 kHz，多普勒一阶变化率kd=4 930.00 Hz/s，多普勒二阶变化率kf=2 008.98 Hz/s2，这与之前的设定值误差保持在2%左右，因此所提算法能够准确估计高动态BOC信号的动态参数。

图5 FRFT输出峰值谱

5.2 相关函数振幅谱

图6是对存在二阶多普勒变化率的BOC信号分别采用本文所提算法、DPF算法和PMF-FFT算法进行捕获得到的相关函数归一化振幅值。图 6(a)中，所提算法在码相位为596位置的主峰幅值明显高于其他位置；图6(b)中，DPF算法的振幅谱能够在相同的码相位处得到一个最大幅值，但有多个位置的幅值与最大幅值非常接近，这就造成接收机无法正确锁定主峰位置；图6(c)中，PMF-FFT算法输出最大幅值所对应的码相位位置错误。通过图6对比分析，可以证明本文所提算法的有效性，也使判决器更容易检测到主峰位置，且检测概率更高。

图6 相关函数归一化振幅谱

5.3 检测概率

对不同动态参数下的 BOC信号采用本文所提算法、DPF算法和PMF-FFT算法得到的检测概率如图7所示。从图7可以看到，对于存在二阶多普勒变化率的信号，所提算法比 DPF算法提高了3.5 dB；对于只有一阶多普勒变化率的信号，所提算法比 DPF算法提高了 4.5 dB。而传统 PMF-FFT算法在大动态参数条件下已经无法实现信号的捕获。

图7 本文所提算法、DPF算法和PMF-FFT算法的检测概率

5.4 平均捕获时间

捕获时间的参数设置为搜索单元数q=2 048，虚警概率Pfa=0.01，积分时间τd=1 ms，惩罚因子K=1。从图8可以看出，对于存在二阶多普勒变化率的信号，当SNR=-5.5 dB时，本文所提算法的平均捕获时间降到最低，大约是DPF算法平均捕获时间的；对于只有一阶多普勒变化率的信号，当SNR=-13 dB时，所提算法的平均捕获时间降到最低，大约是DPF算法平均捕获时间的。在低信噪比下，所提算法对多普勒频偏、多普勒变化率等高动态参数的估计更加准确，残余频偏较小甚至为0。这使捕获检测函数引入的误差较小，更有利于码相位的搜索。但是随着信噪比的升高，DPF算法的平均捕获时间要低于所提算法，这是由于此时多普勒频偏估计更准确，平均捕获时间主要取决于算法的计算时间，DPF算法估计时间为16 ms，FRFT估计时间为51.2 ms[12]，即DPF算法的计算复杂度要低于FRFT算法。而接收到的BOC信号的信噪比往往都是较低的，所以对于算法来说在低信噪比条件下的性能表现更加重要。

5.5 计算复杂度

FRFT算法和 PMF-FFT算法的计算量可以表示为[13]

图8 本文所提算法和DPF算法的平均捕获时间

其中，pu是 FRFT阶数搜索的上限值，pl是 FRFT阶数搜索的下限值；；τ1和τ2分别是FRFT和PMF-FFT算法全部估计所需要的时间；T1和T2分别是FRFT和PMF-FFT算法单个输出的时间是FRFT和FFT算法输出的比值。而进行一次瞬时互相关运算需要的计算量为

对于存在二阶多普勒变化率的 BOC信号，所提算法需要进行一次瞬时互相关运算，DPF算法需要进行2次瞬时互相关运算。在实际应用中，当时延长度为，长度为L的BOC序列经过一次瞬时互相关后，其序列长度变为，再进行FRFT，此时；第二次时延长度变为，进行第二次瞬时互相关，最后进行 PMF-FFT，此时则所提算法和DPF算法的总计算量为

可见，当K12=1时，即2种算法输出信号长度相同时，所提算法的计算量比DPF算法要高，这是由于DPF算法的硬件需求更高，在增加匹配滤波器个数的情况下使计算量降低。

6 结束语

对存在二阶多普勒变化率的高动态BOC信号，本文引入DPT算法进行先定阶后降阶的处理，然后基于 FRFT算法估计其多普勒频偏和码相位。所提算法的检测概率、平均捕获时间和计算复杂度都进行了分析，并在仿真中与DPF、PMF-FFT算法进行了性能对比。对于存在二阶多普勒变化率的高动态BOC信号，所提捕获算法的检测概率比DPF算法提高了3.5 dB，并在低信噪比情况平均捕获时间较少。本文提出的一整套捕获方案不但解决了传统捕获算法无法估算二阶多普勒变化率的问题，而且提高了系统的抗噪性能，非常适用于高动态BOC信号。