MATLAB实现二维变分模态分解(VMD)用于地震数据去噪

变分模态分解(Variational Mode Decomposition, VMD)是一种自适应信号分解方法，可以有效地将信号分解为多个本征模态函数(IMF)。对于二维数据(如地震数据)，VMD可以用于去噪和特征提取。

1. 二维VMD算法原理

二维VMD通过解决以下约束变分问题来分解信号：

min_{u_k,ω_k} { ∑_k ||∇[ (δ(t) + j/πt) * u_k(t) ] e^{-jω_kt} ||^2_2 }
s.t. ∑_k u_k = f

其中：

• u_k是第k个模态函数
• ω_k是第k个模态的中心频率
• f是输入信号

2. MATLAB实现二维VMD

二维VMD的核心代码：

function [u, u_hat, omega] = VMD2D(f, alpha, tau, K, DC, init, tol)
    % 二维变分模态分解
    % 输入:
    %   f - 待分解的二维信号
    %   alpha - 带宽约束参数
    %   tau - 时间步长(噪声容忍)
    %   K - 分解的模态数
    %   DC - 是否包含直流分量(1包含，0不包含)
    %   init - 初始化方式 (1=频谱域, 2=时域均匀分布)
    %   tol - 收敛容差
    % 输出:
    %   u - 分解后的模态(K×M×N数组)
    %   u_hat - 模态的频谱
    %   omega - 中心频率
    
    % 参数检查与默认值设置
    if nargin < 7, tol = 1e-6; end
    if nargin < 6, init = 1; end
    if nargin < 5, DC = 0; end
    
    % 信号尺寸
    [M, N] = size(f);
    
    % 频率网格
    [fx, fy] = meshgrid(1:N, 1:M);
    fx = fx - floor(N/2) - 1;
    fy = fy - floor(M/2) - 1;
    
    % 频率坐标
    omega_x = 2*pi*fx/N;
    omega_y = 2*pi*fy/M;
    
    % 对称扩展
    f_hat = fftshift(fft2(f));
    
    % 初始化
    u_hat = zeros(K, M, N);
    if init == 1
        % 基于频谱的初始化
        for k = 1:K
            u_hat(k,:,:) = f_hat/K;
        end
    else
        % 均匀分布初始化
        for k = 1:K
            u_hat(k,:,:) = rand(M,N).*f_hat;
        end
    end
    
    % 中心频率初始化
    omega = zeros(K, 2);
    for k = 1:K
        omega(k,:) = [0, 0];
    end
    
    % 拉格朗日乘子初始化
    lambda_hat = zeros(M, N);
    
    % 迭代
    n = 0;
    u_diff = tol + eps;
    
    while (u_diff > tol && n < 500)
        n = n + 1;
        
        % 更新模态
        for k = 1:K
            % 计算梯度
            gradient = (omega_x - omega(k,1)).^2 + (omega_y - omega(k,2)).^2;
            
            % 更新频谱
            sum_uk = squeeze(sum(u_hat([1:k-1 k+1:K],:,:), 1));
            u_hat(k,:,:) = (f_hat - sum_uk + lambda_hat/2) ./ (1 + alpha*gradient);
            
            % 更新中心频率
            if DC == 0 || k > 1
                intensity = abs(squeeze(u_hat(k,:,:))).^2;
                omega(k,1) = sum(sum(omega_x .* intensity)) / sum(sum(intensity));
                omega(k,2) = sum(sum(omega_y .* intensity)) / sum(sum(intensity));
            end
        end
        
        % 更新拉格朗日乘子
        sum_uk_hat = squeeze(sum(u_hat, 1));
        lambda_hat = lambda_hat + tau*(f_hat - sum_uk_hat);
        
        % 计算收敛条件
        if n > 1
            u_diff = norm(squeeze(sum(u_hat,1)) - sum_uk_prev, 'fro') / norm(sum_uk_prev, 'fro');
        end
        sum_uk_prev = sum_uk_hat;
    end
    
    % 逆变换得到时域模态
    u = zeros(K, M, N);
    for k = 1:K
        u(k,:,:) = real(ifft2(ifftshift(squeeze(u_hat(k,:,:)))));
    end
end

参考代码 对二维数据进行变分模态分解，可用于地震数据去噪等

3. 地震数据去噪应用

% 加载地震数据
load('seismic_data.mat'); % 假设数据存储在变量seismic_data中

% 数据预处理
data = double(seismic_data);
[M, N] = size(data);

% 参数设置
alpha = 2000;   % 带宽约束
tau = 0.1;      % 噪声容忍
K = 4;          % 模态数
tol = 1e-6;     % 收敛容差

% 执行二维VMD分解
[u, ~, ~] = VMD2D(data, alpha, tau, K, 0, 1, tol);

% 显示分解结果
figure;
for k = 1:K
    subplot(2,2,k);
    imagesc(squeeze(u(k,:,:)));
    title(['IMF ', num2str(k)]);
    colormap('gray');
    axis off;
end

% 去噪处理(假设噪声主要在高阶IMF中)
denoised_data = squeeze(sum(u(1:K-1,:,:), 1));

% 显示去噪结果
figure;
subplot(1,2,1);
imagesc(data);
title('原始地震数据');
colormap('gray');
axis off;

subplot(1,2,2);
imagesc(denoised_data);
title('去噪后数据');
colormap('gray');
axis off;

4. 性能评估

% 计算信噪比(SNR)和峰值信噪比(PSNR)
function [snr, psnr] = evaluate_denoising(original, noisy, denoised)
    % 计算噪声水平
    noise = original - noisy;
    noise_power = mean(noise(:).^2);
    
    % 计算去噪后的噪声
    residual = original - denoised;
    residual_power = mean(residual(:).^2);
    
    % 计算SNR和PSNR
    signal_power = mean(original(:).^2);
    snr = 10*log10(signal_power / noise_power);
    psnr_noisy = 10*log10(max(original(:))^2 / noise_power);
    psnr_denoised = 10*log10(max(original(:))^2 / residual_power);
    
    fprintf('原始噪声SNR: %.2f dB\n', snr);
    fprintf('去噪后PSNR: %.2f dB (原始PSNR: %.2f dB)\n', psnr_denoised, psnr_noisy);
end

% 使用示例
% [snr, psnr] = evaluate_denoising(clean_data, noisy_data, denoised_data);

5. 参数优化建议

1. 模态数K：通常选择3-6个模态，过多会导致过分解
2. 带宽约束α：控制模态带宽，值越大模态越紧凑
3. 噪声容忍τ：影响收敛速度，通常设为0.1-0.5
4. 收敛容差tol：通常设为1e-6到1e-5

可以通过网格搜索或优化算法寻找最佳参数组合：

% 参数优化示例
alpha_values = [1000, 2000, 5000];
K_values = [3, 4, 5];
best_psnr = 0;

for alpha = alpha_values
    for K = K_values
        [u, ~, ~] = VMD2D(data, alpha, 0.1, K, 0, 1, 1e-6);
        denoised = squeeze(sum(u(1:K-1,:,:), 1));
        [~, current_psnr] = evaluate_denoising(clean_data, data, denoised);
        
        if current_psnr > best_psnr
            best_psnr = current_psnr;
            best_alpha = alpha;
            best_K = K;
        end
    end
end

fprintf('最佳参数: alpha=%.0f, K=%d, PSNR=%.2f dB\n', best_alpha, best_K, best_psnr);

对于大规模地震数据，可以考虑以下优化：

1. 分块处理大数据
2. 使用并行计算(parfor)
3. 结合GPU加速