二维经验模态分解(BEMD)算法详解与MATLAB实现

二维经验模态分解(BEMD)是近年来新兴的一种适用于非线性、非平稳数据的自适应数据分析方法，它是一维经验模态分解(EMD)在二维空间的扩展。BEMD能够自适应地将复杂图像分解为不同尺度的本征模态函数(BIMFs)和残余分量，在图像处理、模式识别和特征提取等领域有广泛应用。

算法原理

BEMD的核心思想是通过迭代筛选过程将图像分解为多个不同尺度的本征模态函数(BIMFs)和一个残余分量：

  1. 初始化：令残余信号\(R_0(x,y) = I(x,y)\)（原始图像）
  2. 提取第k个BIMF：
    a. 初始化\(h_0(x,y) = R_{k-1}(x,y)\)
    b. 识别\(h_j(x,y)\)的所有局部极大值和局部极小值点
    c. 构造上包络面\(U_j(x,y)\)和下包络面\(L_j(x,y)\)
    d. 计算平均包络面：\(m_j(x,y) = \frac{U_j(x,y) + L_j(x,y)}{2}\)
    e. 更新：\(h_{j+1}(x,y) = h_j(x,y) - m_j(x,y)\)
    f. 重复b-e直到\(h_{j+1}\)满足BIMF条件
    g. 设置\(BIMF_k = h_{j+1}\)
  3. 更新残余：\(R_k(x,y) = R_{k-1}(x,y) - BIMF_k\)
  4. 重复步骤2直到残余满足停止条件

MATLAB实现

function [bimfs, residual] = bemd(img, num_bimfs, max_iter, tol)
% 二维经验模态分解(BEMD)实现
% 输入:
%   img - 输入图像(M×N矩阵)
%   num_bimfs - 提取的BIMF数量
%   max_iter - 每个BIMF的最大筛选迭代次数
%   tol - 筛选停止容差
% 输出:
%   bimfs - 分解得到的BIMFs(三维矩阵M×N×K)
%   residual - 残余分量

% 参数检查与初始化
if nargin < 2, num_bimfs = 4; end
if nargin < 3, max_iter = 20; end
if nargin < 4, tol = 0.05; end

[M, N] = size(img);
residual = double(img);
bimfs = zeros(M, N, num_bimfs);

% 循环提取每个BIMF
for k = 1:num_bimfs
    fprintf('正在提取BIMF %d/%d...\n', k, num_bimfs);
    h = residual;
    
    % 筛选过程
    for iter = 1:max_iter
        % 检测局部极值点
        [max_pts, min_pts] = find_extrema(h);
        
        % 构造上包络面(极大值点插值)
        if size(max_pts, 1) < 4
            warning('极值点不足，终止筛选');
            break;
        end
        upper_env = scatteredInterpolant(max_pts(:,2), max_pts(:,1), max_pts(:,3), 'natural', 'linear');
        
        % 构造下包络面(极小值点插值)
        if size(min_pts, 1) < 4
            warning('极值点不足，终止筛选');
            break;
        end
        lower_env = scatteredInterpolant(min_pts(:,2), min_pts(:,1), min_pts(:,3), 'natural', 'linear');
        
        % 创建网格坐标
        [X, Y] = meshgrid(1:N, 1:M);
        
        % 计算包络面
        U = upper_env(X, Y);
        L = lower_env(X, Y);
        
        % 处理边界点(防止NaN)
        U(isnan(U)) = h(isnan(U));
        L(isnan(L)) = h(isnan(L));
        
        % 计算平均包络
        m = (U + L) / 2;
        
        % 更新h
        h_prev = h;
        h = h - m;
        
        % 计算停止准则
        sd = sum((h_prev(:) - h(:)).^2) / sum(h_prev(:).^2);
        if sd < tol
            fprintf('  筛选迭代 %d: SD = %.4f < %.4f, 停止筛选\n', iter, sd, tol);
            break;
        end
        
        % 显示当前筛选信息
        if mod(iter, 5) == 0
            fprintf('  筛选迭代 %d: SD = %.4f\n', iter, sd);
        end
    end
    
    % 存储当前BIMF
    bimfs(:, :, k) = h;
    
    % 更新残余
    residual = residual - h;
    
    % 如果残余变化很小，停止分解
    if std(residual(:)) < 0.05 * std(img(:))
        fprintf('残余变化很小，提前终止分解\n');
        break;
    end
