CN110377872B - 一种基于通用计算显卡的多普勒数据处理方法 - Google Patents

Abstract

本发明提供了一种基于通用计算显卡的多普勒数据处理方法，步骤包括：对连续电波数据流进行分段处理获得电波信号数据段，建立随时间变化的电波信号数据模型；根据电波信号数据模型构建用于微分演化算法的目标函数；设置微分演化算法的初始参数、控制参数以及运行时自动参数；利用通用计算显卡的多线程并行加速对各个电波信号数据段进行微分演化算法处理。该数据处理方法在数据拟合时使用了微分演化算法，因此采用通用计算显卡技术来加速数据处理，拟合得到的相位表达式为泰勒多项式，可以进一步计算得到瞬时相位、频率、频率一阶导数以及不同积分尺度的总相位，这些观测量可以方便地应用于掩星、重力场等行星无线电科学研究中。

Description

一种基于通用计算显卡的多普勒数据处理方法

技术领域

本发明涉及一种多普勒数据处理方法，尤其是一种基于通用计算显卡的多普勒数据处理方法。

背景技术

在深空探测任务中，通常利用无线电手段来探测感兴趣的物理现象，这种研究一般称之为无线电科学由于飞行器和台站的运动和引力环境差异以及信号链路传播空间的介电特质、引力特性，电磁波的物理参数会发生变化，包括频率、幅度、极化等参数的变化，通过电磁波参数变化的监测，可以反演飞行器的运动状态和信号链路传播空间的介电参数、引力参数，这种研究手段被广泛地应用于行星科学、空间物理学以及基础物理学等研究领域。在行星科学研究中，电磁波的相位和幅度变化通常被用来研究行星大气组分、电离层分布以及行星重力场。在大气和电离层的研究中，通过检测电磁波相位和幅度的瞬时变化来反演行星大气和电离层参数以及行星际电离层分布；在行星重力场研究中，为了提高测量精度，一般采用积分多普勒或者总相位作为观测量来解算重力场。这些研究都以高精度的电磁波测量精度作为基础。

在多普勒基带数据处理中，所需计算量甚为庞大，block和RSR接收机采用了专用电路(Application Specific Integrated Circuit，ASIC)板卡来计算相位，通过相位计数算法实现多普勒数据处理。随着通用计算设备计算能力的提高，也有一些通用计算设备被用于多普勒数据的处理，如基于gnuradio平台的USRP(Universal Software RadioPeripheral)板卡，这种板卡通过FFT算法实现相位追踪。基于通用计算显卡的相位计算方法，这种方法通过泰勒多项式拟合基带相位，从而给出数据块中心邻域内解析形式的相位表达式，泰勒多项式拟合使用了微分演化算法来计算多项式系数，对基带数据拟合所需计算甚为庞大。

通过数据拟合得到的相位和频率的表达形式为分段泰勒多项式，通过调整数据块的长度和多项式阶数，相位拟合整体截断误差将保持在噪声水平σ_T(积分时间T对应数据块长度，通常取为2秒)。在行星掩星研究中通常需要高分辨率的多普勒观测，采样频率一般为100Hz，多普勒的提取采用的是FFT算法，多普勒的计算精度为σ_0.01。一般情况下，积分时间越长得到的多普勒观测量精度越高，因此在掩星观测中利用FFT算法并非最优方法。在卡西尼的观测测试中，在频率快速变化情况下(df/dt>1500rad/s)，硬件多普勒设备会发生信号失锁的问题。

发明内容

本发明要解决的技术问题是：在采用微分演化算法后对原始跟踪数据的相位拟合所需计算量很大，需要采用较快的数据加速处理技术。

为了解决上述技术问题，本发明提供了一种基于通用计算显卡的多普勒数据处理方法，包括如下步骤：

步骤1，对连续电波数据流进行分段处理获得电波信号数据段，建立随时间变化的电波信号数据模型；

步骤2，根据电波信号数据模型构建用于微分演化算法的目标函数；

步骤3，设置微分演化算法的初始参数、控制参数以及运行时自动参数；

步骤4，利用通用计算显卡的多线程并行加速对各个电波信号数据段进行微分演化算法处理。

进一步，步骤1中，建立随时间变化的电波信号数据模型的具体步骤为：

首先，将电波信号数据段的通过有限项泰勒多项式展开为：

式(1)中，t为数据时间点，t∈[-T/2,T/2]，T为电波信号数据段的长度，电波信号数据段为中心为时间零点；

