合成孔径雷达（SAR）成像：线性调频（LFM）信号原理与MATLAB仿真

线性调频（LFM）脉冲是合成孔径雷达（SAR）实现高分辨率成像的核心技术。

一、 LFM信号：SAR高分辨率的数学基础

1.1 为什么需要LFM信号？

SAR通过发射大时宽-带宽积信号解决经典雷达的“分辨率-探测距离”矛盾：

• 大时宽：提高平均功率，增加探测距离
• 大带宽：通过脉冲压缩获得高距离分辨率

1.2 LFM信号的数学表达

发射信号（基带表示）：
$s_t(t) = \text{rect}\left(\frac{t}{T_p}\right) \exp\left(j\pi K t^2\right)$

接收信号（经过目标反射）：
$s_r(t) = A \cdot \text{rect}\left(\frac{t-\tau}{T_p}\right) \exp\left[j\pi K (t-\tau)^2\right] \exp(-j2\pi f_c\tau)$

其中：

• $T_p$：脉冲宽度
• $K = B/T_p$：调频率（带宽/脉宽）
• $B$：信号带宽
• $\tau = 2R/c$：目标双程时延
• $f_c$：载波频率
• $\text{rect}(u) = 1$ 当 $|u| \leq 0.5$

1.3 脉冲压缩原理

通过匹配滤波将宽脉冲压缩为窄脉冲：
$s_{pc}(t) = s_r(t) \otimes h_{mf}(t)$
其中匹配滤波器 $h_{mf}(t) = s_t^*(-t)$

压缩后脉冲宽度：
$\Delta \tau = \frac{1}{B}$

距离分辨率：
$\rho_r = \frac{c}{2B}$

---

二、 MATLAB仿真：从LFM信号生成到SAR成像

2.1 LFM信号生成与脉冲压缩仿真

%% SAR LFM信号生成与脉冲压缩仿真
clear; clc; close all;

%% 1. 参数设置
c = 3e8;                % 光速 (m/s)
fc = 10e9;              % 载频 10GHz (X波段)
lambda = c/fc;          % 波长

% LFM信号参数
Tp = 10e-6;             % 脉冲宽度 10us
Br = 100e6;             % 带宽 100MHz
Kr = Br/Tp;             % 调频率
Fs = 2*Br;              % 采样率 (满足Nyquist)
Ts = 1/Fs;              % 采样间隔

% 目标参数
R0 = 10000;             % 目标距离 10km
RCS = 1;                % 目标雷达截面积

%% 2. 生成LFM发射信号
t = -Tp/2 : Ts : Tp/2 - Ts;    % 快时间轴
N = length(t);

% 发射信号 (基带)
st = exp(1j*pi*Kr*t.^2);       % 线性调频信号
st_envelope = ones(size(t));   % 矩形包络

% 添加载频 (实际发射信号)
st_rf = real(st .* exp(1j*2*pi*fc*t));

figure('Position', [100, 100, 1200, 800]);

% 时域波形
subplot(3,3,1);
plot(t*1e6, real(st), 'b', 'LineWidth', 1.5);
xlabel('时间 (\mus)'); ylabel('幅度');
title('LFM基带信号 (实部)'); grid on;

% 瞬时频率
subplot(3,3,2);
inst_freq = Kr*t/(2*pi);  % 瞬时频率 = d(phase)/dt / (2π)
plot(t*1e6, inst_freq/1e6, 'r', 'LineWidth', 1.5);
xlabel('时间 (\mus)'); ylabel('频率 (MHz)');
title('LFM信号瞬时频率'); grid on;

% 频谱
subplot(3,3,3);
Nfft = 2^nextpow2(N);
freq_axis = (-Nfft/2:Nfft/2-1)*Fs/Nfft;
ST = fftshift(fft(st, Nfft));
plot(freq_axis/1e6, 20*log10(abs(ST)/max(abs(ST))), 'g', 'LineWidth', 1.5);
xlabel('频率 (MHz)'); ylabel('幅度 (dB)');
title('LFM信号频谱'); grid on; xlim([-Br/2e6, Br/2e6]);

%% 3. 模拟目标回波
tau0 = 2*R0/c;          % 目标双程时延
t_delay = t + tau0;     % 延迟后的时间轴

