基于多载波宽带OCDM的太赫兹汽车雷达

Bhattacharjee S, Mishra K V, Annavajjala R, et al. Multi-carrier wideband OCDM-based THz automotive radar[C]//ICASSP 2023-2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2023: 1-5.

1. 引言与研究背景

自动驾驶是汽车工业的大趋势之一，大多数汽车制造商已经在商用车辆中引入了各种级别的自动驾驶功能。虽然相机、雷达、激光雷达和超声波等多种传感器被用于实现车辆自主性，但雷达因其成本低廉、全天候工作的特性而备受青睐。然而，目前工作在24 GHz和77 GHz的毫米波汽车雷达仅有几GHz的带宽，无法实现光学传感器的高分辨率成像。因此，在太赫兹（THz）频段感知汽车环境正逐渐成为研究热点。

高频工作具有多重优势：系统尺寸更小、可实现更多通道从而获得更高角分辨率、可用带宽更大从而产生更高的距离分辨率。目前，低THz频率如0.15 THz和0.3 THz频段分别提供6 GHz和16 GHz的非授权连续带宽。这些宽频带使汽车雷达能够实现类似激光雷达的成像能力。

然而，在100 GHz以上工作的汽车雷达的主要缺点是由于大气吸收和衰减增加导致的高传播损耗。更重要的是，低THz频谱表现出距离相关的频谱窗口。虽然在1米以下的距离，整个频段可以被视为带宽约为1 THz的单一传输窗口，但在更高距离处，由于分子吸收增加，存在多个宽度为数十或数百GHz的传输窗口。事实上，每个传输窗口的带宽随传输距离而收缩，当距离从1米增加到10米时，由于高吸收峰，带宽减少了一个数量级。因此，在THz频段存在一个关键权衡：在高带宽下操作雷达以提高距离分辨率，与维持足够的最大可探测距离之间的平衡。

2. 系统模型与信道特性

2.1 双基地雷达系统配置

我们考虑一个双基地汽车系统，其中发射（Tx）车辆发射的信号被$P$个感兴趣的目标反射，然后被接收（Rx）车辆捕获。雷达场景包含$P$个遵循Swerling-0目标模型的非起伏点目标。目标位置相对于双基地雷达在Tx信号的照射时间内线性变化：

$$r_p(t) = r_p + v_p t$$

其中$r_p$是$t = 0$时的初始距离，$v_p$是常数径向速度。

2.2 THz信道建模

发射机发出的信号通过THz频段的频率选择性时变信道，其冲激响应表达为：

$$h(t, \tau) = \sum_{p=0}^{P-1} h_p \exp(j2\pi f_{D_p}t)\delta(\tau - \tau_p)$$

这里：

• $h_p$是第$p$个点目标的复散射系数
• $\tau_p = \tau_p^{(1)} + \tau_p^{(2)}$是总时延，由Tx到目标和目标到Rx两段路径组成
• 双基地距离$r_p = c\tau_p$，其中$c$是光速
• 多普勒频移$f_{Dp} = f{Dp}^{(1)} + f{D_p}^{(2)}$由目标的线性运动引起
• 目标速度$vp = c\frac{f{D_p}}{f_c}$，$f_c$是工作频率

THz频段信道的路径损耗特性尤为重要，从第$p$个目标反射造成的接收信号功率路径增益由以下路径损耗表征：

$$PL_{LoS}(f_c, r_p) = \left(\frac{4\pi f_c r_p}{c}\right)^2 e^{k_{abs}(f_c)r_p}$$

其中$k_{abs}$是介质的频率相关吸收系数。这个表达式清楚地显示了THz频段的两个关键特征：自由空间路径损耗和指数形式的分子吸收损耗。

3. OCDM波形设计与调制

3.1 多载波系统架构

