QSEIS2025：从模型到位移波形#

此教程使用仓库提供的 AK135 弹性参数、一个震源深度、一个接收深度和三个距离。默认串行计算，不需要 MPI，也不依赖本机 test/ 目录。

1. 运行位移示例#

先完成安装，进入仓库根目录并激活环境：

python examples/qseis2025.py --regional --output-dir examples/output/qseis2025-regional

Windows 非交互调用可使用：

conda run -n pygrnwang python examples/qseis2025.py --regional --output-dir examples/output/qseis2025-regional

脚本自动写出模型、预处理输入、写入 64 秒震源样本、运行 QSEIS2025、转换结果，并用 output_type="disp" 读取合成位移。示例使用 10 km 深度震源、地表接收、300/600/900 km 距离、4 秒采样间隔（Nyquist 频率 0.125 Hz），并启用平地球变换。底层格林函数库计算 4092 秒、1024 个采样点；合成后保存发震后 0–1020 秒（包含两端，共 256 点）的数组和绘图。求解器运行需要几分钟。模型取前 24 个数值行，常数 Qp=600、Qs=300 是教程设置，并不代表完整 AK135-F 衰减模型。

2. 查看结果#

输出目录中包含 disp.npz、disp.png、source_time_function.json、summary.json、模型文件与 library/。NPZ 保存 0–1020 秒的波形、时间、距离、分量标签及单位，位移数组形状为 (3, 3, 256)；摘要记录环境、运行时间、输出形状和文件大小。library/ 保留完整的 1024 点原生计算结果。

QSEIS2025 小算例在 300、600、900 公里处的三分量位移波形。

10 km 深度震源在 300、600、900 km 处的合成位移，发震后 0–1020 秒；精确参数以共用脚本为准。#

此示例的 rotate=True 对应 东、北、上 三分量。图中位移单位为米，震源矩为脚本中指定的 \(10^{15}\) N m。图的横轴使用示例设定的时间零点；改变约化时间或平移选项后，应重新核对横轴含义。

震源采用 wavelet_type=0 和 1024 个自定义矩率样本：64 秒归一化 sin² 脉冲，预先补偿 QSEIS 的数值阻尼，使实际 STF 面积为 1、质心为 32 秒，与 SPGRN2020 一致。格林函数库保存速率核，读取程序在 output_type="disp" 时积分一次。

3. 计算应变与应力#

在另一个目录中启用完整的教程输出：

python examples/qseis2025.py --regional --observables all --output-dir examples/output/qseis2025-regional-tensors

新增 strain 和 stress 的数组与图，脚本直接向读取程序请求位移、应变和应力，均保存 0–1020 秒；应变与应力数组形状为 (3, 6, 256)。启用地理旋转时，对称张量排列为 [EE, EN, EU, NN, NU, UU]，U 表示向上，也对应代码中的 Z。应变无量纲，应力单位为 Pa。不同后端未旋转张量的排列并不统一，详见科学约定。

4. 复用已完成的计算#

python examples/qseis2025.py --regional --output-dir examples/output/qseis2025-regional --reuse

复用运行会另写 summary-reuse.json，保留首次计算的摘要。仅对同一模型、网格和输出设置复用。改变计算参数时使用新目录，以免将旧库误认为新参数的结果。

模型取前 24 个数值行形成分层半空间，与完整球形地球是不同近似，不能只凭相同距离判断结果应完全相等。参见远距离教程和后端对比。

更短的入门计算#

不加 --regional 时，同一脚本计算 30/60/90 km、0.5 秒采样的轻量算例：原生库为 127.5 秒、256 点，保存 0–100 秒（201 点），几秒内完成：

python examples/qseis2025.py --output-dir examples/output/qseis2025

完整脚本#

