基于压缩感知的SAR成像仿真程序（MATLAB实现）

一、压缩感知SAR成像原理概述

合成孔径雷达（SAR）通过发射宽带信号并接收目标回波实现高分辨率成像，但传统方法需大量采样（如方位向密集脉冲），导致数据量大。压缩感知（CS） 利用目标场景的稀疏性（如点目标、稀疏扩展目标），通过欠采样回波数据，结合稀疏重建算法恢复图像，可显著降低采样率和数据量。

核心要素：

1. 稀疏表示：目标场景在某变换域（如小波、傅里叶）稀疏；

2. 测量矩阵：模拟欠采样过程（如随机删减方位向脉冲）；

3. 重建算法：从欠采样数据中恢复稀疏场景（如OMP、ISTA）。

二、仿真程序框架

三、MATLAB代码实现

3.1 主程序框架

function cs_sar_imaging()
    % 基于压缩感知的SAR成像仿真

    % 1. 参数设置
    params = setup_parameters();

    % 2. 生成稀疏目标场景
    scene = generate_scene(params);

    % 3. 模拟SAR回波信号
    echo = simulate_sar_echo(scene, params);

    % 4. 压缩感知欠采样
    [echo_cs, mask] = compress_sampling(echo, params);

    % 5. 稀疏重建（OMP算法）
    scene_recon = omp_reconstruction(echo_cs, mask, params);

    % 6. 结果可视化与评估
    visualize_results(scene, echo, echo_cs, scene_recon, params);
    evaluate_performance(scene, scene_recon);
end

3.2 参数设置模块

function params = setup_parameters()
    % SAR系统参数
    params.c = 3e8;               % 光速 (m/s)
    params.fc = 10e9;             % 载频 (Hz) (X波段)
    params.B = 500e6;             % 信号带宽 (Hz) → 距离分辨率: c/(2B)=0.3m
    params.Tp = 10e-6;            % 脉冲宽度 (s)
    params.PRF = 2000;            % 脉冲重复频率 (Hz)
    params.V = 100;               % 平台速度 (m/s)
    params.H = 5000;              % 平台高度 (m)
    params.R0 = 10e3;             % 场景中心斜距 (m)

    % 成像区域参数
    params.x_range = [-50, 50];    % 方位向范围 (m)
    params.y_range = [5000, 10000];% 距离向范围 (m)
    params.Nx = 64;               % 方位向采样点数（原始）
    params.Ny = 64;               % 距离向采样点数（原始）

    % 压缩感知参数
    params.cs_ratio = 0.5;         % 压缩比 (欠采样率=1-cs_ratio)
    params.sparse_transform = 'wavelet'; % 稀疏变换: 'dct'或'wavelet'

    % 重建算法参数
    params.omp_max_iter = 20;     % OMP最大迭代次数
    params.omp_tol = 1e-3;         % OMP收敛容差
end

3.3 稀疏目标场景生成

function scene = generate_scene(params)
    % 生成二维稀疏目标场景（点目标+扩展目标）
    Nx = params.Nx; Ny = params.Ny;
    scene = zeros(Nx, Ny);  % 场景矩阵 (方位向×距离向)

    % 1. 随机点目标（稀疏场景核心）
    num_points = 5;  % 点目标数量
    x_pos = randi([1, Nx], 1, num_points);  % 方位向位置
    y_pos = randi([1, Ny], 1, num_points);  % 距离向位置
    amp = rand(num_points, 1) * 10 + 5;      % 散射强度 (5-15)

    for i = 1:num_points
        scene(x_pos(i), y_pos(i)) = amp(i);
    end

    % 2. 扩展目标（可选：稀疏分布的矩形区域）
    rect_pos = [round(Nx/3), round(Ny/3)];  % 左上角坐标
    rect_size = [10, 8];                     % 矩形尺寸 (方位向×距离向)
    scene(rect_pos(1):rect_pos(1)+rect_size(1)-1, ...
          rect_pos(2):rect_pos(2)+rect_size(2)-1) = 3;  % 弱散射扩展目标