[图1描述] ：MCW-OCDM雷达发射机将带宽$B$划分为$K$个子带（不一定相等）进行多路复用。来自载荷数据帧的输入被映射到$K$个独立的OCDM调制器块中。图中显示了系统架构，其中每个调制器处理特定带宽$Bi$的信号，然后通过上变频到相应的中心频率$f{c_i}$并组合输出。

发射波形是OCDM信号，它在相同的时间段和带宽内复用一组啁啾信号。总带宽$B$被划分为$K$个子带，第$i$个子带在中心频率$f_{c_i}$处跨越带宽$B_i$。来自载荷数据帧的输入被映射到$K$个独立的OCDM调制器块。

3.2 OCDM信号数学表示

第$i$个子带中的Tx帧由$N_i$个时域符号组成，通过使用数据比特调制$M_i$个子啁啾的相位和幅度获得，占据总带宽$B_i = M_i\Delta f_i$，其中$\Delta f_i$是每个啁啾的带宽，$M_i$是偶数正整数。基带Tx OCDM信号在第$i$个调制器输出处为：

$$\mathbf{S}_i = \boldsymbol{\Phi}_{M_i}^H \mathbf{X}_i$$

其中：

$$\boldsymbol{\Phi}_{M_i} = \frac{1}{\sqrt{M_i}} \boldsymbol{\Theta}_1 \mathbf{F}_{M_i} \boldsymbol{\Theta}_2$$

$$\boldsymbol{\Theta}_1 = \text{diag}\{\Theta_{1,0}; \cdots; \Theta_{1,M_i-1}\}, \quad \Theta_{1,u} = e^{-j\frac{\pi}{4}} e^{j\frac{\pi}{M_i}u^2}$$

$$\boldsymbol{\Theta}_2 = \text{diag}\{\Theta_{2,0}; \cdots; \Theta_{2,M_i-1}\}, \quad \Theta_{2,v} = e^{j\frac{\pi}{M_i}v^2}$$

这里$\boldsymbol{\Phi}_{M_i}^H \in \mathbb{C}^{M_i \times M_i}$表示阶数为$M_i$的逆离散Fresnel变换（IDFnT），$\mathbf{X}_i \in \mathbb{C}^{M_i \times N_i}$是数据符号矩阵。

3.3 循环矩阵性质与特征分解

由于矩阵$\boldsymbol{\Phi}_{M_i}$是循环矩阵，利用特征分解性质，OCDM信号可以重写为：

$$\mathbf{S}_i = \mathbf{F}_{M_i}^H \boldsymbol{\Gamma}^H \mathbf{F}_{M_i} \mathbf{X}_i = \mathbf{F}_{M_i}^H \mathbf{Z}_i$$

其中$\boldsymbol{\Gamma}^H = \mathbf{F}_{Mi} \boldsymbol{\Phi}{Mi}^H \mathbf{F}{Mi}^H$，矩阵$\mathbf{G} \triangleq \boldsymbol{\Gamma}^H \mathbf{F}{M_i} \in \mathbb{C}^{M_i \times M_i}$将输入数据符号$\mathbf{X}_i$变换为缩放的频域符号$\mathbf{Z}_i \triangleq \mathbf{G}\mathbf{X}_i$。$\boldsymbol{\Gamma} \in \mathbb{C}^{M_i \times Mi}$是对角矩阵，其第$m$个对角元素$\Gamma(m)$是$\boldsymbol{\Phi}{M_i}$的第$m$个特征值，对应于根Zadoff-Chu序列：

$$\Gamma(m) = e^{-j\frac{\pi}{M_i}m^2}, \quad \forall m, \quad M_i \equiv 0 \pmod{2}$$

时域OCDM信号可以写为：

$$s_i(t) = \sum_{n=0}^{N_i-1} \sum_{m=0}^{M_i-1} [X_i]_{m,n} e^{j\frac{\pi}{4}} e^{-j\pi M_i \frac{\left(t-nT_i-\frac{mT_i}{M_i}\right)^2}{T_i^2}} e^{j2\pi f_{c_i}t} \text{rect}(t-nT_i)$$

