1.   引言

近年来，当在传统电磁波上加载轨道角动量调制时，可在空间中叠加形成涡旋电磁波，利用不同形式天线可产生不同特性的轨道角动量电磁波，从而适用于不同的应用需求^[1-3]。特别地，在阵列雷达成像技术中，利用涡旋电磁波照射，同一距离单元内不同方位向的目标所在位置的辐射场相位波前存在一定的差异性，从而使雷达接收到的目标散射回波携带有更多的目标信息，进而通过方位向的信息解耦最终可实现目标凝视成像^[4-6]。对于涡旋电磁波照射下的目标成像，利用照射波束内波前的差异性，并且基于不同本征模态的正交性，通过多模态的照射可有望提升对目标的分辨性能^[7]。在电磁涡旋成像技术的基础上，结合合成孔径雷达(Synthetic Aperture Radar, SAR)成像或逆合成孔径雷达(Inverse Synthetic Aperture Radar, ISAR)成像特性，可进一步实现对目标或场景的三维(Three-Dimensional, 3D)成像^[8,9]。

在太赫兹频段，由于雷达信号载频高，极易发射大带宽信号，实现较高的距离分辨率^[10,11]，另外，相比于微波雷达，由于太赫兹波的波长短，在合成孔径成像或者逆合成孔径成像模式下可以获得更高的横向分辨率^[12-14]。文献[9]将太赫兹波的优势与涡旋电磁波相结合，建立了基于太赫兹电磁涡旋ISAR的新体制成像模型，实现了对目标的高分辨三维成像，然而在已有的成像算法中，俯仰向信息需要通过极坐标下三角函数和方位向的插值得到，这是一个间接求解的过程，俯仰维成像分辨率较低。

为进一步提高目标俯仰维的成像分辨率，本文在前期工作的基础上^[9]，提出了一种基于稀疏贝叶斯学习(Sparse Bayesian Learning, SBL)的图像重建方法，可直接在直角坐标系下求解目标的三维散射分布函数。利用稀疏贝叶斯学习方法进行图像重建时，通过稀疏表示模型建立观测矢量与待重建信号之间的关系，进而直接进行求解，求解过程能够自动学习重建模型中的未知参数，且人工参数设置对重建结果的影响较少^[15]。文中首先推导建立了电磁涡旋三维成像的回波模型；其次，构建了目标三维成像的稀疏表示模型，对SBL图像重建过程进行了详细介绍；最后，通过设置不同的成像实验场景，对不同成像方法的成像性能进行了对比分析，并对不同信噪比(Signal-to-Noise Ratio, SNR)条件下SBL方法的重构性能开展了仿真实验。

2.   电磁涡旋三维成像模型与方法

2.1   电磁涡旋三维成像模型

在传统ISAR成像模型下，利用带宽信号和雷达与目标之间的相对运动可实现距离-方位二维成像。另外，在涡旋电磁波的照射下，结合带宽信号可获得距离-方位角的二维分辨能力。因此，本文将二者结合，建立了基于电磁涡旋ISAR的三维成像模型，如图1所示，图中均匀圆阵表示雷达，通过在线性调频发射信号上加载轨道角动量调制而产生涡旋电磁波，假设目标上一理想散射点P的坐标可以表示为$(x,y,z)$，且目标绕Y轴旋转，旋转方向如图中箭头所示，转角为$\theta \in [ - \Delta \theta /2,{\text{ }}\Delta \theta /2]$, $\Delta \theta$表示目标旋转的最大角度，图中方位角$\varphi$为理想散射点P在$XOY$平面内的投影与X 轴的夹角。

图  1  基于电磁涡旋ISAR的三维成像几何

Figure  1.  Sketch map of the 3D imaging geometry based on electromagnetic vortex ISAR

下载: 全尺寸图片 幻灯片

一般来说，当目标由大量散射点构成，雷达接收到的目标总回波可表示为式(1)积分形式^[9]

其中，$f(x,y,z)$表示目标三维散射分布函数，$k$为波数，$l$为轨道角动量模式数。

当雷达和目标之间转动角度较小时，对式(1)所示雷达回波分别在频率维$k$和方位维$\theta$进行二维傅里叶逆变换，再在拓扑荷域$l$进行一维傅里叶变换，即可得到目标的散射分布函数

