一种波动方程正演的瑞利面波频散响应计算方法及其装置

摘要

本发明公开了一种波动方程正演的瑞利面波频散响应计算方法及其装置，通过使用高精度的高阶有限差分、合理设置自由地表边界条件、有效压制模型外边界的伪反射、使用Radon变换法来计算面波的相速度谱，即瑞利面波在频率域中的频散响应，本发明基于各向异性弹性波动方程通过时域有限差分正演模拟来计算瑞利面波频散响应特征，其地质模型可以是各向异性、非水平层状介质，自由地表可以是起伏自由地表，可以用来研究复杂介质的瑞利面波频散响应特征。

说明书

技术领域

本发明涉及地球物理勘探中的地震波数值模拟领域，特别涉及一种波动方 程正演的瑞利面波频散响应计算方法及其装置。

背景技术

弹性波理论证明，当介质为半无限弹性介质时，在自由地表和弹性介质分 界面上将出现一种波，这种沿界面在弹性介质内部传播的不均匀波就是瑞利波。 瑞利波由瑞利于1887年首先指出其存在而得名。瑞利波是纵波与横波在弹性介 质分界面附近干涉的结果。在表层附近，质点的运动轨迹为椭圆。在均匀介质 中，瑞利波的传播速度大约是横波传播速度的0.92倍。瑞利波在弹性介质中的 传播深度大约为1个波长。瑞利面波的传播速度与频率、介质参数的分布有密 切关系：在均匀半空间介质中，瑞利面波的传播速度与频率呈线性关系，无频 散现象；当地下为多于两层的层状介质时，瑞利面波的传播速度与频率不再保 持线性关系，呈现出显著的频散特征，而且呈现多组模式，即含有多条频散曲 线。因此，频散响应能够一定程度上反映地下介质的分布特性，研究瑞利面波 的特征最直观的方式就是研究瑞利面波的频散响应。已知地下介质参数的分布 研究频散响应特征就是瑞利面波的正演问题，通过频散响应特征研究地下介质 参数的分布就是瑞利面波的反演问题。在水平层状介质模型中，瑞利面波的频 散响应特征表现为多条不同模式的频散曲线。利用瑞利面波的频散响应，反演 地下或地球内部的结构，这在工程物探和天然地震领域得到了广泛的应用。

在探测深度介于天然地震和工程物探之间的石油地震勘探中，瑞利面波是 一种很常见的干扰波，在地震剖面上常常表现为强能量、低频低速特性，在以 反射波勘探为主的地震勘探中，一般在地震资料处理当作噪声干扰被压制或消 除。在复杂的山地、戈壁、沙漠等地区，强能量的面波往往掩盖掉了近地表的 大量反射波信息，这对地震资料的静校正和后续处理非常不利。实际上，石油 地震勘探中的瑞利面波也和折射波、反射波一样包含有地下介质的信息，充分 利用瑞利面波信息有可能为研究表层结构开辟一条新的途径。另外，基于模型 的瑞利面波波动方程正演也可以用来消除或压制地震记录的面波干扰。因此， 在石油地震勘探中，将瑞利面波当作一种有用信息，研究瑞利面波的特性，特 别是研究瑞利面波的频散响应特征具有十分重要的意义。

在实现本发明的过程中，发明人发现现有技术至少存在以下缺点：

目前国内外对瑞利面波频散响应特征的研究都是基于水平层状介质中瑞利 面波频散方程，该方程基于弹性波动方程推导，并假设地下介质由多层水平层 状介质组成，自由地表为水平自由地表。该方法计算频散曲线快速，常被用来 做面波正演和面波反演。但是，该方法计算的频散曲线只能反映相位变化，不 能反映振幅强弱的变化；不能解决地下介质为各向异性、粘弹性、双相介质等 更复杂介质；不能解决自由地表为起伏地表的情形；无法研究非水平层状介质 问题。弹性波动方程理论上可以模拟任意复杂介质和复杂地形中的地震波传播， 例如基于波动方程有限差分正演来研究纵波和横波等体波的传播规律，其相应 的应用如地震叠前逆时偏移、全波形反演、基于模型正演去多次波等已经在石 油地球物理勘探领域发挥了重要作用。由于瑞利面波相对于体波传播得更慢、 其能量在深度方向上衰减得很快、其激发对自由地表要求严格等因素，基于波 动方程使用有限差分正演来研究面波要比研究体波复杂得多。只有高精度模拟 了面波传播，才可以利用面波记录来研究复杂介质中面波的频散响应特征。

