云计算百科
云计算领域专业知识百科平台

SAR成像原理和实现(含python代码)

目录

  • SAR成像基础
    • 简介
    • 基本工作流程
  • 线性调频信号与脉冲压缩
    • 距离向高分辨的关键
    • 方位向处理
  • SAR回波信号模型
    • 点目标回波
  • 距离徙动
    • 距离徙动(Range Cell Migration, RCM)
      • 距离徙动分量
    • RCM校正
  • 成像算法简介
    • 距离多普勒算法 (RDA)
      • 算法流程
    • Chirp Scaling 算法 (CSA)
      • 算法流程
    • ωK 算法 (波数域算法)
      • 算法流程
    • 后向投影算法 (BP)
    • 算法对比
    • 总结
  • 仿真
  • python代码(RDA算法)
  • 参考
    • 参考资料:

SAR成像基础

简介

SAR成像,全称为合成孔径雷达成像(Synthetic Aperture Radar)。SAR利用雷达平台(飞机、卫星、无人机等)的运动,在不同位置连续发射和接收信号。通过信号处理算法,把这些不同位置采集的回波“合成”出一个等效的超大天线(孔径),从而大幅提高方位向分辨率,达到与光学图像相当甚至更高的分辨率。 基本总结 通过距离压缩 + 方位压缩等信号处理算法,形成二维高分辨率图像。

  • 合成孔径雷达 (Synthetic Aperture Radar) 通过雷达平台运动合成虚拟大天线
  • 全天时、全天候获取高分辨率二维图像
  • 距离向分辨率:依靠发射大带宽信号(如线性调频)实现
  • 方位向分辨率:依靠合成孔径长度获得,与距离无关

基本工作流程

  • 距离向(垂直于飞行方向):雷达沿轨迹运动,周期性调频脉冲,通过回波的时间延迟和频率变化,分辨出目标远近,这叫脉冲压缩。
  • 方位向(沿飞行轨迹方向):利用平台运动,把同一目标在不同位置的回波进行匹配处理,等效形成长天线,这叫合成孔径处理。 SAR 成像工作示意图如下所示:

请添加图片描述

线性调频信号与脉冲压缩

距离向高分辨的关键

  • 发射信号:线性调频(LFM)脉冲

    s

    (

    τ

    )

    =

    rect

    (

    τ

    T

    p

    )

    exp

    (

    j

    2

    π

    f

    c

    τ

    +

    j

    π

    K

    r

    τ

    2

    )

    s(\\tau) = \\text{rect}\\left(\\frac{\\tau}{T_p}\\right) \\exp\\left(j2\\pi f_c\\tau + j\\pi K_r\\tau^2\\right)

    s(τ)=rect(Tpτ)exp(j2πfcτ+jπKrτ2)

    T

    p

    T_p

    Tp为脉冲宽度 ,

    f

    c

    f_c

    fc为载频,

    K

    r

    K_r

    Kr为距离向调频率,

    τ

    \\tau

    τ 为距离向脉冲时间, 通常取脉冲中心为参考原点.

  • 接收回波通过匹配滤波(脉冲压缩)实现距离向聚焦
  • 距离分辨率:

    δ

    r

    =

    c

    2

    B

    \\delta_r = \\frac{c}{2B}

    δr=2Bc

方位向处理

  • 方位向回波亦为近似线性调频信号
  • 利用多普勒历史进行方位压缩,获得方位分辨率

    δ

    a

    =

    L

    2

    \\delta_a = \\frac{L}{2}

    δa=2L

SAR回波信号模型

点目标回波

  • 经正交解调后,基带回波为二维信号:

    s

    0

    (

    τ

    ,

    η

    )

    =

    A

    w

    r

    (

    τ

    2

    R

    (

    η

    )

    c

    )

    w

    a

    (

    η

    η

    c

    )

    s_0(\\tau, \\eta) = A\\cdot w_r\\left(\\tau – \\frac{2R(\\eta)}{c}\\right) w_a(\\eta – \\eta_c)

    s0(τ,η)=Awr(τc2R(η))wa(ηηc)

    ×

    exp

    (

    j

    4

    π

    f

    c

    R

    (

    η

    )

    c

    )

    exp

    (

    j

    π

    K

    r

    (

    τ

    2

    R

    (

    η

    )

    c

    )

    2

    )

    \\times \\exp\\left(-j\\frac{4\\pi f_c R(\\eta)}{c}\\right) \\exp\\left(j\\pi K_r \\left(\\tau – \\frac{2R(\\eta)}{c}\\right)^2\\right)

    ×exp(jc4πfcR(η))exp(jπKr(τc2R(η))2)

  • 瞬时斜距:

    R

    (

    η

    )

    =

    R

    0

    2

    +

    v

    2

    η

    2

    R(\\eta) = \\sqrt{R_0^2 + v^2\\eta^2}

    R(η)=R02+v2η2

    (正侧视简化为双曲线)

