基于新离散公式和OpenMP优化的有限差分声波数值模拟

第３９卷　第４期２０２０年１１月 世　界　地　质ＧＬＯＢＡＬＧＥＯＬＯＧＹ
Ｖｏｌ ３９　Ｎｏ ４
Ｎｏｖ ２０２０ 文章编号：１００４
５５８９（２０２０）０４
０８９６
０９
基于新离散公式和ＯｐｅｎＭＰ优化的有限差分声波数值模拟
郑如秋，王波涛，冯永照，余卫江
中海油田服务股份有限公司特普公司，广东湛江５４２０５１
摘要：有限差分法算法是声波数值模拟算法中最广泛的数值计算方法。

该方法具有计算速度快、占用
内存相对较小、易于编程实现及模拟精度高等优点。

规则网格下的有限差分算法相比于近年来提出的新的有限差分算法计算效率最高，消耗内存最少。

然而，引入ＰＭＬ边界，导致常规网格的有限差分离散公式复杂，计算过程需要对计算区域和边界区域进行判断，导致计算效率低。

针对这个问题，笔者新推导的离散公式，形式简单，整个模拟区域计算代码一致，可以很好地解决这个问题。

ＯｐｅｎＭＰ并行算法，语言简洁和可移植性高，通过结合ＯｐｅｎＭＰ并行算法，对模拟算法进一步优化，可以较大提高数值模拟的计算效率。

关键词：有限差分；规则网格；离散公式；Ｏ
ｐｅｎＭＰ并行算法；计算效率中图分类号：Ｐ６３１ 文献标识码：Ａ
ｄｏｉ：１０ ３９６９／ｊ ｉｓｓｎ １００４ ５５８９ ２０２０ ０４ ０１５收稿日期：２０２０ ０３ ２７；改回日期：２０２０ ０８ ２８
基金项目：国家科技重大专项项目（２０１６ＺＸ０５０２６ ００２ ００４）
Ｆｉｎｉｔｅｄｉｆｆｅｒｅｎｃｅａｃｏｕｓｔｉｃｎｕｍｅｒｉｃａｌｓｉｍｕｌａｔｉｏｎｂａｓｅｄ
ｏｎｎｅｗｄｉｓｃｒｅｔｅｆｏｒｍｕｌａａｎｄＯｐｅｎＭＰｏｐｔｉｍｉｚａｔｉｏｎ
ＺＨＥＮＧＲｕ ｑｉｕ，ＷＡＮＧＢｏ ｔａｏ，ＦＥＮＧＹｏｎｇ ｚｈａｏ，ＹＵＷｅｉ ｊｉａｎｇ
ＤａｔａＰｒｏｃｅｓｓｉｎｇＣｏｍｐａｎｙ，ＣｈｉｎａＯｉｌｆｉｅｌｄＳｅｒｖｉｃｅｓＬｉｍｉｔｅｄ ＣＯＳＬ，Ｚｈａｎｊｉａｎｇ５４２０５１，Ｇｕａｎｇｄｏｎｇ，Ｃｈｉｎａ
Ａｂｓｔｒａｃｔ：Ｆｉｎｉｔｅｄｉｆｆｅｒｅｎｃｅｍｅｔｈｏｄｉｓｔｈｅｍｏｓｔｗｉｄｅｌｙｕｓｅｄｍｅｔｈｏｄｉｎａｃｏｕｓｔｉｃｎｕｍｅｒｉｃａｌｓｉｍｕｌａｔｉｏｎ，ｗｈｉｃｈｈａｓｔｈｅａｄｖａｎｔａｇｅｓｏｆｆａｓｔｃｏｍｐｕｔａｔｉｏｎ，ｒｅｌａｔｉｖｅｌｙｓｍａｌｌｍｅｍｏｒｙｃｏｎｓｕｍｐｔｉｏｎ，ｅａｓｙｐｒｏｇｒａｍｍｉｎｇａｎｄｈｉｇｈｓｉｍｕｌａｔｉ｀ｏｎａｃｃｕｒａｃｙ．Ｃｏｍｐａｒｅｄｗｉｔｈｓｏｍｅｎｅｗｉｍｐｒｏｖｅｄｆｉｎｉｔｅｄｉｆｆｅｒｅｎｃｅａｌｇｏｒｉｔｈｍｓｉｎｒｅｃｅｎｔｙｅａｒｓ，ｆｉｎｉｔｅｄｉｆｆｅｒｅｎｃｅａｌｇｏ ｒｉｔｈｍｂａｓｅｄｏｎｒｅｇｕｌａｒｇｒｉｄｉｓｔｈｅｍｏｓｔｅｆｆｉｃｉｅｎｔｏｎｅｗｉｔｈｔｈｅｌｏｗｅｓｔｍｅｍｏｒｙｃｏｎｓｕｍｐｔｉｏｎ．Ｈｏｗｅｖｅｒ，ｔｈｅｔｒａｄｉｔｉｏｎａｌｄｉｓｃｒｅｔｅｆｏｒｍｕｌａｂｅｃｏｍｅｓｃｏｍｐｌｉｃａｔｅｄｓｉｎｃｅｔｈｅｉｎｔｒｏｄｕｃｔｉｏｎｏｆＰＭＬｂｏｕｎｄａｒｙ．Ｔｈｅｃａｌｃｕｌａｔｉｏｎｐｒｏｃｅｓｓｎｅｅｄｓｔｏｄｉｓ ｔｉｎｇｕｉｓｈｔｈｅｃａｌｃｕｌａｔｉｏｎａｒｅａｆｒｏｍｂｏｕｎｄａｒｙａｒｅａ，ｗｈｉｃｈｒｅｓｕｌｔｓｉｎｌｏｗｃａｌｃｕｌａｔｉｏｎｅｆｆｉｃｉｅｎｃｙ．Ｔｏｓｏｌｖｅｔｈｉｓｐｒｏｂｌｅｍ，ａｎｅｗｄｉｓｃｒｅｔｅｆｏｒｍｕｌａｉｓｄｅｒｉｖｅｄｂｙｔｈｅａｕｔｈｏｒｓ，ｗｈｉｃｈｉｓｓｉｍｐｌｅｉｎｆｏｒｍａｎｄｃｏｎｓｉｓｔｅｎｔｉｎｔｈｅｗｈｏｌｅｓｉｍｕｌａｔｉｏｎａｒｅａ．ＯｐｅｎＭＰｐａｒａｌｌｅｌａｌｇｏｒｉｔｈｍｕｓｅｓｓｉｍｐｌｅｌａｎｇｕａｇｅａｎｄｉｓｗｉｔｈｈｉｇｈｐｏｒｔａｂｉｌｉｔｙ．ＢｙｃｏｍｂｉｎｉｎｇｗｉｔｈＯｐｅｎＭＰｐａｒａｌｌｅｌａｌｇｏｒｉｔｈｍａｎｄｆｕｒｔｈｅｒｏｐｔｉｍｉｚｉｎｇｔｈｅｓｉｍｕｌａｔｉｏｎａｌｇｏｒｉｔｈｍ，ｔｈｅｃｏｍｐｕｔａｔｉｏｎａｌｅｆｆｉｃｉｅｎｃｙｏｆｎｕｍｅｒｉｃａｌｓｉｍｕｌａ ｔｉｏｎｉｓｇｒｅａｔｌｙｉｍｐｒｏｖｅｄ．
Ｋｅｙｗｏｒｄｓ：ｆｉｎｉｔｅｄｉｆｆｅｒｅｎｃｅ；ｒｅｇｕｌａｒｇｒｉｄ；ｄｉｓｃｒｅｔｅｆｏｒｍｕｌａ；ＯｐｅｎＭＰｐａｒａｌｌｅｌｐｒｏｇｒａｍ；ｃｏｍｐｕｔａｔｉｏｎａｌｅｆｆｉ ｃｉｅｎｃｙ
０　引言
有限差分算法是偏微分方程的主要数值解法之
一，是最早的数值模拟方法。

