CN114039823A

CN114039823A - 一种低信噪比下高动态多普勒频偏及其加速度估计方法

Info

Publication number: CN114039823A
Application number: CN202111304237.5A
Authority: CN
Inventors: 刘慧�; 张兆维
Original assignee: Jinling Institute of Technology
Current assignee: Nanjing Xiongjun Software Technology Co ltd
Priority date: 2021-11-05
Filing date: 2021-11-05
Publication date: 2022-02-11
Anticipated expiration: 2041-11-05
Also published as: CN114039823B

Abstract

本发明公开了一种低信噪比下高动态多普勒频偏及其加速度估计方法，包括：依次接收连续的TL个码元信号；对每一个码元信号进行快速傅里叶变换；每L个连续码元信号进行非相干累积，得到T个累积结果；从每一个累积结果中选出Nc个元素；计算聚类终止参数；对TNc个元素进行凝聚层次聚类；计算每一个簇的大小，并寻找最大的簇；根据最大簇中元素序号的平均值，得到多普勒频偏的估计结果；将最大簇的元素划分到T个子簇；计算非空子簇中元素序号的平均值；根据所有非空子簇，得到多普勒频偏加速度的估计结果。

Description

一种低信噪比下高动态多普勒频偏及其加速度估计方法

技术领域

本发明属于空间通信领域，尤其涉及一种低信噪比下高动态多普勒频偏及其加速度估计方法。

背景技术

随着空间技术的快速发展，空间通信载体(火箭、卫星和航天器等)早已由近地空间扩张到深空领域。它们的飞行距离越来越远，且飞行速度也越来越快。在信号传输过程中，空间通信通常面临以下两个问题：(1)通信距离(几十，甚至几百亿公里)的增大会造成很大的路径损耗，从而严重降低接收信号的信噪比；(2)通信载体的高速(超过第三宇宙速度16.7km/s)和大加速度(可达几百g，g＝9.8m/s²)导致信号载波频率发生很高动态的多普勒频偏及加速度。为了恢复信号，接收端的载波同步必须补偿多普勒频偏。补偿多普勒频偏的关键就在于如何准确地估计多普勒频偏及其加速度。因此，低信噪比下高动态多普勒频偏及其加速度估计就成为载波频率同步的核心技术之一。

当前，空间通信的载波多普勒频偏及其加速度估计方法主要分为以下两类：滑动相关和快速傅里叶变换(fast Fourier transform，FFT)。滑动相关方法首先将频率搜索范围和加速度搜索范围划分为若干个小间隔进行二维搜索，然后通过不断地调整本地频率并与接收信号进行相关处理，直至相关值超过检测门限。但是，在多普勒频偏高动态变化的情况下，滑动相关方法需要搜索的次数过多，从而导致估计时间很长。FFT方法首先将接收到的码元信号从时域转换到频域，然后在频域上进行最大似然搜索得到频域谱峰，是当前应用最为广泛的方法。但是，在低信噪比条件下，FFT对单个码元信号的处理不能提供足够的输出信噪比，使得信号元素谱线低于最大噪声谱线，无法完成估计。为此，FFT方法通常需要长时间信号累积以提高信噪比。另外，当通信载体在高速运动，尤其是存在很大加速度的情形下，不同码元信号上的多普勒频偏会发生很大的漂移，使得他们的频域谱线不在同一个点上。此时，多个信号的累积能量会扩散到多个FFT点，导致严重的能量扩散问题。现有研究表明，多普勒频偏的检测概率是信号累积周期长度的凸函数并且存在唯一的最优累积周期长度，而且加速度越大，最优累积周期长度越小，检测概率也会越低。这说明，在大加速度下，更长时间的累积周期不仅不能提高检测概率，反而会降低检测概率。除此之外，FFT方法缺乏时间维度的信息，并不能估计多普勒频偏的加速度。

通过以上的分析，可以发现：一方面，低信噪比要求长时间信号累积以提高信号能量；另一方面，大加速度引起的能量扩散问题严重限制了信号累积效果。因此，低信噪比下的高动态多普勒频偏及其加速度估计就成为一个极富挑战性的技术难题。

发明内容

发明目的：本发明所要解决的技术问题是针对现有技术的不足，提供一种低信噪比下高动态多普勒频偏及其加速度估计方法，包括如下步骤：

