BEMD图像分解实战：MATLAB实现原理与参数调优指南

简介：这是一份基于MATLAB实现的二维经验模态分解（BEMD）算法资源包，面向数值算法与人工智能领域的研究者、工程师及高校师生，适用于非线性、非平稳信号及多变量耦合数据的分析与处理，在环境、金融、生物医学等场景均有应用。压缩包共39个文件，以34个m脚本为主体，涵盖BEMD核心分解函数、曲面插值、极值点提取、包络构造等关键模块；另有4张lena测试图像与1份license说明文档，便于直接运行验证算法效果。整体约221KB，结构紧凑，已有274人学习下载。通过研读运行这套代码，读者可掌握从信号预处理、希尔伯特变换、IMF分解到耦合分析、信号重构与可视化的完整流程，理解二维分解与传统EMD的差异，还可将其中sift、gridfit、extrema等基础工具直接复用，灵活迁移至其他数据分析、模式识别任务中。

1. BEMD是什么：一个把图像拆成“纹理层+残差”的MATLAB工具箱入口

如果你手头有一张带噪声的纹理图像，想把它拆成一层一层的“细节”和一块“基底”，BEMD（二维经验模态分解）可能是你最快的答案。它不是深度神经网络，也不依赖训练数据，而是靠数值筛选把图像分解成若干个二维固有模态函数（BIMF）和残差。这个名叫“BEMD.zip”的MATLAB压缩包，常见于数值算法课程和人工智能大作业里，本质是一套用插值和筛选实现的非线性自适应分解工具。适合谁？打算做图像去噪、纹理特征提取、多聚焦图像融合，又不想从头写有限差分或小波栈的工程师。它会告诉你什么是“数据驱动的分解”，也可能让你踩进边界延拓和计算量的坑。

2. BEMD的数值算法原理：从EMD到二维筛选，为什么用插值而不是神经网络

2.1 一维EMD筛选的回顾：包络、均值、停止准则

BEMD的前身是Huang等人在1998年提出的经验模态分解（EMD）。给定一个一维信号，EMD把它拆成有限个本征模态函数（IMF），每个IMF代表信号中的一个振荡模式。筛选过程是递归的：先找到信号所有局部极大值和局部极小值，用三次样条分别插值出上包络和下包络，计算两条包络的平均值m，然后让原始信号减去m，得到候选IMF。如果候选不满足停止准则，就以它为输入继续筛，直到满足为止。之后从剩余信号中提取下一个IMF，直到剩余信号变成单调趋势或很小。

这里有两个细节值得记牢。第一，极值点不是简单的“找一个点”，而是要处理连续相等值的平台区。如果信号在一段时间内完全相同，那么平台两端都可以视为极值点，但也有实现取平台中心点。选错了，包络会出现台阶，IMF里会出现伪频率。第二，端点效应特别头疼。三次样条插值在端点处没有数据约束，包络往往会翘头或下沉，导致分解出来的前几个IMF在端点区域失真。我通常会对信号做两端的镜像延拓，延拓长度约为信号长度的5%，筛完再裁掉。

停止准则最常见的是一维的标准差准则：

SD = sum((h_{k-1} - h_k).^2) / sum(h_{k-1}.^2)

当SD小于某个阈值时停止，典型的阈值是0.2到0.3。但这是Huang经典论文里的建议，实际工程中要按噪声水平调整。如果信号本身很干净，SD可以设小一点；如果噪声大，SD太小会把噪声也筛成一个“IMF”，结果看着就像把噪声单独抽出来了，其实那并不是你想要的。

下面是MATLAB示意代码，展示一维EMD的单次筛选核心循环，但注意它省略了端点延拓和平台区处理，真实实现里这些不能省：

