SEISAMUSE

import os
import h5py
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from datetime import datetime, timedelta
from matplotlib.dates import DateFormatter


def parse_filename_to_start_time(filename):
    """从文件名推断数据起始时间"""
    date_str = filename.split('_')[0]  # 提取"20130512"
    return datetime.strptime(date_str, "%Y%m%d")  # 转为 datetime


def load_cross_spec_from_files(hdf5_dir, start_date, end_date):
    """读取目录下所有落在给定时间范围的互相关谱"""
    files = sorted([f for f in os.listdir(hdf5_dir) if f.endswith("_cross_spec.hdf5")])
    #files = sorted([f for f in os.listdir(hdf5_dir) if f.endswith("_power_spec.hdf5")])
    time_list = []
    spec_list = []
    for file in files:
        try:
            base_time = parse_filename_to_start_time(file)
        except Exception as e:
            print(f"跳过无法解析时间的文件：{file}")
            continue

        if not (start_date <= base_time <= end_date):
            continue

        file_path = os.path.join(hdf5_dir, file)
        with h5py.File(file_path, 'r') as f:
            for i in range(1,24):
                ds_name = f'window_{i}'
                if ds_name in f:
                    data = f[ds_name][:]
                    current_time = base_time + timedelta(hours=i)
                    time_list.append(current_time)
                    spec_list.append(data)
    return time_list, spec_list


def plot_spectrogram(time_list, spec_list, freqs=None, title="Cross Spectrogram", save_path=None):
    """绘制时间-频率图像"""
    if not time_list or not spec_list:
        print("没有数据可供绘图。")
        return

    spec_array = np.array(spec_list).T  # shape: freq x time
    #spec_array = spec_array - np.mean(spec_array, axis=1, keepdims=True)
    n_freq = spec_array.shape[0]

    if freqs is None:
        freqs = np.linspace(0.005, 0.1, n_freq)

    fig, ax = plt.subplots(figsize=(14, 6))
    time_nums = mdates.date2num(time_list)

    cax = ax.pcolormesh(time_nums, freqs, spec_array, shading='auto', cmap='inferno')
    #cax = ax.pcolormesh(time_nums, freqs, spec_array, shading='auto' )
    fig.colorbar(cax, label='Cross Spectral Amplitude')
    ax.set_title(title)
    ax.set_yscale('log')  # 设置纵轴为对数刻度
    ax.set_xlabel("Time (UTC)")
    ax.set_ylabel("Frequency (Hz)")
    ax.xaxis.set_major_formatter(DateFormatter("%Y-%m-%d\n%H:%M"))
    plt.xticks(rotation=30)
    plt.tight_layout()

    if save_path:
        plt.savefig(save_path, dpi=300)
        print(f"图像已保存到：{save_path}")

    plt.show()

# 设置参数
hdf5_dir = "output_1_0"  # HDF5 文件夹路径
start_date = datetime(2006, 5, 1)
end_date = datetime(2006, 5,31)
# 加载数据
time_list, spec_list = load_cross_spec_from_files(hdf5_dir, start_date, end_date)
# 绘图并保存为 PNG（可选）
plot_spectrogram(time_list, spec_list,
                 title=f"Cross Spectral Density ({start_date.date()} ~ {end_date.date()})",
                 save_path="cross_spectrogram_2006_1_0.png")

学习匹配场

发表于 2025-06-27 本文字数： 2.8k 阅读时长 ≈ 3 分钟

今天学习匹配场处理（Matched Field Processing, MFP），见这里。

阅读全文 »

有趣的jupyter-notebook插件

发表于 2025-06-21 分类于 python 本文字数： 675 阅读时长 ≈ 1 分钟

有趣的jupyter notebook插件:

插件名称	简介	官方链接
Jupyter Contrib NBe	核心扩展包，集成50+插件（代码补全、变量监控、执行时间显示等）	GitHub
Hinterland	实时代码补全工具，支持 Python/R/Julia，减少拼写错误	文档
Table of Contents (2)	动态生成目录，支持标题跳转，适合长文档导航	配置页
Variable Inspector	侧边栏实时显示变量类型/大小/值，替代频繁 `print(type)` 操作	说明
ExecuteTime	记录每个 Cell 的执行时间和完成时间，优化性能分析	详情
Autopep8	一键格式化代码为 PEP8 标准，提升可读性	PyPI
jupyterthemes	界面主题美化，支持暗黑模式/护眼配色（如 `monokai`, `solarized`）	GitHub
Notify	内核空闲时发送浏览器通知，适合长时间任务（如模型训练）	文档
jupyter Widgets	创建交互控件（滑块/下拉菜单），将静态图表转为动态仪表盘	官网
Voilà	将 Notebook 转为独立 Web 应用，隐藏代码仅保留交互结果	文档
RISE	实时代码幻灯片工具，用 Markdown 标题分页，支持演示中修改参数	GitHub

