尧图网页设计 
简介二维经验模态分解BEMD是一种面向非线性、非平稳二维数据的多尺度分析技术在图像去噪、增强与特征提取中应用广泛。该压缩包提供了完整的 BEMD 图像分解实现共二十二个文件包含十七张仿真结果图、四个 MATLAB 源文件与一份算法说明文档。其中源文件覆盖包络构造、均值面计算、IMF 迭代分离等核心步骤仿真图直观展示原始图像、各阶固有模态函数、包络及残差重构图。整体大小约 493KB轻便易用。代码基于 MATLAB 编写结构清晰变量命名规范便于研究者直接运行或二次调整。资源适合图像处理、信号分析方向的研究生、科研人员或工程师能帮助快速理解 BEMD 原理并复现分解流程对于医学影像、遥感图像或纹理特征分析场景也可利用源码中的参数调节机制优化分解效果。说明文档对迭代停止条件、包络插值等细节作了较清晰交代进一步降低了上手门槛。目前已有 140 人学习使用表明其在教学和科研中具备一定参考价值。1. 为什么二维经验模态分解 BEMD 值得自己跑一遍二维经验模态分解BEMD是把一维 EMD 的自适应分解思想搬到图像上的直接产物它不预设小波基或傅里叶基而是靠图像自身的局部极值逐层筛分出频率从高到低的二维本征模态函数BIMF最后剩下一张残差。对做图像处理的人来说BEMD 的核心价值是数据驱动——同一套代码能同时处理光学纹理、医学影像和遥感图分解结果天然是一组从细节到背景的图层去噪、融合、增强都能直接复用。这篇文章把 BEMD 中文资料里经常被跳过的实现细节补齐二维极值怎么找、包络面怎么插值、停止准则设多少、边界怎么处理并用 Python 给出一套能直接改的最小实现。2. BEMD 的核心机制包络、筛分与停止准则2.1 从一维 EMD 到二维极值点从序列变成散点一维 EMD 的筛分过程可以用一句话概括把信号的极大值点用三次样条连成上包络极小值点连成下包络取均值后从原信号中减去反复迭代直到剩余分量满足本征模态函数IMF的两个条件。对一幅 M×N 的图像来说极值点不再是按时间排序的序列而是平面上没有天然邻接关系的散点上下包络也从两条曲线变成两张曲面。这个从曲线到曲面的转变是 BEMD 实现里所有复杂度上升的根源下面所有环节都在围绕这一件事转。因为要找的是曲面BEMD 每一步筛分都要回答三个子问题如何定义并定位二维局部极值如何用稀疏的极值散点插值出完整的包络面如何判定一次筛分可以结束。这三个问题分别对应形态学极值检测、散点插值、停止准则也是调参时唯一需要关心的三个位置。2.2 筛分过程的三个关键环节2.2.1 用形态学滤波定位二维极值最稳妥的二维极值定义是邻域比较某个像素的灰度值等于它在 k×k 邻域内的最大值它就是局部极大点等于最小值则是局部极小点。用scipy.ndimage的maximum_filter和minimum_filter一次就能拿到两个掩膜比手写双层循环快两个数量级而且modereflect能顺带处理图像边界处的滤波填充。import numpy as np from scipy.ndimage import maximum_filter, minimum_filter def find_extrema(img, win3): 返回极值点坐标注意返回顺序是 x, y 数组 max_f maximum_filter(img, sizewin, modereflect) min_f minimum_filter(img, sizewin, modereflect) max_pts (img max_f) min_pts (img min_f) my, mx np.where(max_pts) ny, nx np.where(min_pts) return mx, my, nx, ny窗口win在这里决定了极值的最小尺度取 3 时单像素噪声就可能成为极值点取 5 或 7 时更平滑的局部峰值才会被选中。这个参数直接影响第一层 BIMF 吸收的是纯噪声还是细小纹理实际调参时一般从 3 起步观察第一层结果再决定是否加大。2.2.2 稀疏极值插值成包络面极值点在整个图像上通常只占几百到几千个而上、下包络必须返回完整的 M×N 曲面中间的空洞全靠插值填补。我一般用scipy.interpolate.griddatapoints传极大点的 (x, y) 坐标values传这些点的灰度值method选cubic。它的底层是 Delaunay 三角剖分加 Clough-Tocher 三次插值散点分布不均匀时也不会出现线性插值那种明显的棱线。另一种常见做法是径向基函数RBF薄板样条曲面更光滑但求解规模随极值点数近似三次方增长图像稍大就等不起。有一个必须处理的细节griddata对凸包外的查询点返回nan。筛分循环里任何一个nan都会在下一轮迭代中扩散到整张图所以插值结果必须带一个fill_value兜底通常取当前极值点灰度的均值。2.2.3 IMF 判定与二维停止准则一维 IMF 的经典定义要求极值点数量与过零点数量之差的绝对值不超过 1且上下包络均值处处为零。二维没有直接的过零点概念工程实现普遍简化成两个近似条件包络均值趋近于零、极值点数量随分解层数单调下降。前者用 SD 准则度量后者靠观察每一层 BIMF 的极值点数量来判断是否还能继续分解。SD 准则计算的是连续两次筛分结果之间的相对变化量SD Σ(h_prev - h_cur)² / Σ(h_prev)²一维 EMD 的经验阈值是 0.20.3到了二维包络插值本身有数值误差SD 想稳定收敛到 0.2 以下往往要多跑十几轮所以我把阈值放宽到 0.30.5同时把单次筛分的最大迭代次数限制在 2050两个条件哪个先满足都退出。这个取舍在后续参数表里还会展开。2.3 分解结果的结构BIMF 的物理含义分解完成后结果看起来像一组从细到粗的边缘图层。BIMF1 集中了灰度变化最剧烈的高频成分通常是噪声、细微纹理、细胞膜这类精细结构越往后层的 BIMF 携带的纹理尺度越大最后的残差约等于图像的低频背景或光照场。对做图像增强的人来说残差就是可以直接做直方图均衡的光照层对做去噪的人来说BIMF1 往往是默认要被丢弃的对象。需要注意BEMD 的分解结果没有唯一解同一幅图极值窗口、插值方法和停止准则不同分解出的 BIMF 频带就不同。这是自适应分解方法的固有属性不是程序 bug对比不同实现的效果时一定要先固定住这三个参数。3. 用 Python 从零实现一个最小可用的 BEMD3.1 包络插值函数与筛分循环下面这套代码是 BEMD 的最小骨架极值检测复用上一章的find_extremaenvelope负责把稀疏极值点插值成整张包络面sift_one_bimf负责完整走完一次筛分。把这三段拼起来就是一个可运行的分解器。from scipy.interpolate import griddata def envelope(img, px, py, pv): 由稀疏极值点插值出完整包络面。 px, py 是极值点坐标pv 是这些点上的灰度值 h, w img.shape if len(pv) 6: # 极值点太少插值会严重失真退化为常数面 return np.full_like(img, pv.mean()) yy, xx np.mgrid[0:h, 0:w] return griddata((px, py), pv, (xx, yy), methodcubic, fill_valuepv.mean()) def sift_one_bimf(img, win3, max_iter25, sd_goal0.3): 对输入图像做一次完整筛分返回一个 BIMF h img.copy() for _ in range(max_iter): mx, my, nx, ny find_extrema(h, win) up envelope(h, mx, my, h[my, mx]) low envelope(h, nx, ny, h[ny, nx]) mean_env (up low) / 2.0 cand h - mean_env # SD 用相邻两次筛分结果的相对变化衡量收敛程度 sd np.sum((cand - h) ** 2) / (np.sum(h ** 2) 1e-12) h cand if sd sd_goal: break return henvelope里先检查极值点数量少于 6 个直接返回常数面这是因为 Clough-Tocher 插值至少需要 4 个非共线的点才能构建三角剖分数量太少时插值结果会出现大范围的畸形起伏。sd计算时给分母加了一个1e-12的 epsilon防止图像全为零灰度时除零。每次迭代都重新调用find_extrema因为减去包络均值后极值点位置和数量都会变化不能复用上一轮的检测结果。3.2 完整分解主流程与停止条件有了单个 BIMF 的筛分函数主流程就是一个层层剥离的循环先对残差做筛分提取 BIMF再把这个 BIMF 从残差中减掉然后对新的残差做下一轮。什么时候停除了用户指定层数还有一个硬性条件残差里找不到足够的极值点。def bemd_decompose(img, n_imfs3, win3, max_iter25, sd_goal0.3): BEMD 主流程返回 imfs 列表和残差 imfs [] residue img.copy() for k in range(n_imfs): mx, my, nx, ny find_extrema(residue, win) if len(mx) 2 or len(nx) 2: print(f第 {k1} 层极值不足提前停止) break bimf sift_one_bimf(residue, win, max_iter, sd_goal) imfs.append(bimf) residue residue - bimf return imfs, residue分解是否继续的判定依据是极大点和极小点是否同时存在。每一轮筛分都会消耗掉一部分局部极值所以残差会变得越来越平滑极值点越来越少。当某一层连 2 个极大点都找不到时插值已经无法构造出有意义的包络面继续分解只会产出数值噪声。3.3 在合成图像上运行并验证重建scipy.misc.face在较新的 SciPy 版本里被移除了直接用numpy生成一张多尺度合成图更省事而且能确定性地知道哪一层该出现什么结构yy, xx np.mgrid[0:256, 0:256] test_img (128 60 * np.sin(xx / 8) * np.cos(yy / 12) # 大尺度纹理 30 * np.sin(xx / 2) # 中等尺度纹理 10 * np.random.randn(256, 256)) # 加性噪声 imfs, residue bemd_decompose(test_img, n_imfs3, win3, max_iter20, sd_goal0.3) recon sum(imfs, residue) print(重建误差:, np.sqrt(np.mean((recon - test_img) ** 2)))sum(imfs, residue)是 Python 内建求和的第二个参数用法效果等价于residue imfs[0] imfs[1] ...。如果重建误差在1e-10量级说明分解和重构过程是自洽的误差只来源于浮点舍入如果误差明显偏大大概率是筛分循环里h的更新顺序写错了。跑通这一步之后就可以用真实图像替代合成图观察每一层 BIMF 的纹理尺度是否按预期从细到粗排布。4. BEMD 参数怎么调、边界怎么处理、分解质量怎么判断4.1 四个必调参数及推荐取值BEMD 需要人工干预的参数就四个极值窗口win、单次筛分最大迭代数max_iter、SD 目标值sd_goal、分解层数n_imfs。其余如插值方法选定后一般不再改动。参数推荐值控制什么调大/调小的后果win37极值检测邻域决定最细纹理尺度调大BIMF1 吞掉更大结构调小单像素噪声成为极值第一层被噪声主导max_iter2050每次筛分的轮数上限太小BIMF 里残留包络均值模态不纯太大过筛且耗时成倍增加sd_goal0.30.5筛分收敛阈值越小越严格但二维插值误差导致实际收益有限建议不要低于 0.2n_imfs35分解层数越多残差越干净但后层极值稀疏包络插值失真加剧win和sd_goal是耦合的win加大后极值点变少包络面更平缓SD 收敛更快此时sd_goal可以适当收紧到 0.3反过来win3时极值点多每一轮筛分的改动量大sd_goal0.5更实际。我一般先用win3, sd_goal0.3, max_iter20跑一遍看 BIMF1 是否混入了明显的中频结构再决定是否把win升到 5。4.2 边界效应与反射填充griddata的cubic插值只能在 Delaunay 三角剖分的凸包范围内计算凸包外的点全部由fill_value兜底。这导致分解结果在图像四角和边缘出现一整片均值平坦区下一轮筛分在这个区域里找不到任何极值点边界失真会一层层向内侵蚀。处理这个问题的常见做法是分解前对图像做反射填充分解完再裁掉填充区域pad 16 padded np.pad(test_img, pad, modereflect) pimfs, presidue bemd_decompose(padded, n_imfs3, win3, max_iter20, sd_goal0.3) imfs [pimfs[k][pad:-pad, pad:-pad] for k in range(len(pimfs))] residue presidue[pad:-pad, pad:-pad]我倾向于用反射填充而不是常数或对称填充因为常数填充会在边界制造一条明显的灰度跳变直接变成 BIMF1 里的假纹理对称填充在图像边界本身不连续时会让边界两侧的极值点被认为相邻产生虚假的包络振荡。反射填充相当于让图像在边界处光滑延拓对极值分布的影响最小。填充宽度取 1020 像素即可超过极值窗口 3 倍就足够了。4.3 平坦区域与极值板结灰度完全相同的平坦区域会给极值检测带来一个隐蔽问题整块区域都等于邻域最大值find_extrema会把每一个像素都标记为极大点极值点数量瞬间爆炸插值结果也毫无意义。这种情况在医学影像和遥感图的均匀背景区非常常见处理手段是对极值掩膜做连通域收缩只保留每个平板区域的中心点from scipy.ndimage import label def thin_plateau(mask): 把连成片的极值掩膜收缩成每块区域的中心点 lab, n label(mask) centers [] for i in range(1, n 1): ys, xs np.where(lab i) centers.append((int(ys.mean()), int(xs.mean()))) arr np.array(centers) if arr.size 0: return np.array([], dtypeint), np.array([], dtypeint) return arr[:, 1], arr[:, 0] # 返回 x, ylabel把连通的极值区域赋予同一个编号取每块区域的质心作为代表极值点。这个函数替换掉find_extrema里np.where那一步之后极值板结问题就消失了。实际处理中还要注意质心的坐标是浮点均值转整型可能出现两个相邻区域质心跳到同一个像素的情况筛分循环里对极值点做一次去重更稳妥。4.4 用两个硬指标判断分解质量分解质量不能只看重建误差还要看层间泄漏和极值递减规律。两个指标我每跑一组参数都会算一是重建误差 RMSE验证分解自洽二是正交性指标 OI衡量各层 BIMF 之间的信息重叠程度。OI 的定义是不同层逐像素乘积的绝对值之和除以分解总和能量的平方工程上小于 0.05 说明分解质量不错超过 0.1 就需要检查参数。此外还有一个经验规律每层 BIMF 的极值点数应该逐层明显下降如果 BIMF2 的极值点数和 BIMF1 差不多说明win太小或sd_goal太松第一层没有把该吸收的高频成分吸收干净。计数方法直接复用find_extrema(imf, win)的返回长度即可。5. BEMD 图像去噪与纹理分离一个直接能用的场景5.1 噪声落在哪一层取决于窗口和纹理尺度加性高斯噪声的尺度是单像素级在win3时绝大部分噪声极值会在第一轮就被吸收进 BIMF1。如果图像本身带有比噪声更细的纹理两者会一起出现在 BIMF1 里此时直接丢弃整层会连纹理一起抹掉。区分方法是看 BIMF1 的均值纯噪声层的均值接近零纹理层的均值往往有明显的空间结构。处理这类情况时对 BIMF1 再做一次阈值收缩而不是整层丢弃效果通常更好。5.2 去噪重建与正交性验证代码def orthogonality_index(imfs, residue): OI 越小说明各层信息重叠越少 terms imfs [residue] recon sum(imfs, residue) total 0.0 for i in range(len(terms)): for j in range(i 1, len(terms)): total np.abs(np.sum(terms[i] * terms[j])) return total * 2.0 / (np.sum(recon ** 2) 1e-12) imfs, residue bemd_decompose(noisy_img, n_imfs3, win3, max_iter20, sd_goal0.3) denoised residue sum(imfs[1:]) # 丢弃 BIMF1保留其余层 oi orthogonality_index(imfs, residue) print(正交性指标 OI:, oi)丢弃 BIMF1 后直接用sum(imfs[1:])叠加剩余层和残差得到的就是去噪结果。OI 的值在 0.01 到 0.05 之间属于正常范围如果大于 0.1说明相邻 BIMF 之间出现了频带交叠主要原因是sd_goal偏大导致筛分提前退出模态没有充分分离把阈值从 0.5 收到 0.3 重跑即可。5.3 一套可以直接复制的参数组合实际项目中我用的默认组合是win3, n_imfs3, max_iter20, sd_goal0.3配合 16 像素的反射填充。对 512×512 的灰度图这个配置单次完整分解在普通笔记本上大约几十秒OI 通常在 0.03 以下。如果第一层纹理混入严重把win升到 5 并把sd_goal收到 0.25代价是分解时间增加约一倍。如果想要更强的去噪效果可以对 BIMF1 做软阈值收缩后再重建imfs[0] * (np.abs(imfs[0]) t)其中t取 BIMF1 标准差的 1.5 倍这比整层丢弃更柔和也保留了噪声与纹理混叠时的高频细节。本文还有配套的精品资源点击获取 BEMD二维经验模态分解的Python实现与参数调优指南
简介二维经验模态分解BEMD是一种面向非线性、非平稳二维数据的多尺度分析技术在图像去噪、增强与特征提取中应用广泛。该压缩包提供了完整的 BEMD 图像分解实现共二十二个文件包含十七张仿真结果图、四个 MATLAB 源文件与一份算法说明文档。其中源文件覆盖包络构造、均值面计算、IMF 迭代分离等核心步骤仿真图直观展示原始图像、各阶固有模态函数、包络及残差重构图。整体大小约 493KB轻便易用。代码基于 MATLAB 编写结构清晰变量命名规范便于研究者直接运行或二次调整。资源适合图像处理、信号分析方向的研究生、科研人员或工程师能帮助快速理解 BEMD 原理并复现分解流程对于医学影像、遥感图像或纹理特征分析场景也可利用源码中的参数调节机制优化分解效果。说明文档对迭代停止条件、包络插值等细节作了较清晰交代进一步降低了上手门槛。目前已有 140 人学习使用表明其在教学和科研中具备一定参考价值。1. 为什么二维经验模态分解 BEMD 值得自己跑一遍二维经验模态分解BEMD是把一维 EMD 的自适应分解思想搬到图像上的直接产物它不预设小波基或傅里叶基而是靠图像自身的局部极值逐层筛分出频率从高到低的二维本征模态函数BIMF最后剩下一张残差。对做图像处理的人来说BEMD 的核心价值是数据驱动——同一套代码能同时处理光学纹理、医学影像和遥感图分解结果天然是一组从细节到背景的图层去噪、融合、增强都能直接复用。这篇文章把 BEMD 中文资料里经常被跳过的实现细节补齐二维极值怎么找、包络面怎么插值、停止准则设多少、边界怎么处理并用 Python 给出一套能直接改的最小实现。2. BEMD 的核心机制包络、筛分与停止准则2.1 从一维 EMD 到二维极值点从序列变成散点一维 EMD 的筛分过程可以用一句话概括把信号的极大值点用三次样条连成上包络极小值点连成下包络取均值后从原信号中减去反复迭代直到剩余分量满足本征模态函数IMF的两个条件。对一幅 M×N 的图像来说极值点不再是按时间排序的序列而是平面上没有天然邻接关系的散点上下包络也从两条曲线变成两张曲面。这个从曲线到曲面的转变是 BEMD 实现里所有复杂度上升的根源下面所有环节都在围绕这一件事转。因为要找的是曲面BEMD 每一步筛分都要回答三个子问题如何定义并定位二维局部极值如何用稀疏的极值散点插值出完整的包络面如何判定一次筛分可以结束。这三个问题分别对应形态学极值检测、散点插值、停止准则也是调参时唯一需要关心的三个位置。2.2 筛分过程的三个关键环节2.2.1 用形态学滤波定位二维极值最稳妥的二维极值定义是邻域比较某个像素的灰度值等于它在 k×k 邻域内的最大值它就是局部极大点等于最小值则是局部极小点。用scipy.ndimage的maximum_filter和minimum_filter一次就能拿到两个掩膜比手写双层循环快两个数量级而且modereflect能顺带处理图像边界处的滤波填充。import numpy as np from scipy.ndimage import maximum_filter, minimum_filter def find_extrema(img, win3): 返回极值点坐标注意返回顺序是 x, y 数组 max_f maximum_filter(img, sizewin, modereflect) min_f minimum_filter(img, sizewin, modereflect) max_pts (img max_f) min_pts (img min_f) my, mx np.where(max_pts) ny, nx np.where(min_pts) return mx, my, nx, ny窗口win在这里决定了极值的最小尺度取 3 时单像素噪声就可能成为极值点取 5 或 7 时更平滑的局部峰值才会被选中。这个参数直接影响第一层 BIMF 吸收的是纯噪声还是细小纹理实际调参时一般从 3 起步观察第一层结果再决定是否加大。2.2.2 稀疏极值插值成包络面极值点在整个图像上通常只占几百到几千个而上、下包络必须返回完整的 M×N 曲面中间的空洞全靠插值填补。我一般用scipy.interpolate.griddatapoints传极大点的 (x, y) 坐标values传这些点的灰度值method选cubic。它的底层是 Delaunay 三角剖分加 Clough-Tocher 三次插值散点分布不均匀时也不会出现线性插值那种明显的棱线。另一种常见做法是径向基函数RBF薄板样条曲面更光滑但求解规模随极值点数近似三次方增长图像稍大就等不起。有一个必须处理的细节griddata对凸包外的查询点返回nan。筛分循环里任何一个nan都会在下一轮迭代中扩散到整张图所以插值结果必须带一个fill_value兜底通常取当前极值点灰度的均值。2.2.3 IMF 判定与二维停止准则一维 IMF 的经典定义要求极值点数量与过零点数量之差的绝对值不超过 1且上下包络均值处处为零。二维没有直接的过零点概念工程实现普遍简化成两个近似条件包络均值趋近于零、极值点数量随分解层数单调下降。前者用 SD 准则度量后者靠观察每一层 BIMF 的极值点数量来判断是否还能继续分解。SD 准则计算的是连续两次筛分结果之间的相对变化量SD Σ(h_prev - h_cur)² / Σ(h_prev)²一维 EMD 的经验阈值是 0.20.3到了二维包络插值本身有数值误差SD 想稳定收敛到 0.2 以下往往要多跑十几轮所以我把阈值放宽到 0.30.5同时把单次筛分的最大迭代次数限制在 2050两个条件哪个先满足都退出。这个取舍在后续参数表里还会展开。2.3 分解结果的结构BIMF 的物理含义分解完成后结果看起来像一组从细到粗的边缘图层。BIMF1 集中了灰度变化最剧烈的高频成分通常是噪声、细微纹理、细胞膜这类精细结构越往后层的 BIMF 携带的纹理尺度越大最后的残差约等于图像的低频背景或光照场。对做图像增强的人来说残差就是可以直接做直方图均衡的光照层对做去噪的人来说BIMF1 往往是默认要被丢弃的对象。需要注意BEMD 的分解结果没有唯一解同一幅图极值窗口、插值方法和停止准则不同分解出的 BIMF 频带就不同。这是自适应分解方法的固有属性不是程序 bug对比不同实现的效果时一定要先固定住这三个参数。3. 用 Python 从零实现一个最小可用的 BEMD3.1 包络插值函数与筛分循环下面这套代码是 BEMD 的最小骨架极值检测复用上一章的find_extremaenvelope负责把稀疏极值点插值成整张包络面sift_one_bimf负责完整走完一次筛分。把这三段拼起来就是一个可运行的分解器。from scipy.interpolate import griddata def envelope(img, px, py, pv): 由稀疏极值点插值出完整包络面。 px, py 是极值点坐标pv 是这些点上的灰度值 h, w img.shape if len(pv) 6: # 极值点太少插值会严重失真退化为常数面 return np.full_like(img, pv.mean()) yy, xx np.mgrid[0:h, 0:w] return griddata((px, py), pv, (xx, yy), methodcubic, fill_valuepv.mean()) def sift_one_bimf(img, win3, max_iter25, sd_goal0.3): 对输入图像做一次完整筛分返回一个 BIMF h img.copy() for _ in range(max_iter): mx, my, nx, ny find_extrema(h, win) up envelope(h, mx, my, h[my, mx]) low envelope(h, nx, ny, h[ny, nx]) mean_env (up low) / 2.0 cand h - mean_env # SD 用相邻两次筛分结果的相对变化衡量收敛程度 sd np.sum((cand - h) ** 2) / (np.sum(h ** 2) 1e-12) h cand if sd sd_goal: break return henvelope里先检查极值点数量少于 6 个直接返回常数面这是因为 Clough-Tocher 插值至少需要 4 个非共线的点才能构建三角剖分数量太少时插值结果会出现大范围的畸形起伏。sd计算时给分母加了一个1e-12的 epsilon防止图像全为零灰度时除零。每次迭代都重新调用find_extrema因为减去包络均值后极值点位置和数量都会变化不能复用上一轮的检测结果。3.2 完整分解主流程与停止条件有了单个 BIMF 的筛分函数主流程就是一个层层剥离的循环先对残差做筛分提取 BIMF再把这个 BIMF 从残差中减掉然后对新的残差做下一轮。什么时候停除了用户指定层数还有一个硬性条件残差里找不到足够的极值点。def bemd_decompose(img, n_imfs3, win3, max_iter25, sd_goal0.3): BEMD 主流程返回 imfs 列表和残差 imfs [] residue img.copy() for k in range(n_imfs): mx, my, nx, ny find_extrema(residue, win) if len(mx) 2 or len(nx) 2: print(f第 {k1} 层极值不足提前停止) break bimf sift_one_bimf(residue, win, max_iter, sd_goal) imfs.append(bimf) residue residue - bimf return imfs, residue分解是否继续的判定依据是极大点和极小点是否同时存在。每一轮筛分都会消耗掉一部分局部极值所以残差会变得越来越平滑极值点越来越少。当某一层连 2 个极大点都找不到时插值已经无法构造出有意义的包络面继续分解只会产出数值噪声。3.3 在合成图像上运行并验证重建scipy.misc.face在较新的 SciPy 版本里被移除了直接用numpy生成一张多尺度合成图更省事而且能确定性地知道哪一层该出现什么结构yy, xx np.mgrid[0:256, 0:256] test_img (128 60 * np.sin(xx / 8) * np.cos(yy / 12) # 大尺度纹理 30 * np.sin(xx / 2) # 中等尺度纹理 10 * np.random.randn(256, 256)) # 加性噪声 imfs, residue bemd_decompose(test_img, n_imfs3, win3, max_iter20, sd_goal0.3) recon sum(imfs, residue) print(重建误差:, np.sqrt(np.mean((recon - test_img) ** 2)))sum(imfs, residue)是 Python 内建求和的第二个参数用法效果等价于residue imfs[0] imfs[1] ...。如果重建误差在1e-10量级说明分解和重构过程是自洽的误差只来源于浮点舍入如果误差明显偏大大概率是筛分循环里h的更新顺序写错了。跑通这一步之后就可以用真实图像替代合成图观察每一层 BIMF 的纹理尺度是否按预期从细到粗排布。4. BEMD 参数怎么调、边界怎么处理、分解质量怎么判断4.1 四个必调参数及推荐取值BEMD 需要人工干预的参数就四个极值窗口win、单次筛分最大迭代数max_iter、SD 目标值sd_goal、分解层数n_imfs。其余如插值方法选定后一般不再改动。参数推荐值控制什么调大/调小的后果win37极值检测邻域决定最细纹理尺度调大BIMF1 吞掉更大结构调小单像素噪声成为极值第一层被噪声主导max_iter2050每次筛分的轮数上限太小BIMF 里残留包络均值模态不纯太大过筛且耗时成倍增加sd_goal0.30.5筛分收敛阈值越小越严格但二维插值误差导致实际收益有限建议不要低于 0.2n_imfs35分解层数越多残差越干净但后层极值稀疏包络插值失真加剧win和sd_goal是耦合的win加大后极值点变少包络面更平缓SD 收敛更快此时sd_goal可以适当收紧到 0.3反过来win3时极值点多每一轮筛分的改动量大sd_goal0.5更实际。我一般先用win3, sd_goal0.3, max_iter20跑一遍看 BIMF1 是否混入了明显的中频结构再决定是否把win升到 5。4.2 边界效应与反射填充griddata的cubic插值只能在 Delaunay 三角剖分的凸包范围内计算凸包外的点全部由fill_value兜底。这导致分解结果在图像四角和边缘出现一整片均值平坦区下一轮筛分在这个区域里找不到任何极值点边界失真会一层层向内侵蚀。处理这个问题的常见做法是分解前对图像做反射填充分解完再裁掉填充区域pad 16 padded np.pad(test_img, pad, modereflect) pimfs, presidue bemd_decompose(padded, n_imfs3, win3, max_iter20, sd_goal0.3) imfs [pimfs[k][pad:-pad, pad:-pad] for k in range(len(pimfs))] residue presidue[pad:-pad, pad:-pad]我倾向于用反射填充而不是常数或对称填充因为常数填充会在边界制造一条明显的灰度跳变直接变成 BIMF1 里的假纹理对称填充在图像边界本身不连续时会让边界两侧的极值点被认为相邻产生虚假的包络振荡。反射填充相当于让图像在边界处光滑延拓对极值分布的影响最小。填充宽度取 1020 像素即可超过极值窗口 3 倍就足够了。4.3 平坦区域与极值板结灰度完全相同的平坦区域会给极值检测带来一个隐蔽问题整块区域都等于邻域最大值find_extrema会把每一个像素都标记为极大点极值点数量瞬间爆炸插值结果也毫无意义。这种情况在医学影像和遥感图的均匀背景区非常常见处理手段是对极值掩膜做连通域收缩只保留每个平板区域的中心点from scipy.ndimage import label def thin_plateau(mask): 把连成片的极值掩膜收缩成每块区域的中心点 lab, n label(mask) centers [] for i in range(1, n 1): ys, xs np.where(lab i) centers.append((int(ys.mean()), int(xs.mean()))) arr np.array(centers) if arr.size 0: return np.array([], dtypeint), np.array([], dtypeint) return arr[:, 1], arr[:, 0] # 返回 x, ylabel把连通的极值区域赋予同一个编号取每块区域的质心作为代表极值点。这个函数替换掉find_extrema里np.where那一步之后极值板结问题就消失了。实际处理中还要注意质心的坐标是浮点均值转整型可能出现两个相邻区域质心跳到同一个像素的情况筛分循环里对极值点做一次去重更稳妥。4.4 用两个硬指标判断分解质量分解质量不能只看重建误差还要看层间泄漏和极值递减规律。两个指标我每跑一组参数都会算一是重建误差 RMSE验证分解自洽二是正交性指标 OI衡量各层 BIMF 之间的信息重叠程度。OI 的定义是不同层逐像素乘积的绝对值之和除以分解总和能量的平方工程上小于 0.05 说明分解质量不错超过 0.1 就需要检查参数。此外还有一个经验规律每层 BIMF 的极值点数应该逐层明显下降如果 BIMF2 的极值点数和 BIMF1 差不多说明win太小或sd_goal太松第一层没有把该吸收的高频成分吸收干净。计数方法直接复用find_extrema(imf, win)的返回长度即可。5. BEMD 图像去噪与纹理分离一个直接能用的场景5.1 噪声落在哪一层取决于窗口和纹理尺度加性高斯噪声的尺度是单像素级在win3时绝大部分噪声极值会在第一轮就被吸收进 BIMF1。如果图像本身带有比噪声更细的纹理两者会一起出现在 BIMF1 里此时直接丢弃整层会连纹理一起抹掉。区分方法是看 BIMF1 的均值纯噪声层的均值接近零纹理层的均值往往有明显的空间结构。处理这类情况时对 BIMF1 再做一次阈值收缩而不是整层丢弃效果通常更好。5.2 去噪重建与正交性验证代码def orthogonality_index(imfs, residue): OI 越小说明各层信息重叠越少 terms imfs [residue] recon sum(imfs, residue) total 0.0 for i in range(len(terms)): for j in range(i 1, len(terms)): total np.abs(np.sum(terms[i] * terms[j])) return total * 2.0 / (np.sum(recon ** 2) 1e-12) imfs, residue bemd_decompose(noisy_img, n_imfs3, win3, max_iter20, sd_goal0.3) denoised residue sum(imfs[1:]) # 丢弃 BIMF1保留其余层 oi orthogonality_index(imfs, residue) print(正交性指标 OI:, oi)丢弃 BIMF1 后直接用sum(imfs[1:])叠加剩余层和残差得到的就是去噪结果。OI 的值在 0.01 到 0.05 之间属于正常范围如果大于 0.1说明相邻 BIMF 之间出现了频带交叠主要原因是sd_goal偏大导致筛分提前退出模态没有充分分离把阈值从 0.5 收到 0.3 重跑即可。5.3 一套可以直接复制的参数组合实际项目中我用的默认组合是win3, n_imfs3, max_iter20, sd_goal0.3配合 16 像素的反射填充。对 512×512 的灰度图这个配置单次完整分解在普通笔记本上大约几十秒OI 通常在 0.03 以下。如果第一层纹理混入严重把win升到 5 并把sd_goal收到 0.25代价是分解时间增加约一倍。如果想要更强的去噪效果可以对 BIMF1 做软阈值收缩后再重建imfs[0] * (np.abs(imfs[0]) t)其中t取 BIMF1 标准差的 1.5 倍这比整层丢弃更柔和也保留了噪声与纹理混叠时的高频细节。本文还有配套的精品资源点击获取