【翻译】被动雷达信号处理

这是一篇翻译，为了我自己看着方便。原文在此。

杂波消除

这是我在系列文章中关于被动雷达信号处理算法的第一篇。如果你还没有看过，先去主被动雷达页面了解一下这个项目是关于什么的。这篇文章比上一篇更数学化一些，但我会尽量保持直观的解释。

被动雷达的一个主要挑战是从噪声中提取微弱的目标回波。这种噪声包括直接路径信号泄漏到观测通道以及从静止物体（如山脉和建筑物）的强反射信号。在雷达术语中，这些信号统称为“杂波”。对于被动雷达，杂波信号通常比目标回波强数十倍或数百倍。因此，需要一种杂波消除算法来去除这些信号，以便能够观察到目标回波。

在深入推导之前，这里有一个概念验证的快速示例。下面的动画展示了如何使用杂波消除算法来提取原本完全被噪声淹没的目标信号（使用真实的被动雷达数据）。应用该算法可使噪声降低近40dB！从噪声中出现的尖峰就是目标回波。

为了推导杂波消除算法，需要对被动雷达接收到的信号性质做一些假设。我们将参考通道接收到的信号表示为 $s_r[n]$ ，观测通道接收到的信号表示为 $s_o[n]$ 。接下来，我们可以将观测通道信号建模为三个部分的总和：所需的目标回波信号 $s_{tar}[n]$ 、杂波信号 $s_{cl}[n]$ 和附加噪声项 $\nu[n]$ 。

\[ \begin{equation} \displaystyle s_{o}[n] = s_{tar}[n] + s_{cl}[n] + \nu[n] \end{equation} \]

为了简化起见，我们假设参考信号是发射信号的完美副本。然后，由于杂波信号由发射信号的直接路径副本和来自静止物体的各种多路径反射组成，因此很自然地将其建模为 $ M $ 个缩放和延迟的参考信号副本的总和。数字 $ M $ 被选择为具有显著杂波返回的最大延迟，这会根据具体情况而变化。

\[ \begin{equation} \displaystyle s_{cl}[n] = \sum_{k=0}^{M – 1} b_k s_{r}[n-k] \end{equation} \]

需要注意的是，我们还不知道杂波信号中的系数 $ b_k $ ，因为它们取决于雷达环境中存在的散射物体。然而，如果我们能够使用一些额外信息来找到杂波系数的适当估计值 $\hat{b}_k$ ，则可以利用它们来重建杂波信号的估计值。然后，可以从观测通道信号中减去估计的杂波信号，以获得目标回波信号的估计值。

\[ \begin{equation} \displaystyle \hat{s}_{tar}[n] = s_o[n] – \hat{s}_{cl}[n] = s_o[n] – \sum_{k=0}^{M – 1} \hat{b}_k s_{r}[n-k] \end{equation} \]

因此，杂波消除算法的主要目标可以通过找到杂波系数的估计值来实现。这是一个自适应滤波问题。系数是根据输入数据选择的，以便获得所需的输出。下面的框图展示了这一过程。

那么，自适应滤波算法如何确定最佳滤波系数呢？乍一看，这似乎是不可能的，因为滤波器看起来像是在试图预测自己的输出。为了解释这一点，我们将从构建目标回波的信号模型开始。

类似于将杂波信号建模为参考信号的函数，来自 $ N $ 个不同目标的回波信号可以表示为 $ N $ 个缩放、延迟和多普勒频移的参考信号副本的总和。这里， $ w_k $ 、 $ d_k $ 和 $ f_k $ 分别是第 $ k $ 个目标回波的幅度、延迟和多普勒频移， $ F_s $ 是采样率。

\[ \begin{equation} \displaystyle s_{tar}[n] = \sum_{k=0}^{N-1} w_{k} s_{r}[n-d_k] \exp \left(\frac{j2\pi n f_k}{F_s} \right) \end{equation} \]

这个方程看起来很像杂波信号的方程，只是现在总和中的每个元素都乘以了一个多普勒频移项 $\exp(j2\pi n f_k / F_s)$ 。这给了我们一个关于如何将回波信号与杂波信号分开的提示。首先，回忆不同频率的复指数函数彼此正交。这意味着由于多普勒频移，只要多普勒频移不为零且信号在足够长的时间间隔内测量，目标回波就与参考信号线性无关。相反，我们假设杂波信号是延迟参考信号副本的线性组合。换句话说，如果我们能够找到观测通道信号中与延迟参考信号线性相关的所有部分并减去它们，剩下的将只是回波信号。

