摘要——本文针对空间有色噪声下双站多输入多输出 (MIMO) 雷达角度估计应用中的范德蒙约束 CANDECOMP/PARAFAC (CP) 张量分解问题进行了研究。

利用有色噪声的时间不相关特性，提出了一种基于时间平滑互相关方法的新型去噪方案。随后，在将平滑互相关矩阵重排为四阶张量后，构建了一个范德蒙约束 CP 张量分解模型，该模型充分利用了阵列测量的多维结构以及因子矩阵的范德蒙结构。

为了求解该模型，开发了一种高效的约束交替最小二乘 (ALS) 算法以分解范德蒙因子矩阵。最后，利用范德蒙因子矩阵的生成元获得了出发角 (DOD) 和到达角 (DOA) 的联合估计。

与现有的先进方法相比，所提方法通过联合使用时间平滑互相关去噪操作并在张量分解中强制利用范德蒙结构信息，能够实现更优的角度估计性能。数值仿真结果验证了本文算法的有效性和性能提升。

• 时间平滑互相关去噪方案。 由于脉冲数有限，现有的时间互相关算法无法完全消除有色噪声的影响。为了进一步抑制空间有色噪声，我们将匹配滤波器的输出矩阵沿时间维度划分为若干个重叠的子矩阵，计算它们的互相关矩阵，然后对这些互相关矩阵求平均，形成一个时间平滑互相关矩阵。在脉冲数有限的情况下，时间互相关矩阵的平滑估计提供了对有色噪声更好的鲁棒性。

• 面向 MIMO 雷达的范德蒙约束 CP 张量分解模型。 为了充分利用接收信号的多维结构，我们通过重排无噪声互协方差矩阵构建了一个四阶张量，然后建立了一个考虑因子矩阵中潜在范德蒙结构的约束 CP 张量分解模型。

• 基于新型范德蒙因子校正策略的约束 ALS。 我们提出了一种基于约束 ALS 的高效算法来求解范德蒙约束张量分解模型。其基本思想是将约束 ALS 问题分解为两个连续的 LS 子问题，其中第一个子问题求解因子矩阵中没有范德蒙结构约束的原始 LS 问题，而第二个子问题在范德蒙约束集中搜索第一个子问题因子矩阵的解。
  
  在第二个子问题中，我们建立了范德蒙向量的相位与其生成元之间的 LS 拟合关系，从而为范德蒙因子校正产生更好的抗噪声干扰能力。这种新的范德蒙因子校正策略属于 RecALS 族，可以被视为 [36] 中提出的方案的一种替代方案。

• 所提算法的优势。 借助时间平滑互相关和约束 ALS，所提方法表现出对空间有色噪声的鲁棒性，并实现了范德蒙约束因子矩阵的非盲估计（non-blind estimation of Vandermonde constrained factor matrices）。仿真结果表明，在存在空间有色噪声的情况下，所提方法比现有方法能产生更好的角度估计性能。

符号说明：标量、向量、矩阵和张量分别可以用小写字母、加粗小写字母、加粗大写字母和花体大写字母表示，例如 $x$、$\mathbf{x}$、$\mathbf{X}$ 和 $\mathcal{X}$。转置、复共轭、复共轭转置、逆、伪逆、Frobenius 范数和 2-范数分别表示为 $(\cdot {)}^T$、$(\cdot {)}^{\ast }$、$(\cdot {)}^H$、$(\cdot {)}^{-1}$、$(\cdot {)}^{†}$、$\Vert \cdot {\Vert }_F$ 和 $\Vert \cdot {\Vert }_2$。向量 $\mathbf{x}$ 的第 $j$ 个元素表示为 $\mathbf{x}(j)$，矩阵 $\mathbf{X}$ 的第 $(i,j)$ 个元素表示为 $\mathbf{X}(i,j)$。我们分别用 $\mathbf{X}(i,:)$ 和 $\mathbf{X}(:,j)$ 来表示矩阵 $\mathbf{X}$ 的第 $i$ 行和第 $j$ 列。三阶张量 $\mathcal{X}\in {\mathbb{C}}^{I_1\times I_2\times I_3}$ 的水平切片、侧面切片和正面切片分别表示为 $\mathcal{X}(i_1,:,:)\in {\mathbb{C}}^{I_2\times I_3}$、$\mathcal{X}(:,i_2,:)\in {\mathbb{C}}^{I_1\times I_3}$ 和 $\mathcal{X}(:,:,i_3)\in {\mathbb{C}}^{I_1\times I_2}$。算子 E { ⋅ } E\{\cdot\} E{⋅} 代表期望，符号 ${\mathbf{I}}_N$ 代表 $N\times N$ 单位矩阵。符号 $\circ$、$\otimes$ 和 $\odot$ 分别代表向量外积、Kronecker 积和 Khatri-Rao 积。最后，符号 $\text{diag}(\cdot )$、$\text{vec}(\cdot )$、$\text{phase}(\cdot )$ 和 $\arcsin (\cdot )$ 分别表示对角算子、向量化操作、向量的相位和反正弦函数。

2. Prerequisites and tensor signal model

2.1. Tensor prerequisites

以下张量基础知识对于更好地理解我们的工作至关重要。综述文章 [38,39] 提供了关于张量及其分解的更详细信息。

