简介：一份基于枝切法的相位解包裹 MATLAB 代码，面向干涉测量、SAR 成像、光学计量与电子全息等领域的科研人员和工程师，解决因噪声、断点或欠采样导致的相位跳变恢复难题。代码完整实现经典枝切法流程：首先通过残差检测定位正负极点，再依据邻近原则生成枝切线加以连接，随后以洪水填充策略绕开枝切线进行相位展开，兼顾精度与效率。相比最小二乘法，枝切法能有效抑制噪声传播，保持局部相位细节，适合复杂场形解包裹。同时附带椭圆裁剪模块，可自动屏蔽背景与无效区域，主程序集成各步骤，便于一键测试与二次开发。资源共 14 个文件，以 5 个 .m 源文件为核心，配合 9 幅 .bmp 相位图样例，压缩包大小约 962KB，目录清晰、注释明确，适合算法验证、课程实验及工程落地参考。已有 242 人学习下载，是掌握枝切法原理与实现细节的实用资料。

1. 枝切法解包裹的 Matlab 实现：找对障碍，比积分更重要

做白光干涉测量或者 InSAR 形变提取的时候，相位解包裹始终是最玄学的一环：包裹相位图看着一片平滑，直接对相邻像素做差分再累加，结果解出来的相位却东一块西一块飘着斜条纹，怎么也拼不回连续曲面。真正翻车的点不在积分本身，而在于积分路径撞上了那些残差点——相位场上不可积的地方。这套基于枝切法的 matlab 解包裹代码，核心就是先找残差、用枝切线搭桥、再绕开枝切线做路径积分。它适合正在做干涉测量、全息、结构光三维重建的开发者，拿回去能直接跑，也能按自己的数据改掩膜和连接策略。

2. 算法原理与代码结构：残差、枝切线、路径积分三位一体

2.1 核心流程：为什么解包裹要先找「坏点」

从数学上看，解包裹要做的就是把每个像素的相位加上适当的2π整数倍，恢复连续的真实相位。最直觉的做法是逐行扫描：从左到右相邻像素做差，差值超过π就补一个2π。这个方法在无噪声、无遮挡的理想数据上没问题，可实际采集的相位里经常有噪声、阴影、低相干区域，这些地方相邻像素的差分折叠到(-π, π]之后累加一圈并不为零，形成了一个「残差」。

残差相当于积分路径上的一个旋涡：你从左边跨过它和从右边跨过它，得到的积分结果会差2π。如果不做处理，解包裹结果就会沿着某条路径突然跳变出一整块偏移，看起来像相位面被撕裂了。

枝切法的思路是：把正负残差用「枝切线」连接起来，让每一条枝切线上的电荷（正负残差个数）相互抵消，然后规划一条不穿越任何枝切线的积分路径。只要路径不跨过这些障碍，积分结果就和路径无关，解包裹就稳定了。

这套代码在 Matlab 里的执行顺序是：读入包裹相位wrapped，生成掩膜mask，计算残差图residue，生成枝切线branch_map，最后在掩膜内从种子点开始洪水填充积分。整个数据流是单向的，非常适合按模块调试。

2.2 工程目录与数据流：这份代码包怎么组织

代码包按「一主流程、四子模块、一数据目录」的方式组织，主脚本只管串流程，算法全部下沉到utils下的独立函数，这样你要替换枝切线策略或者换积分方式，都不需要动主流程。

我一般建议先单独跑02_calc_residue.m看残差分布，再跑03_build_branchcut.m看枝切线画在哪，最后才整体跑完。这样每一步都能看到中间状态，查起来也快。

2.3 相位输入与掩膜：从复数干涉图到可用相位

枝切法要求输入是单通道的包裹相位，单位是弧度，范围在(-π, π]。很多第一次用的人直接拿强度图（条纹图）往里塞，结果残差图一团乱。正确做法是先对复数干涉图取幅角：

