经典文献翻译：合成孔径雷达快速反投影算法(Fast Backprojection Algorithm for Synthetic Aperture Radar)

排版不方便，图片不好上传，Word版翻译及原文下载地址：//download.csdn.net/download/u010112513/12576783

Ali F. Yegulalp

翻译：Ivan

摘要：介绍了一种新的合成孔径雷达（SAR）时域反投影算法。该算法将重现由标准的逐像素反投影算法生成的N×N像素图像，其运行速度大约快√N  倍。快速反投影保留了标准反投影的优点：任何飞行路径的完美运动补偿、低副瓣水平、无限场景大小、任意宽频带和积分角的完美聚焦，以及严格的本地处理（即，在采集脉冲时可以进行处理，而无需沿轨迹缓冲或拐弯）。新算法还可以将正在进行中的图像存储在磁盘上（而不是内存中），只需稍微降低处理速度。

1 简介

合成孔径雷达成像算法有两大类：基于FFTs的方法和反投影（时域）方法。这两个类都有各自的优点，并且都在使用中。

基于FFT的方法倾向于最小化处理负载，但并不普遍适用。距离多普勒和极坐标格式算法依赖于在各种限制（大图像尺寸、带宽或积分角等）下分解的几何近似。更复杂的w-k算法在几何上是精确的，但它仍然假设飞机飞行的轨迹是完全直线和均匀的。偏离轨道的运动误差只能以近似的方式进行补偿，这会导致在积分角较大和/或运动误差较大时聚焦不良。飞机沿轨道速度的变化要求数据重新采样到均匀的间距。对于大多数FFT方法来说，另一个问题是它们需要在频域内对数据进行插值；插值误差会导致伪影扩散到整个图像上。chirp scaling算法不需要频域插值，但它依赖于Stolt变换的近似，该变换对于低频下的高分辨率图像会失效。困难也可能与算法实现有关，特别是对于实时处理场景。为了计算沿航迹的fft，在开始处理之前必须获取和缓冲全孔径的脉冲。二维fft需要转角操作，这可能成为并行计算体系结构的处理瓶颈。

先前的每一个问题都可以通过使用反投影来消除，代价是大大增加了处理时间。本文描述了一种新的算法，称为“快速反投影算法”，它在运行速度更快的同时，模仿了标准反投影的结果和优点。该算法是在美国国防部高级研究计划局赞助的麻省理工学院林肯实验室（MIT Lincoln Laboratory）FOPEN项目中开发、实现和使用的。在处理CARABAS 2 VHF合成孔径雷达数据时，我们通常获得大约20倍的速度提高。

论文的其余部分格式如下：我们首先回顾标准的反投影。回顾也有助于建立我们的基本假设和符号。我们将描述标准反投影的合理实现，并分析运算量。然后论文转向快速反投影。在简要介绍了快速反投影的工作原理后，给出了基于图像的基本采样原理的简单演示。在这个演示之后，将对采样要求进行严格分析，并描述如何实现该算法。分析了快速反投影的运算量，并与标准反投影进行了比较。最后，利用仿真数据对标准反投影和快速反投影进行了并行比较。

2 标准BP算法

1. 数据模型

要进行反投影，我们需要雷达运动历史和一组距离压缩脉冲。用矢量  表示雷达位置，参数  表示路径长度。在位置  处，接收距离压缩数据  。我们做了一个标准的近似，即当脉冲被发送和接收时，雷达不会移动很长的距离。

例如，考虑理想点目标位于位置  的情况。假设雷达工作在低频vmin和高频vmax之间，数据在谱域具有平坦的幅度和线性相位。在时间域中，数据表示为

 

 

1. 反投影公式

反投影的第一步是对时域数据  应用斜坡滤波器。在谱域中，滤波操作可以表示为

 

位置  处的图像像素的值是通过对滤波的时域数据进行积分来给出的：

 

实际上，我们的数据是离散采样的，因此我们必须用和来代替积分：

 

注意，离散情况中的孔径点  不需要以规则的方式采样。只要孔径采样率高于奈奎斯特速率，离散型就近似于积分型。连续和离散之间的对应关系是相当明显的，因此我们在后面简单地使用连续体表示法。

反投影方程适用于任意形状的轨迹，原理上可得到三维图像。对于强非线性轨迹，例如圆，实际上可以在所有三个空间方向上获得合理的分辨率。对于直线轨迹，三维图像函数  沿环绕飞行轨迹的圆完全恒定。在这种情况下，通常将图像点  限制为单个水平面，我们称之为“图像平面”。如果已知地面地形，也可以计算真实地面上的  。虽然直线情况是最容易讨论的，但本文提出的方法基本上可以应用于任何类型的飞行轨迹。

