姚姚，奚先，唐文榜，许云，黄延诂，

姚姚，奚先，唐文榜，许云，黄延诂，_大人地球物理科学

姚　姚¹，奚　先¹，唐文榜²，许　云³，黄延诂¹， 吴俊峰⁴，秦　臻⁵，周玉冰⁶，董桥梁⁷，董雪华⁸

（1.中国地质大学地球物理与空间信息学院，武汉，430074；

2.中国石化集团公司勘探开发研究院，北京100083；

3.中国地质大学地球物理与空间信息学院，北京，100083；

4.中国石化集团公司江汉油田分公司勘探开发研究院，潜江，433124；

5.中国石油集团公司勘探开发研究院，北京，100083；

6.中国石油集团公司海外分公司，北京，100083；

7.CGG公司中国分部，北京，100083；

8.中国石油集团公司勘探开发研究院西北分院，兰州，730020）

摘　要：地震正演数值模拟无论在地震波场的理论研究中，还是在地震勘探、开发的实际应用中都具有重要的意义。系统地介绍了作者多年来对常规非均匀介质地震数值模拟方法的改进及应用情况，包括一维非线性褶积模型正演模拟、二维射线追踪数值模拟边缘绕射波计算、非均匀介质随机介质模型数学建模、弹性波波动方程数值模拟计算和在碳酸盐岩孔洞型储层预测中的应用、细胞自动机方法地震波数值模拟等。

关键词：地震正演，数值模拟，非均匀介质

前　言

地震数值模拟是根据野外地质体的形状和介质参数建立起数学模型，用数值求解的方法将野外地震勘探过程在数字计算机中重演的一种地震正演模拟方法。

地震数值模拟的方法很多。从模型的维数上说，有一维模拟、二维模拟和三维模拟等；从模拟的地质对象种类上说，有非均匀介质模拟、各向异性介质模拟、裂隙介质模拟等；从数值计算方法上说，有合成地震记录计算、射线追踪计算、波动方程数值解等，每一类还包含有多种不同的方法。

作者多年来对常规非均匀介质地震数值模拟方法进行了系统的改进及应用，从一维非线性褶积模型正演模拟^［1～3］、二维射线追踪数值模拟边缘绕射波计算^［4～7_］、非均匀随机介质模型数学建模^［8～16］、弹性波波动方程数值模拟计算和在碳酸盐岩孔洞型储层预测中的应用，到细胞自动机方法地震波数值模拟等。本文将对比作系统的介绍，由于受篇幅所限，有关细节请读者参考本文后面所附的参考文献。

1　一维非线性褶积模型正演模拟

最简单的地震数值模拟是形成一维合成地震记录。常规方法是用地震子波与由声测井资料得到的反射系数序列褶积产生一维合成地震记录。即：

写成Z变换的形式为：

形成这种合成地震记录的地质模型是一维线性褶积模型。它是一种理想的层状构造模型，称为古皮特模型。它假设：①传播介质由大量不同弹性性质的“薄层”构成，薄层假设为等时厚，每层的时间厚度均为半个采样间隔；②弹性振动以平面纵波形式垂直入射，即对于模型各层均是法线入射；③每一界面的反射子波形状都相同，即忽略介质的吸收作用；④地震波由上方入射到界面上时完全反射没有透射作用，而由下方入射到界面上时完全透射没有反射作用。

由此可见，上述层状介质线性褶积模型中的界面只是一种自然界中并不存在的理想数学界面，它忽略了实际介质中总是存在着的透射波和多次反射波。当一维合成地震记录考虑透射波和多次反射波的影响时，层状构造模型就不再是线性褶积模型了，而是一种非线性褶积模型。

考虑层状介质的每一层中应该既有上行波，又有下行波。如图1所示，第n层介质顶和介质底上、下行波的Z变换关系为：

第n层介质底界面即第n界面上面有第n层底的上、下行波，下面有第n＋1层顶的上、下行波，如图2所示。方向各不相同的4个波通过上方入射反、透射系数r′_n和t′_n以及下方入射反、透射系数r_n和t_n相互联系，其Z变换形式为：

将式（3）代入式（4），可以得到第n＋1层顶的上、下行波与第n层顶的上、下行波之间的关系：