步骤1，依次接收TL个连续码元信号；

步骤2，对每一个码元信号进行快速傅里叶变换；

步骤3，每L个连续码元信号进行非相干累积；

步骤4，从每个累积结果中选取Nc个元素；

步骤5，计算聚类终止参数；

步骤6，对TNc个元素进行凝聚层次聚类；

步骤7，计算每一个簇的大小，并寻找最大簇；

步骤8，根据最大簇中元素的序号平均值，得到多普勒频偏的估计结果；

步骤9，将最大簇的元素划分到T个子簇；

步骤10，计算非空子簇中元素序号的平均值；

步骤11，根据所有非空子簇，得到多普勒频偏加速度的估计结果。

步骤1包括：信号接收端依次接收TL个连续码元信号{x^[0],...,x^[TL-1]}，其中第i个信号

表示为：

式中，i＝0,...,TL-1，

表示第i个码元信号x^[i]的第n个元素，a是码元振幅，b是比特速率为r_b的二进制数据信息，j表示复数，f_s是采样频率，f_d,0是码元上的多普勒初始频偏，f_a是多普勒频偏加速度，

是初始相位，w_n是均值为0、方差为σ²的高斯白噪声；N＝f_s/r_b是一个码元信号长度。

步骤2包括：对每一个长度为N的码元信号进行N_f点的快速傅里叶变换，得到TL个频域向量，其中，N_f≥N，第i个码元信号x^[i]对应的N_f点频域向量I^[i]为：

I^[i]＝fft(x^[i]) (2)

式中，fft(·)表示N_f点的复数快速傅里叶变换。

步骤3包括：每L个连续码元信号进行非相干累积，得到T个累积结果

其中，第t个累积结果

如下所示：

式中，

表示累积结果

中的第i个元素，|·|表示取绝对值。

步骤4包括：从每一个累积结果中选取N_c个元素，其中，第t个累积结果

对应的元素集合为

如下所示：

式中，i＝0,...,N_c-1，1≤N_c＜＜N_f，

表示元素集合

的第i个元素，

表示

中等于

的元素，s(i)表示

在

中的序号。

步骤5包括：计算聚类终止参数δ，如下所示：

式中，f_a,max表示多普勒频偏加速度f_a的最大值。

步骤6包括：对

中的TN_c个元素进行凝聚层次聚类，具体包括：

步骤6-1，计算任意两个元素

和

之间的距离

如下所示：

式中，t＝0,1,...,T-1，t'＝0,1,...,T-1，i＝0,1,...,N_c-1，i'＝0,1,...,N_c-1；

步骤6-2，每一个元素自成一簇，构成TN_c个簇

如下所示：

式中，

表示第tN_c+i个簇，t＝0,...,T-1，i＝0,...,N_c-1，

表示簇的总个数；

步骤6-3，计算簇间距离矩阵M，如下所示：

式中，NULL表示空值(不作任何运算)，M_i,k表示距离矩阵M中的第i行第k列元素，x1表示簇C_i中的元素，x2表示簇C_k中的元素，N_c,i和N_c,k分别表示簇C_i和C_k的元素个数；

步骤6-4，计算簇间距离矩阵M中的最小元素

如下所示：

式中，i^*<k^*；

步骤6-5，判断

是否满足：

如果满足，进入步骤6-6，否则进入步骤6-10；

步骤6-6，合并两个簇

和

如下所示：

步骤6-7，删除簇

并更新簇的序号，如下所示：

步骤6-8，更新簇的总个数

如下所示：

步骤6-9，更新簇间距离矩阵M，返回步骤6-3；

步骤6-10，)得到聚类结果

其中，C_i表示聚类结果中的第i个簇，

步骤7包括：计算每一个簇C_i的大小S_i，如下所示：

然后，根据

寻找最大的簇

其中，i_max表示最大簇的序号。

步骤8包括：计算最大簇

中所有元素序号的平均值

如下所示：

然后，根据

得到多普勒频偏f_d的估计结果

如下所示：

步骤9包括：将最大簇

中所有元素划分到T个子簇

其中第t个子簇

是由既属于

又属于

的所有元素组成，如下所示：

步骤10包括：如果子簇

非空，计算其元素序号的平均值

如下所示：

式中，

表示非空子簇

