由地球重力场模型计算扰动位、重力异常、高程异常等重力场参量的方法称为球谐综合法。当采用高阶次地球重力场模型串行计算大规模格网时，球谐综合法的计算工作量大，耗时长。为加快数据处理速度，不少学者从算法上进行了优化，在一定程度上提高了计算效率^[1-5]。并行计算是高性能计算领域使用最为广泛的一种技术^[6]，它的基本思想是采用多个处理器来协同求解同一个问题。利用并行计算技术突破大规模重力数据处理的瓶颈、提高计算效率已逐渐成为一项共识^[7]。典型的并行编程环境主要有MPI(message passing interface)、OpenMP(open multi-processing)、HPF(high performance Fortran)、CUDA (compute unified device architecture)等^[8-9]，其中MPI、OpenMP是物理大地测量领域较为常用的两种并行模型^[10-16]。

MPI、OpenMP分别为基于消息传递和共享内存的并行编程模型^[17]，其中OpenMP具有易编程、内存开销小等特点，但仅能在1个计算节点内并行；MPI采取分布式内存并行，计算作业可以在1个或者多个节点内并行，适用于个人计算机、服务器、计算机集群、超算等；而MPI具有移植性强、拓展性好、通信函数丰富等优点，是目前可信赖的高效率大规模并行计算平台^{[16, 18]}。因此，本文将适用范围更广、并行潜能更大的MPI运用于大规模格网重力场参量的球谐综合计算中。

1 基本原理

地球重力场一般可通过扰动重力场进行研究，地球外部任意一点的扰动位T的球谐级数可表示为：

$ T\left( {r,\theta ,\lambda } \right) = \frac{{GM}}{r}{\left( {\frac{R}{r}} \right)^n}\sum\limits_{m = 0}^n {\left( {\bar C_{nm}^*\cos \left( {m\lambda } \right) + \bar S_{nm}^ * \sin \left( {m\lambda } \right)} \right)} {\bar P_{nm}}\left( {\cos \theta } \right) $

(1)

式中，$\bar P_{nm}^{}\left( {{\rm{cos}}\theta } \right) $为完全规格化的缔合勒让德函数；其余变量含义可参考文献[1]。常用的缔合勒让德函数计算方法主要有标准前向列推法、标准前向行推法、Belikov列推法和跨阶次递推法，其中跨阶次递推法计算量大、精度高、稳定性好^{[4, 19-20]}，故本文采用跨阶次递推法计算缔合勒让德函数：

$ \begin{array}{*{20}{l}} \begin{array}{l} {{\bar P}_{nm}}\left( {{\rm{cos}}\theta } \right) = {\alpha _{nm}}{{\bar P}_{n - 2, m}}\left( {{\rm{cos}}\theta } \right) + \\ \;\;\;\;\;\;{\beta _{nm}}{{\bar P}_{n - 2, m - 2}}\left( {{\rm{cos}}\theta } \right) - {\gamma _{nm}}{{\bar P}_{n, m - 2}}\left( {{\rm{cos}}\theta } \right), {\rm{}} \end{array}\\ {\;\;\;\;\;\;\forall n > m \ge 2} \end{array} $

(2)

式中, α_nm、β_nm、γ_nm为递推系数，具体计算公式参见文献[3]。

对式(1)中的扰动位T进行简单的数学变换可以得到重力异常、高程异常、垂线偏差等地球重力场参量。以重力异常Δg为例，将正常椭球面做球近似假设，则T与Δg之间的关系式可表示为：

$ {\rm{\Delta }}g = - \frac{{\partial T}}{{\partial r}} - \frac{2}{R}T $

(3)

将式(1)代入式(3)，可得到Δg的球谐展开式：

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\Delta g\left( {r,\theta ,\lambda } \right) = \\ \frac{{GM}}{{{r^2}}}\sum\limits_{n = 2}^N {\left( {n - 1} \right)} {\left( {\frac{R}{r}} \right)^n}\sum\limits_{m = 0}^n {\left( {\bar C_{nm}^*\cos \left( {m\lambda } \right) + } \right.} \\ \;\;\;\;\;\;\;\;\;\;\;\;\left. {\bar S_{nm}^*\sin \left( {m\lambda } \right)} \right){{\bar P}_{nm}}\left( {\cos \theta } \right) \end{array} $

(4)

2 串行算法

方案1  在球谐综合中，勒让德函数$ {\bar P_{nm}}\left( {{\rm{cos}}\theta } \right)$的计算占整个计算量的大部分^[20]。由于$ {\bar P_{nm}}\left( {{\rm{cos}}\theta } \right)$与纬度相关，在利用球谐函数计算格网重力场参量时，为避免$ {\bar P_{nm}}\left( {{\rm{cos}}\theta } \right)$重复计算，一般逐纬度计算，在计算同一纬度的重力场参量时，仅一次计算$ {\bar P_{nm}}\left( {{\rm{cos}}\theta } \right)$。