目前对于面波的频散响应特征研究都是基于水平自由地表、水平层状各向 同性介质，此时瑞利面波的频散方程具有解析解形式；而实际的地下介质绝大 多数为起伏自由地表、非水平层状介质，频散方程很难求出解析解形式。为了 研究复杂介质中瑞利面波的频散响应特征，有必要从最原始的弹性波动方程出 发，来解决上面上述存在的问题。

发明内容

为了实现对复杂介质的瑞利面波频散响应特征的研究，解决地下介质为各 向异性、粘弹性、双相介质等更复杂介质，解决自由地表为起伏地表的情形等 非水平层状介质问题，本发明提供了一种波动方程正演来计算瑞利面波频散响 应的方法，通过使用高精度的高阶有限差分、合理设置自由地表边界条件、有 效压制模型外边界的伪反射、使用相移法来计算面波的相速度谱，即瑞利面波 在频率域中的频散响应。

本发明实施例提供了一种波动方程正演的瑞利面波频散响应计算方法，所 述方法包括以下步骤：

步骤一，根据波动方程建立地下介质的正演模型，并将地下介质的介质参 数转换为波动方程张量矩阵中的弹性系数；

步骤二，对正演模型中的弹性波场分量和弹性系数进行二维网格离散化， 所述弹性波场分量和弹性系数交错地位于网格节点上；根据波速分布的范围 [V_min，V_max]和给定的峰值频率f_peak，设置正演模拟的时间迭代步长dt、离散网格 大小和空间网格离散步长dx、dz；

步骤三，设定自由地表边界条件和不同介质间的边界条件，

所述介质密度为零时，其对应的位移分量为零；

所述弹性系数C₅₅所涉及的其中一个相邻介质的弹性系数为零时，其对应的 剪应力为零；

步骤四，通过弹性波动方程进行时域有限差分正演模拟，在地面上接收包 含面波的地震记录zz(x，t；z＝0)，并将该地震记录变换到频率域

步骤五，对步骤四中得到的频率域进行瑞利面波速度谱计算，得到面波地 震记录在频率域的速度谱s(1:nω，1:nυ)，即瑞利面波在频率域中的频散响应。

进一步，所述步骤一中所建正演模型公式如下所示：

$\left(\begin{array}{c}\rho \frac{{\partial v}_x}{\partial t}=\frac{{\partial \tau }_{\mathrm{xx}}}{\partial x}+\frac{{\partial \tau }_{\mathrm{xz}}}{\partial z}\\ \rho \frac{{\partial v}_z}{\partial t}=\frac{{\partial \tau }_{\mathrm{xz}}}{\partial x}+\frac{{\partial \tau }_{\mathrm{zz}}}{\partial z}\\ \frac{{\partial \tau }_{\mathrm{xx}}}{\partial t}=C_{11}\frac{{\partial v}_x}{\partial x}+C_{33}\frac{{\partial v}_z}{\partial z}\\ \frac{{\partial \tau }_{\mathrm{zz}}}{\partial t}=C_{13}\frac{{\partial v}_x}{\partial x}+C_{33}\frac{{\partial v}_z}{\partial z}\\ \frac{{\partial \tau }_{\mathrm{xz}}}{\partial t}=C_{55}\frac{{\partial v}_x}{\partial z}+C_{55}\frac{{\partial v}_z}{\partial x}\end{array}\right)$

其中，介质参数包括：

V_p0——垂向纵波速度；

V_s0——垂向横波速度；

——介质密度；

和——各向异性参数，

所述弹性系数同介质参数之间的关系为：

$C_{33}={\mathrm{\rho V}}_{p0}^2$

$C_{55}={\mathrm{\rho V}}_{S0}^2$