定义 1. (展开或矩阵化)： 张量通过展开或矩阵化操作展开为矩阵。张量 $\mathcal{X}\in {\mathbb{C}}^{I_1\times I_2\times \cdots \times I_N}$ 的模-$n$ 矩阵化可以表示为 $[\mathcal{X}{]}_n$，并且 $\mathcal{X}$ 的 $(i_1,i_2,\dots ,i_N)$ 元素对应于 $[\mathcal{X}{]}_n$ 的 $(i_n,j)$ 元素，其中 $j=1+{\sum }_{k=1,k\ne n}^N(i_k-1)J_k$ 且 $J_k={\prod }_{m=1,m\ne n}^{k-1}{I}_m$。

定义 2. (模-$n$ 张量-矩阵乘积)： 张量 $\mathcal{X}\in {\mathbb{C}}^{I_1\times I_2\times \cdots \times I_N}$ 与矩阵 $\mathbf{U}\in {\mathbb{C}}^{J\times I_n}$ 的模-$n$ 乘积可以写为 $\mathcal{Y}=\mathcal{X}{\times }_n\mathbf{U}$，其中 $\mathcal{Y}$ 是大小为 $I_1\times \cdots \times I_{n-1}\times J\times I_{n+1}\times \cdots \times I_N$ 的张量。在元素层面上，我们有

$$\mathcal{Y}(i_1,\cdots ,i_{n-1},j,i_{n+1},\cdots ,i_N)=\sum _{i_n=1}^{I_n}\mathcal{X}(i_1,i_2,\cdots ,i_N)\mathbf{U}(j,i_n)\text{\hspace{1em}(1)}$$

张量 $\mathcal{Y}$ 也可以用模-$n$ 展开格式表示为

$$\mathcal{Y}=\mathcal{X}{\times }_n\mathbf{U}\iff [\mathcal{Y}{]}_n=\mathbf{U}[\mathcal{X}{]}_n\text{\hspace{1em}(2)}$$

定义 3. (CANDECOMP/PARAFAC (CP) 分解)： CP 分解的一个目标是生成张量的秩-1 分量之和。对于 $N$ 阶张量 $\mathcal{X}\in {\mathbb{C}}^{I_1\times I_2\times \cdots \times I_N}$，$\mathcal{X}$ 的秩-$P$ CP 分解由下式给出
$$\mathcal{X}=\sum _{p=1}^P{g}_p{\mathbf{a}}_p^{(1)}\circ {\mathbf{a}}_p^{(2)}\cdots \circ {\mathbf{a}}_p^{(N)}=\mathcal{G}{\times }_1{\mathbf{A}}^{(1)}{\times }_2{\mathbf{A}}^{(2)}{\times }_3\cdots {\times }_N{\mathbf{A}}^{(N)}\text{\hspace{1em}(3)}$$
其中 $\mathcal{G}\in {\mathbb{C}}^{P\times P\times \cdots \times P}$ 是一个核心张量，其第 $(p,p,\dots ,p)$ 个元素为 $g_p$，其他位置为零。${\mathbf{A}}^{(n)}=[{\mathbf{a}}_1^{(n)},{\mathbf{a}}_2^{(n)},\dots ,{\mathbf{a}}_P^{(n)}]\in {\mathbb{C}}^{I_n\times P}(n=1,2,\cdots ,N)$ 是第 $n$ 个因子矩阵，且 ${\mathbf{a}}_p^{(n)}\in {\mathbb{C}}^{I_n\times 1}(p=1,2,\dots ,P)$ 是 ${\mathbf{A}}^{(n)}$ 的第 $p$ 个列向量。张量 $\mathcal{X}$ 的模-$n$ 矩阵化可以写为
$$[\mathcal{X}{]}_n\approx {\mathbf{A}}^{(n)}\mathbf{G}({\mathbf{A}}^{(N)}\odot \cdots \odot {\mathbf{A}}^{(n+1)}\odot {\mathbf{A}}^{(n-1)}\odot \cdots \odot {\mathbf{A}}^{(1)}{)}^T\text{\hspace{1em}(4)}$$
其中 $\mathbf{G}=\text{diag}([g_1,g_2,\cdots ,g_P])\in {\mathbb{C}}^{P\times P}$ 是一个对角矩阵。

2.2. Signal model

考虑一个具有 $M$ 个发射天线和 $N$ 个接收天线的双站 MIMO 雷达。接收和发射阵列均为均匀线性阵列 (ULA)。远场存在 $P$ 个不相关目标，其 DOD 和 DOA 分别表示为 $\mathbf{\phi }=[{\phi }_1\dots {\phi }_P{]}^T$ 和 $\mathbf{\theta }=[{\theta }_1\dots {\theta }_P{]}^T$，其中 ${\phi }_p(p=1,\dots ,P)$ 和 ${\theta }_p(p=1,\dots ,P)$ 分别代表第 $p$ 个目标的 DOD 和 DOA。发射和接收天线的位置向量分别表示为 ${\mathbf{d}}_t=[d_{t1}\dots d_{tM}{]}^T$ 和 ${\mathbf{d}}_r=[d_{r1}\dots d_{rN}{]}^T$，其中 $d_{tm}$ 和 $d_{rn}$ 分别表示第 $m$ 个 $(m=1,\dots ,M)$ 发射天线和第 $n$ 个 $(n=1,\dots ,N)$ 接收天线的位置。令 $\mathbf{S}=[{\mathbf{s}}_1,{\mathbf{s}}_2,\dots ,{\mathbf{s}}_M{]}^T\in {\mathbb{C}}^{M\times Q}$ 表示 $M$ 个发射天线具有相同带宽和中心频率的正交窄带波形，其满足

