基于二维地质建模的两种地震数值模拟方法的应用及分析

赵忠泉

（广州海洋地质调查局 广州 510760）

作者简介：赵忠泉，男，（1983—），硕士，主要从事海洋油气资源调查研究工作，E-mail：zzqhello＠163com。

摘要 利用地震数值模拟技术结合实际资料，可以建立各种地质体的地震识别模型，有效地避免地震现象的多解性，从而可以提高解释的精度。本文介绍了二维地质建模的方法流程及两种模拟方法-褶积法和PSPI波动方程法，前者无边界条件约束和频率域中的信号损失，简洁易行，计算稳定，应用广泛，是最早的地震波场模拟方法；后者通过求解波动方程，包含丰富的波场信息，能够充分反映地震波的动力学和运动学特征。实际应用中利用褶积法对三维潮道模型及简化的碳酸盐岩多旋回倾斜薄互层沉积模型进行了模拟；利用零炮检距的频率波数域的波动方程法模拟了生物礁的地震响应，结果对于碳酸盐岩生物礁识别有一定指导意义。

关键词 地质建模 数值模拟 褶积法PSPI法

不同地质体由于其岩性、物性、含油气性、内部结构和岩石组合等的差异，在地震上具有不同的反射特征，包括内部结构、外部形态、振幅、频率等参数。由于地震波在地下地质体中传播的复杂性，加上各种干扰，造成了地震剖面中的各种反射现象存在多解性，大大增加了地震解释的难度。利用地震数值模拟技术结合实际资料，在建立不同地质体的地震识别模型的同时也有效地避免地震现象的多解性，从而可以提高解释的精度。

1 地质建模

地震数值模拟技术的基础是地质地球物理模型的建立，可归结为对地质及地球物理模型结构的数学描述。

二维封闭结构模型用于建立复杂地质模型。二维封闭结构模型就是定义相同地质属性为一独立封闭的地质单元，按照地质属性将地质模型划分成多个独立封闭的地质单元，把所有独立封闭地质单元按照空间分别有序地排列起来，这样组成的集合体就构建了一个二维地质模型。封闭结构模型是以积木方式定义地下地质结构，可以描述非常复杂的地质体。二维封闭结构模型被描述为具有相同地质属性（速度、密度等），并被地层界面、断层界面或模型边界所围成的地质单元的有机组合。对封闭结构模型的描述，实际上就是描述封闭地质单元和封闭地质单元之间的关系，前者包括对封闭地质单元属性和封闭地质单元边界的描述；后者是对地质单元空间关系的描述，也就是描述封闭地质单元边界相接关系及地层属性［1］。

在进行数值模拟过程中，为了验证某些复杂地质体的波场特征，需要绘制多种不同的地质模型，通常可借助常规绘图软件（绘图板、Photoshop、CorelDraw，AutoCAD等）绘制好二维封闭结构面，再根据图像处理中的区域填充算法（种子填充和扫描转换填充），对不同二维封闭结构面进行不同颜色的填充。其中不同颜色代表不同的二维封闭结构面属性（速度、密度等）；合并相同属性的封闭面，形成最终的二维封闭结构模型［1］。为了得到二维封闭结构模型的属性（速度、密度等）模型，需要对二维封闭结构模型的彩色图进行速度像素空间和属性空间转换，根据颜色空间和属性空间的相互映射，就可以得到复杂地质体的属性（速度、密度等）模型，如图1为模型创建流程图。

图1 二维封闭结构模型建立流程图

2 两种数值模拟方法

21 褶积模型

在褶积模型中，我们把地震反射信号s（t）看作是地震子波w（t）与地下反射率r（t）的褶积。地震子波w（t），使用实际地震系统记录到的地下一个单独的反射界面反射的波形（如图2，理想的无噪声褶积过程）。反射率r（t）则代表理想的无噪声地震记录。记录到的地震道s（t）可看作是地震信号w（t） r（t）与可加噪声n（t）之和，因此可以把地震道看作是一种有噪声干扰的，经过了滤波的地下反射率的变形。

在无噪声褶积模型中，我们把地震信号S（t）看作是地震子波w（t）和地下反射系数r（t）的褶积：