$C_{11}=(1+2ϵ){\mathrm{\rho V}}_{p0}^2$

$C_{13}=\rho \sqrt{(V_{P0}^2-V_{s0}^2)\left[(1+2\delta )\right](V_{P0}^2-V_{s0}^2)}-{\mathrm{\rho V}}_{s0}^2$

进一步，所述步骤三中的两种所述边界条件均采用介质平均法计算位于相 邻两离散网格节点中间的弹性系数。

进一步，所述步骤四中对正演模型进行时间偏导数和空间位移偏导数的计 算，所述时间偏导数采用二阶中心差分进行离散，所述空间位移偏导数采用十 二阶中心差分进行离散；通过对正演模型进行计算，分别得出所述弹性波场分 量的二阶时间差分精度和十二阶空间差分精度的数值计算迭代公式。

进一步，所述步骤四中，在正演模型的左边界区域、右边界区域和下边界 区域设置完全匹配层，所述完全匹配层的厚度为8～12层，在对应的一阶空间 偏导数处设置一个辅助变量，用于吸收来自外边界的反射。

进一步，所述步骤五中，根据正演模型介质中面波速度的分布范围确定频 率域中需要分析的单频个数，并对每个频率通过相移法计算得出所对应的速度， 通过对每个频率和速度进行循环，得到面波的地震记录在所述频率域的速度谱 s(1:nω，1:nυ)。

优选地，所述步骤二中，一个波长内至少被离散成8个离散网格，一个地 震子波周期至少被离散为10个时间采样点。

进一步，将频率域作振幅一致性处理，用于对低频能量和高频 能量进行补偿，其实现公式如下：

另外本发明还提供了一种波动方程正演的瑞利面波频散响应计算装置，所 述装置包括：

弹性系数正演模型模块，根据波动方程建立地下介质的正演模型，并将地 下介质的介质参数转换为波动方程张量矩阵中的弹性系数；

正演模型网格离散模块，用于对正演模型中的弹性波场分量和弹性系数进 行二维网格离散化，并设置时间迭代步长和空间离散步长；

边界条件设置模块，用于设定自由地表边界条件和不同介质间的边界条件；

有限差分正演模拟模块，通过弹性波动方程进行时域有限差分正演模拟， 在地面上接收包含面波的地震记录zz(x，t；z＝0)，并进行地震记录的计算；

面波记录变换域模块，用于将地震记录zz(x，t；z＝0)变换到频率域

及面波速度谱计算模块，通过相移法计算频率域面波速度谱，得到面波地 震记录在频率域的速度谱s(1:nω，1:nυ)，即瑞利面波在频率域中的频散响应。

所述边界条件设置模块包括：自由地表边界条件模块、不同介质间边界条 件模块及正演模型边界匹配层模块，所述正演模型边界匹配模块用于在正演模 型的左边界、右边界和下边界设置完全匹配层。

本发明实施例提供的技术方案的有益效果是：

常规方法基于瑞利面波频散方程来计算面波频散响应特征，其地质模型只 能是各向同性、水平层状介质模型，自由地表只能是水平自由地表。本发明基 于各向异性弹性波动方程通过时域有限差分正演模拟来计算瑞利面波频散响应 特征，其地质模型可以是各向异性、非水平层状介质，自由地表可以是起伏自 由地表，因此可以用来研究复杂介质的瑞利面波频散响应特征。针对水平层状 介质模型的频散响应特征，如图6所示，常规方法与本发明的结果保持一致， 说明了本发明的正确性、可信性。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所 需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明 的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下， 还可以根据这些附图获得其他的附图。

图1弹性波波动方程有限差分交错网格离散示意图；

图2本发明瑞利面波频散响应计算流程图；

图3本发明瑞利面波频散响应计算装置；

图4瑞利面波在0.5s处传播的波场快照；

图5瑞利面波在1.0s处传播的波场快照；

图6三层水平层状介质正演合成的面波记录；

图7水平层状介质模型的瑞利面波频散响应；

图8水平层状介质模型的频散曲线比较(虚线为本发明，实线为理论解)；