% 假设 complex_img 是复数干涉图 wrapped = angle(complex_img); % 得到 [-pi, pi] 的包裹相位 mask = (abs(complex_img) > 0.05); % 低幅值区域通常不可信，进入掩膜

这里的mask是 0/1 矩阵，1 表示该像素参与解包裹。掩膜的作用非常大：噪声大、遮挡、阴影区域的相位本来就是乱的，如果让这些像素参与残差计算，会冒出一大片假残差，枝切线的长度会急剧膨胀，影响全局结果。常见做法是用幅值阈值生成初始掩膜，再配合形态学开运算把零散的小孔去掉：

mask = imopen(mask, strel('disk', 5));

在需要更精细控制时，我会额外引入质量图，比如用相位梯度的倒数作为质量指标，把质量过低的像素并入掩膜排除。质量图在第三步选种子、第四步做积分优先级的场景里也会用到。

3. 核心模块实现：残差计算与枝切线生成的完整代码

3.1 残差极性计算：2×2 环路积分

残差的计算是枝切法的基础。对每个 2×2 的像素块，沿着环路的四条边做相位差分并折叠到(-π, π]，再累加。理想平滑区域累加值接近零；如果某个角上存在2π跳变，累加值就会落在+2π或-2π附近，据此判定该点为正残差或负残差。代码如下：

function residue = calc_residue(wrapped, mask) % 输入: % wrapped - MxN double, 包裹相位, 单位弧度 % mask - MxN logical, 1 表示像素有效 % 输出: % residue - MxN int8, +1 / -1 / 0 [M, N] = size(wrapped); residue = zeros(M, N, 'int8'); for j = 1:N-1 for i = 1:M-1 % 四邻域全覆盖才参与计算, 否则跳过 if mask(i,j) && mask(i,j+1) && mask(i+1,j) && mask(i+1,j+1) d1 = wrap_pi(wrapped(i, j+1) - wrapped(i, j)); d2 = wrap_pi(wrapped(i+1, j+1) - wrapped(i, j+1)); d3 = wrap_pi(wrapped(i+1, j ) - wrapped(i+1, j+1)); d4 = wrap_pi(wrapped(i, j ) - wrapped(i+1, j )); s = d1 + d2 + d3 + d4; if s > 2 residue(i,j) = 1; % 正残差 elseif s < -2 residue(i,j) = -1; % 负残差 end end end end end

这里wrap_pi的作用是把任意角度差折叠到主值区间，我常用atan2(sin(x), cos(x))实现，而不是依赖某个工具箱函数，跨机器跑不会报错。判定阈值为什么取2而不是π？因为理论上总和应该是0或±2π，约等于±6.28，取2是给噪声留一点余量；如果你发现残差过多，可以把阈值抬高到3甚至4，过滤掉弱噪声引起的虚假残差。参数说明里还有一点要留意：残差坐标落在 2×2 块的左上角像素上，后续画枝切线时按这个约定对齐，不要混用。

3.2 枝切线生成：最近邻贪心连接与电荷平衡

枝切线的目标是用线段把正负残差连起来，使每一组内正负电荷之和为零。最常见的是贪心最近邻策略：每次取一个残差作为起点，找最近的、能平衡当前组电荷的残差，用直线连接并合并到同一组；如果残差周围找不到可配对的点，就允许把枝切线延伸到边界，因为边界上的积分路径可以绕过去，不会造成闭环误差。

先把残差点坐标和电荷拆出来：

function branch_map = build_branchcut(residue, max_len) % 输入: % residue - MxN int8, 残差极性图 % max_len - 单段枝切线最大长度, 建议 15~30, 防止跨过整幅图 % 输出: % branch_map - MxN logical, true 表示不可穿越 [M, N] = size(residue); [ys, xs] = find(residue ~= 0); charge = residue(residue ~= 0); n = numel(xs); branch_map = false(M, N); used = zeros(n, 1); % 组号, 0 表示未处理 group_charge = []; % 每组的净电荷 % 初始每个残差单独成组 for k = 1:n used(k) = k; group_charge(k, 1) = charge(k); end

