群时延傅里叶分解模型及其估计方法*

肖志斌，黄仰博，唐小妹，刘增军，倪少杰

(国防科技大学 电子科学学院， 湖南 长沙 410073)

大部分卫星导航系统采用了直接序列扩频系统，如美国的GPS、中国的BeiDou、俄罗斯的GLONASS以及欧洲的Galileo等。直接序列扩频系统通过扩展频谱提高信号增益，而扩展频谱后，接收信道带宽将增大，此时接收通道很难保证理想的幅频及群时延特性，导航信号会产生畸变，进而产生测距偏差[1-2]，已成为高精度测距的主要误差源之一。随着导航系统的现代化，大量现代化导航信号采用了二进制偏移载波(Binary Offset Carrier, BOC)[3]调制技术，部分高阶BOC信号带宽很宽，如AltBOC(15,10)[4]、BOC(15,2.5)[5]等，其受幅频及群时延非理想特性的影响更明显。

群时延非理想导致信号畸变，进而产生测距偏差，因此有必要对其进行测量与校准。目前的群时延测量技术主要有两类：矢量网络分析仪法[6]以及基于信号自相关函数频谱分析的群时延估计方法。矢量网络分析仪法主要通过测量通道的相位特性进而得到通道的群时延特性，矢量网络分析仪法会存在测量精度与测量分辨率无法同时兼顾的问题[6]。另外，该方法的实施复杂，需要采用一系列外围设备，如扫频源、矢量网络分析仪等，硬件链路复杂，测试流程也复杂。且该方法一般适合测量输入输出频段一致的单个通道的群时延特性，对于导航接收机，会存在多个变频环节，整机存在多个变频器及多个滤波器，因此测量整机的群时延特性将更为繁杂。

通道群时延特性可以为关于频率的任意连续函数，可以采用函数级数对任意群时延特性进行描述，这样任意群时延特性可用一组函数级数的系数表示，群时延的测量可通过估计这些函数级数的系数来实现，此时预先计算并存储各种群时延分解系数组合下的相关峰波形数据，通过测量信号相关峰并进行相关峰匹配搜索即可对这些系数进行估计。

该方法的可行性主要取决于函数级数的选取，目前比较经典的群时延分解模型是基于泰勒级数展开的分解模型[7]，该模型将任意群时延特性分解为一系列多次项群时延，如1阶群时延、2阶群时延(抛物线群时延)、3阶群时延等，上述模型可以描述实际中的大部分滤波器，比如巴特沃兹、切比雪夫滤波器均表现为线性群时延与二阶群时延的叠加，但该模型很难描述群时延具有三角函数特性的滤波器，如声表面波滤波器。另外该模型分解后的各阶群时延仍然没有明确的物理意义，只是在数学意义上的曲线，不能和相关峰畸变建立直接的联系，因此无法确定函数分解的阶数。

本文通过理论推导分析了三角函数群时延对信号相关峰的影响，发现三角函数群时延对于信号相关峰的影响是产生了一系列多径信号，且三角函数频率越高，多径距离主瓣信号越远，幅度也越低，影响也越小，高频的三角函数群时延对信号相关峰的影响可忽略不计。基于此项分析，可以将任意群时延特性进行傅里叶分解，将其分解为有限阶数余弦型与正弦型群时延的级联，经过本文分析，采用3阶余弦型群时延和正弦型群时延即可表征大部分通道群时延对信号相关峰的影响，这样只需要估计这3阶余弦型群时延和正弦型群时延的系数即可，此时预先计算并存储各种群时延傅里叶分解系数组合下的相关峰波形数据，通过测量信号相关峰并进行相关峰匹配搜索即可对这些系数进行估计。

1 三角函数群时延对信号相关峰的影响

三角函数群时延主要分为余弦型群时延和正弦型群时延两类，余弦型群时延及正弦型群时延可分别表示为：

(1)

(2)

其中，a、b分别为余弦型群时延和正弦型群时延的波动幅度,fp为群时延波动频率。余弦型群时延传输函数可表示为：

Hc(f)=exp[jθcsin(2πTf)]