该方法的基本原理是采用离散后的差分算子取代相应的连续微分算子，具有计算效率高、占用内存相对较小的优势，并且
对于近远场及复杂边界都有广泛的适用性，能够准确地模拟波在各种介质及复杂结构地层中的传播规律，是勘探地震学中应用最广泛的数值计算方法［１］。

目前，国内外学者对基于波动方程的地震波正演模拟技术进行了大量的研究工作，地震波正演模拟技术十分成熟。

但在模拟计算过程中由于传统的弹性波动方程数值模拟存在参数繁多、计算量大、耗时长和占用内存空间大等问题，因而声波方程数值模拟得到广泛应用［２］。

１９７１年，Ｃｌａｅｒｂｏｕｔｅｔａｌ．将有限差分法在波场模拟中的应用做了很大改进，提高了模拟效率［３］；１９７４年，Ａｌｆｏｒｄｅｔａｌ．研究声波方程数值模拟，给出精度与差分阶数和网格步长关系［４］；１９７６年Ｍａｄａｒｉａｇａ首次提出交错网格差分算法［５］；１９８４年，Ｖｉｒｉｅｕｘｅｔａｌ．推导了一阶应力 速度方程有限差分算法，提高了计算精度［６］；１９８６年，Ｄａｂｌａｉｎｅｔａｌ．给出声波方程交错网格算法差分格式［７］；１９９７年，Ｃｌａｙｔｏｎ将吸收边界应用到声波方程数值模拟当中［８］；２００１年，Ｃｏｌｌｉｎｏｅｔａｌ．将ＰＭＬ边界应用到应力 速度方程的数值模拟当中［９］；２００４年，Ｓａｅｎｇｅｒ将交错网格进一步发展，提出了旋转交错网格法，这一方法对复杂介质体的模拟精度有显著提高［１０］；２０００年，董良国等详细给出了基于交错网格的一阶应力 速度的弹性波方程有限差分法［１１］；２００４年，裴正林等进一步研究了空间任意偶数阶精度的交错网格差分法，使其更加适应不同复杂介质的数值模拟［１２］。

通过对有限差分声波波场数值模拟研究的分析发现，基于一阶应力 速度方程的有限差分声波模拟是业界应用最为广泛的算法，很对学者对其进行优化，例如优化差分系数的求解方法，使用旋转网格代替常规网格。

虽然这些方法可以提高模拟精度，不容易出现频散现象，但是计算效率会降低，而在工业生产中计算效率更为重要。

通常情况下可以通过使用高阶有限差分算法来提高计算精度，克服频散现象，因此，常规网格下有限差分算法是保证较高的模拟精度条件下，计算效率最高的算法。

随着计算机水平的发展，利用ＯｐｅｎＭＰ并行计算可以大大提高数值模拟速度，但是常规网格下的有限差分算法的ＰＭＬ边界公式复杂，包含三阶偏导数，需要将计算区域分为有效计算区域和边界区域，不利于实现并行计算，容易产生内存错误。

对此，模仿一阶应力 速度方程ＰＭＬ边界公式格式，改进不含有三阶偏导数倒数的ＰＭＬ边界。

保证了计算区域和边界区域计算代码统一，并行计算编程易于实现，对ＯｐｅｎＭＰ的制导语句参数进行测试优化，进一步提高了计算效率。

１　声波方程有限差分算法原理
利用有限差分法进行数值模拟，对空间采用高阶差分离散，而在时间上只采用二阶中心差分，推导出常规网格条件下计算区域与边界区域一致的离散公式，使得ＯｐｅｎＭＰ并行计算过程中效率提高，所以对于方程离散原理，差分系数计算，频散问题不做详细论述。

标准二维声波方程可以表示为：
２Ｐ
ｘ２
＋
２Ｐ
ｚ２
＝
１
ρＶ２
２Ｐ
ｔ２
（１）
式中：Ｐ＝Ｐ（ｘ，ｚ；ｔ）表示声波波场；Ｖ＝Ｖ（ｘ，ｚ）表示模型速度；ρ＝ρ（ｘ，ｚ）表示模型密度。

相应的ＰＭＬ边界条件下边界区域计算方程是包含三阶偏导数的复杂方程组，以右边界为例［１３］：
（
ｔ
＋ｄ（ｘ））２Ｐ
１
＝Ｖ２
２Ｐ
ｘ２
（
ｔ
＋ｄ（ｘ））３Ｐ
２
＝－Ｖ２ｄ′（ｘ）
Ｐ
ｘ
２Ｐ３＝Ｖ２
２Ｐ
ｚ２
Ｐ＝Ｐ
１
＋Ｐ
２
＋Ｐ
３
（２）
其中：
ｄ（ｘ）＝
３Ｖ
ｍａｘ
２Ｌ
（
ｌ
ｘ
Ｌ
）２ｌｇ１
()
Ｒ
（３）
式中：Ｌ为ＰＭＬ区域宽度；Ｒ为反射系数；ｌ
ｘ
为该点到边界的距离。