南海地质研究2012

式中：s（t）——合成地震记录；

r（t）——反射系数；

w（t）——地震子波。

图2 褶积过程

22 PSPI波动方程法

通过求解波动方程的数值模拟方法，能够充分反映地震波的动力学和运动学特征，波场信息丰富，模拟结果较为准确。这里仅介绍适合横向速度剧烈变化的频率-波数域相移加插值的波场延拓方法［2］。

相位移加插值的波场延拓方法，简称PSPI法，基本思想是在波场向下延拓的每个深度步长Δz之内，将波场的延拓分成两部进行，首先用L个参考速度V1，V2，…VL，将位于深度zi处的波场p（x，zi，ω）延拓到zi＋1＝zi＋Δz处，得到L个参考波场p1（x，zi＋1，ω），p2（x，zi＋1，ω），…，PL（x，zi＋1，ω）。第二步，按实际的偏移速度V（x，z）同参考速度V1，V2…，VL的关系，用波场插值的方法求出zi＋1处的波场p（x，zi＋1，ω），按同样的步骤，可将zi＋1处的波场值p（x，zi＋1，ω）延拓到深度zi＋2，得p（x，zi＋2，ω），直到延拓到最大的深度zmax为止。

对于各向同性介质，取二维标量声波方程作为延拓的基本方程：

南海地质研究2012

式中，p＝p（x，z，t）为二维地震波场值；x，z分别为水平方向和垂直方向坐标轴；t为时间轴；v（x，z）为纵、横向都可变的地震波传播速度。将式（2）分别对x、t作傅氏变换，考虑到并考虑到ə2/əx2和与（-ikx）2和（iw）2的对应关系，可得：

南海地质研究2012

式中， 是p（x，z，t）的二维傅氏变换；v为地震波速度；w为圆频率；kx为水平波数；kz为垂直波数。零炮检距情况下的地震记录模拟只考虑单程波，因此可得到相位移波场延拓公式如下：

南海地质研究2012

式中， （kx，zi，w）为频率波数域波场值；Δz为深度延拓步长；kx为测线方向波数；kz为深度方向波数。式（4）为二维波场正演公式，其延拓方向为由地下向地面延拓；式（5）为二维波场偏移公式，其延拓方向为由地面向地下延拓。

为了适应地下地震波场速度在纵横向均可变的要求，在同一延拓深度内用几个不同地震波速度分别作相移，再用拉格朗日插值公式进行插值，就可求出所有的以不同速度传播的延拓波场值P（x，zi＋1，t），从而近似地解决了横向变速时的波场延拓问题［3］。

3 模拟实例

31 三维潮道数值模拟

运用褶积原理建立了一个简单三维潮道模型，此三维潮道事实上为多个（128）二维剖面排列而成，三维模型的采样点为128×128×128，利用MATLAB实现。选用子波为雷克（Ricker）子波，其公式为：

南海地质研究2012

其中fp为主频。在处理过程中选用主频为fp＝40 Hz、采样间隔2 ms，对称采样点数为24，子波波形如图3。

图3 雷克子波

图4 潮道平面图

图4为潮道平面图，该图仅反映了潮道的平面形态，作为计算机实现三维建模的边界控制，横坐标代表inline线，纵坐标代表xline（crossline）线，图5为三维地质模型示意图，模型较简单，整体由三个水平层叠置而成，在第二层和第三层之间镶嵌了形如图4的潮道，此潮道没有考虑进水方向，根据此地质模型进行计算机地震正演模拟，可得到相应三维地震数据体，从图中可以看到，**虚线（上）和蓝色虚线（下）位置上，分别横跨了三个潮道分支和两个潮道分支，就是说在相应两条虚线位置上的两条测线应该分别有三个和两个潮道显示，提取相应的两条剖面如下图6和图7：

图5 三维地质模型

图6 xline＝100（**虚线）剖面

图7 xline＝100（蓝色虚线）剖面

