圆周扫描地基SAR地形表面提取辅助的三维成像算法

赵子恒 林 赟 王彦平 李 洋 申文杰

（北方工业大学信息学院雷达监测技术实验室，北京 100144）

1 引言

地基合成孔径雷达（Ground-based Synthetic Aperture Radar，GBSAR）作为一种新型高精度测量技术，具有非接触、全天候、设站灵活等优点，能够实时监测场景整体形变情况，已成为国内外各监测领域研究的热点［1-2］。传统GBSAR 系统通常采用直线轨道，该系统在滑坡监测预警和矿山桥梁形变监测等领域的技术适应性已经得到充分的应用检验。然而该系统只能获取二维图像和一维形变，且在地形起伏区域存在叠掩现象。为了解决上述问题，学者们提出了多基线直线轨道GBSAR 和圆周扫描GBSAR（Ground-based Circular Synthetic Aperture Radar，GBCSAR）两种具备三维分辨能力的GBSAR模式［3-5］，国内外研究机构经过大量实验工作后证明了GBCSAR 具备三维分辨能力，且可以有效解决二维成像中的叠掩现象。GBCSAR的合成孔径通过附在旋臂末端上的相位中心旋转生成，在保证系统分辨率的前提下，能够在一次旋转观测中覆盖监测区域的全部场景，实现三维成像。然而GBCSAR 独特的运动轨迹给成像算法的实现增加了难度，目前对于GBCSAR，学者们通常采用后向投影（Back Projection，BP）成像算法。该算法成像精度高，且可以应用于多种架构的雷达而不受天线运动轨迹所限制。文献［5］中提出了GBCSAR 模式下的三维层析BP算法，该算法通过叠加距离向上不同距离处的二维图像实现三维成像［6］。但是该成像算法未曾应用于大场景成像，且面对多目标点的场景时，层析BP算法数据量巨大，存在冗余现象，成像效率较低。在连续监测时数据量更大，无法满足实时形变监测的需求［7-8］。

因此本文提出了一种地形表面BP 成像算法，该算法首先通过层析BP 算法获取场景目标点的幅值三维矩阵，利用幅值矩阵提取出地形有效点的DEM 数据，而后在有效点三维曲面上实现BP 成像，不再借助参考平面投影成像［9］。在连续地形形变监测时，同一位置不同高度存在很多目标点，它们的形变趋势基本一致，因此可以将这些同一位置不同高度的多个目标点用一个有效点代替，这个有效点的成像数据可以有效地反映该位置不同高度处多个目标点的成像数据［10-12］。这样将场景中所有目标点的三维矩阵降维成一个由各个有效点接连的三维曲面，极大地减少了数据量，且三维曲面上的BP成像结果可以反映出监测场景的形变信息。

本文其余部分安排如下：第2部分介绍GBCSAR的几何模型和信号模型。第3部分对算法的原理进行论述，并具体介绍地形表面有效点曲面的获取方法以及如何在有效点曲面上实现BP 成像。第4 部分使用机载三维数据进行仿真实验，对算法的可行性进行验证，并与层析BP 算法进行计算量比较。最后，结论总结。

2 GBCSAR的几何模型和信号模型

GBCSAR 系统在旋臂的末端装有收发天线，该系统通过旋臂的旋转形成二维合成孔径。旋臂在YOZ 平面沿逆时针旋转，方位向和高度向都在YOZ平面内，距离向与旋转平面YOZ垂直。

GBCSAR 的几何模型如图1 所示。相位中心位于旋臂尾端a处，旋臂以坐标系原点O为中心在YOZ 平面旋转，监测区域位于与YOZ 平面垂直方向。在图1 中，a为相位中心位置，θ为相位中心的瞬时旋转角，a坐标随着θ变化。因此a坐标可以表示为(0，rcosθ，rsinθ)，其 中r为旋臂长度。p1，p2，…pn分别代表不同目标点，根据相位中心a坐标(xa，ya，za)和目标点pn坐标(xn，yn，zn)可以得到目标点与相位中心a之间的斜距Rn的表达式

