论文部分内容阅读
摘要：本文应用有限差分法求解声波波动方程，推导出时间域二阶空间域四阶的递推公式，同时对波场数值模拟过程中初始条件的给予、震源子波的选择、边界吸收的实现、稳定性的保证及频散的压制等问题进行了探讨和研究，为应用有限差分法进行波场数值模拟提供理论指导。
关键词：有限差分 地震波场 数值模拟 波动方程 时间切片
中图分类号：TD353.5 文献标识码：A 文章编号：1009-914X（2013）29-583-02
1引言
地震波场数值模拟是地震学领域中一项非常重要的技术方法。它是对特定地质、地球物理问题作适当的抽象，从而形成一个简化的数学模型，采用数值计算方法获取地震响应的过程，其目标就是假定在地下地质体特征已知的情况下，通过计算得到理论时深关系和时间切片图。开展地震波场数值模拟技术研究对于能够提高我们对地震波在地下介质中传播规律的认识，解决实际勘探工作中所面临的各种棘手问题[1]。地震波场数值模拟方法不但在石油、天然气、矿产资源、工程和环境地球物理等地震勘探中得到广泛应用，而且在地震灾害预测、地震区带划分以及地壳构造和地球内部结构研究等天然地震领域也得到普遍应用。
2有限差分原理
有限差分法就是把需要求解的地下地质特征体简化为一个关于速度和密度的模型，将此模型所在的区域划分为有限个差分网格，以离散的空间位置和离散时刻的值代替连续的物理量，以有限差分代替无限微分，以差分方程代替微分方程及其边界条件，以差商代替微商，利用微商与差商之间的近似关系将描述介质传播的微分波动方程转化为差分方程进行求解。
波动方程的差分离散方法有两种，第一种方法是将单变量的二阶波动方程直接转化为时间空间的二阶中心差分进行离散求解；第二种方法是将用位移表示的二阶波动方程转化为用应力和速度表示的一阶方程组，然后对此一阶方程组进行离散求解[47]。
在应用有限差分法的时候要注意选择合适的差分格式将微分方程转化为差分方程和微分方程数值求解二个方面。
3声波方程
在没有弹性横波只有弹性纵波存在时，对上式二边取散度：
利用旋度与散度的关系： ， ，交换上式的微分次序并化简可得：
其中， 表示纵波速度，其表达式为：
用u表示某一时刻t二维空间任意点（x，z）处的位移，则声波方程可以用u表示为：
为雷克子波的时域波形及其振幅谱，因为雷克子波的能量集中在正波峰处，与反射系数褶积后能较好的反应反射问题。
图1 雷克子波（主频25HZ）
用有限差分方式计算所得的地震道，我们看到，在z=150m处，速度由1200m/s变为2000m/s，密度由1.5g/cm3变为2.0g/cm3，反射系数发生较大变化，所以在0.25s处产生了较大的振幅波动；在z=300m处，速度由2000m/s变为1600m/s，密度不变，同样产生波动，但反射系数的变化不如z=150m处，所以振幅变化也小，同时由于速度变小，反射系数应与z=150m处相反，对应的就是在0.4s处的振幅波动；在z=450m，速度由1600m/s变为2500m/s，密度由2.0g/cm3变为2.6g/cm3，反射系数变化较大，也产生了较大的波动，对应的为0.9s处的振幅波动。理论模型和实际情况吻合，计算结果得到了验证。
4程序实现
（1）均匀介质倾斜地层模型
含有倾斜界面的四层地质模型，整个模型横向x=512m，纵向z=512m。当激发点在模型顶界面中心点处（256，0），接收线在z=0处，按照dx=1m，dz=1m进行网格离散化，震源项所用雷克子波的主频为35Hz，记录时间0.6s。
（2）建模方法
为了能将模型离散化，编制了MAKEMODEL程序，可以根据输入的控制节点建立速度和密度模型。
在编程过程中，考虑了弯曲界面的情况，插值方法采用线性插值法和最小曲率法2种方法，以便在不同情况下选择使用。
线性插值是数学、计算机图形学等领域广泛使用的一种简单插值方法，这种方法是通过在数据点之间连线以建立离散化模型的，模型采用直线连接的方式，在表达水平地层是比较合适，但对于弯曲界面等情况则无法使用。
最小曲率法是一种弯曲界面建模时的好方法，其试图在尽可能严格地尊重数据的同时，生成尽可能圆滑的曲面，广泛用于地球科学。用最小曲率法生成的插值面类似于一个通过各个数据值的，具有最小弯曲量的长条形薄弹性片。使用最小曲率法时要涉及到两个参数：最大残差参数和最大循环次数参数，以用来控制最小曲率的收敛。
（3）程序试算
水平四层模型，整个模型横向x=512m，纵向z=512m，分界面分别在z=100m，z=250m，z=400m处；介质速度由上到下分别为：1200m/s，1800 m/s，2200 m/s，2800 m/s；介质密度由上到下分别为：1.6g/m3，2.0 g/m3，2.2 g/m3，2.6 g/m3，具体如所示。当激发点在模型顶界面中心点处（256，0），接收线在z=0处，按照dx=1m，dz=1m进行网格离散化，震源项所用雷克子波的主频为35Hz，记录时間0.6s。
根据理论计算，三层水平大地反射波的时距曲线应该为双曲线，双曲线的最小值对应炮点所在位置，双曲线的曲率和速度有关，速度越大则曲率越小。T1层反射波所形成双曲线的最小值应在：
T2层反射波所形成双曲线的最小值应在：
T3层反射双曲线的最小值应在：
5结论
本文围绕地震波场数值模拟及波场特征分析等问题，在总结分析前人研究成果的基础上，以波动理论为基础，对有限差分法的一些关键性问题进行了探讨和研究，实现了地震波在几种典型地质模型中传播的数值模拟，并对模拟所得的单炮记录和波场快照进行特征分析。研究结果如下：
1.有限差分法在地震波场数值模拟技术中适用性强、算法容易实现、计算速度快、模拟结果与理论吻合良好、所得到得波场信息丰富，拥有广泛的应用前景；
2.二维有限差分法波场数值模拟作为一个好的辅助工具，可以很方便的计算地质体空间位置和形态大小变化对地震记录的影响，得到对实际工作有指导意义的结论；
3.对大倾斜、小断层、起伏界面、溶腔体等模型进行了计算，得到了明确的地震波场的特征，对模拟结果分析可知：在大倾斜地区，应当根据目的层的空间位置和倾角大小设置合理的观测系统；断层面产生的绕射波对于了解地下断层大致结构和特征有重要意义；凸界面的扩散波和凹界面的回转波是识别弯曲界面的特征波，在解决弯曲界面问题时加加以重视；
参考文献：
[1] 牟永光， 裴正林. 三维复杂介质地震数值模拟北京： 石油工业出版社， 2005. 215
[2] 张钋， 刘洪. 射线追踪方法的发展现状. 地球物理学进展， 2000， 15（1）： 36～45
[3] 刘福田等译. 地震学中的射线方法北京： 地质出版社， 1986.
作者简介：
董丽娜（1983-），女，助理工程师，2007年毕业于黑龙江科技大学，获工学学士学位。
关键词：有限差分 地震波场 数值模拟 波动方程 时间切片
中图分类号：TD353.5 文献标识码：A 文章编号：1009-914X（2013）29-583-02
1引言
地震波场数值模拟是地震学领域中一项非常重要的技术方法。它是对特定地质、地球物理问题作适当的抽象，从而形成一个简化的数学模型，采用数值计算方法获取地震响应的过程，其目标就是假定在地下地质体特征已知的情况下，通过计算得到理论时深关系和时间切片图。开展地震波场数值模拟技术研究对于能够提高我们对地震波在地下介质中传播规律的认识，解决实际勘探工作中所面临的各种棘手问题[1]。地震波场数值模拟方法不但在石油、天然气、矿产资源、工程和环境地球物理等地震勘探中得到广泛应用，而且在地震灾害预测、地震区带划分以及地壳构造和地球内部结构研究等天然地震领域也得到普遍应用。
2有限差分原理
有限差分法就是把需要求解的地下地质特征体简化为一个关于速度和密度的模型，将此模型所在的区域划分为有限个差分网格，以离散的空间位置和离散时刻的值代替连续的物理量，以有限差分代替无限微分，以差分方程代替微分方程及其边界条件，以差商代替微商，利用微商与差商之间的近似关系将描述介质传播的微分波动方程转化为差分方程进行求解。
波动方程的差分离散方法有两种，第一种方法是将单变量的二阶波动方程直接转化为时间空间的二阶中心差分进行离散求解；第二种方法是将用位移表示的二阶波动方程转化为用应力和速度表示的一阶方程组，然后对此一阶方程组进行离散求解[47]。
在应用有限差分法的时候要注意选择合适的差分格式将微分方程转化为差分方程和微分方程数值求解二个方面。
3声波方程
在没有弹性横波只有弹性纵波存在时，对上式二边取散度：
利用旋度与散度的关系： ， ，交换上式的微分次序并化简可得：
其中， 表示纵波速度，其表达式为：
用u表示某一时刻t二维空间任意点（x，z）处的位移，则声波方程可以用u表示为：
为雷克子波的时域波形及其振幅谱，因为雷克子波的能量集中在正波峰处，与反射系数褶积后能较好的反应反射问题。
图1 雷克子波（主频25HZ）
用有限差分方式计算所得的地震道，我们看到，在z=150m处，速度由1200m/s变为2000m/s，密度由1.5g/cm3变为2.0g/cm3，反射系数发生较大变化，所以在0.25s处产生了较大的振幅波动；在z=300m处，速度由2000m/s变为1600m/s，密度不变，同样产生波动，但反射系数的变化不如z=150m处，所以振幅变化也小，同时由于速度变小，反射系数应与z=150m处相反，对应的就是在0.4s处的振幅波动；在z=450m，速度由1600m/s变为2500m/s，密度由2.0g/cm3变为2.6g/cm3，反射系数变化较大，也产生了较大的波动，对应的为0.9s处的振幅波动。理论模型和实际情况吻合，计算结果得到了验证。
4程序实现
（1）均匀介质倾斜地层模型
含有倾斜界面的四层地质模型，整个模型横向x=512m，纵向z=512m。当激发点在模型顶界面中心点处（256，0），接收线在z=0处，按照dx=1m，dz=1m进行网格离散化，震源项所用雷克子波的主频为35Hz，记录时间0.6s。
（2）建模方法
为了能将模型离散化，编制了MAKEMODEL程序，可以根据输入的控制节点建立速度和密度模型。
在编程过程中，考虑了弯曲界面的情况，插值方法采用线性插值法和最小曲率法2种方法，以便在不同情况下选择使用。
线性插值是数学、计算机图形学等领域广泛使用的一种简单插值方法，这种方法是通过在数据点之间连线以建立离散化模型的，模型采用直线连接的方式，在表达水平地层是比较合适，但对于弯曲界面等情况则无法使用。
最小曲率法是一种弯曲界面建模时的好方法，其试图在尽可能严格地尊重数据的同时，生成尽可能圆滑的曲面，广泛用于地球科学。用最小曲率法生成的插值面类似于一个通过各个数据值的，具有最小弯曲量的长条形薄弹性片。使用最小曲率法时要涉及到两个参数：最大残差参数和最大循环次数参数，以用来控制最小曲率的收敛。
（3）程序试算
水平四层模型，整个模型横向x=512m，纵向z=512m，分界面分别在z=100m，z=250m，z=400m处；介质速度由上到下分别为：1200m/s，1800 m/s，2200 m/s，2800 m/s；介质密度由上到下分别为：1.6g/m3，2.0 g/m3，2.2 g/m3，2.6 g/m3，具体如所示。当激发点在模型顶界面中心点处（256，0），接收线在z=0处，按照dx=1m，dz=1m进行网格离散化，震源项所用雷克子波的主频为35Hz，记录时間0.6s。
根据理论计算，三层水平大地反射波的时距曲线应该为双曲线，双曲线的最小值对应炮点所在位置，双曲线的曲率和速度有关，速度越大则曲率越小。T1层反射波所形成双曲线的最小值应在：
T2层反射波所形成双曲线的最小值应在：
T3层反射双曲线的最小值应在：
5结论
本文围绕地震波场数值模拟及波场特征分析等问题，在总结分析前人研究成果的基础上，以波动理论为基础，对有限差分法的一些关键性问题进行了探讨和研究，实现了地震波在几种典型地质模型中传播的数值模拟，并对模拟所得的单炮记录和波场快照进行特征分析。研究结果如下：
1.有限差分法在地震波场数值模拟技术中适用性强、算法容易实现、计算速度快、模拟结果与理论吻合良好、所得到得波场信息丰富，拥有广泛的应用前景；
2.二维有限差分法波场数值模拟作为一个好的辅助工具，可以很方便的计算地质体空间位置和形态大小变化对地震记录的影响，得到对实际工作有指导意义的结论；
3.对大倾斜、小断层、起伏界面、溶腔体等模型进行了计算，得到了明确的地震波场的特征，对模拟结果分析可知：在大倾斜地区，应当根据目的层的空间位置和倾角大小设置合理的观测系统；断层面产生的绕射波对于了解地下断层大致结构和特征有重要意义；凸界面的扩散波和凹界面的回转波是识别弯曲界面的特征波，在解决弯曲界面问题时加加以重视；
参考文献：
[1] 牟永光， 裴正林. 三维复杂介质地震数值模拟北京： 石油工业出版社， 2005. 215
[2] 张钋， 刘洪. 射线追踪方法的发展现状. 地球物理学进展， 2000， 15（1）： 36～45
[3] 刘福田等译. 地震学中的射线方法北京： 地质出版社， 1986.
作者简介：
董丽娜（1983-），女，助理工程师，2007年毕业于黑龙江科技大学，获工学学士学位。