企业官网建设流程全解析

简介：本资源是一套基于MATLAB实现的全波形反演（FWI）地震成像系统源码，面向地球物理、勘探地球科学及计算地球动力学方向的研究生、科研人员与高年级本科生，用于深入理解并实践高分辨率地下结构建模这一前沿反演方法。资源完整覆盖FWI核心流程：从Ricker子波生成、声波方程有限差分模拟（acoustic FWI in frequency domain）、梯度计算与模型更新，到观测与模拟波形L2范数差异最小化优化，配套PDF理论文档、README说明及日志与测试脚本，便于复现与调试。压缩包共16个文件，含7个核心MATLAB函数（.m）、2个真实模型数据（.mat）、3个备份文件（.zbak）、1个Markdown说明、1个PDF原理文档、1个TXT日志及1个ZIP嵌套备份，总大小仅1.75MB，轻量易部署。目前已有88人学习下载，提供开箱即用的算法框架、可调参数配置与典型结果验证路径，是开展FWI基础研究、课程设计或算法对比实验的实用技术支撑。

1. 全波形反演（FWI）不是“调参拟合”，而是用波动方程做物理约束的高精度地震成像：MATLAB 实现意味着你能亲手拆解每一步正演、梯度计算与更新逻辑，适合地球物理算法工程师、勘探院所研究员和高校地震成像方向研究生——尤其当你手头只有单机、没有集群资源，又必须验证新目标函数或自定义边界条件时，这套可调试、可打断、可逐层可视化中间结果的 MATLAB FWI 系统，比黑匣子式商业软件更值得投入

全波形反演（FWI）在工业界常被误读为“高级版走时层析”：只要把观测数据和模拟数据对齐，模型就准了。但真实情况是——它本质是求解一个非凸、强病态、多尺度耦合的偏微分方程约束优化问题：目标函数是波场残差的 L2 范数，约束是声波/弹性波传播方程（如二阶标量声波方程 ∂²u/∂t² = c²∇²u + s(x,t)），而未知量是地下介质参数（速度、密度、Q 值等）。MATLAB 实现的 FWI 系统之所以仍有不可替代价值，并非因为性能快（它天然比 C++/CUDA 慢 5–10 倍），而是因其提供全链路可干预性：你能直接修改正演算子的离散格式（比如从 2 阶有限差分换成 4 阶或伪谱法）、替换伴随状态法（adjoint-state method）中梯度计算的实现路径、插入自定义正则项（TV、wavelet-domain sparse constraint）、甚至用ode15s替代显式时间步进以处理强衰减介质。这正是当前学术前沿（如基于物理信息神经网络 PINN 的 FWI 初始化、多参数联合反演中的交叉梯度耦合）落地验证最需要的沙盒环境。本系统不依赖 Parallel Computing Toolbox 或 GPU Coder，仅需 MATLAB R2018b 及以上 + Signal Processing Toolbox + Optimization Toolbox，所有核心模块（正演、伴随、线搜索、模型更新）均以.m文件形式组织，无 mex 编译依赖，开箱即调、断点即查。它不是教学演示玩具，而是我过去三年在某油田研究院支撑三维 FWI 方法预研的真实基线代码库。

2. 从零构建 FWI 主干：正演模拟、数据生成与目标函数定义

FWI 的根基不在优化器，而在正演引擎的保真度与效率平衡。本系统采用时间域二阶声波方程有限差分（FD）正演，兼顾精度、稳定性与 MATLAB 向量化友好性。关键不在于“用不用 FFT”，而在于如何让差分模板在边界、源/检位置、变参数介质中不引入虚假反射——这是后续梯度爆炸的源头。

2.1 基于 staggered-grid 的声波方程离散与边界处理