通过公式（２）可以得出，对于模拟的边界区域要分别判断所属区域，左、右边界使用水平方向的衰减系数ｄ（ｘ），而对于上、下边界需要使用垂直方法的衰减系数ｄ（ｚ），对于四个交点有的学者提出对以对角线进行划分成两部分，上部分使用ｄ（ｚ），下部分使用ｄ（ｘ），但是这些都将导致需要对计算区域进行判定，导致编程实现过程繁琐，不利于并行计算。

通常情况下实现编程计算是通过对方程（１）
７
９
８
第４期 郑如秋，等：基于新离散公式和ＯｐｅｎＭＰ优化的有限差分声波数值模拟
进行离散可以得到空间２Ｎ阶精度的计算区域的波动方程有限差分格式。

Ｐｋ＋１
ｉ，ｊ
＝２Ｐｋ
ｉ，ｊ
－Ｐｋ－１
ｉ，ｊ
＋ρＶ２
Δｔ２
｛１Δｘ２∑Ｎ
ｎ＝１
Ｃｎ［Ｐｋｉ＋ｎ，ｊ－２Ｐｋｉ，ｊ＋Ｐｋｉ－ｎ，ｊ
］　＋１Δｚ２∑Ｎ
ｎ＝１
Ｃｎ［Ｐｋｉ，ｊ＋ｎ－２Ｐｋｉ，ｊ＋Ｐｋ
ｉ
，ｊ－ｎ］｝（４）
式中：Ｐｋ
ｉ，ｊ＝Ｐ
（ｉΔｘ，ｊΔｚ；ｋΔｔ），Δｘ，Δｚ分别表示ｘ方向和ｚ方向的空间采样步长，Δｔ表示时间采样步长，ｉ，ｊ，ｋ分别表示空间和时间的网格数。

而在ＰＭＬ边界区域，由于公式（２）中第二项包含波场的三阶偏导数，需要引入中间变量降阶为二阶偏导数，计算后再计其一阶偏导数，并且公式（２）需要对三个分量方程全部进行离散，使得离散过程十分复杂，增加接近三倍的计算量，会降低并行计算效率，并且，在不同边界需要进行判断，选择使用吸收系数ｄ（ｘ）或ｄ（ｚ）。

Ｃｏｌｌｉｎｏ提出将二阶声波方程转换为一阶应力 速度方程后，再加入ＰＭＬ边界。

这个方法大大降低了直接对二阶声波方程加载ＰＭＬ边界而增加的复杂度。

一阶应力 速度方程加载ＰＭＬ边界条件的控制方程具有如下形式：
ｖｘ
ｔ＋ｄ（ｘ）ｖｘ＝－ Ｐ ｘ ｖ
ｚ ｔ＋ｄ（ｚ）ｖｚ＝－ Ｐ ｚ Ｐｘ ｔ＋ｄ（ｘ）Ｐｘ＝－ρＶ２ ｖｘ ｘ Ｐｚ
ｔ＋ｄ（ｚ）Ｐｚ＝－ρＶ２ ｖｚ ｚＰ＝Ｐｘ＋Ｐ
ｚ
（５）
其中ｖｘ和ｖｚ是质点震动的速度分量，离散公式如下：
　（ｖｘ）ｋ＋１／２ｉ＋１／２，ｊ－（ｖｘ）ｋ－１／２ｉ＋１／２，ｊΔｔ＋ｄ（ｘ）ｉ＋１／２，ｊ（ｖｘ）ｋ＋１／２ｉ＋１／２，ｊ＋（ｖｘ）ｋ－１／２ｉ＋１／２，ｊΔｔ＝－１ρｉ＋１／２，ｊ∑Ｎ
１
ＣｎＰｋｉ＋ｎ，ｊ－Ｐｋｉ－ｎ＋１，ｊ
Δｘ（ｖｚ）ｋ＋１／２ｉ，ｊ＋１／２－（ｖｚ）ｋ－１／２ｉ，ｊ＋１／２Δｔ＋ｄ（ｚ）ｉ＋１／２，ｊ（ｖｚ）ｋ＋１／２ｉ，ｊ＋１／２＋（ｖｚ）ｋ－１／２ｉ，ｊ＋１／２Δｔ＝－１ρｉ＋１／２，ｊ∑Ｎ
１
ＣｎＰｋｉ，ｊ＋ｎ－Ｐｋｉ，ｊ－ｎ＋１
Δｚ（Ｐｘ）ｋ＋１ｉ，ｊ－（Ｐｘ）ｋｉ，ｊΔｔ＋ｄ（ｘ）ｉ，ｊ（Ｐｘ）ｋ＋１ｉ，ｊ＋（Ｐｘ）ｋ＋１ｉ，ｊΔｔ＝－１ρｉ，ｊ∑Ｎ
１
Ｃｎ（ｖｘ）ｋ＋１／２ｉ＋（２ｎ＋１）／２，ｊ－（ｖｘ）ｋ＋１／２ｉ－（２ｎ－１）／２，ｊΔｘ（Ｐｚ）ｋ＋１ｉ，ｊ－（Ｐｚ）ｋｉ，ｊΔｔ＋ｄ（ｚ）ｉ，ｊ（Ｐｚ）ｋ＋１ｉ，ｊ＋（Ｐｚ）ｋ＋１ｉ，ｊΔｔ＝－１ρｉ，ｊ∑Ｎ１Ｃｎ
（ｖｚ）ｋ＋１／２ｉ＋（２ｎ＋１）／２，ｊ－（ｖｚ）ｋ＋１／２ｉ－（２ｎ－１）／２，ｊ
ΔｚＰｋ＋１ｉ，ｊ
＝（Ｐｘ）ｋ＋１ｉ，ｊ
＋（Ｐｚ）
ｋ＋１ｉ，
ｊ
（６）
一阶应力速度方程有限差分算法其相应的ＰＭＬ边界区域离散公式与计算区域离散公式相同，只有阻尼系数在边界区域ｄ
（ｘ）≠０，ｄ（ｚ）≠０，而计算区域阻尼系数为０。

编程过程中不需要对模型进行区域划分，有利于Ｏ
ｐｅｎＭＰ并行计算。

但由于引入中间变量ｖｘ和ｖｚ，相比于常规网线条件下的有限差分算法，增加了大约两倍的计算量。