为了方便计算，式（1）结合相位中心a坐标可得到瞬时斜距Rn表达式

已知了GBCSAR 的几何模型，接下来给出系统信号模型。这里以线性调频信号（LFM）作为发射信号，通常距离向脉冲压缩是在频域执行的，因此直接给出经过频域脉冲压缩后回波的频域表达式

其 中，σ是散射系数，Br是发射信号带宽，f∈[-Br/2，Br/2]，fc是中心频率，c是光速。在频域脉冲压缩后，需要对回波信号进行逆傅里叶变换，由于之后的BP 算法是在时域处理信号，因此这里给出距离脉冲压缩后信号的时域表达式

GBCSAR具备三维分辨能力，它的距离向分辨率为

其中Tr为雷达脉宽，又因为脉宽Tr与带宽Br互为倒数，因此（5）可表示为

由于GBCSAR 的方位向和高度向都在旋转平面YOZ 中，因此在Y方向和Z方向上，GBCSAR 的分辨率相同。这里推导方位向分辨率δy，定义波数K用来表示频率

定义比值αa来表示带宽与频率比，表达式如下

由此可以得到GBCSAR 在Y方向和Z方向上的分辨率

3 GBCSAR地形表面成像算法

对于处在同一(X，Y)坐标上的多个目标点，所在高度Z不同，导致各目标点与雷达之间的斜距Rn和幅值大小f(X，Y，Z) 都不相同。其中幅值f(X，Y，Z)最大的目标点Pn通常在对应(X，Y)坐标上多个目标点中散射系数最大，这个点的形变信息可以基本反映出该(X，Y)坐标上所有目标点的形变趋势。因此将同一(X，Y)坐标上多个目标点降维成一个目标点，即该(X，Y)坐标上的有效点，对这个有效点进行相应的成像就可以有效地实现这个位置的形变监测。而将整个场景所有(X，Y)坐标对应的有效点提取出来后，在有效点曲面上进行BP 成像，即可实现对整片场景的形变监测。本文提出的方法两个主要步骤分别是获取地形有效点所构成的三维曲面和在该曲面上实现BP成像。

本节介绍了如何通过目标点三维坐标获取地形有效点DEM数据和有效点三维曲面，以及如何在有效点三维曲面上BP成像。

3.1 地形表面有效点曲面提取

在采用层析BP 算法获取所有目标点的幅值前需要对来自场景的SAR 回波数据进行预处理，对接收到的实信号进行正交解调得到算法需要的复信号。而后采用层析BP 算法获取所有目标点的幅值，得到目标点幅值矩阵。在连续监测时，只需在开始时使用一次层析BP 算法提取有效点曲面，之后的连续监测都在提取得到的有效点曲面上成像。

提取有效点时需要根据场景大小和目标点的分布对幅值矩阵进行插值或等间隔采样，以确保目标点密集分布时不会产生数据冗余，同时不会出现目标点稀疏分布时提取到的有效点不连续，并且要确保有效点的高度差大于系统理论分辨率δz，否则会出现混叠现象。而后按照目标点的(X，Y)坐标遍历查找对应(X，Y)坐标上的有效点，此时得到的有效点坐标为(X，Y，Z(x，y))，即有效点的高度Z是根据对应(X，Y)坐标一一对应的。当所有(X，Y)遍历结束后，整个场景的有效点会连接成一个三维曲面，这个曲面就是最终BP成像所在的投影曲面。

如图3 所示，右侧的三维曲面就是通过提取地形有效点得到的，需要注意的是这个曲面并不是场景的顶面，而是根据有效点的分布位置连接而成的。相比于左侧的三维立方场景，三维曲面极大地减少了目标点个数，减少了计算量。