用向量空间的概念来重新表述这个想法是很方便的。由于杂波信号由缩放和延迟的参考信号副本组成，因此可以说它占据了一个向量空间，其基向量是延迟的参考信号副本。我们将这个空间称为杂波空间。杂波系数 $ b_k $ 是杂波信号在杂波空间中的坐标。回波信号与每个延迟的参考信号副本线性无关，因此它与杂波空间正交。我们可以通过将观测通道信号投影到杂波空间来找到杂波信号。减去杂波信号后，剩下的只有回波信号。

在被动雷达中可以使用多种不同的自适应滤波算法进行杂波消除。我们使用了一种称为最小二乘滤波器（有时也称为维纳滤波器）的算法。最小二乘滤波器是一种块算法，因为它涉及将输入数据分成许多样本长度的块，并为每个块计算一组滤波器系数。这与基于随机梯度的方法（如最小均方滤波器）形成对比，在后一种方法中，滤波器抽头是使用每个输入样本逐一更新的。这篇论文比较了这两种自适应滤波方案以及几种其他方案在被动雷达杂波消除中的性能。如果你还想了解更多关于自适应滤波的内容，我还强烈推荐西蒙·海金的教科书《自适应滤波理论》。

最小二乘算法的一个好处是，鉴于上述杂波消除过程的向量空间表述，其推导是直观的。首先，设 $ {x_k} $ 是一个向量，包含延迟 $ k $ 个样本的参考信号的 $ L $ 个样本，即 $ s_r[n – k] $ 。对于 $ k $ 在0到 $ M $ 之间的向量 $ {x_k} $ 是杂波空间的基向量。然后，我们可以定义一个矩阵 $ {X} $ ，这是一个 $L\times M$ 的矩阵，其列是基向量 $ {x_k} $ （ $ k $ 在0到 $ M $ 之间）。

接下来，设 $ {} $ 是一个长度为 $ M $ 的任意向量。注意， $ {X} $ 与 $ {} $ 的乘积是一个长度为 $ L $ 的向量。 $ {X} $ 将杂波空间中的任意点（一个长度为 $ M $ 的向量）转换为其在信号空间中的表示（一个长度为 $ L $ 的向量）。然后，设 $ {y} $ 是一个包含观测信号 $ s_o[n] $ 的 $ L $ 个样本的向量。最后，对于任意杂波向量 ${\beta}$ ，设 ${\epsilon}$ 表示 $ {y} $ 和 ${X}{\beta}$ 之间的残差。

\[ \begin{equation} \displaystyle {\epsilon} = {y} – {X}{\beta} \end{equation} \]

或者等价地，

\[ \begin{equation} \displaystyle {y} = {X} {\beta} + {\epsilon} \end{equation} \]

将这个方程与观测通道信号的原始信号模型进行比较，我们看到，如果我们选择正确的杂波向量 $ {}$ ，那么项 ${X} {\beta}$ 就代表杂波信号 $ s_{cl} $ 。 ${\epsilon}$ 则代表剩下的两个项，即 $s_{tar}+\nu$ 。

那么，我们如何找到这个最优值 ${\beta}$ （我们将称之为 ${\hat{\beta}}$ ）呢？我们希望项 ${X} {\beta}$ 能够捕获 ${y}$ 的最大部分，使得在减去它之后剩下的部分最小。因此，最优的杂波向量 ${\hat{\beta}}$ 可以通过最小化 ${\epsilon}$ 的平方模值得到。对应的最小模长向量 ${\hat{\epsilon}}$ 就成为我们对目标信号的最佳估计。由于噪声项 $\nu$ 与参考信号不相关，因此无法去除，但只要它相对于目标信号很小，就可以忽略不计。

这令人兴奋，因为它意味着我们已经将问题简化为普通的最小二乘回归。这种类型问题的解决方案是众所周知的，由下式给出

\[ \begin{equation} \displaystyle {\hat{\beta}} = {(X^H X)^{-1} X^H y} \end{equation} \]