目前关于数值模拟中ＰＭＬ边界加载方式十分成熟，但是绝大部分都是基于一阶应力 速度方程，交错网格框架下的研究。

笔者研究了一阶应力 速度方程Ｐ
ＭＬ边界离散过程，ＳＰＭＬ和ＮＰＭＬ的推导原理，结合参考文献［１４ １５］，引入二阶阻尼系数，给出直接对常规网格下二阶声波方程加载ＰＭＬ边界的方法：
２Ｐｘ ｔ２＋２ｄ（ｘ） Ｐｘ ｔ＋ｄ（ｘ）２Ｐｘ＝ρＶ２ ２
Ｐ
ｘ２ ２
Ｐｚ ｔ２＋２ｄ（ｚ） Ｐｚ ｔ＋ｄ（ｚ）２Ｐｚ＝ρＶ２ ２Ｐ ｚ２Ｐ＝Ｐｘ
＋Ｐ
ｚ
（７）
以公式（７）中第一项为例，详细地展开为离散形式下的有限差分格式。

离散过程中对于 Ｐｘ
ｔ项
和Ｐｘ
项的离散过程有多重选择，例如 Ｐｘ
ｔ项可以选择向前差分或者向后差分，而Ｐｘ项是否要进行
分解为０ ５×［（Ｐｘ）ｋ＋１＋（Ｐｘ）ｋ－１］。

经过数值计算得出以下结果：
８
９８ 世　界　地　质 第３９卷
２Ｐｘ ｔ２＝（Ｐｘ）ｋ＋１＋（Ｐｘ）ｋ－１－２（Ｐｘ
）ｋ
（Δ
ｔ）２
（８）
经过数值模拟测试
Ｐｘ
ｔ使用向前差分公式，边界反射能量最低：
Ｐｘ ｔ＝（Ｐｘ）ｋ＋１－（Ｐｘ
）ｋ
Δ
ｔ（９）
Ｐｘ项经过数值模拟测试，不需要展开为０
５×［（Ｐｘ）ｋ＋１＋（Ｐｘ
）ｋ－１
］，直接保留为原型，此时来自边界反射能量最小。

最终获得边界反射能量最小的常规网格二阶声波方程离散公式为：
（Ｐｘ）ｋ＋１ｉ，ｊ＝１１＋２ｄ（ｘ）Δ
ｔ（［２＋２ｄ（ｘ）Δｔ－ｄ（ｘ）２Δｔ２］（Ｐｘ）ｋｉ，ｊ
－（Ｐｘ）ｋ－１
ｉ，ｊ
＋ρＶ２Δｔ２Δｘ２∑Ｎ
ｎ＝１Ｃｎ［Ｐｋｉ＋ｎ，ｊ－２Ｐｋｉ，ｊ＋Ｐｋ
ｉ－ｎ
，ｊ］）（１０）计算机实现代码：
／ 波动方程计算函数定义 ／ｖｏｉｄｗａｖｅｑｕ（ｉｎｔｉ，ｉｎｔｊ，ｉｎｔｋ，ｄｏｕｂｌｅ ｃｍ，ｄｏｕｂｌｅ（ ｕ１ｔ０）［ＸＮ］，ｄｏｕｂｌｅ（ ｕ２ｔ０）［ＸＮ］，ｄｏｕｂｌｅ（ ｕ
１ｔｐ１）［ＸＮ］，ｄｏｕｂｌｅ（ ｕ
２ｔｐ１）［ＸＮ］，ｄｏｕｂｌｅ（ ｕ１ｔｎ１）［ＸＮ］，ｄｏｕｂｌｅ（ ｕ２ｔｎ１）［ＸＮ］，ｄｏｕｂｌｅ（ ｕ
１ｔｍｐ）［ＸＮ］，ｄｏｕｂｌｅ（
ｕ２ｔｍｐ）［ＸＮ］，ｄｏｕｂｌｅ（ ｕｔ０）［ＸＮ］，ｄｏｕｂｌｅ（ ｕｔｐ１）［ＸＮ］，ｄｏｕｂｌｅ（ ｄｘ）［ＸＮ］，ｄｏｕｂｌｅ（ ｄｚ）［ＸＮ］，ｄｏｕｂｌｅ（ ｖ）［ＸＮ］）｛
　ｉｎｔｓ；
　ｄｏｕｂｌｅｓｕｍｅｑｕ１，ｓｕｍｅｑｕ２；
　ｆｏｒ（ｓ＝０，ｓｕｍｅｑｕ１＝０ ０，ｓｕｍｅｑｕ２＝０ ０；ｓ＜Ｍ；ｓ＋＋）
　｛
　ｓｕｍｅｑｕ１＝ｓｕｍｅｑｕ１＋ｃｍ［ｓ］ （ｖ［ｉ］［ｊ］ ｖ［ｉ］［ｊ］ ｕｔ０［ｉ］［ｊ＋ｓ＋１］＋ｖ［ｉ］［
ｊ］ ｖ［ｉ］［ｊ］ ｕｔ０［ｉ］［ｊ－ｓ－１］－ｖ［ｉ］［ｊ］ ｖ［ｉ］［ｊ］ ２ ｕ
ｔ０［ｉ］［ｊ］）；　ｓｕｍｅｑｕ２＝ｓｕｍｅｑｕ２＋ｃｍ［ｓ］ （ｖ［ｉ］［ｊ］ ｖ［ｉ］［ｊ］ ｕｔ０［ｉ＋ｓ＋１］［ｊ］＋ｖ［ｉ］［
ｊ］ ｖ［ｉ］［ｊ］ ｕｔ０［ｉ－ｓ－１］［ｊ］－ｖ［ｉ］［ｊ］ ｖ［ｉ］［ｊ］ ２ ｕ
ｔ０［ｉ］［ｊ］）；　｝
　ｕ１ｔｍｐ［ｉ］［ｊ］＝ＤＴ Ｄ
Ｔ／ＤＨ／ＤＨ ｓｕｍｅｑｕ１；
　ｕ１ｔｐ１［ｉ］［ｊ］＝１／（１＋２ ｄｘ［ｉ］［ｊ］ Ｄ
Ｔ） （（２－ｄｘ［ｉ］［ｊ］ ｄｘ［ｉ］［ｊ］ ＤＴ ＤＴ－２ ｄ
ｘ［ｉ］［ｊ］ ＤＴ） ｕ１ｔ０［ｉ］［ｊ］－ｕ１ｔｎ１［ｉ］［ｊ］＋ｕ１ｔｍｐ［ｉ］［ｊ］）；
　ｕ２ｔｍｐ［ｉ］［ｊ］＝ＤＴ ＤＴ／ＤＨ／ＤＨ ｓｕｍｅｑｕ２；
　ｕ２ｔｐ１［ｉ］［ｊ］＝１／（１＋２ ｄｚ［ｉ］［ｊ］ ＤＴ） （（２－ｄｚ［ｉ］［ｊ］ ｄｚ［ｉ］［ｊ］ ＤＴ ＤＴ－２ ｄｚ［ｉ］［ｊ］ ＤＴ） ｕ２ｔ０［ｉ］［ｊ］－ｕ２ｔｎ１［ｉ］［ｊ］＋ｕ２ｔｍｐ［ｉ］［ｊ］）；
　ｕｔｐ１［ｉ］［ｊ］＝ｕ１ｔｐ１［ｉ］［ｊ］＋ｕ２ｔｐ１［ｉ］［ｊ］；
｝
推导出来的离散有限差分公式（１０），在计算过程中计算区域与边界区域计算机实现代码统一，不需要对边界进行判断而改变计算公式，更加有利于并行计算。