安装：

1
2
3

pip install jupyter_contrib_nbextensions jupyterthemes voila rise
jupyter contrib nbextension install
jupyter nbextension enable [插件名]  # 启用具体插件

有趣的新闻网站

发表于 2025-06-21 分类于杂本文字数： 2.1k 阅读时长 ≈ 2 分钟

以下是有趣的新闻网站：

AllYouCanRead是一个综合性资讯导航平台，按地域聚合全球六大洲新闻站点，同时以艺术、商业、科技等超 30 个主题精细分类各类网站资源，杂志专区更细分动物、时尚、育儿等数十个品类，通过清晰的结构化导航，为用户提供一站式全球新闻媒体与实用网站的定向检索服务。
Apple News是苹果公司为iOS和macOS设备打造的官方新闻应用，它不仅整合主流新闻媒体的免费内容，还通过订阅服务 Apple News+ 提供大量杂志和深度报道。该平台以用户界面友好、内容精选优质而著称，尤其适合苹果生态中的用户日常阅读使用。借助苹果的算法推荐和编辑精选，用户可以在一个界面中阅读来自多家出版商的新闻、专题和评论。对于追求高品质阅读体验的用户，Apple News 提供了一种简洁但内容丰富的获取方式。
BBC News是英国广播公司旗下的新闻平台，长期以客观、权威的报道风格受到全球用户信赖。它提供来自全球各地的实时新闻、深度分析、专题报道以及视频广播内容，尤其擅长以多语种进行国际传播，包括提供中文页面服务。BBC News 无需订阅即可访问，内容涵盖政治、科技、文化、健康等广泛领域，适合希望获取权威信息、提升英文阅读能力或关注全球动态的广大读者。
CNN作为国际知名新闻平台，以实时追踪全球突发新闻为核心特色，24 小时不间断更新政治、经济、科技、文化及娱乐等领域资讯。其报道兼具速度与深度，既提供俄乌冲突、美国大选等热点事件的现场直击，也通过专家分析、纪录片等形式挖掘新闻背后的背景逻辑。
Flipboard是一个以“杂志式排版”为特色的新闻聚合与社交内容平台，用户可以根据自己的兴趣订阅特定主题、媒体或用户发布的内容源，并将喜欢的内容整理成属于自己的“杂志”。它强调视觉体验和内容个性化，非常适合喜欢图文并茂、轻松滑阅的移动端用户。Flipboard 支持离线阅读和多平台同步，并通过滑动式界面增强交互感，使其成为许多用户在通勤、碎片时间中快速获取资讯的理想选择。
Google News是谷歌推出的一项新闻聚合服务，它通过智能算法从全球各大新闻网站实时抓取内容，为用户推送个性化、主题化的新闻报道。平台聚合同一事件来自不同媒体的视角，帮助读者全面了解事件全貌。其界面简洁、使用免费，同时支持桌面和移动端，是获取国际主流新闻、快速了解时事动态的便捷工具。用户还可以设置感兴趣的主题或地区，实现高度定制化的阅读体验，非常适合需要快速掌握新闻全局的用户。
Newsnow作为一个英国的多元新闻聚合平台，具备诸多显著优势。其新闻内容丰富多元，广泛涵盖英国、世界、商业、娱乐、体育及科技等各个领域，能够满足不同用户的多样化阅读需求。网站的界面设计简洁直观，分类清晰，便于用户快速定位感兴趣的新闻板块。其个性化推荐功能有助于用户发现更多符合个人喜好的内容。
NewsBrief核心优势在于以分钟级速度聚合全球60+语言信源，通过AI驱动的动态聚类技术，将碎片化新闻提炼为可视化事件脉络，其独家”EU Focus”模块深度追踪28个欧盟机构政策动向，并支持跨事件对比分析，为政策研究者、商业机构提供学术级舆情洞察。
Newspaper Map通过地图可视化交互，直观聚合全球报纸资源，支持按地理位置一键直达各地新闻源，突破传统新闻检索的国界与语言限制。
Paperboy在算法推荐新闻泛滥的时代，反其道而行，成为”数字报亭式”的古典新闻聚合平台，满足用户对传统报纸阅读体验的怀旧需求，同时解决地方性信息获取痛点。
PressReader是一个汇集全球报纸和杂志的数字阅读平台，提供来自120多个国家、70多种语言的7000多份出版物原版内容。其最大优势在于内容全面、更新及时，涵盖新闻、财经、科技、时尚等多个领域。平台呈现与实体刊物一致的排版，支持全文搜索、多语种翻译、音频朗读与剪藏分享，提升了阅读的互动性和便捷性。
The New York Times是全球最具影响力的新闻机构之一，提供高质量的新闻报道、评论、深度专题和调查类文章。作为一个付费订阅平台，它以专业记者团队、深入独家的调查报道和精准的国际视角闻名，尤其在政治、经济、国际关系等领域具有广泛权威性。无论是通过网页还是手机 App，纽约时报为读者提供一种深入理解世界的方式，是重度新闻使用者、研究者和政策关注者的重要信息来源。
World News提供全球新闻资讯，整合大量知名媒体的报道，涵盖政治、经济、社会、文化及自然灾害等丰富多样的新闻内容，更新及时，能够快速反映世界各地的最新动态。该网站界面简洁直观，以图文列表形式展示新闻，每条新闻都明确标注了来源媒体和发布时间，增强了新闻的可信度与透明度。
环球网国际新闻以 “全球视野、中国视角” 为核心，兼具时效性、权威性与多元覆盖，既关注国际政治经济格局的宏观变化，也聚焦与中国及华人相关的具体议题，通过结构化板块和图文结合的方式，为读者提供全面、及时的国际资讯。
央视网作为国家级新媒体旗舰平台，既保持主流媒体的权威性，又通过技术创新实现“新闻+政务+服务”深度融合，是媒体融合发展的标杆案例。其核心优势在于依托总台内容资源库构建的全媒体传播体系，以及在无障碍访问等社会责任领域的先行实践。