然后，以n＝3建立电波信号数据模型，通过相位展开式(1)，得到电波信号数据模型为：

式(2)中，c_0～3为泰勒展开系数：

式(2)中，c₄,c₅为信号线性幅度参数。

进一步，步骤2中，根据电波信号数据模型构建用于微分演化算法的目标函数为：

式(4)中，F_obj(c_k)为目标函数，c_k为待拟合参数向量，s_i为电波信号数据段中采样量化数据点，N为采样量化数据点数，σ_i为数据点的标准偏差，s(c_k,t_i)为式(2)构建的电波信号数据模型，当待拟合参数向量c_k满足收敛条件时，式(4)给出噪声功率。

进一步，步骤3中，微分演化算法的初始参数具体设置为：

c₀为电波信号数据段中心的瞬时相位，c₀的估值范围设置为[0,2π]；

c₁为电波信号数据段中心的瞬时频率，抽取电波信号数据段中心的一段数据使用FFT计算其频率，c₁的估值范围设置为FFT计算频率的±10％；

c₂为电波信号数据段中心的瞬时相位的二阶导数，使用FFT计算电波信号数据段中心邻域内的频率序列，利用五点差分方法得到电波信号数据段中心相位的二阶导数φ²(0)，再通过式(3)计算c₂，范围设置为c₂的±10％；

c₃为电波信号数据段的瞬时相位的三阶导数，设置范围为[-50,50]；

c₄为电波信号数据段的幅度，利用FFT估算信号幅度，范围设置为估算值的±10％；

c₅为电波信号数据段内信号的斜率，范围设置为

进一步，步骤3中，微分演化算法的控制参数具体设置为：

种群数量np＝200；参数个数，即问题的维数dim_xc＝6；种群进化策略strategy＝3；迭代显示频率refresh＝5；目标函数期望值vtr＝0；加权因子cr_xc＝0.85；缩放因子f_xc＝0.5；用于进化策略控制和IO控制参数method＝[0,0,0]。

进一步，步骤3中，微分演化算法的运行时自动参数具体设置为：

c₀与初始设置范围一致，设定为[0,2π]；

c₁由上一个电波信号数据段的频率变化信息估算得到，估算值为：

参数范围为：

*号表示前一个数据块；

c₂的估算值为：

参数范围为：

c₄的参数范围为：

c₃和c₅根据试运行的结果设定参数范围，当积分时间小于10秒时，c₃和c₅近似随机分布，在运行中固定搜索范围；当积分时间大于10秒时，c₃和c₅呈现近似线性趋势。

进一步，步骤4中，在进行微分演化算法处理后，需要对参数拟合质量进行判定，具体步骤为：

首先，根据式(1)的有限项泰勒多项式展开形式，电波信号数据段内的相位和频率分别表示为：

则，相邻电波信号数据段边缘的相位和频率为：

式(6)中，n为电波信号数据段序号，根据MEX三程多普勒跟踪过程中的数据处理结果进行判定，若相邻电波信号数据段的相位和频率在误差范围内是连续的，则微分演化算法处理的满足参数拟合质量要求。

本发明的有益效果在于：和传统基于锁相环技术(PLL)的相位计数处理方法不同，这种方法通过泰勒多项式拟合基带相位，给出解析形式的相位表达式，数据拟合使用了微分演化算法(Differentia Evolution algrithm)，由于对原始跟踪数据的相位拟合所需计算量很大，因此用通用计算显卡技术来加速数据处理，拟合得到的相位表达式为泰勒多项式，可以进一步计算得到瞬时相位、频率、频率一阶导数以及不同积分尺度的总相位(totalcount phase)，这些观测量可以方便地应用于掩星、重力场等行星无线电科学研究中。

附图说明

图1为本发明的方法流程图；

图2为本发明的目标函数和相关系数的函数关系示意图。

具体实施方式

下面结合附图对本发明技术方案进行详细说明，但是本发明的保护范围不局限于所述实施例。

实施例1：

如图1所示，本发明公开的基于通用计算显卡的多普勒数据处理方法包括如下步骤：

步骤1，对连续电波数据流进行分段处理获得电波信号数据段，建立随时间变化的电波信号数据模型；

步骤2，根据电波信号数据模型构建用于微分演化算法的目标函数；

步骤3，设置微分演化算法的初始参数、控制参数以及运行时自动参数；

步骤4，利用通用计算显卡的多线程并行加速对各个电波信号数据段进行微分演化算法处理。