图9各向异性介质模型；

图10各向异性介质模型的瑞利面波频散响应；

图11含起伏反射界面的层状介质模型；

图12含起伏反射界面的层状介质模型的面波频散响应；

图13含起伏自由地表的层状介质模型；

图14含起伏自由地表的层状介质模型的瑞利面波频散响应。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明 实施方式作进一步地详细描述。

图2为本发明计算方法的总体流程图，其具体实现原理如下：

在二维各向异性弹性介质中，一阶位移-应力弹性波动方程为：

$\left(\begin{array}{c}\rho \frac{{\partial v}_x}{\partial t}=\frac{{\partial \tau }_{\mathrm{xx}}}{\partial x}+\frac{{\partial \tau }_{\mathrm{xz}}}{\partial z}\\ \rho \frac{{\partial v}_z}{\partial t}=\frac{{\partial \tau }_{\mathrm{xz}}}{\partial x}+\frac{{\partial \tau }_{\mathrm{zz}}}{\partial z}\\ \frac{{\partial \tau }_{\mathrm{xx}}}{\partial t}=C_{11}\frac{{\partial v}_x}{\partial x}+C_{33}\frac{{\partial v}_z}{\partial z}\\ \frac{{\partial \tau }_{\mathrm{zz}}}{\partial t}=C_{13}\frac{{\partial v}_x}{\partial x}+C_{33}\frac{{\partial v}_z}{\partial z}\\ \frac{{\partial \tau }_{\mathrm{xz}}}{\partial t}=C_{55}\frac{{\partial v}_x}{\partial z}+C_{55}\frac{{\partial v}_z}{\partial x}\end{array}\right)$

式中：x为水平方向坐标变量，z为垂直方向坐标变量；ρ＝ρ(x，z)为介质密度， v_x＝v_x(x，z，t)为质点的水平位移，v_z＝v_z(x，z，t)为垂直位移； τ_xx＝τ_xx(x，z，t)，τ_zz＝τ_zz(x，z，t)为正应力，τ_xz＝τ_xz(x，z，t)为剪应力，C₁₁、C₁₃、C₃₃、C₅₅为 介质的弹性系数张量矩阵元素，各向异性介质的速度和各向异性的强度隐含在 弹性参数中，因此需要转换。

根据Thomsen公式，密度、垂向纵波速度V_p0、垂向横波速度V_s0、纵波各向 异性参数ε、各向异性参数δ与张量矩阵弹性系数之间的关系为：

$C_{33}={\mathrm{\rho V}}_{p0}^2$

$C_{55}={\mathrm{\rho V}}_{S0}^2$

$C_{11}=(1+2ϵ){\mathrm{\rho V}}_{p0}^2$

$C_{13}=\rho \sqrt{(V_{P0}^2-V_{s0}^2)\left[(1+2\delta )\right](V_{P0}^2-V_{s0}^2)}-{\mathrm{\rho V}}_{s0}^2$

接下来，采用交错网格离散方案，对弹性波场变量和弹性系数的空间位置 在笛卡尔坐标系中进行网格离散化，并简记：

i·x+Δx，i+0.5·x+0.5Δx

k·z+Δz，k+0.5·z+0.5Δz

式中Δx、Δz分别为水平方向和深度方向的空间网格离散步长。由于弹性系 数离散时位于整网格节点上，考虑到有可能涉及到错半个网格的网格节点上弹 性系数的计算，在这里采用介质平均法来实现：

$\rho (i,k+0.5)=\frac{\rho (i,k)+\rho (i,k+1)}{2}$

$\rho (i+0.5,k)=\frac{\rho (i,k)+\rho (i+1,k)}{2}$

$C_{55}(i+0.5,k+0.5)=\frac{C_{55}(i,k)+C_{55}(i+1,k)+C_{55}(i,k+1)+C_{55}(i+1,k+1)}{4}$

同时对弹性波场变量的迭代更新做以下约束：

如果ρ(i，k+0.5)＝0，则v_z(i，k+0.5)＝0；

如果ρ(i+0.5，k)＝0，则v_x(i+0.5，k)＝0；