上面的初始化比较直接，真正重要的是合并逻辑：每次在所有未配对完整的组里，找一个电荷不为零的组，取该组最近的一个点，去扫描另一组里最近的残差，判断连接后那组的净电荷能否向零靠近。合并后要更新group_charge，并把两点之间的直线像素写入branch_map。这里用到了 Bresenham 直线生成，避免枝切线穿过像素时出现断点：

function line_pixels = bresenham_line(x0, y0, x1, y1) % 从 (x0,y0) 到 (x1,y1) 的像素坐标列表 dx = abs(x1 - x0); dy = -abs(y1 - y0); sx = sign(x1 - x0); sy = sign(y1 - y0); err = dx + dy; line_pixels = [x0, y0]; while ~(x0 == x1 && y0 == y1) e2 = 2 * err; if e2 >= dy err = err + dy; x0 = x0 + sx; end if e2 <= dx err = err + dx; y0 = y0 + sy; end line_pixels(end+1, :) = [x0, y0]; end end

max_len这个参数很关键。最近邻贪心容易犯一个毛病：距离优先，结果可能从图像左上角连到右下角，形成一条超长切口，把高置信区域也隔断了。加一个长度上限，让超长连接改为连到边界，是工程上非常管用的约束。连到边界本身不算错误，边界上的缺口不会形成闭环，解包裹结果依然稳定。

3.3 掩膜与质量信息：如何降低假残差干扰

残差计算里有一个常被忽视的细节：mask边缘的像素参与环路积分时，四邻域全有效才会纳入，这会导致掩膜最外一圈的有效残差值被压掉。我的习惯是先对掩膜做一次imerode(mask, strel('disk', 1))，让残差计算向内收缩一圈，保证参与环路积分的都是有效邻域。

如果质量图可用，常见做法是把质量太差的像素直接并入 mask 排除。质量图我一般用相位梯度模值的倒数构造：

gw = hypot(gradient(wrapped), 0.5); % 梯度模值，加 0.5 防止除零 quality = 1 ./ gw; mask = mask & (quality > quantile(quality(mask), 0.25));

这个quantile的 25% 分位是一个不好拍脑袋定的阈值，先看直方图再选比较稳妥。把低质量区域屏蔽掉之后，残差数量通常会减少一半以上，枝切线数量也跟着下降，解包裹结果会干净很多。

4. 相位解包裹与路径积分：在障碍物下恢复真实相位

4.1 洪水填充：从种子点绕开枝切线做积分

枝切线一旦生成，剩下的工作就是在掩膜内、绕过枝切线做路径积分。枝切法的特色是路径无关性：所有不跨枝切线的路径，积分结果一致，所以可以用最简单的 BFS 填充，不需要质量引导。核心逻辑是维护一个队列，从种子点向四邻域扩散，逐点累加相邻相位差。

function unwrapped = flood_fill(wrapped, mask, branch_map, seed) % 输入: % wrapped - MxN double % mask - MxN logical, 有效区域 % branch_map - MxN logical, true 表示枝切线, 不可跨越 % seed - [row, col] 起始像素 % 输出: % unwrapped - MxN double, 解包裹相位 [M, N] = size(wrapped); unwrapped = zeros(M, N); visited = false(M, N); queue = seed; visited(seed(1), seed(2)) = true; dirs = [-1 0; 1 0; 0 -1; 0 1]; while ~isempty(queue) p = queue(1, :); queue(1, :) = []; for d = 1:size(dirs, 1) ni = p(1) + dirs(d, 1); nj = p(2) + dirs(d, 2); if ni < 1 || nj < 1 || ni > M || nj > N continue; end if visited(ni, nj) || ~mask(ni, nj) || branch_map(ni, nj) continue; end % 相邻相位差折叠到主值区间, 这一步是解包裹的关键 delta = wrap_pi(wrapped(ni, nj) - wrapped(p(1), p(2))); unwrapped(ni, nj) = unwrapped(p(1), p(2)) + delta; visited(ni, nj) = true; queue(end+1, :) = [ni, nj]; end end end

