作者简介: 俞岱(1986-),男,博士研究生,主要从事井中地震资料处理和逆时偏移工作。Email:tablemax47@126.com
随着浅层地裂缝等地质灾害问题精细探测需求的增大,采用初至波速度层析反演技术,可以提供较高精度的浅层速度场的纵横向异常变化,可为识别地裂缝等地质灾害现象提供依据。通常,在采集参数一定的条件下,其反演成像精度受算法、网格类型和剖分尺度的限制,在网格类型不变的条件下,要提高精度就必须加密正反演计算网格,改进算法,从而实现增量而不减速且高精度的效果。文中利用OpenMP的单机多核并行技术,探讨并实现了初至波层析成像并行算法,其中正演部分使用了改进的旅行时线性插值法,针对原有的按列计算,加入了按行计算,充分考虑到来自各方向的射线,使得计算出的旅行时最小,反演部分使用了能迅速稳定收敛的LSQR法以及正则化技术,通过理论和实际资料测试,其收敛速度快,反演结果较好。同时,在加密采样和缩小网格尺度的条件下,与常规串行算法相比,其运算速度和效率有较大提高,与商业软件比较,其运算效率差异不大,但反演结果的精度和可分辨性较好。
In the disaster geological problems such as shallow ground fissure etc. the demand for fine detection has been increasing in recent years. The first-break velocity tomography technique can provide a high accuracy of the vertical and horizontal anomalies of the shallow velocity field, and provide the basis for identifying the disaster geological phenomena such as ground fissures etc. In general, under certain conditions of acquisition parameters the accuracy of inversion imaging is limited by the algorithm, the grid type and the scale of the grid. Therefore, in order to improve the accuracy of the calculation and keep the grid type unchanged, the inversion of the grid need to be increased. On the basis of the improved algorithm, the effect of incremental and non-deceleration and high precision can be realized. In this paper, OpenMP-based stand-alone multi-core parallel technology is used, and the first-break tomography parallel algorithm is implemented. In the forward part, the improved traveltime linear interpolation method is used. The method of column calculation is used, and the method of line calculation is also used. In this case, the rays of each direction are taken into account, so that the minimum traveltime is obtained. In the inversion section, the LSQR method and the regularization technique are used. The results of theoretical and practical data test show that the convergence rate is fast and the inversion result is better. At the same time, compared with the conventional serial algorithm, the speed and efficiency of the algorithm are greatly improved under the condition of added sampling and reduced grid scale. Compared with commercial software,the change of computing efficiency is not obvious, but the accuracy and distinguishability of inversion results are better.
层析成像法是20世纪70年代被引入地球物理学领域的, 经历了几十年的发展。20世纪80年代著名学者Anderson[1]、Daily[2]、Somerstein[3]、Bishop[4]等对该技术进行了大量的研究, 其基本原理已经十分成熟(尤其是地震走时层析成像), 而且在之后近十年的发展中, 其精度和分辨率都有很大的提高, 从二维发展到三维, 从各向同性发展到各向异性, 层析成像一直是地震勘探中一种十分重要的技术方法。地震走时层析成像是通过观测到的天然地震或人工地震走时利用模型重建算法反演地下介质结构、速度及其他弹性参数分布的一种反演方法。目前主要用于反演地下介质的速度分布和速度分界面的形态, 被广泛应用于地球内部结构探测、资源与能源勘探与开发、工程勘查与检测、地质灾害调查等领域。
地裂缝灾害是地质灾害中地面变形灾害的一种, 不仅造成各种建筑的直接破坏, 也引起了一系列的环境问题。现代地裂缝现象在世界许多国家都有发现, 其发生频率和灾害程度逐年加剧。自20世纪80年代后期, 物(化)探作为地裂缝研究的主要手段一直起着非常重要的作用
文中针对浅层地裂缝地震精细探测方法技术中初至波速度层析成像算法精度及运算效率等展开分析, 重点就初至波正演计算中的算法改进和OpenMP单机多核层析并行算法等进行了探讨, 经地裂缝理论模型和实际资料测试, 在加密采样和缩小网格尺度的条件下, 与商业软件比较, 反演结果的精度和可分辨性较好。
旅行时线性插值射线追踪法(linear traveltime interpolation, LTI)是由Asawaka于1993年提出的一种基于线性假设的网格单元扩展射线追踪方法, 其过程分为向前处理和向后处理。在向前处理中考虑来自各个方向的射线, 通过线性插值求出各个节点的旅行时, 在向后处理中从接收点开始利用费马原理对射线的传播路径进行追踪。传统的LTI算法存在问题是只考虑来自震源所在半平面的射线, 所以在向前处理中计算出的节点的旅行时可能不是最小旅行时, 这使得某些情况下正演得到的射线与走时不准确。如图1a所示, S点为震源点, R点为向后处理中的某一个交点, 传统LTI算法对R点最小地震波旅行时的计算只考虑了来自激发点所在左半平面的射线, 如图1b所示, 射线在异常复杂的速度结构当中存在这样的可能:即先绕到R点的右半平面再回折到S点。当出现射线回转的情况时, 用传统的LTI算法不能真正得到R点的最小地震波旅行时, 因而向后追踪得到的射线路径也无法满足费马原理。
针对传统的LTI算法只考虑来自震源所在半平面的射线的问题, 改进方法是在向前处理过程中采用列循环和行循环相结合的方法来计算各单元边界节点的旅行时。下面通过模型试算来验证算法的改进效果, 这里先设置一个3层的模型(图2a), 第1层速度为600 m/s, 第2层速度为1 000 m/s, 第3层速度为1 500 m/s, 模型剖分网格大小为1 m× 1 m; 观测系统:震源位于(0, -10), 接收点位于(0, 0), (4, 0), …, (200, 0), 共51个接收点。传统LTI法追踪的射线路径如图2b。把模型旋转90° , 旋转后的模型如图3a所示, 相应的观测系统变为:震源位于(10, -200), 接收点位于(0, -200), (0, -196), …, (0, 0), 共51个接收点。对旋转后的模型用改进前、后两种射线追踪方法得到的射线路径如图3b、3c所示。对比图2b和图3b、3c中的射线路径和部分接收点旅行时(表1)可知, 旋转后模型用传统LTI法正演只能得到直达波, 这是没有考虑来自各个方向射线的结果, 旋转后模型用改进的LTI法计算得到的射线路径和旅行时与旋转前模型得到的是一致的, 改进后的LTI算法对复杂模型的适用性更强。
![]() | 表1 模型旋转前后部分接收点旅行时计算对比 |
初至波旅行时层析反演, 首先是建立初始模型进行地震初至波走时正演, 此时需要将模型空间进行离散化, 假设每个网格单元内的慢度s为常数, 确定好观测系统之后可以计算每一个炮检对的初至波射线路径以及走时, 将计算的走时与实际观测到的走时做差, 残差为Δ T, 根据广义Radon变换, 假定慢度扰动很小, 慢度模型与扰动慢度模型中射线路径相同, 得到走时残差层析成像公式为AΔ s=Δ T, 其矩阵形式为:
其中:A是M× N维矩阵, 矩阵元素aMN代表第M条射线在第N个网格单元内的长度, Δ sN代表第N个网格单元的慢度扰动量, Δ TM代表第M条射线走时的残差。
传播矩阵A的构建是反演问题的关键, 笔者利用旅行时线性插值法正演可以得到初至波的传播矩阵以及初至波走时。式(1)通常是一个大型稀疏矩阵, 文中选择相对节省内存(利用大型稀疏矩阵的压缩存储技术)与计算量, 且能迅速稳定收敛的LSQR算法对方程进行求解。由于旅行时层析反演的不适定性, 层析问题的解通常含有某些不符合实际的高波数振荡, 因此文中在层析过程中加入正则化并对层析结果进行平滑, 这里正则化主要使用Laplace权值:
选取Laplace权值即中心网格单元的值与周围8个网格的值的总和相等, 与层析反演中最小平方的思想结合, 就相当于在正则化作用的区域范围内慢度值不发生大的变化。
为了测试文中初至波旅行时层析反演程序的正确性, 建立一个起伏地表层状正断层模型进行测试, 模型如图5, 模型信息与观测系统详见表2。
![]() | 表2 理论模型信息与观测系统信息 |
根据此模型和观测系统利用本文介绍的初至波旅行时线性插值法, 正演得到初至波走时, 作为初至波层析反演输入的观测走时。如图6所示为第1炮(坐标(4.00 m, -24.88 m))和第21炮(坐标(243.94 m, -25.67 m))初至波走时曲线对比、第21炮的初至波走时曲线在第37道附近有明显台阶, 这是由于该炮射线穿过断层导致的。建立如图7所示的反演初始速度模型, 模型横向与纵向的网格间距为2 m× 1 m, 从地表线起, 给初始速度740 m/s, 以纵向每5 m增加80 m/s的梯度依次递增直到 2 500 m/s为止。将以上的走时和速度模型作为输入, 利用本文层析反演程序, 得到如图9所示的反演结果, 图8为反演速度收敛曲线, 可见随着迭代次数的增加, 速度差依次递减, 迭代次数为20次时, 反演基本上收敛, 相邻两次的速度差仅为15 m/s。将理论模型的层位投到初至波层析反演的结果上, 第1层的速度分布与层位线对应性很好, 断层在第1层上的位置清晰可见, 第2层的反演结果次之, 断层在第2层上的位置较模糊, 但还可以区分出来, 第3层的速度结果以难区分断面位置, 通过射线分布图可知, 射线最深到达第3层界面附近, 浅层效果好于深层, 由于模型深层和左、右两侧射线较稀, 因此这些地方的反演结果稍差。为了测试加密采样对反演结果的影响, 在前述模型和观测系统基础上, 将道距加密到2 m, 以同样的初始模型又重新反演了一个结果, 如图10所示, 采样加密后的结果断面在浅层的刻画效果以及速度的成层性更好, 成像精度更高。总体而言, 本反演方法的效果较好, 达到了设计要求。
2.2.1 性能评价指标
并行算法的设计主要是为了加快程序的计算速率, 并有多个标准用于评判。通常用加速比和效率来评价并行算法的性能[19]。
1)加速比。多处理器下计算速度提升的倍数就是加速比, 即代码在单处理器下与多处理器下运算时间之比。计算公式为:
式中:SP为加速比, T1是代码在单处理器下运算的时间, TP是指代码在多处理器下运算的时间。
2)效率。效率用于说明并行的效果。计算公式为:
式中:EP为效率, 一般小于1; SP为加速比; P为所使用的线程数。
2.2.2 并行效果分析
利用前述的理论模型资料, 在如表3所示的开发测试平台上进行本文介绍旅行时层析反演程序的并行算法效率评价。选择2 m× 1 m的网格间距, 所用线程数从1~4, 计算出各线程运行时间(为减小其他因素影响, 运算3次求均值)、并行加速比、并行效率, 详见表4。
![]() | 表3 开发测试平台配置 |
![]() | 表4 理论资料运算分析结果 |
从表4中可见, 随着线程数的增加, 程序的运行时间减小, 将串行计算和4线程并行计算的运行时间进行比较, 如图11所示, 对道距4 m资料反演串行计算的运行时间为495.548 s, 而4线程并行计算的时间为249.771 s, 通过并行计算程序的运行时间大幅缩减, 加速比达1.98。如图12所示, 随着线程数的增加, 加速比增加, 效率减小, 加速比并不是按照线性递增, 由于效率随着线程的增加而减小, 导致加速效果有折扣, 整体上来说并行带来的速度提升明显。对比道距4 m时串行和道距2 m时4线程并行的计算时间, 在加密采样点数前提下, 通过并行技术达到了增量而不减速且高精度的效果。
2.2.3 对并行算法效率的影响
超线程技术, 简称HT(Hyper-Threading)技术, 由英特尔公司于2002年提出并研发[20], 现为主流化的技术, 它是利用特殊的硬件指令, 将一个处理器核模拟为两个逻辑内核[21], 让每个单独的处理器都能运用线程级并行计算, 从而减少CPU的闲置时间, 提高CPU的运行速度。英特尔公司表示, 采用超线程技术, 处理器仅增加5%的固件面积, 就可以换来15%到30%的效能提升[20]。在实际使用过程中, 超线程会在某些程度上降低效能, 因为它是两个线程共用单个处理器的资源, 因此当出现资源占用冲突时, 其中一个线程就要临时退出资源, 直到闲置后才能够继续使用[22]。因此, 超线程得出的线程性能并不等价于两颗CPU的性能。
本次测试平台为双核四线程的计算机, 其具有超线程技术。因此有必要测试超线程技术对于本文介绍的并行算法在效率方面的影响。本平台在关闭超线程技术后, 就变为双核二线程计算机, 因此并行运算时最多只进行2线程并行运算。利用前述的道距4 m的理论资料, 在测试平台关闭超线程技术之后对其进行串行和2线程并行运算, 计算出程序运行时间、并行加速比以及并行效率, 详见表5。对比可知:关闭超线程技术后, 2线程并行运算的加速比和效率均高于打开时2线程并行的值, 而略低于4线程并行运行的值, 将关闭超线程下2线程并行的值换算到4线程情况下, 加速比为1.8, 效率为45%, 性能基本与打开超线程下4线程并行时相当, 说明针对计算量大的并行程序时超线程技术优势并不明显, 而且从串行计算的运行时间来说, 关闭超线程技术时运算更快, 建议在运算量较大的工作站或服务器上可将超线程技术关闭, 因为其对提高并行运算效率帮助不大, 而且在关闭超线程技术时处理器进行串行计算可最大发挥运算性能, 毕竟算法程序不太可能做到完全并行计算, 其中会有一些串行计算存在。
![]() | 表5 打开/关闭超线程技术情况下理论资料运算分析结果 |
为了进一步验证本文层析反演的效果, 选取陕西某地实测地裂缝地震勘探资料对程序进行测试。该资料为查明可能存在的地下地裂缝分布以及构造的起伏变化。工区位于关中盆地, 区内浅表层大部分为黄土覆盖, 局部有基岩出露, 浅表层吸收衰减较大且地震地质条件良好。针对工区实际, 通过试验确定二维地震采集参数为:震源为炸药震源, 药量3 kg, 炮孔深度5 m; 激发方式为单边放炮, 测线长 1 305 m, 道距10 m, 变道数接收, 炮数22炮, 采样率0.5 ms, 记录长度为1.0 s。
测线分布如图13中左图所示, 测线较直, 故可利用二维层析反演程序试算。通过对实际资料中初至波走时的拾取(如图13中右图所示), 使用本文的方法和Green Mountain层析软件进行初至波旅行时层析结果比较。图14左边为本文反演方法的速度收敛曲线, 可见其收敛性较好。图14右边为反演的初始模型示意图。
图15为Green Mountain软件的反演结果, 图16为本文方法得到的结果, 根据射线的分布将无射线部分切除, 对反演结果进行解释, 两种反演结果形态大体一致, 本文反演结果的射线分布更广, 有效范围更大, 速度分布成层性更好, 对地裂缝的可分辨性较好。通过初至波层析反演结果可以解释出地下地裂缝分布, 本方法的实际资料应用效果较好。
在本测试平台上将实际资料进行串行和并行(4线程并行)运算时间对比, 串行时间为259.133 s, 并行时间为129.453 s, 说明并行计算时运行时间大幅缩短, 在实际工作中可以有效地节省计算时间, 提高工作效率。
本文介绍的初至波层析成像并行算法, 其中正演部分使用了改进的旅行时线性插值法, 针对原有的按列计算, 加入了按行计算, 充分考虑到来自各方向的射线, 使得计算出的旅行时最小, 反演部分使用了能迅速稳定收敛的LSQR法以及正则化技术。
本方法在地裂缝勘查工作中的应用效果较好, 与专业Green Mountain软件结果对比, 其射线分布更广, 有效范围更大, 速度分布成层性更好, 对地裂缝的可分辨性较好。
通过对测试资料进行并行效率分析可知, 并行算法可以有效的提高程序计算效率, 建议在运算量较大的工作站或服务器上可将超线程技术关闭, 因为其对提高并行运算效率帮助不大。
The authors have declared that no competing interests exist.
[1] |
|
[2] |
|
[3] |
|
[4] |
|
[5] |
|
[6] |
|
[7] |
|
[8] |
|
[9] |
|
[10] |
|
[11] |
|
[12] |
|
[13] |
|
[14] |
|
[15] |
|
[16] |
|
[17] |
|
[18] |
|
[19] |
|
[20] |
|
[21] |
|
[22] |
|