基于稀疏贝叶斯学习的太赫兹电磁涡旋三维成像方法

蒋彦雯 范红旗 李双勋

(国防科技大学电子科学学院 长沙 410073)

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轴旋转，旋转方向如图中箭头所示，转角为θ∈[-Δθ/2,Δθ/2],Δθ表示目标旋转的最大角度，图中方位角φ为理想散射点P在XOY平面内的投影与X轴的夹角。

图1 基于电磁涡旋ISAR的三维成像几何Fig.1 Sketch map of the 3D imaging geometry based on electromagnetic vortex ISAR

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

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

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

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

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

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

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

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

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

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

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

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

图2 基于SBL方法的三维成像几何Fig.2 3D imaging geometry based on the SBL method

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

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

式(10)中忽略了xm,zn这两个重复项，测量矩阵A可进一步表示为

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

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

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

Step 2：在不同轨道角动量模式数l=l1,l2,...,lN，分别对雷达回波进行二维CBP成像，得到目标二维散射分布函数h(x,z,l)；

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

Step 5：采用稀疏贝叶斯学习的方法进行一维重构得到f(xmi,zni,y)；

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

3 仿真结果与分析

3.1 实验场景设置

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

3.2 电磁涡旋三维成像结果

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

图3 目标散射点分布Fig.3 The distribution of point targets

图4 l=0时距离-方位切片及不同幅度阈值划分结果Fig.4 The range-azimuth image atl=0 and the results of different amplitude threshold setting method

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

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

图5 基于SBL方法的三维成像结果Fig.5 3D imaging results based on SBL method

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

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

图6 目标散射点分布及其成像结果对比Fig.6 The distribution of targets and comparison of imaging results

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

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

图7 重构散射系数MSE随信噪比变化情况Fig.7 MSE of reconstructed scattering coefficient as a function of the SNR

4 结束语

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