波动方程求解法正演模拟

Question

如题所述

Answer 1

波动方程求解进行正演模拟常用的方法有有限差分法、傅立叶变换法、克希霍夫积分法以及有限元素法。所有叠后波动方程偏移法均可用于正演模拟而得到自激自收剖面，所不同的是只需要改变延拓方向，由地下向地面延拓，取地面上各点计算值作为计算剖面。

4.4.1 有限差分法二维正演模拟

有限差分法二维正演模拟是利用有限差分法求解二维波动方程定解问题。在计算机上实现这一问题，是在有限区域内进行的，如图4.14所示。

图4.14 有限差分法计算区域示意图

由于取了有限区域，所以在两侧和底界上人为地造成了边界，这个边界叫作计算边界。因为它的存在，波传播碰到边界时就必然产生反射（如图4.14 中所示），此反射造成干扰，影响理论记录的质量。为消除或减弱边界反射效应，需在定解问题中加吸收边界条件。

所考虑的计算区域为D=｛（x，z）|，-α≤x≤α，0≤z≤b，0≤t≤T｝，含吸收边界条件的定解问题为（仅考虑纵波方程）：

地球物理数据处理教程

S（x，z，t）是震源函数，S（x，z，t）=δ（x-x₀）δ（z-z₀）·b（t），b（t）是地震子波，（x₀，z₀）是震源坐标。

关于上面定解问题中的吸收边界条件（4.4.4）、（4.4.5）和（4.4.6）的详细论述请参见本书后面的附录。

用差分法求解上述的定解问题，需建立差分格式（应满足收敛性及稳定性条件，这方面的详细论证请参阅有关文献），把用偏微分方程表示的定解问题转换为用差分方程表示。差分方程的建立过程在此不再重复（与有限差分法波动方程偏移中差分方程建立方法相同），这里直接给出差分方程式：

地球物理数据处理教程

地球物理数据处理教程

右边界条件：

地球物理数据处理教程

j=1，2，…，J；n=1，2，…，N （4.4.8）

左边界：

地球物理数据处理教程

j=1，2，…，J+1；n=1，2，…，N （4.4.9）

底边界：

地球物理数据处理教程

i=1，2，…，I；n=1，2，…，N （4.4.10）

震源子波可取b（t）=

cos2πft。在计算机实现时，首先需将给定的几何界面和物性参数数值化，然后通过迭代计算

（2≤i≤I；2≤j≤J；0≤n≤N），最后输出

（i=1，2，3，…，I；n=0，1，2，…，N）作为地面共炮点地震记录。

4.4.2 傅立叶变换法合成记录

已知纵波波动方程

地球物理数据处理教程

根据傅立叶变换，式（4.4.11）的解可表示为：

地球物理数据处理教程

满足频散关系：

（k_x，k_z，t=0）是t=0时刻波场（初始波场）的傅立叶变换。]]

地球物理数据处理教程

亦即

地球物理数据处理教程

为书写简便，函数u中的变量x、z略写。

显然，如果已知t及t-Δt时刻的波场值，可推得t+Δt时刻的波场值。

由台劳级数，u（t）可表示为

地球物理数据处理教程

当t₀=0时

地球物理数据处理教程

所以

地球物理数据处理教程

取n=3，则

地球物理数据处理教程

对式（4.4.11）两端关于t求导得：

地球物理数据处理教程

式（4.4.11）的初始条件为：

地球物理数据处理教程

当t=0时，由（4.4.11）式有：

地球物理数据处理教程

当t=0时，由（4.4.15）式有：

地球物理数据处理教程

将式（4.4.16）、（4.4.17）、（4.4.18）代入式（4.4.14）得到u（Δt）。再将u（0）、u（Δt）代入式（4.4.13），递推出u（2Δt），如此递推下去，便可求得任意时刻的解。

方程（4.4.13）中的空间导数也可由傅立叶变换得到，即：

地球物理数据处理教程

同样，（4.4.17）式和（4.4.18）式中的空间导数也可由傅立叶变换得到，即：

地球物理数据处理教程

应用递推公式必须满足稳定性及允许误差的条件。应用傅立叶变换法，同样也可把上行波方程来作模型，其方法与上面类似，在此不再赘述。

4.4.3 克希霍夫积分法非自激自收记录正演模拟

由波动理论，克希霍夫积分法公式是惠更斯-菲涅尔原理的定量表达式。地下界面可看作由许多个小面元组成，每个小面元产生的绕射波被地面观测点接收。利用克希霍夫积分公式可以推导出正演模拟公式，二维情况是三维情况的简化，亦即把小面元化为小线元，如图4.15所示。

通过对三维克希霍夫积分公式的简化，忽略一些较小量值的项，引入有关的几何量，可得到二维情况的克希霍夫积分非自激自收记录正演模拟公式

地球物理数据处理教程

式中：（r₁+r₂）/V=t是由震源O点经小线元Δl绕射后传播到接收点G的旅行时间；θ₁和θ₂分别为入射线和绕射线与小面元法线间的夹角；β为绕射线到G点的入射角；r₁是震源O到界面入射波的传播距离；r₂是从界面到接收点G绕射波的传播距离；V是速度；Δl是小线元长度；s（t）是地震子波。有关吸收边界条件的详细推导请见附录C。

图4.15 克希霍夫积分法计算示意图

左为三维情况，右为二维情况