HPSEPS 并行软件包及其在地震层析成像中的应用

1、年伞围高性能计算学术年会并行软件包及其在地震层析成像中的应用赵永华迟学斌刘俊（中国科学院计算机网络信息中心超级计算中心北京（）（：）摘要：基于自主开发的并行软件包阵分解方法，在数据压缩、最小二乘问题、数亨图像处理和特征提取等领域有着广泛的应用，其中的一个典型应用是地震层析成像。自年等人的地震成像的开创性作以来，体波地震成像方法得到了长足的发展，在模型的参数化方面，年刘福出等“在模型中引入了间断面用以描述莫霍面、等速度间断面，这一方法被广泛应用于国内的地震层析成像：作；年赵大鹏等改进了射线追踪方法，使其适用于模型中的非水半间断面，进一步接近了实际情况。在反演方法方面，等最初采用阻尼最小二乘法；随

2、后各种叠代类算法被用于地震层析成像反演，如方法，方法等；年和”提出方法后，该方法被广泛用于地震层析成像反演中。“；年刘福出“引入正交投影算子，将震源参数与速度参数分离，提顺序正交三角化算法，减少了内存占用量和计算量，使得最终的方程可用法求解作为“中科院五信息化专项项目：超级计算环境建设”课题的子课题，五期间超级计算中心同地质与地球物理研究所的科学家们协同合作，完成了“三维地震成像软件包”，求解和解决了地震层析成像中的一些实际问题，得到了一些应用成果）（）所提供的奇异值分解（）并行求解器，给出了该并行求解器在天然地震层析成像的并行反演计算中的应用以及“地球层析成像并行求解平台”的框架模型，在深腾

3、平台上对两个地区的实际地震层析成像数据进行了效率测试，给出了性能测试结果、关键词：地震层析成像分解并行计算转换分解算法（，）：（）：，。为了对这一并行计算半台进一步进行整合和”改善，基于自主开发的软件包”，我们开发了一个：行求解器，通过该并行求解器可更高效地并行求解存在于地震层析成像的分解：，法进行地球层析成像反演为了从整体上对地球层析成像有一个完整的，引言基金项甘基金项：刊家然科学基金项甘（，）；刊家“九匕”重电基础毋充发展规划基矩阵的奇异值分解（）作为一个基本的矩金项目（）并行软件包及其在地震层析成像中的应用描述，下面给出与这一问题有关的主要计算。用于速度罔像重建的层析成像法可归结为一个矩

4、阵方程”】。该矩阵方程对于区内地震为：为矩阵，韵第卅、列向量。通常，在没有确定误差与分辨之问半衡点的情况下，所有“都有可利用的价值，我们一般将各阶解少合并成解矩阵：（）其中，一：，。，（一，：“，）谚。巧：，“占，按照反演理论，凶：，×：则对分辨矩阵有、，其中为分辨矩阵，为矩阵，的前矽组成的子阵。若，矗”，则”，模型的每一参伍，：），。）数都可唯一确定，所得解也完全分辨。各阶解都将对应一个分辨矩阵，被称作分辨阶矩阵式中为模型网格点上的速度扰动量，为相应的速度偏导数矩阵，为研究区域内地震震源位置的扰动量，为相应的震源参数偏导数矩阵，为走时残差。式中变量的下标表示地震序号，。为区内地震总

5、数，引入正交投影算子。使得，可将（）式分解为两个方程分别求解：（）（，一），（）（）。（“，”）（）可知（）其中卢为刀“的对角元索组成的列向量，”，“（）（）其中是，的元氛，对解向量中每个元索的方差有缈（点一）项，直接得到形如（）式的方程：万对于远震数据，采用平滑处理后可忽略震源（）眩却芝的方差矩阵，。有）其中肛口，则由解矩阵，列应卅而矩阵口写成为：却州零勺，十辜（，（）其中。：。上面所有这些构成了地球层析成像反演的主要计算，而所有这些计算中，占据着计算的大部分时间，凼而也是地球层析成像反演的关键。凶此需要有效的算法和并行处理方法。，×，矩阵（通常），而是经过交分解参数分离或半滑处理

6、后的每个地震的编导数矩阵。则通过方法，可将盼解为：奇异值分解（）算法的奇异值形式，（，：，），是甜奇异值，和分别为甜左、有奇异向量组成的正交阵，在计算上通常记，。算法（）最著名并广泛使用的奇异值分解（）算法（）则，（）是（；（）算法（）；该算法由两阶段构成。第一阶段构造变换“（，进而得到，仇，式中）和（一，），使得并行软件包及其在地震层析成像中的应用化。该方法同最重要的不同是其对一个比×××（）更小的矩阵矗（膏）进行转换操作，这样当时，将月二对角化所进行的处理；方法中的右转换要少得多。接下来的问题是如何将矗二对角化。一个显然的方法是将为输入矩阵，采用基于的（；方法