% 回波信号 (忽略幅度衰减和载频)
sr = exp(1j*pi*Kr*(t_delay).^2);

% 添加噪声
SNR_dB = 20;            % 信噪比
signal_power = mean(abs(sr).^2);
noise_power = signal_power / (10^(SNR_dB/10));
noise = sqrt(noise_power/2) * (randn(size(sr)) + 1j*randn(size(sr)));
sr_noisy = sr + noise;

subplot(3,3,4);
plot(t*1e6, real(sr_noisy), 'b', 'LineWidth', 1.5);
xlabel('时间 (\mus)'); ylabel('幅度');
title('含噪回波信号 (实部)'); grid on;

%% 4. 脉冲压缩处理 (频域匹配滤波)
% 匹配滤波器参考函数
t_ref = t;
s_ref = exp(1j*pi*Kr*t_ref.^2);    % 与发射信号共轭
H_ref = conj(fft(s_ref, Nfft));    % 匹配滤波器频域响应

% 回波信号FFT
SR = fft(sr_noisy, Nfft);

% 频域匹配滤波
S_out = SR .* H_ref;

% IFFT得到压缩结果
s_out = ifft(S_out, Nfft);
s_out = s_out(1:N);     % 取有效长度

% 时间轴对应距离轴
range_axis = t * c / 2;  % 距离 = 时间 * 光速 / 2

subplot(3,3,5);
plot(range_axis, abs(s_out)/max(abs(s_out)), 'r', 'LineWidth', 2);
xlabel('距离 (m)'); ylabel('归一化幅度');
title('脉冲压缩结果'); grid on;
xlim([R0-50, R0+50]);

% 绘制对数坐标
subplot(3,3,6);
plot(range_axis, 20*log10(abs(s_out)/max(abs(s_out))), 'r', 'LineWidth', 2);
xlabel('距离 (m)'); ylabel('幅度 (dB)');
title('脉冲压缩结果 (dB)'); grid on;
xlim([R0-50, R0+50]); ylim([-60, 0]);

%% 5. 分析脉冲压缩性能
% 理论分辨率
delta_R = c/(2*Br);
fprintf('理论距离分辨率: %.2f m\n', delta_R);

% 测量主瓣宽度 (3dB宽度)
[peak_val, peak_idx] = max(abs(s_out));
peak_range = range_axis(peak_idx);

% 寻找3dB点
threshold = peak_val/sqrt(2);  % -3dB
left_idx = find(abs(s_out(1:peak_idx)) >= threshold, 1, 'first');
right_idx = find(abs(s_out(peak_idx:end)) >= threshold, 1, 'last') + peak_idx - 1;

if ~isempty(left_idx) && ~isempty(right_idx)
    measured_width = range_axis(right_idx) - range_axis(left_idx);
    fprintf('实测主瓣宽度: %.2f m\n', measured_width);
    
    % 标注3dB宽度
    subplot(3,3,6); hold on;
    plot([range_axis(left_idx), range_axis(right_idx)], ...
         [20*log10(threshold/max(abs(s_out))), 20*log10(threshold/max(abs(s_out)))], ...
         'k--', 'LineWidth', 1.5);
    text(mean([range_axis(left_idx), range_axis(right_idx)]), -5, ...
         sprintf('%.1f m', measured_width), 'HorizontalAlignment', 'center');
end

%% 6. 旁瓣分析
% 计算峰值旁瓣比 (PSLR)
mainlobe_region = (peak_idx-10):(peak_idx+10);
sidelobe_mask = true(size(s_out));
sidelobe_mask(mainlobe_region) = false;

sidelobe_max = max(abs(s_out(sidelobe_mask)));
PSLR_dB = 20*log10(sidelobe_max/peak_val);
fprintf('峰值旁瓣比 (PSLR): %.2f dB\n', PSLR_dB);

% 绘制旁瓣区域
subplot(3,3,7);
sidelobe_range = range_axis;
sidelobe_range(mainlobe_region) = NaN;
plot(sidelobe_range, 20*log10(abs(s_out)/max(abs(s_out))), 'b', 'LineWidth', 1.5);
xlabel('距离 (m)'); ylabel('幅度 (dB)');
title('旁瓣特性'); grid on;
xlim([R0-200, R0+200]); ylim([-60, 0]);

