0 引言

飞行器和弹道导弹的精确制导需要精确计算地球外部扰动引力，利用地球重力场位系数模型可以有效恢复地球外部重力场，但随着阶数的升高，计算效率显著下降，计算耗时呈指数增加^[1-3]。现有的计算方法，如交换求和顺序、快速傅里叶变换(fast Fourier transform，FFT)均是对规则格网分布的扰动引力计算进行加速，很难对不规则飞行轨迹点进行单点快速计算^[4-5]；利用换极法可以对飞行轨迹所在弧段进行换极，使得换极后飞行轨迹位于“新地球”的同一纬度圈上以减少缔合勒让德函数的计算次数，但对不规则飞行轨迹的加速效果较差^[6]；李红伟等^[7]利用MPI(message passing interface)算法对计算扰动引力进行并行加速，加速比达到1.86，但设备利用率低，加速比不高。近些年随着并行技术的发展，异构并行作为更高效的并行技术，已经在各行各业中崭露头角^[8-11]。尤其是NVIDIA (英伟达)公司推出的CUDA(compute unified device architecture)异构运算平台，利用CPU+GPU(显卡)的异构编程模式，并充分使用CPU的调控能力和GPU的计算能力，在社会发展的各个领域获得了广泛应用，并取得了丰硕的成果^[12-15]。因此，本文提出基于向量化原理对原有的扰动引力计算公式进行改进，在保证精度的条件下利用CUDA异构并行算法实现单点扰动引力计算的并行化，从而达到扰动引力快速计算的目的。最后通过仿真实验对其性能进行评估。

1 扰动引力计算公式 1.1 缔合勒让德函数计算

缔合勒让德函数的计算常使用递推法^[16]。主要递推方法有4种：标准前向列推法、标准前向行推法、Belikov递推法和跨阶次递推法。标准前向列推法和标准前向行推法在高阶时难以保证精度，跨阶次递推法与Belikov递推法则可以计算超高阶缔合勒让德函数^[17-19]。因此，本文分别选用跨阶次递推法与Belikov递推法计算缔合勒让德函数。

1.1.1 跨阶次递推法

根据跨阶次递推原理，阶次为n、级次为m的完全规格化缔合勒让德函数P_{n, m}(cos θ)的表达式为：

其中：

式中：t=cos θ；θ表示余纬。在实际应用过程中，由于a₀，b₀，α_nm，β_nm，γ_nm仅与阶数有关而与点位无关，因此可以提前计算出结果，从而作为已知常数备用以提高计算效率。

1.1.2 Belikov递推法

Belikov递推法公式如下：

式中：u=sin θ；表示非完全正常化的勒让德函数。将转化为P_{n, m}(cos θ)：

其中：

式中，仅与阶数有关而与点位无关，可提前计算作为已知常数备用以提高计算效率。

1.2 扰动引力计算

地球重力场模型是用一个截断到有限阶次的引力位球谐函数级数来逼近地球重力场^[4]。以半径为地球平均半径R的球来近似表示地球，可得到扰动引力位级数式如下：

式中：G为万有引力常量；M为地球质量；ρ为轨迹点到地心的向径距离；N为截断阶数；C_{n, m}^*、S_{n, m}为地球重力场模型参数；λ为地形经度。将式(4)分别对3个坐标方向求偏导即可得到扰动引力三分量。在当地东北天坐标系下，T(ρ, θ, λ)三分量公式如下：

式中，φ为地心纬度，φ=90°-θ。

2 扰动引力计算向量化并行算法

利用地球重力场位系数模型计算扰动引力可以划分为两个部分：缔合勒让德函数的递推和球谐级数的运算。缔合勒让德函数递推过程各个阶次计算的前后相关性不利于并行计算。而球谐级数的运算过程内部并无相关性，且包含两重累加运算，随着截断阶数的增加，计算耗时呈指数上涨，在整个扰动引力计算过程中占据主要部分，因此可作基于向量化的并行计算。

2.1 向量化原理及方法