7、。事实上，假如矗的是。少，可以很容易地证明爿昕对应（分解为：×是一个上二对角矩阵。特别地，。使得第角元素以下所有元素为零，“使得第？行中第一个元索左起的所有元索为零。于以上所引入的所有变换是正交的，这使得同具有相同的奇异值、这样，假如，【，是，的分解，那么有（）爿肌将朋的乘积记作，己作这样为了得到，我们必须得到玎朋自乘积、但是，如果不需要显式地得到，则不必对左转换进行任何的累计计算，这样对于的情况，可以极大地减少计算量。将月二对角化的另外一种方法是通过（；变换，这可在一定程度上利用矗中已有的零元索，然而这种方法和上面提到的二对角化方法执行了同样数量的操作（矗）。当需要累计计算左或有转

8、换时，这种方法将需要更多的计算量（岔），主要原凶是该方法在消除每个零元时需要两个旋转，并且这些旋转必须被进行累计乘。但由于转换变换包含了更少的加运算和数组访问，凶此在现代计算机上可能更有竞争力。另外，使用快速（变换也可能在一定程度上改善计算效率。但凶为该方法住并行化实现时不能得到同相比更好的性能，凶此关于曰矩阵的二对角化方法不做进一步讨论。假若我们需要显式地得到，可以列方法的第二阶段进行更进一步的改进。为了显式地得到，我们需要将（，）作用到积阵，从而得到矿在，这些计算可能带来大量的作。凶为单个变换仅影响到白两列，而计算寥“通常需要半均个变换的乘积、凶此，在这样的情况下，首先累计计算所有的“得到

9、其中，而，”，”）第二阶段通过方法将矩阵广对角化。其过程是通过一系迭代使得，（÷÷这里”）（），而和。是由变换乘积得到的交阵，这样（式中的和实质上分别是“，）和“，）的乘积。通常，。的平均迭代次数小于，也就是说，“通常是列角矩阵的一个好的近似形式、改进后的修正算法（当时可以采用两种方法对（；进行改进。在方程（）中，转换和。必须被应用到规模为（）（，一）的子矩阵。因为子矩阵的大部分元索最终将变为零，在情况下，如果能设计一个方法使得“不被施加到子矩阵的次对角部分，则可以节省大量的计算量。这实质上可以通过首先使用变换将矩阵转化为上三角形式：爿÷一个×，矩阵疲为

10、有效，然后在已经收敛到后形成矩阵积在的情况，结合上述两种情况，我们可得到算法给的修算法：（），、来达到上述目的，这里膨黾一个？的上三角矩阵，算法：骆，韵修正算法（）通过一系列列爿的左变换铲厶是一个正交矩阵。然后再将阶的礅日阵二对角并行软件包及其在地震层析成像中的应用将转化为上三角矩阵线性代数子程序集合组成。在原版本基础（瓦一这里矗是？阶上三角矩阵（）通过（；算法得到矗的奇异值分解上，我们对软件包的功能、数据分布方式、并行算法的性能和设计方法等方面进行了扩展和改善。使得该软件成为可以并行求解包括稠密特征问题、稀疏特征问题和奇异值分解等在内的大规模特征问题并行软件包、在中，我们为该软件包对应的并行

11、求解器分别提供了双精度实即任。（）形成：日，数复数和单精度实数复数等不同的代码通过可以并行计算标准厄密（对称）特征问题：和广义厄密（列称）特征问题：凡。最终得到的分解。下面表给卅了两个算法在下面种情况下的渐进计算量列比：需要式得到左奇异向量和有奇异向量其中和是厄密（对称）矩阵，而是正定的，是列应于特征向量的特征值。也能处理与线性代数系统相关的计算，包括的并行分解、的并行分解、正定矩阵的并行分解、并行求解具有多个右端详的上下三角矩阵线性方程组等代码是使用消息传递接口的和言编写，部分代码采用了混合编程方式。可凋用的子程序是和代码程序，软件包的运行要求和库。、：仅需显式得到左奇异向量；：仅需显式得到

12、右奇异向量虬表在，三种情况下和总计算量（）的运行环境是支持消息传递的多种类型体系结构的并行计算机，包括共享存储多处理机、分布存储并行计算机和当前流行的集群系统为了探索在集群环境下的多级并行，我们为标准刘称矩阵特征问题提供了基于自主开发的软件包，我们开发完成了一个并行求解器。通过该并行求解器可以更加高效的并行求解存在于地震层析成像的分解，基于这一并行求解器，结合地质和地球物理研究所已有的计算代码，通过有效整合、代码完善和优化，将形成通用的地球物理反演并行软件包，将建立一个高效的“三维地震层析成像及天然地震数据处理”并行计算半台。地震层析成像的并行求解和应用框架混合并行试验版本，该版本可高效运行在