1. 积分孔径

在应用方程3时，必须对积分孔径作一定的选择（即s的积分极限）。最简单的选择是对每个图像像素  使用整个可用的孔径，通常称为“聚束模式”处理。虽然点响应函数（和分辨率）在空间上不一致，但这种方法可以得到尽可能高的分辨率的图像。

另一个常见的选项是条带模式处理，它对每个像素使用固定的积分角。这种情况下的角度是飞机航向与像素视线之间的角度。如果我们想从  积分到  ，方程3中s的积分极限将由以下关系确定：

 

注意，导数  只是与轨迹  相切的单位向量。  是一个单位向量，这一事实遵循了我们对s的定义，即沿轨迹的路径长度。

由于spotlight与stripmap的问题对本文给出的任何结果都不是关键，因此我们将把讨论限制在spotlight模式的情况下，并将所有内容都以简单的方式概括为stripmap模式的情况。

1. 实现

在讨论反投影时，我们一直在考虑以逐像素的方式进行计算。换言之，我们选择一个像素，在孔径上积分以计算其值，然后移动到下一个像素。在实际中，逐脉工作更有意义。对于每个孔径位置s，我们有一个脉冲  ，它对我们所需图像中的每个像素都有贡献。换句话说，我们是并行计算像素值，而不是串行计算。逐个脉冲工作的优点是我们只需要访问一次每个脉冲，而且在插值方面我们也会受益（如下所述）。

公式2中的滤波步骤可以使用时域FIR滤波器或简单地使用FFT来执行。 等式3中的积分要求我们对t执行  的插值。这样做的合理方法是相对于t上采样  。然后可以对  的上采样版本使用线性或近邻内插法评估积分。上采样步骤可以跳过，但是线性插值或最近邻插值会在采样率接近Nyquist的情况下产生较差的结果。此外，上采样步骤实际上是一个较小的计算负担，因为它是在一维范围数据向量上执行的，而积分步骤是针对二维图像中的每个像素执行的。经验表明，可以升采样至少6倍，而速度几乎没有损失。悲观的观点是，反向投影速度太慢，以至于很难注意到增加采样的额外负担！

最后，应该注意的是，在上采样和插值之前，必须从  中移除载波相位，然后在执行反投影和时恢复它。

1. 计算负荷

想象一下使用标准反投影对Npulse压缩数据形成一个N x N像素的图像。 如果我们使用聚束模式处理，则等式3中的积分将遍及整个可用孔径。 运算数为

 

我们忽略了滤波和上采样所需的运算，因为它们只是对等式6的次要修正。

3 快速BP算法

快速反投影通过将整个合成孔径划分为子孔径而起工作。每个子孔径生成一个图像，并通过对子孔径图像进行相干求和来获得最终图像。关键是子孔径图像在方位向的分辨率非常低；这意味着我们可以在方位向以较粗的像素网格上计算它们。当需要将子孔径图像组合成最终的高分辨率图像时，我们会在方位向内将其上采样到全分辨率。由于子孔径图像在方位向内包含的像素少于最终图像，因此与标准反投影相比，我们需要执行的操作要少得多。

在以下各节中，我们解释了对子孔径使用局部极坐标的必要性，导出采样要求，描述算法的合理实现方式以及估算总操作次数。

1. 子孔径处理

假设合成孔径的总长度为L，那么反投影公式为

 

可以把全孔径分割为Nsub个子孔径  。图像  变为子孔径图像之和 

 

子孔径图像由下式给出

 

第n个子孔径的中心由sn=(n-1/2)l给出。

在公式8中没有近似值。我们只是将公式7中的积分分解为一系列子积分的和。 从这里开始，我们将专注于有效评估单个子孔径图像。

1. 笛卡尔坐标系VS极坐标系

图1是使用子孔径创建的反投影图像。 子孔径的中心位于图像顶部上方的某个位置。 请注意，我们在距离向上具有高分辨率，但在方位向内却具有非常低的分辨率。图像区域有一个很大的角度，所以恒定范围的线很明显是圆弧，而不是直线。图2显示了图1中图像的二维FFT。FFT数据填充了大部分窗口，这表明图像是在Nyquist附近采样的，尽管在方位向分辨率很低。如果我们在图1中使用更粗的像素间距，将导致混叠。

 

图1 子孔径成像

 

图2 图1中图像的二维FFT