在这里插入图片描述

在低斜视角下,

R

(

η

)

R(\\eta)

R(η) 可由菲涅尔近似为:

R

(

η

)

=

R

0

+

(

V

η

)

2

2

R

0

R(\\eta) = R_0 + \\frac{(V\\eta)^2}{2R_0}

R(η)=R0+2R0(Vη)2

距离徙动

距离徙动(Range Cell Migration, RCM)

  • 当雷达平台以恒定速度v飞行时,雷达与目标的瞬时斜距,会随着雷达平台的运动而不断变化。
  • 瞬时斜距:

    R

    (

    η

    )

    =

    R

    0

    2

    +

    v

    2

    η

    2

    R(\\eta) = \\sqrt{R_0^2 + v^2\\eta^2}

    R(η)=R02+v2η2

    (正侧视简化为双曲线)

  • 当这个斜距变化量超过一个距离分辨单元时,目标的回波包络就会在“距离向”(快时间)上跨单元移动,这就是距离徙动。
  • 结果就是:同一个点目标的回波,在多次脉冲间,出现在不同的距离采样点上,仿佛在“远、近”之间迁徙。 请添加图片描述

距离徙动分量

  • 距离走动(Range Walk) 由雷达与目标的相对径向速度引起的线性分量。在大斜视角下非常显著,回波轨迹是一条斜线。
  • 距离弯曲(Range Curvature) 斜距变化中的高阶(主要是二次)分量,使回波轨迹弯曲。正侧视时,走动为零,主要体现为弯曲;高分辨率、长波长或大场景边缘时弯曲更明显。
  • RCM校正

    如果不对距离徙动进行校正,点目标的能量会散布在若干个距离单元上,导致: · 方位向无法有效压缩(匹配滤波失配) · 图像出现严重散焦、几何失真和分辨率下降 · 目标旁瓣升高,信噪比恶化

    因此,距离徙动校正(RCMC) 是SAR成像算法(如距离多普勒算法RDA、Chirp Scaling算法等)中不可缺少的步骤,通常是在距离多普勒域进行距离插值, 可以基于

    (

    s

    i

    n

    c

    )

    (sinc)

    (sinc) 函数进行插值处理. 需要校正的 RCM 为方位频率

    (

    f

    η

    )

    (f\\eta)

    (fη)的函数, 也是

    (

    R

    0

    )

    (R_0)

    (R0)的函数:

    Δ

    R

    (

    f

    η

    )

    =

    λ

    2

    R

    0

    f

    η

    2

    8

    V

    r

    2

    \\Delta R(f_\\eta) = \\frac{\\lambda^2R_0f_\\eta^2}{8V_r^2}

    ΔR(fη)=8Vr2λ2R0fη2. 距离徙动校正因子:

    D

    (

    f

    η

    ,

    V

    r

    )

    =

    1

    c

    2

    f

    η

    2

    4

    V

    r

    2

    f

    0

    2

    =

    1

    λ

    2

    f

    η

    2

    4

    V

    r

    2

    D(f_{\\eta}, V_r) = \\sqrt{1-\\frac{c^2f_{\\eta}^2}{4V_r^2f_0^2}} = \\sqrt{1- \\frac{\\lambda^2f_{\\eta}^2 }{4V_r^2}}

    D(fη,Vr)=14Vr2f02c2fη2

    =14Vr2λ2fη2

    成像算法简介

    距离多普勒算法 (RDA)

    基本原理 距离多普勒算法(Range-Doppler Algorithm,RDA)是合成孔径雷达(SAR)成像中最经典、最成熟的频域成像算法。它的核心思想是将二维的匹配滤波分解为距离向和方位向的级联处理,并在距离-多普勒域通过插值高效校正距离徙动,从而实现高效精确成像。

    算法流程

    • 距离压缩

    雷达接收到的原始回波是目标散射系数与发射信号二维卷积的结果。距离压缩就是把回波信号在快时间(距离向)与发射信号的复共轭进行匹配滤波,将线性调频脉冲压缩成窄脉冲,实现距离向高分辨率。通常这一步骤在距离频域(或时域)完成,之后得到“距离压缩后、方位未压缩”的数据域。

    • 方位向FFT,进入距离-多普勒域

    对距离压缩后的数据沿方位向(慢时间)做快速傅里叶变换(FFT),变换到距离-多普勒域。 在这个域里:横轴是多普勒频率(方位频率),纵轴是距离向时间(或斜距) 此时,由于雷达与目标之间相对运动导致的距离变化,会使同一个点目标的轨迹在距离-多普勒域表现为一条距离徙动曲线。这就是需要校正的“距离徙动”。

    • 距离徙动校正(RCMC)

    距离徙动是SAR成像中特有的现象:雷达与目标的斜距随方位时间变化,导致同一目标的能量在距离压缩后散布在多个距离门上。在距离-多普勒域,同一方位频率处,所有目标的距离徙动量完全相同(只依赖于多普勒频率),这一特性被称为“距离徙动在距离-多普勒域的空不变性”。 距离徙动校正就是在这个域,根据多普勒频率计算出徙动量,然后对距离向数据进行插值平移,将能量搬回正确的距离门。常用的插值方法有sinc插值等。这一步骤是RDA的精髓,将二维耦合的问题解耦,使后续方位压缩可以独立进行。

    • 方位压缩

    经过距离徙动校正后,同一目标的能量已经对齐在一条直线上(同一个距离门),此时方位向等效为一个线性调频信号。方位压缩就是在距离-多普勒域,将方位信号与一个与多普勒频率相关的匹配滤波器相乘(频域匹配滤波),这个匹配滤波器的相位就是方位向调频率决定的二次相位。 相乘后,做方位向逆FFT(IFFT),将数据变回二维时域,就得到了聚焦后的SAR图像。图像中每个像素代表该位置散射系数的估计。

    RDA优点

    • 物理概念清晰,易于实现
    • 计算效率较高(利用FFT)
    • 正侧视及小斜视角下成像质量好

    RDA缺点

    • 大斜视时距离走动显著,需额外处理
    • RCMC 插值运算量较大
    • 假设方位向时域不变频域处理,对轨道非理想情况敏感

    Chirp Scaling 算法 (CSA)

    基本原理

    • 在距离-多普勒域,通过相位相乘(Chirp Scaling)使所有目标的距离徙动曲线形状一致
    • 避免复杂的时域插值,校正过程全通过相位相乘和FFT完成

    算法流程

  • 方位FFT,进入距离多普勒域
  • 乘以Chirp Scaling相位,均衡RCM
  • 距离FFT,在二维频域进行一致RCMC和距离压缩
  • 距离IFFT,接着进行方位压缩及残余相位补偿
  • 方位IFFT 输出图像
  • CSA 优点

    • 无需插值,运算速度更快
    • 对大斜视、宽测绘带适应性强
    • 相位保持精度较高,适合干涉处理

    CSA 局限

    • 需要多次FFT/IFFT,算法结构较复杂
    • 二次距离压缩(SRC)依赖多普勒频率,近似引入一定误差

    ωK 算法 (波数域算法)

    基本原理

    • 基于波方程推导,在二维频域进行精确的匹配滤波
    • 利用 Stolt 插值 完成距离频率轴的映射,实现非均匀距离徙动的统一聚焦

    算法流程

  • 二维FFT 进入波数域
  • 参考函数相乘(二维匹配滤波)
  • Stolt 插值(距离频率映射)
  • 二维IFFT 得到图像
  • ωK算法特点

    • 理论上是精确的线性成像算法(在理想运动下)
    • 可处理超宽测绘带和超大斜视角
    • Stolt 插值计算量较大,工程实现需优化

    后向投影算法 (BP)

    基本原理

    • 对成像网格每个像素,计算平台到该点的双程时延
    • 相干累加相应回波数据,实现像素重建

      I

      (

      x

      ,

      y

      )

      =

      s

      (

      τ

      =

      2

      R

      (

      η

      ;

      x

      ,

      y

      )

      c

      ,

      η

      )

      exp

      (

      j

      4

      π

      f

      c

      R

      (

      η

      ;

      x

      ,

      y

      )

      c

      )

      d

      η

      I(x,y) = \\int s(\\tau=\\frac{2R(\\eta; x,y)}{c}, \\eta) \\, \\exp\\left(j\\frac{4\\pi f_c R(\\eta; x,y)}{c}\\right) d\\eta

      I(x,y)=s(τ=c2R(η;x,y),η)exp(jc4πfcR(η;x,y))dη

    BP算法特点

    • 优点:任意轨迹、非线性孔径均适用,模型最通用,概念简单
    • 缺点:计算复杂度

      O

      (

      N

      3

      )

      O(N^3)

      O(N3),远高于频域算法;需加速技术(如快速BP)

    算法对比

    算法运算复杂度斜视/宽测绘带运动补偿灵活性典型应用
    RDA 低(FFT+插值) 小斜视角 一般 机载正侧视、星载条带
    CSA 中(无插值) 适应大斜视 一般 高斜视星载SAR
    ωK 较高(Stolt插值) 优秀 较弱 超高分辨率、聚束式
    BP 很高 任意 极灵活 机载复杂轨迹、无人机SAR

    总结

    • SAR成像基于匹配滤波与合成孔径原理,算法核心是解决距离徙动校正与方位聚焦
    • RDA 适用于常规正侧视,CSA 增强大斜视适应,ωK 提供理论精确解,BP 最灵活通用
    • 算法选择取决于分辨率要求、平台稳定性、计算资源和实时性需求
    • 未来算法将更智能、更高效,并深度融合运动补偿与优化技术

    仿真

    原始数据如下所示: 在这里插入图片描述 距离徙动校正后: 请添加图片描述 方位向压缩后: 请添加图片描述

    python代码(RDA算法)

    import numpy as np
    import matplotlib.pyplot as plt
    from numpy import pi, exp, sqrt, sinc, sin, cos, tan, arctan, ceil, real, imag, angle
    from numpy.fft import fft, fftshift, ifft
    from tqdm import trange

    plt.rcParams['font.sans-serif'] = ['SimHei']
    plt.rcParams['axes.unicode_minus'] = False

    # ==================== 参数设置 ====================
    # 距离向
    R_eta_c = 10000
    Tr = 2.5e-6
    Kr = 200e12
    alpha_os_r = 1.2
    Nrg = 320

    Bw = abs(Kr) * Tr
    Fr = alpha_os_r * Bw
    d_t_tau = 1 / Fr
    Trg = Nrg / Fr

    # 方位向
    c = 3e8
    Vr = 100
    f0 = 1e9
    Delta_f_dop = 50
    alpha_os_a = 1.25
    Naz = 256
    theta_r_c = 0.0

    lambda0 = c / f0
    t_eta_c = R_eta_c * sin(theta_r_c) / Vr
    f_eta_c = 2 * Vr * sin(theta_r_c) / lambda0
    La = 0.886 * 2 * Vr / Delta_f_dop
    Fa = alpha_os_a * Delta_f_dop
    R0 = R_eta_c * cos(theta_r_c)
    Ka_base = 2 * Vr ** 2 * cos(theta_r_c) ** 2 / (lambda0 * R0)
    theta_bw = 0.886 * lambda0 / La
    Taz = Naz / Fa
    d_t_eta = 1 / Fa

    rho_r = c / (2 * Fr)

    # ==================== 目标设置 ====================
    targets = np.array([[25, 50], [0, 0], [0, 100], [10, 150]], dtype=float)
    targets[:, 0] += R_eta_c

    Tar_t_eta_0 = targets[:, 1] / Vr
    Tar_t_eta_c = (targets[:, 1] targets[:, 0] * tan(theta_r_c)) / Vr

    # ==================== 网格 ====================
    t_tau = np.arange(Trg / 2, Trg / 2, d_t_tau) + 2 * R_eta_c / c
    t_eta = np.arange(Taz / 2, Taz / 2, d_t_eta) + t_eta_c
    r_tau = (t_tau * c / 2) * cos(theta_r_c)

    f_tau = fftshift(np.arange(Fr / 2, Fr / 2, Fr / Nrg))
    f_tau -= np.round(f_tau / Fr) * Fr
    f_eta = fftshift(np.arange(Fa / 2, Fa / 2, Fa / Naz))
    f_eta -= np.round((f_eta f_eta_c) / Fa) * Fa

    t_tauX, t_etaY = np.meshgrid(t_tau, t_eta)
    r_tauX, f_etaY = np.meshgrid(r_tau, f_eta)
    f_tau_X, f_eta_Y = np.meshgrid(f_tau, f_eta)

    # ==================== 生成原始回波 ====================
    st_tt = np.zeros((Naz, Nrg), dtype=complex)

    for i in trange(len(targets), desc="生成回波"):
    R_eta = np.sqrt(targets[i, 0] ** 2 + Vr ** 2 * (t_etaY Tar_t_eta_0[i]) ** 2)
    wr = np.abs(t_tauX 2 * R_eta / c) <= Tr / 2
    wa = sinc(0.886 * arctan(Vr * (t_etaY Tar_t_eta_c[i]) / targets[i, 0]) / theta_bw) ** 2
    phase = 1j * 4 * pi * f0 * R_eta / c + 1j * pi * Kr * (t_tauX 2 * R_eta / c) ** 2
    st_tt += wr * wa * np.exp(phase)

    # ==================== 距离压缩 ====================
    Sf_ft = fft(st_tt, Nrg, axis=1)
    window = fftshift(np.kaiser(Nrg, 2.5)).reshape(1, 1)
    Hrf = (np.abs(f_tau_X) <= Bw / 2) * window * np.exp(1j * pi * f_tau_X ** 2 / Kr)
    srt_tt = ifft(Sf_ft * Hrf, Nrg, axis=1)

    # ==================== 方位向FFT + 距离徙动校正 (RCMC) ====================
    Saf_tf = fft(srt_tt, Naz, axis=0)

    # RCMC (sinc插值)
    RCM = lambda0 ** 2 * r_tauX * f_etaY ** 2 / (8 * Vr ** 2)
    offset = (R0 + RCM R_eta_c) / rho_r

    Srcmf_tf = np.zeros_like(Saf_tf)
    for a in range(Naz):
    for r in range(Nrg):
    off = offset[a, r]
    off_ceil = int(ceil(off))
    frac = round((off_ceil off) * 16)
    idx = (r + off_ceil + np.arange(4, 4)) % Nrg

    if frac == 0:
    Srcmf_tf[a, r] = Saf_tf[a, idx[4]]
    else:
    # 简化的sinc核 (可进一步优化为scipy.interpolate)
    x = np.arange(4, 4) + (off_ceil off)
    hx = np.sinc(x) * np.kaiser(8, 2.5)
    hx /= hx.sum()
    Srcmf_tf[a, r] = np.dot(Saf_tf[a, idx], hx)

    # ==================== 方位压缩 ====================
    Ka = 2 * Vr ** 2 * cos(theta_r_c) ** 2 / (lambda0 * r_tauX)
    Haf = np.exp(1j * pi * f_etaY ** 2 / Ka) * np.exp(1j * 2 * pi * f_etaY * t_eta_c)
    soutt_tt = ifft(Srcmf_tf * Haf, Naz, axis=0)

    # ==================== 绘图 ====================
    def plot_raw():
    fig = plt.figure(figsize=(10, 7.5))
    for i, data in enumerate([real(st_tt), imag(st_tt), np.abs(st_tt), np.angle(st_tt)], 1):
    plt.subplot(2, 2, i)
    plt.pcolor(data, cmap='jet')
    plt.colorbar(shrink=1.0)
    plt.title(['(a)实部', '(b)虚部', '(c)幅度', '(d)相位'][i 1])
    plt.xlabel('距离时域(采样点)')
    plt.ylabel('方位时域(采样点)')
    plt.suptitle('图1,原始数据')

    def plot_rcm():
    fig = plt.figure(figsize=(10, 5))
    plt.subplot(1, 2, 1)
    plt.pcolor(np.abs(fftshift(Saf_tf, 0)), cmap='jet')
    plt.colorbar(shrink=1.0)
    plt.title('距离压缩')
    plt.subplot(1, 2, 2)
    plt.pcolor(np.abs(fftshift(Srcmf_tf, 0)), cmap='jet')
    plt.colorbar(shrink=1.0)
    plt.title('距离徙动校正后')

    def plot_final():
    fig = plt.figure(figsize=(10, 5))
    plt.subplot(121)
    plt.pcolor(real(soutt_tt), cmap='jet')
    plt.colorbar(shrink=1.0)
    plt.gca().invert_yaxis()
    plt.xlabel('距离向(采样点)→')
    plt.ylabel('←方位向(采样点)')
    plt.title('(a)实部')

    plt.subplot(122)
    plt.pcolor(np.abs(soutt_tt), cmap='jet')
    plt.colorbar(shrink=1.0)
    plt.gca().invert_yaxis()
    plt.xlabel('距离向(采样点)→')
    plt.ylabel('←方位向(采样点)')
    plt.title('(b)幅度')
    plt.suptitle('方位压缩后的仿真结果')

    plot_raw()
    plot_rcm()
    plot_final()
    plt.show()

    参考

    参考资料:

    • I. G. Cumming, F. H. Wong, Digital Processing of Synthetic Aperture Radar Data
    • J. C. Curlander, R. N. McDonough, Synthetic Aperture Radar: Systems and Signal Processing
    • 相关论文与开源代码库
    赞(0)
    未经允许不得转载:网硕互联帮助中心 » SAR成像原理和实现(含python代码)
    分享到: 更多 (0)

    评论 抢沙发

    评论前必须登录!