推导的公式（１０）与不含边界的有限差分公式（４）相比，只是增加了阻尼系数的计算，这个阻尼系数在速度模型确定时候即可计算，不需要每次都特殊计算。

２　ＯｐｅｎＭＰ多核并行算法
本文提出基于ＯｐｅｎＭＰ及改进的常规网格有限差分声波数值模拟优化，在第１章中主要通过改进边界条件，从算法原理上是代码简洁为并行计算奠定基础，提高计算效率。

本章主要对于ＯｐｅｎＭＰ制导语句、子句等细节进行调整，可以进一步提高计算效率。

ＯｐｅｎＭＰ标准通过定义编译制导、库例程和环境变量规范的方法，基于ｆｏｒｋ ｊｏｒｎ的并行执行模型，将程序主体划分为并行区和串行区，不同线程之间可以通过共享变量实现数据交换，多核ＣＰＵ进行并行运算时，多个ＣＰＵ可同时对程序各个不相关进程进行运算，从而提高程序的运行效率。

ＯｐｅｎＭＰ可直接在串行代码上通过编写并行制导语句实现程序的并行，以良好的简洁性和可移植性成为并行编程的工业标准。

需要并行化的程序如果比较复杂，需要注意程序内部结构和变量之间的逻辑
９
９８第４期 郑如秋，等：基于新离散公式和ＯｐｅｎＭＰ优化的有限差分声波数值模拟
关系，以防止并行后内存读写冲突。

本章主要分析讨论利用Ｏ
ｐｅｎＭＰ技术进行并行设计的两个关键问题：一是如何将原来的串行程序改为并行程序，避免内存读写冲突；二是如何更为有效地使用制导语句，提高并行效率。

声波波场数值模拟过程当中，任意一点下一时刻的波场值的计算是通过当前时刻的波场值、前一时刻的波场值和当前时刻波场在空间方向的导数获得的，这三个变量之间互不相关，计算过程中不会产生影响，可以共享至全部线路，具有可并行性。

但是控制波场点位的变量ｉ，ｊ只能对应于对应线路计算的波场值，需要进行私有变量处理。

公式（１０）所示，计算（Ｐｘ）ｋ＋１
ｉ，ｊ，则需要（Ｐｘ）ｋｉ，ｊ和（Ｐｘ）ｋ－１ｉ，ｊ，因此需要设置转存（Ｐｘ）ｋ－１
ｉ，ｊ＝（Ｐｘ）ｋｉ，ｊ，（Ｐｘ）ｋｉ，ｊ＝（Ｐｘ）ｋ＋１ｉ，ｊ。

波场转存储过程是十
分消耗计算时间的，并且计算机对于内存的访问相比于对于变量的计算十分慢。

对此进行一些简单优化，尽量减少变量数量和对于内存的访问计算，例如：根据控制方程公式（７）中的第三个方程可以通过对于波场Ｐ使用递归调用。

ＯｐｅｎＭＰ使用Ｃ编辑器的＃ｐｒａｇｍａ扩展机制来定义制导，基本格式为：
＃ｐｒａｇｍａｏｍｐｄｉｒｅｃｔｉｖｅ ｎａｍｅ［ｃｌａｕｓｅ［［，］ｃｌａｕｓｅ］…］ｎｅｗ ｌｉｎｅ
其中，每个制导均以＃ｐｒａｇｍａｏｍｐ开始，ｄｉｒｅｃ ｔｉｖｅ ｎａｍｅ是制导名，［］可以加入相应功能的子句。

ＯｐｅｎＭＰ最基本的单元是ｐａｒａｌｌｅｌ结构，由ｐａｒ ａｌｌｅｌ制导。

使用ＯｐｅｎＭＰ中ｐｒｉｖａｔｅ数据子句解决简单变量的内存读写冲突问题。

ｐｒｉｖａｔｅ子句可将简单变量声明为本线程的私有变量，每个线程都存有变量的副本，其他线程无法访问。

即使在并行区域外有同名的共享变量，此共享变量也不会对并行区域造成影响，且并行区域内部计算也不会改变外部共享的变量值。

因此在并行起始语句中加入ｐｒｉｖａｔｅ（ｉ，ｊ）即可以解决简单的内存读写冲突问题。

ＯｐｅｎＭＰ采用ｆｏｒｋ ｊｏｉｎ（分叉 合并）的并行执行模式，在运行过程中，遇到并行开始代码（由ｐａｒａｌｌｅｌ制导），计算机的主线程则调用多个线程来执行并行任务。

如果没有ｎ
ｏｗａｉｔ子句，所有线程在共享计算结构的结束处隐式同步。

为了提高并行计算效率，子句需要使用ｎｏｗａｉｔ，使参与计算的线程无需同步而继续执行随后的程序段，减少多余的同
步以提高计算效率。

Ｏ
ｐｅｎＭＰ的调度模式使用静态调度ｓ
ｃｈｅｄｕｌｅ（ｓｔａｔｉｃ），同时仍然需要使用ｎｕｍ＿ｔｈｒｅａｄｓ（）子句保证并行区域的占用线程不会影响主线程的计算效率，减少线程调度的时间消耗。