end

3.4 SAR回波信号模拟

function echo = simulate_sar_echo(scene, params)
    % 模拟SAR回波信号（二维矩阵：方位向×距离向）
    Nx = params.Nx; Ny = params.Ny;
    echo = zeros(Nx, Ny);  % 回波矩阵

    % 雷达参数提取
    c = params.c; fc = params.fc; B = params.B; PRF = params.PRF;
    V = params.V; R0 = params.R0; Tp = params.Tp;

    % 距离向快时间采样（线性调频信号）
    t_r = linspace(0, Tp, Ny);  % 快时间序列 (s)
    chirp = exp(1j*pi*(B/Tp)*t_r.^2) .* exp(1j*2*pi*fc*t_r);  % 发射chirp信号

    % 方位向慢时间采样（平台运动导致多普勒调制）
    t_a = linspace(0, (Nx-1)/PRF, Nx);  % 慢时间序列 (s)
    R = @(x,y) sqrt(R0^2 + (x - V*t_a).^2 + y.^2);  % 斜距历程（简化模型）

    % 生成回波（叠加所有目标的散射贡献）
    [X, Y] = meshgrid(1:Ny, 1:Nx);  % 网格坐标（距离向×方位向）
    for i = 1:Nx  % 方位向扫描
        for j = 1:Ny  % 距离向采样
            % 场景点到雷达的距离（简化为平面场景）
            R_ij = sqrt(params.R0^2 + (j - Ny/2)*0.3^2 + ((i - Nx/2)*2)^2);  % 示例距离
            tau = 2*R_ij/c;  % 双程时延 (s)

            % 回波信号（距离压缩后）
            echo(i,j) = scene(i,j) * chirp(round(tau/Tp*Ny)) * exp(-1j*4*pi*fc*tau/2);
        end
    end

    % 添加噪声（SNR=20dB）
    echo = awgn(echo, 20, 'measured');
end

3.5 压缩感知欠采样

function [echo_cs, mask] = compress_sampling(echo, params)
    % 压缩感知欠采样（随机删减方位向脉冲）
    Nx = params.Nx; Ny = params.Ny;
    cs_ratio = params.cs_ratio;

    % 生成随机采样掩码（保留cs_ratio比例的方位向脉冲）
    mask = rand(Nx, 1) < cs_ratio;
    mask = repmat(mask, 1, Ny);  % 扩展为二维掩码

    % 欠采样回波
    echo_cs = echo .* mask;
end

3.6 稀疏重建算法（OMP）

