一种卫星信号载波频率精确估计算法

２０１２年２月１日第３５卷第３期

现代电子技术

Ｍｏｄｅｒｎ　Ｅｌｅｃｔｒｏｎｉｃｓ　Ｔｅｃｈｎｉｑ

ｕｅＦｅｂ．２０１２

Ｖｏｌ．３５Ｎｏ．３

一种卫星信号载波频率精确估计算法

黄富彪，何兵哲

（中国空间技术研究院西安分院，陕西西安　７１００００

）摘　要：针对卫星载波频率高精度估计问题，阐述了现有的频率估计实时处理方法，结合类Ｒｉｆｅ频率修正算法，并考虑卫星信号特点，给出了适用于数字方法实现的卫星信号载波频率高精度估计算法流程。对其进行了仿真实现，其结果表明，该算法能有效提高卫星信号载频估计精度。

关键词：载频估计；ＦＦＴ；ＣＺＴ；Ｒｉｆｅ算法

中图分类号：ＴＮ９１１－３４ 文献标识码：Ａ 文章编号：１００４－３７３Ｘ（２０１２）０３－０１２９－

０３Ａ　ｃａｒｒｉｅｒ　ｆｒｅｑｕｅｎｃｙ　

ａｃｃｕｒａｔｅ　ｅｓｔｉｍａｔｉｏｎ　ａｌｇｏｒｉｔｈｍ　ｆｏｒ　ｓａｔｅｌｌｉｔｅ　ｓｉｇｎａｌＨＵＡＮＧ　Ｆｕ－ｂｉａｏ，ＨＥ　Ｂｉｎｇ

－ｚｈｅ（Ｘｉ’ａｎ　Ｂｒａｎｃｈ，Ｃｈｉｎａ　Ａｃａｄｅｍｙ　

ｏｆ　Ｓｐａｃｅ　Ｔｅｃｈｎｏｌｏｇｙ，Ｘｉ’ａｎ　７１００００，Ｃｈｉｎａ）Ａｂｓｔｒａｃｔ：Ａｉｍｉｎｇ　ａｔ　ｔｈｅ　ｓａｔｅｌｌｉｔｅ　ｓｉｇｎａｌ　ｃａｒｒｉｅｒ　ｆｒｅｑｕｅｎｃｙ　ｈｉｇｈ－ｐｒｅｃｉｓｉｏｎ　ｅｓｔｉｍａｔｉｏｎ，ｔｈｅ　ｅｘｉｓｔｉｎｇ　ｒｅａｌ－ｔｉｍｅ　ｆｒｅｑｕｅｎｃｙ　ｅｓｔｉ－ｍａｔｉｏｎ　ｍｅｔｈｏｄｓ　ａｒｅ　ｐｒｅｓｅｎｔｅｄ．Ｃｏｍｂｉｎｇ　ｗｉｔｈ　Ｒｉｆｅ　ｌｉｋｅ　ｆｒｅｑｕｅｎｃｙ　ｍｏｄｉｆｉｃａｔｏｒｙ　ａｌｇｏｒｉｔｈｍ　ａｎｄ　ｃｏｎｓｉｄｅｒｉｎｇ　ｔｈｅ　ｃｈａｒａｃｔｅｒｉｓｔｉｃｓ　ｏｆｓａｔｅｌｌｉｔｅ　ｓｉｇｎａｌ，ｔｈｅ　ａｌｇｏｒｉｔｈｍ　ｆｌｏｗ　ｃｈａｔ　ｗｈｉｃｈ　ｉｓ　ｓｕｉｔａｂｌｅ　ｆｏｒ　ｓａｔｅｌｌｉｔｅ　ｓｉｇｎａｌ　ｃａｒｒｉｅｒ　ｆｒｅｑｕｅｎｃｙ　ｈｉｇｈ－ｐｒｅｃｉｓｉｏｎ　ｅｓｔｉｍａｔｉｏｎ　ｉｓ　ｐｒｏ－ｐｏｓｅｄ．Ｓｉｍｕｌａｔｉｏｎ　ｒｅｓｕｌｔ　ｓｈｏｗｓ　ｔｈａｔ　ｉｔ　ｃａｎ　ｉｍｐｒｏｖｅ　ｔｈｅ　ｆｒｅｑｕｅｎｃｙ　ｅｓｔｉｍａｔｉｏｎ　ａｃｃｕｒａｃｙ　

ｅｆｆｅｃｔｉｖｅｌｙ．Ｋｅｙｗｏｒｄｓ：ｃａｒｒｉｅｒ　ｆｒｅｑｕｅｎｃｙ　

ｅｓｔｉｍａｔｉｏｎ；ＦＦＴ；ＣＺＴ；Ｒｉｆｅ　ａｌｇｏｒｉｔｈｍ收稿日期：２０１１－０８－

２２０　引　言

航天测控任务中，利用飞行器和观测站之间相对运动引起的载频信号多普勒频移可计算出飞行器的径向速度，

从而获知飞行器的运行姿态。相关数据经过进一步处理后，又能用于精确定轨和地面定位，载频估计精

度的要求甚至会达到毫赫兹量级［１］