=cos[θcsin(2πTf)]+sin[θcsin(2πTf)]

(3)

其中，θc=-Ba,T=1/B，利用第一类贝塞尔函数[8-10]可得到：

(4)

根据余弦型群时延传输函数分解的结果，正弦型群时延传输函数满足：

=exp[jθssin(2πTf)]

(5)

即：

(6)

显然，余弦型群时延τc(f)对信号相关峰的影响是产生了一系列延迟量为kT(k!=0)，幅度为|J-k(θ)| (k!=0)的回波信号；正弦型群时延τs(f)信号相关峰的影响是产生了一系列延迟量为kT(k!=0)，幅度为|J-k(θ)| (k!=0)，相位为kπ/2的回波信号。

2 基于傅里叶级数的通道群时延分解模型

任意通道的传输函数可表示为：

H(f)=A(f)exp[jψ(f)]

(7)

其中，A(f)为幅频响应，ψ(f)为相频响应。群时延的表达式为：

(8)

假设通道的双边带带宽为Bf，对τ(f)在通道带宽范围内进行傅里叶分解，可得到：

(9)

其中，ci与si为傅里叶系数。

根据前一节的分析结果，第i级余弦型群时延与正弦型群时延的传输函数可分别分解为：

(10)

(11)

对于BPSK(n)信号，其码率为fc=1.023nMHz，对应的码率为Tc=1/fc, 通道带宽B取为信号双边带带宽，即B=2fc。显然，当i/B≥2Tc，即i≥2B/fc=4时，第i级余弦型/正弦型群时延将产生一系列与主相关峰相距超过2倍码片的回波信号，且幅度很小，而主相关峰主瓣宽度为两个码片，因此这些回波信号不会对主相关峰产生影响。图1以码率fc=10.23 MHz的扩频信号为例，给出了第4级余弦群时延τ(f)=c·cos(4πf/fc)对信号相关峰的影响，对应群时延波动周期为fc/2=5.115 MHz，c为100 ns。显然，该群时延仅产生了距离主相关峰两倍码片的回波信号，不会影响主相关峰。

图1 余弦型群时延对BPSK(10)信号的相关峰的影响Fig.1 Influence of cosine group delay on the correlation peak of BPSK(10)

虽然第4级(包括第4级)以上的正弦或余弦型群时延产生的回波与信号主相关峰的距离均大于两倍码片，不会对信号的主相关峰产生影响，但是多个余弦或正弦群时延的级联仍然可能产生靠近主相关峰的回波。如第4和第5级余弦/正弦型群时延产生的回波相对于主相关峰的延迟分别为：

dτ4(k)=2kTc

(12)

dτ5(k)=2.5kTc

(13)

显然，第4级余弦型/正弦型群时延产生的第k=-1时的回波与第5级余弦型/正弦型群时延产生的第k=1时的回波进行级联会产生与主相关峰距离为Tc/2的回波，该回波会叠加到主相关峰上，导致主相关峰畸变。

尽管4级(包括第4级)以上余弦型群时延或正弦型群时延的级联会产生接近主相关峰的回波，进而导致主相关峰畸变。但通常情况下，通道滤波器的群时延特性波动较为平缓，具有低通特性(如切比雪夫、巴特沃兹滤波器等)，即对其傅里叶分解后，i越大，群时延波动幅度ci和si越小，相应产生的回波信号幅度也越小，当i足够大时，回波信号的幅度可忽略。例如，对于10.23 MHz码率的扩频信号，若通道滤波器带宽取为20.46 MHz，群时延特性是波动幅度为150 ns的二次曲线。对该通道群时延在信号带宽范围内进行傅里叶分解后得到的各级余弦型群时延的波动幅度以及产生的第1对回波与主相关峰的幅度比如表1所示(通常情况下，根据贝塞尔函数的特性，|k|越大，幅度越小，因此这里仅考虑第1对回波的幅度)。由表1可知，第4级以后的余弦型群时延产生的回波信号幅度相对于主相关峰低37.7 dB以上，基本可以忽略。