function imf = emd_sift_demo(x, maxIter, sdThresh) % 一维EMD单次筛选示例 % 输入：x为信号列向量，maxIter为最大筛选次数，sdThresh为SD停止阈值 % 输出：一个IMF分量 h = x(:); [pks, locsMax] = findpeaks(h); % 局部极大值幅值和位置 [vls, locsMin] = findpeaks(-h); % 局部极小值的相反数及位置 for k = 1:maxIter envUp = spline([1; locsMax], pks, 1:length(h))'; % 上包络 envLo = spline([1; locsMin], -vls, 1:length(h))'; % 下包络 means = (envUp + envLo) / 2; newH = h - means; sd = sum((newH - h).^2) / sum(h.^2); if sd < sdThresh break; end h = newH; end imf = h; end

逻辑说明：每次筛选都减去上下包络的均值，让信号逐渐趋向于局部对称。findpeaks的第一个输出是峰值幅值，第二个是峰值位置，初学者常把顺序搞反，导致spline的节点和值对应不上。参数说明：maxIter通常取10到20，sdThresh取0.2到0.3。若sdThresh太小，噪声会被当成独立IMF分离出来。

2.2 二维化改造：极值点提取、径向基插值与BIMF

把二维图像当作一个“平面信号”，BEMD的筛选流程就变成：寻找图像中的二维局部极大值和极小值点；用这些点拟合上包络面和下包络面；计算包络平均面；从原图中减去；重复筛选。得到的每层结果叫二维固有模态函数（BIMF），最后剩下的是残差。

二维极值点提取的难点在于“局部”和平坦区域。对于像素矩阵，一个极值点要同时和上下左右以及两个对角线共8个邻域比较。如果它严格大于所有8邻域就是局部极大值，严格小于所有8邻域就是局部极小值。但图像里总有成片的平坦区，比如天空、墙面，此时会有大量像素同时“等于”邻域，无法判定极大极小。常见解决方法是先用imregionalmax和imregionalmin找到极大/极小值区域，再用bwmorph(BW,'thin',Inf)把每个区域收缩成一个像素，也可以求区域重心作为极值点坐标。这样得到的极值点数量可控，包络也更平滑。

包络面插值的选择，直接决定BEMD分解质量。MATLAB里常见做法是：

F = scatteredInterpolant(xs, ys, vals, 'natural', 'linear');

其中xs、ys是极值点的横纵坐标，vals是极值点的灰度值，插值方法用'natural'（自然邻域）在平滑性和速度之间平衡。如果你追求更光滑的包络面，可以自己写薄板样条（TPS）插值，但TPS需要求解一个N*N的线性方程组，极值点超过两三千就非常慢，且容易内存不足。所以我在实际项目里，经常先用“抽稀”策略把极值点数量压到2000以内，再做'natural'插值。

筛选停止条件和一维类似，但二维SD的计算要对整幅图像逐像素做：

SD = sum(sum((H_prev - H_current).^2)) / sum(sum(H_prev.^2))

此外还有Cauchy型条件：计算包络平均面的能量与原图能量的比值，小于某个百分比就停止。我个人偏好用能量比值，因为它更直观，也不受边界像素突变影响。下面是二维筛选的一次迭代示意，采用'natural'插值：

function H = bemd_sift2d_demo(img, minDist, sdThresh) % 二维筛选一次迭代的近似示意 % 输入：img为double灰度图像，minDist为极值点最小间距，sdThresh停止阈值 % 输出：一个BIMF分量H I = img; % 用imregionalmax/min找极值区域，再提取中心点 maxReg = imregionalmax(I); minReg = imregionalmin(I); % 简化：直接用find取坐标，实际要处理区域收缩 [my, mx] = find(maxReg); % 如果极值点太多，按minDist抽稀 if ~isempty(my) && length(my) > 1 keep = randperm(length(my), min(length(my), 1000)); % 示例抽稀 mx = mx(keep); my = my(keep); end valMax = I(sub2ind(size(I), my, mx)); Fmax = scatteredInterpolant(mx, my, valMax, 'natural', 'linear'); [xx, yy] = meshgrid(1:size(I,2), 1:size(I,1)); envUp = Fmax(xx, yy); % 同理计算下包络envLo，这里省略 M = (envUp + envLo) / 2; H = I - M; % 停止判断：这里简化为前后差异小于阈值 sd = sum(sum((H - I).^2)) / sum(sum(I.^2)); if sd < sdThresh return; end end