方案2  由式（4）可知，cos(mλ)、sin (mλ)两项与经度相关。为提高计算效率，将cos(mλ)、sin(mλ)预先算好并存放到数组中，再通过调用的方式参与计算。此项措施联合方案1在本文中标记为方案2。

假定全球纬度带条数为M，经度带条数为2M，参考阶次取为N，则这两项的计算次数为${\rm{Num}} = {M^2} \times \left( {2 + N} \right) \times \left( {N - 1} \right) $；而通过数组预存的方式，计算次数仅为${\rm{Nu}}{{\rm{m}}^{\rm{*}}} = 2MN $次，其运算压缩比为$ Z=M(2+N)(N-1)/2N$。

方案3  由文献[3]可知，缔合勒让德函数的3个系数α_nm、β_nm、γ_nm与经纬度无关，只与阶次相关。由于同一纬度所有阶次的${\bar P_{nm}}\left( {{\rm{cos}}\theta } \right) $值通过递推的方式一次性求得，因此这3个系数在同一条纬度带上并没有重复计算，但在不同纬度带上进行了重复计算。为此，在方案2的基础上，进一步对α_nm、β_nm、γ_nm进行数组预存，记为方案3。不难发现，在α_nm、β_nm、γ_nm的处理上，方案3相对于方案2的运算压缩比为Z=M。

3 并行设计与实现

MPI可以通过分解任务实现多进程并行计算。由式（4）可知，在逐纬度计算时，纬度条带之间的计算是相互独立的，利用MPI将大规模格网重力异常按纬度条带分配给各个进程，在计算结束后收集各个进程的结果，避免了额外的通信开销。

在进行任务分解时有两种模式可选，一种是按纬度分块分解，即每一个进程负责处理一块连续的纬度带，如图 1所示。另一种是按纬度循环分解，以进程P₀为例，假设进程个数为n，则进程P₀负责计算的纬度带为W₁，W_1+n…W_1+kn。在任务无法均等分配时，两种模式的进程任务量有所差异，但最高任务量相同。由于MPI的数据收集函数MPI_Gather()是按进程编号由小到大进行收集，即先将P₀中的所有计算结果存放到一维数组中，然后逐个存放其他进程的数据，当采取按纬度分块分解模式时，收集到的结果可直接输出；而按纬度循环分解模式收集的结果需要进行复杂的顺序变换才可以输出，本文采取按纬度分块对任务进行分解。

由于MPI采用分布式内存数据存储模式，随着进程数P增加，各项系数运算压缩比下降Z/P，内存耗用倍增，尤其是方案3耗用内存较多，在一定程度上会影响计算性能，为了获得最优方案组合，对方案1~3均进行并行分析。

4 算例与分析

本实验采用EGM2008重力场模型，取至2 160阶次，并设置了两组算例。算例1计算了全球格网1°×1°的重力异常，算例2计算了跨纬度15°、跨经度360°的5'×5'格网的重力异常(即格网个数为180×4 320)。

本实验计算平台之一为广东工业大学3S技术实验室配置的DELL PowerEdge R730服务器，该服务器配置有2个物理中央处理器（central processing unit，CPU）（E5-2630L，主频1.80 GHz，每个CPU 8核，共16核），32个逻辑CPU，内存32 GB，内存频率2 133 MHz，搭载Red Hat Enterprise Linux server release 6.7操作系统。程序并行环境为MPI 3.0.4，编译器为GCC 4.4.7，编程语言为C++。

为评价并行算法和数组预存再调用方法对计算过程的优化作用，引入加速比S_P：

$ {S_P} = \frac{{T_1^S}}{{T_i^P}} $

(5)

式中，P(P≥1)表示进程个数；$ T_1^S$为方案1的串行执行时间；$ {T_i^P}$为方案i(i =1, 2, 3)采用P个进程时的执行时间。本文中的加速比皆是相对于方案1的串行耗时。

按照上文3种方案，分别实现串行及并行运算。当进程个数相同时，3种方案的CPU使用率一致，皆为3.125P%（P为进程数，3.125P≤100，而3.125%实际上是1个逻辑CPU占总逻辑CPU个数（32）的百分比，这说明了逻辑CPU的启动个数与进程个数相等（在进程个数不超过逻辑CPU个数的情况下）。表 1为两个算例中3种方案对应进程数的内存耗用情况。由表 1可知，在内存充足的情况下，3种方案的内存耗用与进程数P存在良好的线性关系，算例1分别约为180P、193P、268P，算例2分别约为190P、332P、410P。在算例1中，当进程个数增加到180时，3种方案均遭遇内存不足（内存空间31 714 MB），硬盘交换空间分别为3 316 MB、10 010 MB、19 609 MB；在算例2中，方案3在开启90个进程时，硬盘交换空间为6 944 MB，在开启180个进程时，3种方案的硬盘交换空间分别为6 685 MB、29 297 MB、43 008 MB。