通过测试以计算机总线程数为８为例，在并行时尽量不要使用全部线程数量，在并行制导语句中加入ｎｕｍ＿ｔｈｒｅａｄｓ（７），保证并行过程最多使用７个线程，以保证主线程的计算效率。

本文中的有限差分算法，计算每一时刻的波场快照需要通过ｆｏｒ循环结构完成，需要使用ｐａｒａｌｌｅｌｆｏｒ
制导，ＯｐｅｎＭＰ对程序的循环部分进行并行化，要求用于并行的循环部分每次循环互不相关。

二维模拟的计算任务分为三层循环（最外层是时间，内层依次为ｘ
，ｚ方向上的空间循环），由于每个时间切片上的当前循环，依赖于上两个时间切片上的计算结果，即时间循环中各个时间切片存在相关性，因此对于时间不适合使用ＯｐｅｎＭＰ并行而对于ｘ，ｚ方向上的空间循环，具备应用ＯｐｅｎＭＰ并行计算条件，因此每一次空间循环都可以被指定进行并行计算。

使用有限差分算法进行地震波场声波数值模拟区域时，需要在边界加入特殊的边界条件，以消除来自边界的虚假反射问题。

使用上一节推导出的有限差分公式，在空间循环外层加上并行制导指令即可。

并行计算核心代码结构如下所示：
ｆｏｒ（ｋ＝０；ｋ＜Ｋ＿ｒｅｃ；ｋ＋＋）｛　
＃ｐｒａｇｍａｏｍｐｐａｒａｌｌｅｌ　ｎｕｍ＿ｔｈｒｅａｄｓ（７）｛
＃ｐｒａｇｍａｏｍｐｆｏｒｐｒｉｖａｔｅ（ｉ，ｊ）ｓｃｈｅｄｕｌｅ（ｓｔａｔ ｉｃ）　ｎｏｗａｉｔ
　ｆｏｒ（ｉ＝Ｍ；ｉ＜ＺＮ－Ｍ；ｉ＋＋）　｛　ｆｏｒ（ｊ＝Ｍ；ｊ＜ＸＮ－Ｍ；ｊ＋＋）　｛　逐点计算的波动方程计算函数代码 ｝　｝
　波场转存储代码｝　｝
３　数值模拟计算
本章主要为了验证两个问题：①改进的离散的
０
０９ 世　界　地　质 第３９卷
波场模拟公式在任何模型中可以准确模拟波场；②结合ＯｐｅｎＭＰ可以进一步提高计算效率。

目前一阶应力 速度方程时声波模拟的标准方法，为了证明改进的离散公式（１０）可以准确模拟波场，以下的试验将从波场模拟的准确程度和计算效率与其对比。

试验１
：使用一个５００×５００，ｄｘ＝ｄｚ＝１０ｍ，ｄｔ＝８００ｍｓ，速度为４０００ｍ／ｓ的均匀速度模型，震源放置
在中心（２５０，２５０）。

分别使用本文方法进行模拟并与一阶应力
速度方程的模拟结果进行对比。

如图１所示，在ｔ＝０ ５６ｓ，波场未到达边界，模型为均匀模型，因此波场快照表现为圆形，两种方法波形一致。

在ｔ＝０ ８ｓ时，波场触碰到人工边界，均为产生反射现象，说明改进的离散计算公式可以正确模拟波场，并且在边界处无明显反射。

值得注意的是，由于波场传播的控制方程和模拟过程的差
分格式差异，波场在能量振幅存在差异。

ａ．本文方法（ｔ＝０ ５６ｓ）；ｂ．本文方法（ｔ＝０ ８ｓ）；ｃ．一阶应力 速度方程（ｔ＝０ ５６ｓ）；ｄ．一阶应力 速度方程（ｔ＝０ ８ｓ）。

图１　简单模型的波场快照
Ｆｉｇ １　Ｓｎａｐｓｈｏｔｓｏｆｗａｖｅｆｉｅｌｄｏｆａｓｉｍｐｌｅｍｏｄｅｌ
１
０９第４期 郑如秋，等：基于新离散公式和ＯｐｅｎＭＰ优化的有限差分声波数值模拟
为进一步证明推导出的公式（１０）在复杂模型中也可以准确模拟波场，使用Ｓ
ｉｇｓｂｅｅ模型进一步验证。

模型大小为１００１×５０１，ｄｘ＝ｄｚ＝１０ｍ，ｄｔ＝０ ００１ｓ。

速度模型和波场快照如图２所示，ｓｉｇｓｂｅｅ速度模型中包含一个高速盐体，盐体上表面起伏，容易产生绕射波场，将震源放置于模型中间（
５０００，２６００）位置。

理论上，震源上方的崎岖界面会产生一系列绕射波，而盐体的下界面由于界面平缓，界面处速度差异明显，会产生强反射波
和沿着界面滑行的折射波，折射波的波前表现为直线。

选择ｔ＝０ ８ｓ时的波场快照进行展示，两种算法的波场快照中都可以看到震源激发的透射波场波前（绿色箭头所指）、盐体界面产生的反射波（蓝色箭头所指）、折射波（黑色箭头所指）和绕射波（红色箭头所指）。

通过对于ｓｉｇｓｂｅｅ模型进行模拟可以充分说明所推导的离散公式（１０）可以准确
模拟波场。

图２　Ｓｉｇｓｂｅｅ速度模型
Ｆｉｇ ２　Ｓｉｇｓｂｅｅｖｅｌｏｃｉｔｙｍｏｄｅｌ
试验２：
该实验中的程序均为并行程序，使用Ｌｉｎｕｘ系统机群的一个节点进行并行计算效率测试对比。

节点主频２
４ＧＨｚ，十六线程。

分别使用基于一阶应力 速度方程，本文方法公式（１０）和并行代码优化后的程序进行计算效率测试。

测试使用上一个试验Ｓｉｇｓｂｅｅ模型进行模拟，总采样点数设置为１０００，总炮数设置为１炮、２０炮和１００炮。

测试结果如表１所示。

表１　计算效率对比表
Ｔａｂｌｅ１　Ｃｏｍｐａｒｉｓｏｎｏｆｃｏｍｐｕｔａｔｉｏｎｅｆｆｉｃｉｅｎｃｙ
／ｓ
１炮
２０炮１００炮一阶应力—速度方程４２ ２１７６５ ３１３８９３ ９０本文方法公式（１０）３２ ０８５８５ ８４３０３０ ５７并行代码优化
２７ ２１
４８９ ５６
２４００ ３４
从试验结果可以看出，本文改进的有限差分公式（１０）相比于一阶应力 速度方程的计算效率已经有所提高。