注意，在前面的方程中， $ {X^H} $ 表示 $ {X} $ 的共轭转置，因为 $ {X} $ 的元素可以是复数。

在实践中，矩阵 $ {X^H X} $ 有时是病态的（接近奇异），这意味着最小二乘问题的解不是唯一确定的。一种处理方法是使用正则化技术，如提霍诺夫正则化（即 $ L2 $ 正则化，译者注）。使用提霍诺夫正则化的杂波系数向量 $ {} $ 的修改方程如下，其中 $ {} $ 是一个 $ M M $ 的矩阵，称为提霍诺夫矩阵。

\[ \begin{equation} \displaystyle {\hat{\beta}} = ({X^H X} + {\Gamma} {^H} {\Gamma})^{-1}{X^H y} \end{equation} \]

对提霍诺夫矩阵的不同选择会导致不同类型的正则化。一个常见的选择是 ${\Gamma}=\alpha{I}$ ，其中 ${I}$ 是单位矩阵， $$ 是一个常数。这种选择（相当于 $L2 $ 正则化）对解向量的平方模长添加了惩罚，以选择模长较小的解。

一旦获得了杂波系数向量，恢复目标回波信号就变得简单了：

\[ \begin{equation} \hat={y}-{X}{\hat{\beta}} \end{equation} \]

至此，杂波消除过程就完成了！继续阅读这个系列的下一篇，看看回波信号是如何被处理以揭示目标特征的。

距离多普勒

这是我在系列文章中关于被动雷达信号处理算法的第二篇。第一部分在这里。如果你还没有看过，先去主被动雷达页面了解一下这个项目是关于什么的。

在上一篇文章中，我们推导了一种算法，用于从被动雷达的观测通道中去除雷达杂波。理论上，这应该只剩下我们试图分析的目标回波。这些回波能告诉我们关于它们所反射的目标的哪些信息呢？

这里感兴趣的目标特性是每个目标回波的时间延迟和多普勒频移。这些分别与目标的距离和速度有关。这些量之间的精确关系取决于接收机和发射机的几何位置以及目标的位置和方位 —— 目前我们只关注从回波信号中提取延迟和多普勒频移。

为了更具体地说明这一点，考虑我们在上一篇文章中用于回波信号 $ s_{tar}[n] $ 的信号模型。它由参考信号 $ s_r[n] $ 的缩放、延迟和多普勒频移副本组成。这里， $ w_k $ 、 $ d_k $ 和 $ f_k $ 分别是第 $ k $ 个目标回波的幅度、延迟和多普勒频移， $ F_s $ 是采样率。

\[ s_{tar}[n] = \sum_{k=0}^{N-1} w_{k} s_{r}[n-d_k] \exp \left(\frac{j2\pi n f_k}{F_s} \right) \]

目标检测相当于找出回波信号中存在的延迟值和多普勒频移。

这可以通过计算互模糊函数 $\chi_{s_r, s_{tar}}(\tau, f)$ 来实现。可以将互模糊函数视为互相关函数的扩展。对两个相似信号进行互相关可以确定它们之间的延迟。互模糊函数类似，但它还能确定信号之间的频率差异。这是通过利用不同频率的正弦波彼此正交的特性来实现的。

\[ \chi_{s_1, s_2}(\tau, f) = \sum_{n=0}^{N-1} s_1[n]s_2^*[n – \tau]\exp \left(\frac{-j 2 \pi n f}{F_s} \right) \]

如果在一系列延迟和频率偏移值的范围内计算互模糊函数，结果将是一个二维曲面，称为距离 - 多普勒图。目标将出现在该曲面上的 $ $ 和 $ f $ 值处的尖峰，这些值对应于目标回波的延迟和多普勒频移。

值 $ N $ 称为相干处理区间（CPI），必须谨慎选择，因为它限制了系统的多普勒分辨率。一般来说，选择较大的 CPI 意味着可以区分频率差异较小的信号（即多普勒分辨率与 CPI 成反比）。例如，1 秒的 CPI 对应最小多普勒分辨率为 1Hz。

直接计算互模糊函数计算量相当大，特别是如果需要大量的距离和多普勒频段时。幸运的是，有一种近似算法可以显著加快计算速度。其基本思想是使用频域方法，这样可以利用快速傅里叶变换（FFT）算法。在进行 FFT 之前，我们对输入信号进行抽取，因为我们感兴趣的多普勒频移范围远小于输入信号的带宽（小于 500Hz，而输入信号带宽为 200kHz）。