向量化就是将标量运算转换成向量运算的过程，目的是将原有计算改进为向量间的运算，以达到减少循环计算的目的。向量化的过程可以分为两步进行：一是数据预处理过程向量化；二是指令向量化(即计算公式向量化)^[10]。

2.1.1 数据预处理向量化

数据预处理向量化的目的是将计算过程所需数据，通过一定顺序进行向量化排列，以满足向量运算法则。原有数据存储方式为下三角矩阵形式，在并行化过程中若不改变存储方式，则要进行矩阵运算(即高维向量运算)，需要将矩阵空余部分以零元素补齐，这样在并行计算过程中就会有一半的线程在计算零元素，从而影响计算效率。因此，为更好达到向量化运算效果，需将以高维向量形式排列的数据进行降维处理，使之有序地排列至一维向量中。为有效有序地进行数据预处理向量化过程，本文将高维向量数据按照如图 1所示方法进行降维向量化存储。

如图 1所示，将左侧以三角矩阵形式存储的数据，以行为单位，同行数据由左至右依次排列至右侧一维向量中，得到一个长度为J的向量。依据上述向量化排列准则，分别定义C_1×J，S_1×J, V_1×J, K_1×J4个1×J阶行向量：

定义P_J×1，dP_J×1，R_1J×1，R_2J×1，M_J×15个J×1的列向量：

2.1.2 指令向量化

指令的向量化过程是在原串行逻辑函数算法的基础上，利用向量化方法对函数进行向量运算改造，以适用于并行算法，提高计算效率。现有扰动引力计算公式是以串行逻辑运算为基础展开的，计算均为标量运算，每个运算步骤间有先后顺序。扰动引力的计算是通过循环截断阶次n和级次m实现的，并以n、m的递增为先后顺序，每计算一个分量，需要循环(n+1)×(n+2)/2次。假定截断阶次n=2 160，则每计算一个分量需要进行233万次循环。并且，在进行并行计算时，不同线程无法同时对n、m进行赋值，依据原公式难以实现并行化。因此，为实现向量运算并行化，必须进行指令向量化，也就是对原有串行逻辑公式进行向量化改进，实现计算逻辑并行化。

本文通过数据预处理向量化所得向量，通过向量间加减、内积和对应元素相乘等运算，对原有标量计算公式进行向量化改进。将式(6)、(7)中向量代入式(5)，将累加求和过程通过向量运算替代，可得到向量化后的扰动引力计算公式：

式中：“·”表示向量内积；“ ° ”表示同维向量对应元素相乘。

此时，扰动引力计算过程已从传统标量计算转换为向量运算，原本复杂的循环过程被转换为简单的向量运算，实现了指令向量化。扰动引力向量化计算公式相比于标量公式，每一步向量间的运算都包含了多步标量运算，但这些标量运算间无先后顺序，可以同时进行，成为并行计算实现的理论基础。

2.2 并行算法实现

向量化后，扰动引力的并行计算过程可以分为两步: 一是各个向量对应元素的运算；二是向量内部元素的归约求和。在进行对应元素运算时，首先要确定计算单元(即每个CUDA核心所计算的部分)。根据向量化公式，可得计算单元：

式中：i表示第i个CUDA计算核心；M_i，C_i，S_i，P_i，dP_i，R_1i，R_2i分别表示对应向量中第i个元素。这样就将向量中不同位置对应元素的运算分配给了不同的GPU线程。

在完成各个向量对应元素间的计算后，得到3个新向量t_ρ，t_φ，t_λ，然后分别对这3个向量进行向量内部元素的归约求和。传统串行逻辑运算需要进行J-1轮循环才可得到向量归约和。为更高效地计算求和过程，引入如图 2所示阶梯求和方案。

如图 2所示，将待求和向量各个元素相邻两两结合为一组(若元素总数为奇数，则最后一个元素单独为一组)，利用不同GPU线程同时分别计算每一组元素之和；经过一轮计算，得到一个新的向量；再以此向量定为待求向量，重复上述过程；经过int[(J+1)/2]轮计算，即可得到向量归约和。

