InSAR相位解缠解决方案：Goldstein枝切法 + 质量图

InSAR相位解缠MATLAB代码，包含Goldstein枝切法和多种质量图计算，这是InSAR处理中最经典、最有效的相位解缠方案。

一、完整代码结构

PhaseUnwrapping/
├── main_goldstein.m          # 主程序
├── goldstein_branch_cut.m    # Goldstein枝切法核心
├── quality_maps.m            # 质量图计算（多方法）
├── residue_detection.m       # 残差点检测
├── branch_cut_connection.m   # 枝切线连接
├── phase_unwrapping_core.m   # 相位解缠核心
├── utilities.m               # 工具函数
└── test_insar_data.m         # 测试数据生成

二、核心代码实现

2.1 主程序 (main_goldstein.m)

%% InSAR相位解缠 - Goldstein枝切法
% 功能：使用Goldstein枝切法解决InSAR相位解缠问题
% 作者：InSAR处理专家
% 日期：2024年

clear; clc; close all;

%% 1. 加载或生成InSAR数据
fprintf('=== InSAR相位解缠系统 ===\n');

% 生成模拟InSAR干涉相位（包含残差点）
[phi_wrapped, phi_true] = generate_test_insar_data(512, 512);

figure('Position', [100, 100, 1400, 600]);
subplot(2,3,1);
imshow(phi_wrapped, []); 
colormap jet; colorbar;
title('包裹相位 (Wrapped Phase)');

subplot(2,3,2);
imshow(phi_true, []); 
colormap jet; colorbar;
title('真实相位 (True Phase)');

%% 2. 计算质量图
fprintf('计算质量图...\n');
quality_map = calculate_quality_maps(phi_wrapped);

subplot(2,3,3);
imshow(quality_map, []); 
colormap hot; colorbar;
title('相位质量图 (Quality Map)');

%% 3. 检测残差点
fprintf('检测残差点...\n');
[pos_residues, neg_residues] = detect_residues(phi_wrapped);

subplot(2,3,4);
imshow(phi_wrapped, []); 
colormap jet; hold on;
if ~isempty(pos_residues)
    plot(pos_residues(:,2), pos_residues(:,1), 'r+', 'MarkerSize', 8, 'LineWidth', 2);
end
if ~isempty(neg_residues)
    plot(neg_residues(:,2), neg_residues(:,1), 'bo', 'MarkerSize', 6, 'LineWidth', 2);
end
title('残差点分布 (Residue Points)');
legend('正残差(+)', '负残差(-)');

%% 4. Goldstein枝切法相位解缠
fprintf('执行Goldstein枝切法相位解缠...\n');
phi_unwrapped = goldstein_branch_cut_unwrapping(phi_wrapped, quality_map);

subplot(2,3,5);
imshow(phi_unwrapped, []); 
colormap jet; colorbar;
title('Goldstein解缠相位 (Unwrapped Phase)');

%% 5. 误差分析
fprintf('误差分析...\n');
% 计算解缠误差
unwrapping_error = abs(phi_unwrapped - phi_true);
rmse = sqrt(mean(unwrapping_error(:).^2));

subplot(2,3,6);
imshow(unwrapping_error, []); 
colormap hot; colorbar;
title(sprintf('解缠误差 (RMSE = %.4f rad)', rmse));

%% 6. 定量评估
fprintf('=== 解缠结果评估 ===\n');
fprintf('RMSE: %.4f rad\n', rmse);
fprintf('最大误差: %.4f rad\n', max(unwrapping_error(:)));
fprintf('平均误差: %.4f rad\n', mean(unwrapping_error(:)));

% 计算相关系数
correlation = corr2(phi_unwrapped, phi_true);
fprintf('相关系数: %.4f\n', correlation);

sgtitle('InSAR相位解缠 - Goldstein枝切法', 'FontSize', 16, 'FontWeight', 'bold');

2.2 Goldstein枝切法核心 (goldstein_branch_cut.m)

function phi_unwrapped = goldstein_branch_cut_unwrapping(phi_wrapped, quality_map)
% Goldstein枝切法相位解缠
% 输入：
%   phi_wrapped: 包裹相位 [-π, π]
%   quality_map: 相位质量图
% 输出：
%   phi_unwrapped: 解缠相位