%% 7. 加窗处理降低旁瓣
% 使用Hamming窗
window = hamming(N)';
s_ref_windowed = s_ref .* window;
H_ref_windowed = conj(fft(s_ref_windowed, Nfft));

% 加窗后的脉冲压缩
S_out_windowed = SR .* H_ref_windowed;
s_out_windowed = ifft(S_out_windowed, Nfft);
s_out_windowed = s_out_windowed(1:N);

subplot(3,3,8);
plot(range_axis, 20*log10(abs(s_out)/max(abs(s_out))), 'r', 'LineWidth', 1.5); hold on;
plot(range_axis, 20*log10(abs(s_out_windowed)/max(abs(s_out_windowed))), 'b', 'LineWidth', 1.5);
xlabel('距离 (m)'); ylabel('幅度 (dB)');
title('加窗前后对比'); grid on;
legend('不加窗', 'Hamming窗');
xlim([R0-50, R0+50]); ylim([-80, 0]);

% 计算加窗后的PSLR
sidelobe_max_windowed = max(abs(s_out_windowed(sidelobe_mask)));
PSLR_windowed_dB = 20*log10(sidelobe_max_windowed/max(abs(s_out_windowed)));
fprintf('加窗后PSLR: %.2f dB\n', PSLR_windowed_dB);

%% 8. 距离分辨率验证
% 模拟两个相近目标
R1 = 10000;             % 目标1距离
R2 = 10015;             % 目标2距离 (间隔15m)
tau1 = 2*R1/c;
tau2 = 2*R2/c;

% 两个目标的回波
sr1 = exp(1j*pi*Kr*(t + tau1).^2);
sr2 = exp(1j*pi*Kr*(t + tau2).^2);
sr_two = sr1 + sr2 + noise;

% 脉冲压缩
SR_two = fft(sr_two, Nfft);
S_out_two = SR_two .* H_ref;
s_out_two = ifft(S_out_two, Nfft);
s_out_two = s_out_two(1:N);

subplot(3,3,9);
plot(range_axis, abs(s_out_two)/max(abs(s_out_two)), 'm', 'LineWidth', 2);
xlabel('距离 (m)'); ylabel('归一化幅度');
title('两个相近目标的区分'); grid on;
xlim([R1-50, R2+50]);

% 检查是否能分辨
[peaks, locs] = findpeaks(abs(s_out_two), 'MinPeakHeight', 0.3);
if length(peaks) >= 2
    peak_positions = range_axis(locs(1:2));
    separation = abs(diff(peak_positions));
    fprintf('两个目标实际间距: %.1f m\n', R2-R1);
    fprintf('测量间距: %.1f m\n', separation);
    
    if separation > delta_R
        fprintf('结论: 可以分辨 (间距 > 分辨率)\n');
    else
        fprintf('结论: 难以分辨 (间距 < 分辨率)\n');
    end
end

sgtitle('SAR LFM信号脉冲压缩仿真', 'FontSize', 14, 'FontWeight', 'bold');

2.2 SAR成像仿真（简化版距离多普勒算法）

%% 简化SAR成像仿真（距离多普勒算法）
clear; clc; close all;

%% 1. 成像几何与参数
c = 3e8;
fc = 5.3e9;                 % C波段
lambda = c/fc;

% SAR平台参数
V = 100;                    % 平台速度 100 m/s
H = 5000;                   % 飞行高度 5000 m
R0 = 10000;                 % 中心斜距

% 信号参数
Br = 50e6;                  % 距离向带宽
Tp = 5e-6;                  % 脉冲宽度
Kr = Br/Tp;                 % 距离向调频率
PRF = 1000;                 % 脉冲重复频率
Ta = 2;                     % 合成孔径时间

% 目标场景（5个点目标）
targets = [0, 0, R0;        % 场景中心
           -20, 10, R0;     % 左前方
           20, 10, R0;      % 右前方
           -20, -10, R0;    % 左后方
           20, -10, R0];    % 右后方