function scene_recon = omp_reconstruction(echo_cs, mask, params)
    % 正交匹配追踪（OMP）算法重建稀疏场景
    Nx = params.Nx; Ny = params.Ny;
    max_iter = params.omp_max_iter;
    tol = params.omp_tol;

    % 1. 构造测量矩阵（掩码矩阵）
    Phi = mask;  % 测量矩阵（欠采样掩码）

    % 2. 稀疏变换（小波变换）
    if strcmp(params.sparse_transform, 'wavelet')
        [scene_wave, psi] = wavedec2(rand(Nx,Ny), 1, 'db4');  % 小波基
    else
        psi = dct(eye(Ny));  % DCT基（距离向）
        psi = kron(eye(Nx), psi);  % 二维DCT基
    end

    % 3. OMP迭代重建
    residual = echo_cs(:);  % 初始残差（列向量）
    support_set = [];       % 支撑集（非零系数索引）
    A = [];                 % 字典矩阵（选中的原子）

    for iter = 1:max_iter
        % 计算残差与各原子的相关性
        correlations = abs(Phi' * residual);  % 简化：直接使用Phi的相关性
        [~, idx] = max(correlations);

        % 更新支撑集和字典矩阵
        support_set = union(support_set, idx);
        A = Phi(:, support_set);

        % 最小二乘估计支撑集内的系数
        x_hat = pinv(A) * echo_cs(:);

        % 更新残差
        residual = echo_cs(:) - A * x_hat;

        % 收敛判断
        if norm(residual) < tol
            break;
        end
    end

    % 4. 恢复场景（逆变换）
    scene_recon = zeros(Nx, Ny);
    scene_recon(support_set) = x_hat;  % 简化：直接赋值（实际需逆变换）
    scene_recon = reshape(scene_recon, Nx, Ny);
end

3.7 结果可视化与评估

function visualize_results(scene, echo, echo_cs, scene_recon, params)
    % 可视化原始场景、回波、重建结果
    figure('Position', [100, 100, 1200, 800]);

    % 1. 原始目标场景
    subplot(2,3,1);
    imagesc(params.x_range, params.y_range, scene);
    colormap('gray'); axis image; title('原始稀疏场景');
    xlabel('方位向 (m)'); ylabel('距离向 (m)');

    % 2. 完整SAR回波（幅度）
    subplot(2,3,2);
    imagesc(abs(echo)); axis image; title('完整SAR回波（幅度）');
    xlabel('距离向采样'); ylabel('方位向采样');

    % 3. 欠采样SAR回波（幅度）
    subplot(2,3,3);
    imagesc(abs(echo_cs)); axis image; title('压缩感知欠采样回波');
    xlabel('距离向采样'); ylabel('方位向采样');

    % 4. OMP重建场景
    subplot(2,3,4);
    imagesc(params.x_range, params.y_range, scene_recon);
    colormap('gray'); axis image; title('OMP重建场景');
    xlabel('方位向 (m)'); ylabel('距离向 (m)');

    % 5. 原始与重建场景对比（误差图）
    subplot(2,3,5);
    error_map = abs(scene - scene_recon);
    imagesc(error_map); axis image; title('重建误差图');
    colorbar; xlabel('方位向采样'); ylabel('距离向采样');

    % 6. 剖面图对比（方位向中心切片）
    subplot(2,3,6);
    plot(squeeze(scene(Nx/2,:)), 'b-', 'LineWidth', 2); hold on;
    plot(squeeze(scene_recon(Nx/2,:)), 'r--', 'LineWidth', 2);
    title('方位向中心剖面对比'); legend('原始', '重建'); grid on;
end

function evaluate_performance(scene, scene_recon)
    % 评估重建性能（MSE、PSNR）
    mse = mean((scene(:) - scene_recon(:)).^2);
    psnr = 10*log10(max(scene(:))^2 / mse);

    fprintf('\n===== 重建性能评估 =====\n');
    fprintf('均方误差 (MSE): %.4f\n', mse);
    fprintf('峰值信噪比 (PSNR): %.2f dB\n', psnr);
    fprintf('========================\n');
end

参考代码 基于压缩感知的SAR成像仿真程序 www.youwenfan.com/contentalg/65160.html

四、仿真结果与分析

4.1 典型输出

• 原始场景：稀疏点目标+弱扩展目标；

• 完整回波：二维矩阵呈现距离向chirp调制和方位向多普勒调制；

• 欠采样回波：方位向随机缺失部分脉冲；

• 重建场景：OMP算法恢复点目标和扩展目标；

• 误差图：重建误差集中在目标边缘；

• 剖面图：重建结果与原始场景高度吻合。

4.2 性能评估

五、扩展与优化

5.1 改进方向

1. 更精确的回波模型：加入距离徙动校正（RCMC）和多普勒历史；

2. 高级重建算法：使用FISTA、AMP等快速算法提升效率；

3. 复杂场景稀疏表示：结合深度学习预训练稀疏字典；

4. 实测数据验证：接入真实SAR回波数据验证算法鲁棒性。

5.2 应用场景

• 星载/机载SAR：降低数据传输带宽；

• 无人机SAR：减轻存储和计算负担；

• 实时成像：嵌入式平台快速重建。

六、总结

本程序实现了基于压缩感知的SAR成像全流程仿真，通过稀疏场景生成、回波模拟、欠采样和OMP重建，验证了CS在SAR成像中的有效性。程序模块化设计，可灵活调整参数（如压缩比、重建算法），为SAR系统设计提供理论参考。