图1　第n层中上、下行波

图2　第n界面上、下各波间关系

若层状介质有N层，对于整个N层构造可以得到：

式中：T_N＝t₁·t₂…t_N-1·t_N

式中：F_N（Z）与G_N（Z）都是最高阶为N-1的Z的多项式，而F_N（Z^-1）与G_N（Z^-1）则是它们的复共轭多项式。式中f_i和g_i可以用递推的方法由各界面反、透射系数得到。

若一平面纵波垂直入射到N层构造层状介质［D₁（Z）＝S（Z），U_N＋1（Z）＝0］，则地面上接收到的反射波［Y₀（Z）＝U₁（Z）］为：

因为考虑透射波和多次反射波的反射地震记录的Z变换不是地震子波的Z变换与反射系数序列的Z变换多项式相乘，而是与一个不具有线性关系的有理分式相乘，故称这一模型为层状介质非线性模型。

利用多项式除法可以将上述有理分式展开。一般而言，可以得到一个无穷级数，该无穷级数中各项的系数称为拟似反射系数。因此，包含所有一次反射、多次反射和透射影响的合成地震记录的计算应是地震子波与拟似反射系数序列的褶积，这就是一维非线性褶积模型正演模拟。

拟似反射系数序列的Z变换是一个无穷级数反映了实际地层中多次反射是永无止境的。当然，实际计算时不能取无穷多项，只能截断使用。这也符合实际，因为实际地层中经过太多次数的反射，能量已经可以忽略不计了。

图3　用线性褶积模型计算的二道合成地震记录

图4　用非线性褶积模型计算的同样二道合成地震记录

图3和图4是由中原油田两个测井资料得到的反射系数序列和拟似反射系数序列计算的合成地震记录。比较它们可以发现两个明显的不同：一是非线性褶积模型计算的合成地震记录比线性褶积模型计算的合成地震记录要长很多，这是因为多次反射在一次反射结束后还继续存在；二是线性褶积模型计算的合成地震记录的较后部分比非线性褶积模型计算的合成地震记录相同时间部分的值要大，这是因为前者没有考虑透射影响的原因，而且在原始反射系数数值较大的记录（左边一道）上差异更为明显。当然，由于原始反射系数数值普遍都较小，故层间多次反射不太明显，所以两种合成记录的前面部分差异不大。如果使用减去法可以将所有的多次反射都提取出来，进而研究该区多次反射的特征^［3］。

2　二维射线追踪数值模拟边缘绕射波计算

射线追踪数值模拟只能模拟地震波场的运动学特征，但是运算速度快，而且提供的射线图

十分直观，在生产实际中还是很受青睐的。为了解决射线追踪数值模拟无法解决绕射波的问题，发展了边缘绕射理论。

根据边缘绕射理论，入射波到达介质性质发生突变的绕射点时，会向四面八方发出绕射波

射线。注意，与反射波射线代表反射波主能流方向不同，绕射波射线与绕射波主能流方向无关。绕射波射线能直观地给出绕射波的运动学特征。至于绕射波的动力学特征，Berryhill^［4］利用波动方程的克希霍夫积分解计算得到半无限水平界面边缘上的绕射波是地震子波与下述绕射响应算子的褶积。

其中：

而：

式中各量如图8所示，t为时间，它以震源子波起始时间为原点，t_a为绕射波走时，V为介质波速，a为半炮检距，z₀为炮检中点到半平面的法线距离，x₀为炮检中点在半平面（可能在其延展平面）上的投影到中断边缘线的距离,，R₀＝√x²₀＋z²₀为炮检中点到中断边缘线的距离，θ₀为半平面法线与炮检中点到中断边缘线的最短连线间的夹角，U（t）为单位跃变函数，即：

上述公式只对半无限水平界面边缘上的绕射波成立。为了能使用上述公式计算半无限倾斜界面边缘上的绕射波，需要根据波的物理性质引入等效炮点和等效检波点的概念。它们必须满足的条件为：①入射波射线不变；②绕射波射线不变；③等效炮检点的连线与半无限倾斜界面平行；④走时不变。据此，可以导出等效炮点S_e（x_se，z_se）和等效检波点G_e（x_ge，z_ge）应该满足下列联立方程组：