表1 带宽为20.46 MHz，波动幅度为150 ns的二次群时延傅里叶分解后的前6级余弦型群时延的特性Tab.1 The first 6-lever cosine type group delay characteristics after Fourier decomposition of the parabola group delay with a bandwidth of 20.46 MHz and a fluctuation amplitude of 150 ns

上述例子采用了波动幅度为150 ns的抛物线群时延，这种群时延特性的通道滤波器已经非常恶劣了，大部分通道的群时延波动幅度要小于150 ns。因此通常情况下，对于二进制相移键控(Binary Phase Shift Keying, BPSK)信号，若滤波带宽取为信号带宽，可对其在信号带宽内进行傅里叶分解，只需取前3级三角函数群时延即可等效表征该滤波器的群时延特性对信号相关峰的影响。

3 基于傅里叶级数模型的群时延测量方法

根据上述两节的分析，三角函数群时延对于信号相关峰的影响是产生了一系列多径信号，且三角函数频率越高，多径距离主瓣信号越远，幅度也越低，影响也越小，高频的三角函数群时延对信号相关峰的影响可忽略不计，大部分通道的群时延特性采用3级余弦型群时延和正弦型群时延即可等效表征，这样只需要估计各级余弦型和正弦型群时延傅里叶分解系数ci和si即可完成通道等效群时延的估计。这里采用相关峰匹配的方法估计群时延的傅里叶分解系数，具体实现框图如图2所示，采用离线分析计算的方式预先计算好各种群时延傅里叶系数组合下的相关峰数据并存储至本地，然后对经过通道的导航信号进行相关峰测量，利用测量得到的相关峰在本地相关峰数据库中进行搜索并计算相关峰波形匹配度，波形匹配误差最小时对应的系数ci、si即为傅里叶分解系数估计值，相应群时延估计值为：

(14)

图2 基于相关峰匹配搜索的群时延测量框图Fig.2 Block diagram of group delay measurement based on correlation peak matching search

3.1 离线预处理

采用离线分析计算的方式预先计算好各种群时延傅里叶系数组合下的相关峰数据并存储至本地。

根据线性系统理论，经过相位响应为ψ′(f)的通道的扩频信号相关峰可写为：

(15)

其中，Gr(f)为扩频信号功率谱。相应通道群时延特性为：

(16)

按一定步进遍历一定幅度范围内的傅里叶系数，并计算对应的群时延特性τ(f)与对应的相位响应ψ(f)，然后根据式(15)计算相关峰曲线。显然，群时延的各级傅里叶分解系数遍历步进越小，对任意群时延遍历越精细，相关峰采样间隔dt越小，相关峰曲线越精细，通过一阶线性插值可获取各相关峰采样点之间的相关值。

3.2 群时延测量流程

对经过通道的导航信号进行相关峰测量，利用测量得到的相关峰在本地相关峰数据库进行匹配搜索，搜索匹配误差最小的相关峰，具体方法如下：

1)信号相关峰测量。令输入的理想基带信号时域波形为q(t)，其经过通道滤波器后的信号为q′(t)。首先采用经典的早迟环估计器估计q′(t)的时延τ0(早迟间隔取为d0)，对计算输入信号的相关峰如式(17)所示：

(17)

由于实现资源受限，只能得到一定延迟间隔d0τ下的相关峰采样值，相关峰采样范围为正负一个码片。令扩频码片周期为T，延迟间隔d0τ=T/L，则相关峰采样值如式(18)所示(已进行归一化，相关峰最大值为1)：

(18)

2)相关峰匹配搜索。将测量得到的相关峰与本地存储的精细相关峰数据进行匹配误差计算，搜索匹配误差最小的群时延傅里叶系数组合，该组系数即表征了群时延估计值。具体计算流程如下所示(对于任意一组群时延系数组合，假设其离线计算并存储的相关峰为Rlocal(τ)(已进行归一化)，测量得到的相关峰采样值为R(k))：

①首先采用传统早迟环估计器来估计某组本地相关峰的中心位置τ′0，其满足(其中d0为早迟间隔)：

Rlocal(τ′0-d0/2)=Rlocal(τ′0+d0/2)

(19)