end

% 显示分解结果
plot_bemd_results(img, bimfs, residual);
end

function [max_pts, min_pts] = find_extrema(img)
% 寻找图像的局部极大值和极小值
% 返回: max_pts - [行, 列, 值] 的极大值点矩阵
%       min_pts - [行, 列, 值] 的极小值点矩阵

[M, N] = size(img);
max_pts = [];
min_pts = [];

% 使用3x3邻域检测极值点
for i = 2:M-1
    for j = 2:N-1
        % 获取3x3邻域
        region = img(i-1:i+1, j-1:j+1);
        center = region(2, 2);
        
        % 检测极大值(大于所有邻居)
        if center > max(region(:)) && center > mean(region(:)) + 0.1*std(region(:))
            max_pts = [max_pts; i, j, center];
        end
        
        % 检测极小值(小于所有邻居)
        if center < min(region(:)) && center < mean(region(:)) - 0.1*std(region(:))
            min_pts = [min_pts; i, j, center];
        end
    end
end

% 添加边界点作为候选极值点
boundary_pts = [1, 1; 1, N; M, 1; M, N; 
                1, round(N/2); M, round(N/2); 
                round(M/2), 1; round(M/2), N];
for i = 1:size(boundary_pts, 1)
    r = boundary_pts(i, 1);
    c = boundary_pts(i, 2);
    max_pts = [max_pts; r, c, img(r, c)];
    min_pts = [min_pts; r, c, img(r, c)];
end
end

function plot_bemd_results(original, bimfs, residual)
% 可视化BEMD分解结果

num_bimfs = size(bimfs, 3);
figure('Position', [100, 100, 1200, 800], 'Name', '二维经验模态分解(BEMD)结果');

% 显示原始图像
subplot(3, num_bimfs+1, 1);
imshow(original, []);
title('原始图像');
axis image;

% 显示各个BIMF
for k = 1:num_bimfs
    subplot(3, num_bimfs+1, k+1);
    imshow(bimfs(:, :, k), []);
    title(sprintf('BIMF %d', k));
    axis image;
    
    % 显示BIMF的3D表面
    subplot(3, num_bimfs+1, num_bimfs+1 + k+1);
    surf(bimfs(:, :, k), 'EdgeColor', 'none');
    title(sprintf('BIMF %d 表面', k));
    axis tight;
    view(-30, 60);
    camlight;
    lighting gouraud;
end

% 显示残余分量
subplot(3, num_bimfs+1, num_bimfs+2);
imshow(residual, []);
title('残余分量');
axis image;

% 显示残余分量的3D表面
subplot(3, num_bimfs+1, 2*(num_bimfs+1));
surf(residual, 'EdgeColor', 'none');
title('残余分量表面');
axis tight;
view(-30, 60);
camlight;
lighting gouraud;

% 添加算法说明
annotation('textbox', [0.05, 0.01, 0.9, 0.05], 'String', ...
    '二维经验模态分解(BEMD): 一种自适应图像分解方法，适用于非线性、非平稳数据', ...
    'EdgeColor', 'none', 'HorizontalAlignment', 'center', ...
    'FontSize', 12, 'BackgroundColor', [0.9, 0.95, 1]);
end

应用示例：图像分解与特征提取

%% BEMD应用示例：图像分解与特征提取
clear; clc; close all;

% 读取测试图像
img = imread('cameraman.tif');
img = double(img);

% 设置BEMD参数
num_bimfs = 4;    % 提取4个BIMF
max_iter = 50;    % 每个BIMF最大迭代次数
tol = 0.03;       % 筛选停止容差

% 执行BEMD分解
[bimfs, residual] = bemd(img, num_bimfs, max_iter, tol);

%% 应用1：图像去噪
% 添加高斯噪声
noisy_img = img + 30 * randn(size(img));

