2010, Vol. 53
2. 油气资源与探测国家重点实验室, CNPC物探重点实验室, 北京 102249
2. College of Geophysics and Information Engineering, China University of Petroleum, Beijing 102249, China
地震数值模拟是地震勘探和地震学的重要基础，三维地震资料的处理和解释都需要有效的三维正演模型予以验证.在地震勘探领域有两个重要的地震模型已被地球物理领域科研人员用于许多相关的科研项目，其一是Marmousi模型，1988年，由法国石油研究院根据西非安哥拉真实地质情况做出了二维模型地震数据体[1]，其二是美国SEG/EAEG于20世纪90年代末，根据墨西哥湾Sigsbee地区的地质情况，给出了三维盐丘模型的地震数据体[2, 3].我国目前还没有根据实际地质资料给出的完整的三维地震数据体，本文所作的主要研究工作是根据库车地质模型给出了三维正演数据体.在选择正演数值计算方法上，经过前期研究并从计算速度、稳定性能、易于并行等方面考虑，选用了任意差分精细积分方法进行库车地区的地震数值模拟.
任意差分精细积分方法（简称ADPI）是一种半解析的解数理方程的数值方法.对时间域采用局部积分的解析方法，与有限差分对时间和空间进行离散格式相比，在少量增加计算量的情况下可较大地提高计算精度[4].1995年，钟万勰和朱建平[5]用单点精细积分法改善了有限差分法造成的数值不稳定缺点；1999年，强士中等[6]提出了改善热传导方程的时间精度的方法；2003年，贾晓峰等[7]导出了任意差分精细积分法的公式并对地表模型进行正演，提高了差分精度；2007年，王润秋等提出了用精细积分法消除数值计算中有限波场区域边界引起的边界反射[8].
根据三维任意差分精细积分方法计算方便，边界条件易控制的特点，我们采用OpenMP宏指令并行算法理论，进行了大规模的地震波场模拟，共采集了300炮全方位的地震数据体，整个数据体存储量约为2T，本文所做的研究工作既为库车地区石油勘探提供基础地震数据体，同时也对大规模地震波场模拟并行计算效率进行了数值试验.
2 三维数值模拟计算格式 2.1 数值计算公式任意差分精细积分法采用任意差分的方法计算空间偏导数，采用子域精细积分的方法计算时间偏导数.对于三维波动方程[9]：
|
(1) |
首先对方程（1）进行三维空间坐标离散，在任一离散点j（x，y，z）（点序号）处，有：
|
(2) |
将u在j点的x方向邻域作Taylor展开：
|
(3) |
（2）式中，Δxi=xi-xj，i是j的x邻域内的任意点，对（3）式移项并加权得：
|
(4) |
式（4）中αi是加权系数，m原则上可任意选取，一般取m=n.对比（4）式和（2）式得到：
|
(5) |
（5）式是一个关于αi的m维线性方程组，解之可求出加权系数αi，从而得到下列计算格式（其详细推导见文献[10]）
|
(6) |
其中，
实验使用了二阶Enguist[11]同轴近似吸收边界条件，吸收效果不尽人意.因此改进吸收边界条件，与原先相比有了很大程度的提高.
改进分两个方面：
（1）对Enguist方法的一阶旁轴近似进行扩展，在入射波为平面波，且入射角θ已知的情况下，通过在微分方程中边界点处引入系数cosθ，依据改进的微分方程计算边界，可将随角度θ的入射波完全吸收.
（2）对吸收边界上每一个格点，假设在这局部范围内可将入射波近似看成平面波，利用波场导数估算该点附近的入射方向.利用估算得到的入射波角度值θ，应用下列方法计算边界：
|
(7) |
式中，a=10/（3π），b=2/（3π）.
在匀速模型下，中心线左端2/5处放置点震源激发后反射波场强度对比见图 1.
|
图 1 改进吸收边界条件效果图 (a)原方法；(b)改进的方法. Fig. 1 The effort figure of improved absorbing boundary conditions (a) Original method; (b) Improved method. |
从反射波前的峰值来比较，改进型的反射波峰值仅是原方法的20%~30%.
2.3 数值稳定性试验稳定性是衡量地震波数值模拟方法优劣的重要标准之一.在实际计算过程中，减小空间和时间步长Δd和Δt在一定程度上能够提高计算精度.但是这种提高是以牺牲空间资源和效率为代价的，因此要适当选取空间和时间步长.
从数值试验来看，需要根据计算域的尺寸，介质速度特性等设置，依照三维库车速度模型的情况，Δd=16m，V=5200m/s，最后确定Δt=2/3ms=0.667ms.此时计算过程稳定.
3 正演模拟资源占用以及计算速度库车工区原始3D数据块的网格总节点数为1001×813×2001（依次x，y，z方向），文件大小为6.24GB，空间步长分别为16m，16m，4m，此时程序计算需要的内存大小不少于1001×813×2001×4（float型是4个字节）×4（速度场，上一时刻波场、当前时刻波场、下一时刻波场）=26GB，内存要求较大.
为了在单CPU上估算所需的计算时间，将原始速度模型在3个方向上各采样点数为原来的1/4，得到了一个100MB大小的速度数据.用此数据进行ADPI正演模拟得到表 1结果.
|
表 1 单机模拟数据 Table 1 Single machine simulated data |
由表 1可以看出ADPI方法在提高了时间和空间精度的同时，以增加了大概8倍左右的计算时间作为代价.而由此估算得到单核计算原始的速度数据需要的时间为20×64=1280h.这显然是不能接受的.
我们采取了以下措施来解决上面的问题：
（1）用并行工作站作为计算设备：曙光EP-850工作站，64bit WinServer操作系统，64G物理内存，8CPU（Quad-Core AMD Opteron 8354 2.2 GHz），总共32Core.
（2）代码优化和并行化：应用了部分代码优化的准则（例如手工展开循环，内存地址对齐与连续存取）；使用了Intel编译器取代VS2008自带编译器，并设定根据当前处理器类型进行优化（20%性能提升）；使用OpenMP宏指令实现程序并发，其效果优于人为多线程实现，在32Core的机器上经测试能够获得相当于串行程序的20倍加速.
（3）考虑到计算时间与正演分辨率的折中，我们最终采用的模型大小是1001×813×501（z方向采样点数为原来的1/4），空间步长则一致为16m.这样可以将计算量减少为1/4.每一炮的计算时间约为14.6h.真实速度数据正演性能如表 2所示.
|
|
表 2 曙光EP-850并行工作站模拟数据结果 Table 2 The results of Dawning EP-850 parallel data workstation simulation |
为了充分反映库车前陆盆地起伏地表、复杂构造的典型地质特征，使设计的三维地质模型能作为西部地区的一种典型的地质构造，所获得的模拟地震数据能对这一地区的地震资料处理和解释提供帮助.首先对库车地区实际地质地貌以及取得的数据进行分析，最终选择拜城地区典型地震剖面作为背景资料（见图 2），进行关键测线位置的地质模型设计，对地层及其倾角、断层等作了适当的调整（见图 3）.
|
图 2 拜城地区典型地震剖面 (a)、(b)、(c)为该区块南北侧线地震剖面构造解释图；(d)为东西向侧线地震剖面构造解释图. Fig. 2 The typical seismic profile of Baicheng region (a)、(b)、(c) Profiles of this region innorth-south lateral line; (d) Profile of this region ineast-west lateral line. |
|
图 3 关键测线地质模型设计 (a)、(b)、(c)为该区块南北侧线剖面图，图(d)为东西向侧线剖面图. Fig. 3 The design of the critical line of the geological model (a)、(b)、(c) Profiles of this region in north-south lateral line; (d) Profile of this region ineast-west lateral line. |
模型充分反映了库车地区复杂地表和复杂构造特征：大型逆冲推覆构造，大倾角地层，高角度的断层，高陡构造特征.浅层设计高陡构造层，深部设计为平缓构造，不整合面发育.模型浅部有两个逆断层，深部有3个逆断层.地表起伏控制在800m左右；地下构造模型共设计为7层（不包括起伏地表），纵向深度大于8000m，横向分别为13000m、16000m；模型速为3300~5300m/s（见图 4）.
|
图 4 3D地质模型 (a)、(b)、(c)、(d)为不同角度观测三维模型，模型共设计为六层; 东西方向长13000以，南北方向长16000以，深度为8000:m Fig. 4 3D geological model (a)、(b)、(c)、(d) are 3D models viewed from different direction. the model has 6 layers, 13000 m from east to west; 16000 m from north to south; and 8000 m in depth. |
各地层的速度参考实际地层的速度调整得到，其值如表 3.
|
|
表 3 三维模型的设计速度 Table 3 The designed speed of three-dimensional model |
模型速度体积16000m×13000m×8000m，16m为网格划分单位，三维速度模型x、y和z方向上网格划分成1001×813×501，y方向的中间测线作为炮线，地表接收.模拟起始炮点位置在网格点（9，406，0）位置，最后一炮炮点位置（999，406，0），每一炮炮点沿炮线位置逐次移动10个网格点；正演模拟100炮×3条测线，每一炮都是全地表（1001×813点）接收.
三维速度模型：x方向长16km，y方向长13km，z方向长8km.速度模型x、y和z方向都采用16m网格划分，网格点为1001×813×501.对其进行三维地震任意精细积分波场数值模拟，所有炮点在y方向的中间一条测线（第406条测线）上，沿x方向布置.炮点初始位置在（9，406，0），最小炮距为10个网格点，y方向上有813条测线接收.采样间隔2ms.每一炮计算时间大约在14h38min，每一炮的数据存储量为6.06G.图 5为炮点在（489，406，0）位置上的8条测线的全三维单炮记录，其中的每条测线在x方向上相隔50条测线.
|
图 5 全三维单炮记录 Fig. 5 3D single shot records of 8 test lines at the point (489, 406, 0) |
图 6所示为12个不同炮点的第406条测线上的单炮记录.炮点位置在y方向上顺序相隔50条测线.三维库车速度模型经过降采样后的数据是1.51G，模拟每一炮生成的纯波场数据量是6.06G，每一炮的计算时间大约14h38min，在并行工作站昼夜24h运行的情况下，整个前陆盆地库车三维速度模型300炮正演耗时180天左右.
|
图 6 第406条测线(即炮线）上12个不同炮点的单炮记录 Fig. 6 Single shot records of12 different shot points at No. 406 test line (shot line) |
经与地质模型和实际野外采集地震记录比对，三维ADPI正演模拟可以很好显示地质模型中的明显构造以及分层信息，在较复杂区域能够很好描述地下地质构造，振幅、频率、相位等地震主要物理参数可以达到野外防真效果.
6 结论三维任意差分精细积分法计精度高，稳定性好，又能较好地描述非均匀介质的物理特性及灵活处理各类边界条件，是一种有效的计算地震波场并进行地震正演模拟的方法.由于三维勘探数据量巨大，所以正演模拟耗时长.本文采用任意差分精细积分和OpenMP宏指令并行方法正演模拟，大大减少了运算时间，而且精度和稳定性也得以保证.
通过对三维精细积分并行算法的研究，我们得到了以下结论：
（1）三维任意差分精细积分算法在时间域和空间域都采用局部积分，精度高、稳定性强，是一种高效的三维地震正演方法；
（2）采用Intel编译器使大数据在VisualStudio上得以实现，OpenMP宏指令并行提高了计算速度，耗时相比较使用单机计算减少了1/20左右；
（3）三维正演边界处理上采用改进的自适应吸收边界，相比较普通的吸收边界，界面更清楚，反射减弱；
（4）能够较为精确地模拟复杂地下介质构造波场，是一种大规模数据模拟防真的实用正演算法.
| [1] | Grary S. Martin, Robert Wiley, Kurt J. Marfurt.Marmousi2:An elastic upgrade for Marmousi. The Leading Edge , 2006, 25(2): 156-166. |
| [2] | Paul Sava. Subsalt Exploration and Development Imaging:Interpretation, and Drilling-What have we learned. The Leading Edge , 2006, 25(11): 1370-1376. DOI:10.1190/1.2399267 |
| [3] | 刘礼农, 刘洪, 李幼铭. SEG/EAGE盐丘和推覆体模型的波动方程三维叠前深度偏移成像. 地球物理学报 , 2004, 47(2): 312–320. Liu L N, Liu H, Li Y M. Wave-equation 3-D prestack depth migration for the SEG/EAGE salt and overthrust model. Chinese J.Geophys. (in Chinese) , 2004, 47(2): 312-320. |
| [4] | 孙建国. 复杂地表条件下地球物理场数值模拟方法评述. 世界地质 , 2007, 26(3): 345–357. Sun J G. Methods for numerical modeling of geophysical fields under complex topographical conditions:a critical review. Global Geology (in Chinese) , 2007, 26(3): 345-357. |
| [5] | 钟万勰. 子域精细积分及偏微分方程数值解. 计算结构力学及其应用 , 1995, 12(3): 253–260. Zhong W X. Subdomain precise integration method and numerical solution of partial differential equations. Computational Structural Mechanics and Applications (in Chinese) , 1995, 12(3): 253-260. |
| [6] | 强士中, 王孝国, 唐茂林, 等. 任意差分精细积分法及数值稳定性分析. 应用数学与力学 , 1999, 20(3): 256–262. Qiang S Z, Wang X G, Tang M L. On the arbitrary difference precise integration method and its numerical stability. Applied Mathematics and Mechanics (in Chinese) , 1999, 20(3): 256-262. |
| [7] | 贾晓峰, 王润秋, 胡天跃. 求解波动方程的任意差分精细积分法. 中国地震 , 2003(3): 236–242. Jia X F, Wang R Q, Hu T Y. The Arbitrary Difference Precise Integration Method for Solving the Seismic Wave Equation. Earthquake Research in China (in Chinese) , 2003(3): 236-242. |
| [8] | 王润秋, 张明, 张崴. 三维精细积分法地震正演及边界条件. 石油地球物理勘探 , 2007, 42(2): 202–207. Wang R Q, Zhang M, Zhang W. 3-D fine integration seismic forward modeling and boundary condition. Oil Geophysical Prospecting (in Chinese) , 2007, 42(2): 202-207. |
| [9] | 牟永光, 裴振林. 三维复杂介质地震数值模拟. 北京: 石油工业出版社, 2005 . Mou Y G, Pei Z L. Seismic Numerical Modeling for 3-D Complex Media (in Chinese). Beijing: Petroleum Industry Press (in Chinese), 2005 . |
| [10] | Toldi J L, Hale D. Data-dependent absorbing side boundaries. Stanford Exploration Project Report, 1982.111~121 |
| [11] | Clayton R W, Engquist B. Absorbing boundary conditions for wave equation migration. Geophysics , 1980, 45: 895-904. DOI:10.1190/1.1441094 |