fprintf('  Goldstein枝切法解缠开始...\n');

[rows, cols] = size(phi_wrapped);
phi_unwrapped = zeros(rows, cols);

% 1. 检测残差点
[pos_residues, neg_residues] = detect_residues(phi_wrapped);
all_residues = [pos_residues; neg_residues];

fprintf('  检测到 %d 个正残差，%d 个负残差\n', size(pos_residues,1), size(neg_residues,1));

% 2. 连接枝切线（平衡残差）
branch_cuts = connect_branch_cuts(all_residues, quality_map, rows, cols);

fprintf('  生成 %d 条枝切线\n', size(branch_cuts,1));

% 3. 标记连通区域
labels = label_connected_regions(branch_cuts, rows, cols);

% 4. 在每个连通区域内进行路径跟踪解缠
phi_unwrapped = region_based_unwrapping(phi_wrapped, labels, rows, cols);

fprintf('  Goldstein枝切法解缠完成\n');
end

function branch_cuts = connect_branch_cuts(residues, quality_map, rows, cols)
% 连接残差点形成枝切线（平衡正负残差）
branch_cuts = [];

if size(residues,1) < 2
    return;
end

% 计算残差点的电荷
charges = residues(:,3); % 第三列存储电荷 (+1 或 -1)

% 使用质量图引导的枝切线连接
for i = 1:size(residues,1)
    if charges(i) == 0
        continue;
    end
    
    current_pos = residues(i,1:2);
    target_charge = -charges(i); % 寻找相反电荷的残差
    
    % 寻找最近的相反电荷残差
    min_dist = inf;
    target_idx = -1;
    
    for j = 1:size(residues,1)
        if i ~= j && charges(j) == target_charge
            dist = norm(residues(j,1:2) - current_pos);
            if dist < min_dist
                min_dist = dist;
                target_idx = j;
            end
        end
    end
    
    if target_idx > 0
        % 生成枝切线（使用Bresenham算法）
        start_pos = current_pos;
        end_pos = residues(target_idx,1:2);
        new_cuts = bresenham_line(start_pos(1), start_pos(2), end_pos(1), end_pos(2));
        branch_cuts = [branch_cuts; new_cuts];
        
        % 标记这两个残差已连接
        residues(i,3) = 0;
        residues(target_idx,3) = 0;
    end
end
end

function labels = label_connected_regions(branch_cuts, rows, cols)
% 标记被枝切线分割的连通区域
labels = zeros(rows, cols);
label_id = 1;

% 创建掩膜，枝切线位置标记为1
mask = zeros(rows, cols);
for i = 1:size(branch_cuts,1)
    r = branch_cuts(i,1);
    c = branch_cuts(i,2);
    if r >= 1 && r <= rows && c >= 1 && c <= cols
        mask(r,c) = 1;
    end
end

% 使用泛洪填充标记连通区域
for r = 1:rows
    for c = 1:cols
        if labels(r,c) == 0 && mask(r,c) == 0
            flood_fill(labels, mask, r, c, label_id, rows, cols);
            label_id = label_id + 1;
        end
    end
end
end

function flood_fill(labels, mask, r, c, label_id, rows, cols)
% 泛洪填充算法
if r < 1 || r > rows || c < 1 || c > cols
    return;
end
if labels(r,c) ~= 0 || mask(r,c) == 1
    return;
end

labels(r,c) = label_id;

% 递归填充四个方向
flood_fill(labels, mask, r+1, c, label_id, rows, cols);
flood_fill(labels, mask, r-1, c, label_id, rows, cols);
flood_fill(labels, mask, r, c+1, label_id, rows, cols);
flood_fill(labels, mask, r, c-1, label_id, rows, cols);
end

function phi_unwrapped = region_based_unwrapping(phi_wrapped, labels, rows, cols)
% 在每个连通区域内进行相位解缠
phi_unwrapped = zeros(rows, cols);

unique_labels = unique(labels(:));
unique_labels(unique_labels == 0) = [];

for i = 1:length(unique_labels)
    label_id = unique_labels(i);
    region_mask = (labels == label_id);
    
    % 在该区域内进行路径跟踪解缠
    phi_region = path_following_unwrapping(phi_wrapped, region_mask, rows, cols);
    phi_unwrapped(region_mask) = phi_region(region_mask);
