# hk

**Repository Path**: codepanda64/hk

## Basic Information

- **Project Name**: hk
- **Description**: 使用远震P波接收函数，进行HK搜索，得到莫霍面深度和k(vp/vs)
- **Primary Language**: Python
- **License**: Apache-2.0
- **Default Branch**: master
- **Homepage**: None
- **GVP Project**: No

## Statistics

- **Stars**: 1
- **Forks**: 0
- **Created**: 2023-12-21
- **Last Updated**: 2026-06-17

## Categories & Tags

**Categories**: Uncategorized

**Tags**: None

## README

# H-κ Stacking: 接收函数地壳厚度与 Vp/Vs 比值估计

对远震 P 波接收函数应用 H-κ（H-kappa）叠加搜索方法，估计台站下方的地壳厚度（H）和波速比（Vp/Vs, κ）。

---

## 1. 方法原理

### 基本概念

远震 P 波接收函数是径向分量与垂直分量反褶积得到的时间序列，主要包含 P 到 S 的转换波（Ps）及其在 Moho 面与地表之间的多次反射波（PpPs、PsPs+PpSs）。这些相位的走时相对直达 P 波的延迟时间与地壳厚度 H 和波速比 κ = Vp/Vs 之间存在确定关系：

```
t_Ps     = H · ( √(1/β² − p²) − √(1/α² − p²) )
t_PpPs   = H · ( √(1/β² − p²) + √(1/α² − p²) )
t_PsPs+PpSs = 2H · √(1/β² − p²)
```

其中 α = Vp, β = Vs, p 为射线参数。

### 叠加搜索

Zhu & Kanamori (2000) *J. Geophys. Res.* 提出的 H-κ 叠加方法在给定的 H 和 κ 范围内搜索，对每条接收函数在预测走时处取振幅值，并按相位加权叠加：

```
S(H, κ) = w₁ · r(t_Ps) + w₂ · r(t_PpPs) − w₃ · r(t_PsPs+PpSs)
```

其中 w₁, w₂, w₃ 为权重系数（w₁ + w₂ + w₃ = 1），`r(t)` 为接收函数在时间 t 处的振幅值。使 S(H, κ) 最大的 (H, κ) 即为最优估计，其周围的等值线用于估计不确定度。

**参考文献：**
> Zhu, L. & Kanamori, H. (2000). Moho depth variation in southern California from teleseismic receiver functions. *Journal of Geophysical Research: Solid Earth*, 105(B2), 2969–2980.

---

## 2. 安装