$\LaTeX$的数学公式表达

发表于 2025-06-21 分类于 web 本文字数： 2.2k 阅读时长 ≈ 2 分钟

以下是$\LaTeX$的数学算符

说明：公式需包裹在 $...$ （行内公式）或 $$...$$（块级公式）中，以下表格中的 渲染效果 需在支持 LaTeX 的 Markdown 环境中显示（如 Typora、Obsidian 等）。

一、基础运算符号

符号名称	$\LaTeX$ 命令	渲染效果
加号	`a + b`	$a + b$
减号	`a - b`	$a - b$
乘号（叉乘）	`a \times b`	$a \times b$
乘号（点乘）	`a \cdot b`	$a \cdot b$
除号	`a \div b`	$a \div b$
加减号	`a \pm b`	$a \pm b$
减加号	`a \mp b`	$a \mp b$

二、关系运算符

符号名称	$\LaTeX$ 命令	渲染效果
等于	`a = b`	$a = b$
不等于	`a \neq b`	$a \neq b$
约等于	`a \approx b`	$a \approx b$
大于等于	`a \geq b`	$a \geq b$
小于等于	`a \leq b`	$a \leq b$
远大于	`a \gg b`	$a \gg b$
远小于	`a \ll b`	$a \ll b$
正比于	`a \propto b`	$a \propto b$

三、集合运算符

符号名称	$\LaTeX$ 命令	渲染效果
并集	`A \cup B`	$A \cup B$
交集	`A \cap B`	$A \cap B$
属于	`x \in A`	$x \in A$
不属于	`x \notin B`	$x \notin B$
子集	`A \subset B`	$A \subset B$
真子集	`A \subseteq B`	$A \subseteq B$
空集	`\emptyset`	$\emptyset$
实数集	`\mathbb{R}`	$\mathbb{R}$
自然数集	`\mathbb{N}`	$\mathbb{N}$

四、微积分符号

符号名称	$\LaTeX$ 命令	渲染效果
积分	`\int_{a}^{b} f(x) dx`	$\int_{a}^{b} f(x) dx$
偏导数	`\frac{\partial f}{\partial x}`	$\frac{\partial f}{\partial x}$
极限	`\lim_{x \to 0} \frac{\sin x}{x}`	$\lim_{x \to 0} \frac{\sin x}{x}$
求和	`\sum_{i=1}^{n} i^2`	$\sum_{i=1}^{n} i^2$
导数（撇号形式）	`f'(x)`	$f’(x)$
梯度	`\nabla f`	$\nabla f$
二阶导数	`\frac{d^2 y}{dx^2}`	$\frac{d^2 y}{dx^2}$

