2025, Vol. 55
2. 青岛海洋科学与技术国家实验室 区域海洋动力学与数值模拟功能实验室，山东 青岛 266237
毫米波雷达是一种利用毫米波电磁频谱来探测和定位物体的雷达系统。它的核心原理和传统雷达类似，但使用的频率更高，通常在30~300 GHz之间，对应的波长在1~10 mm之间。与可见光和红外传感器相比，毫米波雷达受天气的影响较小，且具有全天时全天候主动工作的优点[1-5]。同时，毫米波雷达的频率高于通常的遥感雷达，比起常用的L、X等波段雷达[6]，能够提供更高的空间分辨率。毫米波雷达天线尺寸较小，使得其更易于集成到小型合成孔径雷达(Synthetic aperture radar，SAR)系统上[7]，从而使其易操作、成本低。由于毫米波具有一定的穿透性，并且对人体无害，近年来在安检、生命体征监测、遥感、无损检测等领域[8-12]被广泛关注。
而应用毫米波SAR实现目标精准快速成像，涉及到成像算法研究。当前，SAR成像算法主流的有两类：一类是时域算法，一类是频域算法。后向投影算法(Back projection algorithm，BPA)作为最早的SAR成像时域算法，起源于计算机断层扫描的成像技术，这种方法最早出现在医学成像领域[13]，后期的时域算法都是基于对BPA的优化。对于频域成像算法而言，虽然不如时域算法压缩性能优异，但是由于其在成像过程中借助了快速傅里叶变换，从而大幅提升了运算效率，例如：距离多普勒算法(Range doppler algorithm，RDA)[14]、距离徙动算法(Range migration algorithm，RMA)[15]等。当前，近场二维SAR成像中使用最多的算法是BPA和RMA。与BPA和RMA相比，RDA中通过距离和方位上的频域操作，实现了高效的模块化处理要求，其成像效率远远高于BPA，同时相比于RMA又具有一维操作的简便性[16]。然而，RDA最初是为了机载或星载等远场成像场景而设计的[17-24]，其近似方位向调频率不随距离的变化而变化，这个近似能够提升计算效率，但是并不适用于近场成像。本文在传统RDA的基础上对其进行优化以适应近场成像条件，对各个距离门的方位向单独进行脉冲压缩，提升了近场条件下的RDA的二维成像质量。
通常，现实中的目标具有x、y、z三个空间维度，二维雷达成像只能给出目标在垂直于雷达视向平面内的散射特征，如果想要获取目标更为完整的三维散射特征信息，需要进行三维成像。毫米波SAR的近场三维成像技术是对二维成像技术的扩展，能提供目标更加丰富和详细的三维空间散射信息。相对于二维SAR成像技术而言，当前，针对毫米波雷达的三维成像算法的研究较少。目前，毫米波雷达在三维成像中最常用的算法是RMA，早在2001年，美国西北太平洋国家实验室的David M S等[25]就通过RMA实现了三维毫米波成像；在2019年，德州大学Muhammet Emin Yanik等用77G毫米波雷达就通过RMA实现三维成像和隐藏目标成像[26]。然而，对于毫米波雷达近场的三维成像技术而言，当前依然缺少针对RDA进行三维成像的研究；此外，大多学者在毫米波雷达三维成像研究过程中，通常的成像对象是跨越较少距离门的小目标或平面目标。但是，对于跨越多个距离门的大目标而言，由于高度向距离徙动的存在和平台运动误差等因素会导致实验数据获取的过程中出现数据错位的情况，进而显著影响成像质量。
针对以上这些问题，本文基于实际实验数据，通过引入二维互相关算法用以矫正数据错位，将已做优化的二维RDA拓展至三维成像，并对比RDA与RMA的三维成像结果。实验结果显示：引入二维互相关算法可有效消除数据错位的影响，显著提升三维成像质量，并且成像效果也优于RMA。
1 算法原理及流程 1.1 RDA原理雷达想要提升距离分辨率，就要发射带宽更大即时宽更窄的脉冲信号，但是窄脉冲意味着发射能量小，导致探测距离短。对于一般的脉冲信号，带宽和时宽互为倒数，即带宽和时宽不能同时增大，为了解决这一问题，SAR一般采用线性调频(Linear frequency modulation, LFM)信号，雷达发射的LFM信号可以表示为
| $ \begin{gather*} \boldsymbol{S}\left(t_{\mathrm{k}}, t_{\mathrm{m}}\right)= \\ \operatorname{rect}\left(\frac{t_{\mathrm{k}}}{T_{\mathrm{P}}}\right) \cdot \mathtt{ω}_{\mathrm{a}}\left(t_{\mathrm{m}}\right) \cdot \exp \left[\mathrm{j} 2 \mathsf{π} \left(f_{\mathrm{c}} t+\frac{1}{2} K_{\mathrm{r}} t_{\mathrm{k}}^{2}\right)\right] 。\end{gather*} $ | (1) |
式中: rect(·)为矩形窗函数；tk为脉冲传播时间，即快时间；TP为脉冲调频时间，即脉宽。ωa(·)为方位包络，tm表示慢时间，其大小为当前脉冲数与脉冲重复时间的乘积。fc为中心频率，t表示全时间，即快时间与慢时间的加和，Kr为距离向调频率。
对于回波信号距离向脉压，可以采取解线频调的方式，解线频调是用一时间固定，而频率、调频率相同的LFM信号作为参考信号，用它和回波作差频处理，差频输出的表达式为[27]
| $ \begin{align*} \boldsymbol{S}_{\mathrm{if}}\left(t_{\mathrm{k}}, t_{\mathrm{m}}\right)= & \sigma_{0} \operatorname{rect}\left(\frac{t_{\mathrm{k}}-\tau}{T_{\mathrm{p}}}\right) \cdot \mathtt{ω}_{\mathrm{a}}\left(t_{\mathrm{m}}\right) \cdot \exp \left(-\mathrm{j} 2 \mathsf{π} K_{\mathrm{r}} t_{\mathrm{k}} \tau\right) \cdot \\ & \exp \left(-\mathrm{j} 2 \mathsf{π} f_{\mathrm{c}} \tau\right) \cdot \exp \left(\mathrm{j} \mathsf{π} K_{\mathrm{r}} \tau^{2}\right) 。\end{align*} $ | (2) |
通过差频处理后，全时间t在回波信号与参考信号的相干检波时消去，σ0为场景中某点的后向散射系数，慢时间tm也在雷达到该点目标的双程时延τ中体现，$\tau=2 R\left(t_{\mathrm{m}}\right) / \mathrm{c}, R\left(t_{\mathrm{m}}\right)$为雷达到目标点的距离。式(2)有三个相位，其中第三个相位称为剩余视频相位(Residual video phase，RVP)，这个相位是解线频调方法特有的，对式(2)的快时间维为做傅里叶变换并补偿掉RVP项，最终距离脉压后的表达式为
| $ \begin{gather*} \boldsymbol{S}_{\mathrm{rc}}\left(f_{\mathrm{k}}, t_{\mathrm{m}}\right)=\sigma_{0} T_{\mathrm{P}} \operatorname{sinc}\left[T_{\mathrm{P}}\left(f_{\mathrm{k}}+K_{\mathrm{r}} \tau\right)\right] \cdot \\ \omega_{\mathrm{a}}\left(t_{\mathrm{m}}\right) \cdot \exp \left(-\mathrm{j} 2 \mathsf{π} f_{\mathrm{c}} \tau\right) 。\end{gather*} $ | (3) |
在正侧视或小斜视角的情况下，R(tm)可近似为
| $ \begin{equation*} R\left(t_{\mathrm{m}}\right) \approx R_{0}+\frac{V_{\mathrm{r}}^{2} t_{\mathrm{m}}^{2}}{2 R_{0}} 。\end{equation*} $ | (4) |
式中：R0为目标点到雷达飞行轨迹的最短距离；Vr为雷达移动速度，将R(tm)代入式(3)可得：
| $ \begin{gather*} \boldsymbol{S}_{\mathrm{rc}}\left(f_{\mathrm{k}}, t_{\mathrm{m}}\right)=\sigma_{0} T_{\mathrm{P}} \operatorname{sinc}\left[T_{\mathrm{P}}\left(f_{\mathrm{k}}+K_{\mathrm{r}} \tau\right)\right] \cdot \mathtt{ω}_{\mathrm{a}}\left(t_{\mathrm{m}}\right) \cdot \\ \exp \left(-\mathrm{j} 4 \mathsf{π} \frac{f_{\mathrm{c}} R_{0}}{c}\right) \cdot \exp \left(-\mathrm{j} \mathsf{π} K_{\mathrm{a}} t_{\mathrm{m}}^{2}\right) 。\end{gather*} $ | (5) |
其中Ka为方位向调频率，可近似表示为
| $ \begin{equation*} K_{\mathrm{a}} \approx \frac{2 V_{\mathrm{r}}^{2}}{\lambda R_{0}} 。\end{equation*} $ | (6) |
式中λ为雷达发射的电磁波的波长，利用驻定相位原理(POSP)，可以得到方位向上的时频关系为
| $ \begin{equation*} f_{\mathrm{m}}=-K_{\mathrm{a}} t_{\mathrm{m}} 。\end{equation*} $ | (7) |
将式(7)代入式(5)中，可得方位向傅里叶变换后的信号为
| $ \begin{gather*} \boldsymbol{S}_{\mathrm{rd}}\left(f_{\mathrm{k}}, f_{\mathrm{m}}\right)=\mathrm{FFT}_{t_{\mathrm{m}}}\left(\boldsymbol{S}_{\mathrm{rc}}\left(f_{\mathrm{k}}, t_{\mathrm{m}}\right)\right)= \\ \sigma_{0} T_{\mathrm{P}} \operatorname{sinc}\left[T_{\mathrm{P}}\left(f_{\mathrm{k}}+K_{\mathrm{r}} \tau\right)\right] \cdot W_{\mathrm{a}}\left(f_{\mathrm{m}}\right) \cdot \\ \exp \left(-\mathrm{j} 4 \mathsf{π} \frac{f_{\mathrm{c}} R_{0}}{c}\right) \cdot \exp \left(\mathrm{j} \mathsf{π} \frac{f_{\mathrm{m}}^{2}}{K_{\mathrm{a}}}\right) 。\end{gather*} $ | (8) |
Wa(fm)为方位天线方向图ωa(tm)的频域形式，二者在形状上一致，在距离多普勒域，位于同一距离门处的不同目标的距离徙动轨迹一致，可以较为高效的进行距离徙动矫正和方位向脉冲压缩。由式(4)可得距离徙动量可以表示为
| $ \begin{equation*} \Delta R\left(f_{\mathrm{m}}\right)=R\left(t_{\mathrm{m}}\right)-R_{0} \approx \frac{V_{\mathrm{r}}^{2} t_{\mathrm{m}}^{2}}{2 R_{0}} 。\end{equation*} $ | (9) |
距离徙动矫正可以基于距离向上的sinc函数插值来实现，在距离徙动矫正后，需要对方位向进行脉冲压缩，由式(8)，其方位向匹配滤波器的表达式可以写为：
| $ \begin{equation*} H_{\mathrm{a}}=\exp \left(-\mathrm{j} \mathsf{π} \frac{f_{\mathrm{m}}^{2}}{K_{\mathrm{a}}}\right) 。\end{equation*} $ | (10) |
将滤波器变换到频域与式(8)相乘，再将方位向逆傅里叶变换到时域可得聚焦后的信号：
| $ \begin{equation*} \boldsymbol{S}_{\mathrm{ac}}\left(f_{\mathrm{k}}, t_{\mathrm{m}}\right)=\operatorname{IFFT}_{t_{\mathrm{m}}}\left(\boldsymbol{S}_{\mathrm{rd}}\left(f_{\mathrm{k}}, f_{\mathrm{m}}\right) \cdot \operatorname{FFT}\left(\boldsymbol{H}_{\mathrm{a}}\right)\right) \text { 。} \end{equation*} $ | (11) |
由式(6)可知方位向调频率Ka为有关于R0的函数，在远场条件下，雷达距成像场景区域内目标点的斜距变化的幅度较小，一般取R0为成像场景中心参考点处到雷达飞行轨迹的最短距离，也就是说R0在整个压缩过程中是一个定值，相应的Ka也是一个定值，由于远场数据量一般较大，这样做可以提升成像效率。但是在近场条件下，距离变化的幅度较大，Ka并不能近似为一个定值，同样的方位向频率轴fm也在变化，如果还是按照远场近似来进行方位向脉冲压缩，其效果较差。
现对传统RDA在近场条件下进行优化，不采用定值R0，而是循环单独对每一个距离门的方位向数据进行距离徙动矫正和脉冲压缩，根据式(6)和式(7)，其滤波器可构建为：
| $ \begin{equation*} \boldsymbol{H}_{\mathrm{a}}=\exp \left(-\mathrm{j} 2 \mathsf{π} \frac{x^{2}}{\lambda y}\right) 。\end{equation*} $ | (12) |
式中：x=Vrtm，为雷达方位坐标，y为每个距离门的距离，可表示为
| $ \begin{equation*} y=i \cdot \Delta y \text { 。} \end{equation*} $ | (13) |
其中Δy=c/2B，为雷达的距离分辨率，B为带宽，i=1, 2, 3…, n表示不同的距离门，最终方位向压缩完成的表达式可由式(11)得出：
| $ \begin{gather*} \boldsymbol{S}_{\mathrm{ac}}\left(f_{\mathrm{k}}, t_{\mathrm{m}} ; y\right)=\sigma_{0} T_{\mathrm{P}} \operatorname{sinc}\left[T_{\mathrm{P}}\left(f_{\mathrm{k}}+K_{\mathrm{r}} \tau\right)\right] \cdot \\ \operatorname{sinc}\left(t_{\mathrm{m}}\right) \cdot \exp \left(-\mathrm{j} 4 \mathsf{π} \frac{f_{\mathrm{c}} y}{c}\right) 。\end{gather*} $ | (14) |
关于距离徙动矫正，由上可将式(9)中的距离徙动量与雷达方位坐标x和距离y之间的关系表示为：
| $ \begin{equation*} \Delta R=\frac{x^{2}}{2 y}。\end{equation*} $ | (15) |
对于本文所涉及到的实验，雷达合成孔径长度为1 m，距离均限制在0.5~3 m之间。距离徙动量的变化如图 1所示，图中z轴代表偏移的距离单元数。
|
图 1 距离徙动量变化图 Fig. 1 Change of the range migration momentum |
由图 1可看出距离徙动量随雷达方位坐标成抛物线式变化，并且随距离的增加而减小，在本文涉及到的实验的成像区域内其距离徙动量均在五个距离分辨单元内，说明其距离徙动效应并不明显，但是为了增加本文算法的实用性，还是选择进行距离徙动矫正，然后再对方位向进行压缩，RDA近场优化后的二维成像流程如图 2所示。
|
图 2 RDA近场优化后二维成像流程图 Fig. 2 Flowchart of the 2D imaging process following the optimization of the RDA under near-field conditions |
原始数据在距离向傅里叶变换之后得到距离脉压完的数据，然后对距离向进行8倍的线性插值以提升最终的显示效果。在划分完成像区域后，需要对方位向进行处理，解线频调体制下基于sinc插值实现的距离徙动矫正最好在方位时域中进行。在矫正完后再对其进行方位向傅里叶变换，从第一个距离门开始，将当前距离门的滤波器傅里叶变换至频域并与该距离门的数据相乘，将相乘后的数据逆傅里叶变换回时域，完成该距离门的方位压缩，在遍历完所有距离门后划分成像区域并输出图像。
根据以上流程，基于毫米波雷达仿真了5个点目标的成像结果，雷达仿真参数如表 1。
|
表 1 雷达仿真参数 Table 1 Radar simulation parameters |
五个点目标与雷达之间的距离在0~20 m之间，雷达合成孔径长度为6 m，优化前的方位向匹配滤波器中的参考斜距选用中心点目标与雷达之间的距离，其优化前后成像结果对比如图 3所示。
|
((a)RDA优化前仿真结果图；(b)RDA优化后仿真结果图。(a)Illustration of the simulation results prior to the optimization of the RDA; (b)Illustration of the simulation results following the optimization of the RDA.) 图 3 RDA优化前后仿真结果图对比 Fig. 3 Comparison of the RDA simulation results before and after optimization |
由图 3(a)可看出除中心点压缩较好外，其余4点在方位向均有不同程度的展宽，这是因为仅以中心点为参考斜距时，只有中心点的回波能与滤波器响应，其余点的能量并不能得到很好的聚焦，反观图 3(b)，优化后的滤波器随距离的变化而变化，对各个距离门单独压缩，五个点目标均能聚焦的很好。
1.3 基于优化后RDA近场三维成像算法 1.3.1 二维互相关算法三维成像本文基于平面扫描模式，即雷达在距离零点处由高度向和方位向组成的平面内扫描。如图 4所示，x、y、z分别代表方位向、距离向以及高度向，雷达在扫描完一个高度层之后间隔一段距离扫描第二个高度层，间隔的距离需要满足高度向的奈奎斯特采样定理以免产生混叠。由于雷达是左右来回扫描，首先需要将偶数行的数据反向，然后将原始的整体二维回波矩阵按高度扫描层数排列成一个三维的矩阵，如图 5所示。
|
图 4 毫米波雷达三维成像几何构型图 Fig. 4 3D imaging geometry of millimeter-wave radar |
|
图 5 三维数据排列示意图 Fig. 5 Schematic diagram of the arrangement of 3D data |
由于平台移动误差以及在排列数据时将雷达上升时接收的回波排列到下一层，这都会导致数据在方位向上的错位。对于平面类型即跨越较少距离门的目标而言，一维互相关算法可以很好的矫正数据错位，但是对于类似飞机模型这种跨越多个距离门的非平面目标，高度向的距离徙动会导致数据会在距离向上也产生错位。在这类情况下，使用一维互相关算法效果较差，甚至会起到反作用。本文在一维互相关的基础上引入了二维互相关算法来同时对每一个高度层距离向和方位向错位的数据进行矫正，以提升其三维成像质量，二维互相关函数可以表示为[28]
| $ \begin{aligned} & ~~~~~~~\rho(u, \nu)= \\ & \frac{\sum\limits_{i=1}^M \sum\limits_{j=1}^N\left[\boldsymbol{f}_{k+1}^{u, v}(i, j)-\overline{\boldsymbol{f}}_{k+1}^{u, v}\right]\left[\boldsymbol{f}_k^{u, v}(i, j)-\overline{\boldsymbol{f}}_k^{u, v}\right]}{\sqrt{\sum\limits_{i=1}^M \sum\limits_{j=1}^N\left[\boldsymbol{f}_{k+1}^{u, v}(i, j)-\overline{\boldsymbol{f}}_{k+1}^{u, v}\right]^2 \sum\limits_{i=1}^M \sum\limits_{j=1}^N\left[\boldsymbol{f}_k^{u, v}(i, j)-\overline{\boldsymbol{f}}_k^{u, v}\right]^2}}。\end{aligned} $ | (16) |
式中：ρ(u, ν)就代表滑动位置(u, v)处的互相关系数；M，N分别为数据矩阵中距离向和方位向像素点数；第k+1层为需要进行数据对齐并挪位的位移层，第k层为参考层，即做互相关的模板。$f_{k+1}^{u, v}(i, j)$ 和$f_{k}^{u, v}(i, j)$分别表示位移层和参考层在滑动位置(u, v)处(i, j)元素的数据幅值，$\overline{\boldsymbol{f}}_{k+1}^{u, v}$ 和$\overline{\boldsymbol{f}}_{k}^{u, v}$分别表示位移层和参考层在滑动位置(u, v)处的均值。其数据对齐的简要原理图如下：
|
图 6 二维互相关算法原理图 Fig. 6 Schematic diagram of a two-dimensional cross-correlation algorithm |
每一个高度层的数据矩阵大小一致，均为M×N，位移层与参考层做二维互相关会得到一个(2M－1)×(2N－1)的矩阵，设为C，定义上述函数ρ(u, ν)表示矩阵C中第u行第v列的元素，即：
| $ \begin{equation*} \rho(u, \nu)=C_{u v} \text { 。} \end{equation*} $ | (17) |
要找到数据对齐所需要的位移量，就需要找到相关系数最大的位置，设使得ρ(u, ν)取最大值的索引为(u*, v*)，其满足：
| $ \left(u^*, v^*\right)=\arg \max\limits_{1 \leqslant u \leqslant 2 M-1,1 \leqslant v \leqslant 2 N-1} \rho(u, v)。$ | (18) |
则行列两个维度的位移量可表示为：
| $ y_{\text {offset }}=u^{*}-M 。$ | (19) |
| $ x_{\text {offset }}=v^{*}-N 。$ | (20) |
行数位移量为yoffset，yoffset大于0时为下移，小于0时为上移，列数位移量为xoffset，xoffset大于0时为右移，小于0时为左移。假设相关系数矩阵中的第1行第1列为做二维互相关时窗口滑动的初始位置，此时位移层的位移量最大，为上移M－1，左移N－1。为提升二维互相关算法的矫正效果，选择将参考层实时进行更新。
由于是按照高度层对齐的，对齐的过程对于高度向距离徙动的矫正也有着很好的作用。最终对齐效果图如图 7所示，对齐前每一层某点目标P的回波既不在同一个方位点也不在同一个距离门，对齐后P的回波位于同一个距离门和同一个方位点处。
|
图 7 数据对齐效果图 Fig. 7 Illustration of the results after data alignment |
优化后的RDA三维成像是对其二维成像的延伸，由式(3)可推导出雷达的三维空间域距离脉压后的信号：
| $ \begin{gather*} \boldsymbol{S}_{\mathrm{rc}}\left(f_{\mathrm{k}}, t_{\mathrm{m}}, t_{\mathrm{h}}\right)=\sigma_{0} T_{\mathrm{P}} \operatorname{sinc}\left[T_{\mathrm{P}}\left(f_{\mathrm{k}}+K_{\mathrm{r}} \tau^{\prime}\right)\right] \cdot \\ \omega_{\mathrm{a}}\left(t_{\mathrm{m}}\right) \cdot \omega_{\mathrm{h}}\left(t_{\mathrm{h}}\right) \cdot \exp \left(-\mathrm{j} 2 \mathsf{π} f_{\mathrm{c}} \tau^{\prime}\right)。\end{gather*} $ | (21) |
式中：th为高度向时间轴；ωh(·)为高度包络；τ′为雷达到三维空间域中某一点的双程时延，即τ′=2R(tm, th)/c，其中与式(4)不同的是，斜距R需要考虑雷达的高度坐标，即：
| $ \begin{gather*} R\left(t_{\mathrm{m}}, t_{\mathrm{h}}\right)=\sqrt{R_{0}^{2}+V_{\mathrm{r}}^{2} t_{\mathrm{m}}^{2}+V_{\mathrm{r}}^{2} t_{\mathrm{h}}^{2}} \approx \\ R_{0}+\frac{V_{\mathrm{r}}^{2} t_{\mathrm{m}}^{2}}{2 R_{0}}+\frac{V_{\mathrm{r}}^{2} t_{\mathrm{h}}^{2}}{2 R_{0}} 。\end{gather*} $ | (22) |
将式(22)代入式(21)，并对方位向和高度向分别做傅里叶变换后可得：
| $ \begin{gather*} \boldsymbol{S}_{\mathrm{rdh}}\left(f_{\mathrm{k}}, f_{\mathrm{m}}, f_{\mathrm{h}}\right)=\mathrm{AT}_{\mathrm{P}} \operatorname{sinc}\left[T_{\mathrm{P}}\left(f_{\mathrm{k}}+K_{\mathrm{r}} \tau^{\prime}\right)\right] \cdot \\ W_{\mathrm{a}}\left(f_{\mathrm{m}}\right) \cdot W_{\mathrm{h}}\left(f_{\mathrm{h}}\right) \cdot \exp \left(-\mathrm{j} 4 \mathsf{π} \frac{f_{\mathrm{c}} R_{0}}{c}\right) \cdot \\ \exp \left(\mathrm{j} \mathsf{π} \frac{f_{\mathrm{m}}^{2}}{K_{\mathrm{a}}}\right) \cdot \exp \left(\mathrm{j} \mathsf{π} \frac{f_{\mathrm{h}}^{2}}{K_{\mathrm{h}}}\right) 。\end{gather*} $ | (23) |
Wh(fh)为ωh(th)的频域形式，二者在形状上一致。Kh为雷达的高度向调频率，表达式同式(6)，由于本文中雷达方位向速度和上升速度大小相等，故Kh和Ka在数值上是一样的。可分别将x=Vrtm，z=Vrth代入式(23)，并将R0替换为y，式(23)后三个相位可表示为：
| $ \mathit{\varPhi}_{1}=\exp \left(-\mathrm{j} 4 \mathsf{π} \frac{f_{\mathrm{c}} y}{c}\right)，$ | (24) |
| $ \mathit{\varPhi}_{2}=\exp \left(\mathrm{j} 2 \mathsf{π} \frac{x^{2}}{\lambda y}\right) ，$ | (25) |
| $ \mathit{\varPhi}_{3}=\exp \left(\mathrm{j} 2 \mathsf{π} \frac{z^{2}}{\lambda y}\right) 。$ | (26) |
Φ1为各距离门目标固有的相位信息，与图像的强度无关，无需去除。正侧视和小斜视角下高度向和方位向之间的耦合效应比较弱，故可以分别对方位向和高度向进行匹配滤波。Φ2可由式(12)的滤波器进行压缩，根据Φ3的表达式，高度向的滤波器表达式可设计为：
| $ \begin{equation*} \boldsymbol{H}_{\mathrm{h}}=\exp \left(-\mathrm{j} 2 \mathsf{π} \frac{z^{2}}{\lambda y}\right) 。\end{equation*} $ | (27) |
与方位向一致，循环单独对每一个距离门的高度向进行压缩，RDA近场优化三维成像流程如图 8。
|
图 8 RDA近场优化后三维成像流程图 Fig. 8 Flowchart of the 3D imaging process following the optimization of the RDA under near-field conditions |
首先对三维数据进行距离向脉压，然后采用二维互相关算法进行数据错位矫正，由于三维成像数据量较大，为了节省计算空间也为了剔除除目标外其他杂波的影响，在数据错位矫正后限制成像区域截取数据，只对有目标的区域进行成像，对截取完后的数据进行方位向的距离徙动矫正，并将方位向和高度向的傅里叶变换到频域，从成像区域第一个距离门开始，对高度向和方位向分别进行匹配滤波，最终均逆傅里叶变换变回至时域，在遍历完所有距离门后输出三维成像结果，为了提升显示效果，最后选取阈值使无目标区域不显示在图像中。
2 RDA近场二维成像及三维成像实验结果 2.1 实验成像系统构成实验器材如图 9所示，主要由滑轨和雷达板组成，其中滑轨的控制台用于设置程序控制电机转动，两个方向的电机分别带动雷达在X方向和Z方向上移动，X方向即方位向，Z方向即高度向，对于二维成像而言，固定高度向，雷达只在方位向上移动，本文实验部分的雷达板采用的是TI公司生产的IWR1642BOOST毫米波雷达开发板，雷达参数设置如表 2。
|
图 9 实验室成像系统实物图 Fig. 9 Physical image of laboratory imaging system |
|
|
表 2 雷达参数设置 Table 2 Radar parameter settings |
实验目标采用飞机模型，滑轨控制雷达从右往左扫描，如图 2所示，扫描速度为0.04 m/s，由于飞机翼展90 cm左右，扫描长度控制在1 m左右，二维成像实验结果对比图如图 10所示。
|
((a)飞机朝向0°实验场景图；(b)飞机朝向0°时优化前近场成像结果；(c)飞机朝向0°时优化后近场成像结果；(d)飞机朝向45°实验场景图；(e)飞机朝向45°时优化前近场成像结果；(f)飞机朝向45°时优化后近场成像结果；(g)飞机朝向90°实验场景图；(h)飞机朝向90°时优化前近场成像结果；(i)飞机朝向90°时优化后近场成像结果；(j)飞机朝向135°实验场景图；(k)飞机朝向135°时优化前近场成像结果；(l)飞机朝向135°时优化后近场成像结果；(m)飞机朝向180°实验场景图；(n)飞机朝向180°时优化前近场成像结果；(o)飞机朝向180°时优化后近场成像结果。(a)Experimental scene diagram of aircraft facing 0°; (b)Near field imaging results before optimization when the aircraft is facing 0°; (c)Near field imaging results after optimization when the aircraft is facing 0°; (d)Experimental scene diagram of aircraft facing 45°; (e)Near field imaging results before optimization when the aircraft is facing 45°; (f)Near field imaging results after optimization when the aircraft is facing 45°; (g)Experimental scene diagram of aircraft facing 90°; (h)Near field imaging results before optimization when the aircraft is facing 90°; (i)Near field imaging results after optimization when the aircraft is facing 90°; (j)Experimental scene diagram of aircraft facing 135°; (k)Near field imaging results before optimization when the aircraft is facing 135°; (l)Near field imaging results after optimization when the aircraft is facing 135°; (m)Experimental scene diagram of aircraft facing 180°; (n)Near field imaging results before optimization when the aircraft is facing 180°; (o)Near field imaging results after optimization when the aircraft is facing 180°.) 图 10 飞机朝向0°、45°、90°、135°、180°时RDA近场优化前后结果对比图 Fig. 10 Comparison of RDA results before and after optimization for aircraft orientations at 0°, 45°, 90°, 135°, and 180° |
图 10展示了飞机分别朝向0°、45°、90°、135°、180°时RDA近场优化前后的结果对比，可以看出RDA经过近场优化后，压缩效果较优化前好。图 10(b)相比于图 10(c)，飞机模型机翼和机尾部分散焦较为严重，机身部分细节丢失；图 10(e)相比于图 10(f)，机翼部分和机身部分的对比较为明显，图 10(f)明显压缩效果更好；对于图 10(h)和图 10(i)，图 10(h)飞机头部和机翼连接部位压缩效果较差，相比而言图 10(i)回波能量更为集中；对比图 10(k)和图 10(l)，图 10(k)成像效果远不如图 10(l)，机身和飞机头部无法分开，甚至无法分辨出目标的大致形状；最后图 10(n)和图 10(o)相比，图 10(n)的图像质量也是不如图 10(o)，机尾部分细节丢失，飞机头部在方位向上的展宽较为严重，机翼的能量也不集中。
为定量化分析近场优化前后二维成像实验结果，引入图像熵(Image entropy，IE)和图像对比度(Image contrast，IC)两个评估指标[29]，IE代表着图像的混乱程度，图像越聚焦反映在IE上其值越小，IC是作为图像聚焦的度量，IC越大代表着大部分后向散射能量越集中，压缩效果越好。因此，IE越小，IC越大，其成像质量越好。IE和IC计算方法如下：
| $ I E=-\sum\limits_{i=1}^{m} \sum\limits_{j=1}^{n} \boldsymbol{D}(i, j) \cdot \ln [\boldsymbol{D}(i, j)], $ | (28) |
| $ I C=\frac{\mathrm{E}\left(|\boldsymbol{I}|^{2}\right)}{[\mathrm{E}(|\boldsymbol{I}|)]^{2}} 。$ | (29) |
式中：m，n分别为成像矩阵中距离向和方位向像素点数；E(·)为数学期望；·为绝对值；I为像素幅值。D(i, j)为散射强度密度，可用以下表达式来计算：
| $ \begin{equation*} \boldsymbol{D}(i, j)=\frac{|\boldsymbol{I}(i, j)|^{2}}{\sum\limits_{i=1}^{m} \sum\limits_{j=1}^{n}|\boldsymbol{I}(i, j)|^{2}} 。\end{equation*} $ | (30) |
定量分析结果如表 3所示。
|
|
表 3 近场优化前后成像结果定量化对比 Table 3 Quantitative comparison of imaging results before and after optimization |
由表 3可以看出，近场优化后二维成像结果图像熵普遍低于优化前，图像对比度优化后普遍高于优化前，这与前文的分析一致，本文所提出的基于传统RDA的近场优化，其二维成像质量比起传统RDA要高。
2.3 优化后近场三维成像实验结果三维成像实验目标采用上文二维成像实验中的飞机模型，成像场景如图 11所示，飞机放在回波较弱的泡沫柱上，成像区域限制在0~1.5 m之间，最终实验结果如图 12所示。
|
图 11 飞机模型三维成像场景图 Fig. 11 3D imaging scene of aircraft model |
|
图 12 经数据错位矫正后RDA近场优化三维成像结果图 Fig. 12 RDA near-field optimized 3D imaging results after data misalignment correction |
图 13展示了未进行数据错位矫正时的三维成像结果图，可以明显看出矫正后较校正前成像效果要好，图 13相较于图 12，机身部位展宽，压缩效果较图 12差，机头部分中飞机的四个螺旋桨连在一起，远没有图 12成像清晰。现将二维互相关算法应用于三维成像中最常用的RMA，用以对比RDA近场优化后和RMA的三维成像结果，如图 14所示。
|
图 13 未进行数据错位矫正三维成像结果图 Fig. 13 3D imaging results without data misalignment correction |
|
((a)RDA近场优化后三维成像结果图；(b)RMA三维成像结果图。(a)3D imaging results of RDA optimization under near-field conditions; (b)3D imaging results of RMA.) 图 14 RDA近场优化后和RMA的三维成像结果对比 Fig. 14 Comparison of 3D imaging results between RDA optimization under near-field conditions and RMA |
图 14(a)展示的是RDA近场优化后的三维成像结果图，图 14(b)展示的是RMA的三维成像结果图，在选取相同阈值的前提下，可以看出本文所提出的RDA近场优化后的三维成像效果优于RMA，首先是机头和机翼部分，图 14(a)相较于图 14(b)，图 14(a)中飞机的四个螺旋桨与机头能明显的区分开，机翼部分相比于图 14(b)在高度向上更为紧凑，说明RDA在高度向上的聚焦效果更好，其次对于离方位轴最近的机尾部分，RMA的成像效果欠佳，分辨能力较弱，混成一团，相比较而言，图 14(a)中机尾部分的成像效果较好，能较为清楚的看到飞机模型两端的水平尾翼和上方的垂直尾翼。
3 结语由于方位向调频率会随着距离的变化而变化，且其在近场条件下变化的更为显著，传统RDA方位向调频率不变的近似在近场条件下并不适用，成像效果较差。本文基于传统RDA在近场条件下进行优化，并结合仿真和实验与传统RDA的成像结果进行对比，结果显示RDA经近场优化后的成像质量要高于传统RDA，实现了RDA在近场条件下的目标精准成像，对于优化方案中由于循环带来的成像时间成本的增加在近场等数据量较小的场景下可以忽略不计。
近场三维成像方面，在平面扫描模式下对跨越多个距离门的非平面目标进行三维成像时，其回波数据受平台移动误差和高度向距离徙动的影响，会在方位向和距离向均产生不同程度的错位，本文通过引入二维互相关算法计算数据在这两个维度的偏移量来进行数据错位矫正，提升了最终的成像质量，由于RDA在近场优化后所表现出来二维成像的优势，本文也将其拓展至三维成像并在相同的条件下与现今使用最多的RMA来进行比较，可看出本文算法的实验结果其成像目标的聚焦情况更好，RMA之所以是使用最多的算法必然有着其无可替代的优势，其无需距离徙动矫正和利用二维快速傅里叶变换至频域来计算的特点使其计算效率远高于其他算法，且RMA在大斜视角下也能够获得较好的成像效果，但是RMA相比于RDA可能更易受平台运动误差的影响，这类影响仅靠二维互相关算法矫正效果并不佳，本文算法相比于RMA可能更适用于该类成像场景，为毫米波雷达近场多维成像提供了另一种思路。
| [1] |
寇蕾蕾, 郜海阳, 林正健, 等. 星载主动遥感测云现状与展望[J]. 遥感学报, 2023, 27(9): 2041-2059. Kou Leilei, Gao Haiyang, Lin Zhengjian, et al. Current status and prospects of spaceborne active remote sensing for cloud observation[J]. Journal of Remote Sensing, 2023, 27(9): 2041-2059. (  0) |
| [2] |
黄岩, 张慧, 兰吕鸿康, 等. 汽车毫米波雷达信号处理技术综述[J]. 雷达学报, 2023, 12(5): 923-970. Huang Yan, Zhang Hui, Lan Lühongkang, et al. A Review of automotive millimeter-wave radar signal processing techniques[J]. Journal of Radar, 2023, 12(5): 923-970. ( 0) |
| [3] |
阚瀛芝. 毫米波近场隐匿目标三维成像技术[D]. 长沙: 国防科学技术大学, 2017. Kan Yingzhi. 3D Imaging Technology for Concealed Targets in Millimeter-Wave Near Field[D]. Changsha: National University of Defense Technology, 2017. ( 0) |
| [4] |
宋少秋, 邢世其, 刘燊文, 等. 近场毫米波三维成像算法[J]. 太赫兹科学与电子信息学报, 2022, 20(12): 1231-1237. Song Shaoqiu, Xing Shiqi, Liu Shenwen, et al. Near-field millimeter-wave 3d imaging algorithm[J]. Journal of Terahertz Science and Electronic Information Technology, 2022, 20(12): 1231-1237. ( 0) |
| [5] |
涂胜龙. 高分辨车载FMCW SAR成像方法研究[D]. 西安: 西安电子科技大学, 2023. Tu Shenglong. Research on High-Resolution Imaging Methods for Vehicle-Mounted FMCW SAR[D]. Xi'an: Xidian University, 2023. ( 0) |
| [6] |
石星. 毫米波雷达的应用和发展[J]. 电讯技术, 2006(1): 1-9. Shi Xing. Applications anddevelopment of millimeter-wave radar[J]. Telecommunication Engineering, 2006(1): 1-9. ( 0) |
| [7] |
Zeng C Z, Zhu W G, Jia X, et al. Sparse aperture ISAR imaging method based on joint constraints of sparsity and low rank[J]. IEEE Transactions on Geoscience and Remote Sensing, 2020, 59(1): 168-181. ( 0) |
| [8] |
李瓒. 行人毫米波安检成像算法研究[D]. 西安: 西安电子科技大学, 2022. Li Zan. Research on Millimeter-Wave Security Imaging Algorithm for Pedestrians[D]. Xi'an: Xidian University, 2022. ( 0) |
| [9] |
Mercuri M, Lorato R I, Liu Y, et al. Vital-sign monitoring and spatial tracking of multiple people using a contactless radar-based sensor[J]. Nature Electronics, 2019, 2(6): 252-262. DOI:10.1038/s41928-019-0258-6 ( 0) |
| [10] |
赵翔, 王威, 李晨洋, 等. 基于毫米波雷达微动信号和脉搏波数据融合的睡眠呼吸暂停低通气综合征筛查技术[J]. 雷达学报(中英文), 2025, 14(1): 102-116. Zhao Xiang, Wang Wei, Li Chenyang, et al. Screening technology for sleep Apnea-hypopnea syndrome based on fusion of millimeter-wave radar micro-motion signal and pulse wave data[J]. Journal of Radar (Chinese and English Edition), 2025, 14(1): 102-116. ( 0) |
| [11] |
冯利鹏, 王辉, 郑世超, 等. FMCW体制毫米波SAR实时信号处理方法研究[J]. 上海航天(中英文), 2022, 39(3): 116-121. Feng Lipeng, Wang Hui, Zheng Shichao, et al. Research on Real-time signal processing method for Millimeter-wave FMCW SAR system[J]. Shanghai Aerospace (Chinese and English Edition), 2022, 39(3): 116-121. ( 0) |
| [12] |
Essen H, Lorenz F, et al. Millimeter wave radar for runway debris detection[C]//Tyrrhenian International Workshop on Digital Communi cations—Enhanced Surveillance of Aircraft and Vehicles. Capri, Italy: IEEE, 2011: 65-68.
( 0) |
| [13] |
Labor T H, Stuart W R. Three methods for reconstructing objects from x-rays: A comparative study[J]. Computer Graphics and Image Processing, 1973, 2(2): 151-178. DOI:10.1016/0146-664X(73)90025-7 ( 0) |
| [14] |
田雪, 梁兴东, 李焱磊, 等. 基于子孔径包络误差校正的SAR高精度运动补偿方法[J]. 雷达学报, 2014, 3(5): 583-590. Tian Xue, Liang Xingdong, Li Yanlei, et al. High-precision motion compensation method for sar based on sub-aperture envelope error correction[J]. Journal of Radar, 2014, 3(5): 583-590. ( 0) |
| [15] |
Yang Y, Liu J, Wang Y, et al. Automotive SAR imaging by range migration algorithm[C]//Proceedings of the 2022 InternationalConfere-nce on Microwave and Millimeter Wave Technology, Harbin, IEEE, 2022: 1-3.
( 0) |
| [16] |
Cumming I G, Wong F H. Digital Signal Processing of Synthetic Aperture Radar Data: Algorithms and Implementation[M]. Beijing: Publishing House of Electronics Industry, 2004: 140-141.
( 0) |
| [17] |
Wu C. A digital system produce imagery from SAR[C]//Proceedings of the AIAA System Design Driven by Sensors Conference. Washington, DC: American Institute of Aeronautics and Astronautics, 1976: 76-968.
( 0) |
| [18] |
Cumming I G, Bennett J R. Digital processing of SEASAT SAR data[C]//Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing. Washington, DC: IEEE, 1979: 710-718.
( 0) |
| [19] |
Bennett J R, Cumming I G. A digital processor for the production of SEASAT synthetic aperture radar imagery[C]//Proceedings of the 9th SURGE Workshop. Washington, DC: NASA, 1979: SP-154.
( 0) |
| [20] |
Wu C, Liu K Y, Jin M J. A modeling and correlation algorithm for spaceborne SAR signals[J]. IEEE Transactions on Aerospace and Electronic Systems, 1982, 18(5): 563-575. ( 0) |
| [21] |
Barber B C. Theory of digital imaging from orbital synthetic aperture radar[J]. International Journal of Remote Sensing, 1985, 6(6): 1009-1057. ( 0) |
| [22] |
Smith A M. A new approach to range Doppler SAR processing[J]. International Journal of Remote Sensing, 1991, 12(2): 235-251. DOI:10.1080/01431169108929650 ( 0) |
| [23] |
Bamler R, Breit H U, Steinbrecher J, et al. Algorithms for X-SAR processing[J]. IEEE Transactions on Geoscience and Remote Sensing, 1993, 31(4): 1589-1592. ( 0) |
| [24] |
Bennett J R, Cumming I G, Deane R A, et al. SEASAT imagery shows St. Lawrence[J]. Aviation Week and Space Technology, 1979, 111(8): 19. ( 0) |
| [25] |
David M S, Douglas L M, Thomas E H. Three-dimensional millimeter-wave imaging for concealed weapon detection[J]. IEEE Tran-sactions on Microwave Theory and Techniques, 2001, 49(9): 1581-1592. DOI:10.1109/22.942570 ( 0) |
| [26] |
Yanik M E, Torlak M. Near-field 2-D SAR imaging by millimeter-wave radar for concealed item detection[C]//Proceedings of the IEEE Radio and Wireless Symposium (RWS). Orlando, FL, USA: IEEE, 2019: 1-4
( 0) |
| [27] |
保铮, 邢孟道, 王彤. 雷达成像技术[M]. 北京: 电子工业出版社, 2005: 5-9. Bao Zheng, Xing Mengdao, Wang Tong. Radar Imaging Technology[M]. Beijing: Publishing House of Electronics Industry, 2005: 5-9. ( 0) |
| [28] |
Lewis J P. Fast normalized cross-correlation[C]//Proceedings of Vision Interface. Quebec City, Canacla: Canaclian Image Processing and Pattern Recognition Society, 1995: 120-123.
( 0) |
| [29] |
郑易颖, 李郝亮, 刘燊文, 等. 极化圆迹SAR自适应子孔径成像方法[J/OL]. 现代雷达, https://kns.cnki.net/kcms/detail1132.1353.tn.20230104.1005.001.html. Zheng Yiying, Li Haoliang, Liu Shenwen, et al. Adaptive sub-aperture imaging method for polarimetric circular-track SAR[J/OL]. Modern Radar, https://kns.cnki.net/kcms/detail1132.1353.tn.20230104.1005.001.html. ( 0) |
2. Laboratory for Regional Oceanography and Numerical Modeling, Qingdao National Laboratory for Marine Science and Technology, Qingdao 266237, China