②对本地存储的相关峰数据Rlocal(τ)进行采样，各采样点对应的延迟序列即为相关峰R(k)各采样点对应的延迟量，采样后得到的本地相关峰序列如式(20)所示：

R′local(k)=Rlocal(τ′0+kT/L-T)

(20)

③计算相关峰匹配误差:

(21)

4 仿真验证

采用码率为10.23 MHz的扩频码信号对本文提出的群时延特性估计方法进行仿真验证，滤波器带宽为20.46 MHz。滤波器的群时延采用3阶余弦型/正弦型群时延等效表示。

仿真验证框图如图3所示，给定某通道滤波器H(f),首先采用本文所提出的基于傅里叶级数分解的群时延估计方法对该滤波器的群时延特性进行测量，得到各级傅里叶分解系数估计值c1,c2,c3，s1,s2,s3，具体步骤如下所示。

步骤1：计算各类傅里叶分解系数组合下的相关峰，并存储至本地，各级傅里叶分解系数的遍历步进取为0.1 ns，遍历范围取为-60～60 ns。

步骤2：计算给定传输函数H(f)通道的滤波器脉冲响应h(t)，生成理想基带信号r(n)，并将其通过脉冲响应为h(t)的滤波器，输出为rI(n)+jrQ(n)。

步骤3：对rI(n)+jrQ(n)进行相关峰测量，并在本地相关峰数据库中进行匹配搜索。

步骤4：相关峰匹配误差最小的系数组合即为群时延傅里叶级数的估计结果c1,c2,c3，s1,s2,s3。

随后根据该估计结果修正原始传输函数H(f)，假设H(f)的群时延特性为τ(f)，则修正后的传输函数H′(f)的群时延响应如式(22)所示：

(22)

采用将理想基带信号通过修正后的滤波器H′(f)，来计算精细相关峰以及不同相关间隔下的时延估计值，并与理想相关峰进行比对。

图3 仿真验证框图Fig.3 Simulation and verification block diagram

仿真采用两个滤波器：一个为给定系数的巴特沃兹滤波器，另一个为某导航接收机实际使用的中频滤波器(其群时延特性采用矢量网络分析仪测量得到)，两个滤波器的群时延特性分别如图4和图5所示。

图4 巴特沃兹仿真滤波器的群时延特性Fig.4 Group delay characteristics of simulation Butterworth filter

图5 实际滤波器的群时延特性Fig.5 Group delay characteristics of actually measured filter

仿真结果如图6～11所示，显然理想信号经过上述两个滤波器后均发生了不同程度的畸变(信号相关峰畸变，不同早迟间隔下的时延估计偏差互不相同)，这主要是由于群时延非理想导致。采用本文方法估计得到的群时延特性对滤波器进行群时延修正，显然经过校准后，相关峰畸变得到修正，其与通过滤波器之前的信号相关峰高度吻合，差异幅度很小(主瓣范围内差异约为主峰幅度的1/10 000量级)，同时不同早迟相关间隔的时延估计偏差在0.1 ns之内。仿真结果说明了本文所提出的群时延估计方法的正确性。

注：1—理想信号相关峰;2—经过滤波器后的相关峰; 3—采用群时延估计结果对滤波器校准后的相关峰。图6 通过实际滤波器之后的信号相关峰、群时延校准后的相关峰与理想相关峰的对比Fig.6 Comparison of the correlation peaks after passing the actual filter, the correlation peak after group delay calibration and the ideal correlation peaks

注：1—理想信号相关峰;2—经过滤波器后的相关峰; 3—采用群时延估计结果对滤波器校准后的相关峰。图7 通过巴特沃兹滤波器之后的信号相关峰、群时延校准后的相关峰与理想相关峰的对比Fig.7 Comparison of the correlation peaks after passing the Butterworth filter, the correlation peak after group delay calibration and the ideal correlation peaks

图8 对通过实际滤波器的信号进行群时延校准后的信号相关峰与理想相关峰的偏差Fig.8 Bias between the ideal correlation peak and the correlation peak after group delay calibration of the signal passing through the actual filter