再在三维数据体中沿水平方向做切片，即提取时间切片。图8为时间切片在地震剖面上的位置示意图，图中五条标示线从上到下依次为白色实线、**虚线、白色实线、红色虚线和白色实线，与之对应的时间分别为70 ms、85 ms、95 ms、99 ms和110 ms（时间范围是0～128 ms），图9～图13为相应切片，从图中可以看出，随着所做切片时间的增大（深度的增加），潮道的展布范围逐渐减小，由于地层是水平层状的，使得时间切片等同于地层切片和沿层切片，其切片效果非常明显，切片中潮道形态得到了很好的展示，但是在多个切片中发现，从可以见到潮道形态一直到潮道消失的时间范围是在70～110 ms之间，而潮道的真实范围是在80～100 ms之间，显然依据切片所圈定的潮道的范围相比真实的范围扩大了，究其原因是由于不管选取哪一种子波，子波都有一定的延续长度和有限频宽，这就限制了合成地震记录本身的分辨率并不能达到等时厚度反射系数序列的分辨率。因此在对实际地震资料进行解释的时候，对地质异常体边界的识别应该考虑地震子波并非脉冲波所带来的影响。

图8 剖面示意图

图9 切片t＝70 ms

图10 切片t＝85 ms

图11 切片t＝95 ms

图12 切片t＝99 ms

图13 切片t＝110 ms

32 薄互层沉积模型

图14为简化的碳酸盐岩多旋回倾斜薄互层沉积模型（Zeng，2003），模型简化是为了更好地突出由岩相控制的波阻抗结构和地震信号之间的相互关系。该模型所有倾斜的倾角都相同，每层都有相同的垂直时间厚度（5 ms或15 m，速度为6000 m/s），泥岩与低孔隙度颗粒灰岩的波阻抗差，以及低孔隙度颗粒灰岩与高孔隙度颗粒灰岩的波阻抗差都相同，所有高孔隙度颗粒灰岩具有相同的深度范围，综合起来形成了一个水平的岩性地层单元。

其时间域地震响应（图15）中，高频情况下（60 Hz雷克子波），地震反射被建设性地调谐到时间地层单元，因此地震同相轴沿着时间地层单元分布（图15a）。当子波频率减到40 Hz时，地震反射对时间地层单元和岩性地层单元都有响应（图15b）。当用30 Hz雷克子波时（图15c），地震同相轴破坏性地调谐到时间地层单元和建设性地调谐到岩性地层单元，因而时间地层单元的反射进一步变弱，地震同相轴被岩相反射所控制［4］。

这个模拟过程强调了了解地质格架和时间地层单元以及岩性地层相带厚度尺度的重要性。时间地层（图15a）和岩性地层（图15c）成像都是有用的，前者用于对比，后者用于粗略的储层评价。然而，这两种响应不能混淆在一起。图15b中的两组相互矛盾的地震同相轴会造成地震假象［4］。

图14 简化的碳酸盐岩多旋回倾斜薄互层沉积模型

33 生物礁数值模拟［5～7］

频率—波数域的相移加插值偏移（PSPI）在每一个深度间隔内使用多个参考速度进行偏移，由多个偏移结果插值生成最终的偏移剖面，所用插值的速度越多，越能反映实际介质的速度变化情况，此方法在成像精度及横向变速适应性上具有很大的优越性，但处理所需的时间稍长，鉴于本文的二维叠后建模对处理时间没有过高要求，因此应用PSPI方法做正演、偏移。

图16为某区块过生物礁的原始地震剖面，图17为根据此剖面建立的生物礁速度模型：模型速度变化范围是5600 m/s到5980 m/s，从图16中可以看出生物礁的底界面清晰可辨，围岩有披覆现象，内部呈杂乱反射。为了检验该地质建模的正确性，先采用PSPI方法对该模型进行了波场正演模拟计算，其模拟剖面如图18所示。由于生物礁埋藏深，生物礁顶底反射的弧度较大，不规则点的绕射波杂乱，因此用图15的速度模型对其进行叠后时间偏移，得到了偏移剖面（图19），横向表示256个地震道，纵向表示零偏移距反射时间，礁体最大时间厚度约40 ms。从图19可以看出，模拟记录中的礁体顶界与原始剖面有一定差距，但是生物礁底界反射和内幕反射以及侧翼反射与原始剖面基本一致，其他的地层界面形态与原始剖面也吻合较好，在一定程度上验证了地质模型的正确性，说明当生物礁与围岩之间存在一定波阻抗差异时，在地震剖面上必然出现异常反映，经过有效的构造和参数反演，能够将其分辨出来。相信通过模型改进以及算法中参数的调整，能够与原始剖面更好地吻合，从而为生物礁的地震解释提供一种有力的验证工具。