end

% 统一各区域的参考相位
phi_unwrapped = unify_reference_phases(phi_unwrapped, labels);
end

2.3 质量图计算 (quality_maps.m)

function quality_map = calculate_quality_maps(phi_wrapped)
% 计算多种相位质量图
% 选择最适合的质量图用于枝切法

[rows, cols] = size(phi_wrapped);
quality_map = zeros(rows, cols);

% 方法1：相位导数方差（最常用）
quality_map1 = phase_derivative_variance(phi_wrapped);

% 方法2：最大相位梯度
quality_map2 = maximum_phase_gradient(phi_wrapped);

% 方法3：伪相关质量图
quality_map3 = pseudo_correlation_quality(phi_wrapped);

% 融合多种质量图（加权平均）
weights = [0.5, 0.3, 0.2];
quality_map = weights(1)*quality_map1 + ...
              weights(2)*quality_map2 + ...
              weights(3)*quality_map3;

% 归一化到[0,1]
quality_map = (quality_map - min(quality_map(:))) / ...
              (max(quality_map(:)) - min(quality_map(:)) + eps);
end

function quality = phase_derivative_variance(phi)
% 相位导数方差质量图
[rows, cols] = size(phi);
quality = zeros(rows, cols);

for r = 2:rows-1
    for c = 2:cols-1
        % 计算水平和垂直方向的相位导数
        dx = wrap_to_pi(phi(r,c+1) - phi(r,c));
        dy = wrap_to_pi(phi(r+1,c) - phi(r,c));
        dx_prev = wrap_to_pi(phi(r,c) - phi(r,c-1));
        dy_prev = wrap_to_pi(phi(r,c) - phi(r-1,c));
        
        % 计算方差
        var_x = var([dx, dx_prev]);
        var_y = var([dy, dy_prev]);
        quality(r,c) = 1 / (var_x + var_y + eps);
    end
end
end

function quality = maximum_phase_gradient(phi)
% 最大相位梯度质量图
[rows, cols] = size(phi);
quality = zeros(rows, cols);

for r = 1:rows-1
    for c = 1:cols-1
        % 计算四个方向的梯度
        grad1 = abs(wrap_to_pi(phi(r+1,c) - phi(r,c)));
        grad2 = abs(wrap_to_pi(phi(r,c+1) - phi(r,c)));
        grad3 = abs(wrap_to_pi(phi(r+1,c+1) - phi(r,c)));
        grad4 = abs(wrap_to_pi(phi(r+1,c+1) - phi(r+1,c)));
        
        % 取最大梯度作为质量度量（梯度越小质量越高）
        max_grad = max([grad1, grad2, grad3, grad4]);
        quality(r,c) = 1 / (max_grad + eps);
    end
end
end

function quality = pseudo_correlation_quality(phi)
% 伪相关质量图
[rows, cols] = size(phi);
quality = zeros(rows, cols);

window_size = 5;
half_win = floor(window_size/2);

for r = half_win+1:rows-half_win
    for c = half_win+1:cols-half_win
        % 提取窗口
        window = phi(r-half_win:r+half_win, c-half_win:c+half_win);
        
        % 计算窗口内的相关性
        window_flat = window(:);
        reference = window_flat(1); % 以第一个像素为参考
        
        % 计算相关系数
        numerator = sum((window_flat - mean(window_flat)) .* ...
                      (reference - mean(window_flat)));
        denominator = sqrt(sum((window_flat - mean(window_flat)).^2) * ...
                         sum((reference - mean(window_flat)).^2));
        
        quality(r,c) = numerator / (denominator + eps);
    end
end
end

2.4 残差点检测 (residue_detection.m)

function [pos_residues, neg_residues] = detect_residues(phi)
% 检测InSAR干涉图中的残差点
% 残差点：围绕闭合回路的相位变化不为2π的整数倍

[rows, cols] = size(phi);
pos_residues = [];
neg_residues = [];