。因此，对载频的精

确估计有着十分重要的意义。

虽然利用Ｍｕｓｉｃ算法、ＡＲ模型算法以及最大似然估计算法等现代谱估计的方法可以对正弦信号频率进行精确估计，但由于算法复杂，计算量大，因而限制了其在具有实时性要求的卫星信号中的应用。得益于数字信号处理技术的发展，基于离散傅里叶变换（ＤＦＴ）的直接谱估计法在卫星信号处理中得到了广泛的应用。如何在保证实时性的前提下进一步提高载频估计精度已成为研究的重点。１　信号频率估计的基本方法

快速傅里叶变换（ＦＦＴ）是ＤＦＴ的一种快速有效算法，在工程实践中得到了普遍的应用。实际中基于ＤＦＴ的直接谱估计算法也都是采用了快速傅里叶变换

来处理的，以减小运算量，提高处理速度。１．１　常规ＦＦＴ方法

如果对Ｎ点数据进行ＦＦＴ运算，将同时在频域得到Ｎ个结果。如果采样频率是ｆｓ，ｋ０为谱线幅值最大值

点，取其中最大谱线值Ｘ（ｋ０）对应频率＾ｆ０＝ｋ０ｆｓ

／Ｎ，则＾ｆ０为频率估计值，频率分辨率为ｆｓ／Ｎ。显然，当实际信号频率接近ＦＦＴ频率分辨率的整数倍频率点时，估计出的频率性能较好。否则，由于ＦＦＴ“栅栏”效应引起频谱泄露，估计出的误差将会很大。在采样率保持不变时，为了减小频率抽样间隔，提高频率分辨率，只能增加抽样点数，这将使运算量大为增加。ＦＦＴ法可对信号进行实时处理，

是一种最常规的频谱分析手段。１．２　ＣＺＴ法［２

］

常规的ＦＦＴ变换实质是在ｚ平面的单位圆周上对ｚ变换进行Ｎ点等间隔抽样，而线性调频ｚ变换（Ｃｈｉｒｐ－ｚ变换，ＣＺＴ）却可对ｚ平面上的单位圆进行局部抽样，即在有限的频率范围内进行频谱分析。ＣＺＴ的定义式如下，对已知的时间序列ｘ（ｎ）（０≤

ｎ≤Ｎ－１

），其ｚ变换为：Ｘ（ｚ）＝

∑Ｎ－１

ｎ＝０

ｘ（

ｎ）ｚ－

ｎ（１

） 如前所述，

常规的ＦＦＴ变换是在ｚ平面的单位圆周上进行Ｎ点均匀取样，

为使ｚ可在ｚ平面上沿更一

般路径取值，令ｚ的取值样点：

ｚｋ＝ＡＷ－ｋ， ｋ＝０，１，２，…，Ｍ－１（２）式中：Ｍ表示复频谱的点数，取值可以不等于Ｎ；Ａ和Ｗ是任意复数，即：Ａ＝Ａ０ｅｊθ０，Ｗ＝Ｗ０ｅ－ｊφ０，其中，Ａ决定谱分析起始点ｚ０的位置；Ｗ０的值决定谱分析路径的盘旋趋势；φ０表示相邻分析点之间的夹角。如果Ｗ０＜１，表示随着ｋ增大，分析点以φ０为步长向外盘旋，而Ｗ０＞１时，向内盘旋。

对ＣＺＴ的计算过程，一般是将其转化为线性卷积，再利用ＦＦＴ来计算。ＣＺＴ算法的关键在于频谱范围的选定，合适的频谱范围保证了频谱估计的正确性，在相对窄的频谱范围内做同等点数的ＣＺＴ计算时得到的分辨率相对提高，测频精度也就相应提高。从ＣＺＴ的定义来看，一旦确定了频谱估计范围，只要Ｍ值取得足够大，估计的频率值似乎就能达到足够理想的精度。但实际上，它的分辨率还受到数据长度、信号截取方式等的限制［３］。

１．３　ＺＦＦＴ细化方法

基于复调制频移的频谱细化分析法，即Ｚｏｏｍ－ＦＦＴ（ＺＦＦＴ）法也是一种有效的分析方法。ＦＦＴ只能分析从零频开始的一个低通频带，且频带越窄分辨率越高，ＺＦＦＴ的原理就是利用复调制技术将感兴趣的那段频谱移到零频附近，再进行常规的ＦＦＴ运算。相较于ＣＺＴ法，ＺＦＦＴ所需的原始数据会比较长，在某些场合如瞬态过程将不再适用。

２　卫星信号频率估计算法

卫星在轨运动中会存在较大的多普勒频移，因而使载频的估计范围比较大。如果采用常规ＦＦＴ进行载频估计，在ＦＦＴ点数受制于硬件条件，分析点数有限的情况下，估计精度并不高。如果采用ＣＺＴ算法，则由于载频动态范围大，无法恰当地选择频谱估计范围，估计精度也不理想。因而，可考虑先采用常规ＦＦＴ粗估，缩小频谱估计范围，再利用ＣＺＴ进行进一步的细化，以达到精确估计卫星信号载频的目的。