%% 2. 生成SAR原始数据
% 时间轴设置
Fs_r = 1.2*Br;              % 距离向采样率
Ts_r = 1/Fs_r;
Nr = round(2*R0/c * Fs_r) + 512;  % 距离向采样点数
tr = (0:Nr-1)*Ts_r - Nr*Ts_r/2;   % 距离向快时间

Na = round(Ta * PRF);       % 方位向采样点数
ta = (0:Na-1)/PRF - Ta/2;   % 方位向慢时间

% 初始化原始数据矩阵
raw_data = zeros(Na, Nr);

fprintf('开始生成SAR原始数据...\n');
for target_idx = 1:size(targets, 1)
    fprintf('  处理目标 %d/%d\n', target_idx, size(targets, 1));
    
    x0 = targets(target_idx, 1);    % 方位向位置
    y0 = targets(target_idx, 2);    % 地面距离向位置
    Rc = targets(target_idx, 3);    % 中心斜距
    
    for i = 1:Na
        % 计算瞬时斜距
        R = sqrt(Rc^2 + (V*ta(i) - x0)^2);
        
        % 计算双程时延
        tau = 2*R/c;
        
        % 生成LFM回波（距离徙动忽略）
        phase = -4*pi/lambda * R + pi*Kr*(tr - tau).^2;
        sr = exp(1j*phase) .* (abs(tr - tau) <= Tp/2);
        
        % 添加到原始数据
        raw_data(i, :) = raw_data(i, :) + sr;
    end
end

% 添加噪声
SNR_dB = 15;
signal_power = mean(abs(raw_data(:)).^2);
noise_power = signal_power / (10^(SNR_dB/10));
noise = sqrt(noise_power/2) * (randn(size(raw_data)) + 1j*randn(size(raw_data)));
raw_data = raw_data + noise;

%% 3. 距离多普勒（RD）算法处理
fprintf('开始RD算法成像处理...\n');

% 步骤1：距离向脉冲压缩
fprintf('  步骤1: 距离向脉冲压缩\n');
% 生成距离向参考函数
t_ref = tr;
s_ref_r = exp(1j*pi*Kr*t_ref.^2) .* (abs(t_ref) <= Tp/2);
H_ref_r = conj(fft(s_ref_r, Nr));

% 距离向压缩（每行做FFT）
data_rc = zeros(size(raw_data));
for i = 1:Na
    row_fft = fft(raw_data(i, :), Nr);
    row_filtered = row_fft .* H_ref_r;
    data_rc(i, :) = ifft(row_filtered);
end

% 步骤2：方位向FFT（转到距离多普勒域）
fprintf('  步骤2: 方位向FFT\n');
data_rd = fft(data_rc, [], 1);

% 步骤3：距离徙动校正（RCMC - 简化版）
fprintf('  步骤3: 距离徙动校正\n');
% 计算多普勒频率轴
fa = (-Na/2:Na/2-1) * PRF/Na;  % 多普勒频率

% 计算距离徙动量
delta_R = lambda^2 * R0 * fa.^2 / (8*V^2);

% 插值校正（简化：只校正参考距离）
data_rcmc = zeros(size(data_rd));
for i = 1:Na
    % 计算该多普勒频率下的徙动量
    shift_samples = round(delta_R(i) / (c/(2*Fs_r)));
    
    if shift_samples > 0
        data_rcmc(i, 1:end-shift_samples) = data_rd(i, shift_samples+1:end);
    elseif shift_samples < 0
        data_rcmc(i, -shift_samples+1:end) = data_rd(i, 1:end+shift_samples);
    else
        data_rcmc(i, :) = data_rd(i, :);
    end
end

% 步骤4：方位向压缩
fprintf('  步骤4: 方位向压缩\n');
% 计算方位向调频率
Ka = 2*V^2 / (lambda * R0);

% 生成方位向参考函数
s_ref_a = exp(1j*pi*Ka*ta.^2);
H_ref_a = conj(fft(s_ref_a, Na));

% 方位向压缩（每列做IFFT）
data_ac = zeros(size(data_rcmc));
for j = 1:Nr
    col_fft = data_rcmc(:, j);
    col_filtered = col_fft .* H_ref_a.';
    data_ac(:, j) = ifft(col_filtered);
end