中英文页面直接引用同一脚本：

  1"""Build QSEIS2025 introductory or regional traces; optionally include tensors."""
  2from pathlib import Path
  3
  4from common import (MECHANISM, MOMENT_NM, REGIONAL_SAMPLING_INTERVAL_S, REGIONAL_STF, finish, parser_for, prepare,
  5                    require_library_settings, save_waveforms)
  6from pygrnwang.create_qseis2025_bulk import (
  7    pre_process_qseis2025, create_grnlib_qseis2025_sequential)
  8from pygrnwang.read_qseis2025 import get_outfile_name_list, seek_qseis2025
  9from source_time_function import prepare_qseis_stf, validate_qseis_stf
 10from spectral_settings import qseis_spectral_settings
 11
 12
 13def main():
 14    parser = parser_for("qseis2025")
 15    parser.add_argument("--regional", action="store_true",
 16                        help="Use 300/600/900 km, a 64 s wavelet and Earth flattening")
 17    parser.add_argument("--point-source", action="store_true",
 18                        help="Regional control run: disable the default Gaussian spatial source smoothing")
 19    parser.set_defaults(output_dir=None)
 20    parser.add_argument("--observables", choices=("disp", "all"), default="disp",
 21                        help="all also computes strain and stress")
 22    args = parser.parse_args()
 23    if args.point_source and not args.regional:
 24        parser.error("--point-source requires --regional")
 25    source_radius_ratio = 0.0 if args.point_source else 0.05
 26    if args.output_dir is None:
 27        directory = "qseis2025-regional" if args.regional else "qseis2025"
 28        if args.point_source:
 29            directory += "-point-source"
 30        args.output_dir = Path(__file__).resolve().parent / "output" / directory
 31    output, library, model, report, started = prepare(args, "QSEIS2025")
 32    dt, window = (REGIONAL_SAMPLING_INTERVAL_S, 4092.0) if args.regional else (0.5, 127.5)
 33    output_end = 1020.0 if args.regional else 100.0
 34    output_samples = int(round(output_end / dt)) + 1  # Include the final sample.
 35    native_samples = int(round(window / dt)) + 1
 36    distances = [300.0, 600.0, 900.0] if args.regional else [30.0, 60.0, 90.0]
 37    wavelet_duration = 16 if args.regional else 4
 38    wavelet_type = 0 if args.regional else 2
 39    source_time_function = None
 40    if not args.reuse:
 41        pre_process_qseis2025(
 42            processes_num=1, path_green=library, event_depth_list=[10.0],
 43            receiver_depth_list=[0.0], dist_range=[distances[0], distances[-1]],
 44            delta_dist=distances[0],
 45            N_each_group=3, time_window=window, sampling_interval=dt, source_radius_ratio=source_radius_ratio,
 46            output_observables=([1, 0, 1, 1, 0] if args.observables == "all"
 47                                else [1, 0, 0, 0, 0]),
 48            wavelet_type=wavelet_type, wavelet_duration=wavelet_duration, time_reduction_velo=0,
 49            flat_earth_transform=args.regional, path_nd=model, earth_model_layer_num=24,
 50        )
 51        if args.regional:
 52            source_time_function = prepare_qseis_stf(
 53                library, duration_s=64.0, samples=1024)
 54        create_grnlib_qseis2025_sequential(library, remove_pd=False)
 55    require_library_settings(
 56        library, event_depth_list=[10.0], receiver_depth_list=[0.0],
 57        grn_dist_range=[distances[0], distances[-1]], grn_delta_dist=distances[0],
 58        sampling_interval=dt, time_window=window, sampling_num=native_samples,
 59        wavelet_type=wavelet_type, wavelet_duration=wavelet_duration, time_reduction_velo=0,
 60        flat_earth_transform=args.regional, earth_model_layer_num=24,
 61        slowness_window=None, wavenumber_sampling_rate=12, anti_alias=0.01,
 62        free_surface=0,
 63    )
 64    native_input = Path(library) / "10.00" / "0.00" / "0_0" / "grn.inp"
 65    lines = native_input.read_text(encoding="utf-8-sig").splitlines()
 66    headings = [i for i, line in enumerate(lines) if "WAVENUMBER INTEGRATION PARAMETERS" in line]
 67    if len(headings) != 1:
 68        raise ValueError("Expected one native wavenumber section")
 69    data = [line.split("#", 1)[0].strip() for line in lines[headings[0] + 1:]]
 70    data = [line for line in data if line]
 71    controls = [float(value) for value in data[1].split()]
 72    if controls != [1e-6, source_radius_ratio]:
 73        raise ValueError("Existing native spatial-source settings differ; use a fresh --output-dir")
 74    if args.regional and args.reuse:
 75        source_time_function = validate_qseis_stf(library)
 76    observables = ("disp", "strain", "stress") if args.observables == "all" else ("disp",)
 77    if args.reuse:
 78        # Backend metadata does not store observable flags. Check the requested
 79        # binary outputs before reading a displacement-only library as tensors.
 80        native_dir = Path(library) / "10.00" / "0.00" / "0_0"
 81        required = []
 82        for observable in observables:
 83            psv, sh = get_outfile_name_list(observable)
 84            required.extend(native_dir / ("grn_%s.bin" % name) for name in psv + sh)
 85        if any(not path.is_file() for path in required):
 86            raise ValueError("Requested outputs are absent. Recalculate in a fresh "
 87                             "--output-dir without --reuse using --observables all.")
 88    for observable in observables:
 89        # Regional type-0 kernels are rates; the reader integrates them before cropping.
 90        arrays = [MOMENT_NM * seek_qseis2025(
 91            path_green=library, event_depth_km=10.0, receiver_depth_km=0.0,
 92            az_deg=30.0, dist_km=distance, focal_mechanism=MECHANISM,
 93            srate=1 / dt, output_type=observable, rotate=True,
 94            before_p=None, shift=False, pad_zeros=False,
 95        )[:, :output_samples] for distance in distances]
 96        labels = ["E", "N", "U"] if observable == "disp" else ["EE", "EN", "EU", "NN", "NU", "UU"]
 97        unit = {"disp": "m", "strain": "1", "stress": "Pa"}[observable]
 98        save_waveforms(output, report, observable, arrays, distances, dt, labels, unit,
 99                       expected_samples=output_samples, time_limits=(0.0, output_end))
100    report.update(sampling_interval_s=dt, max_frequency_hz=0.5 / dt, time_window_s=window,
101                  output_time_range_s=[0.0, output_end], distances_km=distances,
102                  native_samples=native_samples, earth_model_numeric_rows=24,
103                  wavelet_type=wavelet_type, wavelet_duration_samples=wavelet_duration,
104                  wavelet_duration_s=wavelet_duration * dt,
105                  flat_earth_transform=args.regional, regional=args.regional, source_radius_ratio=source_radius_ratio, point_source=args.point_source)
106    if args.regional:
107        report["source_time_function"] = source_time_function
108        report["physical_source_time_function"] = dict(REGIONAL_STF)
109        report["spectral_settings"] = qseis_spectral_settings(library)
110    finish(output, report, started)
111
112
113if __name__ == "__main__":
114    main()

下一步可阅读 QSEIS2025 后端指南、读取与后处理及验证记录。小算例用于学习完整流程；正式研究还需检查空间采样、频率范围和数值收敛。