要利用低方位向的分辨率，我们必须处理分布在极坐标网格上的像素。换句话说，我们应该定位我们的像素，以便它们在半径和角度上以均匀的增量间隔（从子孔径的中心开始测量）。图3显示了与图1相同的数据，但是像素位于极坐标网格上。垂直轴对应距离向坐标，水平轴对应角坐标。注意，在图3中，恒定范围的直线是水平直线。图4显示了极坐标网格图像的二维FFT。频域数据在水平方向上只占很窄的区域，这意味着图像数据沿水平轴高度过采样。快速反投影利用了极坐标的低带宽优势，它在角方向上的步数非常粗。像素数量的急剧减少导致了比标准反投影速度的显著提高。

 

图3 图1中极坐标网格上的图像

 

图4 图3中极坐标网格图像的二维FFT

我们提醒读者不要将本文讨论的极坐标网格与众所周知的SAR图像形成的极坐标算法混淆。快速反投影中使用的极坐标网格是图像域中的极坐标网格，而极坐标算法使用的是频域中的极坐标网格。

1. 极坐标系下的奈奎斯特采样率

在上一节中，我们使用了一个简单的可视化示例来说明使用极坐标的带宽优势。本节分析反投影方程，导出极坐标网格图像Nyquist速率的简单表达式。这些表达式精确地告诉我们极坐标网格子孔径图像中像素的间距。

第n个子孔径的图像由方程9给出。在本节中，我们假设子孔径上的航迹是一条直线；对于足够短的子孔径，这显然是一个合理的假设。此外，我们做这个假设只是为了计算奈奎斯特采样率；在实际的图像形成过程中，我们使用真实的飞行轨迹。真实飞行轨迹在子孔径期间不完全是直线的这一事实意味着奈奎斯特速率可能略高于直线子孔径。如果我们不太接近奈奎斯特的采样率，应该没有问题。

我们的第一步是根据极坐标(r,   )重新表达图像公式。坐标r定义为距子孔径中心的距离：

 

坐标  被定义为飞行轨迹与像素视线之间的角度余弦：

 

对于直线航迹，我们也有有用的标识

 

图像公式中的关键量是图像像素与孔径上的点之间的距离。 我们可以在极坐标中重新表达它，如下所示：

 

极坐标中的子孔径图像由下式给出

 

下一步是将  傅里叶展开

 

我们分析的最后一步是相对于r和  计算  的傅立叶变换：

 

使用固定相位方法可以非常准确地计算出该变换。 固定相条件为：

 

求解上述方程式可知，  为非零时

 

在写最后两个方程式时，我们做了一个额外的近似，即子孔径长度远小于斜距。这个近似值不需要使用，但它适用于几乎所有感兴趣的情况，而且它还使结果特别明显。

  和  的界限可直接转换为r和  的奈奎斯特采样要求：

 

  中所需的样本数量将明显与子孔径长度l成比例。另一方面，r中的采样率仅取决于系统带宽，而不取决于子孔径的大小。

1. 实现

到目前为止，我们主要关注快速反投影背后的数学推理。我们现在将描述算法的合理实现。和标准的反投影一样，我们将逐个脉冲地处理数据：我们读取一个脉冲，计算它对图像的贡献，然后移动到下一个脉冲。方程2中的脉冲滤波可以用与标准反投影相同的方式进行：使用时域FIR滤波器或简单地使用FFTs。

第一步是为所需的输出图像选择图像网格。这仅仅意味着我们必须知道图像中每个像素的三维坐标，并且我们必须分配足够的存储空间（内存或磁盘空间）来容纳每个像素的复数值。如果将像素排列在规则的笛卡尔网格中，则容易指定角像素的位置以及距离向和方位向上的步长的位移矢量。如果像素位于不规则的位置（例如，在遵循真实局部地形的网格上），则需要存储每个像素的三维位置。

接下来，我们必须选择子孔径的大小。 如下一节所示，对于N * N像素的图像，最佳尺寸大约是每个子孔径√N  个脉冲数。 在实践中，通常可以通过少量的反复试验来选择合理的子孔径。

该处理逐子孔径进行。 在每个子孔径的开始处，我们计算以当前子孔径的中间为中心的极坐标网格的像素位置。径向和角度方向上的像素间距被选择为略高于等式11和12中的奈奎斯特速率。如前所述，像素间距在角度方向上非常稀疏。

子孔径的脉冲处理与标准反投影完全相同，只是像素位于我们特殊的极坐标网格上。特别是，我们使用子孔径上的真实运动历史，而不是直线近似。