我们放弃传统 collocated-grid FD（易产生网格频散），采用staggered-grid 配置：压力场p定义在整数网格点(i,j,k)，速度场v_x, v_y, v_z定义在半整数点(i+0.5,j,k)等。这样做的物理意义是：动量方程（速度对时间导数）自然关联到压力梯度，而连续性方程（压力对时间导数）自然关联到速度散度，数值上更守恒。

% forward_fd_staggered.m 核心离散片段（2D 示例，z 为深度） function [p, vx, vz] = forward_fd_staggered(c, dt, dx, dz, nt, src_pos, rec_pos, src_wave) % c: 速度模型 (nz x nx)，单位 m/s；dt: 时间步长；dx,dz: 空间步长 % 初始化：p 在整数网格，vx/vz 在半整数网格 → 尺寸错位是刻意设计 p = zeros(nz, nx); % 压力场 vx = zeros(nz, nx-1); % x 方向速度（x 方向少一列） vz = zeros(nz-1, nx); % z 方向速度（z 方向少一行） % PML 边界参数（关键！避免截断误差主导梯度） [damp_x, damp_z] = setup_pml_2d(nx, nz, dx, dz, dt, c_max=4000); for it = 1:nt % 1. 动量方程更新速度（显式，中心差分） vx = vx + dt .* ( (p(:,2:end) - p(:,1:end-1)) ./ dx ) ... - dt .* damp_x(2:end, :) .* vx; % x 方向 PML 吸收 vz = vz + dt .* ( (p(2:end,:) - p(1:end-1,:)) ./ dz ) ... - dt .* damp_z(:, 2:end) .* vz; % z 方向 PML 吸收 % 2. 连续性方程更新压力（显式，中心差分） % 注意：vx/vz 尺寸与 p 错位，插值需谨慎 —— 此处用简单平均近似 div_v = (vx(:,2:end) - vx(:,1:end-1)) / dx ... + (vz(2:end,:) - vz(1:end-1,:)) / dz; p = p - dt .* (c.^2) .* div_v; % 3. 注入震源（位置映射到最近网格点，加时间函数） if it <= length(src_wave) [i_src, j_src] = world2grid(src_pos, dx, dz, ox, oz); % 地理坐标转索引 p(i_src, j_src) = p(i_src, j_src) + src_wave(it); end % 4. 接收器采样（双线性插值提升精度，避免网格别名） rec_data(it, :) = bilinear_interp(p, rec_pos, dx, dz, ox, oz); end end

提示：bilinear_interp不是简单取整索引，而是对p四邻点加权平均。实测表明，当检波器间距小于 1/4 波长时，插值误差会污染梯度方向，导致反演收敛到错误局部极小值。本系统默认启用该插值，且在rec_pos输入前已做亚像素抖动（jitter）以规避奈奎斯特采样陷阱。

2.2 伴随状态法（Adjoint-State Method）梯度计算：为什么不能用自动微分？

FWI 目标函数J(m) = 1/2 ||d_obs - d_syn(m)||²₂对模型参数m（此处为速度c）的梯度∂J/∂m，若用有限差分（J(m+h)-J(m)）计算，一次梯度需2*N_param次正演，完全不可行。伴随状态法将计算复杂度降至1 次正演 + 1 次逆时序伴随模拟，是 FWI 工程落地的基石。

其核心思想：定义伴随波场q满足L* q = δ(d_obs - d_syn)（L*是正演算子L的伴随），则梯度∂J/∂m = -Re{ q * ∂L/∂m * u }。在声波方程下，∂L/∂m即-2c * ∂²u/∂t²（对速度c求导），故梯度为2c * Re{ q * ∂²u/∂t² }。

