留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

OpenMP并行计算在全球电离层反演中的应用

张强 赵齐乐

张强, 赵齐乐. OpenMP并行计算在全球电离层反演中的应用[J]. 武汉大学学报 ● 信息科学版, 2018, 43(2): 227-233, 240. doi: 10.13203/j.whugis20160065
引用本文: 张强, 赵齐乐. OpenMP并行计算在全球电离层反演中的应用[J]. 武汉大学学报 ● 信息科学版, 2018, 43(2): 227-233, 240. doi: 10.13203/j.whugis20160065
ZHANG Qiang, ZHAO Qile. Application of Parallel Computing with OpenMP in Global Ionosphere Mapping[J]. Geomatics and Information Science of Wuhan University, 2018, 43(2): 227-233, 240. doi: 10.13203/j.whugis20160065
Citation: ZHANG Qiang, ZHAO Qile. Application of Parallel Computing with OpenMP in Global Ionosphere Mapping[J]. Geomatics and Information Science of Wuhan University, 2018, 43(2): 227-233, 240. doi: 10.13203/j.whugis20160065

OpenMP并行计算在全球电离层反演中的应用

doi: 10.13203/j.whugis20160065
基金项目: 

国家高技术研究发展计划 2014AA121501

中央高校基本科研业务费专项 2014618020202

详细信息

Application of Parallel Computing with OpenMP in Global Ionosphere Mapping

Funds: 

The National High-tech R & D Program of China 2014AA121501

the Fundamental Research Funds for the Central Universities 2014618020202

More Information
图(9) / 表(2)
计量
  • 文章访问数:  1299
  • HTML全文浏览量:  101
  • PDF下载量:  372
  • 被引次数: 0
出版历程
  • 收稿日期:  2016-12-28
  • 刊出日期:  2018-02-05

OpenMP并行计算在全球电离层反演中的应用

doi: 10.13203/j.whugis20160065
    基金项目:

    国家高技术研究发展计划 2014AA121501

    中央高校基本科研业务费专项 2014618020202

    作者简介:

    张强, 博士生, 主要从事GNSS电离层及精密定位研究。zhangqiang@whu.edu.cn

    通讯作者: 赵齐乐, 博士, 教授。zhaoql@whu.edu.cn
  • 中图分类号: P352

摘要: 对全球电离层反演数据处理中的计算密集型任务进行分析,针对数据预处理、组建法方程矩阵、参数预消除和法方程矩阵求逆等主要模块设计了基于OpenMP(Open Multi-Processing)的并行计算方案。该实验方案在单台服务器下实施,通过算例验证了本文并行计算方案的有效性和可靠性。实验结果表明:采用并行计算后,全球电离层快速解执行时间只需要10~13 min,计算速度加快了约6倍;最终解执行时间只需要39~47 min,计算速度加快了约5倍。本文全球电离层模型精度约为2.8~3.8 TECU,最终解模型精度相比快速解精度提高了约0.2 TECU,与IGS各个分析中心电离层产品精度基本相当。

English Abstract