本程序使用 [uv](https://docs.astral.sh/uv/) 管理 Python 虚拟环境和依赖。

### 安装 uv

```bash
# macOS / Linux
curl -LsSf https://astral.sh/uv/install.sh | sh

# Windows (PowerShell)
powershell -c "irm https://astral.sh/uv/install.ps1 | iex"
```

### 创建虚拟环境并安装依赖

```bash
# 克隆仓库
git clone https://gitee.com/codepanda64/hk.git
cd hk

# 创建虚拟环境并安装依赖（Python 3.11）
uv sync

# 激活虚拟环境
source .venv/bin/activate      # Linux / macOS
.venv\Scripts\activate          # Windows
```

> 程序依赖：numpy, scipy, matplotlib, pandas, obspy，详见 `pyproject.toml`。

---

## 3. 使用

### 3.1 接收函数数据准备

本程序使用 **SAC 格式** 的接收函数文件，径向分量，建议数据预处理步骤：

1. 从三分量地震数据中截取远震 P 波波形（震中距 30°–90°，震级 ≥ 5.5）
2. 旋转至 Z-R-T 坐标系
3. 在时间域或频率域对径向分量做反褶积（以垂直分量为输入函数），得到 P 波接收函数
4. 人工或自动挑选信噪比高的接收函数，放入每个台站目录下的 `best/` 子目录

**目录结构示例：**

```
RF_results/
├── X1.53110/
│   └── best/
│       ├── 53110y020a.eqr
│       ├── 53110y023a.eqr
│       └── ...
├── X1.53111/
│   └── best/
│       ├── 53111y020a.eqr
│       ├── 53111y023a.eqr
│       └── ...
├── X1.53112/
│   └── best/
│       ├── 53112y014a.eqr
│       └── ...
└── X1.53114/
    └── best/
        ├── 53114y015a.eqr
        └── ...
```

接收函数文件为 SAC 二进制格式，需包含以下头段变量：

| 头段 | 含义 |
|---|---|
| `gcarc` | 震中距（度） |
| `dist` | 震中距（km） |
| `evdp` | 震源深度（km） |
| `b` | 文件起始时间（相对参考时刻, s） |
| `delta` | 采样间隔（s） |
| `npts` | 数据点数 |

### 3.2 配置参数

编辑 `hk.ini` 文件：

```ini
[input]
rf_path = ./RF_results/
stations_name = *
best_dir_name = best
rf_ext = eqr

[hk]
model_name = iasp91
hmin = 30
hmax = 60
kmin = 1.5
kmax = 2.5
num = 200
vp = 6.2
pre_time = 0
post_time = 40
w1 = 0.7
w2 = 0.2
w3 = 0.1

[output]
out_root = ./HK_results/
hK_results_file_name = {station_name}_{model_name}_{hmin}-{hmax}_{kmin}-{kmax}_{num}_{vp}.txt
hk_img_name = {station_name}_{model_name}_{hmin}-{hmax}_{kmin}-{kmax}_{num}_{vp}.png
```

**主要参数说明：**

| 参数 | 含义 | 建议值 |
|---|---|---|
| `stations_name` | 台站代码，多个用逗号分隔，`*` 表示处理所有 | `*` 或 `X1.53110,X1.53111` |
| `model_name` | TauP 走时模型 | `iasp91` |
| `hmin, hmax` | 莫霍面深度搜索范围（km） | 30–60 |
| `kmin, kmax` | Vp/Vs 搜索范围 | 1.5–2.5 |
| `num` | 搜索网格数 | 200 |
| `vp` | 地壳 P 波平均速度（km/s） | 6.2 |
| `pre_time, post_time` | 接收函数分析时间窗（s） | 0–40 |
| `w1, w2, w3` | 各相位权重 | 0.7, 0.2, 0.1 |

### 3.3 运行

```bash
python hk.py
```

### 3.4 输出

#### 文本结果

输出目录 `HK_results/` 下每个台站生成一个文本文件，格式为：

```
台站名, best_H, best_H_error, best_kappa, best_kappa_error
```

示例（`X1.53111_iasp91_30-60_1.5-2.5_30_6.2.txt`）：

```
X1.53111, 32.07, 2.20, 1.81, 0.10
```

本目录 4 个台站的处理结果汇总：

| 台站 | 接收函数数 | 莫霍面深度 H (km) | Vp/Vs (κ) |
|---|---|---|---|
| X1.53110 | 51 | 34.14 ± 2.38 | 1.707 ± 0.080 |
| X1.53111 | 37 | 32.07 ± 2.20 | 1.810 ± 0.102 |
| X1.53112 | 57 | 34.14 ± 2.56 | 1.741 ± 0.084 |
| X1.53114 | 40 | 36.21 ± 3.10 | 1.776 ± 0.089 |

#### 图形输出

每个台站生成一张 2×2 的 H-κ 图像（`*.png`），包含：
- 左上：Ps 相叠加结果
- 右上：PpPs 相叠加结果
- 左下：PsPs+PpSs 相叠加结果
- 右下：加权叠加结果（H-κ 图），标记最优解及误差等值线

---

## 4. 目录结构

```
hk/
├── hk.py              # H-κ 叠加搜索主程序
├── main.py            # 入口占位
├── hk.ini             # 配置文件
├── pyproject.toml     # 项目依赖
├── uv.lock            # uv 锁定文件
├── environment.yml    # conda 环境（备选）
├── .python-version    # Python 版本
├── README.md          # 本文件
├── RF_results/        # 接收函数数据目录
│   ├── X1.53110/best/
│   ├── X1.53111/best/
│   ├── X1.53112/best/
│   └── X1.53114/best/
└── HK_results/        # 结果输出目录
    ├── *.txt
    └── *.png
```