时域二维声波方程全波形反演（FWI）MATLAB程序

一、全波形反演原理概述

全波形反演（Full Waveform Inversion, FWI）是一种利用地震波场的全部信息反演地下介质参数的高精度成像技术。其核心是通过最小化模拟数据与观测数据的差异，迭代更新地下介质参数（如速度模型）。

声波方程（二维）：

其中 $p$为压力场，$c(x)$为声波速度，$∇2$为拉普拉斯算子，$f$为震源函数。

反演目标：

其中 $p_{cal}$为模拟数据，$p_{obs}$为观测数据。

二、MATLAB程序实现

1. 主程序框架

% 时域二维声波方程全波形反演
clear; clc; close all;

%% 参数设置
% 模型参数
nx = 100; ny = 100;        % 网格点数
dx = 10; dy = 10;          % 网格间距 (m)
dz = 10;                   % 深度间距 (m)
c_true = initialModel(nx, ny); % 真实速度模型
c_init = smoothModel(c_true);  % 初始速度模型

% 时间参数
dt = 0.001;                % 时间步长 (s)
nt = 2000;                 % 时间步数
t = (0:nt-1)*dt;           % 时间轴

% 震源参数
f0 = 25;                   % 主频 (Hz)
src_x = nx/2; src_y = 5;    % 震源位置
src_func = ricker(f0, t);  % Ricker子波

% 接收器参数
rec_x = 1:nx; rec_y = ny/2*ones(1,nx); % 接收器排列
nrec = length(rec_x);       % 接收器数量

% 反演参数
niter = 20;                % 反演迭代次数
alpha = 0.0001;             % 步长
c = c_init;                 % 当前速度模型

%% 生成观测数据（模拟）
d_obs = forwardModel(c_true, src_x, src_y, src_func, rec_x, rec_y, dx, dy, dt, nt);

%% 全波形反演迭代
for iter = 1:niter
    % 正演模拟
    d_cal = forwardModel(c, src_x, src_y, src_func, rec_x, rec_y, dx, dy, dt, nt);
    
    % 计算残差
    residual = d_cal - d_obs;
    
    % 计算梯度
    grad = gradientCalculation(c, residual, src_func, rec_x, rec_y, dx, dy, dt, nt);
    
    % 预处理（高斯平滑）
    grad_smooth = smoothGradient(grad, 3);
    
    % 更新模型 (梯度下降法)
    c_new = c - alpha * grad_smooth ./ (smoothModel(c) + 1e-6);
    
    % 更新模型并约束速度范围
    c = constrainVelocity(c_new, 1500, 4500);
    
    % 显示进度
    misfit(iter) = 0.5 * norm(residual(:))^2;
    fprintf('Iteration %d: Misfit = %.4f\n', iter, misfit(iter));
    
    % 可视化
    if mod(iter, 5) == 0
        visualizeModels(c_true, c, misfit(iter));
    end
end

%% 最终结果可视化
visualizeFinalResults(c_true, c, d_obs, d_cal, misfit);

2. 声波方程正演模拟

function p = forwardModel(c, src_x, src_y, src_func, rec_x, rec_y, dx, dy, dt, nt)
    % 声波方程正演模拟 (二阶时间差分)
    % 输入:
    %   c: 速度模型 (nx×ny)
    %   src_x, src_y: 震源位置
    %   src_func: 震源时间序列
    %   rec_x, rec_y: 接收器位置
    %   dx, dy: 空间步长
    %   dt: 时间步长
    %   nt: 时间步数
    % 输出:
    %   p: 接收器记录的波形 (nt×nrec)
    
    [nx, ny] = size(c);
    p = zeros(nx, ny, 3); % 压力场 (存储三个时间层)
    p_rec = zeros(nt, length(rec_x)); % 接收器记录
    
    % 二阶空间差分系数
    coeff_x = [1, -2, 1]/(dx^2);
    coeff_y = [1, -2, 1]/(dy^2);
    
    % 时间步进
    for it = 1:nt
        % 更新压力场 (中心差分)
        p_new = zeros(nx, ny);
        for ix = 2:nx-1
            for iy = 2:ny-1
                laplacian = coeff_x(1)*p(ix-1,iy,:) + coeff_x(2)*p(ix,iy,:) + coeff_x(3)*p(ix+1,iy,:) + ...
                           coeff_y(1)*p(ix,iy-1,:) + coeff_y(2)*p(ix,iy,:) + coeff_y(3)*p(ix,iy+1,:);
                p_new(ix,iy) = 2*p(ix,iy,2) - p(ix,iy,1) + (dt^2)*(laplacian/c(ix,iy)^2);
            end
        end
        
        % 添加震源项 (位于src_x, src_y)
        p_new(src_x, src_y) = p_new(src_x, src_y) + dt^2 * src_func(it);
        
        % 边界吸收条件 (简单实现)
        p_new(1,:) = 0; p_new(end,:) = 0;
        p_new(:,1) = 0; p_new(:,end) = 0;
        
        % 更新时间层
        p(:,:,1) = p(:,:,2); % t-1
        p(:,:,2) = p(:,:,3); % t
        p(:,:,3) = p_new;    % t+1
        
        % 记录接收器数据
        for ir = 1:length(rec_x)
            p_rec(it, ir) = p(rec_x(ir), rec_y(ir), 3);
        end
    end
    
    p = p_rec;