表 2为算例1、算例2中3种方案串行及并行的运算时长。由表 2可知，在串行计算中，通过将系数预存调用的方式，可大幅度提高计算效率，方案2、3相对方案1在算例1中的加速比分别为1.58、1.70，在算例2中分别为1.44、1.47。并行计算中，3种方案均在进程数为30时获得最短计算时长，在算例1中由19 260 s分别减少至1 169 s、1 024 s、978 s，最大加速比分别为16.48、18.81、19.69；在算例2中由274 533 s分别减少至14 226 s、12 695 s、12 402 s，最大加速比分别为19.30、21.63、22.14。

为了直观地比较3种方案的并行性能，图 2展示了3种方案的并行加速比。由图 2可知，在算例1和算例2中，各方案的加速比在进程数为30时达到峰值。随着进程数进一步增加，加速比开始缓慢下降，原因在于实验服务器共拥有32个逻辑CPU，而每个逻辑CPU同一时刻只能运行1个进程，当进程个数超出逻辑CPU的个数时，进程之间会出现替换使用逻辑CPU的情况，此时继续增加进程并不能进一步缩短计算耗时；相反，计算资源的过度分散将会降低计算性能。另外，进程数越多，并行的开销时间也越多，反而会拖慢计算速度，因此并行进程的最优个数为逻辑CPU的个数。

由图 2可知，整体上方案3的加速比优于方案2，方案1最差。但当进程数增加到90或180以后，方案1的加速比优于方案2和方案3，其原因一是方案2和方案3的预存系数运算压缩比随着进程数增加而减小（每个进程需要对方案2或方案3所述的预存系数进行独立预存）；二是结合表 1可知，当进程数增加到90或180后，方案2、方案3面临内存严重不足问题，进程部分内存被释放到swap分区中，频繁的页面调度和交换增加了磁盘负载，大大降低了计算效率，而方案1的swap分区交换内存相对于方案2、方案3则要小得多，故其加速比、相对效率反而优于方案2、方案3。

为进一步验证上述结论，本文也在中山大学超算平台“天河二号”上进行了实验。“天河二号”每个节点配置两颗主频为2.2 GHz的Intel Xeon E5-2692 V2系列CPU，共24核，48个逻辑CPU，内存为128 GB，内存频率为1 333 MHz，并配有3个Xeon Phi 57核心的协处理器（运算加速卡）。程序并行环境为MPI 3.1，编译器为GCC 4.4.7，编程语言为C++。值得注意的是，“天河二号”可以根据进程个数动态分配节点，并始终满足分配的逻辑CPU个数大于等于进程个数。

由表 3可知，由于“天河二号”节点的计算资源充足，在两个算例中，3种方案的计算耗时均随着进程数的增加一直减少。图 3中算例1中3种方案之间的差距很小，方案2、方案3略优于方案1；算例2中方案1则明显要优于方案2和方案3，其加速比约为方案2和方案3的两倍。这是由于“天河二号”每一节点的计算速度快(达422.4 Gflops)，并且配有运算加速卡，其运算的瓶颈并不在于CPU，因此通过预存数组（即减少CPU运算工作量）所减少的计算耗时很小；相反，随着格网分辨率的提高，预存数组急剧增大，庞大的数组无法在高速缓存器中完全驻留，CPU在调用数组时需要频繁访问内存。

由上文可知，“天河二号”内存频率与CPU主频的比值为0.606。因此在算例2中，其访问内存增加的时间要远多于CPU减少的时间；而在DELL服务器上，内存频率与CPU主频的比值为1.185。两个算例中预存数组在DELL服务器上均取得了良好效果。

5 结语

本文以2 160阶次的EGM2008模型为例，基于MPI实现了按纬度带分解并行计算格网重力异常，通过对MPI并行、MPI结合数组预存调用等方案的运行效率进行分析比较，得出如下结论。

相比串行计算，MPI并行计算能充分利用多核计算机资源，大幅度缩短计算时长；增加进程数可以进一步分解计算任务量，但由于1个逻辑CPU同一时刻只能运行1个进程，当进程个数超出逻辑CPU个数时，各进程需要替换使用逻辑CPU，故在此基础上继续增加进程已失去意义，相反，进程数越多，并行的开销时间也越多，反而会增加整体计算耗时。因此建议开启的进程个数取为逻辑CPU的个数。此外，应尽量避免内存耗用超出计算平台内存。

对cos(mλ)、sin (mλ)以及缔合勒让德函数系数α_nm、β_nm、γ_nm采取数组预存调用的方式可有效减少CPU的运算工作量，但会增加内存的访问次数，建议在CPU性能一般而内存性能相对较好的计算平台上采用，CPU与内存的性能衡量可参照文中的两个计算平台。

本文中提出的数组预存调用方式以及逐纬度分块并行计算模式适用于采用球谐综合法计算扰动位、重力扰动、垂线偏差、高程异常、大地水准面高等重力场参量。