高性能法方程求逆方法：GNSS基准站网数据处理方法与应用

在基于最小二乘方法的参数估计中，需要逐历元读取观测数据并叠加法方程，至观测时段历元结束时形成最终法方程并加以求逆，得到待估参数的最小二乘估值。因此，法方程解算(主要是法方程求逆过程)是最小二乘参数估计中非常重要的环节，也是影响数据处理效率的关键因素之一。为了实现大规模基准站网的高性能数据处理，需要解决两个关键问题：在特定维数下提高法方程解算效率和减小法方程的维数。后者将在第4章中进行详细研究，本节详细探讨研究提高法方程解算效率的方法。

1.高性能线性代数函数库

法方程求逆过程是一个标准的线性代数问题，可以用任意标准方法进行解算。在最小二乘估计方法中，法方程系数矩阵为一正定对称矩阵，因此高斯消元法、LU分解法、Bunch-Kaufman选主元法以及乔里斯基分解法均可用于法方程求解。从原理上说，不同方法求解效率并不一致。对于大规模基准站网的高性能数据处理，需要找出与最小二乘估计相关的最优法方程求逆算法。

此外，对于法方程求逆问题，不仅运用不同的方法求逆效率有差异，即便同一个方法、不同编程实现方式的解算效率也可能相差较大。因此，对于非数值计算的专业人员而言，面对这样一个标准的数值计算问题，要使自行编写的函数具有较高的计算机数值解算效率，需要投入精力与时间进行反复调试。

具体来说，对于高斯消元法这样一个相对简单的算法，在编程时主要问题是如何提高内存读/写的效率。目前计算机处理器具有分级的存取速度设计，如L1缓存、L2缓存等，其中可以直接被处理器使用的部分缓存速度最快(通常来说只有几百KB)，而普通内存虽然容量更大，但速度慢了几个量级。因此，如果计算操作可以在处理器缓存中完成，那么执行效率较高。但处理器缓存毕竟容量有限，为了避免溢出常常需要在处理器缓存和内存之间交换数据。这些操作会大大地影响计算效率，因为处理器的计算速度比内存读取速度快几百，甚至上千倍。可见，如何尽量减少对内存等慢速存储器的访问，决定了数值计算的效率。同时，现有的高性能线性代数函数库，包括BLAS、LAPACK、MKL等，其运算效率已经被很多学者加以验证(Anderson，et al.，2002；Demmel，1989)。因此，使用这些已有的函数库，可以将精力放在与GNSS相关的算法研究上，而不必在与计算机底层硬件相关的计算效率问题上花很多时间。

(1)BLAS

BLAS(Basic Linear Algebra Subroutines，基本线性代数子程序库)是一个高质量的基本向量、矩阵运算子程序库。最早是由FORTRAN语言编写，后来又发展了其他语言接口，包括C语言、JAVA等。BLAS仅涉及最基本的向量和矩阵运算，因此可将高性能数值计算程序的开发同特定硬件平台的性能优化区分出来。开发者只需将所涉及的计算过程转换为基本的矩阵、向量运算并调用相应的BLAS子程序即可，而不必考虑与计算机硬件相关的性能优化问题。这一模式提高了高性能线性代数函数库的开发效率，如LAPACK、ScaLAPACK等都是基于BLAS而设计的(张林波，2006)。

BLAS从结构上可以分为三个层次：Level 1 BLAS、Level 2 BLAS和Level 3 BLAS。其中Level 1 BLAS涉及向量与向量、向量与标量之间的运算，Level 2 BLAS涉及向量和矩阵之间的运算，Level 3 BLAS则涉及矩阵之间的运算。一般来说，层次越高的BLAS子程序效率越高(张林波，2006)。因此，应尽可能地使用Level 3 BLAS的子程序。

(2)LAPCAK