式中：（x_s，z_s）、（x_g，z_g）和（x_d，z_d）分别是真实炮点S、检波点G和绕射点D的坐标，而（n_x，n_z）则是绕射点上界面法线的方向余弦。求解上述方程组就可以得到等效炮点和等效检波点（图5），进而计算半无限倾斜界面边缘上的绕射波。

图5　等效炮点和等效检波点示意图

图6　尖灭结构的等效界面示意图

前面讨论的半无限界面比较简单，当地下存在如图6和图7所示的楔形结构（如尖灭、不整合等）时，绕射响应算子的计算就困难了（图8）。如图7所示，由于界面R₂和R₃以上的介质不是均匀的，故不能直接使用式（10）来计算它们的绕射响应算子。这个困难同样可以基于波的物理性质引入等效界面加以解决，此时其上部的界面视为不存在。等效界面必须满足的条件为：①绕射点位置不变；②入射波射线不变；③绕射波射线不变；④反射系数不变。如图6所示，利用统一的波的运动学边界公式，可以导出等效界面的单位法线向量为：

图7　角度不整合模型

图8　计算绕射响应算子时各量的示意图

式中：ξ₁和η₂分别为界面1的入射波和出射波的单位射线向量。η₂由下列诸式确定：其等效反射系数为：

式中：cosα₁＝ξ₁·n₁；cosα₂＝η₁·n₂；cosα₃＝ξ₂·n₁。注意，α₁、α₂和α₃必须要取0～之间的值。

再者，实际地层界面是不满足半无限平面这一条件的，但是许多物理现象表明：绕射波基本上只与绕射点附近的界面性质有关。绕射波的形成可以用边缘处第一菲涅尔带的缺损来解释，而射线法是在高频条件下的一种近似，高频条件下的第一菲涅尔带很小，其缺损只受边缘附近界面性质的影响。因此，可以用实际地层界面在绕射点处的切平面来近似计算。这就是弯曲界面边缘绕射波计算的等效平面。

至此，就在实际上完全解决了二维射线追踪中边缘绕射波计算的全部问题。下面给出部分计算的结果。图9是用包括边缘绕射波的射线追踪数值模拟方法计算出来的半无限水平界面（图8）上的自激自收剖面（a）和炮集记录（b）。图10是用射线追踪数值模拟方法计算出来的角度不整合模型（图7）上的自激自收剖面，其中（a）为没有考虑边缘绕射波的结果，（b）为考虑了边缘绕射波的结果。从这些记录剖面上可以清楚地看到绕射波的影响。

图9　半无限水平界面上自激自收剖面（a）和炮集记录（b）

图10　角度不整合模型上自激自收剖面

（a）未加绕射波；（b）加了绕射波

3　弹性波波动方程数值模拟与应用

弹性波波动方程数值模拟包括地质模型的数学建模和波动方程的数值计算两部分内容。一般来说，越是存在油气藏的地方，地下介质越是复杂，其数学建模越是困难。利用随机介质理论能方便地实现复杂非均匀介质的数学建模。波动方程的数值计算方法有很多种，目前实际工作中使用最多的是有限差分法。

3.1　非均匀随机介质模型数学建模^［8～16］

在复杂的非均匀介质中存在大量分布极不规则的异常。对于这些异常，用确定性方法描述它们十分困难，而使用统计学方法则可以只用少数几个简单的参数（如空间自相关函数、自相关长度和均值等）灵活、方便、完整地描述它们。用统计学方法描述介质的非均匀性所形成的非均匀介质模型就是随机介质模型。

各向同性弹性介质由其密度ρ和拉梅参数λ、μ所确定。在随机介质中，介质参数可以表示为：

式中：ρ₀、λ₀、μ₀为背景介质参数，假设为常数或随空间坐标（x，z）缓慢变化；δρ、δλ、δμ为加在上述背景上的非均匀扰动量，假设其为具有零均值、一定方差及某一自相关函数的空间平稳随机过程。

假定空间随机扰动σ＝σ（x，z）为二阶平稳随机过程，且均值为零，方差为ε²，协方差函数为C（x，z），从而自相关函数为：