3.2 有效点三维曲面BP成像

与传统的二维BP算法和层析三维BP算法不同的是，传统二维BP算法设定像素点都在同一高度平面上，通常选择地平面。将整个场景中某一个高度的回波数据都投影在一个平面，因此传统二维BP算法只能实现单一高度的成像，无法获取整个场景信息。而层析三维BP 算法是将各个高度的回波数据分别投影在各个平面，需要在三维空间中每个高度上都设定等量的像素点，因此计算量庞大，无法满足连续监测的需求。层析BP算法示意图如图4所示。

本文提出的地形表面有效点成像则是依据有效点的高度在对应高度设定像素点，没有有效点的高度也就没有像素点。最终像素点会连接成有起伏的三维曲面，并在该三维曲面上进行投影成像。地形表面有效点成像示意图如下图5所示。

三维曲面上的BP 成像首先需要借助距离脉冲压缩后的回波获得目标点幅值矩阵，借助幅值矩阵获取有效点坐标，然后根据有效点坐标和DEM数据构建成像网格。计算成像网格中每个像素点距离相位中心a的瞬时斜距，得到每个像素点的双程延时以提取对应像素点的幅值，最后叠加各个方位向在各像素点的幅值得到SAR 图像。具体实现步骤如下：

步骤1：根据系统理论分辨率和提取到的三维曲面划分成像网格；

步骤2：计算各个像素点距相位中心a当前位置的双程延时；

步骤3：通过插值从雷达回波中依据双程延时提取当前方位向上像素点的回波值；

步骤4：对得到的回波值进行相位补偿得到对应方位向上的子图像；

步骤5：遍历各个方位向，重复步骤3、4，叠加各个方位向的子图像；

步骤6：得到最终的SAR图像。

根据式（6）和（9）得到GBCSAR 的理论分辨率，X和Y方向设定的网格间隔应略小于δx和δy，在提取有效点时已经保证了有效点在Z方向的间隔大于理论分辨率δz。X和Y方向设定的成像网格个数根据场景大小设定，最终X和Y方向的成像网格范围包含场景所有目标点即可，成像网格在Z方向的分布和有效点在Z方向的分布一致。设定好的像素点坐标可以表示为(xi，yj，z(xi，yj))，其中i，j表示XOY 平面中第i行第j列的像素点，像素点所在的高度z(xi，y j)是从有效点DEM 中遍历出的，最终像素点都是分布在提取出的有效点附近。

根据相位中心a的坐标(0，ya，za)和像素点的坐标(xi，yj，z(xi，yj))计算出各个像素点到相位中心的距离

像素点到相位中心的双程延时表示为：

由于回波数据是离散地存储在二维矩阵中，回波值在矩阵中对应的都是整数索引，而双程延时在回波矩阵中对应的回波值索引大多数不是整数，想要得到离散信号非整数点的值就需要对信号进行插值。因此先对回波数据插值，再利用延时和回波的对应关系从回波数据中取得对应方位向的回波值，本文采用的是sinc插值，插值核为8，插值公式为：

插值后得到各像素点的回波值需要进行相位补偿，以实现各方位向的聚焦。将回波值与对应像素点的相位补偿因子相乘即可实现补偿：

其中R(i，j)为当前像素点与相位中心的距离。最终得到像素点P处的幅值

3.3 计算量分析

层析BP 算法逐个平面进行插值和相位补偿，成像结果准确且立体，但在大场景下高度向Z的采样点数可能很大，因此需要操作的点数庞大，计算效率较低。本文提出的地形表面有效点成像实现了降维成像，减少了像素点个数和插值及补偿的次数，能够实现更高效的成像。假设算法中一个像素点的插值和补偿操作耗时t1，一次有效点DEM 提取操作耗时t2。

（1）层析BP 算法：设定X，Y，Z方向像素点个数分别为Nx，Ny，Nz，算法耗时为

