基于同心方形网格插值处理的柱面SAR成像算法

丁 丽 何华港 王 韬 储得苗

①(上海理工大学光电信息与计算机工程学院 上海 200093)

②(上海理工大学生物医学工程研究所 上海 200093)

1 引言

柱面毫米波合成孔径雷达(Cylindrical millimeter-wave Synthetic Aperture Radar, CSAR)采用圆迹扫描实现了目标场景360°无盲区检测[1]，相比于传统直线轨迹SAR, CSAR具有高空间覆盖和空间分辨率、高信噪比等优势，在安检成像、目标检测等领域得到了广泛的研究[2-6]。目前，CSAR成像算法主要基于傅里叶变换理论，需要通过2维插值解决柱面扫描引入的方位维和距离维向波数域非均匀填充的问题，但是这两个维度的强耦合性，采用2维逐点遍历插值的算法复杂度高，导致整个成像算法实时性差。

文献[7]首次提出基于波前重建的CSAR成像算法，为后续柱面成像算法的发展奠定了理论基础。文献[8]提出的CSAR成像算法利用球面波分解理论，通过傅里叶变换、相位补偿、2维插值等频域处理方式得到了目标3维图像。文献[9]提出了改进的3维CSAR成像算法，通过对高度维进行运动补偿，消除了阵列天线在高度维运动造成的运动误差，提高了成像性能。文献[10]提出了一种降维时域相关法用于CSAR，在距离域和方位域采用相干累加，使3维时域相关降维至2维时域相关，但是该算法需要对空间目标逐点进行相干累加成像，成像效率低，无法满足实时成像的需求。文献[11]在聚束SAR成像处理中利用Chirp-Z变换(Chirp-Z Transform, CZT)实现方位维和距离维重采样替代2维插值，提高了算法效率。但是CSAR由于采用圆迹扫描使得径向波数域是非均匀采样，而CZT变换具有从均匀数据重采样到均匀数据点的特点，使得无法使用CZT处理CSAR中的2维插值。因此，对于CSAR中波数域方位维与距离维的2维非均匀性和强耦合性造成成像实现效率低、难以开展并行处理的问题仍需进一步研究。

为解决方位维和距离维的强耦合性，本文提出一种基于同心方形网格处理的插值方法。基于柱面成像中波数域数据分布特点的分析，首先，在每一个高度维下对方位维和距离维波数域数据在极坐标系下进行径向扩充补零，得到膨胀的虚拟同心圆环带，实现波数域频带利用率最大化。接着，将补0后分布在极坐标系下的数据沿着径向1维插值到同心方形网格中。在进行径向1维插值之后，原始波数域方位维和距离维的2维强耦合已经可以在直角坐标系下进行分解，因此沿着波数域方位维×距离维横截面方形的对角线进行区间划分，把数据划分为两个正交不交叠的基本垂直点数据区间和基本水平点数据区间，至此实现对方位维和距离维的解耦。再分别进行方位维和距离维的独立1维插值，使数据重采样到均匀网格中，避免了传统2维插值在遍历2维网格时所耗的时间复杂度和在待插值点查找局部邻域所需的计算复杂度等难题，大大减小插值的整体运算量。最后，将两部分数据进行拼接得到直角坐标系下均匀分布的波数域数据，利用快速傅里叶逆变换重构目标图像。实验证明，改进算法插值处理速度较传统算法提高了7倍，符合算法复杂度理论推导，证明所提算法能大幅度提高成像效率。

2 柱面SAR成像算法

图1(a)为柱面3维成像系统模型。以场景中心O为原点建立直角坐标系，X轴、Y轴和Z轴分别表示为方位维、距离维和高度维。天线阵列的高度为Lz，以R为半径进行圆周扫描，θ表示某时刻天线旋转方位角，且θ ∈[π,-π]。设天线阵列任意采样的位置坐标为(x′,y′,z′)=(Rcosθ,Rsinθ,z′)，某一散射点目标的位置坐标为(x,y,z)=(rcosφ,rsinφ,z)，对应的后向散射系数为σ(x,y,z)。