可以根据实际情况选择自相关函数，并由此构造相应的随机介质模型。地震勘探中，随机介质模型的自相关函数通常可选择高斯型、指数型或冯·卡尔曼型的函数。它们分别具有不同的特点，适应不同的地质情况。

例如，使用指数型椭圆自相关函数：

式中：a、b分别是介质在x方向和z方向上的自相关长度。图11给出了指数型随机介质模型不同自相关长度所对应的随机介质（P波速度）模型，其中均值都为3 500m/s，标准差为ε＝0.1＝10%（这时其扰动范围基本上不超过20%）。从图中可以看出自相关长度a、b分别描述了随机介质在水平（x）方向及深度（z）方向上非均匀异常的平均尺度。图11所示的随机介质模型与地下取得的岩心图非常相似，这说明可以用指数型随机介质模型很好地模拟实际地质情况。当a＝∞时二维随机介质退化为一维随机介质模型［图1（f）］，它与陆相薄互层剖面图十分相似，这说明用该随机介质模型可以很好地描述陆相薄互层的非均匀性。

除了上述几种随机介质模型之外，又发展了混合型、分形型等多种随机介质模型，还发展了非平稳和多尺度等更复杂的随机介质模型。

图11　具有不同自相关长度的6个椭圆型自相关函数随机介质模型

3.2　弹性波波动方程数值模拟计算^［8～16］

在各向同性、完全弹性的非均匀介质中，弹性波波动方程为如下的方程组：

式中：τ_xx＝τ_xx（x，z，t）、τ_zz＝τ_zz（x，z，t）、τ_xz＝τ_xz（x，z，t）是应力张量；ρ＝ρ（x，z）是密度；［v_x＝v_x（x，z，t），v_z＝v_z（x，z，t）］是速度向量；λ＝λ（x，z）和μ＝μ（x，z）是拉梅系数。

用交错网格差分格式的有限差分法求解上述方程组。该差分格式可以精确、稳定地应用于任意复杂变化的非均匀介质模型。当一个波长中的网格点数多于10个时，网格色散与网格各向异性均可忽略不计，但在计算时均需要使用吸收边界条件。图12给出了对3个不同的介质模型、取两种不同频率的地震子波所进行的数值模拟结果：6个弹性波垂直分量自激自收记录。

图12　用不同频率地震子波对不同自相关长度随机介质模型数值模拟的自激自收记录

左边2图对应于均匀介质模型，其弹性参数为v_p＝3 500m/s，v_s＝2 121m/s，ρ＝2 600 kg/m³；中间2图对应于a＝b＝1的随机介质模型，右边2图对应于a＝b＝20的随机介质模型，它们的背景参数值与均匀介质相同，扰动标准差ε＝10%。上面3图是使用主频为150Hz（视波长为23.3m）的高频地震子波计算出来的，下面3图是使用频率为40Hz的低频地震子波计算出来的。

观察高频地震子波的模拟结果：图12a.2与图12a.1无明显的差异，说明当异常尺度（a＝b＝1）远小于地震波长（λ＝23.3m）时，随机异常不会十分明显地影响地震波的传播；而图12a.3与图12a.1有明显的差异，说明当异常尺度（a＝b＝20）接近地震波长（λ＝23.3m）时，随机异常会十分明显地影响地震波的传播。

观察低频地震子波的模拟结果：图12b.1～图12b.3的3个图基本相似。比较图12b.3与图12a.3可以看出，若异常尺度（a＝b＝20）不变而降低子波的频率（使地震波长远大于异常尺度），则随机异常对地震波传播的影响便相应地减小。

3.3　在碳酸盐岩孔洞型储层预测中的应用^{［17～21］}

为了回答能否利用常规频带（20～40Hz）地震资料直接检测我国西部地区埋藏较深（一般在5 000m以上）的下奥陶统碳酸盐岩储层中高度数米至数十米的岩溶洞穴的问题，也为了解决定量预测碳酸盐岩孔洞型储层的孔洞宽度问题，开展了系统的弹性波波动方程数值模拟研究。