其中$\text{rect}(t) \triangleq \begin{cases} 1 & 0 \leq t \leq T \ 0 & \text{otherwise} \end{cases}$，$T_i = \frac{1}{\Delta f_i}$是OCDM符号持续时间。

4. 雷达信号处理与目标参数估计

4.1 接收信号模型

[图2描述] ：接收处理需要通过为每个子带分配最优权重来组合估计值。图中展示了$K$个解调器分别处理中心频率$f_{c1}$到$f{cK}$的信号，每个解调器输出进入相应的感知处理器，产生距离和速度估计$\hat{r}{pi}, \hat{v}_{pi}$，最后通过最优组合模块得到最终估计$\hat{r}_p, \hat{v}_p$。

Rx在频率$f_{c_i}$处接收到经过双扩展THz雷达信道的雷达回波，表征为目标的时延和多普勒频移的和：

$$y_i^{rad}(t) = \sqrt{PL_{LoS}(f_{c_i}, r_p)} \sum_{p=1}^P h_p s_i(t-\tau_p) e^{j2\pi f_{c_i}\vartheta_p(t-\tau_p)} + w_i(t)$$

其中$\vartheta_p = \frac{v_p}{c}$是归一化速度，$w_i(t) \sim \mathcal{CN}(0, \epsilon_i^2)$表示加性高斯白噪声。

4.2 离散时间处理

经过下变频后，信号在$t = nT_i + m\frac{T_i}{M_i}$处采样得到：