3 实验计算分析

为验证本文所提方法的可靠性和有效性，利用模拟飞行轨迹点进行实验。分别使用两种硬件环境：硬件环境1，便携式笔记本电脑，处理器为Intel^®Core^TM i7-4700MQ CPU @ 2.40 GHz；GPU为NVIDIA GTX 765 M；内存8.00 GB。硬件环境2，台式高性能工作站，处理器为Intel^®Core^TM i9-9900K CPU @ 3.60 GHz；GPU为NVIDIA TITAN V；内存32.00 GB。软件环境均为64位Windows 10专业版操作系统，VS2013 + CUDA10.0平台。实验采用超高阶地球重力场模型EIGEN6C4^[20]，拓展阶次2 160阶，空间分辨率5′。

3.1 精度分析

为验证本文所提方法精度的可靠性，在CPU端使用传统重力场模型串行计算方法，计算共计1 000个飞行轨迹点的扰动引力三分量，作为“真值”。利用CPU+GPU设备，基于上述向量化并行计算公式，计算相同飞行轨迹点扰动引力三分量。将其结果分别与“真值”进行对比，对比结果如表 1所示。

由表 1可以看出，基于向量化的并行算法相较于传统串行算法计算扰动引力三分量的最大误差绝对值仅为10^-6 mGal，平均误差和均方根误差均小于10^-12 mGal，充分证明了本文所提向量化并行算法的精度可靠性。

3.2 效率分析

利用跨阶次递推法计算缔合勒让德函数，不同实验环境下计算不同阶次扰动引力单点计算耗时、加速比和并行比例(可并行部分耗时与总耗时之比)，统计至表 2。

由表 2可以看出，随着阶数的增加，利用跨阶次递推法计算扰动引力的串行算法耗时迅速增加，当截断阶数达到2 160阶时，硬件环境1的单点计算耗时超过1.18 s，硬件环境2也达到0.68 s。在利用向量化公式进行并行计算改进后，并行比例普遍高于92%，充分利用了计算设备的计算能力；且随着阶数的增加，并行比例略有提高。不同阶数扰动引力计算效率普遍提高了一个量级，随着截断阶数的增加，加速比在不断提高，在普通便携式计算机上加速比就可以达到8.11以上，最高可达9.22。而在工作站上加速比普遍高于11.42，最高可达13.20，计算效率提升显著。

利用Belikov递推法计算缔合勒让德函数，不同实验环境下计算不同阶次扰动引力单点计算耗时、加速比和并行比例分别统计至表 3。

将不同实验环境下，依据跨阶次递推法和Belikov递推法进行不同阶次扰动引力计算耗时绘制图 3、4。

由表 2、3和图 3、4可以看出，基于Belikov递推法计算缔合勒让德函数并利用向量化公式进行并行计算，并行比例普遍高于89%，但低于利用跨阶次递推法。不同阶数扰动引力计算在普通便携式计算机上加速比可以达到6以上，最高可达7.26。而在工作站上加速比普遍高于7，最高可达8.99，加速效果不如利用跨阶次递推法明显。且计算相同阶数缔合勒让德函数时，跨阶次递推法相较Belikov递推法耗时较少，并随着阶数的增加，耗时差距逐渐增大。

4 结论

1) 本文对利用地球重力场位系数模型计算扰动引力的效率问题进行了研究。兼顾所用设备效能和计算需要，提出了基于向量化的扰动引力矢量快速并行算法，并推导出利用地球重力场模型计算扰动引力三分量的向量化表达公式。

2) 相比于原有串行逻辑公式，向量化公式减少了计算过程中的累加循环，且符合主流并行计算设备的计算逻辑，更有利于对扰动引力矢量计算进行并行化。

3) 利用CUDA异构并行算法对上述扰动引力三分量向量化计算公式进行并行计算实验，充分验证了本文所提方法的有效性。结果表明，利用跨阶次递推法进行单点计算加速比可达8以上，最高可达13.20；利用Belikov递推法进行单点计算加速比可达6以上，最高可达8.99。