图15 图14模型时间域地震响应

4 结论

地震数值模拟（正演）技术基于地球物理模型的建立，运用概念二维封闭结构地质模型的建立方法，得到复杂地质体的数学模型，结合各种算法对其进行模拟从而可以验证相应地质体的地震波场特征；结合实际资料建立不同地质体的地震识别模型，可以有效地减少地震现象的多解性，从而提高解释的精度；褶积法无边界条件约束和频率域中的信号损失，简洁易行，计算稳定，应用广泛，本文用此方法模拟的伪三维潮道模型及倾斜薄互层模型取得了较好的效果；通过求解波动方程的数值模拟方法，包含丰富的波场信息，能够充分反映地震波的动力学和运动学特征，PSPI波场沿拓方法为其中之一，利用正演与偏移相结合的流程模拟了生物礁的地震响应特征，检验解释成果的正确性，为生物礁的地震解释提供了一种有力的检验工具。

图16 原始剖面

图17 生物礁地质速度模型（256×256）

图18 正演记录（子波主频30Hz）

图19 偏移剖面（子波主频30Hz）

参考文献

［1］刘远志碳酸盐岩地震相分析与数值模拟［D］成都：成都理工大学，2009

［2］韩建彦复杂地质体地震正演与偏移［D］成都：成都理工大学，2008

［3］贺振华，王才经等反射地震资料偏移处理与反演方法［M］重庆：重庆大学出版社，1989

［4］Zeng Hongliu ＆Kerans，CSeismic frequency control on carbonate seismic stratigraphy；a case study of the Kingdom Abosequence，West Texas，American Association of Petroleum Geologists Bulletin，200387，273～293

［5］贺振华，黄德济，文晓涛，等碳酸盐岩礁滩储层多尺度高精度地震识别技术［R］成都：成都理工大学地球探测与信息技术教育部重点实验室，2009

［6］熊晓军，贺振华，黄德济生物礁地震响应特征的数值模拟［J］石油学报，2009，30（1）：7～65

［7］熊忠，贺振华，黄德济生物礁储层的地震数值模拟与响应特征分析［J］石油天然气学报，2008，30（1）：75～78

The application and analysis of two kinds of seismic numerical simulation method based on the2D-geological modeling

Zhao Zhongquan

（Guangzhou Marine Geological Survey，Guangzhou，510760）

Abstract：Pick to using seismic numerical simulation technology combined with the actual seismicdata，we can build all kinds of seismic recognition model of geologic body and effectively avoidthe multiple solutions of seismic phenomenon，which can improve the precision of the explana-tionThis paper describes the method of the process of 2D geological modeling and two simulationmethods，seismic convolution method and PSPI wave equation method，the former has no bounda-ry condition and the signal loss in frequency domain，is concise and easy，it can be calculatedsteadily and be applied widely，is the earliest simulation method in seismic wave field，the latterbased on the wave equation，it contains the rich information in wave field，can fully reflect thedynamics and kinematics characteristics of seismic waveIn the practical application，we use theconvolution model in 3D-tidal channel model and the multi-cyclic simplified deposition model oftilt thin interbed layer of carbonate；We simulate the seismic response of reefs using the method ofzero-offset wave equation in frequency and wave number domain，it is confirmed that the resulthas definite significance for the identify of the reef

Key words：Geological modeling Numerical simulation Convolution PSPI method

这是要借书看的，比如：（第一本也是讲matlab的）

《数学建模及其基础知识详解》 王文波 武汉大学出版社

《MATLAB在数学建模中的应用》 卓金武 北京航空航天大学出版社

Matlab软件及教程百度网盘免费下载

链接：>pwd=ffh6 提取码：ffh6