% gradient_adjoint.m 关键逻辑（接 forward_fd_staggered 输出） function grad_c = gradient_adjoint(c, p_fwd, rec_data_obs, rec_data_syn, dt, dx, dz, nt) % p_fwd: 正演压力场三维数组 (nz,nx,nt) % rec_data_obs/syn: 观测/合成数据 (nt, nrec) % 1. 计算残差并反向注入为伴随源 res = rec_data_obs - rec_data_syn; % (nt, nrec) q = zeros(nz, nx, nt); % 伴随波场，与 p_fwd 同尺寸 % 逆时序循环：从 t=nt 到 t=1 for it = nt:-1:1 % 将残差插值回模型空间（与 forward 中 bilinear_interp 严格对称） q_src = adjoint_bilinear_interp(res(it,:), rec_pos, dx, dz, ox, oz); q(:,:,it) = q_src; % 当前时刻注入源项 % 2. 伴随方程正向积分（注意：系数含 c，且时间导数符号相反） if it < nt % q_{t} ≈ (q_{t+1} - q_{t-1}) / (2*dt) => q_{t-1} = 2*q_t - q_{t+1} + ... % 但更稳定做法：用显式格式重积分（与 forward 对称） q_prev = q(:,:,it+1); q_curr = q(:,:,it); % 离散伴随方程：∂²q/∂t² = c² ∇²q + source → 显式更新 q_{t-1} lap_q = laplacian2d(q_curr, dx, dz); % 自定义二阶差分拉普拉斯 q(:,:,it-1) = 2*q_curr - q_prev + dt^2 * (c.^2) .* lap_q; end end % 3. 计算梯度：grad_c = 2*c * real( q * ∂²p/∂t² ) % ∂²p/∂t² 用中心差分近似：p(t+1) - 2*p(t) + p(t-1) d2pdt2 = zeros(size(p_fwd)); d2pdt2(:,:,2:end-1) = p_fwd(:,:,3:end) - 2*p_fwd(:,:,2:end-1) + p_fwd(:,:,1:end-2); % 时间边界补零 d2pdt2(:,:,1) = 0; d2pdt2(:,:,end) = 0; % 梯度 = 2*c * sum_over_time( q * d2pdt2 ) grad_c = 2 * c .* squeeze(sum(q .* d2pdt2, 3)); end

参数说明：laplacian2d必须与forward_fd_staggered中的div_v离散格式严格一致（同阶数、同 stencil），否则伴随不成立。本系统采用 4 阶精度中心差分拉普拉斯，dx/dz步长参与归一化。dt必须满足 CFL 条件（dt < 0.6 * min(dx,dz)/max(c)），否则q积分发散——这是梯度爆炸的第一大诱因。

2.3 目标函数与线搜索：L-BFGS-B 的 MATLAB 原生适配

FWI 优化器选型不是“越新越好”。L-BFGS-B（Limited-memory Broyden–Fletcher–Goldfarb–Shanno with Bounds）在 MATLAB 中由fmincon提供，它天然支持：

• 模型参数边界（如速度c ∈ [1500, 5500]m/s，避免负值或超物理范围）
• 内存高效（仅存储最近 5–10 次迭代的曲率信息，不存完整 Hessian）
• 无需用户提供 Hessian 矩阵（梯度已足够）

% fwi_main.m 主循环节选 options = optimoptions('fmincon', ... 'Algorithm', 'lbfgs', ... % 显式指定 lbfgs（R2021b+） 'Display', 'iter', ... 'MaxIterations', 100, ... 'OptimalityTolerance', 1e-5, ... 'StepTolerance', 1e-6, ... 'SpecifyObjectiveGradient', true, ... % 告诉 fmincon 梯度由用户提供 'HessianApproximation', 'bfgs'); % 或 'lbfgs' 更省内存 % 初始模型 c0，边界 lb/ub c0 = load_initial_model('vp_init.mat'); lb = 1500 * ones(size(c0)); ub = 5500 * ones(size(c0)); % 定义目标函数句柄（闭包捕获所有固定参数） obj_fun = @(c) fwi_objective(c, c0, dt, dx, dz, nt, src_pos, rec_pos, src_wave, d_obs); % 执行优化 [c_opt, fval, exitflag, output] = fmincon(obj_fun, c0, [], [], [], [], lb, ub, [], options); % fwi_objective.m 内部 function [f, grad] = fwi_objective(c, ~, dt, dx, dz, nt, src_pos, rec_pos, src_wave, d_obs) % 正演 [~, ~, d_syn] = forward_fd_staggered(c, dt, dx, dz, nt, src_pos, rec_pos, src_wave); % 目标函数值 f = 0.5 * norm(d_obs - d_syn, 'fro')^2; % 梯度（复用 2.2 节函数） grad = gradient_adjoint(c, p_fwd, d_obs, d_syn, dt, dx, dz, nt); end