图1 柱面3维成像系统模型

在实际系统中，CSAR回波信号是体目标在3维空间内散射系数的体积分，信号的幅度随着距离的衰减忽略不计，因为它对图像的聚焦产生很小的影响，所以在(θ,w,z)域中接收到目标的回波信号为[12]

其中，k为波数，令Φ为相位

对式(1)两边做高度维和角度维的傅里叶变换得

其中，kθ和kz分别为波数域角度维分量和高度维分量。根据驻定相位法，对相位Φ沿z′方向求导取0点得

根据式(4)，令ε为天线到目标视线的俯仰角，则可以得到

因此利用正切与正弦三角函数关系解得驻相点z′

将式( 6 ) 代入式( 3 ) 中得Φ=-kr同理可求θ的驻相点，对相位Φ沿θ方向进行求导并令其为0可得

根据式(7)，将CSAR沿着高度维投影到XOY平面上，其投影后的空间几何关系如图1(b)所示，其中，天线到目标直线分别与天线到原点直线之间的夹角，和与目标到原点直线之间的夹角记为α和β，且有θ=π-α-β+φ。利用图1(b)中所示三角形的几何关系可得

根据三角函数变换求得驻相点θ

将式(9)代入到Φ中得

令Φ=Φ1+Φ2，其中

将式(11)代入式(2)中得

式(11)表明相位量Φ1与积分变量无关，故对式(12)两边乘以exp(-jΦ1)进行相位补偿得

对式(13)两边进行kθ维傅里叶逆变换得

令式(14)中被积函数的相位为

根据驻定相位法，将式(15)在kθ维进行求导且令其为0，得到kθ=rkrsin(θ-φ)，代入式(14)

其中，kx是方位向波数，ky是距离维波数

式(17)表明波数域数据在直角坐标系下的方位维kx和距离维ky是非均匀分布的[14]，将直角坐标系下非均匀分布的S6(θ,w,kz)通过方位维和距离维的插值，得到均匀分布的频率数据

其中，interpkx,ky{·}表示插值函数，和分别表示插值后均匀采样的波数域方位维分量和距离维分量，对式(18)进行3维傅里叶逆变换可得到目标重建图像

3 同心方形网格插值方法

根据式(16)，给定任一高度维波数，CSAR波数域是一个在极坐标系下表现出为kr非均匀分布、kθ为均匀分布的圆环。通过插值从非均匀极坐标系获得直角坐标系下均匀分布的波数域值，获得式(18)所需的呈现同心方形环带的波数域填充样式，如图2所示。尽管如此，式(17)表明波数域数据在方位维和距离维满足存在强耦合性，传统插值需要通过2维逐点遍历实现，即每一个待插值点在2维空间寻找邻域点，运算复杂度高。因此，本文利用同心方形网格的分布对称性与CSAR波数域在极坐标下的采样特点，提出一种波数域数据的2维解耦方法，将2维逐点遍历插值分解为两个独立的1维数据插值处理，大大减小插值所耗的运算量，完成2维插值的高效实现。改进后的插值步骤如图3所示，分为径向kr数据补0、径向插值、平面数据分区、波数域方位向与距离向独立插值、数据拼接5个步骤来实现。

图2 柱面SAR成像数据分布图

步骤1 径向补0如图3(a)-图3(c)所示，得到膨胀的虚拟同心圆环带。

所需的同心方形带填充样式如图3(a)中实线框定的区域所示，能完全覆盖原始带宽引入的同心圆环波数域填充区域(实心环带所覆盖的区域)，这里对S5(θ,ω,kz)频谱数据在径向进行最内环和最外环的补0膨胀，找到能够覆盖同心方形区域的虚拟同心圆环带。kr_min和kr_max表示虚拟同心圆环带的内外圆半径kr_min，由图3(a)和图3(b)可以得到kr_min=其中min(kr)和max(kr)为初始圆环的内外圆半径。根据原始环带波数域径向分量kr的采样规则进行补0，如图3(c)所示。该步骤保证波数域频带利用率最大化，与传统直接插值方法相比，所提方法考虑到径向非均匀采样的特点先通过补0操作确定覆盖所需同心方形环带的环带填充区域可以减少区间的逻辑计算和索引查找等运算量。