其中，$f(x,z,\varphi )$为方位-距离-方位角上的目标三维散射分布函数。因此，在已知目标方位$x$、距离$z$、方位角$\varphi$信息的情况下，要获得直角坐标系下的目标三维散射分布函数，即$f(x,y,z)$，首先需要根据方位$x$、方位角$\varphi$的值计算得到俯仰向Y 的分布，计算表达式为$y = x \cdot \tan \varphi$，然后在Y 轴上将计算得到的非均匀网格插值为均匀网格，最终可得$f(x,y,z)$，实现目标三维成像^[9]。

2.2   基于稀疏贝叶斯学习的成像方法

根据2.1节中建立的成像模型，式(2)中的成像过程可改写为

首先，对雷达回波采用卷积逆投影(Convolution Back-Projection, CBP)的方法进行成像，令$G(k,\theta ,l) =$$k \cdot {S_r}(k,\theta ,l)$，且$g(w,\theta ,l)$为$G(k,\theta ,l)$的一维傅里叶逆变换。因此，式(3)中的第1步积分可表示为

其中，$h(x,z,l)$为不同轨道角动量模式数时目标二维散射分布函数，在不同$\theta$值计算式(4)中对应的$- x\sin \theta + z\cos \theta$，然后将$g(w,\theta ,l)$沿不同的$\theta$进行叠加，即可计算得到$h(x,z,l)$，这就是CBP成像算法的基本步骤^[16]。

接下来，将式(4)代入式(3)，可以得到

式(5)完全符合傅里叶变换形式，因此，利用傅里叶变换与逆变换相对应的性质可知

将$\varphi = \arctan (y/x)$代入式(6)等号右侧的积分表达式中，可得

下面，再对X方向和Z方向的二维成像网格进行离散化，将相应方向的成像区域分别划分为M, N 个网格，在每一个离散值${x_m}$, ${z_n}$对应的目标二维散射分布函数$h(x,z,l)$处，采用SBL方法进行一维重构，成像几何如图2所示，将式(7)进一步改写为

图  2  基于SBL方法的三维成像几何

Figure  2.  3D imaging geometry based on the SBL method

下载: 全尺寸图片 幻灯片

根据稀疏贝叶斯恢复的基本原理，在式(8)的基础上构建如下信号模型

其中，${\boldsymbol{H}}$表示观测矢量，${\boldsymbol{A}}$为对应的测量矩阵，${\boldsymbol{y}}$为待求解的俯仰向Y 的散射系数矢量，${\boldsymbol{n}}$表示噪声。假设雷达发射信号的轨道角动量调制模式数的采样点数为L，俯仰向的离散网格数为P，将式(9)改写为

式(10)中忽略了${x_m}$, ${z_n}$这两个重复项，测量矩阵${\boldsymbol{A}}$可进一步表示为

在每一个${x_m}$, ${z_n}$处，重复进行式(10)中的稀疏贝叶斯重构，最终可直接得到直角坐标系下的目标三维散射分布函数$f(x,y,z)$。

在上述重构过程中，${x_m}$, ${z_n}$各有总共M, N 个取值，这样就需要完成M×N次贝叶斯重构，当M, N 取值较大时，重构过程计算复杂度高且算法运行时间较长。在实际成像场景中，目标散射点通常是成稀疏分布的，并不是所有俯仰向都存在目标散射点，因此，为降低计算量，只需要在有目标散射点的距离-方位切片上进行一维重构。当$l = 0$时，雷达发射信号仅为线性调频信号，此时，$h(x,z,0)$为传统ISAR成像结果，即目标在$XOZ$平面上的投影，当$h({x_m},{z_n},0)$的幅度很小时，本文认为在${x_m}$, ${z_n}$所在的距离-方位切片上不存在目标散射点。因此，本文对实际成像过程作进一步优化以节约计算成本，与文献[17]中不同的是，为避免距离-方位切片上强散射点的影响，本文提出分区域幅度阈值设置方法，首先在距离-方位切片上寻找局部最大值点$h({x_m},{z_n},0)$，以各局部最大值点$h({x_m},{z_n},0)$为中心划分$W \times W$个分辨单元为选定区域$\Omega$, $W$通常选择为1/2距离-方位向点扩散函数主瓣宽度所占的分辨单元个数，设定幅度阈值为$\eta \cdot h({x_m},{z_n},0)$。当该区域内$\mathop {h(x,z,0)}\limits_{x,z \in \Omega } \ge \eta \cdot h({x_m},{z_n},0)$时，记为$h({x_{{m_i}}},{z_{{n_i}}},0)$，该距离-方位切片包含目标散射点，将$h({x_{{m_i}}},{z_{{n_i}}},l)$作为观测矢量，构建测量矩阵，求解$f({x_{{m_i}}},{z_{{n_i}}},y)$；当$\mathop {h(x,z,0)}\limits_{x,z \in \Omega } < \eta \cdot h({x_m},{z_n},0)$时，不包含目标散射点，该距离-方位分辨单元内俯仰向散射强度均设置为0。最终，优化后的SBL成像处理流程如下。