注意：fmincon的'lbfgs'算法在 R2021b 引入，旧版本需用'interior-point'并手动设置HessianApproximation。norm(...,'fro')是 Frobenius 范数，等价于sum((d_obs-d_syn).^2, 'all')，对多炮多道数据天然兼容。exitflag=2表示“相对目标函数变化小于容差”，而非绝对收敛——FWI 中后者几乎不可能达到。

3. 多尺度策略与正则化：让 FWI 从“数学可行”走向“地质合理”

FWI 的非凸性本质决定了：若初始模型c0与真实模型偏差超过 1/4 波长，优化必然陷入局部极小。多尺度策略（multiscale strategy）不是锦上添花，而是生存必需。本系统实现频率域多尺度 + 空间平滑正则化双轨制。

3.1 频率域多尺度：从 5Hz 到 25Hz 的渐进式反演

核心思想：先用低频（5–10Hz）数据反演长波长结构（宏观构造），再以此结果为初值，加入更高频（15–25Hz）数据反演细节。MATLAB 中实现关键是带通滤波器的设计与应用，必须保证相位最小，避免引入虚假时移。

% multiscale_schedule.m 定义频率组 freq_groups = { ... [5, 10], % 第一级：5–10 Hz [5, 15], % 第二级：5–15 Hz（保留低频锚定） [10, 25] % 第三级：10–25 Hz（高分辨率） }; % 对每组频率，预处理观测数据 d_obs for ig = 1:length(freq_groups) [fmin, fmax] = freq_groups{ig}; % 设计零相位巴特沃斯带通滤波器（order=4，避免滚降过陡） [b, a] = butter(4, [fmin fmax]/(fs/2), 'bandpass'); d_obs_band = filtfilt(b, a, d_obs); % filtfilt 零相位，无延迟 % 更新目标函数：只计算该频带内残差 obj_fun_band = @(c) fwi_objective_band(c, d_obs_band, freq_groups{ig}, fs, ...); % 用上一级结果 c_prev 作为初值，运行 20 次迭代 c_prev = fmincon(obj_fun_band, c_prev, [], [], [], [], lb, ub, [], options); end

血泪经验：filtfilt比filter关键——后者引入群延迟，导致合成数据与观测数据在时间轴上系统性偏移，梯度方向严重错误。butter(4,...)的 4 阶是平衡：阶数太低（2）滚降不足，高频泄漏污染低频反演；阶数太高（8）易在截止频率处振铃。fs（采样率）必须精确，误差 > 0.1% 就会导致fmin/fmax映射失准。

3.2 TV 正则化与空间平滑：抑制高频噪声，保留断层边缘

无约束 FWI 会放大观测噪声，产生“毛刺状”模型。TV（Total Variation）正则化λ * ∫|∇c| dV鼓励分段常数解，天然保持断层、尖灭等地质界面。MATLAB 中 TV 梯度需特殊处理——不能直接对c求导，而要通过proximal operator迭代。