五、希腊字母

小写字母	$\LaTeX$ 命令	渲染效果	大写字母	LaTeX 命令	渲染效果
α (alpha)	`\alpha`	$\alpha$	Γ (Gamma)	`\Gamma`	$\Gamma$
β (beta)	`\beta`	$\beta$	Δ (Delta)	`\Delta`	$\Delta$
θ (theta)	`\theta`	$\theta$	Θ (Theta)	`\Theta`	$\Theta$
π (pi)	`\pi`	$\pi$	Π (Pi)	`\Pi`	$\Pi$
σ (sigma)	`\sigma`	$\sigma$	Σ (Sigma)	`\Sigma`	$\Sigma$

六、箭头符号

符号名称	$\LaTeX$ 命令	渲染效果
右箭头	`\rightarrow`	$\rightarrow$
左箭头	`\leftarrow`	$\leftarrow$
双向箭头	`\leftrightarrow`	$\leftrightarrow$
蕴含符号	`\Rightarrow`	$\Rightarrow$
等价符号	`\Leftrightarrow`	$\Leftrightarrow$
映射箭头	`\mapsto`	$\mapsto$

七、括号与定界符

符号名称	$\LaTeX$ 命令	渲染效果
圆括号（自适应）	`\left( \frac{a}{b} \right)`	$\left( \frac{a}{b} \right)$
方括号	`\left[ x \right]`	$\left[ x \right]$
花括号	`\left\{ x \right\}`	${ x }$
绝对值	`\lvert x \rvert`	$\lvert x \rvert$
范数	`\lVert \mathbf{v} \rVert`	$\lVert \mathbf{v} \rVert$

八、矩阵环境

$$ 
\begin{pmatrix}  % 圆括号矩阵
a & b \\
c & d 
\end{pmatrix}
\quad
\begin{bmatrix}  % 方括号矩阵
a & b \\
c & d 
\end{bmatrix}
\quad
\begin{vmatrix}  % 行列式
a & b \\
c & d 
\end{vmatrix}
$$

渲染效果：
$$
\begin{pmatrix} % 圆括号矩阵
a & b \
c & d
\end{pmatrix}
\quad
\begin{bmatrix} % 方括号矩阵
a & b \
c & d
\end{bmatrix}
\quad
\begin{vmatrix} % 行列式
a & b \
c & d
\end{vmatrix}
$$

看样子花括号和矩阵都不太行啊。

PYTHON脚本练习（一）

发表于 2025-06-21 更新于 2025-06-30 分类于 work 本文字数： 4.7k 阅读时长 ≈ 4 分钟

以下是python脚本练习1，功能包括：

遍历目录events_20250619下所有子目录中以bhz.SAC_rm结尾的SAC文件；
对这些数据进行窄带滤波，宽度为中心频率（周期分之一）的$\pm$5mHz，滤波器为4个极点0相位的Butterworth，滤波周期为arange(25,145,10)；
计算窄带滤波后的每个周期的信噪比。信噪比定义为信号窗口内，波形包络的最大值比上噪声窗口的均方根。信号窗口定义为2.5-5km/s的到时。噪声窗定义为信号末端之后1000秒开始的1000秒长度的窗口。计算的信噪比写入到user0；
将处理后的数据写到新的文件夹bp_sac中，文件名命名为z.year.jday.00.STA.bhz.period，仅保留信噪比大于3的数据。
采用并行处理(8个cpu)。
统计每个周期信噪比大于3的波形数据。
统计每个周期信噪比大于3的波形的平均信噪比。

将统计结构写入csv，并画出统计结果。

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from obspy import read
from obspy.signal.filter import envelope
from concurrent.futures import ThreadPoolExecutor, as_completed
from collections import defaultdict

# 参数设置 
input_root = "events_20250619"
output_root = "bp_sac"
center_periods = np.arange(25, 150, 10)
bandwidth = 0.010
vmin, vmax = 2.5, 5.0
noise_offset, noise_len = 1000, 1000
snr_threshold = 3.0
max_workers = 8  # 根据 CPU 核心数调整

# 存储每个周期的 SNR 值
period_snr_map = defaultdict(list)