Step 1：设置目标散射点，根据式(1)生成三维采样下的雷达回波${S_r}(k,\theta ,l)$；

Step 2：在不同轨道角动量模式数$l = {l_1},{l_2}, \cdots ,l{}_N$，分别对雷达回波进行二维CBP成像，得到目标二维散射分布函数$h(x,z,l)$；

Step 3：在距离-方位切片$h(x,z,0)$内，寻找所有局部最大值$h({x_m},{z_n},0)$；

Step 4：以$h({x_m},{z_n},0)$为中心分别划分区域$\Omega$，为不同区域分别设置幅度阈值$\eta \cdot h({x_m},{z_n},0)$，当$\mathop {h(x,z,0)}\limits_{x,z \in \Omega } < \eta \cdot h({x_m},{z_n},0)$，设置$\mathop {f(x,z,y)}\limits_{x,z \in \Omega } = 0$；当$\mathop {h(x,z,0)}\limits_{x,z \in \Omega } \ge \eta \cdot h({x_m},{z_n},0)$，将$h({x_{{m_i}}},{z_{{n_i}}},l)$作为观测矢量，根据式(11)构建测量矩阵${\boldsymbol{A}}$，建立信号模型${{{\boldsymbol{H}} = {\boldsymbol{A}}}} \cdot {{{\boldsymbol{y}} + {\boldsymbol{n}}}}$；

Step 5：采用稀疏贝叶斯学习的方法进行一维重构得到$f({x_{{m_i}}},{z_{{n_i}}},y)$；

Step 6：计算得到的所有离散值$f({x_m},{z_n},y)$，最终得到目标三维成像结果。

3.   仿真结果与分析

3.1   实验场景设置

仿真实验中，假设雷达发射信号中心频率为330 GHz，信号带宽为20 GHz，频率采样间隔为0.1 GHz。圆形阵列的阵元半径设置为500$\lambda$，最大轨道角动量模式数取${l_{\max }}$=30，相应的轨道角动量模式数的变化范围为[–30, 30]，成像时设置方位向最大转角为$\Delta \theta$=5°，对应转角$\theta \in [ -2{\text{.}}{{\text{5}}^{\text{o}}}{\text{, 2}}{\text{.}}{{\text{5}}^{\text{o}}}{\text{]}}$，仿真时设置方位向成像场景宽度为0.6 m，对应满足不混叠距离的角度采样点数可设置为126，设置各散射点的散射强度均为1，其几何位置关系如图3所示。

图  3  目标散射点分布

Figure  3.  The distribution of point targets

下载: 全尺寸图片 幻灯片

3.2   电磁涡旋三维成像结果

成像仿真中，本文的稀疏贝叶斯求解方法采用变分贝叶斯推断(Variational Sparse Bayesian Inference, VSBI)^[18]的方法，首先对图3所示的目标散射点根据式(1)生成雷达回波，然后再采用本文提出的SBL方法进行三维成像，成像过程中，设置X, Y, Z 3个方向的离散网格数均为$M = N = P = 601$。

图4(a)为$l = 0$时的归一化的距离-方位切片$h(x,z,0)$，采用分区域幅度阈值设置方法时，根据点扩散函数的主瓣宽度，设置$\eta = 0.7$，得到如图4(b)所示的分区域幅度阈值划分结果。区域划分时，将$l = 0$时的归一化的距离-方位切片视为一幅图像，采用图像处理中寻找局部最大值的经典八连通方法进行处理^[19]，首先会找到散射点P₁所在分辨单元的局部最大值点，以该点幅度值的$\eta$倍为阈值，根据距离分辨率和方位分辨率与离散网格数计算选取$W{\text{ = }}40$，对选定区域$\Omega$内的分辨单元作幅度阈值划分，保证了后续稀疏贝叶斯重构方法的准确性。若采用文献[17]中的幅度阈值方法，为所有$h(x,z,0)$设置单一幅度阈值，结合图3和图4(a)可以看出，除散射点P₁之外，距离-方位平面上的其他散射点均由2个目标散射点投影形成，使得散射点P₁的相对散射强度较弱，仅为0.35(远低于1)，若设置阈值$\eta = 0.7$，最终三维成像结果中将不包含散射点P₁，而若将幅度阈值降低，设置为0.3，得到如图4(c)所示的划分结果。对比图4(b)和图4(c)可知，分区域幅度阈值设置方法得到的各散射点所在区域分辨单元个数是大致相同的，而图4(c)中虽然能够得到距离-方位切片上5个散射点所在区域，但P₁点所在的分辨单元个数明显少于其他散射点，这样会造成P₁点回波能量的大量损耗，进而影响后续重构结果。另外，对比图4(b)和图4(c)中选定区域分辨单元总数，图4(c)中除P₁外其他4个散射点由阈值0.3得到的区域面积过大，图4(c)的分辨单元总数约超过图4(b)中一倍，极大地增加了后续稀疏重构的计算量。