本发明通用计算显卡采用的是双K80显卡，双K80显卡的性能参数为：带宽768×2bit，核心数4992×2个，标称双精度计算能力3728/Gflops，计算速度约为GTX580的20倍，能够完成数据的实时处理；微分演化算法是成熟的遗传(GA,Genetic Algorithm)算法，和通常的遗传算法一样，最初的种群是在参数空间内随机选取的，在迭代过程中，通过父代微分交叉优化，得到子代参数矢量，微分演化算法的实现可以分为以下步骤：

(1)给出参数搜索区域范围、迭代次数限制和收敛精度阈值；

(2)在参数空间产生随机种群，计算个体适应度；

(3)通过加权微分交叉操作产生子代个体；

(4)收敛性判断。

由于在微分演化算法的迭代处理过程存在因果性，迭代过程自身无法实现并行化，因此在进行相位拟合计算中本发明采用了细粒度并行(Fine-grained parallelism)的方案实现计算加速，即仅仅在计算目标函数(objective function)时采用GPU多线程并行加速。

GPU全称为图形处理单元(Graphic Processing Unit)，与CPU的设计理念不同(处理逻辑操作和数值计算)，早期的GPU是专门为处理图形任务而设计的芯片，GPU可以完成图形图像复杂的渲染、纹理以及光线追踪等操作，被广泛地应用于游戏和动画制作领域。随着GPU的性能增强，相应的软件开发包(Software Development Kit，SDK)CUDA，通过CUDA软件包，使用者可以将矩阵和矢量的操作转换为GPU擅长的纹理及渲染操作，从而实现数值计算的目的。第一个专门应用于数值计算的GPU芯片Tesla C870以及相应的应用包CUDA1.0，早期的CUDA应用包较为简单，只有单精度和整数的矢量和矩阵操作，不过在数值计算性能上已经超越了同时期的CPU。随着显卡性能的不断提升，NVIDIA公司发布了专门应用于数值计算的通用计算显卡芯片，如Tesla、Fermi、Kepler、Maxwell等系列，由于GPGPU的计算性能大大超过CPU的性能，这些专用显卡被广泛应用于超级计算平台，实验证实GPU的计算性能远远超过CPU的计算性能。

进一步，步骤1中，建立随时间变化的电波信号数据模型的具体步骤为：

首先，将电波信号数据段的通过有限项泰勒多项式展开为：

式(1)中，t为数据时间点，t∈[-T/2,T/2]，T为电波信号数据段的长度，电波信号数据段为中心为时间零点；

然后，以n＝3建立电波信号数据模型，通过相位展开式(1)，得到电波信号数据模型为：

多项式展开阶数n与被观测对象的运动状态有关，n＝3在大多数情况下可以满足数据处理精，式(2)中，c_0～3为泰勒展开系数：

式(2)中，c₄,c₅为信号线性幅度参数，当电波信号数据段较短(1～5秒)时，信号幅度变化很小，c₅的不确定度通常较大，但并不影响相位的估算精度，在实际的数据试验中，发现引入c₅可以得到较好的收敛效果。

进一步，步骤2中，根据电波信号数据模型构建用于微分演化算法的目标函数为：

式(4)中，F_obj(c_k)为目标函数，c_k为待拟合参数向量，s_i为电波信号数据段中采样量化数据点，N为采样量化数据点数，σ_i为数据点的标准偏差，s(c_k,t_i)为式(2)构建的电波信号数据模型，当待拟合参数向量c_k满足收敛条件时，式(4)给出噪声功率，在相关系数的零点位置处，目标函数迅速收敛，如图2所示。

进一步，在信号模型参数估算环节需要设置模型参数的搜索范围，将最优解矢量包含在搜索范围内，由于微分演化算法为全域最优化算法，理论上可以将搜索范围尽量放大以包含最优解矢量，但过大的搜索范围相应地需要扩大种群数目，影响计算效率，因此合理地设置搜索范围是有必要的。搜索范围的设置包括初始参数设置和运行时自动参数设置，以式(4)的目标函数模型为例分开描述。初始参数设置需要预先估计参数值，由式(3)可知c_k参数和电波信号数据段中心相位的导数有关，因此通过估算相位中心的导数可以得到各个参数的预估值，于是步骤3中，微分演化算法的初始参数具体设置为：

c₀为电波信号数据段中心的瞬时相位，c₀的估值范围设置为[0,2π]；

c₁为电波信号数据段中心的瞬时频率，抽取电波信号数据段中心的一段数据使用FFT计算其频率，c₁的估值范围设置为FFT计算频率的±10％；