如果C₅₅(i，k)、C₅₅(i+1，k)、C₅₅(i，k+1)、C₅₅(i+1，k+1)中只要其中任何一个为零， 那么C₅₅(i+0.5，k+0.5)强制为零，τ_xz(i+0.5，k+0.5)＝0。

上述介质平均法和特殊规定就能够实现起伏自由地表边界条件和两种介质 分界面处的边界条件。

在数值实现弹性波动方程时，涉及到时间偏导数和空间偏导数的计算。时 间偏导数用二阶中心差分进行离散，

$\frac{\partial f}{\partial t}\approx \frac{f(t+\mathrm{\Delta t})-f(t+\mathrm{\Delta t})}{\mathrm{\Delta t}}$

其中f代表弹性波场分量v_x，v_z，τ_xx，τ_zz，τ_xz。

空间偏导数用十二阶中心差分进行离散，

$\frac{\partial f}{\partial h}\approx \frac{1}{\mathrm{\Delta h}}\underset{m=1}{\overset{6}{\Sigma }}{D}_6^m\left\{f\right[h+\frac{\mathrm{\Delta h}}{2}(2m-1)]-f[h-\frac{\mathrm{\Delta h}}{2}(2m-1)\left]\right\}=L_h\left[f\right]$

其中，h代表x或z，L_h[f]表示数值离散后的f对h的一阶导数。差分系数 计算公式为

$D_6^m=\frac{{(-1)}^{m+1}\underset{i=1,i\ne n}{\overset{6}{\Pi }}{(2i-1)}^2}{(2m-1)\underset{i=1}{\overset{5}{\Pi }}[{(2m-1)}^2-{(2i-1)}^2]},(m=\mathrm{1,2,3,4,5,6})$

记fⁿ(i，k)＝f(nΔt，iΔx，kΔz)＝f(t；x，z)，以正演模型公式中第三个方 程的数值离散为例，可表示为

${\tau }_{\mathrm{xx}}^{n+1}(i,k)={\tau }_{\mathrm{xx}}^n(i,k)+\frac{C_{11}(i,k)}{\mathrm{\Delta t}}{L}_x\left[v_x^{n+0.5}(i,k)\right]+\frac{C_{13}(i,k)}{\mathrm{\Delta t}}{L}_z\left[v_z^{n+0.5}(i,k)\right]$

同理可得正演模型公式中其它4个方程的二阶时间差分精度、十二阶空间 差分精度的数值计算迭代公式。

在模型的左边界区域、右边界区域和下边界区域设置完全匹配层吸收来自 外边界的反射，其规律是：为了吸收某一方向的伪反射，在对应的一阶空间偏 导数处引入一个辅助变量。以正演模型公式中的第三个方程所在模型的左边界 为例，在x轴设置完全匹配层的厚度为10层离散网格，并在关于x的空间偏导 数项L_x[v_x(i，k)]处引入一个辅助变量有：

${\tau }_{\mathrm{xx}}^{n+1}(i,k)={\tau }_{\mathrm{xx}}^n(i,k)+\frac{C_{11}(i,k)}{\mathrm{\Delta t}}\left\{L_x\right[v_x^{n+0.5}(i,k)]-E_{\mathrm{xx}}^n(i,k)\}+\frac{C_{13}(i,k)}{\mathrm{\Delta t}}{L}_z\left[v_z^{n+0.5}(i,k)\right]$

(i＝0，9，k＝0，nz)

其中下一时刻的计算可以使用下面公式来进行迭代更新

$E_{\mathrm{xx}}^{n+1}(i,k)=e^{-{\sigma }_i\mathrm{\Delta t}}{E}_{\mathrm{xx}}^n(i,k)+(1-e^{-{\sigma }_i\mathrm{\Delta t}})L_x\left[v_x^{0.5}(i,k)\right]$

(i＝0，9，k＝0，nz)

其中${\sigma }_i=0.14*{\left(\frac{i+1}{10}\right)}^2,(i=\mathrm{0,1,2},·,9)$