这里每一行的注释都在讲一个关键点：delta先把相邻差折叠到2π主值区间再累加，等价于「补上合适的 2π 整数倍」。branch_map(ni, nj)的判断放在visited之后也很重要——如果某段枝切线只有 1 像素宽，从另一侧绕过来访问时仍然能正常通过，不会被堵死。种子点的选择一般取掩膜中心附近质量最高的像素，避免从边界噪音区域开始扩散导致误差传播到一个大扇区。

4.2 积分路径的效果评估：用相邻跳变占比校验

解包裹跑完，不能只看图「像不像」。我一般跑一个客观指标：统计结果里相邻像素相位差超过π的像素占比。理想情况下解包裹相位应该是连续的，这个比例应该非常低；如果占比很高，说明大概率有未连接的残差或者枝切线穿透了有效区。

function jump_ratio = phase_jump_ratio(phase) % 计算相邻相位差超过 pi 的像素占比, 用于质量评估 dx = abs(diff(phase, 1, 2)); dy = abs(diff(phase, 1, 1)); jumps = sum(dx(:) > pi) + sum(dy(:) > pi); total = (size(phase,1)-1) * size(phase,2) + ... size(phase,1) * (size(phase,2)-1); jump_ratio = jumps / total; end

这个指标特别适合在换参数之后做横向对比：同样一张图，枝切线max_len从 10 改成 50，jump_ratio是变好还是变差，一测就知道，不用每版都肉眼盯着看。配合可视化脚本把残差图、枝切线、解包裹图三张图并排输出，问题出一眼能定位。

4.3 可视化对比：哪些区域最容易被解坏

可视化我倾向于用subplot排四张图：原始包裹相位、残差极性图（spy或者imagesc标红）、枝切线覆盖图、解包裹相位图。残差图和枝切线图叠在一起看，能快速发现两类问题：一类是残差密度过高的区域，说明掩膜没挡住噪声；另一类是枝切线横穿高置信区，说明贪心连接策略需要加max_len限制。

figure; subplot(2,2,1); imagesc(wrapped); title('wrapped phase'); axis image; colorbar; subplot(2,2,2); imagesc(residue); title('residue (+1/-1)'); axis image; colorbar; subplot(2,2,3); imagesc(branch_map); title('branch cuts'); axis image; colorbar; subplot(2,2,4); imagesc(unwrapped); title('unwrapped'); axis image; colorbar; colormap jet;

实际调参时，我通常先看第二张子图的残差有没有聚集在某个物理边界上，再看第三张的枝切线长度是否可控。如果残差图上有大片同号残差，说明掩膜阈值太松，噪声相位被当成了信号；如果枝切线又长又横穿图像，把max_len调小，让长线改连边界，通常能明显降低jump_ratio。

5. 避坑指南：跑枝切法时最常见的 5 个问题

5.1 残差图全是噪点，密集到看不见结构

现象：calc_residue.m跑完，残差图上密密麻麻全是正负残差，完全看不出真实相位缺陷的位置。

原因：输入相位噪声太大，或者掩膜没有排除低幅值、阴影区域。相位噪声一多，2×2 环路积分处处不闭合，残差分不出真假。

解决：先用medfilt2(wrapped, [3 3])或者imgaussfilt(wrapped, 0.8)对包裹相位做一次轻度平滑，再重新计算残差。掩膜侧再用幅值阈值加形态学开运算过滤掉孤立小区域，残差数量通常能下降 70% 以上。如果还嫌多，把残差判定阈值从2提到3。

5.2 枝切线过长，解包裹结果出现长条状偏移

现象：枝切线很长，横跨图像中间，解包裹结果在枝切线两侧的高度明显差了一截，像地形图上出现一道断崖。