步骤2 径向插值，如图3(d)-图3(f)所示，确定所需同心方形带内波数域径向的均匀采样。

为了最大化利用有效频谱数据，需要将插值前数据即空心点插值到正方实心点上。如图3(e)所示，由于kr是非均匀的，假设原始径向kr的采样个数为NR，先将kr均匀化

图3(e)为某一原始数据点位置与待插数据点位置的几何关系，可以得到krv(n)的表达式，径向插值后数据如图3(f)所示

步骤3kx-ky平面正交分区 ，如图3(g)-图3(h)所示。

经过步骤2，图3(f)的同心方形带填充呈现出分维度分块的均匀性。根据正方形对角线进行2维的正交区间划分，获得两组数据，每组数据包含两个两两不交叠相对的数据块，分别记为基本垂直点数据和基本水平点数据。图3(g)为基本垂直点数据Sv，该区间数据在波数距离维ky是均匀分布的，但是波数方位维kx为非均匀；图3(h)为基本水平点数据SH，该区间数据在波数方位维kx是均匀分布的，但是波数距离维ky为非均匀。根据波数域分布特点，可以对2维插值进行分区降维处理，实现对kx和ky维度的解耦。

步骤4 两个维度独立的1维插值，如图3(i)-图3(j)所示。

根据分区，将基本垂直点数据Sv沿着方位维利用3次样条插值法对kx进行内插，基本水平点数据SH沿着距离维利用3次样条插值法对ky进行内插。注意到在以正方形对角线划分数据时，对于正方形上半部分数据，从上往下每行数据点的个数依次减少2个，下半部分则反之。所以对于基本垂直点数据Sv的上半部分进行内插时，从上往下每行插值点的个数依次减少2个，下半部分则相反，如图3(i)所示。空心点为插值前的数据，实心点代表插值后的数据。对基本水平点数据SH亦是如此，基本水平点数据SH重采样后数据如图3(j)实心点所示。

步骤5 数据拼接，如图3(k)所示。

将基本垂直点数据Sv和基本水平点数据SH重采样后的数据点进行拼接得到完整的直角坐标系下分布的均匀数据，拼接后完整的数据如图3(k)所示，其中实心点为插值后的有效数据，空心圆点为置0数据点。最后将处于直角坐标系下均匀分布的数据进行3维傅里叶逆变换即可得到重建图像。

4 算法运算量分析

本节对比所提插值方法与逐点遍历插值方法的运算量，以算法所需的浮点运算量(FLoating point Operations Per Second, FLOPS)来进行衡量[10]。令Nker是1维插值核长度，角度维度采样点数为Nθ，距离维采样点数Nf，高度维采样点数为Nz，插值后均匀网格的方位维和距离维大小分别为Nx和Ny，即Nx=Ny。

2维逐点遍历(细胞单元法)插值所需运算量为

同心方形网格插值所需运算量为

分析所提插值算法与传统算法的运算量，定义运算量比γ=C2/C1，将式(22)和式(23)代入可得

Nker一般取8,16或32，这里取Nker=8，且由图3(i)可知，Nf ＜Nx/2，则

CSAR场景下通常取Nθ ＜4Nx，所以γ ＜1/6。可知所提插值处理方法与传统插值处理方法运算量比小于1/6，说明所提插值方法运算量小，实际成像中耗时少于传统插值处理方法，理论上提升的速度最少可达到6倍以上。

5 实验验证

本节通过数值仿真和实测数据对比分析本文所提插值方法与传统2维逐点遍历插值方法的成像时间与成像效果。实验参数见表1。成像处理的软件平台为MATLAB R2018a，为了更好地对比分析两种插值方法成像质量和成像时间，成像处理过程中不使用任何滤波处理和并行计算。

