单台三分量背景噪声数据如何分析极化与信号源背方位角？

当然是DOP-E啦(基于论文 Schimmel et al., 2011; Berbellini et al., 2019)。

总体流程概述

  1. 数据预处理（去趋势、去仪器、对齐）
  2. S-transform 时频分解
  3. 构造谱矩阵与时域平滑
  4. 本征分解求极化参数（DOP、planarity、semimajor axis）
  5. 筛选高质量段（DOP > 阈值 & planarity 接近水平）
  6. 计算背方位角（BAZ）并消除 180° 二义性
  7. 输出可视化与统计（时频图、直方图、极坐标图）
  8. 验证与质量控制（QA）
  9. 批量与性能优化

数据预处理

目标

确保三分量（Z/N/E）记录在同一时间基线，且仅保留稳定的噪声段。

步骤

  1. 去仪器响应

  2. 去仪器响应，转换到速度记录。

  3. 去趋势 / 去均值

  4. trace.detrend("linear"); trace.detrend("demean")

  5. 剔除事件或脉冲干扰

  6. 可通过短时能量（STA/LTA）检测排除瞬变信号（例如地震）。

  7. 时间对齐与补零

  8. 三分量必须同起止时刻；短 gap 可补零，长 gap 建议跳过。

S-transform

定义

$$ S(t,f)=\int x(\tau)\frac{|f|}{\sqrt{2\pi}}e^{-(t-\tau)^2f^2/2}e^{-i2\pi f\tau}d\tau $$

参数与理由

参数 推荐值 理由
频率点数 n_freqs 60 覆盖 0.01–0.06 Hz 区间，低频对数分辨率更合理
窗宽 σ 1/(2πf) 论文定义，时间-频率能量守恒
调整系数 k 1.0–1.5 增大可提高时域平滑性
降采样步长 step 5 s 提高计算效率，仍能解析低频信号

说明：
σ = 1/(2πf) 表示随频率增加窗宽减小，确保多分辨率特性。
|f|/√(2π) 为归一化系数，维持 Parseval 能量一致性。

谱矩阵构造与平滑

谱矩阵定义

$$ S_{ij}(t,f) = \widetilde{X_i}(t,f)\widetilde{X_j}^*(t,f) $$ 其中 $ \widetilde{X_i}(t,f) $ 为分量 i 的 S-transform 复系数。

平滑

对时间方向做高斯平滑： $$ \sigma_{\text{tf}} = \text{tf_window_periods} \times T \times f_s $$

参数 推荐值 理由
tf_window_periods 3.0 论文示例值（2s ≈ 3T），平滑噪声但保留相干波
DOP 窗长度 4 × T 保证至少 4 周期稳定极化，论文建议值

本征分解与极化属性

对每个 (t,f) 做 3×3 复矩阵本征分解： $$ S = V \Lambda V^H,\quad \Lambda = \operatorname{diag}(\lambda_1,\lambda_2,\lambda_3) $$

DOP（Degree of Polarization）

$$ \mathrm{DOP} = \frac{\lambda_1 - \lambda_2}{\lambda_1 + \lambda_2 + \lambda_3} $$

Planarity

通过主、次特征向量叉积求平面法线，与竖直方向夹角 α： - 若 α ≈ 90° → 水平 planarity（Rayleigh 波特征）
- 若 α ≈ 0° → 垂直平面（Love 波或噪声）

参数 推荐值 理由
dop_thresh 0.8 论文常用 0.75–0.85；取 0.8 折中稳健
alpha_min 60° 仅保留水平 planarity

背方位角（BAZ）计算

公式

$$ \mathrm{BAZ} = \operatorname{atan2}(E, N) $$ 以度为单位，0° = 北，顺时针增加。

处理细节

  • 使用主特征向量的水平分量 (N,E)。
  • 实部幅值过小则跳过或使用相位差法。
  • 输出角度范围 [0, 360°)。

二义性处理

  • 单台存在 ±180° 二义性。
  • 假设 retrograde（瑞利波基模），可选反向修正：若竖直与水平分量相位差 ~180°，翻转 BAZ。

筛选与聚合输出

筛选规则

DOP >= 0.8
planarity_angle > 60°

脚本

DOPE