设计了一系列高度相同、宽度不同的孔洞模型（一个模型上只有一个孔洞），进行弹性波波动方程数值模拟计算，对计算出来的各模型垂直分量地震记录进行动校、叠加等处理得到叠加剖面，从叠加剖面各叠加道上提取孔洞绕射波的最大振幅，制作出各模型（即各孔洞）绕射波振幅的变化曲线（图13中不同颜色的曲线就对应着不同洞宽的孔洞模型），将同一洞高、不同洞宽的孔洞模型绕射波振幅变化曲线放在一起就组成了不同洞宽的孔洞模型绕射波振幅对比图（图13中的每一个小图均是这样的对比图）。由图中的每一个小图都可以清楚地看出，随着孔洞宽度由小变大，绕射波的振幅由弱变强。

图13　高度相同、宽度不同孔洞的绕射波垂直分量振幅对比图

（a）洞高2m；（b）洞高3m；（c）洞高4m；（d）洞高6m；（e）洞高8m；（f）洞高10m；（g）洞高15m；（h）洞高30m

将与孔洞中心点相对应的那一道的绕射波振幅最大值（以宽度600m孔洞的振幅为1）相对于孔洞宽度画一个图，图14为高度10m孔洞的图。由图可见，开始时孔洞宽度变大振幅值迅速增大；到70m时增幅减慢，并逐渐趋于平稳；孔洞宽度继续增大时振幅趋于极限1。称宽度很大的孔洞为宽度无限的孔洞，它们实际上相当于薄层。

图14　高度＝10m、不同宽度孔洞绕射波垂直分量振幅变化图

一个宽度无限孔洞的顶、底面反射波会叠合在一起，形成复合反射波。反射波振幅与用下式计算的顶底面复合反射系数R成正比。

式中：h为孔洞高度，r为孔洞顶、底面反射系数，λ为孔洞中传播的地震波波长。

根据式（17）计算了频率为25Hz时不同高度孔洞的复合反射系数R（表1）。

表1　f＝25Hz时，不同高度孔洞的复合反射系数R

由表1可见，宽度无限的孔洞能够形成很强的复合反射波。在频率为25Hz时，10m高孔洞的复合反射波相当于碳酸盐岩与砂泥岩界面（反射系数为0.2）反射波振幅的4.5倍以上，几乎近似于绝对反射；即使1m高的孔洞（相当于1/60波长），其复合反射波的振幅也相当于灰岩与砂泥岩界面反射波的2倍以上。

当然，碳酸盐岩中一般的孔洞都是宽度有限的。为了计算宽度有限孔洞的绕射波振幅大小，以宽度无限孔洞的振幅为1，将图14中的高度固定、振幅随宽度变化的曲线作为宽度振幅因子曲线，利用式（17）求出宽度无限孔洞的复合反射系数，再乘上相应的宽度振幅因子值计算出宽度有限孔洞的复合反射系数，估算出其反射（此时应称其为绕射波）的强弱。例如，在频率为25Hz时，10m高宽度无限孔洞的复合反射系数为0.968，宽度18m孔洞的宽度振幅因子值为0.21，则其复合反射系数为0.203。以上计算均是在孔洞中充填流体的情况下得出。

对于可检测性的一种合理认识是：一个地质体所形成的地震波，当其振幅不低于地震剖面上它所出现时段的背景振幅时，就可以被识别，因而该地质体就可以用地震方法检测出来。在研究孔洞储层时，由于它的顶、底面反射系数远远大于0.2（超常反射系数），即使很小的洞也能形成很强的波；而其周围是准均质介质，所形成的是杂乱无章的幅度较低的散射波。因此，在5 000m以下的深层，在常规地震频带内可以从地震剖面上识别出孔洞形成的波，并对它进行研究。

根据这一认识，假设背景水平R_s＝0.2，可以得出充填流体的高度10m、宽度大于18m的孔洞即可被检测到，而孔洞中有其他充填物时的可检测极限宽度也可以使用同样的方法得到。

利用上面计算的宽度振幅因子曲线也可以定量预测碳酸盐岩孔洞型储层的孔洞宽度：当确定了地震剖面上某一绕射波为孔洞型储层所产生的之后，提取该绕射波的最大振幅并进行一定的归一化处理，然后根据处理后的最大振幅在宽度振幅因子曲线上的位置，就可以得到相应的孔洞宽度。

4　细胞自动机方法地震波数值模拟^［22］