$$[Y_i^{rad}]_{m,n} = \sum_{p=1}^P \tilde{h}_{pi} e^{j2\pi\frac{f_{c_i}}{c}v_p\left(nT_i+m\frac{T_i}{M_i}\right)} \sum_{m'=0}^{M_i-1} [X_i]_{n,m'} e^{j\frac{\pi}{4}} \times e^{-j\pi\frac{M_i}{T_i^2}\left[\left(m-m'\right)\frac{T_i}{M_i}-\tau_p\right]^2} + [W_i]_{m,n}$$

其中$\tilde{h}{pi} = \sqrt{PL{LoS}(f_{c_i}, r_p)} h_p$。

接下来应用DFnT观察跨啁啾的雷达回波：

$$[Y_i^{rad}]_{m,n} = \frac{1}{M_i} \sum_{l=0}^{M_i-1} [Y_i^{rad}]_{l,n} \exp\left(-j\frac{\pi}{4}\right) \exp\left[j\frac{\pi}{M_i}(m-l)^2\right]$$

由于假设最大多普勒频移远小于子载波间隔（$f_{D_i}^{max} \ll \Delta f_i$），近似可得：

$$[Y_i^{rad}]_{m,n} \approx \sum_{p=1}^P [X_i]_{m,n} \tilde{h}_{pi} e^{-j\pi M_i\left(\frac{\tau_p}{T_i}\right)^2} e^{j2\pi(n\vartheta_p f_{c_i}T_i - m\tau_p\Delta f_i)} + [W_i]_{m,n}$$

4.3 最大似然估计

由于载荷数据$\mathbf{X}_i$在雷达接收器已知，通过逐元素除法去除它们，得到雷达观测值：

$$[Z_i^{rad}]_{m,n} = \sum_{p=1}^P \tilde{h}_{pi} e^{j2\pi(n\vartheta_p f_{c_i}T_i - m\tau_p\Delta f_i)} + [W_i]_{m,n}$$

假设目标数量已通过假设检验确定，我们寻找目标参数$\thetai = [\theta{1i}, \cdots, \theta{Pi}]^T$的最大似然估计，其中$\theta{pi} = (r{pi}, v{pi})$。简化的对数似然函数为：

$$\mathcal{L}(Z_i^{rad}; \theta_{pi}) = 2\tilde{h}_{pi}\mathfrak{R}\left[\sum_{m,n}[Z_i^{rad}]_{m,n}e^{-j2\pi n\vartheta_p f_{c_i}T_i}e^{j2\pi m\tau_p\Delta f_i}\right] - \tilde{h}_{pi}^2$$

这是一个二维复周期图。通过量化频率并沿所需维度进行FFT，离散化的对数似然函数为：

$$\mathcal{L}_q(Z_i^{rad}; m', n') = 2\tilde{h}_{pi}\mathfrak{R}\left[\sum_{m=0}^{M_{Per}-1}\sum_{n=0}^{N_{Per}-1}[Z_i^{rad}]_{m,n}e^{-j2\pi\frac{nn'}{N_{Per}}}e^{j2\pi\frac{mm'}{M_{Per}}}\right]$$

其中$m' := \tau_p\Delta f_i$和$n' := \varthetap f{c_i}Ti$是离散化频率，搜索网格为$m' = 0, \cdots, M{Per}-1$和$n' = -\frac{N{Per}}{2}, \cdots, \frac{N{Per}}{2}-1$，且$M_{Per} > Mi$，$N{Per} > N_i$实现过采样。

5. 最优加权组合方案

5.1 组合估计理论

在我们的最优加权组合（OWC）方案中，来自$K$个感知处理器的估计需要组合以获得目标距离和速度的最终估计$\theta_p = (r_p, v_p)$。目标在解调器块$i$处的估计参数可以在小误差假设下通过一阶泰勒级数近似线性化：

$$\hat{\zeta} = \mathbf{1}_K\zeta + \mathbf{e}_\zeta \in \mathbb{C}^K$$

其中$\hat{\zeta} = [\hat{\zeta}_1, \cdots, \hat{\zeta}_K]^T$，$\zeta_i \in {\hat{r}_i, \hat{v}i}$是SP $i$处的估计参数，$\zeta \in {r, v}$是真实目标参数，$\mathbf{e}\zeta = [e_{\zeta1}, \cdots, e{\zeta_K}]^T$是估计误差。

5.2 克拉美-罗下界分析

在高信噪比下，基于DFT的估计器达到克拉美-罗下界（CRLB）。因此，可以将高信噪比下$e_\zeta$的方差近似为：

$$\sigma_{r_i}^2 \approx \frac{6\epsilon_i^2}{(2\pi)^2 M_i N_i(N_i^2-1)|\tilde{h}_i|^2 P_{avg}}\left(\frac{c}{\Delta f_i}\right)^2$$

$$\sigma_{v_i}^2 \approx \frac{6\epsilon_i^2}{(2\pi)^2 M_i N_i(M_i^2-1)|\tilde{h}_i|^2 P_{avg}}\left(\frac{c}{T_i f_{c_i}}\right)^2$$

5.3 最优权重推导

命题1：定义$T(\hat{\zeta}) = \boldsymbol{\beta}^T\hat{\zeta}$；则$T(\hat{\zeta})$是估计$\zeta$的充分统计量，其中$\boldsymbol{\beta} = [\beta_1, \cdots, \beta_K]^T \in \mathbb{R}^K$表示线性组合器的权重。

定理2：最优组合方案是对$K$个感知处理器得到的估计值进行线性加权组合，使估计误差最小化：

$$\min_{\boldsymbol{\beta}} \mathbb{E}\left[|\boldsymbol{\beta}^H\hat{\boldsymbol{\zeta}} - \zeta|^2\right] \quad \text{s.t.} \quad \mathbf{1}_K^H\boldsymbol{\beta} = 1$$

最优组合权重为：

$$\beta_k = \frac{\sigma_{\zeta_k}^{-2}}{\sum_{i=1}^K \sigma_{\zeta_i}^{-2}}, \quad \forall k \in \{1, \cdots, K\}$$

6. 数值实验结果分析

6.1 实验参数设置

在所有实验中，总带宽$B$设置为1.4 THz。基于可用的距离相关传输窗口，将总带宽划分为子带$B_i$，每个跨越1 GHz。每个子带考虑相等数量的啁啾和OCDM符号，即$M = 256$和$N = 256$。因此，每个子带的两个啁啾之间的间隔$\Delta f_i$为3.9 MHz，这完全在THz信道的相干带宽内。OCDM帧时间$T$为0.25 μs。将参考目标放置在距雷达收发器不同距离处，速度为23 m/s。

6.2 实验结果分析

[图3(a)描述] ：频率相关THz路径损耗对距离估计的影响与信噪比的关系（$r = 0.1$ m）。图中比较了三条曲线：SP-K（最后一个子带）、SP-1（第一个子带）和OWC（最优加权组合）。可以看到，在RMSE为$10^{-3}$时，OWC相比单个SP估计有约6 dB的性能提升。OWC曲线始终低于其他两条，表明组合估计的优越性。

[图3(b)描述] ：距离和频率相关THz路径损耗对距离估计的影响与信噪比的关系。图中展示了不同目标距离（0.1 m、1 m、10 m、50 m）下的性能。可以观察到，随着距离增加，所有曲线都向上移动（RMSE增大），但OWC始终保持最佳性能。在较近距离（0.1 m）时，OWC能达到亚毫米级精度。

[图3(c)描述] ：速度估计与信噪比的关系，展示了不同双基地目标距离$r$的影响。图中包含了$r = 50$ m、10 m、1 m、0.1 m等多个距离的曲线，以及OWC和无组合（No combining）的对比。OWC在所有距离下都显著优于单独估计，改善约6 dB。随着目标距离增加，可用的带宽窗口数量减少，路径损耗变得更严重，导致速度估计误差增加。

实验结果表明，我们提出的OWC方案实现了亚毫米级的感知精度，在合理的接收信噪比水平（> 0 dB）下表现优异。通过分配最优权重，组合器能够有效应对THz频段的路径损耗，比单个感知处理器的估计性能显著提升。这种改进在较低距离时更为明显，因为此时有更多的感知处理器可用，能够实现更好的目标定位。

7. 结论

本文提出了一个创新的多载波宽带OCDM框架，成功解决了THz频段雷达目标参数估计中频率和距离相关路径损耗的挑战。我们开发了一种新颖的多级感知算法，有效利用来自不同THz子带的数据帧。基于参数估计的克拉美-罗下界，推导了跨不同感知处理器的最优组合权重。

数值实验充分证明，使用最优加权组合器处理来自不同THz传输窗口的雷达回波，通过优先考虑经历较低吸收损耗的子带的更准确估计，能够显著提高估计精度。相比单个传输窗口的处理，我们的系统实现了三个数量级的距离估计精度改进，达到了亚毫米级的精度。此外，在不同目标距离下，速度估计也得到了显著改善。

附录：数学推导

A. 离散Fresnel变换的循环性质证明

离散Fresnel变换矩阵$\boldsymbol{\Phi}_{M_i}$的定义为：

$$\boldsymbol{\Phi}_{M_i} = \frac{1}{\sqrt{M_i}} \boldsymbol{\Theta}_1 \mathbf{F}_{M_i} \boldsymbol{\Theta}_2$$

其中$\mathbf{F}_{M_i}$是标准DFT矩阵，其元素为：

$$[\mathbf{F}_{M_i}]_{u,v} = \frac{1}{\sqrt{M_i}} e^{j\frac{2\pi}{M_i}uv}$$

相位矩阵定义为：

$$\boldsymbol{\Theta}_1 = \text{diag}\left\{e^{-j\frac{\pi}{4}} e^{j\frac{\pi}{M_i}u^2}\right\}_{u=0}^{M_i-1}$$

$$\boldsymbol{\Theta}_2 = \text{diag}\left\{e^{j\frac{\pi}{M_i}v^2}\right\}_{v=0}^{M_i-1}$$

由于$\boldsymbol{\Phi}_{M_i}$是循环矩阵，其特征值可以通过DFT对角化得到。设$\boldsymbol{\Gamma}$为特征值对角矩阵，则：

$$\boldsymbol{\Phi}_{M_i} = \mathbf{F}_{M_i}^H \boldsymbol{\Gamma} \mathbf{F}_{M_i}$$

通过直接计算可得第$m$个特征值：

$$\Gamma(m) = e^{-j\frac{\pi}{M_i}m^2}, \quad m = 0, 1, \cdots, M_i-1$$

这对应于根Zadoff-Chu序列的特征结构。

B. 最优权重的拉格朗日乘数法推导

优化问题为：

$$\min_{\boldsymbol{\beta}} \mathbb{E}\left[|\boldsymbol{\beta}^H\hat{\boldsymbol{\zeta}} - \zeta|^2\right] \quad \text{s.t.} \quad \mathbf{1}_K^H\boldsymbol{\beta} = 1$$

将目标函数展开：

$$\mathbb{E}\left[|\boldsymbol{\beta}^H\hat{\boldsymbol{\zeta}} - \zeta|^2\right] = \mathbb{E}\left[|\boldsymbol{\beta}^H(\mathbf{1}_K\zeta + \mathbf{e}_\zeta) - \zeta|^2\right]$$

由于约束条件$\mathbf{1}_K^H\boldsymbol{\beta} = 1$，上式简化为：

$$= \mathbb{E}\left[|\boldsymbol{\beta}^H\mathbf{e}_\zeta|^2\right] = \boldsymbol{\beta}^H\mathbf{R}_{e_\zeta}\boldsymbol{\beta}$$

其中$\mathbf{R}{e\zeta} = \mathbb{E}[\mathbf{e}\zeta\mathbf{e}\zeta^H]$是误差协方差矩阵，其对角元素为$\sigma_{\zeta_i}^2$。

构造拉格朗日函数：

$$L = \boldsymbol{\beta}^H\mathbf{R}_{e_\zeta}\boldsymbol{\beta} + \lambda(\mathbf{1}_K^H\boldsymbol{\beta} - 1)$$

对$\boldsymbol{\beta}$求导并令其为零：

$$\frac{\partial L}{\partial \boldsymbol{\beta}} = 2\mathbf{R}_{e_\zeta}\boldsymbol{\beta} + \lambda\mathbf{1}_K = 0$$

解得：

$$\boldsymbol{\beta} = -\frac{\lambda}{2}\mathbf{R}_{e_\zeta}^{-1}\mathbf{1}_K$$

代入约束条件：

$$\mathbf{1}_K^H\boldsymbol{\beta} = -\frac{\lambda}{2}\mathbf{1}_K^H\mathbf{R}_{e_\zeta}^{-1}\mathbf{1}_K = 1$$

因此：

$$\lambda = -\frac{2}{\mathbf{1}_K^H\mathbf{R}_{e_\zeta}^{-1}\mathbf{1}_K}$$

最终得到最优权重：

$$\boldsymbol{\beta} = \frac{\mathbf{R}_{e_\zeta}^{-1}\mathbf{1}_K}{\mathbf{1}_K^H\mathbf{R}_{e_\zeta}^{-1}\mathbf{1}_K}$$

由于$\mathbf{R}{e\zeta}$是对角矩阵，第$k$个权重为：

$$\beta_k = \frac{\sigma_{\zeta_k}^{-2}}{\sum_{i=1}^K \sigma_{\zeta_i}^{-2}}$$

C. 高信噪比下的CRLB推导

对于参数$\theta = [\tau, \vartheta]^T$的Fisher信息矩阵为：

$$\mathbf{J}(\theta) = \mathbb{E}\left[\left(\frac{\partial \ln p(\mathbf{Z};\theta)}{\partial \theta}\right)\left(\frac{\partial \ln p(\mathbf{Z};\theta)}{\partial \theta}\right)^T\right]$$

在高信噪比下，对于时延估计：

$$J_{\tau\tau} = \frac{2|\tilde{h}|^2}{\epsilon^2} \sum_{m,n} \left(\frac{\partial \phi_{m,n}}{\partial \tau}\right)^2$$

其中$\phi_{m,n} = 2\pi(n\vartheta f_c T - m\tau\Delta f)$。

计算导数：

$$\frac{\partial \phi_{m,n}}{\partial \tau} = -2\pi m\Delta f$$

因此：

$$J_{\tau\tau} = \frac{2|\tilde{h}|^2}{\epsilon^2} \sum_{m=0}^{M-1} \sum_{n=0}^{N-1} (2\pi m\Delta f)^2$$

$$= \frac{2|\tilde{h}|^2}{\epsilon^2} N \cdot (2\pi\Delta f)^2 \sum_{m=0}^{M-1} m^2$$

利用求和公式$\sum_{m=0}^{M-1} m^2 = \frac{M(M-1)(2M-1)}{6}$，在$M \gg 1$时近似为$\frac{M^3}{3}$：

$$J_{\tau\tau} \approx \frac{2|\tilde{h}|^2 N (2\pi\Delta f)^2 M^3}{3\epsilon^2}$$

CRLB为Fisher信息的倒数：

$$\sigma_\tau^2 \geq \frac{1}{J_{\tau\tau}} = \frac{3\epsilon^2}{2|\tilde{h}|^2 N (2\pi\Delta f)^2 M^3}$$

转换到距离域（$r = c\tau$）：

$$\sigma_r^2 = c^2\sigma_\tau^2 \approx \frac{6\epsilon^2}{(2\pi)^2 MN(N^2-1)|\tilde{h}|^2 P_{avg}}\left(\frac{c}{\Delta f}\right)^2$$

类似地可推导速度估计的CRLB。

D. 二维周期图的FFT实现

二维复周期图：

$$\mathcal{L}(Z_i^{rad}; \theta_{pi}) = 2\tilde{h}_{pi}\mathfrak{R}\left[\sum_{m,n}[Z_i^{rad}]_{m,n}e^{-j2\pi n\vartheta_p f_{c_i}T_i}e^{j2\pi m\tau_p\Delta f_i}\right] - \tilde{h}_{pi}^2$$

令$\omega_1 = \tau_p\Delta f_i$，$\omega_2 = \varthetap f{c_i}T_i$，则：

$$\mathcal{L}(Z_i^{rad}; \omega_1, \omega_2) = 2\tilde{h}_{pi}\mathfrak{R}\left[\sum_{m,n}[Z_i^{rad}]_{m,n}e^{j2\pi(m\omega_1 - n\omega_2)}\right] - \tilde{h}_{pi}^2$$

这正是二维离散傅里叶变换的形式。通过在频率域量化并使用FFT算法，可以高效计算：

$$\hat{\mathcal{L}}_q(m', n') = \text{FFT2D}\{[Z_i^{rad}]_{m,n}\}$$

其中$(m', n')$是离散频率索引，通过过采样因子$\frac{M_{Per}}{Mi}$和$\frac{N{Per}}{N_i}$提高频率分辨率。

峰值检测得到：

$$[\hat{m}', \hat{n}'] = \arg\max_{m',n'} |\hat{\mathcal{L}}_q(m', n')|$$

从而恢复目标参数：

$$\hat{\tau}_{pi} = \frac{\hat{m}'}{M_{Per}\Delta f_i}, \quad \hat{\vartheta}_{pi} = \frac{\hat{n}'}{N_{Per}f_{c_i}T_i}$$

最终得到距离和速度估计：

$$\hat{r}_{pi} = c\hat{\tau}_{pi}, \quad \hat{v}_{pi} = c\hat{\vartheta}_{pi}$$