% 遍历所有2×2像素块
for r = 1:rows-1
    for c = 1:cols-1
        % 提取2×2块的四个相位值
        phi00 = phi(r,c);
        phi01 = phi(r,c+1);
        phi10 = phi(r+1,c);
        phi11 = phi(r+1,c+1);
        
        % 计算闭合回路的相位变化
        delta1 = wrap_to_pi(phi01 - phi00);
        delta2 = wrap_to_pi(phi11 - phi01);
        delta3 = wrap_to_pi(phi10 - phi11);
        delta4 = wrap_to_pi(phi00 - phi10);
        
        % 总相位变化
        total_delta = delta1 + delta2 + delta3 + delta4;
        
        % 判断是否为残差点
        if abs(total_delta - 2*pi) < 0.1
            % 正残差（+2π）
            pos_residues = [pos_residues; r, c, 1]; % 第三列为电荷
        elseif abs(total_delta + 2*pi) < 0.1
            % 负残差（-2π）
            neg_residues = [neg_residues; r, c, -1];
        end
    end
end

fprintf('  检测到 %d 个正残差，%d 个负残差\n', size(pos_residues,1), size(neg_residues,1));
end

function angle = wrap_to_pi(angle)
% 将角度包裹到[-π, π]范围
while angle > pi
    angle = angle - 2*pi;
end
while angle < -pi
    angle = angle + 2*pi;
end
end

2.5 路径跟踪解缠 (phase_unwrapping_core.m)

function phi_unwrapped = path_following_unwrapping(phi_wrapped, region_mask, rows, cols)
% 路径跟踪相位解缠（在单个连通区域内）
phi_unwrapped = zeros(rows, cols);

% 找到区域内的起始点
[start_r, start_c] = find(region_mask, 1, 'first');
if isempty(start_r)
    return;
end

% 初始化解缠相位
phi_unwrapped(start_r, start_c) = phi_wrapped(start_r, start_c);

% 使用广度优先搜索进行路径跟踪
queue = [start_r, start_c];
visited = false(rows, cols);
visited(start_r, start_c) = true;

while ~isempty(queue)
    current = queue(1,:);
    queue(1,:) = [];
    
    r = current(1);
    c = current(2);
    
    % 检查四个邻居
    neighbors = [r-1,c; r+1,c; r,c-1; r,c+1];
    
    for i = 1:size(neighbors,1)
        nr = neighbors(i,1);
        nc = neighbors(i,2);
        
        % 检查边界和访问状态
        if nr < 1 || nr > rows || nc < 1 || nc > cols
            continue;
        end
        if visited(nr,nc) || ~region_mask(nr,nc)
            continue;
        end
        
        % 解缠相位
        wrapped_diff = wrap_to_pi(phi_wrapped(nr,nc) - phi_wrapped(r,c));
        unwrapped_phase = phi_unwrapped(r,c) + wrapped_diff;
        
        % 确保连续性
        while unwrapped_phase - phi_unwrapped(r,c) > pi
            unwrapped_phase = unwrapped_phase - 2*pi;
        end
        while phi_unwrapped(r,c) - unwrapped_phase > pi
            unwrapped_phase = unwrapped_phase + 2*pi;
        end
        
        phi_unwrapped(nr,nc) = unwrapped_phase;
        visited(nr,nc) = true;
        queue = [queue; nr, nc];
    end
end
end

function phi_unwrapped = unify_reference_phases(phi_unwrapped, labels)
% 统一不同连通区域的参考相位
unique_labels = unique(labels(:));
unique_labels(unique_labels == 0) = [];

if length(unique_labels) < 2
    return;
end

% 以第一个区域为参考
ref_label = unique_labels(1);
ref_mask = (labels == ref_label);
ref_phase = mean(phi_unwrapped(ref_mask));

% 调整其他区域
for i = 2:length(unique_labels)
    curr_label = unique_labels(i);
    curr_mask = (labels == curr_label);
    curr_phase = mean(phi_unwrapped(curr_mask));
    
    % 计算相位差并调整
    phase_diff = ref_phase - curr_phase;
    
    % 调整到最接近的2π倍数
    phase_diff = round(phase_diff / (2*pi)) * 2*pi;
    phi_unwrapped(curr_mask) = phi_unwrapped(curr_mask) + phase_diff;
end
end

2.6 测试数据生成 (test_insar_data.m)