在水平自由地表z＝0处记录每个采样步长内的波场，最终我们获得了面波 正演记录τ_zz(1:nt，1:nx)。

将该记录变换到频率域得到由于使用的子波为带限雷克子波， 在面波记录的低频部分或高频部分能量比较弱，需要对低频能量或高频能量进 行补偿，即振幅一致性处理，实现公式为：

确定正演模型介质中面波速度的分布范围上限υ_min和下限υ_max，然后等分成 nυ份，则第iυ等份处的速度变量为

$\upsilon ={\upsilon }_{\min }+\frac{\mathrm{i\upsilon }}{\mathrm{n\upsilon }-1}({\upsilon }_{\max }-{\upsilon }_{\min }),(\mathrm{i\upsilon }=\mathrm{0,1},·,\mathrm{n\upsilon }-1)$

然后，对于每个频率，使用Radon变换法计算速度谱：

其中x＝(ix-1)dx，(ix＝0，1，·，nx-1)

对频率和速度进行循环，最终得到面波记录在频率域的速度谱s(1:nω，1:nυ)。 在该速度谱中的振幅谱中，观测能量团强弱分布，拾取能量团的峰 值曲线，就可以得到面波的频散曲线

首先设计模型。该模型含有多种介质，每种介质包含垂向纵波速度V_p0、垂 向横波速度V_s0和密度、各向异性参数和等五个参数。为了正演计算方便， 将这五个参数转换为波动方程张量矩阵中的弹性系数。根据速度分布的范围 [V_min，V_max]和给定的峰值频率f_peak，来设置正演模拟的时间迭代步长dt和空间网 格离散步长dx、dz。一个波长内至少被离散成8个网格，同时一个地震子波周 期至少得离散为10个时间采样点。根据模型的大小，以及设定的时间采样步长 dt和空间采样步长dx、dz，对模型进行网格离散化，离散后的网格大小为nx*nz。

然后基于弹性波动方程使用有限差分进行正演模拟，并记录面波。主要涉 及四方面的内容：一是弹性波场变量和弹性参数在离散网格上的分布；二是自 由地表边界条件的实现；三是有限差分的精度；四是吸收边界条件的实现。

如图1所示，弹性波场变量交错地分布于离散网格上。其中正应力分量xx 和zz位于实线相交的网格节点上，剪应力分量xz位于虚线相交的网格节点上， 位移分量x和z位于实线和虚线相交的网格节点上。弹性系数、C₁₁、C₁₃、C₃₃和C₅₅皆位于实线相交的网格节点上，即与正应力在相同的网格节点上。剪应力 分量和位移分量所在位置如果涉及到弹性参数，就使用介质平均法来计算弹性 参数。

采用介质平均法，使得剪应力波场分量位于不同介质的分界面处。同时规 定凡是密度为零的地方，其对应的位移分量也为零；弹性系数C₅₅涉及的相邻介 质只要其中一个为零，其对应的剪应力也为零。起伏自由地表处的边界条件和 不同介质分界处的边界条件得以自动实现。

为了保持计算效率和计算精度，时间偏导数的计算采用二阶中心差分近似， 空间偏导数的计算采用十二阶中心差分近似。

为了压制有限边界模型在有限边界处引起的虚假发射，在模型的左边界、 右边界和下边界设置完全匹配层来消除虚假反射。一般完全匹配层的厚度设置 为10层。

通过弹性波动方程有限差分正演模拟，在地面上接收到了包含面波的地震 记录zz(x，t；z＝0)。接下来，将该地震记录变换到频率域，对单个频率成分在水平 方向上做归一化处理，即对振幅做归一化处理。

确定正演模型介质中面波速度的分布范围上限υ_min和下限υ_max，然后等分成 nυ份；确定需要分析的单频个数。对于每个频率，使用相移法计算速度谱。对 每个频率和速度进行循环，最终得到面波记录在频率域的速度谱s(1:nω，1:nυ)。

本发明还提供了一种波动方程正演的瑞利面波频散响应计算装置，所述装 置包括：

弹性系数正演模型模块，根据波动方程建立地下介质的正演模型，并将地 下介质的介质参数转换为波动方程张量矩阵中的弹性系数；