LAPACK(Linear Algebra PACKage，线性代数软件包)是由美国Argonne国家实验室、Courant研究院和NAG公司联合开发完成的一个大型线性代数函数库，以FORTRAN语言编写。第一个版本发布于1992年2月(Anderson，et al.，1999；Anderson，et al.，1995)。LAPACK包含了求解科学与工程计算中最常见的线性代数数值计算问题，如线性方程组求解、矩阵奇异值求解、特征值求解等，还可以实现各种矩阵分解、矩阵求逆、条件数估计等相关计算问题。LAPACK通过对矩阵进行分块，将许多操作转换为矩阵运算，调用Level 3 BLAS的高效子程序来完成。这样有助于提高处理器缓存的命中率，从而进一步提升了运算效率。

LAPACK是一个开放源代码的工程，包含了LAPACK所有源代码的程序安装包可以从网址http：//www.netlib.org/lapack/免费获得，其软件包结构如图3-21所示(张林波，2006)。

(3)MKL

图3-21　LAPACK软件包结构图

英特尔数学核心函数库(Intel Math Kernel Library，MKL)是英特尔公司为科学和工程计算设计的数学函数库，并针对目前和将来的英特尔公司的处理器(包括奔腾系列、至强系列、安腾系列、酷睿系列等英特尔公司全系列处理器)进行了优化，特别针对英特尔至强、酷睿系列多核处理器进行了多线程的性能优化(Gustafson，et al.，2007；Basics，2005)。MKL中包含了BLAS库和LAPACK库，并已集成在英特尔FORTRAN编译器ifort中，因此使用时可直接调用LAPACK子函数，而不必显式链接LAPACK库。

为了研究分析适用于最小二乘估计方法的最优法方程解算方法，分别利用高斯消元法、LU分解法、Bunch-Kaufman选主元法和乔里斯基分解法来解算法方程，并采用2009年5月2日的IGS跟踪站网数据(采样率为300s)，分别在50站、100站、150站和200站的不同测站数目情况下进行解算。所用数据如图3-19所示，解算时的计算机软、硬件平台见表3-11。采用上述四种不同方法解算法方程的效率见表3-12。其中，高斯消元法是自行编写的子程序，LU分解法采用了LAPACK/MKL库中的DGETRF、DGETRI子程序，Bunch-Kaufman选主元法采用了LAPACK/MKL库中的DSYTRF、DSYTRI子程序，乔里斯基分解法采用了LAPACK/MKL库中的DPOTRF、DPOTRI子程序。

表3-12　不同方法法方程求解效率(单位：秒)

由表3-12可知，对于最小二乘估计法方程的解算，高斯消元法所用时间较多，接近LU分解法所用时间的20倍。LU分解法、Bunch-Kaufman选主元法和乔里斯基分解法所用时间基本在同一量级，而乔里斯基分解法所用时间最少，大约为LU分解法所用时间的一半。可见，上述情况一方面说明了自行编写的法方程求解函数与专业函数库相比效率过低，另一方面也说明了高斯消元法的计算量较大而导致与LU分解法等方法相比解算时间过长。

上述四种不同测站数目情况下法方程的维数如图3-22所示。

图3-22　不同测站情况下法方程维数示意图

结合表3-12和图3-22，可以看到，对于法方程解算，乔里斯基分解法的效率较高，在200个测站、法方程维数高达12219的情况下，其解算时间694s，仅为同等情况下高斯消元法的2.4%。原因在于其与最小二乘估计法方程特征的匹配。最小二乘估计法方程的系数矩阵，是一个实正定对称矩阵，乔里斯基分解法可以充分利用其正定对称特性，而LU分解法(适用于所有可逆实矩阵)和Bunch-Kaufman选主元法(适用于所有可逆实对称)则不能充分利用这一特征，因此乔里斯基分解法的解算效率最高。

因此，对于大规模GNSS基准站网高性能数据处理中的法方程解算，可以得到如下结论：

①高斯消元法由于计算量较大，并不适用于大规模基准站网的高性能数据处理；

②对于法方程解算这样的标准线性代数问题，建议采用BLAS/LAPACK/MKL等高性能线性代数函数库以提升解算效率；