张强, 赵齐乐. OpenMP并行计算在全球电离层反演中的应用[J]. 武汉大学学报 ● 信息科学版, 2018, 43(2): 227-233, 240. doi: 10.13203/j.whugis20160065
引用本文: 张强, 赵齐乐. OpenMP并行计算在全球电离层反演中的应用[J]. 武汉大学学报 ● 信息科学版, 2018, 43(2): 227-233, 240. doi: 10.13203/j.whugis20160065
ZHANG Qiang, ZHAO Qile. Application of Parallel Computing with OpenMP in Global Ionosphere Mapping[J]. Geomatics and Information Science of Wuhan University, 2018, 43(2): 227-233, 240. doi: 10.13203/j.whugis20160065
Citation: ZHANG Qiang, ZHAO Qile. Application of Parallel Computing with OpenMP in Global Ionosphere Mapping[J]. Geomatics and Information Science of Wuhan University, 2018, 43(2): 227-233, 240. doi: 10.13203/j.whugis20160065
  • 由于GNSS信号穿过电离层区域会产生与频率相关的色散效应,利用GNSS双频观测数据可以反演全球电离层总电子含量(total electron content, TEC)、监测电离层变化活动以及分析电离层的时空变化规律[1]。1998年,IGS成立了电离层工作组(Ionosphere Working Group, IWG),负责全球电离层TEC建模与监测,并发布了日常的全球电离层快速产品和最终产品[2],国内部分单位也于近些年开始发布全球电离层产品[3-6]。随着GNSS的发展,IGS全球分布测站已达到400多个,且可以提供包括GPS、GLONASS、BDS、Galileo这4个系统在内的观测数据。此外,掩星观测数据、DORIS观测数据、雷达测高数据也为反演全球电离层TEC和电子密度提供了丰富的数据源[7]。更多的GNSS观测数据以及多种观测技术的数据融合,一方面有助于提高全球电离层模型的时空分辨率,另一方面却大量增加了模型复杂度和计算时间,严重降低了全球电离层建模的效率。对于全球电离层模型的优化测试、长时间段(10 a)的数据再处理以及实时电离层建模,数据处理效率是必须考虑的因素之一。

    OpenMP(Open Multi-Processing)并行计算作为一种高效的计算方法在大地测量领域得到广泛应用,显著提高了地球重力场反演[8-10]和大规模GNSS网平差数据处理[11-12]的计算效率,但是在全球电离层建模方面的应用较少。文献[5]初步分析了在测站数量为376个的条件下,利用OpenMP解算1天的GPS观测数据来反演全球电离层,其结果表明, 8核服务器下单天解的执行时间仍然需要60 min左右。文献[6]基于Bernese软件,将全球200多个测站分为11个子网进行处理,其结果表明, 4核服务器下快速解执行时间不分网(串行处理)需要150 min左右,分网(并行处理)需要40 min左右。当前的并行方法没有针对数据处理中的计算任务进行研究,计算速度仍有较大的提升空间。为了进一步提高全球电离层建模的数据处理效率,本文对数据处理中的计算密集型任务进行了深入分析,并提出了基于OpenMP的全球电离层反演并行计算方案。

    • 根据GNSS双频观测数据反演全球电离层TEC的基本表达式为:

      $$ \begin{align} {\bar P}^k_{4,i} &= \xi F\left( z \right)\sum\limits_{n = 0}^{{n_{\max}}} {\sum\limits_{m = 0}^n [{{\tilde P}_{\rm nm}}\left( {\sin\beta } \right)({{\tilde C}_{\rm nm}}\cos\left( {ms} \right) + } \\ &{{\tilde S}_{ nm}}\sin\left( {ms} \right))] + c\left( {{\rm DCB}_{{f_1},{f_2},i} + {\rm DCB}^k_{{f_1},{f_2}}} \right) \end{align} $$ (1)

      式中,β表示电离层穿刺点的地磁纬度; s表示电离层穿刺点的日固经度; $ \tilde P_{nm}$表示正交化缔合勒让德函数; nmax表示球谐函数的最大阶次; $\tilde C $nm和$\tilde S $nm表示球谐函数的待估参数; DCBf1, f2k表示卫星k的差分码偏差; DCBf1, f2, i表示测站接收机i的差分码偏差;ξ=40.28(1/f12-1/f22);F(z)表示卫星在测站处与天顶距z相关的投影函数;P4, ik表示由载波相位平滑伪距得到的非差电离层观测组合。

      在全球电离层建模中,为了分离卫星和接收机的差分码偏差,各电离层分析中心通常采用卫星的差分码偏差“零均值”作为基准[13]。为了更准确地描述电离层TEC在时间尺度上的变化特性,通常在一天内估计多组球谐系数,一般将时段长度设置为2 h,即每天估计13组球谐系数。本文在多组系数之间采用分段线性(piece wise linear, PWL)处理策略[13]。由于IGS观测站在全球的分布不均匀,全球电离层模型存在部分区域TEC的解算异常问题,本文在不同时段球谐系数之间添加随机游走约束[13]

      本文采用附有限制条件的间接平差同时估计球谐函数模型系数和卫星与接收机的差分码偏差,参数估计与精度评估公式为:

      $$ {\left\{ \begin{array}{l} \hat X = {\left( {{\boldsymbol{N}_{\rm BB}} + {\boldsymbol{N}_{\rm DCB}} + {\boldsymbol{N}_{\rm PWL}} + {\boldsymbol{N}_{\rm REL}}} \right)^{ - 1}} \\ \qquad \times \left( {{\boldsymbol{W}_{\rm BB}} + {\boldsymbol{W}_{\rm DCB}} + {\boldsymbol{W}_{\rm PWL}} + {\boldsymbol{W}_{\rm REL}}} \right)\\ {\boldsymbol{D}_{{\rm{\hat X\hat X}}}} = {\hat \sigma }^2_0{({\boldsymbol{N}_{\rm BB}} + {\boldsymbol{N}_{\rm DCB}} + {\boldsymbol{N}_{\rm PWL}} + {\boldsymbol{N}_{\rm REL}})^{ - 1}} \end{array} \right.} $$ (2)

      式中,NBBNDCBNPWLNREL分别表示法方程矩阵和以上3种约束条件组成的约束矩阵; WBBWDCBWPWLWREL分别表示法方程的右项式; $\hat \sigma_0 $表示验后单位权中误差。

      基于球谐函数模型的待估参数的个数为:

      $$ {u = {t_n}{{({n_{\max}} + 1)}^2} + {n_{\rm site}} + {n_{\rm sat}}} $$ (3)

      式中,tn表示时段个数; nmax表示最大阶次; nsite表示接收机的差分码偏差个数; nsat表示卫星的差分码偏差个数。

      本文中球谐函数的阶数为15阶,快速解(单天解)的时段数为13,最终解(3天解)的时段数为37,卫星和接收机的差分码偏差每天估计一组。当GPS观测站为280个、观测到32颗卫星,GLONASS观测站为200个、观测到24颗卫星时,全球电离层快速解和最终解的待估参数的个数如表 1所示。

      表 1  全球电离层快速解和最终解的参数个数

      Table 1.  Parameters Number for GIM Rapid and Final Solutions

      观测系统 快速解 最终解
      GPS 3 640 9 784
      GPS/GLONASS 3 864 10 008
    • 基于地基GNSS的全球电离层建模主要分为以下3个数据处理流程。

      1) 数据下载:从CDDIS网站下载获取全球分布的280多个观测站数据;

      2) 数据预处理:根据原始观测数据提取电离层穿刺点信息,主要包括粗差剔除、周跳探测、载波相位平滑伪距,以及电离层穿刺点地理经纬度、卫星高度角、卫星方位角等信息的解算;

      3) 全球电离层模型系数和卫星与接收机差分码偏差的解算:主要包括组建法方程矩阵、参数预消除、法方程矩阵求逆、高维矩阵乘法,以及参数精度评估与残差编辑等。

      根据全球电离层建模数据处理流程的特点,课题组基于C++语言自主开发了一套GNSS电离层监测与分析软件(GNSS Ionosphere Monitoring and Analysis Software, GIMAS),该软件基于Visual Studio编译器平台开发。在引入OpenMP并行计算前,首先基于Visual Studio的性能分析工具对串行代码进行性能监测,主要统计相关函数调用次数和执行时间,然后逐步优化主要耗时部分的代码,如减少递归、减少内存分配与释放操作等。在串行优化的基础上,对影响计算速度的关键算法部分引入并行计算。

    • OpenMP属于共享内存编程模型的多线程程序设计技术。通过对特定模型算法的分析和并行计算方案的设计,在源代码中加入指导性注释、编译制导指令#pragma omp parallel行来指明程序的并发属性,并在必要之处加入同步互斥以及通信[14]。OpenMP可支持的编程语言包括C/C++和Fortran,支持OpenMP的商用编译器包括Visual Studio和Intel Compiler等,以及开放源码的GCC、Omni、OMPi编译器等。

      OpenMP是基于线程的并行编程模型,遵循fork-join执行模式,其运行过程如图 1所示。进程中的主线程负责执行串行代码,直到进入并行域开始多线程并行执行并行代码;主线程在进入并行域时创建或唤醒一个从线程队列(fork),从线程与主线程一起并行执行并行域代码;当并行域代码执行完之后,从线程被推出或挂起,只有主线程继续执行串行域代码(join)。运行在不同处理器上的线程之间可通过共享变量来实现数据的交换。

      图  1  OpenMP Fork-join并行执行模型

      Figure 1.  OpenMP Fork-join Parallel Model

      并行的性能主要受机器性能和算法性能的影响。机器性能指的是CPU和存储器的一些基本性能指标以及并行通信速度等。算法性能主要体现在算法的加速比和效率上。并行系统的加速比是指对于一个给定的计算荷载,并行算法的执行速度相对于串行算法的执行速度加快了多少倍,可以使用阿姆达尔定律(Amdahl)计算[5]加速比。并行系统的效率计算公式为[8]

      $$ {{S_p} = 1/\left( {1 - \alpha + \alpha /p} \right)} $$ (4)
      $$ {E_p} = {S_p}/p \times 100\% $$ (5)

      式中,α表示并行计算部分所占比例; p表示线程数。当采用p个线程,且α=100%时,Sp的理论值为pEp的理论值为1。

    • 为了验证本文提出的全球电离层建模并行计算方案,选取2015年年积日001-007全球分布的约280个IGS观测站GNSS观测数据对全球电离层TEC进行建模,并对建模执行时间和模型精度进行分析。

      本文实验平台为单台服务器,基本配置如下:操作系统为OpenSUSE12.3;内存大小为256 GB;CPU型号为Intel Xeon(R) E5-2667,CPU主频为2.90 GHz;物理CPU个数为2个,每个CPU核心数为6个,由于开启了Intel超线程技术,逻辑CPU个数为24个。

    • 本文将需要下载的280多个IGS观测站数据分为7个子网,采用多进程模式同时下载。在特定网速环境下,执行时间约为采用单进程模式的1/7。数据预处理执行时间主要集中于粗差剔除和周跳探测部分,传统的方法是逐测站进行数据预处理。由于各个测站观测数据相互独立,引入并行计算后,可以同时对多个测站进行数据预处理操作。

    • 组建法方程矩阵可以归纳为运算BTB,其中B为设计矩阵,B的维数可以记为n×tn为观测值个数,t为参数个数。BTB的计算复杂度为O(t2n),求逆算法的复杂度为O(t3),当n显著大于t时,组建法方程矩阵比矩阵求逆更耗时。对于全球电离层建模快速解,其观测值个数约为100万左右;对于最终解,其观测值个数约为300万左右。本文采用法方程矩阵叠加的方法直接组建法方程矩阵,可以跳过组建设计矩阵的过程且可以节省存储空间,其计算公式为:

      $$ \left\{ \begin{array}{l} \boldsymbol{B} = \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{B}_1}}&{{\boldsymbol{B}_2}}& \cdots &{{\boldsymbol{B}_n}} \end{array}} \right]\\ {\boldsymbol{N}_{BB}} = {\boldsymbol{B}^{\rm{T}}}\boldsymbol{B} = \sum\limits_{i = 1}^n {\boldsymbol{B}^{\rm{T}}_i{\boldsymbol{B}_i}} = \sum\limits_{i = 1}^n {{\boldsymbol{N}_i}} \end{array} \right. $$ (6)

      由式(6)可知,根据观测值组建法方程矩阵时,各观测值之间相互独立,因此可以采用并行计算,如图 2所示。图中Nit×t的矩阵,在采用多线程并行计算时需要保证足够的内存容量。在法方程矩阵叠加的过程中,每个线程都需要和NBB累加并更新NBB,OpenMP采用了数据规约(reduction)的方式或者互斥锁(omp_set_lock)的机制来控制线程之间的冲突和竞争。在组建法方程矩阵的过程中,为了尽量减小计算通信量的计算荷载,采用基于测站的粗粒度[8]并行算法相较于基于观测值的细粒度并行算法更优。

      图  2  组建法方程矩阵的并行计算

      Figure 2.  Parallel Computing for NBB Stacking

      图 3表示不同线程数下组建法方程矩阵的平均执行时间。从图中可以看到,在串行模式下,GPS快速解和最终解执行时间分别需要12.5 min和52.8 min,加入了GLONASS观测数据后,其执行时间分别增加了5 min和9 min。在采用并行计算后,执行时间显著降低,且随着线程数的增加逐步降低,在采用24个线程时,快速解只需要约1.9 min左右,最终解只需要12 min。

      图  3  组建法方程矩阵的执行时间

      Figure 3.  Computing Time for NBB Stacking

    • 在形成法方程矩阵后,由表 1可知,快速解约为3 800维,而最终解则高达10 000维,对于下一步的法方程矩阵求逆无疑是一个时间瓶颈。在估计多种类型的未知参数时,参数预消除(parameter pre-elimination, PPE)和回代求解是一种常用的策略[15]。全球电离层最终解选取中间一天的结果作为输出,本文针对最终解预消除前后一天的参数,降低法方程矩阵的维度,参数预消除后,最终解法方程矩阵维度大小与快速解一致。本文暂不预消除快速解的参数,从而避免了参数回代求解过程,而且可以节约存储空间。

      对于全球电离层模型的未知参数,假设X1, X2, X3, X4, S分别表示时段参数和全局参数,则法方程为:

      $$ \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{N}_{{X_1}{X_1}}}}&{{\boldsymbol{N}_{{X_1}{X_2}}}}&\boldsymbol{O}&\boldsymbol{O}&{{\boldsymbol{N}_{{X_1}S}}}\\ {{\boldsymbol{N}_{{X_2}{X_1}}}}&{{\boldsymbol{N}_{{X_2}{X_2}}}}&{{\boldsymbol{N}_{{X_2}{X_3}}}}&\boldsymbol{O}&{{\boldsymbol{N}_{{X_2}S}}}\\ \boldsymbol{O}&{{\boldsymbol{N}_{{X_3}{X_2}}}}&{{\boldsymbol{N}_{{X_3}{X_3}}}}&{{\boldsymbol{N}_{{X_3}{X_4}}}}&{{\boldsymbol{N}_{{X_3}S}}}\\ \boldsymbol{O}&\boldsymbol{O}&{{\boldsymbol{N}_{{X_4}{X_3}}}}&{{\boldsymbol{N}_{{X_4}{X_4}}}}&{{\boldsymbol{N}_{{X_4}S}}}\\ {{\boldsymbol{N}_{S{X_1}}}}&{{\boldsymbol{N}_{S{X_2}}}}&{{\boldsymbol{N}_{S{X_3}}}}&{{\boldsymbol{N}_{S{X_4}}}}&{{\boldsymbol{N}_{SS}}} \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {{\boldsymbol{X}_1}}\\ {{\boldsymbol{X}_2}}\\ {{\boldsymbol{X}_3}}\\ {{\boldsymbol{X}_4}}\\ \boldsymbol{S} \end{array}} \right] = \left[ {\begin{array}{*{20}{l}} {{\boldsymbol{b}_{{X_1}}}}\\ {{\boldsymbol{b}_{{X_2}}}}\\ {{\boldsymbol{b}_{{X_3}}}}\\ {{\boldsymbol{b}_{{X_4}}}}\\ {{\boldsymbol{b}_S}} \end{array}} \right] $$ (7)

      消除不感兴趣的时段参数X1X4后的法方程为:

      $$ \left[ {\begin{array}{*{20}{l}} {{\boldsymbol{\tilde N}_{{X_2}{X_2}}}}&{{\boldsymbol{\tilde N}_{{X_2}{X_3}}}}&{{\boldsymbol{\tilde N}_{{X_2}S}}}\\ {{\boldsymbol{\tilde N}_{{X_3}{X_2}}}}&{{\boldsymbol{\tilde N}_{{X_3}{X_3}}}}&{{\boldsymbol{\tilde N}_{{X_3}S}}}\\ {{\boldsymbol{\tilde N}_{S{X_2}}}}&{{\boldsymbol{\tilde N}_{S{X_3}}}}&{{\boldsymbol{\tilde N}_{SS}}} \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {{X_2}}\\ {{X_3}}\\ S \end{array}} \right] = \left[ {\begin{array}{*{20}{l}} {{\boldsymbol{\tilde b}_{{X_2}}}}\\ {{\boldsymbol{\tilde b}_{{X_3}}}}\\ {{\boldsymbol{\tilde b}_S}} \end{array}} \right] $$ (8)

      式中,

      $$ \left\{ \begin{array}{l} {\boldsymbol{\tilde N}_{{X_2}{X_2}}} = {\boldsymbol{N}_{{X_2}{X_2}}} - {\boldsymbol{N}_{{X_2}{X_1}}}N^{ - 1}_{{X_1}{X_1}}{\boldsymbol{N}_{{X_1}{X_2}}}\\ {\boldsymbol{\tilde N}_{{X_2}{X_3}}} = {\boldsymbol{N}_{{X_2}{X_3}}}\\ {\boldsymbol{\tilde N}_{{X_2}S}} = {\boldsymbol{N}_{{X_2}S}} - {\boldsymbol{N}_{{X_2}{X_1}}}N^{ - 1}_{{X_1}{X_1}}{\boldsymbol{N}_{{X_1}S}}\\ {\boldsymbol{\tilde N}_{{X_3}{X_2}}} = {\boldsymbol{N}_{{X_3}{X_2}}}\\ {\boldsymbol{\tilde N}_{{X_3}{X_3}}} = {\boldsymbol{N}_{{X_3}{X_3}}} - {\boldsymbol{N}_{{X_3}{X_4}}}N^{ - 1}_{{X_4}{X_4}}{\boldsymbol{N}_{{X_4}{X_3}}}\\ {\boldsymbol{\tilde N}_{{X_3}S}} = {\boldsymbol{N}_{{X_3}S}} - {\boldsymbol{N}_{{X_3}{X_4}}}N^{ - 1}_{{X_4}{X_4}}{\boldsymbol{N}_{{X_4}S}}\\ {\boldsymbol{\tilde N}_{S{X_2}}} = {\boldsymbol{N}_{S{X_2}}} - {\boldsymbol{N}_{S{X_1}}}N^{ - 1}_{{X_1}{X_1}}{\boldsymbol{N}_{{X_1}{X_2}}}\\ {\boldsymbol{N}_{S{X_3}}} = {\boldsymbol{N}_{S{X_3}}} - {\boldsymbol{N}_{S{X_4}}}N^{ - 1}_{{X_4}{X_4}}{\boldsymbol{N}_{{X_4}{X_3}}}\\ {\boldsymbol{\tilde N}_{SS}} = {\rm{ }}{\boldsymbol{N}_{SS}} - {\boldsymbol{N}_{S{X_1}}}N^{ - 1}_{{X_1}{X_1}}{\boldsymbol{N}_{{X_1}S}} - \\ \qquad {\boldsymbol{N}_{S{X_4}}}N^{ - 1}_{{X_4}{X_4}}{\boldsymbol{N}_{{X_4}S}}\\ {\boldsymbol{\tilde b}_{{X_2}}} = {\boldsymbol{b}_{{X_2}}} - {\boldsymbol{N}_{{X_2}{X_1}}}N^{ - 1}_{{X_1}{X_1}}{\boldsymbol{b}_{{X_1}}}\\ {\boldsymbol{\tilde b}_{{X_3}}} = {\boldsymbol{b}_{{X_3}}} - {\boldsymbol{N}_{{X_3}{X_4}}}N^{ - 1}_{{X_4}{X_4}}{\boldsymbol{b}_{{X_4}}}\\ {\boldsymbol{\tilde b}_S} = {\boldsymbol{b}_S} - {\boldsymbol{N}_{S{X_1}}}N^{ - 1}_{{X_1}{X_1}}{\boldsymbol{b}_{{X_1}}} - {\boldsymbol{N}_{S{X_4}}}N^{ - 1}_{{X_4}{X_4}}{\boldsymbol{b}_{{X_4}}} \end{array} \right. $$ (9)

      由式(8)和式(9)可知,采用参数预消除法可以消除最终解前后一天的时段参数,最终解法方程矩阵由10 000维降低为3 800维。参数预消除后的法方程矩阵各部分同样采用并行计算。

    • 基于最小二乘原理建立全球电离层模型的法方程矩阵为对称正定阵,求解法方程可以采用高斯-约当(Gauss-Jordan)算法、乔列斯基(Cholesky)分解算法、改进的Cholesky分解算法[9, 16]。采用Cholesky分解算法可以显著降低参数求解的时间[16],但是在全球电离层建模中需要解算待估参数的方差-协方差阵来评估模型系数的精度,因此,严密的数据分析中矩阵求逆运算不可避免。本文采用列主元Gauss-Jordan算法对法方程矩阵进行求逆运算,并在法方程矩阵求逆中引入OpenMP并行计算,如图 4所示。

      图  4  法方程矩阵求逆的并行计算

      Figure 4.  Parallel Computing for NBB Inversion

      图 5表示最终解参数预消除前后法方程矩阵求逆的执行时间。从图中可以看到,参数预消除前,最终解法方程矩阵求逆在单线程串行模式下的执行时间约为35 min,在多线程并行模式下的执行时间约为10 min左右,其加速比约为3.5。采用参数预消除策略后,最终解法方程矩阵求逆执行时间在单线程串行模式下需要3.8 min,在多线程并行模式下仅仅需要0.8 min。采用参数预消除法,最终解法方程矩阵执行时间提高了约90%。

      图  5  最终解法方程矩阵求逆的执行时间

      Figure 5.  Computing Time for NBB Inversion of Final Solutions

      图 6表示不同线程数下法方程矩阵求逆的平均执行时间。从图中可以看到,在串行模式下,快速解和最终解执行时间分别需要1.8 min和3.8 min。在采用参数预消除和并行计算后,快速解与最终解执行时间基本相当,约为0.8 min。比较图 5中参数预消除前的最终解法方程矩阵求逆时间与图 6中快速解法方程矩阵求逆时间可知,最终解执行时间约为快速解的18倍,与求逆算法的复杂度的理论值基本一致。比较图 3图 6可以发现,组建法方程矩阵执行时间显著大于法方程矩阵求逆执行时间,与算法复杂度的理论值基本一致。

      图  6  法方程矩阵求逆的执行时间

      Figure 6.  Computing Time for NBB Inversion

    • §3.2~§3.5节详细分析了全球电离层建模数据处理中主要模块的并行计算执行时间。图 7表示不同线程下,全球电离层建模的平均执行时间。表 2分别列出了线程个数为1、2、3、4、6、9、12、18和24的情况下的平均建模时间。需要指出的是,本文全球电离层建模数据处理中包括了3次参数估计和残差编辑。

      图  7  不同线程数下全球电离层建模执行时间

      Figure 7.  GIM Parallel Computing Time with Different Threads

      表 2  不同线程个数下的全球电离层建模执行时间/min

      Table 2.  GIM Parallel Computing Time with Different Threads/min

      线程个数 1 2 3 4 6 9 12 18 24
      GPS_Rapid 59.20 33.50 24.28 19.62 15.17 12.20 10.82 9.92 10.45
      GNSS_Rapid 77.10 41.50 29.58 23.98 18.55 14.98 13.90 12.62 12.48
      GPS_Final 192.10 106.17 79.47 67.90 53.17 43.90 42.02 38.68 42.26
      GNSS_Final 221.12 126.55 95.08 81.42 66.42 55.88 50.38 46.93 49.08

      图 7表 2中可以看到,采用OpenMP并行计算后,全球电离层建模执行时间随着线程个数的增加显著降低,而后趋于稳定。在单线程串行模式下,快速解执行时间约60~77 min,最终解执行时间约为192~221 min,约为快速解的3倍。采用并行计算后,快速解和最终解在线程个数增加到12时,执行时间基本趋于稳定,快速解只需要约10~13 min,最终解只需要约39~47 min,最终解约为快速解的3.7倍。最终解相较于快速解,观测数据量约为后者的3倍,显著增加了组建法方程矩阵的执行时间,如图 3所示。由于最终解采用了参数预消除策略,法方程矩阵求逆的时间显著降低,进一步采用并行计算后,执行时间基本和快速解一致,如图 5图 6所示。在采用并行计算后,全球电离层建模的快速解和最终解的执行时间的差异,主要取决于组建法方程矩阵的时间差异,即观测数据量的差异。

      若考虑所有的计算机性能因素对并行计算的影响,则加速比的实际值[8]可以表示为Sp=T1/Tp,其中T1表示串行计算的执行时间,Tp表示p个线程个数的并行计算的执行时间。图 8表示本文计算的全球电离层建模在不同线程数下的并行加速比和并行效率与理论值的比较。本文中全球电离层建模并行代码约占所有代码的92%。图 8中黑色圆点画线表示α=92%时并行加速比与并行效率的理论值。

      图  8  全球电离层建模并行计算加速比和并行效率

      Figure 8.  GIM Parallel Speedup and Efficiency

      图 8中可以看到,并行加速比随着线程个数的增多逐渐增大,当线程个数达到12个以后,基本趋于稳定。采用并行计算后,快速解的加速比约为6,最终解的加速比约为5。当线程个数为2时,并行效率约为90%,当线程个数增加到24个时,并行效率约为20%。

      图 8中可以看到,随着线程个数的增多,实际加速比逐渐偏离理论值。从图 5表 2中也可以看到,当线程个数达到一定数量后,执行时间反而略有升高。其主要原因在于并行计算的实际执行时间为计算时间、并行开销时间和线程间相互通信时间的综合。增加更多的线程会增加并行开销时间,且OpenMP为共享内存的多线程程序设计技术,程序中通常包括临界区,需要采用互斥锁来确保计算的正确性,从而增加了线程间的通信和等待时间。因此,对于一定计算规模的并行计算,并非线程个数越多越好,需要通过实验选取最适合计算任务的线程个数。

      最后,将本文解算的全球电离层模型与IGS最终电离层模型做了对比分析,模型精度如图 9所示。从图 9中可以看到,本文解算的模型精度相较于IGS最终电离层模型精度约为2.8~3.8TECU,与IGS各个分析中心电离层产品精度基本相当。最终解模型精度相比快速解模型精度提高了约0.2 TECU,原因在于最终解消除了第1个时段和第13个时段的边际误差。此外,从图中可以看到,加入GLONASS观测数据后,模型精度略有提高。

      图  9  全球电离层模型精度

      Figure 9.  GIM Global RMS w.r.t IGS Products

    • 为了提高全球电离层反演的数据处理效率,基于自主开发的GNSS电离层监测与分析软件(GIMAS),针对全球电离层反演中数据预处理、组建法方程矩阵、参数预消除和法方程矩阵求逆等主要模块设计了基于OpenMP的并行计算优化方案。实验结果表明:全球电离层模型快速解执行时间由60~77 min降低为10~13 min,计算速度提高了约6倍;最终解执行时间由192~221 min降低为39~47 min,计算速度提高了约5倍。本文全球电离层模型精度约为2.8~3.8 TECU,最终解模型精度相比快速解精度提高了约0.2 TECU;增加GLONASS观测数据后,模型精度略有提高。

      本文的实验平台为高性能服务器,在并行计算设计中并未充分考虑内存消耗的问题。假如在内存有限的PC机平台上采用全球电离层建模最终解并行方案,可能会导致内存溢出。笔者将在后续工作中进一步优化适合并行计算的矩阵存储等问题,以及结合MKL、LAPACK等并行函数库进一步优化高维矩阵运算。

参考文献 (16)

目录

    /

    返回文章
    返回