【程序分享】三维到一维速度模型转换程序

Python >= 3.7numpy >= 1.19.0pandas >= 1.2.0scipy >= 1.6.0matplotlib >= 3.3.0

# 下载程序文件wget https://github.com/zhangsuxiang/3D-to-1D-Velocity-Model-Converter
# 安装依赖包pip install numpy pandas scipy matplotlib
# 可选：在虚拟环境中安装python -m venv velo_envsource velo_env/bin/activate  # Windows系统: velo_env\Scripts\activatepip install numpy pandas scipy matplotlib

# 测试导入python -c "from velocity_3d_to_1d import Velocity3Dto1D; print('安装成功！')"

from velocity_3d_to_1d import Velocity3Dto1D
# 加载三维速度模型converter = Velocity3Dto1D('velocity_3d.txt', velocity_type='Vp')
# 提取一维模型depths, velocities = converter.extract_1d_velocity(method='weighted_mean')
# 插值到1km间隔depths_new, vel_new = converter.interpolate_1d_model(depth_interval=1.0)
# 保存结果converter.save_1d_model('velocity_1d.dat', depths_new, vel_new, format='hypo71')

经度  纬度  深度(km)  速度(km/s或比值)

102.5  25.3  0.0   5.80102.5  25.3  5.0   6.10102.5  25.3  10.0  6.35...

# 自动检测分隔符converter = Velocity3Dto1D('data.txt', velocity_type='Vp')
# 指定速度类型converter = Velocity3Dto1D('data.csv', velocity_type='Vs')

import numpy as np
# 创建或加载形状为(n, 4)的数组data_array = np.loadtxt('data.txt')converter = Velocity3Dto1D(data_array=data_array, velocity_type='VpVs')

depths, velocities = converter.extract_1d_velocity(method='mean')

# 按纬度加权，考虑地球球面效应depths, velocities = converter.extract_1d_velocity(    method='weighted_mean',    weight_by_area=True)

# 从特定区域提取depths, velocities = converter.extract_1d_velocity(    method='weighted_mean',    region=(经度最小值, 经度最大值, 纬度最小值, 纬度最大值))

# 线性插值（适用于有明显界面的模型）depths_linear, vel_linear = converter.interpolate_1d_model(    depth_interval=0.5,    kind='linear')
# 三次样条插值（适用于平滑模型）depths_cubic, vel_cubic = converter.interpolate_1d_model(    depth_interval=1.0,    kind='cubic')
# Akima插值（保形插值，防止过冲）depths_akima, vel_akima = converter.interpolate_1d_model(    depth_interval=2.0,    kind='akima')

# 执行质量检查qc_report = converter.quality_check(depths, velocities)
# 查看报告详情print(f"状态: {qc_report['status']}")print(f"警告: {qc_report['warnings']}")print(f"信息: {qc_report['info']}")

# HYPO71 velocity model - Vp# Depth(km)  Vp(km/s)   0.00   5.800   1.00   5.850   2.00   5.900

* HYPOINVERSE VELOCITY MODEL - Vp  5.80    0.00  5.85    1.00  5.90    2.00

VELEST velocity model - Vp   3  5.800  5.850  5.900  0.00  1.00  2.00

# 1D Velocity Model - Vp# Depth(km)  Vp(km/s)0.0005.8001.0005.8502.0005.900

import numpy as npimport matplotlib.pyplot as pltfrom velocity_3d_to_1d import Velocity3Dto1D
# 步骤1：加载P波速度模型print("加载三维P波速度模型...")converter_vp = Velocity3Dto1D('vp_model_3d.txt', velocity_type='Vp')
# 步骤2：定义感兴趣区域region = (100.0, 105.0, 20.0, 25.0)  # 经度最小, 经度最大, 纬度最小, 纬度最大
# 步骤3：提取一维模型print("提取一维模型...")depths_vp, velocities_vp = converter_vp.extract_1d_velocity(    method='weighted_mean',    region=region,    weight_by_area=True)
# 步骤4：插值到所需分辨率print("插值模型...")depths_interp, vel_interp = converter_vp.interpolate_1d_model(    depth_interval=0.5,  # 500m间隔    depth_max=100.0,     # 最大深度100km    kind='cubic')
# 步骤5：质量检查print("执行质量检查...")qc_report = converter_vp.quality_check(depths_interp, vel_interp)
if qc_report['status'] == 'PASS':    print(" 模型通过所有质量检查")else:    print(" 模型存在警告:")    for warning in qc_report['warnings']:        print(f"  - {warning}")
# 步骤6：保存多种格式print("保存模型...")converter_vp.save_1d_model('vp_hypo71.dat', depths_interp, vel_interp, format='hypo71')converter_vp.save_1d_model('vp_velest.dat', depths_interp, vel_interp, format='velest')
# 步骤7：可视化结果print("创建可视化图表...")fig = converter_vp.plot_comparison(    (depths_vp, velocities_vp),    (depths_interp, vel_interp),    save_figure='vp_model_comparison.png')plt.show()
# 步骤8：生成统计报告converter_vp.plot_statistics(save_figure='vp_statistics.png')

# 加载Vp和Vs模型converter_vp = Velocity3Dto1D('vp_3d.txt', velocity_type='Vp')converter_vs = Velocity3Dto1D('vs_3d.txt', velocity_type='Vs')
# 使用相同参数提取一维模型depths_vp, vel_vp = converter_vp.extract_1d_velocity(method='weighted_mean')depths_vs, vel_vs = converter_vs.extract_1d_velocity(method='weighted_mean')
# 计算Vp/Vs比值vpvs_ratio = vel_vp / vel_vs
# 创建Vp/Vs模型vpvs_data = np.column_stack([    np.full(len(depths_vp), 105.0),  # 虚拟经度    np.full(len(depths_vp), 25.0),   # 虚拟纬度    depths_vp,    vpvs_ratio])
converter_vpvs = Velocity3Dto1D(data_array=vpvs_data, velocity_type='VpVs')
# Vp/Vs质量检查qc_vpvs = converter_vpvs.quality_check(depths_vp, vpvs_ratio)

# 定义多个区域regions = {    '北部': (100, 102, 25, 27),    '南部': (100, 102, 23, 25),    '东部': (102, 104, 24, 26),    '西部': (98, 100, 24, 26)}
# 处理每个区域results = {}for region_name, bounds in regions.items():    print(f"处理区域: {region_name}")
    depths, velocities = converter.extract_1d_velocity(        method='weighted_mean',        region=bounds    )
    # 插值    depths_i, vel_i = converter.interpolate_1d_model(depth_interval=1.0)
    # 保存    converter.save_1d_model(        f'velocity_1d_{region_name}.dat',        depths_i, vel_i    )
    results[region_name] = (depths_i, vel_i)
# 对比各区域plt.figure(figsize=(10, 8))for region_name, (depths, vels) in results.items():    plt.plot(vels, depths, label=region_name)plt.xlabel('速度 (km/s)')plt.ylabel('深度 (km)')plt.title('区域一维速度模型对比')plt.legend()plt.gca().invert_yaxis()plt.grid(True, alpha=0.3)plt.savefig('区域对比.png')