中元素的个数；

步骤11包括：根据所有非空子簇，得到多普勒频偏加速度f_a的估计结果

如下所示：

式中，

表示非空子簇中元素序号平均值

的平均值，

表示非空子簇序号t的平均值，n_c,max表示非空子簇的总个数。

有益效果：本发明将多个信号元素聚为一个簇，能够聚集已经扩散的信号能量，能够提高多普勒频偏的估计性能。同时，本发明将最大簇划分为多个子簇，能够提供时间维信息，从而得到多普勒频偏加速度的估计结果。

附图说明

下面结合附图和具体实施方式对本发明做更进一步的具体说明，本发明的上述和/或其他方面的优点将会变得更加清楚。

图1是本发明流程图。

图2是本发明方法下多普勒频偏的检测概率。

图3是本发明方法下多普勒频偏加速度的估计结果。

具体实施方式

本发明提供一种低信噪比下高动态多普勒频偏及其加速度估计方法，首先，利用聚类思想将多个累积周期中的信号能量聚为一个簇，寻找最大簇作为多普勒频偏的估计结果，然后，将最大簇划分为多个子簇，得到多普勒频偏加速度的估计结果。

本发明具体包括如下步骤：

步骤1：信号接收端依次接收TL个连续码元信号{x^[0],...,x^[TL-1]}，其中第i个信号

表示为：

式中，i＝0,...,TL-1，a是码元振幅，b是比特速率为r_b的二进制数据信息，f_s是采样频率，f_d是码元上的多普勒初始频偏，f_a是多普勒频偏加速度，

是初始相位，w_n是均值为0、方差为σ²的高斯白噪声，N＝f_s/r_b是一个码元信号长度。

步骤2：对每一个长度为N的码元信号进行N_f(N_f≥N)点的快速傅里叶变换，得到TL个频域向量，其中，第i个码元信号x^[i]对应的N_f点频域向量I^[i]为：

I^[i]＝fft(x^[i]) (2)

式中，fft(·)表示N_f点的复数快速傅里叶变换。

步骤3：每L个连续码元信号进行非相干累积，得到T个累积结果

其中，第t个累积结果

如下所示：

式中，|·|表示取绝对值。

步骤4：从每一个累积结果中选取N_c个元素，其中，第t个累积结果

对应的元素集合为

如下所示：

式中，i＝0,...,N_c-1，1≤N_c＜＜N_f，s(i)表示

在

中的序号。

步骤5：计算聚类终止参数δ，如下所示：

式中，f_a,max表示多普勒频偏加速度f_a的最大值。

步骤6：对

中的TN_c个元素进行凝聚层次聚类。

6-1)计算任意两个元素

和

之间的距离

如下所示：

6-2)每一个元素自成一簇，构成TN_c个簇

如下所示：

式中，t＝0,...,T-1，i＝0,...,N_c-1，

表示簇的总个数。

6-3)计算簇间距离矩阵M，如下所示：

式中，N_c,i和N_c,k分别表示簇C_i和C_k的元素个数；

6-4)计算簇间距离矩阵M中的最小元素，如下所示：

式中，i^*<k^*。

6-5)判断

是否满足

如果满足，进入下一步6-6)，否则进入6-10)。

6-6)合并两个簇

和

如下所示：

6-7)删除簇

并更新簇的序号，如下所示：

6-8)更新簇的总个数

如下所示：

6-9)更新簇间距离矩阵M，返回6-3)。

6-10)得到聚类结果

步骤7：计算每一个簇C_i的大小S_i，如下所示：

然后，根据

寻找最大的簇

其中，i_max表示最大簇的序号。

步骤8：计算最大簇

中所有元素序号的平均值

如下所示：

然后，根据

得到多普勒频偏f_d的估计结果

如下所示：

步骤9：将最大簇

中所有元素划分到T个子簇

其中第t个子簇

是由既属于

又属于

的所有元素组成，如下所示：

步骤10：如果子簇

非空，计算其元素序号的平均值

如下所示：

式中，

表示非空子簇

中元素的个数。

步骤11：根据所有非空子簇，得到多普勒频偏加速度f_a的估计结果

如下所示：

式中，

和

n_c,max表示非空子簇的总个数。

实施例

本发明的实施例如下：