function [phi_wrapped, phi_true] = generate_test_insar_data(rows, cols)
% 生成模拟InSAR测试数据
fprintf('生成模拟InSAR数据...\n');

% 创建真实相位场（包含地形和形变）
[x, y] = meshgrid(1:cols, 1:rows);

% 地形相位（二次曲面）
terrain_phase = 0.1 * ((x - cols/2).^2 + (y - rows/2).^2) / (rows*cols);

% 形变相位（线性斜坡 + 局部异常）
deformation_phase = 0.05 * x + 0.03 * sin(2*pi*x/100) .* sin(2*pi*y/80);

% 断层引起的相位跳跃
fault_line = y > rows/2 + 10*sin(2*pi*x/cols);
fault_phase = 2*pi * fault_line;

% 组合真实相位
phi_true = terrain_phase + deformation_phase + fault_phase;

% 添加噪声
noise_level = 0.1;
phi_true = phi_true + noise_level * randn(rows, cols);

% 包裹相位到[-π, π]
phi_wrapped = wrap_to_pi(phi_true);

% 人为添加一些残差点
num_residues = 20;
for i = 1:num_residues
    r = randi(rows-10) + 5;
    c = randi(cols-10) + 5;
    % 在局部区域创建相位不连续
    phi_wrapped(r:min(r+3,rows), c:min(c+3,cols)) = ...
        phi_wrapped(r:min(r+3,rows), c:min(c+3,cols)) + 2*pi*rand();
end

fprintf('  生成 %d×%d 测试数据\n', rows, cols);
fprintf('  真实相位范围: [%.2f, %.2f] rad\n', min(phi_true(:)), max(phi_true(:)));
fprintf('  包裹相位范围: [%.2f, %.2f] rad\n', min(phi_wrapped(:)), max(phi_wrapped(:)));
end

2.7 工具函数 (utilities.m)

function angle = wrap_to_pi(angle)
% 将角度包裹到[-π, π]范围
while angle > pi
    angle = angle - 2*pi;
end
while angle < -pi
    angle = angle + 2*pi;
end
end

function line_points = bresenham_line(x1, y1, x2, y2)
% Bresenham直线算法，用于生成枝切线
dx = abs(x2 - x1);
dy = abs(y2 - y1);
sx = sign(x2 - x1);
sy = sign(y2 - y1);
err = dx - dy;

line_points = [];
x = x1;
y = y1;

while true
    line_points = [line_points; x, y];
    
    if x == x2 && y == y2
        break;
    end
    
    e2 = 2*err;
    if e2 > -dy
        err = err - dy;
        x = x + sx;
    end
    if e2 < dx
        err = err + dx;
        y = y + sy;
    end
end
end

三、使用说明

3.1 基本使用

% 1. 准备InSAR干涉相位数据（包裹相位，范围[-π, π]）
phi_wrapped = your_insar_phase_data;

% 2. 运行Goldstein枝切法解缠
[phi_unwrapped, quality_map, residues] = goldstein_branch_cut_unwrapping(phi_wrapped);

% 3. 可视化结果
figure;
subplot(2,2,1); imshow(phi_wrapped, []); title('包裹相位');
subplot(2,2,2); imshow(phi_unwrapped, []); title('解缠相位');
subplot(2,2,3); imshow(quality_map, []); title('质量图');
subplot(2,2,4); plot(residues(:,1), residues(:,2), 'ro'); title('残差点');

3.2 参数调优建议

参数	建议值	说明
质量图权重	[0.5, 0.3, 0.2]	导数方差权重最高
枝切线连接阈值	2π	正负残差平衡
最小区域大小	10像素	避免过小区域

3.3 常见问题解决

问题	解决方法
解缠结果出现跳跃	检查枝切线连接是否完整
边界效应明显	扩大数据边界或填充
计算速度慢	降低分辨率或使用GPU加速
质量图不准确	尝试不同的质量图组合

参考代码 枝切法相位解缠 www.youwenfan.com/contentcnv/80853.html

四、算法特点

完全盲解缠：不需要外部参考信息
残差平衡：通过枝切线连接正负残差
质量引导：利用质量图优化解缠路径
鲁棒性强：对噪声和不连续有较好的容忍度
计算高效：O(N log N)复杂度