尧图网站建设 
简介本资源是一套面向地球物理建模与反演研究者的Python实践教程聚焦pysit地震成像工具的正演模拟、反演优化及MATLAB数据互通适用于具备基础Python和数值计算能力的科研人员与高年级本科生。压缩包共305个文件以179个Python脚本含模型构建、源接收器配置、正反演核心逻辑为主体辅以58个reStructuredText文档提供API说明与流程注解、14张可视化结果图如波场快照、模型更新对比及少量C/C头文件支撑底层求解器调用整体体积仅1.8MB轻量易部署。已有294人学习下载资源结构清晰包含完整可运行示例如pysit-example-master目录、详细配置文件pysit.cfg、make.bat及部署手册.docx覆盖从环境安装、参数设置、迭代反演到.mat格式结果导出的全流程特别适合开展地震成像算法验证、跨平台Python-MATLAB联合分析或课程实验复现。1. 用 PySIT 做地震波正演模拟、反演成像并把模型/数据存成.mat文件这不是 MATLAB 专属流程你手头有一套地下介质参数比如速度场想看看地震波在其中怎么传播——这是正演反过来你有一组实际采集的地震记录想反推地下结构长什么样——这是反演。传统上这类计算常被默认绑定在 MATLAB 生态里尤其当结果要交给地质解释人员时.mat文件几乎是交付标配。但现实是团队里 Python 工程师越来越多MATLAB 许可成本高、部署难、CI/CD 集成卡顿。PySIT 正是为这种场景设计的——它用纯 Python 实现了有限差分波场模拟、全波形反演FWI、最小二乘反演等核心算法且原生支持 NumPy天然兼容 SciPy、Matplotlib更重要的是它不依赖 MATLAB 运行时却能无缝对接.mat文件的读写。本文面向已掌握 Python 基础NumPy/Pandas、了解波动方程基本概念的地球物理建模者或算法工程师不讲泛泛而谈的“Python 多好”只聚焦如何用 PySIT 跑通一个完整正反演闭环并确保中间模型、观测数据、反演结果全部按地质行业通行标准存为.mat文件——不是用scipy.io.savemat粗暴打包而是结构清晰、字段可读、MATLAB 端开箱即用。2. 安装 PySIT 及配套生态避开 conda-forge 的旧版本陷阱用源码编译保障反演稳定性PySIT 并未发布到 PyPI 主索引官方推荐安装方式是通过 GitHub 源码构建。直接pip install pysit会失败而conda install -c conda-forge pysit虽然能装上但截至 2024 年中conda-forge 通道中的最新版仍是 0.4.0发布于 2021 年缺失对scipy1.10的适配且反演模块在多线程下存在内存泄漏风险。生产环境必须使用当前主干main branch代码。2.1 依赖预检与环境隔离PySIT 重度依赖scipy用于稀疏矩阵求解器、numba加速波场传播内核、h5py可选用于大体积数据存储。注意numba必须与llvmlite版本严格匹配否则编译失败。建议新建独立虚拟环境python -m venv pysit_env source pysit_env/bin/activate # Linux/macOS # pysit_env\Scripts\activate.bat # Windows提示不要用--system-site-packages。PySIT 对scipy.sparse.linalg的调用路径敏感混用系统级 scipy 极易触发AttributeError: module object has no attribute cg类错误。2.2 源码编译安装含关键补丁执行以下命令拉取最新代码并安装git clone https://github.com/pysit/pysit.git cd pysit # 应用社区维护的关键补丁修复 FWI 中梯度归一化失效问题 curl -sL https://patch-diff.githubusercontent.com/raw/pysit/pysit/pull/187.patch | git apply pip install -e .[dev]该命令中-e表示开发模式安装后续修改源码可即时生效[dev]安装额外测试依赖如pytest。验证安装是否成功import pysit print(pysit.__version__) # 输出应为类似 0.5.0.dev0 from pysit import DevitoDomain print(PySIT 安装就绪)若报错ModuleNotFoundError: No module named devito说明未启用 Devito 后端PySIT 支持多种求解器后端。此时需单独安装 Devito仅当需要高阶精度或复杂边界条件时才必需pip install devito4.7.3 # 固定版本避免与 PySIT 0.5.x 不兼容2.3.mat文件支持scipy.io是唯一可靠选择PySIT 本身不处理文件 I/O.mat读写完全交由scipy.io。注意scipy1.9.0开始默认使用 MATLAB v7.3 格式HDF5 封装该格式支持大于 2GB 的数组且 MATLAB R2016b 全面兼容。但旧版scipy1.8.0默认生成 v7 格式无法保存超过 2GB 的三维速度模型。务必检查版本python -c import scipy; print(scipy.__version__) # 必须 ≥ 1.9.0若版本过低升级pip install --upgrade scipy注意scipy.io.savemat生成的.mat文件在 MATLAB 中加载后变量名即为字典 key。例如savemat(model.mat, {vel: vel_array, mesh: mesh_dict})MATLAB 中load(model.mat)后直接可用vel和mesh变量——这是地质解释流程中最友好的交互方式。3. 构建正演模型从网格定义、源接收器布设到波场快照保存正演是反演的基础。PySIT 使用pysit.modeling.Modeling类封装波场模拟逻辑。本节以二维声波方程为例构建一个含 3 层介质的合成模型并生成 100 个时间步的波场快照。3.1 定义计算域与离散化参数import numpy as np from pysit import ConstantDensityAcousticWave, CartesianDomain, PointSource, PointReceiver from pysit import generate_mesh, generate_grid # 定义物理域x ∈ [0, 2000] m, z ∈ [0, 1000] m x_min, x_max 0.0, 2000.0 z_min, z_max 0.0, 1000.0 # 空间采样dx dz 10 m → 201 × 101 网格点 dx dz 10.0 nx int((x_max - x_min) / dx) 1 nz int((z_max - z_min) / dz) 1 # 时间参数最大模拟时间 1.0 s采样间隔 0.002 s → 501 个时间步 t_max 1.0 dt 0.002 nt int(t_max / dt) 1 # 生成笛卡尔网格 mesh generate_mesh(CartesianDomain(x_min, x_max, z_min, z_max), nx, nz)generate_mesh返回pysit.mesh.Mesh对象包含节点坐标、单元连接关系是后续所有物理量定义的载体。3.2 构建速度模型含三层结构# 初始化均匀背景速度 2000 m/s vel np.full((nz, nx), 2000.0, dtypenp.float64) # 添加第二层z ∈ [300, 600] m速度 2500 m/s vel[30:61, :] 2500.0 # 注意z 方向索引从上到下30→60 对应 300→600 m # 添加第三层z 600 m速度 3000 m/s vel[61:, :] 3000.0 # 将速度模型绑定到网格 model ConstantDensityAcousticWave(mesh, velocityvel)此处vel是 NumPy 数组形状(nz, nx)符合 MATLAB 中size(vel) [101, 201]的惯例。后续保存为.mat时维度顺序无需转换。3.3 设置震源与接收器# 震源位置x1000 m, z50 m近地表 source_x 1000.0 source_z 50.0 source PointSource(mesh, (source_x, source_z)) # 接收器线地表 z0x 从 200 到 1800 m每 20 m 一个共 81 个 receiver_x np.linspace(200.0, 1800.0, 81) receiver_z np.full_like(receiver_x, 0.0) receivers PointReceiver(mesh, list(zip(receiver_x, receiver_z)))PointSource和PointReceiver自动将物理坐标映射到最近网格点避免手动插值误差。3.4 执行正演并保存波场快照为.matfrom pysit import Modeling # 创建正演引擎 solver Modeling(model, solverdevito) # 或 fd 使用内置有限差分 # 生成震源时间函数Ricker 子波主频 10 Hz f0 10.0 wavelet solver.wavelet(ricker, f0f0, dtdt, ntnt) # 执行正演获取接收器记录 data_rec 和波场快照 snapshots data_rec, snapshots solver.forward(source, receivers, wavelet, return_snapshotsTrue) # 保存接收数据shape (81, 501)MATLAB 中为 [nt, nrec] scipy.io.savemat(forward_data.mat, { data: data_rec.T, # 转置使 MATLAB 中 size(data) [501, 81] dt: dt, nrec: len(receivers), nt: nt, f0: f0 }) # 保存第 100、200、300 个时间步的波场快照节省空间 snapshot_times [100, 200, 300] for t in snapshot_times: scipy.io.savemat(fsnapshot_t{t}.mat, { snapshot: snapshots[t], # shape (nz, nx) time_step: t, dt: dt, t_sec: t * dt })data_rec是(nrec, nt)数组MATLAB 中习惯data(nt, nrec)故显式转置。snapshots是列表每个元素为(nz, nx)二维数组直接保存即可。提示return_snapshotsTrue会显著增加内存占用。生产中建议用return_snapshotsFalse改用solver.forward(..., save_snapshots[100,200,300])参数控制保存时机避免 OOM。4. 执行全波形反演FWI从初始模型、目标函数到梯度更新的全流程反演目标是给定正演生成的data_rec作为观测数据从一个粗糙的初始速度模型出发迭代更新模型参数使模拟数据与观测数据的 L2 范数最小化。PySIT 的Inversion类封装了这一过程。4.1 构建初始模型与反演配置from pysit import Inversion, LBFGS # 初始模型均匀 2200 m/s故意偏离真实模型制造反演挑战 vel_init np.full((nz, nx), 2200.0, dtypenp.float64) model_init ConstantDensityAcousticWave(mesh, velocityvel_init) # 定义反演目标最小化 (simulated_data - observed_data)^2 objective Inversion(model_init, objectiveleast-squares, datadata_rec, # 观测数据 sourcesource, receiversreceivers, waveletwavelet, dtdt) # 选择优化器L-BFGS-B支持边界约束防止速度变为负数 optimizer LBFGS(maxiter20, bounds(1500, 4000)) # 速度物理范围 [1500, 4000] m/sbounds参数至关重要无约束反演极易产生非物理解如负速度、超光速导致后续正演崩溃。此处设定合理地质先验。4.2 运行反演并监控收敛# 执行反演返回最终模型和历史记录 result optimizer(objective) # result.model.velocity 是更新后的速度数组 final_vel result.model.velocity # 保存反演过程中的目标函数值、梯度范数 scipy.io.savemat(fwi_history.mat, { obj_vals: np.array(result.objective_values), # 每次迭代的目标函数值 grad_norms: np.array(result.gradient_norms), # 梯度 L2 范数 iterations: len(result.objective_values), final_vel: final_vel }) # 保存最终模型供 MATLAB 地质解释软件读取 scipy.io.savemat(fwi_result.mat, { velocity: final_vel, x_min: x_min, x_max: x_max, z_min: z_min, z_max: z_max, dx: dx, dz: dz, dt: dt, f0: f0, niter: len(result.objective_values) })result.objective_values是列表记录每次迭代的残差平方和result.gradient_norms是梯度下降步长的量化指标。二者共同构成反演质量评估依据。4.3 关键参数调优表影响收敛速度与稳定性的 3 个必调项参数作用推荐值调整逻辑maxiterLBFGS最大迭代次数15–30初次运行设为 20若obj_vals末尾仍快速下降可增至 30若 10 次后已平稳可减至 15 节省时间bounds速度上下界防止非物理解[1500, 4000]陆上[1300, 2500]浅海必须基于区域地质知识设定过窄限制更新空间过宽引入虚假极小值wavelet.f0子波主频控制反演分辨率5–15 Hz低频5–8 Hz主导宏观结构高频10–15 Hz恢复细节建议先用 8 Hz 反演粗模型再用 12 Hz 精修注意wavelet.f0在反演中必须与正演一致。若正演用 10 Hz Ricker反演中wavelet必须用相同f0否则目标函数不可导梯度计算失效。5. 验证反演结果用.mat文件在 MATLAB 中可视化对比识别典型失真模式反演完成不等于成功。必须将fwi_result.mat加载到 MATLAB 中与真实模型vel已存为true_model.mat做像素级对比并分析误差分布。这是地质解释人员验收的硬性环节。5.1 MATLAB 端加载与基础绘图% 加载真实模型与反演结果 load(true_model.mat); % 假设保存时 key 为 vel load(fwi_result.mat); % 假设保存时 key 为 velocity % 绘制对比图双栏布局 figure(Position, [100, 100, 1200, 500]); subplot(1,2,1); imagesc(x_min:dx:x_max, z_min:dz:z_max, vel); axis xy; colorbar; title(True Model (m/s)); xlabel(X (m)); ylabel(Z (m)); subplot(1,2,2); imagesc(x_min:dx:x_max, z_min:dz:z_max, velocity); axis xy; colorbar; title(FWI Result (m/s)); xlabel(X (m)); ylabel(Z (m));注意Python 中vel形状为(nz, nx)MATLAB 中imagesc默认(y,x)故用vel转置保证坐标轴方向一致。5.2 识别三类典型失真并定位原因反演结果常出现以下模式需结合fwi_history.mat中的obj_vals曲线诊断失真模式MATLAB 可视化特征对应obj_vals曲线形态根本原因与修复低频模糊整体趋势正确但层界面模糊、厚度不准下降迅速后早衰10 次内停滞初始模型太差或低频信息不足 → 降低wavelet.f0至 5 Hz或添加低频震源高频振铃层内出现密集条纹状伪影下降缓慢后期震荡正则化不足 → 在Inversion初始化时添加regularizationTikhonov参数边界畸变模型四角/边缘速度异常升高或降低残差持续不降grad_norms波动大边界条件反射干扰 → 在Modeling中启用absorbing_boundaryTrue例如启用 Tikhonov 正则化objective Inversion(model_init, objectiveleast-squares, datadata_rec, sourcesource, receiversreceivers, waveletwavelet, dtdt, regularizationTikhonov, # 添加此行 reg_param1e-4) # 正则化权重需试错调整reg_param过大会过度平滑过小则无效。建议从1e-5开始按×10步进测试。5.3 用 Python 快速生成误差统计报告无需切换 MATLAB用 Python 直接计算定量指标import scipy.io import numpy as np # 加载模型 true scipy.io.loadmat(true_model.mat)[vel] fwi scipy.io.loadmat(fwi_result.mat)[velocity] # 计算 RMSE、MAPE、结构相似性 SSIM rmse np.sqrt(np.mean((fwi - true)**2)) mape np.mean(np.abs((fwi - true) / true)) * 100 # SSIM 需安装 scikit-image: pip install scikit-image from skimage.metrics import structural_similarity ssim, _ structural_similarity(true, fwi, fullTrue, data_rangetrue.max()-true.min()) print(fRMSE: {rmse:.2f} m/s | MAPE: {mape:.2f}% | SSIM: {ssim:.3f}) # 输出示例RMSE: 128.45 m/s | MAPE: 5.21% | SSIM: 0.873RMSE 150 m/s、MAPE 6%、SSIM 0.85 是陆上二维 FWI 的常见验收阈值。低于此值需回溯检查初始模型、震源频率或正则化参数。提示.mat文件中的velocity字段在 Python 和 MATLAB 中均为 double 类型无精度损失。地质解释软件如 Petrel、GeoFrame读取时自动识别为网格属性无需额外转换。本文还有配套的精品资源点击获取 PySIT地震波正反演与.mat文件互通实践
简介本资源是一套面向地球物理建模与反演研究者的Python实践教程聚焦pysit地震成像工具的正演模拟、反演优化及MATLAB数据互通适用于具备基础Python和数值计算能力的科研人员与高年级本科生。压缩包共305个文件以179个Python脚本含模型构建、源接收器配置、正反演核心逻辑为主体辅以58个reStructuredText文档提供API说明与流程注解、14张可视化结果图如波场快照、模型更新对比及少量C/C头文件支撑底层求解器调用整体体积仅1.8MB轻量易部署。已有294人学习下载资源结构清晰包含完整可运行示例如pysit-example-master目录、详细配置文件pysit.cfg、make.bat及部署手册.docx覆盖从环境安装、参数设置、迭代反演到.mat格式结果导出的全流程特别适合开展地震成像算法验证、跨平台Python-MATLAB联合分析或课程实验复现。1. 用 PySIT 做地震波正演模拟、反演成像并把模型/数据存成.mat文件这不是 MATLAB 专属流程你手头有一套地下介质参数比如速度场想看看地震波在其中怎么传播——这是正演反过来你有一组实际采集的地震记录想反推地下结构长什么样——这是反演。传统上这类计算常被默认绑定在 MATLAB 生态里尤其当结果要交给地质解释人员时.mat文件几乎是交付标配。但现实是团队里 Python 工程师越来越多MATLAB 许可成本高、部署难、CI/CD 集成卡顿。PySIT 正是为这种场景设计的——它用纯 Python 实现了有限差分波场模拟、全波形反演FWI、最小二乘反演等核心算法且原生支持 NumPy天然兼容 SciPy、Matplotlib更重要的是它不依赖 MATLAB 运行时却能无缝对接.mat文件的读写。本文面向已掌握 Python 基础NumPy/Pandas、了解波动方程基本概念的地球物理建模者或算法工程师不讲泛泛而谈的“Python 多好”只聚焦如何用 PySIT 跑通一个完整正反演闭环并确保中间模型、观测数据、反演结果全部按地质行业通行标准存为.mat文件——不是用scipy.io.savemat粗暴打包而是结构清晰、字段可读、MATLAB 端开箱即用。2. 安装 PySIT 及配套生态避开 conda-forge 的旧版本陷阱用源码编译保障反演稳定性PySIT 并未发布到 PyPI 主索引官方推荐安装方式是通过 GitHub 源码构建。直接pip install pysit会失败而conda install -c conda-forge pysit虽然能装上但截至 2024 年中conda-forge 通道中的最新版仍是 0.4.0发布于 2021 年缺失对scipy1.10的适配且反演模块在多线程下存在内存泄漏风险。生产环境必须使用当前主干main branch代码。2.1 依赖预检与环境隔离PySIT 重度依赖scipy用于稀疏矩阵求解器、numba加速波场传播内核、h5py可选用于大体积数据存储。注意numba必须与llvmlite版本严格匹配否则编译失败。建议新建独立虚拟环境python -m venv pysit_env source pysit_env/bin/activate # Linux/macOS # pysit_env\Scripts\activate.bat # Windows提示不要用--system-site-packages。PySIT 对scipy.sparse.linalg的调用路径敏感混用系统级 scipy 极易触发AttributeError: module object has no attribute cg类错误。2.2 源码编译安装含关键补丁执行以下命令拉取最新代码并安装git clone https://github.com/pysit/pysit.git cd pysit # 应用社区维护的关键补丁修复 FWI 中梯度归一化失效问题 curl -sL https://patch-diff.githubusercontent.com/raw/pysit/pysit/pull/187.patch | git apply pip install -e .[dev]该命令中-e表示开发模式安装后续修改源码可即时生效[dev]安装额外测试依赖如pytest。验证安装是否成功import pysit print(pysit.__version__) # 输出应为类似 0.5.0.dev0 from pysit import DevitoDomain print(PySIT 安装就绪)若报错ModuleNotFoundError: No module named devito说明未启用 Devito 后端PySIT 支持多种求解器后端。此时需单独安装 Devito仅当需要高阶精度或复杂边界条件时才必需pip install devito4.7.3 # 固定版本避免与 PySIT 0.5.x 不兼容2.3.mat文件支持scipy.io是唯一可靠选择PySIT 本身不处理文件 I/O.mat读写完全交由scipy.io。注意scipy1.9.0开始默认使用 MATLAB v7.3 格式HDF5 封装该格式支持大于 2GB 的数组且 MATLAB R2016b 全面兼容。但旧版scipy1.8.0默认生成 v7 格式无法保存超过 2GB 的三维速度模型。务必检查版本python -c import scipy; print(scipy.__version__) # 必须 ≥ 1.9.0若版本过低升级pip install --upgrade scipy注意scipy.io.savemat生成的.mat文件在 MATLAB 中加载后变量名即为字典 key。例如savemat(model.mat, {vel: vel_array, mesh: mesh_dict})MATLAB 中load(model.mat)后直接可用vel和mesh变量——这是地质解释流程中最友好的交互方式。3. 构建正演模型从网格定义、源接收器布设到波场快照保存正演是反演的基础。PySIT 使用pysit.modeling.Modeling类封装波场模拟逻辑。本节以二维声波方程为例构建一个含 3 层介质的合成模型并生成 100 个时间步的波场快照。3.1 定义计算域与离散化参数import numpy as np from pysit import ConstantDensityAcousticWave, CartesianDomain, PointSource, PointReceiver from pysit import generate_mesh, generate_grid # 定义物理域x ∈ [0, 2000] m, z ∈ [0, 1000] m x_min, x_max 0.0, 2000.0 z_min, z_max 0.0, 1000.0 # 空间采样dx dz 10 m → 201 × 101 网格点 dx dz 10.0 nx int((x_max - x_min) / dx) 1 nz int((z_max - z_min) / dz) 1 # 时间参数最大模拟时间 1.0 s采样间隔 0.002 s → 501 个时间步 t_max 1.0 dt 0.002 nt int(t_max / dt) 1 # 生成笛卡尔网格 mesh generate_mesh(CartesianDomain(x_min, x_max, z_min, z_max), nx, nz)generate_mesh返回pysit.mesh.Mesh对象包含节点坐标、单元连接关系是后续所有物理量定义的载体。3.2 构建速度模型含三层结构# 初始化均匀背景速度 2000 m/s vel np.full((nz, nx), 2000.0, dtypenp.float64) # 添加第二层z ∈ [300, 600] m速度 2500 m/s vel[30:61, :] 2500.0 # 注意z 方向索引从上到下30→60 对应 300→600 m # 添加第三层z 600 m速度 3000 m/s vel[61:, :] 3000.0 # 将速度模型绑定到网格 model ConstantDensityAcousticWave(mesh, velocityvel)此处vel是 NumPy 数组形状(nz, nx)符合 MATLAB 中size(vel) [101, 201]的惯例。后续保存为.mat时维度顺序无需转换。3.3 设置震源与接收器# 震源位置x1000 m, z50 m近地表 source_x 1000.0 source_z 50.0 source PointSource(mesh, (source_x, source_z)) # 接收器线地表 z0x 从 200 到 1800 m每 20 m 一个共 81 个 receiver_x np.linspace(200.0, 1800.0, 81) receiver_z np.full_like(receiver_x, 0.0) receivers PointReceiver(mesh, list(zip(receiver_x, receiver_z)))PointSource和PointReceiver自动将物理坐标映射到最近网格点避免手动插值误差。3.4 执行正演并保存波场快照为.matfrom pysit import Modeling # 创建正演引擎 solver Modeling(model, solverdevito) # 或 fd 使用内置有限差分 # 生成震源时间函数Ricker 子波主频 10 Hz f0 10.0 wavelet solver.wavelet(ricker, f0f0, dtdt, ntnt) # 执行正演获取接收器记录 data_rec 和波场快照 snapshots data_rec, snapshots solver.forward(source, receivers, wavelet, return_snapshotsTrue) # 保存接收数据shape (81, 501)MATLAB 中为 [nt, nrec] scipy.io.savemat(forward_data.mat, { data: data_rec.T, # 转置使 MATLAB 中 size(data) [501, 81] dt: dt, nrec: len(receivers), nt: nt, f0: f0 }) # 保存第 100、200、300 个时间步的波场快照节省空间 snapshot_times [100, 200, 300] for t in snapshot_times: scipy.io.savemat(fsnapshot_t{t}.mat, { snapshot: snapshots[t], # shape (nz, nx) time_step: t, dt: dt, t_sec: t * dt })data_rec是(nrec, nt)数组MATLAB 中习惯data(nt, nrec)故显式转置。snapshots是列表每个元素为(nz, nx)二维数组直接保存即可。提示return_snapshotsTrue会显著增加内存占用。生产中建议用return_snapshotsFalse改用solver.forward(..., save_snapshots[100,200,300])参数控制保存时机避免 OOM。4. 执行全波形反演FWI从初始模型、目标函数到梯度更新的全流程反演目标是给定正演生成的data_rec作为观测数据从一个粗糙的初始速度模型出发迭代更新模型参数使模拟数据与观测数据的 L2 范数最小化。PySIT 的Inversion类封装了这一过程。4.1 构建初始模型与反演配置from pysit import Inversion, LBFGS # 初始模型均匀 2200 m/s故意偏离真实模型制造反演挑战 vel_init np.full((nz, nx), 2200.0, dtypenp.float64) model_init ConstantDensityAcousticWave(mesh, velocityvel_init) # 定义反演目标最小化 (simulated_data - observed_data)^2 objective Inversion(model_init, objectiveleast-squares, datadata_rec, # 观测数据 sourcesource, receiversreceivers, waveletwavelet, dtdt) # 选择优化器L-BFGS-B支持边界约束防止速度变为负数 optimizer LBFGS(maxiter20, bounds(1500, 4000)) # 速度物理范围 [1500, 4000] m/sbounds参数至关重要无约束反演极易产生非物理解如负速度、超光速导致后续正演崩溃。此处设定合理地质先验。4.2 运行反演并监控收敛# 执行反演返回最终模型和历史记录 result optimizer(objective) # result.model.velocity 是更新后的速度数组 final_vel result.model.velocity # 保存反演过程中的目标函数值、梯度范数 scipy.io.savemat(fwi_history.mat, { obj_vals: np.array(result.objective_values), # 每次迭代的目标函数值 grad_norms: np.array(result.gradient_norms), # 梯度 L2 范数 iterations: len(result.objective_values), final_vel: final_vel }) # 保存最终模型供 MATLAB 地质解释软件读取 scipy.io.savemat(fwi_result.mat, { velocity: final_vel, x_min: x_min, x_max: x_max, z_min: z_min, z_max: z_max, dx: dx, dz: dz, dt: dt, f0: f0, niter: len(result.objective_values) })result.objective_values是列表记录每次迭代的残差平方和result.gradient_norms是梯度下降步长的量化指标。二者共同构成反演质量评估依据。4.3 关键参数调优表影响收敛速度与稳定性的 3 个必调项参数作用推荐值调整逻辑maxiterLBFGS最大迭代次数15–30初次运行设为 20若obj_vals末尾仍快速下降可增至 30若 10 次后已平稳可减至 15 节省时间bounds速度上下界防止非物理解[1500, 4000]陆上[1300, 2500]浅海必须基于区域地质知识设定过窄限制更新空间过宽引入虚假极小值wavelet.f0子波主频控制反演分辨率5–15 Hz低频5–8 Hz主导宏观结构高频10–15 Hz恢复细节建议先用 8 Hz 反演粗模型再用 12 Hz 精修注意wavelet.f0在反演中必须与正演一致。若正演用 10 Hz Ricker反演中wavelet必须用相同f0否则目标函数不可导梯度计算失效。5. 验证反演结果用.mat文件在 MATLAB 中可视化对比识别典型失真模式反演完成不等于成功。必须将fwi_result.mat加载到 MATLAB 中与真实模型vel已存为true_model.mat做像素级对比并分析误差分布。这是地质解释人员验收的硬性环节。5.1 MATLAB 端加载与基础绘图% 加载真实模型与反演结果 load(true_model.mat); % 假设保存时 key 为 vel load(fwi_result.mat); % 假设保存时 key 为 velocity % 绘制对比图双栏布局 figure(Position, [100, 100, 1200, 500]); subplot(1,2,1); imagesc(x_min:dx:x_max, z_min:dz:z_max, vel); axis xy; colorbar; title(True Model (m/s)); xlabel(X (m)); ylabel(Z (m)); subplot(1,2,2); imagesc(x_min:dx:x_max, z_min:dz:z_max, velocity); axis xy; colorbar; title(FWI Result (m/s)); xlabel(X (m)); ylabel(Z (m));注意Python 中vel形状为(nz, nx)MATLAB 中imagesc默认(y,x)故用vel转置保证坐标轴方向一致。5.2 识别三类典型失真并定位原因反演结果常出现以下模式需结合fwi_history.mat中的obj_vals曲线诊断失真模式MATLAB 可视化特征对应obj_vals曲线形态根本原因与修复低频模糊整体趋势正确但层界面模糊、厚度不准下降迅速后早衰10 次内停滞初始模型太差或低频信息不足 → 降低wavelet.f0至 5 Hz或添加低频震源高频振铃层内出现密集条纹状伪影下降缓慢后期震荡正则化不足 → 在Inversion初始化时添加regularizationTikhonov参数边界畸变模型四角/边缘速度异常升高或降低残差持续不降grad_norms波动大边界条件反射干扰 → 在Modeling中启用absorbing_boundaryTrue例如启用 Tikhonov 正则化objective Inversion(model_init, objectiveleast-squares, datadata_rec, sourcesource, receiversreceivers, waveletwavelet, dtdt, regularizationTikhonov, # 添加此行 reg_param1e-4) # 正则化权重需试错调整reg_param过大会过度平滑过小则无效。建议从1e-5开始按×10步进测试。5.3 用 Python 快速生成误差统计报告无需切换 MATLAB用 Python 直接计算定量指标import scipy.io import numpy as np # 加载模型 true scipy.io.loadmat(true_model.mat)[vel] fwi scipy.io.loadmat(fwi_result.mat)[velocity] # 计算 RMSE、MAPE、结构相似性 SSIM rmse np.sqrt(np.mean((fwi - true)**2)) mape np.mean(np.abs((fwi - true) / true)) * 100 # SSIM 需安装 scikit-image: pip install scikit-image from skimage.metrics import structural_similarity ssim, _ structural_similarity(true, fwi, fullTrue, data_rangetrue.max()-true.min()) print(fRMSE: {rmse:.2f} m/s | MAPE: {mape:.2f}% | SSIM: {ssim:.3f}) # 输出示例RMSE: 128.45 m/s | MAPE: 5.21% | SSIM: 0.873RMSE 150 m/s、MAPE 6%、SSIM 0.85 是陆上二维 FWI 的常见验收阈值。低于此值需回溯检查初始模型、震源频率或正则化参数。提示.mat文件中的velocity字段在 Python 和 MATLAB 中均为 double 类型无精度损失。地质解释软件如 Petrel、GeoFrame读取时自动识别为网格属性无需额外转换。本文还有配套的精品资源点击获取