图9 对通过巴特沃兹滤波器的信号进行群时延校准后的信号相关峰与理想相关峰偏差Fig.9 Bias between the ideal correlation peak and the correlation peak after group delay calibration of the signal passing through the Butterworth filter

图10 对通过实际滤波器信号进行群时延校准后，不同相关间隔下的时延估计偏差Fig.10 After the group delay calibration of the signal passing through the actual filter, the delay estimation bias under different correlation intervals

图11 对通过巴特沃兹滤波器信号进行群时延校准后，不同相关间隔下的时延估计偏差Fig.11 After the group delay calibration of the signal passing through the Butterworth filter, the delay estimation bias under different correlation intervals

5 实测验证

采用湖南卫导公司的NSS8000系列导航信号源与实际GNSS接收机对本文提出的群时延特性估计方法进行实测验证，导航信号源播发单颗卫星高载噪比(大于50 dB-Hz)的B2a信号(尽可能消除互相关和测量噪声的影响)，GNSS接收机B2a信号数字接收通道的线性相位有限长单位冲激响应(Finite Impulse Response, FIR)滤波器带宽为25 MHz(接收机设计参数)。

具体的测试步骤如下所示。

步骤1：GNSS接收机对导航信号源发送过来的B2a信号进行接收处理，并测量相关峰(相关峰积累时长为1 s)，目前该接收机相关峰采集点数为32，间隔为1/16码片。

步骤2：计算各类傅里叶分解系数组合下的相关峰(相关峰滤波带宽为25 MHz，具体滤波器取为接收机B2a信号数字接收通道的线性相位FIR滤波器)，并存储至本地，各级傅里叶分解系数的遍历步进取为0.1 ns，遍历范围取为-60～60 ns。

步骤3：根据GNSS接收机测量得到的相关峰，按照本文方法在本地相关峰数据库中进行匹配搜索，相关峰匹配误差最小的系数组合即为群时延傅里叶级数的估计结果c1,c2,c3,s1,s2,s3，对应的通道群时延特性估计值如下(其中带宽B为25 MHz)：

(23)

步骤5：群时延补偿后，GNSS接收机再次进行相关峰测量，比较补偿前后的相关峰测量结果，并计算补偿前后不同早迟间隔下的测距偏差。

测试结果如图12～14所示,其中图12与图13为GNSS接收机监控软件所存储相关峰测量数据的图示，图14为群时延校准前后不同相关间隔下的时延估计值，显然，补偿前信号相关峰存在畸变，不同早迟间隔下的时延估计偏差高达1.7 ns，经群时延估计补偿后，相关峰畸变得到修正，不同早迟相关间隔的时延估计偏差在0.1 ns内。实测结果验证了本文所提出的群时延估计方法的正确性。

图12 群时延补偿之前的B2a信号相关峰Fig.12 B2a correlation peak before the group delay compensation

图13 群时延补偿后的B2a信号相关峰Fig.13 B2a correlation peak after the group delay compensation

图14 群时延校准前后的不同相关间隔下的时延估计偏差Fig.14 Delay estimation bias in different intervals before and after calibration of group delay estimation

上述测量的是导航信号源和接收机级联的群时延特性，导航信号源为标准设备，通常需给出其群时延特性，测量结果减去导航信号源的群时延特性即为GNSS接收机的群时延特性。

6 结论

本文通过理论推导分析了三角函数群时延对信号相关峰的影响，发现三角函数群时延对于信号相关峰的影响是产生了一系列多径信号，且三角函数频率越高，多径距离主瓣信号越远，幅度也越低，影响也越小，高频的三角函数群时延对信号相关峰的影响可忽略不计，从而可将任意群时延特性分解为有限阶数的三角函数群时延的叠加。基于此分析，通过测量信号相关峰，并与不同三角函数群时延组合下的相关峰进行匹配搜索即可估计群时延的三角函数分解级数，从而估计得到相关峰影响等效的群时延特性，仿真结果验证了本文方法的准确性。相对于传统基于矢量网络分析仪的群时延测量方法，本文方法可实现软件自动化测试，且不需要搭建额外复杂的硬件链路，测量效率更高，具有较强的工程应用价值。