c₂为电波信号数据段中心的瞬时相位的二阶导数，使用FFT计算电波信号数据段中心邻域内的频率序列，利用五点差分方法得到电波信号数据段中心相位的二阶导数φ²(0)，再通过式(3)计算c₂，范围设置为c₂的±10％；

c₃为电波信号数据段的瞬时相位的三阶导数，设置范围为[-50,50]；

c₄为电波信号数据段的幅度，利用FFT估算信号幅度，范围设置为估算值的±10％；

c₅为电波信号数据段内信号的斜率，范围设置为

进一步，步骤3中，微分演化算法的控制参数具体设置为：

种群数量np＝200；参数个数，即问题的维数dim_xc＝6；种群进化策略strategy＝3；迭代显示频率refresh＝5；目标函数期望值vtr＝0；加权因子cr_xc＝0.85；缩放因子f_xc＝0.5；用于进化策略控制和IO控制参数method＝[0,0,0]。

进一步，步骤3中，微分演化算法的运行时自动参数设置就是系统对相位的自动追踪，追踪过程用到了参数的外推法(extrapolation)，具体设置为：

c₀与初始设置范围一致，设定为[0,2π]；

c₁由上一个电波信号数据段的频率变化信息估算得到，估算值为：

参数范围为：

*号表示前一个数据块；

c₂的估算值为：

参数范围为：

c₄的参数范围为：

在跟踪弧段处于轨道近点处，由于频率的变化率变化很快，上面c₁的估值范围可能无法覆盖真值，导致搜索失败，系统会对数据块控制参数自动重置，扩大参数覆盖范围重新处理该电波信号数据段。c₃和c₅根据试运行的结果设定参数范围，当积分时间小于10秒时，c₃和c₅近似随机分布，在运行中固定搜索范围；当积分时间大于10秒时，c₃和c₅呈现近似线性趋势。

进一步，步骤4中，通用计算显卡具有多线程并行加速数据处理功能，因此采用通用计算显卡的多线程并行加速对各个电波信号数据段进行微分演化算法处理是直接可行的，在进行微分演化算法处理后需要对参数拟合质量进行判定，参数拟合的质量有一个简单的判定方法，那就是相邻数据块边缘的相位(2π取模)和频率是否连续，于是在进行微分演化算法处理后，需要对参数拟合质量进行判定，具体步骤为：

首先，根据式(1)的有限项泰勒多项式展开形式，电波信号数据段内的相位和频率分别表示为：

则，相邻电波信号数据段边缘的相位和频率为：

式(6)中，n为电波信号数据段序号，根据MEX三程多普勒跟踪过程中的数据处理结果进行判定，MEX的多普勒数据处理的精度约为50mrad(1σ，1秒积分长度)，若相邻电波信号数据段的相位和频率在误差范围内是连续的，则微分演化算法处理的满足参数拟合质量要求。

在比对试验中，通常锁相环方法在追踪相位时一般都是基于FFT算法，通过频谱的极值点追踪相位，当频率变化太大、频谱展宽时会发生相位失锁的情况。对于这种情况原理上可以通过动态调整接收机本振(NCO)频率来解决，即让接收机的本振频率变化率和信号的频率变化率相当，但是目前接收机还不具有动态调节本振的功能，另外通过动态小波分析算法也可以实现快速变化相位的追踪。在卡西尼的观测实验中，有两次碰到了这种情况，一次是2018/6/10的轨道近点跟踪弧段，另外一次是卡西尼的坠落弧段。使用了PRSR(Planetary Radio science Receiver)接收机自带的硬件多普勒处理机和通用显卡处理设备对数据进行了处理，发现硬件多普勒处理机在频率变化率大于1500rad/s时相位失锁，无法正常处理数据，而通用显卡处理设备则可以正常处理数据，得到X波段卡西尼的多普勒处理精度为20mrad/s|_1s、4mrad/s|_30s。

本发明给出的基于通用计算显卡(GPGPU)技术的深空飞行器多普勒数据处理方法。本发明的处理方法中用到的GPU加速技术和微分演化算法，给出了基于泰勒多项式拟合的相位处理方法，通过泰勒多项式系数可以进一步合成瞬时多普勒、积分多普勒、总相位和视向加速度观测量，由于采用了拟合的方法，瞬时多普勒和视向加速度具有积分观测量的精度水平。