以空间飞行器测控通信为例，信号载波频率f_c＝30GHz，飞行器运动速度和加速度的最大值分别为18×10³m/s和-120g(g＝9.8m/s²)。根据运动速度和加速度可以得到多普勒f_d,0及其加速度f_a分别为1.8MHz和-118kHz/s，多普勒频偏加速度最大值f_a,max为125kHz/s，码元二进制数据速率为5kb/s，采样频率f_s＝5.12MHz，码元信号初始相位

多普勒频偏估计可容忍误差门限θ＝10×10³Hz，FFT点数N_f＝1024，累积周期个数T＝3，累积周期长度L＝250，元素个数N_c＝20，信噪比范围SNR＝a²/σ²∈[-9,-2]dB。

步骤1：信号接收端依次接收TL＝750个连续码元信号{x^[0],...,x^[749]}。

步骤2：对每一个长度为N＝1024的码元信号进行N_f＝1024点的快速傅里叶变换，得到750个频域向量{I^[0],...,I^[749]}。

步骤3：每250个连续码元信号进行非相干累积，得到3个累积结果

其中，第t个累积结果

如下所示：

三个累积周期内的信号元素位置分别为(360,359)、(359,358)和(358,357)，总共六个信号元素。

步骤4：从每一个累积结果中选取20个元素，其中，第t个累积结果

对应的元素集合为

如下所示：

三个累积周期内只有(360)、(358)和(358,357)共四个信号元素被选取出来，属于较大元素。

步骤5：计算聚类终止参数δ，如下所示：

步骤6：对

中的TN_c＝60个元素进行凝聚层次聚类。

6-1)计算任意两个元素

和

之间的距离

如下所示：

6-2)每一个元素自成一簇，构成TN_c个簇

如下所示：

6-3)计算簇间距离矩阵M，如下所示：

6-4)计算簇间距离矩阵M中的最小元素，如下所示：

式中，i^*<k^*。

6-5)判断

是否满足

如果满足，进入下一步6-6)，否则进入6-10)。

6-6)合并两个簇

和

如下所示：

6-7)删除簇

并更新簇的序号，如下所示：

6-8)更新簇的总个数

如下所示：

6-9)更新簇间距离矩阵M，返回6-3)。

6-10)得到聚类结果

聚类结果中，三个累积周期的四个信号元素(360)、(358)和(358,357)会聚为一个簇，而其余的噪声元素会聚为很多簇，因为噪声元素的位置是随机的。

步骤7：计算每一个簇C_i的大小S_i，如下所示：

然后，根据

寻找最大的簇

聚类结果的所有簇中，四个信号元素(360)、(358)和(358,357)构成的簇成为最大的簇，相当于起到了能量累积的作用。

步骤8：计算最大簇

中所有元素序号的平均值

如下所示：

然后，根据

得到多普勒频偏f_d的估计结果

如下所示：

根据四个信号元素(360)、(358)和(358,357)的平均值358.25，得出多普勒频偏为

步骤9：将最大簇

中所有元素划分到T＝3个子簇

其中第t个子簇

是由既属于

又属于

的所有元素组成，如下所示：

第一个子簇为

第二个子簇为

第三个子簇为

步骤10：如果子簇

非空，计算其元素序号的平均值

如下所示：

三个子簇的元素序号的平均值分别为360、358、357.5。

如下所示：

根据三个子簇

和

得到多普勒频偏加速度的估计结果

从附图2可以看出，短时傅里叶变换方法通常要求较高的信噪比条件，其检测概率随信噪比的下降而急剧下降。快速傅里叶变换方法的检测概率较低且不能估计多普勒频偏加速度。搜索范围校正方法可以通过缩小搜索范围来提高检测概率。本发明方法将多个累积周期的信号能量重新聚集起来，提高了信号累积增益，从而显著地提高多普勒频偏检测概率。从附图3可以看出，本发明方法将最大簇中元素划分到多个子簇，从而估计多普勒频偏加速度。

在空间飞行器测控通信中，多普勒频偏f_d及其加速度f_a的估计结果

和

为信号发送端和接收端之间载波频率同步和解调提供重要依据，从而为飞行器测信息的可靠传输奠定基础。

本发明提供了一种低信噪比下高动态多普勒频偏及其加速度估计方法，具体实现该技术方案的方法和途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。