$${\mathbf{s}}_{m_1}^H{\mathbf{s}}_{m_2}=\{\begin{array}{l}0,m_1\ne m_2\\ 1,m_1=m_2\end{array}\text{\hspace{1em}(5)}$$

其中 $Q$ 表示一个脉冲周期内的采样数。假设一个相干处理间隔 (CPI) 包含 $L$ 个脉冲。通过收集 $L$ 个脉冲周期的所有信号，我们可以得到接收信号张量 $\mathcal{X}\in {\mathbb{C}}^{N\times Q\times L}$ [32,40]。

那么，第 $l$ 个脉冲的接收信号可以建模为

$$\mathcal{X}(:,:,l)={\mathbf{A}}_r({\mathbf{d}}_r,\mathbf{\theta })\text{diag}({\mathbf{c}}_l){\mathbf{A}}_t^T({\mathbf{d}}_t,\mathbf{\phi })\mathbf{S}+\mathcal{W}(:,:,l)\text{\hspace{1em}(6)}$$

其中 $\mathcal{X}(:,:,l)\in {\mathbb{C}}^{N\times Q}$ 和 $\mathcal{W}(:,:,l)\in {\mathbb{C}}^{N\times Q}$ 分别表示接收信号张量 $\mathcal{X}$ 和噪声张量 $\mathcal{W}$ 的第 $l$ 个正面切片。${\mathbf{c}}_l=[{\xi }_1{e}^{j2\pi l{f}_1/f_s},{\xi }_2{e}^{j2\pi l{f}_2/f_s},\dots ,{\xi }_P{e}^{j2\pi l{f}_P/f_s}{]}^T\in {\mathbb{C}}^{P\times 1}$ 表示第 $l$ 个脉冲的反射系数向量，其中 ${\xi }_p$ 是第 $p$ 个目标的反射系数，$f_p$ 是第 $p$ 个目标的多普勒频率，$f_s$ 表示脉冲重复频率。矩阵 ${\mathbf{A}}_t({\mathbf{d}}_t,\mathbf{\phi })=[{\mathbf{a}}_t({\mathbf{d}}_t,{\phi }_1),{\mathbf{a}}_t({\mathbf{d}}_t,{\phi }_2),\dots ,{\mathbf{a}}_t({\mathbf{d}}_t,{\phi }_P)]\in {\mathbb{C}}^{M\times P}$ 和 ${\mathbf{A}}_r({\mathbf{d}}_r,\mathbf{\theta })=[{\mathbf{a}}_r({\mathbf{d}}_r,{\theta }_1),{\mathbf{a}}_r({\mathbf{d}}_r,{\theta }_2),\dots ,{\mathbf{a}}_r({\mathbf{d}}_r,{\theta }_P)]\in {\mathbb{C}}^{N\times P}$ 分别是 $M\times P$ 发射导向矩阵和 $N\times P$ 接收导向矩阵。第 $p$ 个发射导向向量和接收导向向量可以描述为 ${\mathbf{a}}_t({\mathbf{d}}_t,{\phi }_p)=[e^{-j\pi d_{t1}\sin {\phi }_p},e^{-j\pi d_{t2}\sin {\phi }_p},\dots ,e^{-j\pi d_{tM}\sin {\phi }_p}{]}^T$ 和 ${\mathbf{a}}_r({\mathbf{d}}_r,{\theta }_p)=[e^{-j\pi d_{r1}\sin {\theta }_p},e^{-j\pi d_{r2}\sin {\theta }_p},\dots ,e^{-j\pi d_{rN}\sin {\theta }_p}{]}^T$。

通过将阵列中的第一个天线设置为参考阵元（即第一个阵元），我们可以将发射和接收天线的位置向量重写如下：

$${\mathbf{d}}_t={\alpha }_1[01\dots M-1{]}^T+{\gamma }_1\text{\hspace{1em}(7)}$$

$${\mathbf{d}}_r={\alpha }_2[01\dots N-1{]}^T+{\gamma }_2\text{\hspace{1em}(8)}$$

其中 ${\alpha }_1$ 和 ${\alpha }_2$ 分别以波长为单位表示发射和接收阵列间距，${\gamma }_1$ 和 ${\gamma }_2$ 是任意偏移量。作为一个经典选择，可以设定 $({\alpha }_1,{\gamma }_1)=({\alpha }_2,{\gamma }_2)=(1,0)$，这对应于发射和接收阵列均为半波长天线间距的 ULA。在这种情况下，阵列导向矩阵 ${\mathbf{A}}_t({\mathbf{d}}_t,\mathbf{\phi })$ 和 ${\mathbf{A}}_r({\mathbf{d}}_r,\mathbf{\theta })$ 具有范德蒙结构，其生成向量 ${\mathbf{\omega }}_t$ 和 ${\mathbf{\omega }}_r$ 描述为

$${\mathbf{\omega }}_t=[e^{-j\pi u_1},\dots ,e^{-j\pi u_P}{]}^T\text{\hspace{1em}(9)}$$

$${\mathbf{\omega }}_r=[e^{-j\pi v_1},\dots ,e^{-j\pi v_P}{]}^T\text{\hspace{1em}(10)}$$

其中 $u_p=\sin {\phi }_p$ 且 $v_p=\sin {\theta }_p$。因此，${\mathbf{A}}_t({\mathbf{d}}_t,\mathbf{\phi })$ 和 ${\mathbf{A}}_r({\mathbf{d}}_r,\mathbf{\theta })$ 可以重写为 [41]

$${\mathbf{A}}_t({\mathbf{d}}_t,\mathbf{\phi })=[{\mathbf{\omega }}_t^0{\mathbf{\omega }}_t^1\dots {\mathbf{\omega }}_t^{(M-1)}{]}^T\text{\hspace{1em}(11)}$$

3. The proposed algorithm

3.2. Joint DOD and DOA estimation

利用定义 3 中的 CP 分解模型，可以将 (21) 中的无噪声平滑互相关矩阵 $\mathbf{R}$ 重排为一个四阶张量 $\mathcal{R}\in {\mathbb{C}}^{M\times N\times M\times N}$ [30]，其表达式为$$\mathcal{R}={\mathcal{R}}_C{\times }_1{\mathbf{A}}_t{\times }_2{\mathbf{A}}_r{\times }_3{\mathbf{A}}_t^{\ast }{\times }_4{\mathbf{A}}_r^{\ast }\text{\hspace{1em}(22)}$$其中 ${\mathcal{R}}_C$ 表示一个核心张量，其第 $(p,p,p,p)$ 个 $(p=1,2,\cdots ,P)$ 元素等于 ${\lambda }_p$，其余位置为零。为了估计 DOD 和 DOA，我们需要从张量 $\mathcal{R}$ 中估计因子矩阵 ${\mathbf{A}}_t,{\mathbf{A}}_r,{\mathbf{A}}_t^{\ast }$ 和 ${\mathbf{A}}_r^{\ast }$。为此，利用最小二乘准则，该问题可以表示为 [38]$$\underset{{\mathbf{A}}_t,{\mathbf{A}}_r,{\mathbf{A}}_t^{\ast },{\mathbf{A}}_r^{\ast }}{\min }\Vert \mathcal{R}-{\mathcal{R}}_C{\times }_1{\mathbf{A}}_t{\times }_2{\mathbf{A}}_r{\times }_3{\mathbf{A}}_t^{\ast }{\times }_4{\mathbf{A}}_r^{\ast }{\Vert }_F\text{\hspace{1em}(23)}$$问题 (23) 通常被称为四阶张量 $\mathcal{R}$ 的 CP 分解模型。显然，该张量分解模型 (23) 是针对具有任意因子矩阵的张量开发的，并且可以通过盲 ALS 技术求解，该技术在没有 Vandermonde 结构约束的情况下迭代地交替更新因子矩阵。由 (11) 和 (12) 可知，在采用 ULA 的双基地 MIMO 雷达中，(22) 中的因子矩阵 ${\mathbf{A}}_t,{\mathbf{A}}_r,{\mathbf{A}}_t^{\ast }$ 和 ${\mathbf{A}}_r^{\ast }$ 具有 Vandermonde 结构，其列采取指数形式。通过将 Vandermonde 结构与常规 CP 分解模型相结合，我们提出了 Vandermonde 结构因子矩阵约束下的 CP 分解模型，其可以表示为$$\underset{\begin{array}{c}{\mathbf{A}}_t\in {\mathcal{V}}_{\phi },{\mathbf{A}}_r\in {\mathcal{V}}_{\theta }\\ {\mathbf{A}}_t^{\ast }\in {\mathcal{V}}_{\phi }^{\ast },{\mathbf{A}}_r^{\ast }\in {\mathcal{V}}_{\theta }^{\ast }\end{array}}{\min }\Vert \mathcal{R}-{\mathcal{R}}_C{\times }_1{\mathbf{A}}_t{\times }_2{\mathbf{A}}_r{\times }_3{\mathbf{A}}_t^{\ast }{\times }_4{\mathbf{A}}_r^{\ast }{\Vert }_F\text{\hspace{1em}(24)}$$其中 ${\mathcal{V}}_{\phi }$ 和 ${\mathcal{V}}_{\phi }^{\ast }$ 分别是由尺寸为 $M\times P$ 的 Vandermonde 因子矩阵组成的约束集。Vandermonde 矩阵的生成器向量可以分别表示为$${\mathbf{\omega }}_{\phi }=[{\omega }_{{\phi }_1},{\omega }_{{\phi }_2},\cdots ,{\omega }_{{\phi }_P}{]}^T\text{\hspace{1em}(25)}$$和$${\mathbf{\omega }}_{\phi }^{\ast }=[{\omega }_{{\phi }_1}^{\ast },{\omega }_{{\phi }_2}^{\ast },\cdots ,{\omega }_{{\phi }_P}^{\ast }{]}^T\text{\hspace{1em}(26)}$$

其中 ${\omega }_{{\phi }_p}=e^{-j\pi \sin {\phi }_p}$ 和 ${\omega }_{{\phi }_p}^{\ast }=e^{j\pi \sin {\phi }_p}$ 分别是生成元，且 ${\phi }_p,p=1,\dots ,P$ 是 $P$ 个目标的 DOD。利用生成元 ${\omega }_{{\phi }_p}$ 和 ${\omega }_{{\phi }_p}^{\ast }$，可以恢复 Vandermonde 矩阵 ${\mathbf{A}}_t$ 和 ${\mathbf{A}}_t^{\ast }$ 的第 $p$ 列。类似地，集合 ${\mathcal{V}}_{\theta }$ 和 ${\mathcal{V}}_{\theta }^{\ast }$ 中的矩阵是具有以下生成元向量的 $N\times P$ Vandermonde 矩阵$${\mathbf{\omega }}_{\theta }=[{\omega }_{{\theta }_1},{\omega }_{{\theta }_2},\cdots ,{\omega }_{{\theta }_P}{]}^T\text{\hspace{1em}(27)}$$和$${\mathbf{\omega }}_{\theta }^{\ast }=[{\omega }_{{\theta }_1}^{\ast },{\omega }_{{\theta }_2}^{\ast },\cdots ,{\omega }_{{\theta }_P}^{\ast }{]}^T\text{\hspace{1em}(28)}$$其中 ${\omega }_{{\theta }_p}=e^{-j\pi \sin {\theta }_p}$ 和 ${\omega }_{{\theta }_p}^{\ast }=e^{j\pi \sin {\theta }_p}$ 分别是生成元，且 ${\theta }_p,p=1,\dots ,P$ 是 $P$ 个目标的 DOA。根据定义 3 中 CP 分解的性质，四阶张量 $\mathcal{R}$ 可以重排为以下模-$n$ ($n=1,2,3,4$) 展开，即

$$\{\begin{array}{l}[\mathcal{R}{]}_1={\mathbf{A}}_t{\mathbf{R}}_C({\mathbf{A}}_r^{\ast }\odot {\mathbf{A}}_t^{\ast }\odot {\mathbf{A}}_r{)}^T\\ [\mathcal{R}{]}_2={\mathbf{A}}_r{\mathbf{R}}_C({\mathbf{A}}_r^{\ast }\odot {\mathbf{A}}_t^{\ast }\odot {\mathbf{A}}_t{)}^T\\ [\mathcal{R}{]}_3={\mathbf{A}}_t^{\ast }{\mathbf{R}}_C({\mathbf{A}}_r^{\ast }\odot {\mathbf{A}}_r\odot {\mathbf{A}}_t{)}^T\\ [\mathcal{R}{]}_4={\mathbf{A}}_r^{\ast }{\mathbf{R}}_C({\mathbf{A}}_t^{\ast }\odot {\mathbf{A}}_r\odot {\mathbf{A}}_t{)}^T\end{array}\text{\hspace{1em}(29)}$$

可以通过应用 ALS 算法对问题 (24) 进行优化，该算法在每次迭代中通过最小二乘 (LS) 方法更新 ${\mathbf{A}}_t,{\mathbf{A}}_r,{\mathbf{A}}_t^{\ast }$ 和 ${\mathbf{A}}_r^{\ast }$ 中的一个分量，同时固定其他分量，然后重复此过程直到达到收敛。根据 ALS 算法的思想，Vandermonde 约束问题 (24) 可以简化为四个线性最小二乘问题以及约束条件，得到

$$\{\begin{array}{l}\underset{{\mathbf{A}}_t\in {\mathcal{V}}_{\phi },{\mathbf{R}}_C}{\min }\Vert [\mathcal{R}{]}_1-{\mathbf{A}}_t{\mathbf{R}}_C(({\mathbf{A}}_r^{\ast }{)}^k\odot ({\mathbf{A}}_t^{\ast }{)}^k\odot ({\mathbf{A}}_r{)}^k{)}^T{\Vert }_F\\ \underset{{\mathbf{A}}_r\in {\mathcal{V}}_{\theta },{\mathbf{R}}_C}{\min }\Vert [\mathcal{R}{]}_2-{\mathbf{A}}_r{\mathbf{R}}_C(({\mathbf{A}}_r^{\ast }{)}^k\odot ({\mathbf{A}}_t^{\ast }{)}^{k+1}\odot ({\mathbf{A}}_t{)}^{k+1}{)}^T{\Vert }_F\\ \underset{{\mathbf{A}}_t^{\ast }\in {\mathcal{V}}_{\phi }^{\ast },{\mathbf{R}}_C}{\min }\Vert [\mathcal{R}{]}_3-{\mathbf{A}}_t^{\ast }{\mathbf{R}}_C(({\mathbf{A}}_r^{\ast }{)}^k\odot ({\mathbf{A}}_r{)}^{k+1}\odot ({\mathbf{A}}_t{)}^{k+1}{)}^T{\Vert }_F\\ \underset{{\mathbf{A}}_r^{\ast }\in {\mathcal{V}}_{\theta }^{\ast },{\mathbf{R}}_C}{\min }\Vert [\mathcal{R}{]}_4-{\mathbf{A}}_r^{\ast }{\mathbf{R}}_C(({\mathbf{A}}_t^{\ast }{)}^{k+1}\odot ({\mathbf{A}}_r{)}^{k+1}\odot ({\mathbf{A}}_t{)}^{k+1}{)}^T{\Vert }_F\end{array}\text{\hspace{1em}(30)}$$

由于因子矩阵 ${\mathbf{A}}_t$ 和 ${\mathbf{A}}_t^{\ast }$ 是共轭倒数对（conjugate reciprocal pairs），我们在每次迭代中只需更新 ${\mathbf{A}}_t$ 和 ${\mathbf{A}}_t^{\ast }$ 中的一个。同样，因子矩阵 ${\mathbf{A}}_r$ 和 ${\mathbf{A}}_r^{\ast }$ 的表现也类似。因此，(30) 可以简化为$$\underset{{\mathbf{A}}_t\in {\mathcal{V}}_{\phi },{\mathbf{R}}_C}{\min }\Vert [\mathcal{R}{]}_1-{\mathbf{A}}_t{\mathbf{R}}_C(({\mathbf{A}}_r^k{)}^{\ast }\odot ({\mathbf{A}}_t^k{)}^{\ast }\odot {\mathbf{A}}_r^k{)}^T{\Vert }_F\text{\hspace{1em}(31)}$$$$\underset{{\mathbf{A}}_r\in {\mathcal{V}}_{\theta },{\mathbf{R}}_C}{\min }\Vert [\mathcal{R}{]}_2-{\mathbf{A}}_r{\mathbf{R}}_C(({\mathbf{A}}_r^k{)}^{\ast }\odot ({\mathbf{A}}_t^{k+1}{)}^{\ast }\odot {\mathbf{A}}_t^{k+1}{)}^T{\Vert }_F\text{\hspace{1em}(32)}$$然后，提出了一种约束迭代 ALS 来依次最小化 (31) 和 (32)，如下所示。

(a) Updating ${\mathbf{A}}_t^{k+1}$：

我们在第 $k$ 次迭代中通过求解优化子问题 (31) 来更新变量 ${\mathbf{A}}_t^{k+1}$，这是一个典型的约束最小二乘 (LS) 问题。本质上，这个约束 LS 问题可以分为两个连续的 LS 子问题 [42]。对于第一个子问题，我们求解无约束的原始 LS 问题。对于第二个子问题，我们在约束集中找到近似于第一个子问题最优解的点。因此，求解问题 (31) 的方法涉及两个阶段。

• 在第一阶段，我们求解无约束的 LS 问题
  $$\underset{{\mathbf{U}}_t}{\min }\Vert [\mathcal{R}{]}_1-{\mathbf{U}}_t(({\mathbf{A}}_r^k{)}^{\ast }\odot ({\mathbf{A}}_t^k{)}^{\ast }\odot {\mathbf{A}}_r^k{)}^T{\Vert }_F\text{\hspace{1em}(33)}$$
  其中 ${\mathbf{U}}_t={\mathbf{A}}_t{\mathbf{R}}_C$。我们可以看到，在这个阶段，因子矩阵上的 Vandermonde 结构约束被移除了。矩阵 ${\mathbf{U}}_t$ 的解可以获得为
  $${\mathbf{U}}_t^{k+1}=[\mathcal{R}{]}_1{\left((({\mathbf{A}}_r^k{)}^{\ast }\odot ({\mathbf{A}}_t^k{)}^{\ast }\odot {\mathbf{A}}_r^k{)}^T\right)}^{†}\text{\hspace{1em}(34)}$$

• 在第二阶段，简单的约束 LS 问题（其中约束集 ${\mathcal{V}}_{\phi }$ 实际上包含因子矩阵的先验结构信息）可以写为
  $$\begin{array}{rl}& \underset{{\mathbf{A}}_t\in {\mathcal{V}}_{\phi },{\mathbf{R}}_C}{\min }\Vert {\mathbf{U}}_t^{k+1}-{\mathbf{A}}_t{\mathbf{R}}_C{\Vert }_F\\ & =\underset{{\mathbf{A}}_t(:,p)\in {\mathcal{V}}_{{\phi }_p},{\lambda }_p}{\min }\sum _{p=1}^P\Vert {\mathbf{U}}_t^{k+1}(:,p)-{\mathbf{A}}_t(:,p){\lambda }_p{\Vert }_F\end{array}\text{\hspace{1em}(35)}$$
  其中 ${\lambda }_p$ 是对角矩阵 ${\mathbf{R}}_C$ 的第 $p$ 个对角元素，且 ${\mathcal{V}}_{{\phi }_p}$ 是由 $M\times 1$ Vandermonde 向量组成的约束集，其生成元为 ${\omega }_{{\phi }_p}=e^{-j\pi u_p}$ 和 $u_p=\sin {\phi }_p$。很容易证明问题 (35) 等价于独立的最小化问题
  $$\underset{{\mathbf{A}}_t(:,p)\in {\mathcal{V}}_{{\phi }_p},{\lambda }_p}{\min }\Vert {\mathbf{U}}_t^{k+1}(:,p)-{\mathbf{A}}_t(:,p){\lambda }_p{\Vert }_F,p=1,2,\cdots ,P\text{\hspace{1em}(36)}$$
  

正确的公式应该是：$${\mathbf{A}}_t(:,p)=[1,e^{-j\pi u_p},\cdots ,e^{-j\pi (M-1)u_p}{]}^T$$

向量 ${\mathbf{A}}_t(:,p)$ 是 Vandermonde 矩阵 ${\mathbf{A}}_t$ 的第 $p$ 列，且 ${\mathbf{A}}_t(:,p)$ 的形式给出为 ${\mathbf{A}}_t(:,p)=[1,e^{-j\pi u_1},\cdots ,e^{-j\pi (M-1)u_p}{]}^T$。显然，${\mathbf{A}}_t(:,p)$ 的估计等价于 $u_p$ 的估计，这可以通过 LS 拟合原理 [43] 来求解。构造以下矩阵和向量如下：

$${\mathbf{P}}_t={\left[\begin{array}{cccc}1& 1& \cdots & 1\\ 0& \pi & \cdots & (M-1)\pi \end{array}\right]}^T\text{\hspace{1em}(37)}$$

$$\begin{array}{rl}{\mathbf{h}}_{tp}^{k+1}& =-phase({\mathbf{U}}_t^{k+1}(:,p))\\ & =[{\eta }_t,{\eta }_t+\pi u_p,\cdots ,{\eta }_t+\pi (M-1)u_p{]}^T+{\mathbf{\epsilon }}_t\end{array}\text{\hspace{1em}(38)}$$

其中 ${\eta }_t$ 是相位偏移，${\mathbf{\epsilon }}_t$ 表示相位误差向量。为了降低对误差的敏感度，我们构建针对 $u_p$ 的 LS 拟合问题如下

$$\underset{{\mathbf{c}}_{tp}}{\min }\Vert {\mathbf{P}}_t{\mathbf{c}}_{tp}-{\mathbf{h}}_{tp}^{k+1}{\Vert }_2\text{\hspace{1em}(39)}$$

其中 ${\mathbf{c}}_{tp}\in {\mathbb{C}}^{2\times 1}$ 表示估计向量，且 ${\mathbf{c}}_{tp}$ 的第二个元素是 $u_p$ 的估计值。(39) 中 ${\mathbf{c}}_{tp}$ 的解可以由下式给出

$${\hat{\mathbf{c}}}_{tp}={\mathbf{P}}_t^{†}{\mathbf{h}}_{tp}^{k+1}\text{\hspace{1em}(40)}$$

因此，$u_p$ 的估计值可以通过 ${\hat{u}}_p={\hat{\mathbf{c}}}_{tp}(2)$ 获得，其可用于估计 Vandermonde 向量 ${\hat{\mathbf{A}}}_t^{k+1}(:,p)=[1,e^{-j\pi {\hat{u}}_p},\cdots ,e^{-j\pi (M-1){\hat{u}}_p}{]}^T$。因此，因子矩阵 ${\hat{\mathbf{A}}}_t^{k+1}$ 可以逐列进行估计。根据 (36)，${\lambda }_p$ 的解可以获得为

$${\hat{\lambda }}_p^{k+1}=e^{-j{\hat{c}}_{tp}(1)}\cdot \Vert {\mathbf{U}}_t^{k+1}(:,p){\Vert }_2/\Vert {\mathbf{A}}_t^{k+1}(:,p){\Vert }_2\text{\hspace{1em}(41)}$$

估计矩阵 ${\mathbf{U}}_t^{k+1},{\mathbf{U}}_r^{k+1}$ 与 ${\mathbf{A}}_t^{k+1},{\mathbf{A}}_r^{k+1}$ 之间的关系可以表述为

$${\mathbf{U}}_t^{k+1}={\mathbf{A}}_t^{k+1}\mathbf{\Pi }{\mathbf{\Delta }}_t+{\mathbf{N}}_t\text{\hspace{1em}(49)}$$

$${\mathbf{U}}_r^{k+1}={\mathbf{A}}_r^{k+1}\mathbf{\Pi }{\mathbf{\Delta }}_r+{\mathbf{N}}_r\text{\hspace{1em}(50)}$$

其中 $\mathbf{\Pi }\in {\mathbb{C}}^{P\times P}$ 是一个排列矩阵。${\mathbf{\Delta }}_t$ 和 ${\mathbf{\Delta }}_r$ 是 $P\times P$ 对角矩阵，其对角线元素表示缩放效应。${\mathbf{N}}_t$ 和 ${\mathbf{N}}_r$ 代表拟合误差。由于矩阵 ${\mathbf{U}}_t^{k+1}$ 和 ${\mathbf{U}}_r^{k+1}$ 具有相同的列模糊性，且 ${\hat{u}}_p$ 和 ${\hat{v}}_p,p=1,2,\dots P$ 是分别从 ${\mathbf{U}}_t^{k+1}$ 和 ${\mathbf{U}}_r^{k+1}$ 中逐列估计得到的，因此 ${\hat{u}}_p$ 和 ${\hat{v}}_p$ 的独立估计值会自动针对同一目标进行配对。回想一下，${\hat{\mathbf{A}}}_t^{k+1}$ 和 ${\hat{\mathbf{A}}}_r^{k+1}$ 是 Vandermonde 矩阵，并分别由其生成元 ${\omega }_{{\phi }_p}=e^{-j\pi {\hat{u}}_p}$ 和 ${\omega }_{{\theta }_p}=e^{-j\pi {\hat{v}}_p},p=1,2,\dots P$ 重构，这使得 ${\hat{\mathbf{A}}}_t^{k+1}$ 和 ${\hat{\mathbf{A}}}_r^{k+1}$ 也具有相同的列模糊性。

当达到最大迭代次数或两次连续迭代之间的张量分解损失差小于预定义的小值时，迭代将停止，即

$$(\Vert \mathcal{R}-{\hat{\mathcal{R}}}^{k+1}{\Vert }_F-\Vert \mathcal{R}-{\hat{\mathcal{R}}}^k{\Vert }_F)/\Vert \mathcal{R}{\Vert }_F\le 10^{-8},\text{其中}{\hat{\mathcal{R}}}^k={\hat{\mathcal{R}}}_C^k{\times }_1{\mathbf{A}}_t^k{\times }_2{\mathbf{A}}_r^k{\times }_3({\mathbf{A}}_t^k{)}^{\ast }{\times }_4({\mathbf{A}}_r^k{)}^{\ast }$$

并且 ${\hat{\mathcal{R}}}_C^k$ 表示核心张量，其第 $(p,p,p,p)$ 个 $(p=1,2,\cdots ,P)$ 元素等于 ${\lambda }_p^k$，其余位置为零。在本文中，采用了 ESPRIT 方法 [44] 算法来初始化因子矩阵以加速收敛。在从最后一次迭代获得 ${\hat{u}}_p$ 和 ${\hat{v}}_p$ 的估计值之后，DOD 和 DOA 导出为

$${\hat{\phi }}_p=\arcsin ({\hat{u}}_p)\text{\hspace{1em}(51)}$$$${\hat{\theta }}_p=\arcsin ({\hat{v}}_p)\text{\hspace{1em}(52)}$$

我们将所提方法的实现步骤总结在算法 1 中。

5. Simulation results

在本节中，给出了在空间色噪声下评估所提角度估计算法性能的仿真结果，并与 Liao 的方法 [27]、Wen 的方法 [29]、Cai 的方法 [30]、VTT-ReCALS [37]、ESPRIT 算法 [44] 和 ND-ESPRIT [46] 的性能进行了比较。

由于 VTT-ReCALS 和 ND-ESPRIT 方法都不适用于色噪声情况，因此在利用 VTT-ReCALS 或 ND-ESPRIT 方法进行角度估计之前，采用了所提的时间平滑互相关去噪方法来消除色噪声的影响。

还提供了双基地 MIMO 雷达中 DOD 和 DOA 估计的克拉美-罗界 (CRB) [43]。我们考虑配置有 10 个发射阵元和 12 个接收阵元的双基地 MIMO 雷达。发射和接收阵列均使用半波长间距的 ULA。第 $m$ 个发射波形 $s_m$ 由 $s_m=(1+j)/\sqrt{2}{\mathbf{h}}_m$ 产生，其中 ${\mathbf{h}}_m$ 表示 $Q\times Q$ Hadamard 矩阵的第 $m$ 行，且 $Q=256$。考虑 $P=3$ 个不相关信源，位于角度 $({\phi }_1,{\theta }_1)=(-10^{\circ },-15^{\circ }),({\phi }_2,{\theta }_2)=(8^{\circ },-5^{\circ })$ 和 $({\phi }_3,{\theta }_3)=(0^{\circ },12^{\circ })$，多普勒频率为 { f p } p = 1 3 = ( 200 , 400 , 800 ) Hz \{f_p\}_{p=1}^3 = (200, 400, 800) \text{Hz} {fp​}p=13​=(200,400,800)Hz。空间色噪声协方差矩阵 $\mathbf{D}$ 中的第 $(p,q)$ 个元素遵循生成表达式 [30] $\mathbf{D}(p,q)=0.9^{|p-q|}{e}^{j\pi (p-q)/2}$。信噪比 (SNR) 定义为 $SNR=10\lg (\Vert \mathcal{X}-\mathcal{W}{\Vert }_F^2/\Vert \mathcal{W}{\Vert }_F^2)\text{dB}$。使用两个指标评估性能，即均方根误差 (RMSE) 和成功检测概率 (probability of the successful detection，PSD)。

• RMSE 定义为
  $$RMSE=\sqrt{\frac{1}{2M_{T}P}\sum _{p=1}^P\sum _{i=1}^{M_T}\left[({\hat{\phi }}_{p,i}-{\phi }_p{)}^2+({\hat{\theta }}_{p,i}-{\theta }_p{)}^2\right]}\text{\hspace{1em}(55)}$$
  其中 $({\hat{\phi }}_{p,i},{\hat{\theta }}_{p,i})$ 是第 $i$ 次蒙特卡洛试验中第 $p$ 个目标的估计 DOD 和 DOA 对，$M_T$ 表示蒙特卡洛试验的次数。对于每次仿真，执行 200 次蒙特卡洛试验。

• PSD 的定义是 $PSD=T/200\times 100\%$，其中 $T$ 是所有估计角度的绝对误差小于 $0.2^{\circ }$ 的成功试验次数。

实验 1： 在本实验中，实验研究了重叠参数 $\tilde{L}$ 对所提算法的影响。图 1 描绘了在不同的 $\tilde{L}$ 值下 RMSE 性能随 SNR 变化的情况，其中收集了 $L=100$ 个脉冲，且 $\tilde{L}$ 从 1 变化到 5。当 $\tilde{L}=1$ 时，用于消除色噪声的所提时间平滑互相关方案退化为 [29-31] 中提出的时间互相关方案。从图 1 中观察到一个有趣的现象，在低 SNR 区域，将 $\tilde{L}$ 增加到 3 显著提高了角度估计性能，而将 $\tilde{L}$ 增加到超过 3 可能不会改善甚至会恶化角度估计性能。此外，在高 SNR 区域，增加 $\tilde{L}$ 可能会略微恶化 RMSE。在没有关于 SNR 的先验信息可用的条件下，在以下对所提方法的仿真中选择了折衷值 $\tilde{L}=3$，这在整个 SNR 范围内实现了可接受的 RMSE。基于在不同脉冲数下的大量实证仿真，建议选择重叠参数值为 $\tilde{L}=(0.02\sim 0.03)\times L$，这将提供良好的性能。