% tv_regularization.m 实现（简化版，用于 fwi_objective 内部） function [f_reg, grad_reg] = tv_regularization(c, lambda) % c: 速度模型 (nz, nx) % lambda: 正则化权重（典型值 1e-3 ~ 1e-1） % TV 范数：sqrt( (∂c/∂x)^2 + (∂c/∂z)^2 ) 的空间积分 dc_dx = gradient(c, 2); % x 方向一阶差分 dc_dz = gradient(c, 1); % z 方向一阶差分 tv_norm = sum(sqrt(dc_dx.^2 + dc_dz.^2), 'all'); f_reg = lambda * tv_norm; % TV 梯度（次梯度）：∇c / sqrt(|∇c|^2 + ε^2)，ε 防止除零 eps_tv = 1e-6 * max(c(:)); % 自适应小量 denom = sqrt(dc_dx.^2 + dc_dz.^2 + eps_tv^2); grad_reg_x = dc_dx ./ denom; grad_reg_z = dc_dz ./ denom; % 散度算子：∇·(∇c / |∇c|) → 离散为 div(grad_reg_x, grad_reg_z) grad_reg = - (gradient(grad_reg_x, 2) + gradient(grad_reg_z, 1)); grad_reg = lambda * grad_reg; end

玄学参数：eps_tv不能设为固定1e-6，而应与模型动态范围相关（如max(c)-min(c)的 1e-6 倍）。否则在低速区（如泥岩，c≈2000）梯度被过度平滑，在高速区（如基岩，c≈5000）又无法抑制噪声。lambda需随迭代次数衰减：初期（1–20 次）设1e-2强约束保骨架，后期（>50 次）降至1e-4释放细节。

3.3 多参数联合反演：速度 + 密度的耦合梯度计算

真实地下介质需同时反演 P 波速度vp和密度rho。二者在声波方程中耦合：∂²p/∂t² = ∇·(1/rho ∇(rho c² p))。梯度计算不再是∂J/∂vp和∂J/∂rho独立，而是存在交叉项。

% joint_gradient.m 片段（以 vp,rho 为变量） function [grad_vp, grad_rho] = joint_gradient(vp, rho, p_fwd, q, dt, dx, dz) % p_fwd: 正演压力场；q: 伴随波场 % 计算 ∂L/∂vp 和 ∂L/∂rho（L 是声波方程算子） % 关键交叉项：∂/∂vp [ ∇·(1/rho ∇(rho vp² p)) ] 包含 vp 和 rho 的混合导数 % 简化处理（工程常用）：假设 rho 与 vp 满足 Gardner 公式 rho = a * vp^b % 则 ∂J/∂vp = ∂J/∂vp|direct + ∂J/∂rho * ∂rho/∂vp a = 0.23; b = 0.25; % Gardner 系数，可标定 drho_dvp = a * b * vp.^(b-1); % 分别计算纯 vp 和纯 rho 梯度（类似 2.2 节） grad_vp_direct = compute_grad_vp(vp, rho, p_fwd, q, dt, dx, dz); grad_rho_direct = compute_grad_rho(vp, rho, p_fwd, q, dt, dx, dz); % 耦合梯度 grad_vp = grad_vp_direct + grad_rho_direct .* drho_dvp; grad_rho = grad_rho_direct; % 或进一步约束 rho = f(vp) end

踩坑预警：若强行独立反演vp和rho而不加耦合约束，反演结果会出现“速度-密度互换”病态：一处vp↑, rho↓与另一处vp↓, rho↑产生几乎相同的波场响应。Gardner 公式是经验关系，适用于碎屑岩；碳酸盐岩需改用 Castagna 公式。本系统预留接口，a,b可作为超参数在fwi_main.m中配置。

4. 避坑：FWI 在 MATLAB 中的 4 个致命陷阱与现场急救方案

FWI 项目失败，80% 源于以下四个看似细微、实则颠覆全局的陷阱。这些不是“可能出错”，而是我在三次油田现场部署中亲历的翻车现场，附带可立即执行的诊断命令。

4.1 现象：梯度grad_c出现大面积 NaN 或 Inf，fmincon直接崩溃

原因：正演中c出现零或负值 →c²为负 → 波动方程系数为负 → 数值不稳定 →p爆炸 →q积分发散 →q * ∂²p/∂t²产生 NaN
解决：