③鉴于最小二乘估计法方程的正定对称特性，建议在求解法方程时使用乔里斯基分解法以提高解算效率。

2.基于OPENMP的高性能并行计算(https://www.xing528.com)

上一节研究了大规模GNSS基准站网数据处理中的法方程解算方法与算法实现的优化问题。然而，数据处理仅占用了计算机处理器的一个线程(thread)，即所谓的单线程解算。目前的计算机技术已经进入了多核处理器时代，Intel和AMD公司面向个人消费者推出了众多款式的多核桌面/移动处理器。个人桌面计算机的处理器已经普遍采用双核处理器，部分已经使用了高端的四核甚至六核处理器；个人移动(笔记本)计算机也已经普遍使用了双核处理器，部分使用了四核处理器。特别是Intel公司在推出的酷睿2系列桌面/移动处理器中，使用了“超线程”(Hyper-Threading)技术，可将一个物理内核扩展为两个线程使用，使得双核处理器可以使用四个线程，而四核处理器可以使用八个线程，大大扩展了多核处理器的应用。事实上，目前上述两家公司已经停止了单核处理器的生产，市场上可以购买到的个人桌面/移动计算机的处理器已经全部是多核处理器。因此，在进行大规模GNSS基准站网的高性能数据处理时，可考虑同时使用多个线程进行并行计算，以提升解算效率。

一般来说，并行计算是由并行计算机来实现的。目前主要的并行计算机包括多计算机系统(Multi-Computer)和集中式多处理器系统(Centralized Multi-Processor)。多计算机系统是由多台单核或多核计算机通过网络互联组成的并行计算机系统，在多个计算机上的处理器之间通过传递消息(Message)来通信。目前国际上流行的超级计算机系统大多是按照这一标准组建的，由于构建多计算机系统在实现上较为复杂且成本较高，本节内容暂不涉及。集中式多处理器系统同时也称为对称多处理器系统(Symmetrical Multi-Processor，SMP)，是集成度更高、更加紧密的并行计算机系统，在该系统中所有处理器共享物理内存，并通过共享内存来实现多个处理器之间的通信和同步。本节内容所涉及的多线程并行计算，就是集中式多处理器系统的一个实现方式，在单台计算机中同时使用多个物理内核或线程，以达到并行计算的目的。

从某种意义上说，目前所有具备多核处理器的个人计算机均可视作集中式多处理器并行计算机。并行算法在并行计算机上需要通过并行编程来实现，当前使用较多的并行编程环境主要有消息传递编程接口(Message Passing Interface，MPI)以及共享存储编程接口(Open Multi-Processing，OPENMP)两种。

(1)消息传递编程接口MPI

MPI是目前国际流行的用于并行计算的消息传递标准，被大多数商用并行计算系统采用。MPI 1.0版于1994年5月推出，2009年4月更新，免费的MPI库MPICH可从Argonne国家实验室网站(hhtp：//www-unix.mcs.anl.gov/mpi/mpich)获取。MPI提供了C/C＋＋、FORTRAN 77/90四种编程语言的接口。

MPI定义了一组实现消息传递模型的函数，并集成了不同消息传递函数库中业已被证明为最有效的部分，形成了功能丰富的并行函数库，具有良好的可移植性、扩展性、易用性，并具有完备的异步通信功能。MPI的发展，为并行计算软件产业的发展打下了坚实的基础。

虽然MPI已被证明为并行计算使用最广泛的编程接口之一，但它的优势主要体现在由多个节点(计算机)组成的大型多计算机系统上。对于单台多核计算机而言，使用MPI需要为每个物理处理器内核或线程单独分配独占的内存资源如图3-23所示(Quinn，2003)，实现起来较为繁琐。同时由于单台计算机的物理内核数目并不多(不超过10个)，MPI的优势并不能充分体现。更为重要的是，使用MPI进行并行编程时，需要对已有的串行(单线程)代码进行较多的修改，因而在单台计算机上不建议采用MPI进行并行编程计算。

图3-23　消息传递模型示意图

(2)共享存储编程接口OPENMP

与MPI不同，OPENMP是专为共享存储环境下进行并行计算而设计的编程接口，即所有处理器内核/线程共用一个本地内存资源，而MPI标准下每个处理器/内核/线程有自己的内存并只能读取本地内存上的数据。共享存储模型示意图如图3-24所示。

图3-24　共享存储模型示意图

OPENMP支持FORTRAN语言的第一个标准公布于1997年10月，支持C/C＋＋语言的标准于1998年11月发布，目前OPENMP可支持包括FORTRAN 77/90，C/C＋＋语言。

基于OPENMP的共享存储并行模式称为“fork/join”模式。如图3-25所示，当程序开始执行时，只有一个主线程存在。主线程按串行模式执行程序代码，当进入并行计算区域时，主线程派生出若干线程并协同工作(即fork)。在并行代码结束时，派生的线程销毁或挂起，程序回到单一的主线程中(即Join)，直至程序执行完毕。

基于OPENMP的并行计算与基于MPI的并行计算的区别在于MPI中所有进程存活于整个程序的执行过程之中。而在OPENMP中，在程序开始与结束时仅有唯一的主线程，程序在执行过程之中线程数会根据情况而动态变化。在基于OPENMP的并行计算算法实现时，可以逐步从程序最费时的模块开始进行并行化改造，然后依次对其他模块进行评估并根据情况进行改造。这样，可以在不大规模改动原有代码的基础上完成串行程序代码到并行代码的改造，节约了编程时间。

图3-25　共享存储模式示意图

MPI和OPENMP是常用的两种并行编程环境。它们都可以应用在多计算机系统上，其中MPI适用于多台计算机组成的集群系统中，而OPENMP由于在多个计算机内核/线程间共享本地内存。因此，并不适合这种无共享内存的多节点集群系统。MPI和OPENMP之间的比较见表3-13。

表3-13　OPENMP与MPI的比较

由表3-13可见，OPENMP与MPI两者各有特点，MPI更适合于多台计算机节点组成的计算机集群系统，而OPENMP更适用于多个计算机内核/线程的单台计算机(节点)。因此，在由多个具有多处理器/内核的计算机节点组成的集群系统中，可以采取OPENMP和MPI混合编程的模式来实现并行计算，即在集群的每个节点上创建一个MPI进程，而在单个节点内使用OPENMP实现共享存储并行计算。一般来说，使用OPENMP和MPI混合编程的效率比仅用MPI编写程序的要高。

综上所述，由于本节内容的研究重点是个人计算机平台上的大规模GNSS基准站网数据处理，因此暂不涉及基于MPI的并行计算。我们在笔记本计算机上实现了基于LAPACK/MKL的OPENMP并行化计算，所用的计算机软、硬件平台见表3-11。

在将3.5.2节中的四种最小二乘法方程求解方法进行并行化改造之后，采用2009年5月2日的IGS跟踪站数据进行了算例分析，所用测站数目分别为50站、100站、150站与200站，采样率为300s。测站分布如图3-19所示。3.5.2中通过算例分析，得出对于最小二乘估计而言，利用乔里斯基分解法进行法方程求解效率较高，乔里斯基分解法在并行化之后的解算效率如图3-26所示，解算时同时使用了四个线程进行并行计算。

图3-26　基于OPENMP的乔里斯基分解法解算时间(单位：秒)

结合表3-12和图3-26可知，在实现了并行化计算后，解算效率提升超过400%，在200个测站、法方程维数为12219的情况下，乔里斯基分解法求解法方程仅用130秒，仅为串行计算时694秒解算时间的18%，可见并行计算对解算效率的巨大影响。

综上所述，在大规模GNSS基准站网法方程解算中，可以采用基于OPENMP的并行化乔里斯基分解方法，充分利用多核/多线程处理器的优势进行法方程求逆解算，从而大幅度提高个人计算机平台的解算效率。