一旦我们处理完当前子孔径中的脉冲，我们必须将结果从极坐标网格传输回最终输出图像网格，这需要几个步骤。首先，我们从极坐标图像数据中去除载波相位。接下来，我们对两个方向的极坐标图像数据进行上采样。然后，我们通过在上采样的极坐标图像数据中找到对应的点来计算输出图像中每个像素的贡献。二维线性插值用于近似估计上采样极坐标网格上像素之间的极坐标图像数据。线性插值后，载波相位恢复。从极坐标网格转换为输出网格的精度完全取决于上采样因子。 我们上采样越多，结果越准确。实际上，上采样系数在2到8之间可提供速度和精度的合理组合。

将极坐标图像数据传输到输出图像网格后，我们移至下一个子孔径并计算新的极坐标网格。 我们将继续处理子孔径，直到使用完所有脉冲。由于每个子孔径仅更新输出图像数据一次，因此将输出图像存储在磁盘上而不是存储在内存中并非没有道理。此过程允许以相对较少的内存创建非常大的图像，而速度仅会受到轻微影响。极坐标栅格图像必须存储在内存中，但是它们需要的空间要少得多（即使在上采样之后），因为它们在角度方向上很稀疏。

1. 计算负荷和最佳子孔径

略

4 FBP与BP比较

为了补充上一节中的理论分析，我们现在提供标准和快速反投影的量化并行比较。为了验证电目标脉冲响应，我们使用合成数据模拟理想的点散射体（即频率和孔径的均匀响应）。数据在频率和孔径上都有一个Kaiser窗口。 如果没有窗口，则标准和快速反投影图像之间的差异将隐藏在脉冲响应旁瓣中。

基本仿真参数如下：

带宽：20~90MHz

入射角：-45~45度

像素大小：距离1.5m，方位1.0m

图像尺寸：2000*2000像素

图像中心：4000m斜距

标准反投影软件对这些数据运行了两次。 在第一次中，我们对距离脉冲进行了6倍的上采样； 此运行用于确定标准反向注入的运行速度。在第二轮中，将距离脉冲上采样12倍； 来自这次运行的图像被用作参考图像，以测试快速反投影图像的质量。快速反投影软件也运行了多次，使用1到8之间的上采样因子。选择不同子孔径大小以在每种情况下优化处理速度。

 

图5 标准反投影，上采样系数=为12

 

图6 快速反投影，上采样因子= 1

 

图7 快速反投影，上采样系数= 8

 

图8 快速反投影残差，上采样系数= 8

1. 成像质量

图5显示了使用标准反投影创建的图像的200×200像素部分。图像被标准化，使得峰值像素为0db。图像中显示宽动态范围（80分贝），以强调远离主峰的处理副瓣。在频率和孔径上使用Kaiser窗口来抑制旁瓣。

图6和图7显示了使用快速反投影创建的图像，上采样因子为1和8。再次，我们在图像中使用宽动态范围，以使处理伪影明显。注意当上采样因子从1增加到8时的改进。

为了数值评估快速反投影图像的质量，我们将每个图像与图5中的标准反投影图像进行了比较。通过从标准反投影图像中相干地减去每个快速反投影图像，进行比较。图8显示了上采样因子为8的图像减去后剩余像素功率的散点图。水平轴表示标准反投影图像中的像素功率，垂直轴表示相干相减后对应像素的功率。对角线表示误差与正确像素值一样大的边界。

1. 处理速度

测量给定反投影算法速度的一种合理方法是以每秒脉冲像素为单位。换句话说，将N个脉冲处理成M像素图像需要多长时间？对于标准反投影，答案仅取决于乘积NM。因此，我们以每秒数百万个脉冲像素（Mpp/s）的速度来表示我们的速度测量。

标准反投影以0.71mpp/s的速度运行。对于快速反投影，处理速度取决于上采样因子。下表显示了处理速度（以Mpp/s为单位）、相对于标准反投影的速度以及最大相干残差（dB）：

 

最大残差总是出现在主叶峰值。例如，在上采样因子为4的情况下，最大剩余功率为-37.5dB，这转化为主瓣峰值处像素复值的1.3%的误差。远离主叶的残差总是小得多。

5 结论

我们提出了一种快速反投影算法，该算法可以复制标准反投影的结果，但是在实际情况下的运行速度要快10到30倍。通过改变上采样因子，用户可以根据特定应用的需求设定速度和准确性。该算法可以应用于来自任何一种SAR的数据，并且效果良好，特别是对于宽带、大观测角或具有长相干孔径且具有明显偏离轨迹运动的传感器特别有吸引力。这也是一种处理在非线性轨迹上收集的数据的方法。