• 预防：在fwi_objective开头强制裁剪c：c = max(c, 1500);（下界取地表水层速度）
• 诊断：在forward_fd_staggered.m结尾加assert(all(c(:) > 0), 'c has non-positive value!');
• 急救：若已发生，加载c模型文件，运行c(c<=0) = 1500; save c_fixed.mat c;后重启反演

4.2 现象：反演收敛极慢，fval下降停滞在 1e-1，gradNorm始终 > 1e-2

原因：PML 吸收层失效，边界反射混入有效信号 →d_syn与d_obs残差含强相干噪声 → 梯度指向错误方向
解决：

• 验证 PML：单独运行forward_fd_staggered，输入一个脉冲源，imshow(p(:,:,end))查看边界是否干净。若见强反射，调大damp_x/damp_z幅度（原公式乘 1.5）
• 参数修正：setup_pml_2d中alpha_max从2*pi*f_max改为3*pi*f_max，kappa_max从1改为3
• 终极手段：在d_obs和d_syn中，对每一道数据d(t)执行d = d .* hamming(length(d))'（汉宁窗），压制边界反射能量

4.3 现象：反演结果出现“条纹状”伪影，平行于网格方向，且随迭代次数增强

原因：正演与伴随的离散格式不对称（如正演用 4 阶 FD，伴随用 2 阶），破坏伴随原理⟨Lδm, d⟩ = ⟨δm, L* d⟩→ 梯度不准确
解决：

• 强制统一：检查forward_fd_staggered.m中div_v计算与gradient_adjoint.m中laplacian2d是否同阶同 stencil。本系统要求二者均为 4 阶：

  % 4阶拉普拉斯（正确） lap_q = (-1/12)*q(i-2,j) + (4/3)*q(i-1,j) - (5/2)*q(i,j) + (4/3)*q(i+1,j) - (1/12)*q(i+2,j) ... + (-1/12)*q(i,j-2) + (4/3)*q(i,j-1) - (5/2)*q(i,j) + (4/3)*q(i,j+1) - (1/12)*q(i,j+2);

• 验证对称性：写测试脚本，随机生成δm，计算Lδm和L*(Lδm)，检查norm(L*(Lδm) - δm)是否 < 1e-10

4.4 现象：多尺度反演中，加入高频后fval突然增大 10 倍，模型c出现剧烈震荡

原因：高频数据信噪比（SNR）过低，d_obs中噪声被当作有效信号拟合 → 过拟合
解决：

• SNR 量化：对每炮数据，计算snr_db = 20*log10(norm(d_signal)/norm(d_noise))，其中d_noise取首 50ms（震源激发前）
• 阈值过滤：若snr_db < 10，该炮数据权重设为 0：weight_shot = (snr_db > 10);
• 数据加权：在fwi_objective中，残差改为sum(weight_shot .* (d_obs - d_syn).^2)
• 后悔药：保存每级反演的c模型（save c_scale1.mat c;），若高频失败，回退到上一级c_scale2.mat并改用fmincon的'HessianMultiplyFcn'加入更强正则

5. 验证与可视化：用三类图证明你的 FWI 结果可信

FWI 成果不能只看fval下降曲线。必须通过物理一致性、数据拟合度、地质合理性三重验证。MATLAB 的强可视化能力在此刻转化为生产力。

5.1 波场残差时空图：诊断拟合质量的黄金标准

这不是简单的imagesc(d_obs - d_syn)，而是按炮、按时间窗、按频率成分的精细化分析。

% residual_analysis.m figure('Position', [100,100,1200,800]); subplot(2,2,1); imagesc(d_obs(1:500,:)); title('Observed (first 500ms)'); axis xy; subplot(2,2,2); imagesc(d_syn(1:500,:)); title('Synthetic (first 500ms)'); axis xy; subplot(2,2,3); imagesc((d_obs - d_syn)(1:500,:)); title('Residual (first 500ms)'); axis xy; subplot(2,2,4); % 计算残差 RMS 随时间变化 rms_res = sqrt(mean((d_obs - d_syn).^2, 2)); plot(rms_res); title('RMS Residual vs Time'); xlabel('Time Sample'); ylabel('RMS');