注意：这个演示代码里我用了randperm抽稀，这是非常不可取的，它会让分解结果不固定。实际应该用空间距离聚类或者网格化取局部最大值。我只是为了展示结构。后面避坑章节会专门说这个。真正的BEMD库不会这么写，但很多粗糙的BEMD.zip里还真有随机抽稀，这就是为什么同一幅图跑两次结果不一样。

这个过程的复杂度主要在于插值。极值点数量是N，'natural'插值构建三角剖分的时间接近O(N log N)，之后在每个网格点求插值还要遍历邻域。所以一幅512x512的图像，极值点如果不过滤，很容易突破5000个，一次筛选可能就要几十秒，整幅分解几十个筛选循环，跑下来10分钟很正常。

2.3 BEMD和人工智能的关系：当预处理特征提取器用，还是当网络层用？

BEMD近几年在人工智能相关项目里出现，多数是以“特征工程模块”身份。图像先做BEMD，把得到的高频BIMF和残差当作多通道输入，再接CNN或SVM。这样做的好处是，BEMD能自适应地抓住图像中不同尺度的振荡模式，这些模式比原始像素更有物理意义。比如在医学图像分类中，BIMF1的纹理细节、残差的亮度分布，分别对应不同诊断特征。

也有少数人尝试把筛选过程改造成可微算子，放进神经网络里。难点在于极值点搜索和插值梯度回传。极值搜索本质是max操作，不可微，插值倒是可以线性化。所以常见方案是：固定极值点集合（比如每层采样网格），只让包络插值参与梯度计算。这写起来非常复杂，而且BEMD本身迭代次数不固定，网络训练时容易不稳。因此我一般不推荐把BEMD嵌入到网络内部，而是离线算好特征存盘。

在手头没有GPU的机器上，BEMD+传统分类器是一个性价比很高的方案。你可以先把训练图像库里的所有图像分解，提取每层BIMF的统计量，存成特征矩阵，然后训练一个SVM。整个过程不需要深度网络，代码量小，还能作为“人工智能大作业”的亮点——把数值算法和AI结合，比单纯调一个CNN有说服力。

下表简单对比BEMD和常见分解方法：

从这个表能看出来，BEMD的最大特色就是完全自适应和局部性，这正好是很多AI任务需要的。但成本高，所以实际使用时要控制图像尺寸和极值点数量。

3. 在MATLAB里跑通BEMD.zip：最小复现代码与文件结构

3.1 拿到压缩包先做三件事：检查函数依赖、测试入口、看数据格式

从网上下载的BEMD.zip打开后，别急着运行，先看看里面有什么。常见做法是解压后用MATLAB的dir命令遍历所有.m文件：