长期以来，地震波数值模拟计算方法都是以由牛顿力学导出的波动方程为基础进行研究的（射线追踪表面上看没有直接对波动方程进行计算；实际上，射线追踪的基本方程——程函方程是波动方程的高频近似，故射线追踪的理论基础也是波动方程）。目前，勘探的对象越来越复杂，而在复杂介质条件下的波动方程推导十分困难，有些还根本不可能。因此，需要建立和发展更灵活、更容易适应复杂地质条件的新的数值模拟方法。

近年来，一些非常规的地震波数值模拟方法出现了，细胞自动机方法就是这样的一种新方法。它是一种在微观级上进行简单计算，在宏观级上模拟复杂物理现象的方法。由于细胞自动机是在微观领域中进行的计算，故具有潜在的模拟任意复杂介质的能力；而且它的计算十分简单（无算术运算），无舍入误差，无条件稳定，还非常适合于并行计算，极有发展前途。

用细胞自动机进行模拟，首先需要将地下研究区域划分为由一系列格点（细胞）所组成的网格，每个格点上的状态（数）为有限的；然后将时间离散化，进行时间演化。此时，每个格点每个时刻的状态由其相邻格点前一时刻的状态所决定，决定过程由一些简单的局部相互作用规则所控制。如果规则中考虑了物理过程的微观机制，就能模拟复杂的物理过程。

由此可知，与有限差分法、有限元法等那些只作空间和时间离散的方法不同，细胞自动机方法不仅离散了时间和空间，而且也离散了状态变量本身。这些离散的状态变量的变化不是由宏观的微分方程所决定的，而是由微观空间的局部相互作用规则所控制的。

细胞自动机的另一特点是其统计平均特性的利用。它的输出是对大量离散的状态变量进行宏观统计平均后得到的。虽然一个细胞自动机的微观量可能是离散的、近乎随机的，但进行统计平均后的宏观量则是连续的、确定性的。这种宏观与微观的统一正是现代统计物理学的基本思想。因此，细胞自动机是一种基于统计物理学的新方法。

用于地震波场数值模拟时，地下研究区域首先需要划分为各种类型的网格，例如长方形网格或三角形网格。如图15所示，先将地下研究区域划分为三角形网格，具有单位速度、单位质量的粒子分布在格点上，在三角形网格中，每个粒子有6个可能的移动方向；当然，在演化过程中，每个粒子只能沿6个可能方向的一个方向移动。时间演化的局部相互作用规则十分简单：每经过一个时间步长，粒子沿规定的方向前进到下一个格点；如果在该格点上没有碰到其他粒子，则下一个时间步长继续沿原方向运动且速度不变，再经过一个时间步长到达再下一个格点；若碰到了其他粒子，则按“碰撞”规则改变运动方向。“碰撞”规则与网格形状有关。例如，在三角形网格中的几个“碰撞”规则如图16所示：如果两个粒子正面碰撞，则每个粒子都要改变60°方向运动，向顺时针方向改变60°和向逆时针方向改变60°的几率是相同的，都是1/2［图16（a）］；如果3个粒子碰撞，则每个粒子都原路返回［图16（b）］。当然，由于三角形网格中粒子碰撞的可能情况比较多，仅仅这两条规则是不够的，其他规则在此就不一一述说了。

图15　三角形网格

图16　“碰撞”规则示意图

模拟开始之前，首先将随机分布的大量背景粒子安置在网格上，每个格点上分布的粒子数和每个粒子的运行方向都是充分随机的。另外，在一个很小的“源”区内模拟波源，即全部安置最多数量（6个）的粒子。然后，就可以按细胞自动机的时间演化方式自动进行演化了。经一定时间的演化后，进行统计平均即可以得到地下介质中传播的地震波场或地面记录到的地震记录。图17就是通过细胞自动机方法数值模拟出来的地震记录和地下波场的快照。由于试验中为了节省计算时间，划分的网格较大，统计所使用的格点较少，所以模拟出来的地震记录和波场快照显得比较粗糙，似乎存在许多计算引起的干扰，这些问题完全可以通过缩小网格、增多格点数的方式来解决。

图17　利用细胞自动机数值模拟出来的地震记录（a）和波场快照（b）

由细胞自动机进一步发展了格子气自动机、Bolzmann机等更复杂的数值模拟方法，可以适应更为复杂的介质。