% 最终图像
image_data = fftshift(abs(data_ac));

%% 4. 结果显示
figure('Position', [100, 100, 1400, 600]);

% 原始数据
subplot(2,3,1);
imagesc(tr*c/2, ta*V, 20*log10(abs(raw_data)/max(abs(raw_data(:)))));
xlabel('距离 (m)'); ylabel('方位位置 (m)');
title('原始数据 (幅度)'); colorbar; colormap('jet');
caxis([-40, 0]); axis tight;

% 距离压缩后
subplot(2,3,2);
imagesc(tr*c/2, ta*V, 20*log10(abs(data_rc)/max(abs(data_rc(:)))));
xlabel('距离 (m)'); ylabel('方位位置 (m)');
title('距离压缩后'); colorbar; colormap('jet');
caxis([-40, 0]); axis tight;

% 距离多普勒域
subplot(2,3,3);
imagesc(tr*c/2, fa, 20*log10(abs(data_rd)/max(abs(data_rd(:)))));
xlabel('距离 (m)'); ylabel('多普勒频率 (Hz)');
title('距离多普勒域'); colorbar; colormap('jet');
caxis([-40, 0]); axis tight;

% 距离徙动校正后
subplot(2,3,4);
imagesc(tr*c/2, fa, 20*log10(abs(data_rcmc)/max(abs(data_rcmc(:)))));
xlabel('距离 (m)'); ylabel('多普勒频率 (Hz)');
title('RCMC后'); colorbar; colormap('jet');
caxis([-40, 0]); axis tight;

% 最终成像结果
subplot(2,3,5);
imagesc(tr*c/2, ta*V, image_data);
xlabel('距离 (m)'); ylabel('方位位置 (m)');
title('最终SAR图像'); colorbar; colormap('gray');
axis tight;

% 三维显示
subplot(2,3,6);
[X, Y] = meshgrid(tr*c/2, ta*V);
surf(X, Y, image_data, 'EdgeColor', 'none');
xlabel('距离 (m)'); ylabel('方位 (m)'); zlabel('幅度');
title('SAR图像3D显示'); colormap('jet');
view(45, 30); axis tight;

sgtitle('SAR成像仿真 - 距离多普勒算法', 'FontSize', 14, 'FontWeight', 'bold');

%% 5. 性能指标计算
fprintf('\n========== SAR成像性能指标 ==========\n');

% 距离分辨率
delta_r = c/(2*Br);
fprintf('距离分辨率: %.2f m\n', delta_r);

% 方位分辨率
delta_a = V/PRF;  % 简化计算
fprintf('方位分辨率: %.2f m\n', delta_a);

% 计算目标位置精度
[~, center_idx] = min(abs(tr*c/2 - R0));
[~, center_a_idx] = min(abs(ta*V));

% 提取中心目标响应
target_response = image_data(center_a_idx-50:center_a_idx+50, center_idx-50:center_idx+50);

% 计算峰值旁瓣比
[peak_val, peak_loc] = max(target_response(:));
[peak_row, peak_col] = ind2sub(size(target_response), peak_loc);

% 创建主瓣掩膜
mainlobe_mask = false(size(target_response));
mainlobe_radius = 5;
for i = 1:size(target_response,1)
    for j = 1:size(target_response,2)
        if sqrt((i-peak_row)^2 + (j-peak_col)^2) <= mainlobe_radius
            mainlobe_mask(i,j) = true;
        end
    end
end

% 计算旁瓣
sidelobe_values = target_response(~mainlobe_mask);
PSLR = 20*log10(max(sidelobe_values)/peak_val);
fprintf('峰值旁瓣比 (PSLR): %.2f dB\n', PSLR);

% 积分旁瓣比
mainlobe_power = sum(target_response(mainlobe_mask).^2);
sidelobe_power = sum(sidelobe_values.^2);
ISLR = 10*log10(sidelobe_power/mainlobe_power);
fprintf('积分旁瓣比 (ISLR): %.2f dB\n', ISLR);

---

三、 关键理论与工程要点

3.1 LFM信号的核心优势