为了推导快速互模糊算法，我们可以重新索引求和式，将其分成两个部分，其中 $ N = ML $ ：

\[ \chi(\tau, f) = \sum_{k=0}^{M-1} \left[ \sum_{l=0}^{L-1} s_1(l + kM)s_2^*(l + kM – \tau)\exp \left(\frac{j 2 \pi l f}{F_s}\right) \right] \exp \left(\frac{-j 2 \pi kM f}{F_s} \right) \]

我们选择合适的 $ L $ 和 $ M $ ，使得 $ f_{max}L << F_s $ ，其中 $ f_{max} $ 是我们感兴趣的最大多普勒频移。在这种条件下，项 $ ( j 2 l f / F_s) $ ，因此可以忽略。接下来，我们定义函数 $ g(k, ) $ ：

\[ g(k, \tau) = \sum_{l=0}^{L-1} s_1(l + kM)s_2^*(l + kM – \tau) \]

然后，互模糊函数的近似表达式为

\[ \chi(\tau, f) \approx \sum_{k=0}^{M-1} g(k, \tau) \exp \left(\frac{-j 2 \pi kM f}{F_s} \right) \]

这其实就是 $ g(k, ) $ 的离散傅里叶变换！我们喜欢离散傅里叶变换，因为它们可以用 FFT 算法高效计算。快速互模糊算法的框图如下所示。

根据计算的距离和多普勒频段的数量，这种方法比直接计算算法快大约两倍！

卡尔曼滤波

在我的第一篇被动雷达文章中，我抱怨不得不通过在几百张图像中点击斑点来手动提取大量被动雷达测量数据。幸运的是，我现在找到了一种更好的方法，正如这个视频所展示的：

这个跟踪器使用了卡尔曼滤波器，这是一种广为人知的算法，已广泛应用于雷达和其他运动跟踪领域。我不会在这里详细讲解卡尔曼滤波器的全部数学细节（维基百科对此有相当简洁的介绍），但我会讲解卡尔曼滤波算法在被动雷达情况下的动态模型。接下来，我会展示一些实验结果，并比较不同参数的卡尔曼滤波器的性能。所有相关代码都可以在这里找到。

卡尔曼滤波器基础

卡尔曼滤波器通过维护一个被观测系统的内部动态模型来处理测量噪声问题。在每个时间步，滤波器使用该模型根据系统的前一个状态预测下一个状态，并生成对该预测的不确定性。接下来，获得系统的新测量值，该测量值也有其自身的不确定性。考虑到这些不确定性，卡尔曼滤波器使用预测值和测量值的加权平均值来估计系统的真正状态。

动态模型的推导

创建系统动态模型的第一步是定义一个状态向量 $ {x}(k) $ ，它指定了时间 $ k $ 时系统的状态。对于被动雷达目标，我们将使用一个包含4个元素的状态向量，其中包括目标的距离 $ r(k) $ 和多普勒频移 $ f(k) $ 以及这些值的时间导数 $ (k) $ 和 $ (k) $ 。

\[ {x}(k) = [r(k) \; \dot{r}(k) \; f(k) \; \dot{f}(k)]^T \]

现在，我们需要定义转换规则，以指定系统如何从一个状态移动到下一个状态。一个简单的选择是下面显示的“恒定速度”模型。位置值 $ r(k) $ 和 $ f(k) $ 根据它们各自的导数进行更新，而导数保持不变。

\[ r(k+1) = r(k) + \dot{r}(k)\Delta t \]

\[ \dot{r}(k+1) = \dot{r}(k) \]

\[ f(k+1) = f(k) + \dot{f}(k)\Delta t \]

\[ \dot{f}(k+1) = \dot{f}(k) \]

通过定义一个状态转移矩阵 $ {F} $ ，我们可以将状态更新方程重写为矩阵乘法的形式。

\[ {x}(k+1) = {F}{x}(k), \quad \text{其中} \quad {F} = \begin{bmatrix} 1 & \Delta t & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & \Delta t \\ 0 & 0 & 0 & 1 \end{bmatrix} \]

