反射率法弹性波正演与Tau-P域合成记录模拟

资 源 简 介

本项目旨在基于经典反射率法（Reflectivity Method）理论，开发一套用于模拟水平层状介质全波场响应的地震正演软件。系统首先读取用户定义的层状地质模型参数，包括每一层的纵波速度、横波速度、密度、层厚度以及介质的品质因子（Q值），以此构建弹性介质模型。核心算法采用Kennett矢量递推矩阵方法，在频率-慢度（f-p）域内严格求解弹性波动方程，计算反射系数与透射系数矩阵；该方法能够精确模拟地下介质中的一次反射波、层间多次波、自由表面多次波以及P-S转换波等各种类型的波场，主要用于高精度的AVO（振幅随偏移距变化）分析和薄层调谐效应研究。在获得频率-慢度域的波场响应后，系统将其与用户指定的震源子波（如雷克子波）频谱进行卷积。最后，依据需求对频域数据执行一维反傅里叶变换（从角频率ω变换到截距时间τ），从而直接生成并输出τ-p域（Tau-P domain）的合成地震记录。该系统不仅能直观展示平面波在层状介质中的传播特征，也能为地震数据处理中的Tau-P反褶积和速度分析提供标准的理论数据模型。

详 情 说 明

反射率法地震波场正演与Tau-P域变换模拟系统

项目简介

本项目是一套基于MATLAB开发的地震波场正演模拟软件。系统基于经典反射率法（Reflectivity Method）理论，采用高效的Kennett矢量递推矩阵算法，专门用于模拟水平层状弹性介质中的全波场响应。

该系统在频率-慢度（f-p）域内严格求解弹性波动方程，不仅考虑了纵波（P波）和横波（S波）的传播，还集成了介质的品质因子（Q值）以模拟衰减效应。通过在频率-慢度域进行震源子波卷积和一维反傅里叶变换，系统能够直接生成高精度的Tau-P域（截距时间-射线参数域）合成地震记录，为AVO分析、薄层调谐效应研究以及速度分析提供标准的理论数据模型。

核心功能

• 层状弹性介质建模：支持定义多层地质模型，参数包含纵波速度、横波速度、密度、层厚度及Q值（Qp, Qs）。
• 全波场响应计算：基于Kennett反射率法，精确计算包含一次波、层间多次波及转换波的系统响应。
• 介质衰减模拟：通过引入复速度模型，模拟地震波在非完全弹性介质中的能量衰减与频散。
• 复频率技术应用：利用复频率参数（Complex Frequency）有效压制时间域混叠效应（Time Aliasing）。
• Tau-P域合成记录：直接从f-p域经过反傅里叶变换生成Tau-P域地震记录，直观展示平面波传播特征。
• 多维度可视化：提供地层模型、震源子波、f-p域振幅谱及Tau-P域地震剖面的综合展示。

系统要求与使用方法

• 运行环境：MATLAB (推荐 R2018b 及以上版本)
• 工具箱依赖：信号处理工具箱 (Signal Processing Toolbox) - 用于FFT/IFFT及基础信号操作。
• 使用方式：直接运行主脚本即可启动模拟全过程，程序会自动完成模型构建、波场计算及图形绘制。

算法原理与代码实现细节

本项目的主程序逻辑严密，严格按照地震正演的物理过程编写，主要包含以下五个核心处理步骤：

1. 模型构建与参数定义

程序首先定义了一个典型的4层地质模型（包含基底半空间）。每一层的物性参数被存储在一个矩阵中，包含纵波速度(Vp)、横波速度(Vs)、密度(Rho)、层厚(H)以及品质因子(Qp, Qs)。同时，定义了模拟的关键参数，包括时间采样率、最大记录时间、震源主频，以及Tau-P变换的核心参数——慢度（Ray Parameter, p）的范围和采样间隔。

2. 数据准备与坐标轴初始化

在进入核心计算前，程序自动计算并生成了后续运算所需的基础坐标轴：

• 时间轴与频率轴：根据最大记录时间和采样间隔，确定FFT点数（确保为偶数），并生成对应的时间序列和频率序列。
• 慢度轴：构建从0到最大慢度的线性向量，用于扫描不同的射线参数。
• 震源子波：通过雷克子波（Ricker Wavelet）发生器生成时域子波，并利用FFT将其转换至频率域，为频域卷积做准备。

3. 核心计算：Kennett反射率法模拟

这是系统的计算核心。程序通过双重循环结构（或向量化处理）进行波均场的求解：

• 复速度计算：为了引入介质的衰减特性，程序依据各层的Q值，利用常Q模型的近似公式将弹性速度转换为复速度。
• 慢度循环：程序遍历每一个慢度值 $p$（模拟不同入射角的平面波）。
• 频率循环与Kennett递推：在每一个慢度下，程序对有效频带内的频率进行扫描。为防止时间混叠，频率变量被转化为复频率（引入衰减因子 $sigma$）。
• 递归求解：调用各向同性介质的Kennett矢量递推算法，从底层半空间向上递推至地表，计算总反射系数矩阵。程序目前主要提取并存储垂直分量（对应P-P反射及转换波投影）的频率响应。

4. 合成地震记录生成

完成频率-慢度（f-p）域的响应矩阵计算后，程序执行以下步骤生成最终记录：

• 频域卷积：将计算得到的介质冲激响应谱与震源子波频谱相乘。
• 频谱重构：为了进行反傅里叶变换，程序利用共轭对称性，将单边频谱扩展为包含负频率的双边全频谱。
• Tau-P变换：对频率轴执行一维反傅里叶变换（IFFT），将数据从角频率 $omega$ 变换到截距时间 $tau$，从而得到 $tau-p$ 域数据。
• 增益恢复：消除在计算过程中为压制混叠而引入的复频率衰减项（乘以 $exp(sigma t)$），恢复真实的波场振幅。

5. 结果可视化

系统最终生成一个包含四个子图的综合图表：

• 地层速度模型：展示深度与层速度的阶梯状分布。
• f-p 域波场响应谱：以热图形式展示不同慢度和频率下的波场能量分布，色调反映振幅强弱。
• 震源子波：显示所使用的雷克子波时域波形。
• Tau-P 域合成地震记录：最终的输出结果，展示不同慢度下的地震记录（波形摆），采用红蓝地震色标，并进行了自动增益控制以便于观察弱反射信号。

关键参数说明

在代码中主要涉及以下关键物理参数的配置：

• sim_params.dt: 时间采样间隔，决定了模拟的最高频率（奈奎斯特频率）。
• sim_params.p_max: 最大慢度，决定了模拟的最大入射角度范围。
• sim_params.sigma: 复频率衰减因子，用于控制时间域折叠噪声的压制程度。
• model_data: 矩阵形式的地质模型，每行代表一层，能够灵活修改以适应不同的地质目标。