而对于并行代码优化后的就算效率进一步提高，计算时间只需要一阶应力 速度方程的６０％。

４　结论
（１）规则网格下的有限差分算法，适当提高差分阶数，既可以满足模拟精度要求，相比于一阶应力 速度方程的计算效率更高。

通过改进边界条件公式，推导出新的离散公式，相比于传统二阶声波方程的ＰＭＬ边界控制方程，边界计算区域中不再包含波场的三阶偏导数和阻尼系数的一阶导数。

并且，新的离散公式在计算区域与边界区域公式一
２
０９ 世　界　地　质 第３９卷
ａ．本文方法；ｂ．一阶应力 速度方程。

图３　ｔ＝０ ８ｓ时的波场快照Ｆｉｇ ３　Ｓｎａｐｓｈｏｔｓａｔｔ＝０ ８ｓ
致，计算机实现过程中不需要加入条件语句判断。

相比于一阶应力 速度方程，减少中间变量引入，模拟的控制方程只需要３个，可以提高计算效率。

（２）通过结合ＯｐｅｎＭＰ并行计算，可以进一步提高计算效率。

并行过程容易产生变量的内存问题，串行代码转换为并行代码后对程序变量优化，减少内存访问。

通过ＯｐｅｎＭＰ制导语句中使用适当的参数，可以进一步加快计算效率。