如上所述，尽管参照特定的优选实施例已经表示和表述了本发明，但其不得解释为对本发明自身的限制。在不脱离所附权利要求定义的本发明的精神和范围前提下，可对其在形式上和细节上作出各种变化。

Claims (6)

1.一种基于通用计算显卡的多普勒数据处理方法，其特征在于，包括如下步骤：

步骤1，对连续电波数据流进行分段处理获得电波信号数据段，建立随时间变化的电波信号数据模型；

步骤2，根据电波信号数据模型构建用于微分演化算法的目标函数；

步骤3，设置微分演化算法的初始参数、控制参数以及运行时自动参数；

步骤4，利用通用计算显卡的多线程并行加速对各个电波信号数据段进行微分演化算法处理；

步骤1中，建立随时间变化的电波信号数据模型的具体步骤为：

首先，将电波信号数据段的通过有限项泰勒多项式展开为：

式(1)中，t为数据时间点，t∈[-T/2,T/2]，T为电波信号数据段的长度，电波信号数据段为中心为时间零点；

然后，以n＝3建立电波信号数据模型，通过相位展开式(1)，得到电波信号数据模型为：

式(2)中，c_0～3为泰勒展开系数：

式(2)中，c₄,c₅为信号线性幅度参数。

2.根据权利要求1所述的基于通用计算显卡的多普勒数据处理方法，其特征在于，步骤2中，根据电波信号数据模型构建用于微分演化算法的目标函数为：

式(4)中，F_obj(c_k)为目标函数，c_k为待拟合参数向量，s_i为电波信号数据段中采样量化数据点，N为采样量化数据点数，σ_i为数据点的标准偏差，s(c_k,t_i)为式(2)构建的电波信号数据模型，当待拟合参数向量c_k满足收敛条件时，式(4)给出噪声功率。

3.根据权利要求2所述的基于通用计算显卡的多普勒数据处理方法，其特征在于，步骤3中，微分演化算法的初始参数具体设置为：

c₀为电波信号数据段中心的瞬时相位，c₀的估值范围设置为[0,2π]；

c₁为电波信号数据段中心的瞬时频率，抽取电波信号数据段中心的一段数据使用FFT计算其频率，c₁的估值范围设置为FFT计算频率的±10％；

c₂为电波信号数据段中心的瞬时相位的二阶导数，使用FFT计算电波信号数据段中心邻域内的频率序列，利用五点差分方法得到电波信号数据段中心相位的二阶导数φ²(0)，再通过式(3)计算c₂，范围设置为c₂的±10％；

c₃为电波信号数据段的瞬时相位的三阶导数，设置范围为[-50,50]；

c₄为电波信号数据段的幅度，利用FFT估算信号幅度，范围设置为估算值的±10％；

c₅为电波信号数据段内信号的斜率，范围设置为

4.根据权利要求2所述的基于通用计算显卡的多普勒数据处理方法，其特征在于，步骤3中，微分演化算法的控制参数具体设置为：

种群数量np＝200；参数个数，即问题的维数dim_xc＝6；种群进化策略strategy＝3；迭代显示频率refresh＝5；目标函数期望值vtr＝0；加权因子cr_xc＝0.85；缩放因子f_xc＝0.5；用于进化策略控制和IO控制参数method＝[0,0,0]。

5.根据权利要求2所述的基于通用计算显卡的多普勒数据处理方法，其特征在于，步骤3中，微分演化算法的运行时自动参数具体设置为：

c₀与初始设置范围一致，设定为[0,2π]；

c₁由上一个电波信号数据段的频率变化信息估算得到，估算值为：

参数范围为：

*号表示前一个数据块；

c₂的估算值为：

参数范围为：

c₄的参数范围为：

c₃和c₅根据试运行的结果设定参数范围，当积分时间小于10秒时，c₃和c₅近似随机分布，在运行中固定搜索范围；当积分时间大于10秒时，c₃和c₅呈现近似线性趋势。

6.根据权利要求1所述的基于通用计算显卡的多普勒数据处理方法，其特征在于，步骤4中，在进行微分演化算法处理后，需要对参数拟合质量进行判定，具体步骤为：

首先，根据式(1)的有限项泰勒多项式展开形式，电波信号数据段内的相位和频率分别表示为：

则，相邻电波信号数据段边缘的相位和频率为：

式(6)中，n为电波信号数据段序号，根据MEX三程多普勒跟踪过程中的数据处理结果进行判定，若相邻电波信号数据段的相位和频率在误差范围内是连续的，则微分演化算法处理的满足参数拟合质量要求。

Images