关键判据：

• Residual图应呈“白噪声”状，无清晰相干事件（如未建模的多次波、面波）
• RMS Residual曲线应在t=0（震源时刻）达峰值，之后指数衰减；若在t=300ms后仍持平，说明深层反射未拟合
• 若Residual中出现与Observed形状相似但符号相反的斑块，说明模型低估了某处速度

5.2 梯度能量分布图：验证反演聚焦性

梯度grad_c的空间分布揭示反演“注意力”所在。理想情况下，能量应集中在构造变化剧烈区（断层、尖灭），而非均匀散布。

% gradient_energy.m grad_abs = abs(grad_c); grad_norm = grad_abs / max(grad_abs(:)); % 归一化到 [0,1] figure; subplot(1,2,1); imagesc(grad_abs); title('Gradient Magnitude'); colorbar; subplot(1,2,2); % 计算梯度能量集中度：90% 能量覆盖面积占比 grad_sorted = sort(grad_norm(:), 'descend'); cum_energy = cumsum(grad_sorted); idx_90 = find(cum_energy >= 0.9, 1, 'first'); area_90 = idx_90 / numel(grad_norm); imagesc(grad_norm > grad_sorted(idx_90)); title(sprintf('Top 90%% Energy Area: %.1f%%', area_90*100)); colorbar;

地质解读：

• area_90 < 15%：梯度高度聚焦，反演有效（如断层两侧）
• area_90 > 40%：梯度弥散，说明数据信息不足或初始模型偏差太大，需加强低频或增加正则
• 若Gradient Magnitude图中出现与已知井位vp_log位置强相关的亮点，是重大利好信号

5.3 井约束剖面对比：用实测数据一票否决

最终交付物必须与测井vp曲线对齐。MATLAB 中实现“模型-井”空间匹配是核心技能。

% well_tie.m % 加载测井数据（depth_m, vp_log） load('well_vp.mat'); % depth_m: (n,1), vp_log: (n,1) % 将井坐标 (x_well, z_well) 映射到模型网格 [i_well, j_well] = world2grid([z_well, x_well], dx, dz, oz, ox); % 注意 z/x 顺序 % 提取模型上该列的 vp_profile vp_model = c(round(i_well), round(j_well)); % 插值更准，但此为快速验证 % 绘制对比 figure; plot(vp_log, depth_m, 'r', 'LineWidth', 2); hold on; plot(vp_model, depth_m, 'b--', 'LineWidth', 2); xlabel('Velocity (m/s)'); ylabel('Depth (m)'); legend('Well Log', 'FWI Model'); grid on; title(sprintf('Well Tie at X=%.0f m', x_well));

验收红线：

• 在目的层（如储层顶底）深度误差< 5m
• 速度绝对误差< 150 m/s（碳酸盐岩可放宽至 200）
• 若vp_model在浅层（<100m）系统性偏低，说明近地表模型不准，需单独反演近地表
• 若vp_model在深层（>2000m）持续高于测井，检查c上边界约束ub是否设得太松

6. 进阶技巧：用 MATLAB OOP 架构重构 FWI 系统，实现算法热插拔与多场景复用

当你的 FWI 系统要支撑“海上宽频带数据”、“陆上微震监测”、“CO₂ 封存时移成像”多个项目时，面向过程的.m文件堆砌会迅速失控。MATLAB 的 OOP（Object-Oriented Programming）不是炫技，而是工程化刚需——它让你把正演引擎、梯度计算器、优化器、数据加载器解耦为可独立测试、可继承扩展的类。

6.1 FWI 系统的四层 OOP 架构设计

本文还有配套的精品资源，点击获取