表1 实验参数

5.1 单点目标仿真验证

仿真实验中设置的点目标为距离场景圆心0.01 m、角度方向为20°、高度方向为0.01 m的单点目标，则点目标在直角坐标系下坐标为(0.009 4, 0.003 4,0.01)，其中X表示方位维，Y表示距离维，Z表示高度维。图4给出了传统2维逐点遍历插值算法(传统算法)与本文算法成像的剖面图。

图4 点目标成像剖面图

从3个方向剖面图可以看出，两种算法成像结果位置位置基本一致。为了定量评估点目标成像性能，表2给出了点目标在3个维度的分辨率、峰值旁瓣比和积分旁瓣比，其运算量参数对应于表1。由表2可知，改进算法成像3个维度的分辨率与传统算法成像的分辨率基本一致，各项测量参数非常接近，说明两种算法点目标成像质量相差不大，并且本文算法插值处理时间只有传统遍历插值处理时间的14%，提升速度达到7倍，说明了本文算法大大提高了成像的速度，验证了算法的高效性。

表2 点目标成像质量参数

5.2 多点目标仿真

以距离场景圆心0.1 m、高度维间隔为0.01 m构建柱形目标点阵进行仿真。图5给出了本文算法及传统算法成像3个维度的成像投影图。从投影图可以看出，成像位置与仿真设置的位置一致，且设置的点目标均能清晰分辨出，验证了算法成像的精确性。

图5 多点目标成像

5.3 实测数据验证

实验采用矢量网络分析仪、喇叭天线、高精度旋转步进电机和竖直向步进电机搭建实验平台如图6(a)所示。图6(b)为测试目标铝箔。

图6 柱面成像实验成像场景

图7给出了目标物品用两种算法成像的3个维度投影图，从投影图可以看出成像结果整体均匀程度和聚焦质量良好，成像目标清晰。对于无参考图像的质量评价，采用图像熵和图像信噪比来说明两种算法的成像效果[15]。传统算法3个维度投影图的图像熵分别为5.6, 2.2, 2.1，改进算法3个维度投影图的图像熵分别为3.6, 2.1, 1.1，由此看出改进算法的图像熵更小，图像聚焦效果较传统算法好；传统算法3维投影图的图像信噪比分别为21.6, 21.3, 19.3，改进算法3维投影图的图像信噪比为22.4, 23.8,24.2，通过信噪比可知改进算法的图像信噪比更大，图像所含噪声较少。

图7 实测数据成像图

表3给出了两种算法的运算量和成像时的插值处理时间，可以看出本文算法的运算量与传统算法相比减少了一个数量级，而且本文算法的插值处理时间约为传统2维逐点遍历插值处理时间的14%，提升速度达到7倍，符合算法复杂度理论分析结果。虽然表3中运算量减小倍数10和实际处理时间提速比7不完全相等，是因为在计算算法所需的浮点运算量时，忽略逻辑运算、索引查找等运算。并且在实际实现时，所耗时间还和代码的编写风格以及成像处理的软件平台MATLAB R2018a的运行模式有关。综合来说，本文改进的算法插值部分的运算量大幅度减少，插值处理时间短，在保证成像质量的基础上，大幅提高了成像算法的效率。

表3 理论运算量与插值处理时间

6 结论

本文对CSAR成像算法插值处理进行研究，针对在大角度下方位向与距离向波数域耦合性强，传统2维逐点遍历插值处理运算量大和耗时长的问题，提出了基于同心方形网格插值处理的方法。所提方法实现了方位向与距离向波数域解耦合，只涉及两个独立的1维插值处理，易于后期硬件实现时的并行化处理。实验结果表明，在保证成像质量的基础上，本文算法的插值处理时间约为传统2维逐点遍历插值处理时间的14%，提升速度达到7倍。因此，本文提出的基于同心方形网格插值处理的方法能够有效解决CSAR 3维成像在大角度下方位向与距离向波数域耦合性强、2维插值处理运算量大导致耗时长的问题。