另外，考虑到ＦＦＴ“栅栏”效应，信号的实际频率通常位于ＦＦＴ主瓣内两条最大谱线之间，因而可借助次大谱线与最大谱线的幅度比值来估计信号的实际频率在两条谱线之间的位置，即基于ＦＦＴ幅度比值的频率插值方法，也称Ｒｉｆｅ算法。

下面介绍一种基于ＣＺＴ的改进Ｒｉｆｅ算法［４］。该算法采用类似Ｒｉｆｅ法对ＣＺＴ算法估计出的频率进行进一步的修正。

设待估信号序列为：

ｘ［ｎ］＝Ａ０ｃｏｓ（Ω０ｎ＋φ０）＋ｒ［ｎ］，　０≤ｎ≤Ｎ－１

（３）式中：Ａ０和φ０分别为正弦信号的幅度和初相位；Ω０＝

２πｆ０／ｆｓ，ｆ０为待估信号频率，ｆｓ为信号的采样率；ｎ为采样点数；ｒ［ｎ］为实随机高斯白噪声，方差为δ２ｒ。

ＣＺＴ参数设置如下：Ａ＝ｅｊ　ｆ１２π／ｆｓ，Ｗ＝ｅｘｐ［ｊ（ｆ２－ｆ１）２π／（Ｍｆｓ）］，ｆｗ＝ｆ２－ｆ１。式中：ｆ１＜ｆ０＜ｆ２，即关于待估频率所处频率范围的先验知识。Ｍ为ＣＺＴ变

换的点数，可任意选定。

设ｋ０为谱线幅值最大值点，频率分辨率Δｆ＝

ｆｗ／Ｍ，频率粗估值为ｆ１＋ｋ０Δｆ，则以ｋ０为中心，正弦信号的ＣＺＴ线谱模型为：

Ｘ（ｋ）ＣＺＴ＝

Ａ０

２

Ｎｅｘｐ　ｊφ０－π（ｋ－δ）

Ｎ

Ｍ

ｆｗ

ｆ

［］

｛｝ｓ·

ｓｉｎ［π（ｋ－δ）（Ｎ／Ｍ）（ｆｗ／ｆｓ）］

π（ｋ－δ）（Ｎ／Ｍ）（ｆｗ／ｆｓ）

＋Ｒ（ｋ）ＣＺＴ（４）式中：ｋ为偏离ｋ０谱线的序号值；δ为正弦信号真实频率与粗估频率间的差值大于Δｆ的倍数；δ∈

［－０．５，０．５］；Ｒ（ｋ）

ＣＺＴ

为ｒ［ｎ］的ＣＺＴ谱，也是一个高斯白噪声过程。显然正弦信号频率：

ｆ０＝ｆ１＋（ｋ０＋δ）Δｆ（５） 由线谱模型，仿照Ｒｉｆｅ幅度比插值方法，在不考虑噪声的理想情况下，由泰勒公式，略去高次项可得：

Ａ±１＝Ｘ（±１）ＣＺＴ／Ｘ（０）ＣＺＴ ≈

１－（１／３！）［π（±１－δ）（Ｎ／Ｗ）（ｆｗ／ｆｓ）］２

１－（１／３！）［πδ（Ｎ／Ｍ）（ｆｗ／ｆｓ）］２

（６）式中·代表ＣＺＴ相应谱线的幅值。

由上式可得δ的估计式：

＾δ＝±１－１－（１－Ａ±１）（１－Ｃ＋ＣＡ±１

槡）

１－Ａ±１

（７）式中：Ｃ＝３！／［π（Ｎ／Ｗ）（ｆｗ／ｆｓ）］２。当Ｘ（－１）ＣＺＴ／Ｘ（１）ＣＺＴ＜１时，上述两式中的“±”取正号，否则取负号。

综合以上分析，可得到卫星信号载频精确估计的过程如下：首先对采样所得数字信号进行Ｎ点ＦＦＴ分析，判断出待测频率的ＦＦＴ主瓣位置，缩小频率估计范围；再在主瓣范围内进行Ｍ点ＣＺＴ细化，以获取更高精度的频率估计值；最后，采用类Ｒｉｆｅ算法进行进一步的修正，得到最终载频估计值。整个算法流程如图１所示。

采用ＦＦＴ与ＣＺＴ联合测频算法，再经类Ｒｉｆｅ算法修正后，在相同的运算量下，载频估计精度将远高于常规ＦＦＴ估计算法。对于常规ＦＦＴ算法，考虑到每一级蝶形运算需进行Ｎ／２次复数乘运算（Ｎ为ＦＦＴ点数）。因此，为了获得ｆｓ／（Ｎ·Ｍ）的分辨率，需进行

（Ｎ·Ｍ／２）ｌｏｇ

２

（Ｎ·Ｍ）次复数乘运算。采用ＦＦＴ和ＣＺＴ综合估计算法后，可先进行Ｎ点的ＦＦＴ运算，再进行Ｍ点的ＣＺＴ运算，则需要进行（Ｎ／２）ｌｏｇ　Ｎ＋

０

３

１现代电子技术２０１２年第３５卷