MATLAB是世氏美国MathWorks公司出品的商业数学软件，用于数据分友闹析、无线通信、深度学习、图像处理与计算机视觉、信号处理、量化金融与风险管理、机器人，控制系统等领域。链接包含各版本Matlab软件及相关基础和进阶视频教程及资料，涉及统计，信搜告散号处理，图像处理、量化等方向。

根据ode45()函数对微分方程的形式要求, y=ill(t,x)实际上是 y'=ill(t,x), y必须是列向量

y(1) = dy1/dt, y(2)=dy2/dt

故这一行就是 dy1/dt = a y1y2-by1, dy2/dt= -ay1y2

正是sir模型方程(y1=i, y2=s)

08年建模A题吧，用我的程序，以相图中椭圆最长弦为圆心，程序如下：

I = imread('clip_image','jpg'); % 读入

J = edge(I,'prewitt'); % 边缘提取

K=bwlabel(J,8); %连同区域分区

%共5个园，可以分成五个区，各○边缘点数量为Xi

X1=find(K==1);

X2=find(K==2);

X3=find(K==3);

X4=find(K==4);

X5=find(K==5);

%各○边缘点坐标

%第2个点

for i=1:length(X1)

x1(i)=ceil(X1(i)/281); %横坐标

y1(i)=X1(i)-(x1(i)-1)281; %纵坐标

end

R=0;

for i=1:length(X1)

for j=i+1:length(X1)

r=sqrt((x1(i)-x1(j))^2+(y1(i)-y1(j))^2);

if r>R

R=r;m=i;n=j; %记录坐标点

end

end

end

d1=[(x1(m)+x1(n))/2;(y1(m)+y1(n))/2] ; %%中心坐标

%第2个点

for i=1:length(X2)

x2(i)=ceil(X2(i)/281); %横坐标

y2(i)=X2(i)-(x2(i)-1)281; %纵坐标

end

R=0;

for i=1:length(X2)

for j=i+1:length(X2)

r=sqrt((x2(i)-x2(j))^2+(y2(i)-y2(j))^2);

if r>R

R=r;m=i;n=j; %记录坐标点

end

end

end

d2=[(x2(m)+x2(n))/2;(y2(m)+y2(n))/2] ; %%中心坐标

%第3个点

for i=1:length(X3)

x3(i)=ceil(X3(i)/281); %横坐标

y3(i)=X3(i)-(x3(i)-1)281; %纵坐标

end

R=0;

for i=1:length(X3)

for j=i+1:length(X3)

r=sqrt((x3(i)-x3(j))^2+(y3(i)-y3(j))^2);

if r>R

R=r;m=i;n=j; %记录坐标点

end

end

end

d3=[(x3(m)+x3(n))/2;(y3(m)+y3(n))/2] ; %%中心坐标

%第4个点

for i=1:length(X4)

x4(i)=ceil(X4(i)/281); %横坐标

y4(i)=X4(i)-(x4(i)-1)281; %纵坐标

end

R=0;

for i=1:length(X4)

for j=i+1:length(X4)

r=sqrt((x4(i)-x4(j))^2+(y4(i)-y4(j))^2);

if r>R

R=r;m=i;n=j; %记录坐标点

end

end

end

R4=R;

d4=[(x4(m)+x4(n))/2;(y4(m)+y4(n))/2] ; %%中心坐标

%第5个点

for i=1:length(X5)

x5(i)=ceil(X5(i)/281); %横坐标

y5(i)=X5(i)-(x5(i)-1)281; %纵坐标

end

R=0;

for i=1:length(X5)

for j=i+1:length(X5)

r=sqrt((x5(i)-x5(j))^2+(y5(i)-y5(j))^2);

if r>R

R=r;m=i;n=j; %记录坐标点

end

end

end

R5=R;

d5=[(x5(m)+x5(n))/2;(y5(m)+y5(n))/2] ; %%中心坐标

Di=[d1,d2,d3,d4,d5] %坐标中心在左上角

以上就是关于基于二维地质建模的两种地震数值模拟方法的应用及分析全部的内容，包括:基于二维地质建模的两种地震数值模拟方法的应用及分析、如何用MATLAB进行数学建模、数学建模软件MATLAB教程在哪能找到要基础的等相关内容解答，如果想了解更多相关内容，可以关注我们，你们的支持是我们更新的动力！