正演模型网格离散模块，用于对正演模型中的弹性波场分量和弹性系数进 行二维网格离散化，并设置时间迭代步长和空间离散步长；

边界条件设置模块，用于设定自由地表边界条件和不同介质间的边界条件；

有限差分正演模拟模块，通过弹性波动方程进行时域有限差分正演模拟， 在地面上接收包含面波的地震记录zz(x，t；z＝0)，并进行地震记录的计算；

面波记录变换域模块，用于将地震记录zz(x，t；z＝0)变换到频率域

及面波速度谱计算模块，通过相移法计算频率域面波速度谱，得到面波地 震记录在频率域的速度谱s(1:nω，1:nυ)，即瑞利面波在频率域中的频散响应。

所述边界条件设置模块包括：自由地表边界条件模块、不同介质间边界条 件模块及正演模型边界匹配层模块，所述正演模型边界匹配模块用于在正演模 型的左边界、右边界和下边界设置完全匹配层。

另外如图3所示，本发明还提供了一种波动方程正演的瑞利面波频散响应 计算装置，包括：

弹性系数正演模型模块，根据波动方程建立地下介质的正演模型，并将地 下介质的介质参数转换为波动方程张量矩阵中的弹性系数；

正演模型网格离散模块，用于对正演模型中的弹性波场分量和弹性系数进 行二维网格离散化，并设置时间迭代步长和空间离散步长；

边界条件设置模块，用于设定自由地表边界条件和不同介质间的边界条件；

弹性波时域有限差分正演模拟模块，通过弹性波动方程进行时域有限差分 正演模拟，在地面上接收包含面波的地震记录zz(x，t；z＝0)，并进行地震记录 的计算；

面波记录变换域模块，用于将地震记录zz(x，t；z＝0)变换到频率域

及面波速度谱计算模块，通过相移法计算频率域面波速度谱，得到面波地 震记录在频率域的速度谱s(1:nω，1:nυ)，即瑞利面波在频率域中的频散响应。

其中边界条件设置模块包括：自由地表边界条件模块、不同介质间边界条 件模块及正演模型边界匹配层模块，所述正演模型边界匹配模块用于在正演模 型的左边界、右边界和下边界设置完全匹配层。

实施例1水平层状介质水平自由地表的面波频散响应

设计一个含3层介质的水平层状介质模型，自由地表为水平自由地表。模 型大小为为562米×100米，被离散为750*200的网格，其中网格水平间距为 0.75米，网格垂直间距为0.5米；第一层垂向纵波速度为650m/s，横波速度为 200m/s，密度为1.82×103Kg/m³，2个各向异性参数皆设置为0。在水平自由地 表处垂直震源激发，水平自由地表接收，主频为30Hz。图4为瑞利面波在0.5s 处传播的波场快照，图5为瑞利面波在1.0s处的波场快照。图6是正演模拟合 成的面波记录。图7是该面波记录在频率域的频散响应。图8是常规方法和本 发明计算的频散曲线比较，实现是常规方法，虚线是本发明。

实施例2各向异性水平层状介质的面波频散响应

设计一个如图9所示的含3层介质的各向异性水平层状介质模型，自由地 表为水平自由地表。与实施例1中弹性系数不同的是，各向异性参数ε为0.1、 各向异性参数δ为0.5。网格大小等其它参数的设置同实施例1。图10为各向 异性介质模型的瑞利面波频散响应。

实施例3含起伏反射界面层状介质模型的面波频散响应

设计一个如图11所示的含起伏反射界面的层状介质模型，各向同性介质， 自由地表为水平自由地表。共3层介质，弹性系数的具体设置同实施例1设置。 图12是含起伏反射界面的层状介质模型的面波频散响应。

实施例4含起伏自由地表层状介质模型的面波频散响应

设计一个如图13所示的含起伏自由地表的层状介质模型，各向同性介质， 自由地表为起伏自由地表。共3层介质，弹性系数的具体设置同实施例1设置。 图14是含起伏自由地表的层状介质模型的面波频散响应。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的 精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的 保护范围之内。