def process_file(filepath, root):
    results = []
    try:
        st = read(filepath)
        tr = st[0]
        sac = tr.stats.sac
    except Exception as e:
        print(f"跳过 {filepath}: {e}")
        return results

    if not hasattr(sac, "o") or not hasattr(sac, "dist"):
        return results

    dist = sac.dist
    starttime = tr.stats.starttime
    win_start = starttime + dist / vmax
    win_end = starttime + dist / vmin
    t = tr.times()
    abs_time = np.array([starttime + float(tt) for tt in t])

    for period in center_periods:
        fc = 1.0 / period
        fmin, fmax = fc - bandwidth / 2, fc + bandwidth / 2
        if fmin <= 0:
            continue
        tr_filt = tr.copy()
        tr_filt.detrend("demean")
        tr_filt.taper(0.05)
        tr_filt.filter("bandpass", freqmin=fmin, freqmax=fmax, corners=4, zerophase=True)
        env = envelope(tr_filt.data)
        idx_sig = np.where((abs_time >= win_start) & (abs_time <= win_end))[0]
        if len(idx_sig) == 0:
            continue
        signal_max = np.max(env[idx_sig])
        noise_start = win_end + noise_offset
        noise_end = noise_start + noise_len
        idx_noise = np.where((abs_time >= noise_start) & (abs_time <= noise_end))[0]
        if len(idx_noise) == 0:
            continue
        noise_rms = np.sqrt(np.mean(env[idx_noise] ** 2))
        snr = signal_max / noise_rms if noise_rms > 0 else 0

        if snr >= snr_threshold:
            tr_filt.stats.sac.user0 = snr
            year = tr.stats.starttime.year
            jday = tr.stats.starttime.julday
            station = tr.stats.station.lower()
            channel = tr.stats.channel.lower()
            outname = f"z.{year}.{jday:03d}.00.{station}.{channel}.{period:03d}"

            rel_dir = os.path.relpath(root, input_root)
            out_dir = os.path.join(output_root, rel_dir)
            os.makedirs(out_dir, exist_ok=True)
            outpath = os.path.join(out_dir, outname)
            tr_filt.write(outpath, format="SAC")

            results.append((period, snr))
            print(f"✔ 保存: {outname}, SNR={snr:.2f}")

    return results
# ===== 收集所有文件路径 =====
all_files = []
for root, dirs, files in os.walk(input_root):
    for file in files:
        if file.endswith("bhz.SAC_rm"):
            all_files.append((os.path.join(root, file), root))

print(f"📁 待处理文件数: {len(all_files)}")

# ===== 并行处理 =====
with ThreadPoolExecutor(max_workers=max_workers) as executor:
    future_to_file = {executor.submit(process_file, fpath, root): (fpath, root) for fpath, root in all_files}
    for future in as_completed(future_to_file):
        try:
            result = future.result()
            for period, snr in result:
                period_snr_map[period].append(snr)
        except Exception as e:
            fpath, _ = future_to_file[future]
            print(f"❌ 文件出错 {fpath}: {e}")

# ===== 统计与可视化 =====
periods = sorted(period_snr_map.keys())
counts = [len(period_snr_map[p]) for p in periods]
means = [np.mean(period_snr_map[p]) if len(period_snr_map[p]) > 0 else 0 for p in periods]

df = pd.DataFrame({
    "Period(s)": periods,
    "Count(SNR>3)": counts,
    "Mean_SNR(SNR>3)": means
})
df.to_csv("snr_stats.csv", index=False)

# === 可视化 ===
plt.figure(figsize=(10, 6))
plt.bar(periods, counts, width=4, color='skyblue', edgecolor='black')
plt.xlabel("Period (s)")
plt.ylabel("Count of SNR > 3")
plt.title("Number of Traces with SNR > 3 per Period")
plt.grid(True, linestyle="--", alpha=0.5)
plt.tight_layout()
plt.savefig("snr_count_bar.png", dpi=150)

plt.figure(figsize=(10, 6))
plt.plot(periods, means, marker='o', linestyle='-', color='orange')
plt.xlabel("Period (s)")
plt.ylabel("Mean SNR (SNR > 3)")
plt.title("Mean SNR of Traces with SNR > 3 per Period")
plt.grid(True, linestyle="--", alpha=0.5)
plt.tight_layout()
plt.savefig("snr_mean_line.png", dpi=150)

print("🎉 并行处理完成，统计结果写入 snr_stats.csv")