5　结论与认识

地震正演数值模拟在地震波场的理论研究和地震勘探、开发的实际中都具有重要的意义，也得到了广泛的应用。长期以来，地震正演数值模拟得到了长足的发展，考虑的对象越来越复杂，研究的问题越来越深入。我们介绍的发展包括从线性到非线性，从只考虑几何地震学到考虑物理地震学，从声波模拟到弹性波模拟，以及数学建模和实际应用等工作。从模拟和应用的结果来看，确实达到了改进和深入的目的。但是，目前地震正演数值模拟的发展仍然不能满足地震理论研究和实际工作的需要，还需要建立和发展更灵活、更容易适应复杂地质条件的新的数值模拟方法。细胞自动机、格子气自动机和Bolzmann机等代表了一类新的数值模拟方法，具有较强的生命力。

参考文献

［1］许云.地震层状构造理论与应用［M］.北京：石油工业出版社，1981.

［2］姚姚.层状介质非线性模型反演问题［C］.地震勘探学术交流会论文集.南京：江苏科学技术出版社，1981：54-57.

［3］姚姚，蔡其新，张果，董桥梁.微屈多次波正演模拟及其在中原油田深层地震资料解释中的应用［J］.物探与化探，2004，28（5）：443-446.

［4］Berryhill J R.Diffraction response for nonzero separation of source and receiver［J］.Geophysics，1977，42（5）：1 158-1 176.

［5］黄延诂，姚姚.二维层状介质中绕射波的计算［C］.岩性物探学术讨论会论文集.成都：四川科技出版社，1988：34-37.

［6］黄延诂，姚姚.二维层状介质中衍射波的计算［J］.地球物理学报，1989，32（1）：76-88.

［7］Huang Yanhu，Yao Yao.Calculation of Diffracted Waves in 2-D layered media［J］.ACTA Geophysics Sinica，1989，32（2）：313-329.

［8］奚先，姚姚.随机介质模型理论概述［J］.地学探索，2000，2（4）：125-128.

［9］姚姚，奚先.浅层地震随机介质理论及发展［M］.工程地球物理学进展.武汉：武汉水利电力大学出版社，2000，52-60..

［10］奚先，姚姚.二维随机介质及波动方程正演模拟［J］.石油地球物理勘探，2001，36（5）：546-552.

［11］姚姚.随机介质模型正演模拟［M］.中国东部地区物探技术研讨会论文集［J］，2001：167-174.

［12］姚姚，奚先.随机介质模型正演模拟及其地震波场分析［J］.石油物探，2002，41（1）：31-36.

［13］奚先，姚姚.随机介质模型的模拟与混合型随机介质［J］.地球科学——中国地质大学学报，2002，27（1）：67-71.

［14］姚姚，奚先.区域多尺度随机介质模型及其波场分析［J］.石油物探，2004，43（1）：1-7.

［15］奚先，姚姚.非平稳随机介质模型［J］.石油地球物理勘探，2005，40（1）：71-75.

［16］Yao Yao，Xi Xian.Review on reflected seismic wave in random media［J］.Journal of Geophysics and Engineering，2004，1（2）：147-152.

［17］姚姚.非平稳随机介质模型与洞缝型油气藏弹性波正演模拟［J］.石油地球物理勘探，2002，37（专刊）：168-171.

［18］姚姚，唐文榜.深层碳酸盐岩岩溶风化壳洞缝型油气藏可检测性理论研究［J］.石油地球物理勘探，2003，38（6）：623-629.

［19］Yao Yao，Sa Liming，Wang Shangxu.Research on the Seismic Wave Field of Karst Cavern Reservoirs near Deep Carbonate Weathered Crusts［J］.Applied Geophysics，2005，2（2）：947-102.

［20］姚姚.菲涅耳带与洞缝型油气藏地震波场［J］.石油物探，2005，44（5）：491-494.

［21］吴俊峰，姚姚.碳酸盐岩特殊孔洞型构造地震响应特征分析［J］.石油地球物理勘探，2007，42（2）：180-185.

［22］姚姚.格子Boltzmann方法地震波场模拟［J］.地球科学——中国地质大学学报，1997，22（6）：54-60.