end

3. 梯度计算（伴随状态法）

function grad = gradientCalculation(c, residual, src_func, rec_x, rec_y, dx, dy, dt, nt)
    % 使用伴随状态法计算梯度
    % 输入:
    %   c: 当前速度模型
    %   residual: 残差 (nt×nrec)
    %   其他参数同forwardModel
    % 输出:
    %   grad: 梯度 (nx×ny)
    
    [nx, ny] = size(c);
    grad = zeros(nx, ny);
    
    % 正向传播存储压力场
    [p_forward, ~] = forwardStorage(c, src_func, dx, dy, dt, nt);
    
    % 反向传播伴随场
    lambda = zeros(nx, ny, 3); % 伴随场 (存储三个时间层)
    
    % 初始化伴随源 (接收器残差)
    adjoint_src = zeros(nx, ny, nt);
    for ir = 1:length(rec_x)
        adjoint_src(rec_x(ir), rec_y(ir), :) = residual(:, ir);
    end
    
    % 时间反向传播
    for it = nt:-1:1
        % 更新伴随场
        lambda_new = zeros(nx, ny);
        for ix = 2:nx-1
            for iy = 2:ny-1
                laplacian = (lambda(ix-1,iy,2) - 2*lambda(ix,iy,2) + lambda(ix+1,iy,2))/(dx^2) + ...
                            (lambda(ix,iy-1,2) - 2*lambda(ix,iy,2) + lambda(ix,iy+1,2))/(dy^2);
                lambda_new(ix,iy) = 2*lambda(ix,iy,2) - lambda(ix,iy,1) + (dt^2)*laplacian;
            end
        end
        
        % 添加伴随源 (负号来自伴随算子)
        lambda_new = lambda_new - adjoint_src(:,:,it);
        
        % 边界条件
        lambda_new(1,:) = 0; lambda_new(end,:) = 0;
        lambda_new(:,1) = 0; lambda_new(:,end) = 0;
        
        % 更新时间层
        lambda(:,:,1) = lambda(:,:,2);
        lambda(:,:,2) = lambda(:,:,3);
        lambda(:,:,3) = lambda_new;
        
        % 累加梯度 (在时间步it)
        if it > 1
            grad = grad + (2*p_forward(:,:,it).*lambda(:,:,3))/c.^3 * dt^2;
        end
    end
end

4. 辅助函数

% Ricker子波生成
function w = ricker(f0, t)
    w = (1 - 2*(pi*f0*t).^2) .* exp(-(pi*f0*t).^2);
end

% 初始模型生成 (层状模型)
function model = initialModel(nx, ny)
    model = 2000 * ones(nx, ny); % 背景速度 2000 m/s
    model(40:60, :) = 2500;      % 高速层
    model(70:80, 30:70) = 3000;  % 高速体
    model = model + 100*randn(nx, ny); % 添加随机扰动
end

% 模型平滑
function model_smooth = smoothModel(model)
    model_smooth = imgaussfilt(model, 2);
end

% 梯度平滑
function grad_smooth = smoothGradient(grad, sigma)
    grad_smooth = imgaussfilt(grad, sigma);
end

% 速度约束
function c_constrained = constrainVelocity(c, vmin, vmax)
    c_constrained = max(min(c, vmax), vmin);
end