cd('D:\workspace\BEMD\'); fileList = dir('**/*.m'); % 递归找所有m文件 for i = 1:length(fileList) disp(fileList(i).name); end

这段代码会列出压缩包里的所有MATLAB脚本和函数。你要关注三件事：有没有demo或test脚本；有没有主函数，比如bemd2d.m或bemd_main.m；有没有明显缺失的依赖。很多压缩包是把论文源码的一部分上传了，少一个私有函数，比如find_extrema.m或interp_surface.m，直接跑主函数会报“Undefined function”。这时你先在整个文件列表里搜一下，如果确实没有，那就要去报错提示的调用链里逐个补。

第二步是测试入口。如果能看到demo.m或example.m，直接在MATLAB里运行它。但运行之前，建议先看开头几行，确认它读的数据文件是否存在。有时候demo里写的图片路径是./data/lena.bmp，但压缩包里根本没有data文件夹。我会写一个快速检查文件存在的代码：

% 检查关键文件是否存在 if exist('demo.m', 'file') == 2 disp('找到demo入口'); end if exist('lena.bmp', 'file') == 2 img = imread('lena.bmp'); disp(size(img)); else warning('缺少演示图像，请先准备一张灰度图'); end

第三步是看数据格式。BEMD一般要求输入是double类型、灰度图。很多代码会默认图像范围是[0,1]，如果直接传uint8的[0,255]，插值结果会是两三百的样子，后面计算能量就爆炸。我习惯统一做一次归一化：

if nargin == 0 img = imread('cameraman.tif'); end if isinteger(img) img = double(img) / double(intmax(class(img))); end

3.2 最小运行代码：读图、分解、显示BIMF

假设我们手头的BEMD.zip里有个函数叫bemd2d，它的接口大致是：

function [bimfs, residue] = bemd2d(img, layers, sd, minDist)

参数说明：img是MxN的double灰度图像；layers是分解层数；sd是停止阈值；minDist是极值点最小间距。这个接口不一定是你下载到的那个，但比较通用。最小运行代码如下：

% 读取并归一化图像 img = imread('cameraman.tif'); img = double(img) / 255; % 设置参数 layers = 3; sd = 0.3; minDist = 3; % 调用分解 [bimfs, residue] = bemd2d(img, layers, sd, minDist); % 可视化第一层和残差 figure('Name', 'BEMD分解结果'); subplot(1,3,1); imshow(bimfs{1}); title('BIMF1'); subplot(1,3,2); imshow(bimfs{2}); title('BIMF2'); subplot(1,3,3); imshow(residue); title('残差');

如果bimfs是三维数组，比如bimfs(:,:,1)，那么你必须按实际返回值方式调整。最好在调用前用which bemd2d和open bemd2d看一下函数定义。

这段代码的逻辑说明：归一化是必须的，因为BEMD的插值过程依赖数值尺度，如果输入是[0,255]的uint8，包络面的拟合矩阵会有较大的条件数，可能在筛选过程中产生NaN。而分解层数设为3是一般经验值：第一层是高频细节，第二层是中频结构，第三层基本就是低频轮廓，残差就是全局亮度趋势。如果你只做去噪，往往只要BIMF1就够了；做融合可能需要前两层；做分类可能前三层加残差。

3.3 把BIMF转成特征向量：供后续聚类/分类使用

对于人工智能项目，BEMD分解后的BIMF不是直接拍平送给分类器的，因为像素维度太高。常见做法是对每层计算一组统计量，得到固定长度的特征向量。下面这段代码从bimfs元胞数组中提取特征：

% 输入bimfs为1*K元胞数组，每层是M*N矩阵 numLayers = length(bimfs); feat = zeros(1, numLayers * 5); for i = 1:numLayers b = bimfs{i}; b = mat2gray(b); % 缩放到0-1，便于计算熵和互不相关 feat(i*5-4) = mean2(b); feat(i*5-3) = std2(b); feat(i*5-2) = sum(b(:).^2) / numel(b); % 平均能量 feat(i*5-1) = entropy(b); % 灰度共生矩阵对比度（需要影像处理工具箱） glcm = graycomatrix(b, 'Offset', [0 1; -1 1; -1 0; -1 -1]); stats = graycoprops(glcm, {'Contrast', 'Correlation'}); feat(i*5) = mean(stats.Contrast); end

这段代码里有几个点要注意：mat2gray是对数据进行线性归一化，让熵函数能够正确统计灰度分布；graycomatrix需要Image Processing Toolbox；graycoprops返回的是结构体数组，所以用mean取平均。如果没装工具箱，可以把熵和对比度换成简单的能量梯度：

gx = diff(b, 1, 2); gy = diff(b, 1, 1); feat(i*5-1) = mean(gx(:).^2) + mean(gy(:).^2);

这样整个特征提取就不依赖任何工具箱了。然后按行堆叠所有样本的特征，生成训练矩阵：

% 假设X_train是N*D矩阵，N是样本数，D是特征维度 X_train = [X_train; feat];

逻辑说明：BIMF特征向量的长度等于分解层数乘以每层特征数，如果分解层数波动（比如不同图像筛出的层数不一样），特征矩阵对不齐。解决办法是固定分解层数，让BEMD总是输出前K层加残差，不足时补零，多余时丢弃。很多MATLAB图像处理大作业里，大家就是用这种方式把非定长输出变成定长特征。

4. BEMD的三个必调参数：筛选次数、停止条件、插值点数

4.1 筛选次数（sifting iterations）直接决定BIMF的“纯度”

筛选次数是BEMD中最容易被忽视的参数。很多代码里它叫maxIter或numSift，默认值从5到20不等。筛选次数越多，每个BIMF越接近“局部对称平衡”的理想状态，但同时也意味着计算量线性增加。更重要的是，筛选次数过多会发生“过筛选”，把噪声从残差中吸入BIMF1。实测经验：对一般图像，10次够用，20次是极限；如果超过30次，BIMF1会全是椒盐状噪声。

判断筛选次数是否合适的办法是看相邻两次筛选结果的差异。你可以在BEMD循环里打点，或者自己写一个小脚本：

% 假设bemd2d内部允许回调iter，或者你自己改写筛选循环 for iter = [5 10 15 20 25] [C, ~] = bemd2d(img, 3, 0.3, 3, iter); % 假设第五个参数是maxIter energy = sum(C{1}(:).^2) / sum(img(:).^2); fprintf('iter=%d, BIMF1能量占比=%.4f\n', iter, energy); end

如果energy在iter=15之后突然升高，说明过筛开始。取拐点前的值比较安全。注意这个前提是BEMD函数真的接受maxIter参数，有些实现写死了，你需要改代码。我通常保留这个接口，因为这是超参搜索的第一步。

4.2 停止条件（SD阈值）调不好就是欠筛或过筛

SD阈值表示两次筛选之间的归一化差。阈值越小，筛选越苛刻，所需迭代次数越多。二维下，由于像素数量大，SD的数值容易偏小，如果沿用一维的0.3，往往还没筛完就停了，导致BIMF1里混着大尺度结构。我一般对图像设阈值0.4~0.6。不过这也和图像对比度有关：如果图像已经归一化到[0,1]，阈值0.5比较稳。

你可以用下面的脚本快速验证不同SD对残差的影响：

sds = [0.2 0.3 0.4 0.5 0.6]; for i = 1:length(sds) [~, res] = bemd2d(img, 3, sds(i), 3); resEnergy = sum(res(:).^2) / sum(img(:).^2); fprintf('SD=%.1f, 残差能量占比=%.3f\n', sds(i), resEnergy); end

一般来说，残差能量占比应该占原图的50%~80%左右。如果你看到SD=0.5时残差能量高达95%，说明分解深度远远不够——你可能需要增加分解层数，而不是继续调小SD。如果SD=0.3时残差能量只有20%，说明过分解了，许多低频信息都被抽到BIMF里去了。

4.3 插值点数与边界延拓：暗藏玄机的计算量黑洞

插值点数指的是参与构造包络面的极值点个数。这个参数通常不直接暴露，而是通过minDist或“最大极值点数”间接控制。我见过一个BEMD实现，内部用find(imregionalmax(img))直接找极值，一张512x512的图能找到6000多个极大值点。对这些点做scatteredInterpolant的'natural'插值，一次包络面就要构造并解一个三角剖分，内存开销几GB。如果你发现MATLAB内存飙升，八成是极值点太多。

设置minDist可以强制极值点之间至少相隔几个像素。在MATLAB里，一种实现方式是对二值极大值图做形态学腐蚀：

% 限制极值点最小间距 maxReg = imregionalmax(img); maxReg = imdilate(maxReg, strel('disk', minDist)) & maxReg; % 再用bwmorph提取中心 maxReg = bwmorph(maxReg, 'shrink', inf); [ys, xs] = find(maxReg);

这样可以把相邻极值点合并成一个。注意，imregionalmax会返回平台区所有像素为1，bwmorph shrink会把每一块连续区域收缩为一个像素，这是稳妥的做法。如果你想让极值点数量进一步减少，可以对xs, ys做一个网格统计，每个网格只留最大灰度极值点。

边界延拓则是另一个决定结果质量的参数。BEMD在图像边缘处没有相邻像素，包络面会在边缘处翘起，导致分解出的BIMF在边缘有伪振荡。常见解决方法是：在分解前对图像做对称延拓，扩展宽度为8~16像素，分解后裁掉。下面是一个延拓与裁剪的封装：

function [bimfs, res] = bemd2d_with_extend(img, layers, sd, minDist, ext) % 对图像四周做镜像延拓，然后调用bemd2d，最后裁剪 extImg = padarray(img, [ext ext], 'symmetric'); [bimfs, res] = bemd2d(extImg, layers, sd, minDist); for k = 1:length(bimfs) bimfs{k} = bimfs{k}(ext+1:end-ext, ext+1:end-ext); end res = res(ext+1:end-ext, ext+1:end-ext); end

这段代码里padarray的'symmetric'是镜像延拓，比'replicate'（重复边缘）更自然，能有效抑制边缘振铃。延拓宽度ext建议设为图像最长边的2%~4%，比如512图取16。

需要特别提醒：有些BEMD源码内部使用griddata('v4')，但MATLAB新版已经移除'v4'方法，导致直接报错。遇到这种情况，把griddata(x,y,v,xi,yi,'v4')改成scatteredInterpolant(x,y,v,'natural')后再插值，输出会有微小不同，但至少能跑。

5. BEMD实战避坑：MATLAB里最常见的5个翻车现场

5.1 现象：分解出的BIMF全是噪声，没有纹理层

你可能看到BIMF1是一幅几乎全是细碎斑点的图，看不出任何边缘或纹理。原因通常是筛选次数过多或者插值方法不当。特别是当代码用了线性插值'linear'时，包络面不光滑，局部均值面到处是折痕，筛出来的残余就是高频伪影。解决方法是改用'natural'插值，同时把最大筛选次数降到10左右。如果图像本身已经有噪声，先做一次轻度高斯滤波再分解，比如img = imgaussfilt(img,0.8)，否则BEMD会把噪声当作最高频振荡模式独立筛出来。

5.2 现象：运行时间爆炸，一幅512x512图像跑十分钟

这是BEMD最劝退的问题。原因除了极值点太多外，还可能是代码里用了双层循环遍历每个像素去判断极值，而不是用imregionalmax向量化。解决思路是先剖开看看哪个函数耗时。在MATLAB中运行：

profile on bemd2d(img, 3, 0.3, 3); profile viewer

如果scatteredInterpolant占了大头，说明极值点数量失控。先把minDist设到4或5，让极值点数量减少到原来的四分之一，运行时间能大幅缩短。如果瓶颈是循环，那只能改写向量化。另外，千万不要在筛选循环内部用imshow或disp打印中间过程，这会拖慢几十倍。

5.3 现象：边界处出现强烈的“振铃”或伪分量

边界振铃是插值包络在图像边缘没有足够极值点约束造成的。常见出现在BIMF1的上下左右边缘，有像是波浪的亮暗条纹。解决方法是先延拓再分解，延拓后裁掉。如果延拓后还有振铃，可以尝试对图像边缘做加权淡入淡出，或者将图像分块分解后做重叠融合。但重叠融合比较复杂，我的经验是先加镜像延拓，90%的振铃都能被压下去。

5.4 现象：不同MATLAB版本结果不一致

同一个BEMD.zip，在MATLAB R2016a和R2024a上跑出来的BIMF不一样，这是正常现象。原因在于griddata和scatteredInterpolant底层算法在不同版本中改动过，三角剖分和插值结果会有细微差异。另一个原因是旧版randperm行为可能不同，如果代码内部有随机抽样，结果天然不稳定。解决方法是固定插值方法为'natural'，固定随机种子，并且把关键参数写进注释。如果要求可复现，最好在同一版本下运行，并在结果保存时记录MATLAB版本和参数。

5.5 现象：内存不足（Out of Memory）在二维筛选中高发

当极值点数量达到几千甚至一万，scatteredInterpolant创建时可能直接抛出“Requested array exceeds maximum array size preference”。原因是'natural'插值需要构建Delaunay三角剖分，复杂度O(N log N)，内存需求随N快速上涨。解决方法：一是减少极值点数量，上面说过；二是用分块分解；三是如果图像太大，先imresize到256x256再分解。注意，BEMD本身对分辨率很敏感，强行缩放也会改变分解层次，但在内存受限的机器上这是必要的妥协。

6. 让BEMD真正工作起来：图像融合场景的完整验证与进阶技巧

最后这一章，我拿一个具体场景把整篇文章串起来：多聚焦图像融合。假设你有两张同一场景的照片，一张左边清晰、右边模糊，另一张右边清晰、左边模糊。经典做法是用拉普拉斯金字塔，但我们这里用BEMD来验证分解结果有没有意义。

融合规则很简单：两幅图像分别做BEMD，得到各自的前两层BIMF和残差。高频部分（BIMF1和BIMF2）用局部能量或梯度作为清晰度度量，取清晰度大的像素作为融合结果；残差部分直接平均，因为没有“清晰不清晰”之分，残差代表背景强度。实现代码如下：

% imgA、imgB是已配准的双灰度图，尺寸相同 [C_A, r_A] = bemd2d(imgA, 2, 0.4, 4); [C_B, r_B] = bemd2d(imgB, 2, 0.4, 4); % 初始融合结果用残差平均 fused = (r_A + r_B) / 2; % 逐层融合高频BIMF for k = 1:2 % 用水平Sobel梯度幅值代表清晰度 ga = abs(imfilter(C_A{k}, fspecial('sobel'))); gb = abs(imfilter(C_B{k}, fspecial('sobel'))); w = ga > gb; fused = fused + w .* C_A{k} + (1-w) .* C_B{k}; end

这段代码里我简化了清晰度衡量，只用了水平Sobel。实际可以用imgradient算多方向梯度，效果更稳：

ga = imgradient(C_A{k}); gb = imgradient(C_B{k});

验证方式也很重要，光看图片觉得“好像更清晰”不够。可以计算融合图像与原两张清晰部分的梯度能量对比，或者用量化指标：信息熵、标准差、平均梯度。MATLAB一行就能算：

metric = std2(fused) + entropy(mat2gray(fused)) + mean2(abs(imfilter(fused, fspecial('sobel'))));

如果metric明显高于两张原图各自的对应值，说明融合确实保留了更多细节。注意，BEMD对两幅图分解时用的是同一组参数，否则BIMF层数可能不同，导致C_A{1}和C_B{1}尺度不匹配。如果出现这个问题，强制设置layers=2，而且必须保证两幅图的尺寸完全一样。

最后分享一个我的个人习惯：每次跑BEMD之前，我都会把输入图像缩放到最大边不超过512像素。不是算法不行，而是工程上值得。BEMD的筛选过程是确定性的，但计算代价随图像面积近似平方增长，缩小四倍分辨率，时间能快十倍以上。先在小图上把参数调好，再放大到最终尺寸，能少熬夜。如果你一开始就在1024图上慢慢等，十有八九会等到怀疑人生。希望这个血泪经验能帮到你。

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