原因：贪心最近邻策略只按距离连接，不区分高低质量区域，把一个残差拉到了很远的配对点，切口切过了大片高质量区域。

解决：给build_branchcut加上max_len参数，超过阈值的连接改为连到边界；或者在距离里叠加质量惩罚项，让枝切线尽量走在低质量区域。

5.3 解包裹结果在边框附近出现同心圆状伪影

现象：结果图边缘有一圈圈阶梯状条纹，越往外越明显，好像相位被多绕了几圈。

原因：残差计算只覆盖了掩膜内部的(M-1)×(N-1)个环路，掩膜最外一圈的残差没有被识别，导致枝切线在边界处没有闭合，积分路径在边界附近绕出了额外的一圈。

解决：对掩膜先做padarray(mask, [1 1], 'replicate')，包裹相位同样扩展一圈，再算残差；或者干脆在掩膜外预留一圈 1 像素的安全边界，保证最小环路都落在有效掩膜之内。

5.4 换了数据源之后，整套代码跑出来的全是 NaN

现象：新的数据加载进来，flood_fill返回的矩阵里出现成片NaN，可视化直接空白。

原因：常见两个——掩膜把种子点所在区域置为 0，种子点一开始就无效；或者包裹相位里有NaN没有被提前处理，wrap_pi传播了NaN。

解决：在跑主流程前加一个断言，检查种子点处mask(seed) == 1；再对输入做一次wrapped(isnan(wrapped)) = 0; mask(isnan(mask)) = 0;或干脆报错提示，让数据问题第一时间暴露。

5.5 相邻差直接算出来是 3π，补值后越补越乱

现象：有人图省事不折叠相位差，直接delta = wrapped(ni,nj) - wrapped(...)然后累加，结果越界位置越跳越大。

原因：没有理解包裹相位的差分必须折叠到(-π, π]才能反映「最小相位差」。原值差分在跳变边缘会出现> π或<-π，累加积分完全失真，翻车属于必然。

解决：所有相邻差一律走wrap_pi，也就是atan2(sin(diff), cos(diff))。这个函数是整个路径积分的灵魂，我在封装时特意把它独立成utils/wrap_pi.m，每次调用都不会忘。

6. 二次开发思路：把固定策略改成质量引导和自适应枝切线

6.1 改造点一：BFS 换成质量优先队列

枝切法本身不挑积分路径，但你可以把第四章的flood_fill改成质量引导版本：不再用先进先出队列，而是每次从边界里挑质量值最高的像素展开。这样噪声区域的误差就不会沿着 BFS 的「波纹状」扩散到整个画面，很适合掩膜复杂、低质量区域比较多的数据。改动很小，把queue换成按quality排序的优先队列即可，delta和避障逻辑完全不动。

6.2 改造点二：最近邻连接换成最小生成树

贪心最近邻容易拉出长枝切，工程上另一个常用方案是用最小生成树把残差连接起来。Matlab 里可以先构造graph对象，再用minspantree找最小生成树，最后用bresenham_line把每条边画到branch_map上。连接权重可以用欧氏距离乘以质量惩罚项，这样枝切线既短又好走。这个改法对残差多而杂的数据提升很明显。

6.3 改造点三：自适应max_len代替单一阈值

把固定的max_len改成根据残差点密度自动估计——比如取残差平均间距的 1.5 倍作为阈值。密度高时阈值小、连线短，密度低时阈值大、允许连到边界。这样能减少大量人工调参工序，数据源变化时也不至于立刻崩。

多年以后我每次跑枝切法解包裹，都会强制自己先看一眼残差图和jump_ratio两个中间结果再谈最终输出。这个习惯来自一次给某课题组的干涉数据跑解包，跳过检查直接出图，结果整片结果图里藏着一道贯穿性错位，差点把分析结论带偏。从那以后，中间产物检查成了我的固定动作。希望帮到你，少走一圈弯路。

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