参考文献：
［１］裴俊勇．基于新差分结构的时空域有限差分逆时偏移研究：硕士学位论文［Ｄ］．成都：成都理工大学，
２０１８．
ＰＥＩＪｕｎ ｙｏｎｇ．Ｔｈｅｒｅｓｅａｒｃｈｏｆｔｉｍｅ ｓｐａｃｅｄｏｍａｉｎｆｉｎｉｔｅ
ｄｉｆｆｅｒｅｎｃｅｂａｓｅｄｏｎｎｅｗｓｔｅｎｃｉｌｆｏｒｒｅｖｅｒｓｅｔｉｍｅｍｉｇｒａ
ｔｉｏｎ：ｍａｓｔｅｒ ｓｄｅｇｒｅｅｔｈｅｓｉｓ［Ｄ］．Ｃｈｅｎｇｄｕ：Ｃｈｅｎｇｄｕ
ＵｎｉｖｅｒｓｉｔｙｏｆＴｅｃｈｎｏｌｏｇｙ，２０１８．
３
０
９
第４期 郑如秋，等：基于新离散公式和ＯｐｅｎＭＰ优化的有限差分声波数值模拟
［２］张千祥．纯Ｐ波各向异性介质波场数值模拟及逆时偏移研究：硕士学位论文［Ｄ］．长春：吉林大学，
２０１６．
ＺＨＡＮＧＱｉａｎ ｘｉａｎｇ．Ｒｅｓｅａｒｃｈｏｎｎｕｍｅｒｉｃａｌｓｉｍｕｌａｔｉｏｎ
ｆｏｒｐｕｒｅＰ ｗａｖｅｉｎａｎｉｓｏｔｒｏｐｉｃｍｅｄｉｕｍａｎｄｒｅｖｅｒｓｅｔｉｍｅ
ｍｉｇｒａｔｉｏｎ：ｍａｓｔｅｒ ｓｄｅｇｒｅｅｔｈｅｓｉｓ［Ｄ］．Ｃｈａｎｇｃｈｕｎ：Ｊｉｌｉｎ
Ｕｎｉｖｅｒｓｉｔｙ，２０１６．
［３］Ｃｌａｅｒｂｏｕｔ，ＪｏｎＦ．Ｔｏｗａｒｄａｕｎｉｆｉｅｄｔｈｅｏｒｙｏｆｒｅｆｌｅｃｔｏｒｍａｐｐｉｎｇ［Ｊ］．Ｇｅｏｐｈｙｓｉｃｓ，１９７１，３６（３）：４６７ ４８１．［４］ＡｌｆｏｒｄＲＭ，ＫｅｌｌｙＫＲ，ＢｏｏｒｅＤＭ．Ａｃｃｕｒａｃｙｏｆｆｉｎｉｔｅ ｄｉｆｆｅｒｅｎｃｅｍｏｄｅｌｉｎｇｏｆｔｈｅａｃｏｕｓｔｉｃｗａｖｅｅｑｕａｔｉｏｎ［Ｊ］．
Ｇｅｏｐｈｙｓｉｃｓ，１９７４，３９（６）：８３４ ８４２．
［５］ＭａｄａｒｉａｇａＲ．Ｄｙｎａｍｉｃｓｏｆａｎｅｘｐａｎｄｉｎｇｃｉｒｃｕｌａｒｆａｕｌｔ［Ｊ］．ＧｅｏｌｏｇｉｃａｌＳｏｃｉｅｔｙｏｆＡｍｅｒｉｃａＢｕｌｌｅｔｉｎ，１９７６，６６
（３）：６３９ ６６６．
［６］ＶｉｒｉｅｕｘＪ．ＳＨ ｗａｖｅｐｒｏｐａｇａｔｉｏｎｉｎｈｅｔｅｒｏｇｅｎｅｏｕｓｍｅｄｉａ：ｖｅｌｏｃｉｔｙ ｓｔｒｅｓｓｆｉｎｉｔｅ ｄｉｆｆｅｒｅｎｃｅｍｅｔｈｏｄ［Ｊ］．Ｅｘｐｌｏｒａｔｉｏｎ
Ｇｅｏｐｈｙｓｉｃｓ，１９８４，１５（４）：２６５ ２６５．
［７］ＤａｂｌａｉｎＭＡ．Ｔｈｅａｐｐｌｉｃａｔｉｏｎｏｆｈｉｇｈ ｄｉｆｆｅｒｅｎｃｉｎｇｔｏｔｈｅｓｃａｌａｒｗａｖｅｅｑｕａｔｉｏｎ［Ｊ］．Ｇｅｏｐｈｙｓｉｃｓ，１９８６，５１（１）：
５４ ６６．
［８］ＣｌａｙｔｏｎＲ，ＥｎｇｑｕｉｓｔＢ．Ａｂｓｏｒｂｉｎｇｂｏｕｎｄａｒｙｃｏｎｄｉｔｉｏｎｓｆｏｒａｃｏｕｓｔｉｃａｎｄｅｌａｓｔｉｃｅｑｕａｔｉｏｎｓ［Ｊ］．Ｂｕｌｌｅｔｉｎｏｆｔｈｅ
ＳｅｉｓｍｏｌｏｇｉｃａｌＳｏｃｉｅｔｙｏｆＡｍｅｒｉｃａ，１９７７，６７（６）：
１５２９ １５４０．
［９］ＣｏｌｌｉｎｏＦ，ＴｓｏｇｋａＣ．Ａｐｐｌｉｃａｔｉｏｎｏｆｔｈｅｐｅｒｆｅｃｔｌｙｍａｔｃｈｅｄａｂｓｏｒｂｉｎｇｌａｙｅｒｍｏｄｅｌｔｏｔｈｅｌｉｎｅａｒｅｌａｓｔｏｄｙｎａｍｉｃｐｒｏｂｌｅｍ
ｉｎａｎｉｓｏｔｒｏｐｉｃｈｅｔｅｒｏｇｅｎｅｏｕｓｍｅｄｉａ［Ｊ］．Ｇｅｏｐｈｙｓｉｃｓ，
２００１，６６（１）：２９４ ３０７．
［１０］ＳａｅｎｇｅｒＥＨ，ＢｏｈｌｅｎＴ．Ｆｉｎｉｔｅ ｄｉｆｆｅｒｅｎｃｅｍｏｄｅｌｉｎｇｏｆｖｉｓｃｏｅｌａｓｔｉｃａｎｄａｎｉｓｏｔｒｏｐｉｃｗａｖｅｐｒｏｐａｇａｔｉｏｎｕｓｉｎｇｔｈｅ
ｒｏｔａｔｅｄｓｔａｇｇｅｒｅｄｇｒｉｄ［Ｊ］．ＳＥＧＴｅｃｈｎｉｃａｌＰｒｏｇｒａｍＥｘ
ｐａｎｄｅｄＡｂｓｔｒａｃｔｓ，１９９９，２１（１）：５８３ ５９１．
［１１］董良国，马在田，曹景忠，等．一阶弹性波方程交错网格高阶差分解法［Ｊ］．地球物理学报，２０００，
４３（３）：４１１ ４１９．
ＤＯＮＧＬｉａｎｇ ｇｕｏ，ＭＡＺａｉ ｔｉａｎ，ＣＡＯＪｉｎｇ ｚｈｏｎｇ，ｅｔａｌ．
Ａｓｔａｇｇｅｒｅｄ ｇｒｉｄｈｉｇｈ ｏｒｄｅｒｄｉｆｆｅｒｅｎｃｅｍｅｔｈｏｄｏｆｏｎｅ ｏｒ
ｄｅｒｅｌａｓｔｉｃｗａｖｅｅｑｕａｔｉｏｎ［Ｊ］．ＣｈｉｎｅｓｅＪｏｕｒｎａｌｏｆＧｅｏ
ｐｈｙｓｉｃｓ，２０００，４３（３）：４１１ ４１９．
［１２］裴正林，牟永光．地震波传播数值模拟［Ｊ］．地球物理学进展，２００４，１９（４）：２１３ ２２１．
ＰＥＩＺｈｅｎｇ ｌｉｎ，ＭＯＵＹｏｎｇ ｇｕａｎｇ．Ｎｕｍｅｒｉｃａｌｓｉｍｕｌａｔｉｏｎ
ｏｆｓｅｉｓｍｉｃｗａｖｅｐｒｏｐａｇａｔｉｏｎ［Ｊ］．ＰｒｏｇｒｅｓｓｉｎＧｅｏｐｈｙｓ
ｉｃｓ，２００４，１９（４）：２１３ ２２１．
［１３］ＫｏｍａｔｉｔｓｃｈＤ，ＴｒｏｍｐＪ．Ａｐｅｒｆｅｃｔｌｙｍａｔｃｈｅｄｌａｙｅｒａｂ ｓｏｒｂｉｎｇｂｏｕｎｄａｒｙｃｏｎｄｉｔｉｏｎｆｏｒｔｈｅｓｅｃｏｎｄ ｏｒｄｅｒｓｅｉｓｍｉｃ
ｗａｖｅｅｑｕａｔｉｏｎ［Ｊ］．ＧｅｏｐｈｙｓｉｃａｌＪｏｕｒｎａｌｏｆｔｈｅＲｏｙａｌ
ＡｓｔｒｏｎｏｍｉｃａｌＳｏｃｉｅｔｙ，２００３，１５４（１）：１４６ １５３．
［１４］ＷａｎｇＴ，ＴａｎｇＸ．Ｆｉｎｉｔｅ ｄｉｆｆｅｒｅｎｃｅｍｏｄｅｌｉｎｇｏｆｅｌａｓｔｉｃｗａｖｅｐｒｏｐａｇａｔｉｏｎ：ａｎｏｎｓｐｌｉｔｔｉｎｇｐｅｒｆｅｃｔｌｙｍａｔｃｈｅｄｌａｙｅｒ
ａｐｐｒｏａｃｈ［Ｊ］．Ｇｅｏｐｈｙｓｉｃｓ，２００３，６８（５）：１７４９
１７５５．
［１５］赵茂强，吴国忱．基于高阶有限差分纵横波分解的弹性波数值模拟［Ｊ］．物探化探计算技术，２０１０，
３２（５）：４８７ ４９４．
ＺＨＡＯＭａｏ ｑｉａｎｇ，ＷＵＧｕｏ ｃｈｅｎ．Ｈｉｇｈｏｒｄｅｒｆｉｎｉｔｅ
ｄｉｆｆｅｒｅｎｃｅｎｕｍｅｒｉｃａｌｓｉｍｕｌａｔｉｏｎｏｆｅｌａｓｔｉｃｗａｖｅｗｉｔｈｓｅｐ
ａｒａｔｅｄＰ ｗａｖｅｓａｎｄＳ ｗａｖｅｓ［Ｊ］．ＣｏｍｐｕｔｉｎｇＴｅｃｈ
ｎｉｑｕｅｓｆｏｒＧｅｏｐｈｙｓｉｃａｌａｎｄＧｅｏｃｈｅｍｉｃａｌＥｘｐｌｏｒａｔｉｏｎ，
２０１０，３２（５）：４８７ ４９４．
４
０
９ 世　界　地　质 第３９卷。