十年匠心定制 · 商业建站与技术教学双线并行 咨询热线:400-886-1026 service@lmnt.cn
ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

BEMD二维EMD分解:原理、工程选型与图像去噪实践

BEMD二维EMD分解:原理、工程选型与图像去噪实践 简介这是一份用于二维经验模态分解BEMD的MATLAB实现资源面向需要分析非线性、非平稳二维数据的科研人员与工程师。压缩包内含主程序bemdtry.m、核心分解函数TwoD_EMD.m以及sift_bicubic.m、findextm.m等辅助脚本覆盖图像读入、极值点提取、包络构造到IMF分解的完整流程同时提供lena.jpg标准测试图便于快速验证算法效果。该工具将经典EMD方法扩展到二维空间可将图像分解为多个内在模式函数IMF与残余分量有助于提取局部特征、去除噪声或揭示多尺度结构。资源共13个文件以MATLAB源码为主附带示例图像整体仅53KB轻量易用。已有533人学习下载适合图像处理、医学成像、气候遥感等领域的二维信号分析需求。通过阅读源码与运行示例可深入理解BEMD的sifting迭代机制并基于现有代码扩展至自定义数据集。1. BEMD 把二维 EMD 分解从理念变成了能跑的代码BEMDBidimensional Empirical Mode Decomposition解决的是图像没法直接用一维 EMD 分解的问题一张灰度图里的噪声、纹理和光照基底叠在同一层逐行套用一维 EMD 会把模态拆得行与行错位边缘还生出一串伪条纹。二维 EMD 分解把极值点检测和包络插值整个搬到平面上筛出的每一层 BIMF 都覆盖全图频率由高到低逐层剥离。bemds.rar 里装的就是一套常见的 MATLAB 实现标题里的可完成就是指它解压后能直接吃进灰度矩阵、吐出模态序列。下面讲讲筛选原理、包络插值选型以及真正把它们跑顺的三个关键位置——这些点定了BEMD 才从演示脚本变成你手上去噪和纹理分析的常规工具。2. BEMD 的筛选骨架与三个必须做对的工程选型2.1 一次筛选迭代里到底发生了什么这里的 BEMD 指的是 Bidimensional EMD处理的是灰度图像矩阵不是处理复数信号的 bivariate EMD别被同名缩写带偏。筛选sifting是 EMD 的核心二维版本只是让所有操作都作用在 I(x,y) 上。每一步做四件事在图像上检测局部极大值和局部极小值用这些散点插值出上包络 EU 和下包络 EL算出包络均值 M(EUEL)/2让当前信号减掉均值得到细节 hI−M。重复直到 h 满足本征模态条件拿到一个 BIMF再对残差做同样的流程。一维 EMD 用极值点与过零点数量至多差一判断 IMF二维图像没有自然的过零点定义所以工程实现基本都用相邻两次筛选的能量归一化差值 SD 来停。SD 的公式和一段可独立运行的筛选循环如下function h sift_2d(h0, sd_tol, max_iter) h h0; new_h h0; [H, W] size(h0); [X, Y] meshgrid(1:W, 1:H); for iter 1:max_iter [max_y, max_x, max_v] local_extrema(h, max); [min_y, min_x, min_v] local_extrema(h, min); if numel(max_v) 4 || numel(min_v) 4 break; % 极值点太少包络插值会退化 end up griddata(max_x, max_y, max_v, X, Y, cubic); low griddata(min_x, min_y, min_v, X, Y, cubic); up(isnan(up)) h(isnan(up)); low(isnan(low)) h(isnan(low)); mean_env (up low) / 2; new_h h - mean_env; sd sum((new_h(:) - h(:)).^2) / (sum(h(:).^2) eps); if sd sd_tol break; end h new_h; end h new_h; end这段代码只做一个模态的筛选。外层还需要对残差 rr−BIMF 循环调用 sift_2d直到残差极值点数量掉到个位数或达到预分层数。逻辑上的要点是griddata 的 cubic 插值在极值点凸包外返回 NaN这里用当前像素值直接回填避免空值把下一轮筛选冲乱SD 分母加 eps 是防止全零区域除零。参数上sd_tol 取 0.2~0.3 是 EMD 文献里的常见区间小于 0.2 会过度筛选把噪声细节硬压进低频模态max_iter 设 50~100超过这个数还收敛不了多半是某个像素在包络均值附近震荡直接取当前 new_h 退出比继续空转更好。2.2 极值点检测邻域窗口在控制分解分辨率local_extrema 的实现通常依赖图像处理工具箱的 imregionalmax / imregionalmin默认用 8 邻域判定function [yy, xx, vv] local_extrema(A, type) if strcmp(type, max) bw imregionalmax(A); % 8 邻域局部极大 else bw imregionalmin(A); % 8 邻域局部极小 end [yy, xx] find(bw); vv A(bw); end3×3 邻域会让每一点像素级噪声都参与建包络结果第一个 BIMF 几乎就是高频噪声本身。想让它承载真实纹理就把邻域放大到 5×5 或 7×7或者先做一次中值滤波再找极值点。邻域越大极值点越稀疏包络越光滑插值越快代价是小于邻域尺寸的细节直接漏进残差。对 512×512 的图我一般从 3×3 起步如果第一层和噪声高度相关而纹理留在第二层就把邻域改成 5×5 重新分解。包络、停止条件、邻域三个选型是互相咬合的建议一开始就按下面的基准值试别同时动三个参数参数常见取值调大的后果调小的后果极值点邻域3×3 ~ 7×7包络光滑、细节混入残差噪声主导 BIMF1包络插值cubic / v4 / natural更慢、包络更连续线性插值快但包络有折痕SD 阈值0.2 ~ 0.3提前停止、模态混叠过度筛选、模态振荡最大迭代数50 ~ 100耗时翻倍、可能过筛模态不稳定2.3 包络插值BEMD 慢和边角空值的根源二维包络插值是把稀疏极值点拟合成整张曲面这是 BEMD 和一维 EMD 在性能上的分水岭。一维三次样条的开销是 O(N)二维每次插值都要重建 Delaunay 三角网griddata 的 cubic 和 v4 在 1024×1024 的图上单次调用就是秒级。分解 4 层、每层 20 次筛选、每次两次插值就是上百次全图插值这是 bemds 跑起来合理慢的根本原因不是死循环。比慢更隐蔽的是边角空值。cubic 和 natural 只在极值点凸包内有值凸包外尤其是四角返回 NaN。回填原像素是保住边界的应急手段但回填区没有经过包络均值修正边缘会留下一条很窄的高频带。对边缘敏感的任务先做对称镜像延展再分解、完事裁回I_pad padarray(I, [8 8], symmetric); [imf_pad, res_pad] bemds(I_pad, maxIMF, 4, SD, 0.25); imf imf_pad(9:end-8, 9:end-8, :); res res_pad(9:end-8, 9:end-8);延展量取 8~16 个像素要求大于插值方法能拿到的有效半径。镜像延展比补零稳补零会在边界造出一圈新极值点让边界附近的包络被假的低洼拽下去。3. 在 MATLAB 里确认 bemds 主函数并跑通第一次二维 EMD 分解3.1 解压后的入口确认bemds.rar 解压出来文件结构各家略有差异但一般都能看到主函数bemds.m 或 bemd.m、插值子函数、示例图像和 demo 脚本。先别双击 demo 就等结果先挂路径、确认哪个文件是入口addpath(/your/path/bemds); which bemd which bemds help bemdwhich 返回空白说明包内没有这个名字就依次试 bemd、bemds、bemd2d。help 的输出重点看两处输入是二维 double 矩阵还是 cell输出是 M×N×K 三维数组还是 cell 数组。这决定了后面怎么按层取值。多数实现收 M×N 的 double 灰度图uint8 读进来要先转 double否则包络减法在整数域溢出分解结果完全不可信。imregionalmax 依赖图像处理工具箱没有工具箱时要用 imdilate 比较或手写 3×3 扫描替代逻辑相同但慢一些。3.2 最小调用代码与重建检查% demo_bemd.m —— 第一次跑通的最小流程 I im2double(imread(grains.png)); if size(I, 3) 1 I rgb2gray(I); end [IMF, RES] bemds(I, maxIMF, 4, SD, 0.25); recon sum(IMF, 3) RES; fprintf(最大重建误差: %g\n, max(abs(recon(:) - I(:)))); figure; subplot(2, 3, 1); imshow(I); title(原始); for k 1:4 subplot(2, 3, k 1); imshow(IMF(:, :, k), []); title(sprintf(BIMF %d, k)); end函数签名里的参数名以你手上包的 help 输出为准我写的 maxIMF、SD 是这类工具包最常见的叫法。重建误差是分解正确性的第一道判据误差超过 1e-10 级别优先怀疑输出格式。有的实现把最后一层和残差合在一起返回这时 sum(IMF,3) 会把残差也算进去重建误差恒等于零但对不上层数。imshow 第二参用 [] 自动拉伸对比度低对比度的模态不会显示成一片灰。3.3 三个必调参数和它们对应的故障信号参数作用故障信号maxIMF 分解层数控制分出几个 BIMF残差里还有明显纹理说明层数不够SD 阈值筛选停止的灵敏度模态出现周期性振荡说明阈值偏小极值点邻域决定极值点密度第一层全是椒盐噪声说明邻域偏小分解层数从 3 到 5 起步判断是否加层看残差残差还有肉眼可见的亮暗起伏就加残差已是平滑背景渐变就停再加就是把基底拆碎。第一次运行把 maxIMF 压到 2先确认单层耗时再放开层数。提示BEMD 耗时随图像尺寸超线性增长1024×1024 的 4 层分解在普通台式机上经常要几分钟到十几分钟。先用 2 层验证流程再跑全量能省掉大量把参数调错后白等的时长。4. BEMD 的边界效应、插值退化与三个提速手段4.1 边角空值与回填的两难第 2 章的回填方案能保住边界不出现 NaN但回填像素没有参与包络修正所以首尾几行几列的第一模态总是偏亮或偏暗。判断边界是不是坏了直接把 IMF(:,:,1) 的四周像素和中心区域对比如果边缘行有规律性的宽条纹说明回填区把假极值带进来了。另一个常见表现是第一模态的四个角出现明显的X形亮度。这时候先差分化看边界行std(IMF(1:8,:,1), 0, 2)如果远大于内部行的方差就该上镜像延展而不是继续调 SD。4.2 耗时到底耗在哪里用 profile on 跑一次 2 层分解耗时几乎全部落在 griddata 上。插值之外imregionalmax 的计算量只是线性扫描占比不到 5%。所以优化方向只有一个减少插值调用次数和降低单次插值成本。减少次数靠调大 SD 阈值、调大极值点邻域降低成本靠换插值方法。常见做法是探测阶段用 linear速度比 cubic 快一个量级参数定案后用 cubic 或 v4 出最终结果。注意 v4 是全图径向基插值内存占用随极值点数平方上涨极值点过万时很容易把内存吃满遇到卡死先看是不是它。提速手段做法代价线性插值探测参数探测用 linear定案换 cubic探测阶段的包络不光滑裁块探测256×256 定参全图只跑一次边缘处的参数可能有偏差极值点抽稀最小距离 4px 合并极值点细纹理可能被抹掉4.3 三个能落地的提速手段第一是裁块探测。在中间裁一个 256×256 的区域跑完整流程把 SD 和层数定下来再全图只跑一次。全图跑之前把结果预分配到 .mat防止跑到一半内存不足前功尽弃probe I(129:384, 129:384); [IMF_probe, RES_probe] bemds(probe, maxIMF, 4, SD, 0.3); % ...... 分析探测结果后定案 ...... [IMF, RES] bemds(I, maxIMF, num_layers, SD, sd_final); save(bimf_result.mat, IMF, RES, -v7.3); % 超过 2GB 必须 v7.3第二是极值点抽稀。imregionalmax 之后做一次最小距离过滤距离小于 4 个像素的极值点只保留强度最大的那个极值点数能压掉一半以上插值时间跟着下降。抽稀过度会让细纹理直接消失抽稀窗口不要超过你设定的最小细节尺度。第三是换 scatteredInterpolant。对同一批散点建立三角网后可以反复查询配合 natural 插值比 griddata cubic 快且同样会外插 NaN 到凸包外回填逻辑不用改。这三招叠加1024×1024 的 4 层分解通常能从十几分钟压到四分钟以内。5. 用 BEMD 做图像去噪的验证步骤与判据5.1 丢弃第一层的方式与适用范围BIMF1 集中了最高频能量对平稳噪声占主导的图像直接丢掉第一层再重建是 BEMD 最稳的去噪路径denoised sum(IMF(:, :, 2:end), 3) RES; [gx, gy] gradient(denoised); flat_mag mean(sqrt(gx(:).^2 gy(:).^2)); c12_all corrcoef(IMF(:, :, 1), IMF(:, :, 2)); c12 c12_all(1, 2);逻辑说明denoised 用除第一层外的所有模态加残差重建等价于原图减掉 BIMF1。flat_mag 是重建图梯度模的均值平坦区域占比大的图这个值应该比原图的对应值明显下降c12 是前两层模态的皮尔逊相关系数按经验阈值取 0.1 以下高于这个值说明两层之间还有模态混叠需要把 SD 调大或邻域调小重跑。有干净参考图时直接比 PSNR/SSIM这两个函数在图像处理工具箱里没有就用 mse 转 psnr 自己算。不要条件反射地丢第一层。纹理分析、边缘检测场景里第一层常常是细节的主要载体丢掉后重建图会出现磨皮效果。动手丢之前先看频谱对第一个模态做 fft2把能量环带画出来噪声主导时能量集中在高频圆环纹理主导时能量会沿特定方向延伸二者在频谱图上非常容易区分。5.2 参数收敛的验证组合只凭一张图的视觉效果调参不靠谱建议固定一个三件套判据重建误差在 1e-10 量级验证分解完整IMF1 与 IMF2 相关系数低于 0.1验证模态分离到位残差梯度均值低于任何单个模态的梯度均值验证基底平坦。三个条件同时满足时参数组合基本收敛否则按故障信号表格反向定位是层数、阈值还是邻域的问题。5.3 边界处理是否到位的最后检查最后一招是边界验证。把镜像延展法用在同一张图上跑两遍一遍不延展、一遍延展 16 像素裁回原尺寸后对比两者 IMF1 的边缘两行。如果边缘像素差超过该模态动态范围的 5%说明边界还在污染分解延展量要再加如果差异集中在角点而不是整条边多半是插值方法在凸包顶点附近退化换成 natural 再试。看到 IMF 间相关系数低于 0.1、残差梯度在平坦区域接近零、被丢弃的模态频谱集中在高频环带这组参数就可以固定下来转入你要做的分割、融合或纹理量化任务了。本文还有配套的精品资源点击获取
返回列表