13、集群环境下。该软件包已用于对光电材料（中山大学光电材料与技术国家重点实验室）和电子轨道结构中上万个未知量的全部本征值问题、百万规模“有机共轭体系分子的本征值问题”（中科院化学所）和规模“地震层析成像”中的分解（中科院地质和地球物理研究所）等多个大规模问题的求解中。概述软件包结构软件包为一个层次结构组织方式。（）是我们自主开发的特征问题并行求解器，主要用于高性能求解实对称或厄密矩阵特征问题及其相关问题。由求解矩阵特征系统的在中包括三类特征问题并行求解器：），称稠密矩阵特征问题并行求并行软件包及其在地震层析成像中的应用解器；程序求解相应的特征问题。在这些接口子程序中分别包含着并行求解特征问题所需要

14、的各种矩阵分解、转换和线性代数子程序。另外在软件包中提供了与通信器和矩阵分解有关的子程序。最低层是支持并行求解的消息传递库和数学库。其层次结构如图。）称稀疏矩阵特征问题并行求解器；）的并行求解器、我们分别为上述三类特征求解器提供了相应的代码。应用程序通过调用相应的接口子厄密对称稠密矩阵特符问题求解器数据在处理器上分布方式并行算法的效率较大地依赖数据在处理器问映射（分布）和处理器间通信所使用的算法有关。的负载半衡和通信特性。的移动代价。：这个代价又与数据在处理器上的初始并且数据在处理器的分布方式极大的影响着算法当前在稠密线性代数中被普遍使用的数据分一、，、一一一一一。，一一一一一一矗、厄密对称稀

15、疏矩阵特衍问题求解器求解器分解并行矩阵分解，转换及与特征求解器自关的线陛代数子程序管卵通占、错吴核鲞和些辅助子程序（）图吏观的层次结构矩阵元索，）的分布位置来确定。在这里假设矩阵的元素（，）被映射为处理器（，）的本地矩阵的元索（，）。：这时全局元索（，）被存储为处理器僦卜蚓此本地赫“”一“：啬卜一一一一，元索。冈给了将一个存储块大小为的阶矩阵分布到处理器网格的示意图。其中“，）中的（表示处理器在网格中的逻辑坐标，其列应布是块一循环数据分布。该数掂分布方式支持可扩展的高性能并行稠密线性代数编码，并被、和其它并行稠密线性代数库选为首要的数据分布方式、一个一块循环数据分布由处理器网格（）、局部块大小

16、（×）和一的处理器序号为：脚亡声矗这样每个处理器序号与它在网格中的逻辑坐标间存在着一对一的关系。俺“霉二【，（【。啦呸。如啦，硅如如如氐龟地呸呸氆呜呜髦静，哓篡囊（劾呜。吗。“毡静蝮，吨吗，气呜。吗吒龌喁呜呸呜呜啦。吗。龌龌吃吃，。喁，喁，强喁，呜咏。燃，吗。蝌瓯【喁唾畦÷一。强瓯嘞强罔将一个按划分的阶矩阵（左）影射到×处理器网格卜（一）的示图并行软件包及其在地震层析成像中的应用并行求解模块和应用框架在软件包中，提供了并行求解通过一系列左变换并行求解×矩阵爿的形式，其中交矩阵是向量得到。通过一系列左、存交替的变换并行完成矗的二列角化分解。下面给了基于该求

17、解器的“地球层析成像并行求解半台”的框架模型结构示意图。器，在该求解器中提供了和两个并行模块分别用于矩阵的上三角矩阵转化和矩阵的二对角化过程。当时，冈“地球层析成像并行求解平台”的框架模犁结构示意罔图为环渤海成像远震数据（×性能测试和×）在深腾超级计算机上的并行顺序交三角化和并行求解器的效率测试刘比和优化前后的加速比（可扩展性能）列比。从冈中可以看到基于我们开发并行求解器，并行计算结果得到了更好的性能和加速比，环渤海成像远震层析成像测试中国大陆及周边地区层析成像数据本测试采用中国大陆周边地区地震层析成像实际数据，矩阵方程组大小为行，列。凶数掂量巨大，在深腾超级计算机上至少需个节点（个）才能完成运算，冈所示为个的情况下程序的运行时间。区内地震的分解运算和远震的平滑运算约占总计算量的，罔列和并行效率测试（左）和优化前后的加速比（，）图×规模问题的并行效率测试时间并行软件包及其在地震层析成像中的应用地球层析成像远程可视化冈是基于地震层析成像网格平台得到的（）：可视化求解结果，数据来源于中围大陆周边地区的反演数掂。，（）：（）刘福田，李强，吴华等，用于速度罔像重建的层析成像法地球物理学报：罔中国大陆周边地区地球层析