给定某个初始状态 $ {x}(0) $ ，可以通过应用 $ {F} $ $ k $ 次来获得状态 $ {x}(k) $ 。不难发现，这个特定模型表示一个以恒定速度沿直线在距离 - 多普勒曲面上移动的目标。

虽然恒定速度模型可能对被动雷达跟踪有效，但通过更详细地考虑情况的几何形状，可以推导出更好的模型。需要注意的是，恒定速度模型将目标距离和多普勒频移的动态视为完全独立的：值或它们的导数之间没有相互作用。然而，正如我们将看到的，实际被动雷达目标的距离和多普勒频移值是相关的。

被动雷达几何

上图显示了被动雷达中目标、发射机和接收机之间的几何关系。双基地距离定义为直接信号路径和回波信号路径长度的差。

\[ r = R_T + R_R – L \]

目标回波的多普勒频移是由于回波信号路径长度的变化引起的。下面显示了一个多普勒频移的方程，其中 $ $ 是载波信号的波长。

\[ f = – \frac{1}{\lambda} \frac{d}{dt} [R_T + R_R] \]

由于发射机 - 接收机距离 $ L $ 是恒定的，其导数为零。因此，多普勒频移可以直接用双基地距离的时间导数来表示。

\[ f = – \frac{1}{\lambda} \dot{r} \rightarrow \dot{r} = – \lambda f \]

作为附注，下面显示了一个目标回波多普勒频移关于目标几何形状的方程。这个方程对于预测具有给定位置和速度的目标所产生的多普勒频移是有用的。然而，我们实际上并不需要它来推导被动雷达状态更新模型。

\[ \frac{d}{dt}[R_T + R_R] = -2v \cos(\delta)\cos(\beta / 2) \rightarrow f = \frac{2v}{\lambda}\cos(\delta)\cos(\beta / 2) \]

在我们的新状态更新模型中，我们将保持与恒定速度模型相同的 $ f $ 和 $ $ 更新方程，即模型仍然假设恒定的多普勒速度。

\[ f(k+1) = f(k) + \dot{f}(k)\Delta t \]

\[ \dot{f}(k+1) = \dot{f}(k) \]

然后我们将距离和多普勒频移之间的关系代入距离方程：

\[ r(k+1) = r(k) + \dot{r}(k)\Delta t = r(k) – \lambda f(k) \Delta t \]

\[ \dot{r}(k+1) = – \lambda f(k+1) = – \lambda [f(k) + \dot{f}(k) \Delta t] \]

新的状态更新模型是通过将这些方程转换为矩阵形式获得的，如下所示。

\[ {x}(k+1) = {F}{x}(k), \quad \text{其中} \quad {F} = \begin{bmatrix} 1 & 0 & -\lambda \Delta t & 0 \\ 0 & 0 & -\lambda & -\lambda \Delta t \\ 0 & 0 & 1 & \Delta t \\ 0 & 0 & 0 & 1 \end{bmatrix} \]

下图比较了两种不同的状态更新模型从相同初始状态开始的行为。红色轨迹是通过依次应用恒定速度状态更新矩阵生成的，蓝色轨迹是使用恒定多普勒速度状态更新矩阵生成的。

目标跟踪逻辑

虽然有时可以让卡尔曼滤波器直接在原始输入数据上运行，但通常最好应用某种初步的数据验证。这可以防止滤波器被远离真实目标位置的虚假测量所干扰。一种解决方案是在执行卡尔曼滤波器更新步骤之前应用一个验证门，排除任何超出前一个状态估计值一定半径范围的测量值。

然而，我们只想在确信已经找到目标的情况下应用验证门。否则，我们应该考虑整个距离 - 多普勒图中的测量值，因为目标可能出现在任何位置。

选择验证门的一种简单方法是引入一个目标跟踪逻辑，它遵循类似于下图所示的状态机。在每个时间步，根据新测量值是否可能对应于当前目标来更新跟踪器状态。为了做出这个决定，我们检查新测量值是否落在最后一个状态估计值的一定半径范围内。如果是，我们沿着蓝色箭头获得下一个跟踪器状态；否则，我们沿着红色箭头前进。然后，我们根据跟踪器的状态来决定验证门的宽度。在完全锁定状态下，我们可以应用一个相当严格的验证门；而在未锁定状态下，我们根本不应用验证门。