1. 大时宽-带宽积：

   • 时宽 $T_p$ 决定能量：$E = P_{avg} \cdot T_p$
   • 带宽 $B$ 决定分辨率：$\rho_r = c/(2B)$
   • 通过脉冲压缩同时获得高能量和高分辨率

2. 多普勒容限：

   • LFM信号对多普勒频率不敏感
   • 适用于高速平台下的SAR成像

3.2 实际工程中的挑战与解决方案

% 实际工程中的改进措施示例

% 1. 运动补偿
function compensated_data = motion_compensation(raw_data, platform_trajectory)
    % 基于GPS/INS数据的运动补偿
    % platform_trajectory: [time, x, y, z, vx, vy, vz]
    % 计算每个脉冲的相位补偿项
    % ...
end

% 2. 自聚焦算法（相位误差校正）
function autofocused_data = autofocus(sar_image, method)
    % method: 'PGA' (相位梯度自聚焦) 或 'MapDrift'
    % 迭代估计并校正相位误差
    % ...
end

% 3. 多视处理（降低斑点噪声）
function multilook_image = multilook_processing(sar_image, nlooks)
    % 将全分辨率图像分为多个子视，然后非相干平均
    % nlooks: 视数
    % ...
end

3.3 现代SAR成像算法对比

算法	优点	缺点	适用场景
距离多普勒(RD)	计算效率高，实现简单	大斜视时精度下降	小斜视星载SAR
调频变标(CS)	处理大斜视数据精度高	计算复杂，参数敏感	机载SAR，大斜视
波数域(ωK)	理论上最精确	计算量最大，内存需求高	超高分辨率SAR
后向投影(BP)	最灵活，适用于任意轨迹	计算量极大，实时性差	前视SAR，曲线轨迹

3.4 扩展应用：极化SAR与干涉SAR

% 极化SAR数据处理框架
classdef PolSAR_Processor
    properties
        HH_data    % HH极化数据
        HV_data    % HV极化数据
        VH_data    % VH极化数据
        VV_data    % VV极化数据
    end
    
    methods
        function obj = PolSAR_Processor(hh, hv, vh, vv)
            % 初始化四极化数据
            obj.HH_data = hh;
            obj.HV_data = hv;
            obj.VH_data = vh;
            obj.VV_data = vv;
        end
        
        function [span, entropy] = compute_pol_parameters(obj)
            % 计算极化总功率和熵
            % 构建协方差矩阵
            % 特征值分解
            % 计算熵和散射机制
        end
    end
end

% 干涉SAR处理
function [height_map, coherence] = insar_processing(master_slc, slave_slc, baseline)
    % master_slc, slave_slc: 单视复图像
    % baseline: 空间基线
    
    % 1. 配准
    registered_slave = register_images(slave_slc, master_slc);
    
    % 2. 生成干涉图
    interferogram = master_slc .* conj(registered_slave);
    
    % 3. 去平地效应
    interferogram = remove_flat_earth(interferogram, baseline);
    
    % 4. 相位解缠
    unwrapped_phase = phase_unwrapping(angle(interferogram));
    
    % 5. 高程反演
    height_map = phase_to_height(unwrapped_phase, baseline);
    
    % 6. 计算相干性
    coherence = compute_coherence(master_slc, registered_slave);
end

参考代码 合成孔径雷达成像脉冲调频信号 www.youwenfan.com/contentcsu/65419.html

四、 运行与调试建议

1. 首次运行：

   • 先运行LFM脉冲压缩仿真，理解基本原理
   • 调整带宽 Br 观察分辨率变化
   • 修改信噪比 SNR_dB 观察噪声影响

2. SAR成像仿真：

   • 从简单点目标开始（代码中的5个点）
   • 逐步增加目标数量，观察成像效果
   • 调整平台速度 V 和高度 H，理解几何关系

3. 性能优化：

   % 使用GPU加速（如有NVIDIA GPU）
   if gpuDeviceCount > 0
       raw_data = gpuArray(raw_data);
       % ... 处理过程
       image_data = gather(image_data);
   end
   
   % 使用并行计算
   parfor i = 1:Na
       % 距离向处理
   end
   

4. 实际数据应用：

   • 替换模拟数据为实际SAR原始数据
   • 根据实际系统参数调整仿真参数
   • 添加实际几何校正和辐射定标