大地电磁一维自动反演算法及OPENMP并行效率分析

本科生毕业论文论文题目 ：大地电磁一维自动反演算法及 OPENMP 并行效率分析姓名 ： 向馨梅学号 ： 09023316班级 ： 090233 班年级 ： 2009 级专业 ： 勘查技术与工程学院 ： 核工程与地球物理学院指导教师 ： 陈辉、邓居智完成时间 ： 2013 年 5 月 30 日作者声明本人以信誉郑重声明：所呈交的学位毕业设计（论文） ，是本人在指导教师指导下由本人独立撰写完成的，没有剽窃、抄袭、造假等违反道德、学术规范和其他侵权行为。文中引用他人的文献、数据、图件、资料均已明确标注出，不包含他人成果及为获得东华理工大学或其他教育机构的学位或证书而使用过的材料。对本设计（论文）的研究做出重要贡献的个人和集体，均已在文中以明确方式标明。本毕业设计（论文）引起的法律结果完全由本人承担。本毕业设计（论文）成果归东华理工大学所有。特此声明。毕业设计（论文）作者（签字）：签字日期： 年 月 日本人声明：该学位论文是本人指导学生完成的研究成果，已经审阅过论文的全部内容，并能够保证题目、关键词、摘要部分中英文内容的一致性和准确性。学位论文指导教师签名：签字日期： 年 月 日大地电磁一维自动反演算法及 OPENMP 并行效率分析向馨梅1-D Automatic Inversion of Magnetotelluric And The Analysis Of Parallelization Efficiency based on OPENMP Xiang,Xinmei2013 年 5 月 30 日东华理工大学 2013 届本科生毕业论文 摘要摘 要大地电磁测深反演算法理论日渐丰富和成熟；一维反演理论作为高维反演方法研究的基础，对于其深入研究是很有必要的。但在一维反演中往往需要人工给定地下地质体层数，要求数据解释人员具有很高的专业水平。同时大地电磁多维反演非常耗时，多核计算机的并行计算过能够提高计算速度。因此本文介绍了一维自动反演理论，并以此为例对基于 OPENMP 的多核计算机并行化的可行性效率进行分析，为多维反演提供参考和借鉴。本文首先对 OPENMP 并行库的特点和优势进行描述，在此基础上设计一个简单求解圆周率的程序，通过对比并行计算和串行计算表明：基于 OPENMP 的并行库采用共享内存模式，非常适合当前流行多核计算机计算有效提高计算速度节省计算时间。然后介绍大地电磁一维自动反演基本算法，它首先以 Bostick 反结果作为初始模型，然后自动构建 2,3,4层模型，通过对比最小的拟合差以及有关统计的 F-test 来确定最终模型层数和厚度。本文通过 G,D,A,H 型模型以及实测数据计算表明：该自动反演算法能够很好的匹配地下模型的层状参数，具有较好的实用性。关键词：一维自动反演； 并行计算； OPENMP东华理工大学 2013 届本科生毕业论文 ABSTRACTABSTRACTThe Theory of MT inversion algorithm is increasingly mature and rich;As the base of the multi-dimensional inversion algorithm,the further research of 1D inversion algorithm is necessary.While setting the number of the model layer is usually required in advance,which asks for the analysts with much experience and technical knowledge.As MT multi-dimensional inversion needs much time,and the parrallelization algorithm on the multi-core PC could improve the caculation efficiency.The paper introduced the Theory Of 1D Automatic Inversion,and did some analysis of the parrallelization feasibility of auto inversion in preparation for Multi-inversion.In the part of the test,with an example of caculating Pi,by comparing the serial algorithm with the parrallelization to illustrate the efficiency and extensibility of parrallelization based on OPENMP using shared memory method,which is good suitable with multi-core PC to decease the time of caculation.Then the paper introduced the theory of 1D Automatic Inversion ,which successively construced 2,3,4layer model from an initial model based on Bostick Inversion Theory,determing the final modle with the minimun misfit and a statistial F-test.With G,D,A,H models to explain Automatic Inversion how to work ,it turned out to be perfectly matched with the theory modles.And to some degree,a further test with a set of pratical data illustrated the applicability of the algorithm.Key words：1D Intersion； Conditions Of Parallelization； Openmp东华理工大学 2013 届本科生毕业论文 目录 目 录一、绪 论 .11（一）MT 一维反演理论介绍 .11（二）并行计算在勘探地球物理中的应用 .11（三）本文的主要内容及安排 .22二、并行计算 .33（一）介绍 .33（二）并行计算的理论基础 .33（三）MPI 与 OPENMP.44（四）openmp 并行效率分析 66三、一维 MT 自动反演算法 88（一）介绍 .88（二）一维 MT 自动反演算法 88（三）实例分析 .1313（四）一维自动 MT 反演并行化可行性分析 1818四、总 结 .1919致 谢 .2020参考文献 .2121东华理工大学 2013 届本科生毕业论文 绪论 1一、绪 论近年来，在地球物理勘探领域，并行技术已经成为炙手可热的研究方向。在地震、电磁法、重磁等方面对于并行计算的应用，地球物理研究者们，已经作出了初步的探讨和研究并取得了一系列的成果。电法勘探，作为地球物理勘探的一个重要分支，以其低成本、简便易行的特点，它的发展也受到人们的关注。对于电磁法日渐深入的研究，由于计算量以及计算难度的增大，其算法也在向着并行化的方向发展。本文将探讨，通过接收天然场源平面电磁波信号，研究地下构造的大地电磁法测深（MT）的一维自动反演算法，并对其并行可行性进行分析。（一）MT 一维反演理论介绍MT 一维反演理论发展较早，在早期的地球物理勘探研究中占有很重要的地位。一方面，可以通过一维反演得到的模型来对研究区域的整体的地电分布情况作定性的认识；另一方面，可以为高维反演的研究提供一定的理论基础。 【19】（二）并行计算在勘探地球物理中的应用地球物理的数据处理，往往需要面对海量的数据以及复杂的数值计算。计算机硬件的发展速度在一定时间范围内是有限的。如何在现行条件下以较低的成本实现计算精度和效率的提高，并行计算为数据处理提供了一个可行的研究方向。目前，在地球物理领域使用较多的是基于机群(COW)的 MPI 和基于共享内存的在多核（Muli-core ）PC 机上实 OPENMP。 【1】当前，并行计算在地震数据处理中应用较为广泛 【8】 。因为地震本身就是多道记录，后续的部分处理，各道之间没有逻辑关联，就可以实现并行化给多个线程，最后，再将结果反馈给主线程 【2】 。应用于电磁法的并行 【5】 ，主要针对正演部分的频率并行（Frequency loop parallelization）和反演部分“最小二乘法”拟合时，求解雅可比矩阵的参数并行（Parameter update parallelization） 。在重磁方面 【17】 ，有人在正演部分采用了 GPU/CPU 协同并行技术。由 CUP负责主程序串行主体，GPU 负责线程化任务。其实这有点类似 MPI 和OPENMP 的混合编程模式。通过试算表明，并行计算确实能够节省程序运行所需要的时间，具有可推广东华理工大学 2013 届本科生毕业论文 绪论 2性。现有的并行化，大多是在串行代码的基础上，通过嵌入一些并行语句，将整个程序的一部分并行化 【9】 。这种并行方式，由于受到原始代码串行结构的限制，并行程度并不高；另一方面，其效率能否很大的提高，很大程度上取决于并行计算部分在整个程序中运行时间所占的比例大小。想要实现效率提高最大化，需要从根本的并行算法研究入手。（三）本文的主要内容及安排本文主要探讨了一维自动反演算法及其并行可行性。通过分析研究了程序并行化所需要满足的条件。全文大体分为两个部分，第一部分结合一个小的算例（关于圆周率值的计算） ，来说明 openmp 在并行运算方面的方便和高效性；第二部分介绍了一维自动反演理论，通过四个模型（G 、D 、H、A）和一组实测数据的试算和分析来对该算法进行验证，并对其并行可行性作了一定的分析。东华理工大学 2013 届本科生毕业论文 并行计算 3二、并行计算（一）介绍并行计算（Parallel Computing） ，就是将某些问题（需要多次循环计算，但各计算过程又相对独立，没有逻辑关系的）开辟多个线程，同时运行计算。并行计算的主要目的是快速解决大型且复杂的计算问题。像其他计算问题一样，并行计算依赖于计算机硬件（并行机）和相应的并行计算理论。现在主要的并行机有以下几类：并行向量处理机（PVP） 、对称多处理机（SMP ） 、大规模并行处理机（MPP ） 、工作站机群（COW） 、分布式共享存储处理机（DSM） 。而现在比较通用的并行编程模型有以下几种：共享内存、线程、消息传递、数据并行、混合模型等。（二）并行计算的理论基础1.并行计算机系统(1)共享存储对称多处理机(SMP)(2)分布式存储并行处理机：大规模并行系统(MPP) 机群系统(Cluster)2.并行计算环境并行机系统，还要在相应的并行环境下才可以得以运行。并行程序开发环境主要包括：并行语言、程序优化重构工具、图形可视化工具、编译工具等。PVM、 MPI 是现在较为主流的并行程序开发环境。PVM(Parrallel Virtual Machine)，是一种通过计算机网络连接得以实现的并行系统。PVM 由：环境控制管理部分、和传递消息的子程序库组成。MPI(Massage Passing Interface)，及消息传递界面标准。其特点主要是：任务之间通过调用接收和发送指令来进行数据通信，数据的传输需要“接收”和“发送”的配合。MPI 既可以在分布式系统（比如说，COW）中运行，也可以在共享内存的的系统中运行。OPENMP，是应用编程接口 API，在共享存储体系结构上的一个编程模型。其主要包括编译制导(Compiler Directive) 、运行库例程(Runtime Library)和环境变量(Environment Variables) 。东华理工大学 2013 届本科生毕业论文 并行计算 43.并行算法程序并行化的实现，一方面可以通过将串行算法改写成并行；另一方面，是构造数学模型，形成新的算法。前者的实现较为为容易，但由于受到串行结构的限制，难以取得较高的并行度。而后者的实现目前进展不大，还有待进一步的深入研究。现较为成熟的一些理论和技术 【2】 主要有：NC 理沦、平衡树法(Balanced Tree)、分治策略(Divideand Conque)、切分原理(Partioning)、及流水线技术(Pipelianing)等。4.并行算法的性能评价一个并行算法的优劣主要由以下几个指标来决定：成本（ ） 、)(*)(nptc加速比（ ） 、效率（ ） ，以及可扩展性。其中： 并行运行psTS/PSE/时间， 处理器个数； 串行运行时间， 并行运行时间；S 加速比，P 处)(ns pT理器个数。（三）MPI 与 OPENMP1.MPI数据并行和消息传递是现在最主要的两种并行编程模型。数据并行 【3】 即将相同的操作作用于不同的数据。但其应用范围相对局限，只针对数据并行类问题。消息传递模型，指的是对于各个并行部分通过“消息传递”来交流、协调和控制。消息传递既适用于分布式内存，也适用于共享内存的并行机。消息传递的灵活性和可移植性，让它应用广泛。MPI 就是基于消息传递。MPI 分为两种模式：主从模式和对等模式。MPI 中最重要的两个语句就是，SEND. 和 RECEIVE.两个语句是配对使用的。很重要的就是数据的标识(TAG) 。好比一个人抛球，另外一个人接球。2.OPENMPOPENMP 应用编程接口 API，是在共享存储体系结构上的一个编程模型。包含编译制导(Compiler Directive)、运行库例程(Runtime Library)和环境变量(Environment Variables) 。它支持增量并行化(Incremental Parallelization) 。OPENMP 简单实用，可移植性强，在单机上操作即可，相较于 MPI 更加方便东华理工大学 2013 届本科生毕业论文 并行计算 5易学，具有很强的推广性。OPENMP 是一种基于线程的并行编程模型(Programming Model) ，它的执行模式是 Fork-Jion（图 1-1） 。图 2-1： Fork-Jion 模式下面将介绍一些 OPENMP 常用的语句，并对其并行模式 【18】 进行简要说明。在 fortran 里面，只要设置好环境变量，就可以直接使用下面语句了。在并行没有开启时，代码将按照串行的运行，并行语句将被视为注解。例如：!$OMP PARALLEL!$OMP END PARALLEL主要的语句有以下一些：parallel for, parallel sections;critical, single;barrier, atomic;master, ordered;threadprivate几个常用的库函数：omp_get_num_procs, omp_get_num_threads, omp_get_thread_num, omp_set_num_threads, 运用 OPENMP 编程，首先要选择合适的并行区域。任务可被划分成若干小块（BLOCK） ，可以等量分配，也可以自己根据任务量的大小实行不等量的分配。其次，在运行过程中有一个主线程（MASTER ） ，它对其他线程（SLAVE）起到控制的作用。对于，各个线程中的变量是否共享（ SHARE）也是一个很关键的问题；在定义并行区域时，申明变量是否为哪个线程私有（PRIVATE） ，是很有必要的。这样，各个线程才可以独立运行，互不影响；否则可能出现数据竞争的现象。如果各个线程的任务量不同，通过LOCK，BARRIR 这样一些语句，就可以对一条线程的运行进行控制；使得各FORK JOIN FORK JOIN主 线 程 并 行 域 并 行 域东华理工大学 2013 届本科生毕业论文 并行计算 6个线程之间协调运行，达到某种平衡。程序的运行将更为高效。（四）openmp 并行效率分析 这里，将通过一个关于计算圆周率 Pi 的算例,来对 openmp 的并行效率进行验证。其中，涉及到的主要参数有：nthreads （线程数） ，N（迭代次数） ，error（误差） 。用于该计算的计算机硬件参数：CUP 主频：2GHZ ；内存：2GB；Intel ie5(双核) 。并行算法设计如下：omp_get_num_procsomp_get_thread_numcall omp_set_num_threads ( )!$OMP parallel private ( )!$OMP end parallelEnd下面，将从算法的正确性和有效性两个方面来验证并行化后的算法是否更优。1.有效性该算法的有效性 【1116】 主要由以下几个指标来说明：加速比（ ） 、psTS/效率（ ） 。PSE/下面的表 1-1，表 1-2，表 1-3，分别呈现了在不同线程数下，算法运行时间，加速比，以及效率。表 2-1 表明，在较大迭代次数的情况下，程序运行时间的提高是很显著的；实际上，程序在较小的迭代次数情况下，并行后运算的时间可能出现，比在相同迭代次数的串行运行的时间更长的情况。这也说明了并行程序设计需要考虑的通信开销的问题。表 2-2 表明，随着线程数的增加，程序运行的速度确实有一定程度的提高。但这并不是说明只要增加进程数，就可以提高程序的运行效率；效率的提高是在一定限度之下的，需要根据具体的并行区域算法的特点，来予以分析线程数的开辟的饱和点在哪里。东华理工大学 2013 届本科生毕业论文 并行计算 7表 2-3 表明，每个线程的使用程度是不同的。也是需要根据实际情况来使得开辟的线程的任务完成度相当，达到负载平衡的要求。另外，当算法更为复杂时，每个线程的任务分配可能是不均等的，需要进一步更为细致的划分（细粒度） 。现在比较流行的 mpi 和 openmp，都只是一种“人为”并行的平台，进一步的研究，应该是二者的结合（机群和多核） ，以及更为自动的“智能”的并行方式。表 2-1：不同线程数下的运行时间Time (s)nthreads 1 2 4 8Iter=1E6 0.031 0.015 0.016 0.011Iter=1E7 0.203 0.125 0.183 0.047表 2-2：不同线程数下的加速比S(加速比)nthreads 2 4 8Iter=1E6 2.1 1.9375 2.819Iter=1E7 1.624 1.11 4.319表 2-3：不同线程数下的效率值E(效率)nthreads 2 4 8Iter=1E6 1.05 0.484 0.352Iter=1E7 0.812 0.2775 0.54这里只是通过小的算例来说明“并行” （openmp）在提高程序运行效率方面的成效性。对于地球物理数据处理，并行计算的意义就显而易见了。不过，上述的并行方法，适用于那些数据量较大，循环变量比较单一独立的算法；满足上述条件，使用 mpi 或者 openmp 将原有的串行代码改为并行才具有可行性。东华理工大学 2013 届本科生毕业论文 一维 MT 自动反演算法 8三、一维 MT 自动反演算法（一）介绍大地电磁测深（MT） ，是通过研究天然场源来研究地球内部结构的一种地球物理方法。当交变电磁场地下介质中向下传播时，地面上将接收关于电磁波频率的频率响应信号。不同的频率信号对应着不同深度的信息 【4】 。根据法国学者 Cagniard 的经典理论，当假设天然电磁场以平面波形式垂直入射均匀各向同性的层状介质时，可得到以下一个的解。（1）fZ52其中， 为卡尼亚是电阻率； 为频率；Z 为波阻抗，其定义为：xyxHEZ早期的一维反演理论有 Bostick 反演、高斯-牛顿法、梯度法、广义逆反演、马夸特法、连续介质反演等 【19】 。（二）一维 MT 自动反演算法1.Bostick 反演理论自动反演中将用 Bostick 反演得到的模型，作为初始模型。下面将介绍与Bostick 反演有关的理论基础。从水平层状介质理论公式（2）出发，可推导当 或 0时，可得到（ 3） （4） 【4】 。（2） nnnnhkhknnn eZekZ2121其中， ； 为电导率；h 为电性层厚度； nZ为第 n 层波阻抗。nik/当下半空间介质东华理工大学 2013 届本科生毕业论文 一维 MT 自动反演算法 9， （3） 21S0， 2D （4） 其中，S 下半空间以上多层介质的总纵向电导；D下半空间以上多层介质的总厚度。再根据绘图中观察发现的规律，进一步推导，可以得到 Bostick 反演理论的公式（Niblett，Sayn-Wittgenstein 1960）：（5）)(1)(nnnTmK02)(nnTZ其中， ， k为电导率； n为深度； ，nZk Tddm101log/)(log)(可以由相位估算得到（Weidelt 1972）： （6）/)(4)(T根据 Schmuckers （1970 ）的推导 nk， Z也可以由 计算得到：2cos)(1nnTK 02)(sinn（5） ， （6）中的 即卡尼亚视电阻率，在实际反演中，就是地面实测数据或者正演得到的电阻率值。 都是关于 Tn 的函数，而 Tn 对应着不同深)(,nnZk度，根据各测深点不同深度的频点值，即可反演得到该深度下的电阻率值。2.一维 MT 自动反演算法 【6】 由于反演问题一般都是非线性的，在构建初始模型的时候，将会有很多的限制条件。找到全局最优解是非常困难的。（1）首先，由 , 或 （ 为卡尼亚视电阻率， 为相位；)(nT)(nTZ为波阻抗） ，基于 Bostick 反演理论得到一个初始模型)(nTZ（ 为电阻率， 为深度） 。 （n=1,2,M ）);(,nKn（2）根据（1）中的反演，可近似认为每一层的分界在 的位置，即：nZ东华理工大学 2013 届本科生毕业论文 一维 MT 自动反演算法 10， （n=1,2,M-1 ） （7）011nnZ其中， ， 由方程（5）给出；关于第 M 层，及),.2(1MZn K基底均匀半空间区域（ ）的 由 频点处的的值给定。这样，就可以MT计算 M 层分界面以上各层的厚度以及纵向电导率的值了（ ）：nd,， （n=1,2,.,M-1） （8）nnd1ndK（3）下面的步骤将开始自动反演。首先，由两层模型开始。由 这两个参1,h数来定义这个两层模型：， （9）1Mndh11Mnh此时，将由最小拟合差（misfit） （10）来计算得出上述的参数（ ） 。)2(1)2(,h（10）)()(12 2(2)(21 knnINnknnRkN ITTks ， NnRnnRT12)(/)/(NnInnI12)(/)/(（11）)()(minax)()(InR TZi（12）其中， 为拟合差的方差，它的自由度为 2N-2k+1；21kNs分别为实际计算出的波阻抗 的实部和虚部；)(,nTIR 0)()(ZTiIRTZ分别为 k 层模型下通过最小拟合差计算得到的 的实部和虚部；)(,knI )(kn分别为 的实部和虚部；IRT)(n是决定高斯分布(Guassian distribution)误差棒（error bar）的权重因子。下一步，将构建一个三层模型（ ） 。21321,;,h东华理工大学 2013 届本科生毕业论文 一维 MT 自动反演算法 11， 21Mndh12)(hdhMn， ， 211n2Mn)2(3其中， ；/)(2通过最小拟合差 的取得，可以计算得到 。52Ns )3(2)3(12)3(1,h当计算 K 层模型时，定义它的参数，计算公式如下：， （j=1,2,.,k-2） （13）kjjMndhkjjMnjh,211knk（14）,1Mnkkh)1(k其中， ；最小拟合差 取决于/)(,)(jnj 12kNs的值。)1()(,.2, kkjjh（4）在“自动反演”进入 k+1 层模型的计算之前，有一个关于统计的判断，将会决定程序是否继续往下进行。K 层模型计算得到的值，将被带入下面的式子（15） 。 值的计算，将会对数据的随机性作出检验。我们只需定义计算最小拟s合差的权重因子 。1.0（15） Nn InkRnkk TT1 2)(2)(2 （5）当 ，且模型经过了（7）的测试，就可以将 k-1 层和 k 层计算得到的3，通过下面（16） （F-test ）的计算，来判断增加了一层的模型，是否比先前2的模型有更好的拟合度。（16）)12/(/()12,(213kNkNFk当然，k 层模型一层一层的构建不会无限制地一直进行下去。当 k 达到 10层时，将会结束程序。此时，程序将比较之前计算得到的，在某层模型下的拟合差；取其中的最小值所对应的 k 层模型参数。 （4） 、 （5）两项检测将会跳过。东华理工大学 2013 届本科生毕业论文 一维 MT 自动反演算法 12上述搜索过程的终止，可能是因为新增加的一层没有使得模型与实际数据的拟合度增加，另外一个原因，就是计算的得到的模型参数的分布不满足误差棒（error bar ）的限制范围。（6）如果，最佳模型是由 L（ ）这个限制条件获得，那么，将进行下面10一个检测公式（17） 。该检测的目的是，将相邻层的电阻率之间差异小于给定的某个值（程序中，我们设定 0.2）时，相邻的两层将合并为一层。这样，可k以减少得到的层状模型的层数。（k=1,2,.,L-1） （17） ,max)(1)(Lkkk（7）经过了步骤（7）-（8） ，最终的模型将被选定。这个时候，可以再次通过拟合差的方差计算公式，来得到最终模型的参数值。此时，可以选定更大的迭代次数值和更小的步长值。前文中提到，波阻抗 ，可以通过公式（ 18）和 相互转换；下面将给出，Z,初始数据为 时，上述自动反演过程的公式（19）-（21） ） 。,constTiTZ)(/)(ln)2/1(ln0（18） Nn knnknnkN Tks1 2)(2)(21 )(l/)(2 （19）NnnnNnnnT12)(12)(/)/(,/)/(（20）)/ln(),/l()2/( minax(minax) n（21） 如上文中所述， “自动反演”在读入数据后，将计算得到一个初始模型，该模型基于 Bostick 反演理论。在这个模型反演得到层参数结论之下，程序计算得到 k 层模型下的最优模型层参数。首先，程序将开始于一个两层模型。根据公式（19） ，程序将不断修改一个东华理工大学 2013 届本科生毕业论文 一维 MT 自动反演算法 13两层的初始模型的层参数，计算得到一个最优的两层模型，此时，程序将返回得到 Epsilon 的值；然后，程序将三层，四层，这样逐层增加下去，计算得到相应层数下的最优的层参数值；在进入 k+1 层的计算前，有一个关于统计的判断（F-test） ，如果不能通过该检测，程序将返回到 k-1 层的模型下；当然，在模型满足上述判断的条件下，程序也不能无限制的在层数上增加下去，程序中设定的最大的层数值为 10，当到达 n=10 时，程序将比较之前各层数下计算得到的 Epsilon 的值，取其最小；程序在上述的两种情况下退出后，将进一步地修改得到模型的层参数，使得模型进一步的优化。StartRead data AutoK layer modelInitial modelFinal modeltesF,2Find misfit 2,3.layerAUTOEnd NNY东华理工大学 2013 届本科生毕业论文 一维 MT 自动反演算法 14图 3-1:一维自动反演流程图（三）实例分析1.理论模型反演下面作为反演得理论模型数据都没有加入噪声。由图 2-22-9 可知，反演结果是比较理想的。其中 R-D 曲线加入了 Bostick 反演得到的模型曲线，可与自动反演得到的模型进行对比。总的来说，无论是 R-D 曲线还是 reZ -logT Im Z-logT 曲线，理论模型和反演模型的响应都比较吻合，说明了一维自动反演在一定模型参数下的适用性。在具体反演过程中，由其他层数的模型与理论模型的拟合差大小可以看出，地球物理反演的多解性；以及在一维反演中，水平介质层数的确定对于反演结果的影响也是很大的。在某种层参数之下，三层的模型也可能由反演得到四、五层的模型。Noise Y with greater number of interations and smaller step size东华理工大学 2013 届本科生毕业论文 一维 MT 自动反演算法 15两层模型：设置 D 型模型，第一层的电阻率为 1000 ohm.m，厚度为 1km；第二层的电阻率为 10 ohm.m。 由成果图可分析，反演结果是比较理想的（如图 3-2，图 3-3） 。反演模型与理论模型的电阻率的相对误差约为 0.1e-3；26 个频点的 reZ的平均相对误差为 0.399e-3 ；ImZ 的平均相对误差为： 0.321e-3。 由反演得到的模型为一个两层模型，其层参数为：第一层的电阻率为998.98 ohm.m，厚度为 1km；第二层的电阻率为 9.99 ohm.m。图 3-2：反演模型与 D 型模型的 R-D 曲线图 3-3：反演模型的 reZ-logT ImZ-logT 曲线设置 G 型模型，第一层的电阻率为 10 ohm.m，厚度为 0.5km；第二层的电阻率为 1000 ohm.m。 由成果图可分析，反演结果是比较理想的（如图 3-4，图 3-5） 。反演模型与理论模型的电阻率的相对误差为 0.101e-3；26 个频点的 reZ的平均相对误差为 0.686e-3 ；ImZ 的平均相对误差为： 0.648e-3。 由反演得到的模型为一个两层模型，其层参数为：第一层的电阻率 9.99 ohm.m，厚度为 1km；第二层的电阻率 998.99 ohm.m。东华理工大学 2013 届本科生毕业论文 一维 MT 自动反演算法 16图 3-4：反演模型与 G 型模型的 R-D 曲线图 3-5：反演模型的 reZ-logT ImZ-logT 曲线三层模型：设置 A 型模型，第一层的电阻率为 10 ohm.m，厚度为 1km；第二层的电阻率为 50 ohm.m，厚度为 0.2km；第三层的电阻率为 1000 ohm.m。由成果图可分析，反演结果是比较理想的（如图 3-6，图 3-7） 。反演模型与理论模型的电阻率的相对误差为 0.1e-3；26 个频点的 reZ的平均相对误差为 0.686e-3 ；ImZ 的平均相对误差为： 0.658e-3。由反演得到的模型为一个三层模型，其层参数为：第一层的电阻率 9.99 ohm.m，厚度 0.2km；第二层的电阻率 49.95 ohm.m，厚度 1km；第三层的电阻率 999.01 ohm.m。图 3-6：反演模型与 A 型模型的 R-D 曲线东华理工大学 2013 届本科生毕业论文 一维 MT 自动反演算法 17图 3-7：反演模型的 reZ-logT ImZ-logT 曲线设置 H 型模型，第一层的电阻率为 50 ohm.m，厚度为 0.5km；第二层的电阻率为 10 ohm.m，第二层的厚度为 0.3km；第三层的电阻率为 1000 ohm.m。 由成果图可分析，反演结果是比较理想的（如图 3-8，图 3-9） 。反演模型与理论模型的电阻率的相对误差约为 0.1e-3；26 个频点的 reZ的平均相对误差为 0.616e-3 ；ImZ 的平均相对误差为： 0.591e-3。由反演得到的模型为一个三层模型，其层参数为，第一层的电阻率 49.95 ohm.m，厚度 0.5km；第二层的电阻率 9.99 ohm.m，厚度 0.3km；第三层的电阻率 999.03 ohm.m。图 3-8：反演模型与 H 型模型的 R-D 曲线图 3-9：反演模型的 reZ-logT ImZ-logT 曲线2.实测数据试算下面将对实测数据进行试算。具体的试算过程中，程序进行到计算得出六层模型时，由于无法通过 f-东华理工大学 2013 届本科生毕业论文 一维 MT 自动反演算法 18test，程序退到五层模型下，此时，给以程序更大的迭代次数和更小的精度要求，进一步对五层模型的层参数进行修改，得到最终模型。得到的五层模型的层参数如下，第一层的电阻率 450.01 ohm.m，厚度 0.1km；第二层的电阻率 50 ohm.m，厚度 0.4km；第三层的电阻率 28 ohm.m，厚度 3km；第四层的电阻率 45.01 ohm.m，厚度 7km；第五层的电阻率10007.17 ohm.m。图 3-8：反演模型与 H 型模型的 R-D 曲线图 3-11：实测数据反演得到模型的 reZ-logT ImZ-logT 曲线从图 3-11 可以看出，反演得到的数据与实测数据大体上是比较吻合的，26 个频点的 reZ的平均相对误差为 0.465e-3 ；ImZ 的平均相对误差为 0.33507e-3。表 3-1：反演模型与理论模型的相对误差模 型G 型 D 型 H 型 A 型电阻率(ohm.m) 0.102e-2 0.101e-2 0.99e-3 0.997e-3深度(km) 0 0 0 0相对误差东华理工大学 2013 届本科生毕业论文 一维 MT 自动反演算法 19表 3-2：阻抗数据的实部与虚部的相对误差模 型G 型 D 型 H 型 A 型 实测数据reZ 0.686e-3 0.399e-3 0.616e-3 0.686e-3 0.465e-3ImZ 0.648e-3 0.321e-3 0.591e-3 0.658e-3 0.33507e-3（四）一维自动 MT 反演并行化可行性分析对于 MT 反演的并行化方法主要有两种 【5】 ，一种是针对调用频率的函数的“频率并行” ，另一种是针对最小二乘法方程组求解中，雅可比迭代的“参数并行” 。相对来说， “频率并行”更为简单。在频率域的电磁法中，各个测深点的数据都是关于不同频率（不同深度）的响应；各频点之间相互独立。对于程序中涉及到的不同频率的计算部分（关于 N 循环） ，从 “并行化”要求并行区域的各条线程（或者各个区块）的计算相对独立、没有逻辑关系的条件上来说，都可以将它们并行化。在上述的一维自动反演中，关于频率的循环并不少，但都比较零碎，各个循环的数据量也不大；开辟并行区间所花去的时间，还不如直接串行计算。这涉及到了决定并行计算是否可行的另外一个关键因素通信开销。上述反演过程中，唯一可能“并行化”的环节就是程序中调用Findmisfit 函数的部分（每次调用该函数的迭代次数为四百多次） 。相对误差东华理工大学 2013 届本科生毕业论文 总结 20四、总 结通过试算以及结果分析我们发现，并行算法对于程序效率的提高确实很有帮助。一个并行计算算法的优劣并不完全在于它对于该具体程序的运行时间缩短了多少（当然，运算时间的长短是一个很重要的评价指标） ，还在于它的可扩展性。对于那种变量单一，且各个运算过程相对独立的循环体（Loop）来说，只要迭代次数足够大（或者对于运算精度的要求很高） ，并行对于其效率的提高时显而易见的。 “并行化”需要考虑以下因素：（一）分配给个线程的运算相对独立，前后没有逻辑关系；（二）并行区域中运算的变量是否共享，各线程之间或者各个区块之间是否需要通信（通信开销） ，以及各部分的运算量是否平衡（负载平衡） 。 在一维反演方面，通过试算与分析发现，一维自动反演（1-D Aotumatic invertion）不失为一种比较高效的反演算法。只要给定数据，程序就可以自动地进行下去，直到最终给出较少层数的层状模型。不过，真正意义上的“自动”很难到达。单单依靠计算机对于给定条件语句的判断，有的时候得到的结果可能会偏差很大；毕竟，实际的情况是复杂多变的，真实的状况由多个参量共同决定，程序无法“面面俱到”地考虑所有可能的参量；其只能在一定范围内，选择最为主要的参量和判断准则。由于时间仓促，文章未能完全实现一维自动反演的并行化，在今后可作进一步的探讨和研究。东华理工大学 2013 届本科生毕业论文 致谢 21致 谢首先，感谢陈辉老师在论文过程中给予的指导和鼓励；其次，感谢程励辉师兄和刘俊峰师兄在论文过程中的帮助；最后，还要感谢身边的同学和朋友的支持，和你们的探讨有时会使我茅塞顿开。东华理工大学 2013 届本科生毕业论文 参考文献 22参考文献1 黄易等，并行计算技术及其在勘探地球物理学中的现状与展望J，武汉：中国地质大学地球物理与空间信息学院，2010：642-6492 刘羽等，地球物理数据处理与并行计算J，武汉：中国地质大学地球物理与空间信息学院，2004：412-4173 都志辉，高性能计算并行编程技术-MPI 并行程序设计，北京4 陈乐寿，大地电磁测深资料处理与解释，北京5 Virginie Maris et al. Parallelizing a 3D finite difference MT inversion algorithm on a multicore PC using OpenMJ，Department of Geology and Geophysics, University of Utah，2010：1384-13876 J.T.Weaver，A.K.Agarwal,Automatic. 1-D inversion of magnetotelluric data by the method of medelling,Department of Phisics and Astronomy and School of Earth and Ocean Science,University of Victoria，Victoria，BC，Canada，1992；115-1237 LaustB.Pedersen，Mehran Gharibi. Automatic 1-D invertion of magnetotelluric data:Finding the simplest possible model that fits the data，Uppsala University，Dept.of Geophysics，1999：773-7818 杨进，地震数据处理中的并行计算技术研究，C，湖南：湖南大学，20089 赵连功，刘洪，地震勘探数据资料处理软件集成化研究现状和发展趋势，北京：中国科学院地质地球物理研究所J，2003：202-20810 谢桂生，刘 洪，李幼铭，地震波正演模拟的并行计算初探J，北京：中国科学院地质与地球物理研究所J，2003：243-24611 王宏琳，高绘生，地震并行处理模式与应用框架J，涿州：石油地球物理勘探局，2001：201-20812 张伟，陈晓非，并行计算在计算理论地震图中的应用应用MPIOpenMP 在 SMPcluster 机群采用三维有限差分计算地震图J ，北京：北京大学地球与空间科学学院，2002:46913 李焱，胡祥云等，基于 MPI 的一维大地电磁并行计算研究J，武汉：中国地质大学地球物理与空间信息学院，2010：1612-161614 胡祥云，李焱等，基于 MPI 的二维大地电磁正演的并行计算J，武汉：中国地质大学地球物理与空间信息学院，2010：392-40115 刘鹏茂等，基于 MPI 大地电磁二维正则化反演并行算法研究J，长沙：中南大学地球科学与信息物理学院：241东华理工大学 2013 届本科生毕业论文 参考文献 2316 柳建新，刘鹏茂等，基于 MPI 瞬变电磁测深一维反演并行算法探究J，长沙：中南大学地球科学与信息物理学院，2011：491-49517 陈召曦等，重磁数据三维物性反演方法进展J，北京：地质过程与矿产资源国家重点实验室，中国地质大学，2012：503-51118 张军，基于 OpenMp 的并行计算技术及应用研究 硕士论文，广西：挂林理工大学，201119 冯思臣，一维大地电磁测深反演算法比较研究 硕士论文，四川：成都理工大学，2007