% 正演存储 (用于梯度计算)
function [p_store, p_rec] = forwardStorage(c, src_func, dx, dy, dt, nt)
    [nx, ny] = size(c);
    p = zeros(nx, ny, 3);
    p_store = zeros(nx, ny, nt);
    p_rec = zeros(nt, nx); % 假设接收器在y=1
    
    for it = 1:nt
        p_new = zeros(nx, ny);
        for ix = 2:nx-1
            for iy = 2:ny-1
                laplacian = (p(ix-1,iy,3) - 2*p(ix,iy,3) + p(ix+1,iy,3))/(dx^2) + ...
                            (p(ix,iy-1,3) - 2*p(ix,iy,3) + p(ix,iy+1,3))/(dy^2);
                p_new(ix,iy) = 2*p(ix,iy,2) - p(ix,iy,1) + (dt^2)*laplacian/c(ix,iy)^2;
            end
        end
        
        % 震源 (中心位置)
        src_x = nx/2; src_y = ny/2;
        p_new(src_x, src_y) = p_new(src_x, src_y) + dt^2 * src_func(it);
        
        % 边界条件
        p_new(1,:) = 0; p_new(end,:) = 0;
        p_new(:,1) = 0; p_new(:,end) = 0;
        
        % 更新时间层
        p(:,:,1) = p(:,:,2);
        p(:,:,2) = p(:,:,3);
        p(:,:,3) = p_new;
        
        % 存储当前压力场
        p_store(:,:,it) = p(:,:,3);
        
        % 记录接收器数据
        p_rec(it, :) = p(:,1,3)';
    end
end

5. 可视化函数

function visualizeModels(c_true, c_current, misfit)
    figure(1);
    subplot(1,2,1);
    imagesc(c_true'); colorbar; title('真实速度模型');
    xlabel('X (m)'); ylabel('Y (m)');
    
    subplot(1,2,2);
    imagesc(c_current'); colorbar; 
    title(sprintf('反演速度模型 (Misfit=%.4f)', misfit));
    xlabel('X (m)'); ylabel('Y (m)');
    
    drawnow;
end

function visualizeFinalResults(c_true, c_inv, d_obs, d_cal, misfit)
    % 速度模型对比
    figure(1);
    subplot(1,3,1);
    imagesc(c_true'); colorbar; title('真实速度模型');
    xlabel('X (m)'); ylabel('Y (m)');
    
    subplot(1,3,2);
    imagesc(c_inv'); colorbar; title('反演速度模型');
    xlabel('X (m)'); ylabel('Y (m)');
    
    subplot(1,3,3);
    imagesc(c_true' - c_inv'); colorbar; title('模型差异');
    xlabel('X (m)'); ylabel('Y (m)');
    
    % 数据对比
    figure(2);
    subplot(2,1,1);
    plot(d_obs(:,1), 'b'); hold on;
    plot(d_cal(:,1), 'r');
    legend('观测数据', '反演数据');
    title('接收器1波形对比');
    xlabel('时间 (s)'); ylabel('振幅');
    
    subplot(2,1,2);
    plot(misfit);
    title('目标函数收敛曲线');
    xlabel('迭代次数'); ylabel('Misfit');
    
    % 误差分布
    figure(3);
    imagesc((c_true - c_inv).^2'); colorbar;
    title('反演误差平方分布');
    xlabel('X (m)'); ylabel('Y (m)');
end

参考代码 时域2维声波方程全波形反演程序 www.youwenfan.com/contentcss/51732.html

三、关键技术说明

1. 正演模拟

• 采用二阶时间差分、二阶空间差分的显式有限差分法
• 使用Ricker子波作为震源函数
• 实现简单吸收边界条件

2. 梯度计算

• 使用伴随状态法高效计算梯度
• 避免直接求导的数值不稳定性
• 通过正演存储减少计算量

3. 反演优化

• 采用梯度下降法更新模型
• 添加高斯平滑正则化项
• 约束速度范围保证物性合理

4. 模型构建

• 层状背景模型+高速异常体
• 添加随机扰动模拟实际地质情况
• 多尺度反演策略（未在代码中体现，实际中重要）

四、应用扩展

1. 弹性波FWI：

   % 添加剪切波速度
   vs = 0.5 * vp; % 泊松固体假设
   

2. 各向异性FWI：

   % 添加各向异性参数
   epsilon = 0.1; delta = 0.2;
   

3. 多尺度反演：

   % 从低频到高频逐步反演
   freq_bands = [5, 10, 20, 30]; % Hz
   

4. 粘滞性补偿：

   % 添加品质因子Q
   Q = 50; % 品质因子
   attenuation = exp(-dt*pi*freq/Q);
   

六、总结

本程序实现了时域二维声波方程全波形反演的核心算法：

1. 采用有限差分法求解声波方程
2. 使用伴随状态法计算梯度
3. 实现梯度下降优化算法
4. 包含模型构建、数据模拟和可视化模块