图  4  l =0时距离-方位切片及不同幅度阈值划分结果

Figure  4.  The range-azimuth image at l =0 and the results of different amplitude threshold setting method

下载: 全尺寸图片 幻灯片

在分区域辐射阈值设置方法的基础上，图5为采用本文提出的VSBI方法重构的三维成像结果，为能直接获得更清晰的观测效果，将目标的三维成像结果分别投影到距离、方位和俯仰3个不同坐标平面内，得到不同的二维成像结果。在利用式(11)建立测量矩阵时，各矩阵元素的幅值${x_m}/(x_m^2 + y_p^2)$根据散射点位置的不同而存在一定差异性，从而使得最终图5中的重构得到的各散射点之间的幅度不尽相同。最后，从图5中可以看出，包括散射点P₁在内的所有散射点位置均得到了准确重构，验证了本文所提的成像方法的有效性。

图  5  基于SBL方法的三维成像结果

Figure  5.  3D imaging results based on SBL method

下载: 全尺寸图片 幻灯片

3.3   不同成像方法的对比分析

本文提出的SBL方法对成像效果的提升主要体现在俯仰向Y的重构上，为进一步说明本文方法的优势，下面将其与文献[9]中所提的基于三维快速傅里叶变换(Fast Fourier Transform, FFT)和基于CBP和功率谱密度(Power Spectrum Density, PSD)估计两种成像方法进行对比分析。仿真实验中，在同一距离-方位分辨单元内设置不同俯仰向位置的4个点，散射点的三维坐标分别为Q₁(0.15 m, 0, 0), Q₂(0.15 m, 0.01 m, 0), Q₃(0.15 m, 0.1 m, 0), Q₄(0.15 m, 0.115 m, 0)，散射点分布如图6(a)所示，各散射点X 轴坐标相同，间隔分布在Y 方向。

图  6  目标散射点分布及其成像结果对比

Figure  6.  The distribution of targets and comparison of imaging results

下载: 全尺寸图片 幻灯片

理想条件下，在$x = 0.15\;{\text{m}}$, $z = 0$时的俯仰维剖面图对比如图6(b)所示，对比图中成像结果可知，对于距离很近的Q₁和Q₂两点、Q₃和Q₄两点，文献[9]中提出的基于三维FFT和基于CBP-PSD的两种成像方法均无法在俯仰向实现分辨，而本文提出的SBL方法能够准确地重构出Q₁, Q₂, Q₃和Q₄ 4个散射点的位置，且成像分辨率优于0.01 m。

为进一步衡量SBL成像方法的性能，计算目标散射点的重构误差，首先，定义散射系数的最小均方误差(Mean Square Error, MSE)

其中，$\tilde y$和$\hat y$分别表示散射系数的理论值和重构得到的估计值。在不同信噪比条件下，计算目标重构误差MSE，信噪比从0 dB到20 dB间隔2 dB变化，500次蒙特卡洛仿真条件下的重构误差统计结果如图7所示，可以看出，SBL成像方法的重构性能随信噪比的增大而提升，且SBL成像方法的重构误差显著低于三维FFT方法和基于CBP-PSD的成像方法。

图  7  重构散射系数MSE随信噪比变化情况

Figure  7.  MSE of reconstructed scattering coefficient as a function of the SNR

下载: 全尺寸图片 幻灯片

4.   结束语

本文将太赫兹ISAR与电磁涡旋相结合，建立了基于电磁涡旋ISAR的成像稀疏表示模型，并提出了基于稀疏贝叶斯学习的目标三维重构方法，通过SBL方法可直接由雷达回波重构目标在空间直角坐标系中的三维信息，极大地简化了成像求解过程，提高了成像分辨率。后续研究将针对涡旋电磁波所特有的贝塞尔函数幅度特性并结合实际电磁涡旋雷达实验系统进行深入分析，为新体制雷达三维成像技术的发展提供参考和借鉴。