（2）地形有效点BP 算法：设定提取的有效点共有M个不同的高度，那么像素点也是分布在不同高度的，在M个高度共分布Nx×Ny个像素点，算法耗时为

可以得到层析BP算法和地形有效点BP算法的耗时比为

由于DEM 提取耗时t2非常少，即t2＜＜Nx×Ny×t1，因此耗时比可以近似为层析BP时成像区域高度向的网格点数Nz，在大场景下Nz通常较大。因此可以得出结论，在同一场景且成像网格间隔设定相同的条件下，地形表面有效点BP 算法效率提升明显。

4 验证实验

4.1 成像结果分析

本节我们利用机载三维图像数据进行仿真实验，主要的实验参数如表1所示。

表1 仿真参数Tab.1 Simulation parameters

仿真实验采用的原始数据是929*938*141的三维矩阵，具有目标点的三维坐标以及散射系数。由于数据量太大不方便仿真实验，因此对数据进行等间隔采样，采样后得到92*93*141 的三维矩阵。首先给出三维数据在中间高度的二维图像。

利用层析BP算法获取到目标点的幅值矩阵后，提取出有效点的DEM数据，并根据有效点的在矩阵中的索引从散射系数矩阵中提取出对应的散射系数值。以下是有效点DEM图像和有效点RCS图像。

可以从图6和图7（a）对比发现原始图像和提取出的有效点DEM图像基本一致，有效点的分布反映了整个场景目标点的分布。

在三维数据图像中进行了高度向成像切片分析，选择了目标点高度变化比较明显的位置Y=40 m 处，下面是Y=40 m 处目标点的幅值与高度的关系以及Y=40 m处有效点的高度。

从图8（a）可以看出在X=0 m 到35 m 处最大幅值交替分布在60 m 到70 m 高度，X=35 m 之后最大幅值稳定在70 m 高度处。这和图8（b）中Y=40 m 处的有效点DEM 分布相同，可以证明提取得到的目标点DEM数据准确。

接下来利用地形有效点BP算法在提取出的有效点曲面上进行成像，并将有效点的散射系数附加在对应有效点上生成回波。得到BP成像结果如下。

通过图9可以发现当回波附加上对应的散射系数后，有效点曲面上的BP 成像结果和原始三维数据图像一致，可以看出场景中目标点的分布。在连续形变监测时，根据前后两张成像结果中像素点的幅值变化就可以分析出场景的形变信息。

图10中可以看到有一个相对独立的强散射点，该点周围的散射点幅值较低，不会影响该点的旁瓣值，因此针对这个点进行成像效果分析。

如图11 所示，将实际场景成像图按照5∶1 放大，可以得到场景中位于坐标（521，527）处的独立点目标周围区域展示图。再对该点的成像结果进行一维成像分析，得到对应的峰值旁瓣比和-3 dB处宽度，与理论值进行比较。

图12是该点方位向的剖面图，可以从图中分析方位向成像指标。

经过分析后得到实际峰值旁瓣比为-8.08，-3 dB处宽度为27 m，与理论值基本一致。

如图13所示，独立点目标在距离向上的一维成像结果是sinc 函数，符合圆周扫描地基SAR 的成像模型原理。

综上，本文中对于原始三维机载图像数据的成像仿真实验达到了预期的效果，且方位向和距离向的成像指标都符合成像模型原理。

5 结论

本文提出了一种针对GBCSAR 的地形表面成像算法，通过将场景三维目标点提取为三维曲面，并在该曲面上投影成像，实现了地形表面成像。本文利用BP 算法的优势避免了由于GBCSAR 弧形孔径带来的轨迹误差，并且在曲面投影时省去了平面坐标转换等复杂操作。此外，与传统的层析后向投影算法相比，地形表面后向投影算法通过降维的方式减小了计算量，提高了成像效率，这使得GBCSAR在地形形变连续监测领域具有优势。