基于MATLAB的地震波信号处理与震相拾取系统

资 源 简 介

本项目基于MATLAB环境开发，旨在提供一套完整的地震监测数据处理解决方案。系统主要功能涵盖三个核心维度：首先是信号的时频分析，利用快速傅氏变换（FFT）算法将时域地震波形转换为频域信号，绘制频谱图以分析地震波的主频及频率成分分布，帮助识别地质特征；其次是震相读取（拾取）功能，系统通过特定的信号突变检测算法（如长短时窗比STA/LTA法或AIC准则），辅助用户识别P波和S波的初至时刻，并提供交互式界面供用户手动微调震相位置，确保数据处理的准确性；最后，项目包含了完整的数据预处理模块，支持去噪、滤波、去趋势等操作。为了方便用户上手，项目内不仅包含核心源代码，还配备了详尽的使用说明文档，详细解释了函数接口、参数设置及运行步骤，适用于地震学研究、地球物理勘探数据分析及相关教学实验。

详 情 说 明

基于MATLAB的地震波信号处理与震相拾取系统

项目简介

本项目是一个基于MATLAB环境开发的地震信号处理与分析系统。针对地震监测数据的特性，系统集成了一套完整的数据处理流程，包括模拟信号生成、信号预处理、频域分析以及自动化的震相拾取（Phase Picking）。该项目不仅实现了经典的STA/LTA初拾取算法，还结合了AIC准则进行P波到达时刻的精细化修正，旨在帮助研究人员和学生理解地震信号处理的核心算法与逻辑。

功能特性

本项目在 main.m 中实际实现了以下核心功能：

• 合成地震数据生成：在缺乏外部数据源的情况下，系统能自动构建包含背景噪声、线性漂移、P波和S波（基于Ricker子波）的合成地震记录。
• 信号预处理：包含去趋势（Detrend）、去均值（Demean）以及零相位带通滤波，有效提升信噪比。
• 频谱分析：利用FFT算法将时域信号转换为频域信号，直观展示信号频率分布。
• 混合震相拾取算法：

* STA/LTA：用于捕捉信号能量突变，进行初步事件检测。 * AIC准则：在初步检测的基础上，利用赤池信息量准则（Akaike Information Criterion）精确锁定P波初至时刻。

• S波估算：基于P波位置和能量最大值原理对S波进行简易估算。
• 多维可视化：系统通过分屏图表展示从原始波形到处理结果的全过程，并计算拾取误差。

系统要求

• MATLAB R2016a 或更高版本
• Signal Processing Toolbox（信号处理工具箱）

使用方法

1. 将 main.m 文件保存在MATLAB的工作路径中。
2. 在MATLAB命令窗口输入 main 并回车，或直接在编辑器中点击“运行”。
3. 系统将自动执行所有计算步骤，弹出可视化窗口，并在控制台输出具体的震相拾取时间及误差统计。

---

技术实现细节

本节详细解析 main.m 代码内部的算法逻辑与实现流程：

1. 模拟数据生成 (generate_synthetic_data)

由于未连接外部数据库，系统内部生成时长为60秒、采样率为100Hz的合成信号：

• 噪声模拟：使用高斯白噪声作为背景干扰。
• 震相模拟：通过Ricker子波（雷克子波）分别模拟P波（10秒到达，主频5Hz）和S波（18秒到达，主频2Hz，振幅更大）。
• 漂移模拟：叠加一个线性趋势项，模拟实际地震计可能存在的基线漂移。

2. 数据预处理 (preprocess_signal)

为了消除干扰，代码严格按照以下顺序处理原始信号：

1. 去趋势与去均值：消除信号中的线性漂移和直流分量。
2. 带通滤波：设计4阶Butterworth滤波器，通带范围设定为 0.5Hz - 15Hz，覆盖主要地震波能量频带。
3. 零相位滤波：使用 filtfilt 函数代替标准的 filter，确保滤波后的信号相较于原始信号没有任何相位延迟（Phase Delay），保证震相拾取的时间准确性。

3. 频域分析 (compute_spectrum)

• 对预处理后的数据进行快速傅氏变换（FFT）。
• 计算双边谱并转换为单边振幅谱（Single-sided Amplitude Spectrum）。
• 对频率轴进行归一化处理，以便在可视化中正确显示对应的赫兹（Hz）值。

4. 震相拾取算法逻辑

#### 第一阶段：STA/LTA 触发检测 (run_stalta)

• 算法原理：计算短时窗（STA, 0.5s）平均能量与长时窗（LTA, 10.0s）平均能量的比值。
• 高效实现：代码利用MATLAB的 filter 函数递归计算滑动平均，避免了低效的 for 循环求和，显著提高了计算速度。
• 触发机制：当比值曲线首次超过设定的阈值（TriggerOn = 3.0）时，标记为触发点。

#### 第二阶段：AIC 精细化拾取 (run_aic_picker)

• 目的：STA/LTA往往对初至的拾取存在滞后，AIC用于修正这一误差。
• 逻辑：以STA/LTA触发点为中心，截取前后各2秒（此处硬编码为2.0s）的窗口数据。
• AIC计算：计算窗口内信号分割点前后的方差对数和：$k cdot log(text{var}_1) + (N-k-1) cdot log(text{var}_2)$。AIC曲线的极小值点被认为是信号统计特性发生改变的最可能时刻（即P波初至）。

#### 第三阶段：S波简易拾取

• 代码实现了一种简化的S波搜索逻辑：仅在P波拾取位置的2秒之后开始搜索，寻找信号振幅绝对值最大的点，将其标记为S波到达时刻。

5. 结果可视化 (visualize_results)

系统生成一个包含四个子图的图形窗口：

• 子图1：原始波形（包含噪声和线性漂移）。
• 子图2：预处理后的波形，并叠加了真实的P波和S波到达时间线（绿色和品红色虚线），用于对比。
• 子图3：FFT频域分析图，展示信号的主要频率成分。
• 子图4：综合展示。左轴显示预处理波形，右轴显示STA/LTA比值曲线。图中标注了算法自动计算出的P波（红色实线）和S波（蓝色实线）拾取时刻。

6. 输出结果

控制台将打印以下信息：

• 系统采样率。
• P波拾取结果：AIC修正后的计算值、合成数据的真实值，以及两者之间的绝对误差（单位：秒）。
• S波估算结果：基于最大能量法推断的S波到达时刻。