% 对噪声图像进行BEMD分解
[noisy_bimfs, noisy_residual] = bemd(noisy_img, num_bimfs, max_iter, tol);

% 去噪策略：去除第一个BIMF（通常包含大部分噪声）
denoised_img = noisy_residual;
for k = 2:size(noisy_bimfs, 3)
    denoised_img = denoised_img + noisy_bimfs(:, :, k);
end

% 显示去噪结果
figure('Name', 'BEMD图像去噪', 'Position', [100, 100, 1200, 400]);
subplot(1,3,1); imshow(noisy_img, []); title('带噪声图像');
subplot(1,3,2); imshow(denoised_img, []); title('BEMD去噪结果');
subplot(1,3,3); imshow(img, []); title('原始图像');

% 计算PSNR
psnr_noisy = psnr(uint8(noisy_img), uint8(img));
psnr_denoised = psnr(uint8(denoised_img), uint8(img));
fprintf('去噪效果评估:\n');
fprintf('带噪声图像PSNR: %.2f dB\n', psnr_noisy);
fprintf('BEMD去噪后PSNR: %.2f dB\n', psnr_denoised);

%% 应用2：边缘检测
% 使用第一个BIMF进行边缘检测
edges = edge(bimfs(:, :, 1), 'canny');

% 显示边缘检测结果
figure('Name', '基于BIMF的边缘检测', 'Position', [100, 100, 1000, 400]);
subplot(1,3,1); imshow(img, []); title('原始图像');
subplot(1,3,2); imshow(bimfs(:, :, 1), []); title('BIMF 1');
subplot(1,3,3); imshow(edges, []); title('提取的边缘');

%% 应用3：纹理分析
% 使用不同BIMF的统计特征进行纹理分析
texture_features = zeros(num_bimfs, 4); % [均值, 方差, 熵, 平滑度]

for k = 1:num_bimfs
    bimf = bimfs(:, :, k);
    texture_features(k, 1) = mean(bimf(:));      % 均值
    texture_features(k, 2) = var(bimf(:));       % 方差
    texture_features(k, 3) = entropy(bimf);      % 熵
    texture_features(k, 4) = 1/(1 + var(bimf(:))); % 平滑度
end

% 显示纹理特征
figure('Name', '基于BEMD的纹理特征分析');
bar(texture_features);
title('不同BIMF的纹理特征');
xlabel('BIMF序号');
ylabel('特征值');
legend({'均值', '方差', '熵', '平滑度'}, 'Location', 'best');
grid on;

参考代码 二维经验模态分解(BEMD)是近年来新兴的一种适用于非线性、非平稳数据的自适应数据 www.youwenfan.com/contentcno/79508.html

BEMD算法特点与优势

  1. 自适应性：

    • 不需要预先设定基函数
    • 分解基于数据本身的局部特性
    • 适用于各种复杂的图像模式

  2. 多尺度分析能力：

    • BIMFs按频率从高到低排列
    • 不同BIMF捕获不同尺度的特征
    • 残余分量表示图像的整体趋势

  3. 非线性处理能力：

    • 能够处理非平稳信号
    • 有效分析非线性系统产生的数据
    • 保留信号的瞬时频率特性

  4. 边界处理改进：

    • 镜像延拓减少边界效应
    • 特殊处理边界极值点
    • 自适应插值方法

结论

二维经验模态分解(BEMD)作为一种自适应、多尺度的图像分解方法，在非线性、非平稳数据处理方面展现出独特优势。通过MATLAB实现，我们可以直观地观察BEMD分解过程及其在图像处理中的应用效果。随着算法优化和计算效率的提升，BEMD有望在更多领域发挥重要作用，特别是在需要多尺度特征分析的复杂场景中。

BEMD的核心价值在于它能够根据数据本身的特性进行自适应分解，无需预先假设数据分布或选择基函数，这使得它在处理现实世界中的复杂数据时具有显著优势。随着研究的深入和应用场景的拓展，BEMD将继续发展并完善，成为非线性数据分析的重要工具之一。

posted @ 2025-12-25 17:43  chen_yig  阅读(25)  评论(0)    收藏  举报