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

当然是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 $$

参数与理由

说明：
σ = 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 $$

本征分解与极化属性

对每个 (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 波或噪声）

背方位角（BAZ）计算

公式

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

处理细节

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

二义性处理

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

筛选与聚合输出

筛选规则

DOP >= 0.8
planarity_angle > 60°

脚本

DOPE