一种利用ESPRIT的线阵DOA与极化参数联合估计

A Joint Estimation of DOA and Polarization Parameters for Linear Arrays Based on ESPRIT

摘要

针对线性极化敏感阵列多参数联合估计中计算复杂度高的问题，提出一种基于旋转不变子空间（ESPRIT）算法的二维波达方向（DOA）和极化参数联合估计算法。首先，构建由极化敏感阵列组成的均匀线阵（ULA）模型，通过对接收数据的协方差矩阵进行特征分解，将信号子空间划分成5类子阵；然后，基于旋转不变特性，分别利用其中两组子阵的相位差闭式求解俯仰角，另3组子阵比值的相位差直接推导方位角，实现方向参数的独立解算；最后，完成极化参数的闭式估计。算法避免了传统多维谱峰搜索与迭代操作，运算量显著降低，保证了实时性。实验结果表明，该算法除极化幅度角外，其他参数估计精度均优于降维MUSIC与Ratio-ESPRIT算法，且计算效率显著优于降维MUSIC算法。

Abstract

To address the issue of high computational complexity in joint multi-parameter estimation for linear polarization-sensitive arrays, an ESPRIT-based algorithm is proposed for the joint estimation of two-dimensional direction of arrival (DOA) and polarization parameters. First, a uniform linear array (ULA) model composed of polarization-sensitive arrays is constructed. By performing eigenvalue decomposition on the covariance matrix of the received data, the signal subspace is divided into five sub-arrays. Based on the rotational invariance property, elevation angles are then obtained in closed form using the phase differences from two of the sub-arrays, while azimuth angles are derived directly from the phase differences of the ratios among the other three sub-arrays. Independent estimation of directional parameters is thereby achvieved. Finally, closed-form estimation of the polarization parameters is achieved. Traditional multi-dimensional spectral-peak searches and iterative operations are avoided by using the algorithm, which significantly reduces the computational load and ensures real-time performance. Experimental results show that the estimation accuracies of all parameters are higher by using the proposed algorithm than those of the reduced-dimension MUSIC algorithm and the Ratio-ESPRIT algorithm except that of the polarization amplitude angle, while its computational efficiency is significantly superior to that of the reduced-dimension MUSIC algorithm.

Graphical abstract

关键词

DOA估计 / 极化参数估计 / 极化敏感阵列 / 均匀线阵 / 旋转不变子空间

Key words

DOA estimation / polarization parameters estimation / polarization-sensitive array / uniform linear array (ULA) / estimation of signal parameters via rotational invariance techniques (ESPRIT)

波达方向（Direction of Arrival, DOA）估计是阵列信号处理领域的核心问题，在雷达目标跟踪、无线通信波束成形及电子侦察等场景中具有重要应用价值^[1-4]。随着电磁环境的日益复杂，电磁干扰技术的不断发展，标量阵列因信息维度有限往往易受多径效应、极化干扰等因素的影响，导致分辨力下降和估计误差增大^[5-6]。极化敏感阵列作为一种新型阵列结构，通过将极化信息与空域信息融合，具有较强的空域滤波能力和参数估计性能。目前，极化敏感阵列已广泛应用于通信、无线电和导航领域^[7-8]。

得益于极化敏感阵列的优异性能，一些基于传统标量阵列的DOA估计算法已经应用到极化敏感阵列，衍生出许多DOA与极化参数联合估计的算法^[9-12]。例如，文献[13]基于正交偶极子设计了极化多重信号分类（Multiple Signal Classification, MUSIC）算法；文献[14]基于极化敏感阵列提出了牛顿迭代法的四维联合参数快速估计算法。除了方阵，基于稀疏拉伸式的L型阵列也有良好的参数估计性能^[15]，提高了阵列孔径。这些算法均是针对面阵设计的DOA估计优化算法。但在现代电子战中，此类面阵依赖特定孔径大小和阵元数量来保证角度估计精度，成本和复杂度较高，需要较大资源代价才能获取二维DOA信息，而全极化敏感矢量传感器的一维排布方式将会有更高的灵活性，可显著简化系统。

极化敏感线阵的一维拓扑性导致其获取空间参数维度有限，传统的正交偶极子阵列无法实现二维DOA估计。例如，Wang等^[16]利用旋转不变子空间（Estimation of Signal Parameters via Rotational Invariance Techniques, ESPRIT）算法设计了均匀线阵的DOA和频率联合估计算法；Zhang等^[17]改进了MUSIC算法，通过降维减小复杂度，实现一维DOA与极化参数的联合估计；为避免谱峰搜索，文献[18]提出的模值约束Root-MUSIC算法进一步降低了计算复杂度；为了提高阵列虚拟孔径数量，文献[19]提出了互质极化敏感阵列的降维快速估计算法。Lan等^[20]将降维MUSIC（Reduced-Dimension MUSIC, RD-MUSIC）算法应用到一维极化敏感线阵中，虽改进了传统的四维搜索方法，但其二维搜索带来的复杂度依然很高。尽管这些算法采用的均为极化敏感线阵，但无法同时满足二维DOA信息获取与低复杂度的要求。

针对以上问题，本文提出基于线性极化敏感阵列的二维DOA快速估计算法，旨在通过不同方向的偶极子天线的相位关系实现多参数的解耦运算。近年来，极化敏感阵列的研究热点主要是面阵、稀疏阵等结构^[15,19]，这些算法虽提升了精度和自由度，但需要更大的布阵空间和阵元资源，与线阵的结构特性存在差异。相比之下，针对线性极化敏感阵列的二维DOA估计研究相对有限，其中线阵降维MUSIC算法^[20]与线阵Ratio-ESPRIT算法^[21]是解决该问题的代表性工作，因此本文将同时与这两种算法进行对比。不同于Ratio-ESPRIT算法，本文算法改进阵列流形信号的子空间分解策略，将其划分为5个子阵，俯仰角和方位角依赖线阵的多重旋转不变性，分别从相位关联子阵和轴向偶极子响应矩阵中解算，最终实现极化参数的闭式解估计。相较于基于极化敏感线阵的降维MUSIC算法，本文算法规避了谱峰搜索与迭代优化，降低了算法复杂度。此外，相较于基于交叉偶极子的线阵Ratio-ESPRIT算法^[21]，本文算法在DOA估计精度上显著优于该算法。因此，本文算法在获取二维DOA信息的同时，实现了相关参数的高精度估计。

1 阵列结构与信号模型

接收天线阵列模型如图1所示。该阵列共由L个极化敏感阵元组成，每个极化敏感阵元由3个正交偶极子组成，总共有\mathrm{3}L个元素。每个极化敏感阵元的中心均匀分布在Z轴上，相邻两个阵元之间的间隔为\delta，工作波长为\lambda，为了避免角度模糊问题，\delta满足\delta \le \lambda /\mathrm{2}。

假设共有K个连续波信号入射到阵列上，第k个信号的角度信息为\left({\theta }_k,{\varphi }_k\right)，{\theta }_k是信号的俯仰角，{\varphi }_k是信号的方位角，其中k满足k\le L-\mathrm{1}。给定一个入射到阵列上的横向电磁波，向量{\mathit{e}}_{\mathit{\gamma }}代表信号的入射方向，为了表示信号的极化方式，定义单位向量{\mathit{e}}_{\varphi }、{\mathit{e}}_{\theta }、{\mathit{e}}_{\mathit{\gamma }}按顺序形成的右手坐标系，根据Jones矢量法的描述，该电场矢量\mathit{E}可以表示为

\mathit{E}=E_{\varphi }{\mathit{e}}_{\varphi }+E_{\theta }{\mathit{e}}_{\theta }=\left[\begin{array}{c}E\mathrm{c}\mathrm{o}\mathrm{s}\gamma \\ E\mathrm{s}\mathrm{i}\mathrm{n}\gamma \cdot e^{\mathrm{j}\eta }\end{array}\right]（1）

式中：E_{\varphi }和E_{\theta }分别表示电磁波在{\mathit{e}}_{\varphi }和{\mathit{e}}_{\theta }方向上的幅度分量；\gamma和\eta分别表示电磁波的极化幅度角和极化相位角。

因此，1个入射信号可以由1个振幅和4个角度参数来表征。例如，第k个信号的参数信息为\left({\theta }_k,{\varphi }_k,{\gamma }_k,{\eta }_k,E_k\right)，\left({\theta }_k,{\varphi }_k\right)表示信号的DOA参数；\left({\gamma }_k,{\eta }_k\right)表示信号的极化参数；E_k表示信号的幅度；k=1,2,…,K。

根据几何关系，对于给定的阵元，经过坐标转换，电磁波接收矢量可以表示为

\mathit{E}={\left[\begin{array}{c}{\mathit{e}}_x\\ {\mathit{e}}_y\\ {\mathit{e}}_z\end{array}\right]}^{\mathrm{{\rm T}}}\cdot \left[\begin{array}{cc}-\mathrm{s}\mathrm{i}\mathrm{n}\varphi & \mathrm{c}\mathrm{o}\mathrm{s}\theta \mathrm{c}\mathrm{o}\mathrm{s}\varphi \\ \mathrm{c}\mathrm{o}\mathrm{s}\varphi & \mathrm{c}\mathrm{o}\mathrm{s}\theta \mathrm{s}\mathrm{i}\mathrm{n}\varphi \\ \mathrm{0}& -\mathrm{s}\mathrm{i}\mathrm{n}\theta \end{array}\right]\cdot \left[\begin{array}{c}E\mathrm{c}\mathrm{o}\mathrm{s}\gamma \\ E\mathrm{s}\mathrm{i}\mathrm{n}\gamma \cdot e^{\mathrm{j}\eta }\end{array}\right]=E\left(-\mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\varphi +\mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{c}\mathrm{o}\mathrm{s}\theta \mathrm{c}\mathrm{o}\mathrm{s}\varphi \cdot e^{\mathrm{j}\eta }\right){\mathit{e}}_x+E\left(\mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{c}\mathrm{o}\mathrm{s}\varphi +\mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{c}\mathrm{o}\mathrm{s}\theta \mathrm{s}\mathrm{i}\mathrm{n}\gamma \cdot e^{\mathrm{j}\eta }\right){\mathit{e}}_y+E\left(-\mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\theta \cdot e^{\mathrm{j}\eta }\right){\mathit{e}}_z（2）

根据阵元位置不同，定义空间相位因子：

q=e^{\mathrm{j}\frac{\mathrm{2}\mathrm{\pi }}{\lambda }\delta \mathrm{c}\mathrm{o}\mathrm{s}\theta }（3）

将时间和空间相位因子代入公式，假设每个信号向量中存在噪声电压向量{\mathit{n}}_i\left(t\right)，且{\mathit{n}}_i\left(t\right)是服从均值为零的高斯分布的随机变量,i=\mathrm{1,2},\cdots ,L。因此第i个阵元接收的信号向量为

{\mathit{z}}_i\left(t\right)={\displaystyle \sum _{k=\mathrm{1}}^K}{\mathit{u}}_k{s}_k\left(t\right)q_k^{i-\mathrm{1}}+{\mathit{n}}_i\left(t\right)（4）

其中：

{\mathit{u}}_k=\left[\begin{array}{c}-\mathrm{c}\mathrm{o}\mathrm{s}{\gamma }_k\mathrm{s}\mathrm{i}\mathrm{n}{\varphi }_k+\mathrm{s}\mathrm{i}\mathrm{n}{\gamma }_k\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_k\mathrm{c}\mathrm{o}\mathrm{s}{\varphi }_k\cdot e^{\mathrm{j}{\eta }_k}\\ \mathrm{c}\mathrm{o}\mathrm{s}{\gamma }_k\mathrm{c}\mathrm{o}\mathrm{s}{\varphi }_k+\mathrm{s}\mathrm{i}\mathrm{n}{\gamma }_k\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_k\mathrm{s}\mathrm{i}\mathrm{n}{\varphi }_k\cdot e^{\mathrm{j}{\eta }_k}\\ -\mathrm{s}\mathrm{i}\mathrm{n}{\gamma }_k\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k\cdot e^{\mathrm{j}{\eta }_k}\end{array}\right]（5）

s_k\left(t\right)=E_k\cdot e^{\mathrm{j}\left(\omega t+\phi \right)}（6）

式中：\omega是信号频率；\phi是信号的初始相位；后续E_k的值均取1。

设\mathit{z}\left(t\right)、\mathit{s}\left(t\right)和\mathit{n}\left(t\right)分别是包含接收信号、入射信号和噪声的列向量，即：

\mathit{z}\left(t\right)={\left[\begin{array}{cccc}{\mathit{z}}_{\mathrm{1}}\left(t\right)& {\mathit{z}}_{\mathrm{2}}\left(t\right)& \cdots & {\mathit{z}}_L\left(t\right)\end{array}\right]}^{\mathrm{{\rm T}}}（7）

\mathit{s}\left(t\right)={\left[\begin{array}{cccc}{\mathit{s}}_{\mathrm{1}}\left(t\right)& {\mathit{s}}_{\mathrm{2}}\left(t\right)& \cdots & {\mathit{s}}_K\left(t\right)\end{array}\right]}^{\mathrm{{\rm T}}}（8）

\mathit{n}\left(t\right)={\left[\begin{array}{cccc}{\mathit{n}}_{\mathrm{1}}\left(t\right)& {\mathit{n}}_{\mathrm{2}}\left(t\right)& \cdots & {\mathit{n}}_L\left(t\right)\end{array}\right]}^{\mathrm{{\rm T}}}（9）

接收向量则表示为

\mathit{z}\left(t\right)=\mathit{A}\mathit{s}\left(t\right)+\mathit{n}(t)（10）

式中，\mathit{A}=\left[\begin{array}{cccc}{\mathit{a}}_{\mathrm{1}}& {\mathit{a}}_{\mathrm{2}}& \cdots & {\mathit{a}}_K\end{array}\right]是一个\mathrm{3}L\times K的阵列流型矩阵，每列由下式给出：

{\mathit{a}}_k={\left[\begin{array}{cccc}{\mathit{u}}_k^{\mathrm{{\rm T}}}& {\mathit{u}}_k^{\mathrm{{\rm T}}}{q}_k& \cdots & {\mathit{u}}_k^{\mathrm{{\rm T}}}{q}_k^{L-\mathrm{1}}\end{array}\right]}^{\mathrm{{\rm T}}}（11）

这里假设{\mathit{a}}_k是线性无关的，它们在一个L维空间中定义了一个K维信号子空间。

为了保证算法的可行性，要求\mathit{A}是列满秩的矩阵，因此排除了一些退化情况。假设阵列在N个不同时间对入射信号进行采样，算法的目标是从采样信号{\mathit{z}}_{\mathit{k}}\left(t_n\right)中得到估计量\left({\theta }_k,{\varphi }_k,{\gamma }_k,{\eta }_k\right)。其中，n=\mathrm{1,2},\cdots ,N,k=\mathrm{1,2},\cdots ,K。

2 二维DOA与极化参数联合估计算法

由于入射信号和噪声之间均线性无关，根据ESPRIT算法的原理可知，上述矩阵具有多个不变特性。假设入射信号数量已知，得到接收信号的协方差矩阵：

\mathit{R}=\mathrm{E}\left\{\mathit{z}\left(t\right){\mathit{z}}^{\mathrm{{\rm H}}}\left(t\right)\right\}=\mathit{A}{\mathit{R}}_{\mathit{s}}{\mathit{A}}^{\mathrm{{\rm H}}}+{\mathit{R}}_n={\mathit{R}}_{\mathrm{0}}+{\sigma }^{\mathrm{2}}\mathit{I}（12）

式中：{\mathit{R}}_{\mathit{s}}=\mathrm{E}\left\{\mathit{s}\left(t\right){\mathit{s}}^{\mathrm{H}}\left(t\right)\right\}；{\sigma }^{\mathrm{2}}是高斯白噪声方差；\mathit{I}是单位阵；\mathrm{E}(·)表示期望。

从特征根解思路求解，可以将协方差矩阵\mathit{R}写成：

\mathit{R}={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{3}L}}{\lambda }_i{\mathit{e}}_i{\mathit{e}}_i^{\mathrm{{\rm H}}}=\mathit{U}\mathit{\Lambda }{\mathit{U}}^{\mathrm{{\rm H}}}={\mathit{E}}_s{\mathit{\Lambda }}_s{\mathit{E}}_s^{\mathrm{{\rm H}}}+{\mathit{E}}_n{\mathit{\Lambda }}_n{\mathit{E}}_n^{\mathrm{{\rm H}}}（13）

在阵列的协方差矩阵中，{\lambda }_{\mathrm{1}}+{\sigma }^{\mathrm{2}}\ge {\lambda }_{\mathrm{2}}+{\sigma }^{\mathrm{2}}\ge \cdots

\ge {\lambda }_K+{\sigma }^{\mathrm{2}}>{\sigma }^{\mathrm{2}}=\cdots ={\sigma }^{\mathrm{2}}是\mathit{R}的特征值；{\mathit{e}}_{\mathrm{1}},{\mathit{e}}_{\mathrm{2}},\cdots ,{\mathit{e}}_K,

{\mathit{e}}_{K+\mathrm{1}},\cdots ,{\mathit{e}}_{\mathrm{3}L}是对应的特征向量；\mathit{\Lambda }、{\mathit{\Lambda }}_s、{\mathit{\Lambda }}_n均是对角矩阵；\mathit{U}是特征向量矩阵。根据式（13）可以证得：

\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\left({\mathit{E}}_s\right)=\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\left(\mathit{A}\right)（14）

式中：{\mathit{E}}_s=\left[\begin{array}{cccc}{\mathit{e}}_{\mathrm{1}}& {\mathit{e}}_{\mathrm{2}}& \cdots & {\mathit{e}}_K\end{array}\right]；\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}(·)表示张成子空间。因此，存在非奇异矩阵\mathit{T}，满足：

{\mathit{E}}_s=\mathit{A}\mathit{T}（15）

将阵列流型矩阵\mathit{A}进行拆分，假设{\mathit{A}}_{q\mathrm{1}}和{\mathit{A}}_{q\mathrm{2}}均是\mathrm{3}\left(L-\mathrm{1}\right)\times K的矩阵，且分别由\mathit{A}的前\mathrm{3}\left(L-\mathrm{1}\right)行和后\mathrm{3}\left(L-\mathrm{1}\right)行组成；假设{\mathit{A}}_x、{\mathit{A}}_y和{\mathit{A}}_z均是L\times K的矩阵，且分别由\mathit{A}的(\mathrm{3}i-\mathrm{2}）行，(\mathrm{3}i-\mathrm{1}）行和\mathrm{3}i行组成，即分别由平行于X轴、Y轴和Z轴的偶极子阵列流型组成。根据ESPRIT的旋转不变性，得到等式关系：

\left\{\begin{array}{l}{\mathit{E}}_{q\mathrm{1}}={\mathit{A}}_{q\mathrm{1}}\mathit{T};\\ {\mathit{E}}_{q\mathrm{2}}={\mathit{A}}_{q\mathrm{2}}\mathit{T}.\end{array}\right.（16）

其中：

\left\{\begin{array}{l}{\mathit{E}}_{q\mathrm{2}}={\mathit{E}}_{q\mathrm{1}}{\mathit{\Psi }}_{\mathrm{1}};\\ {\mathit{A}}_{q\mathrm{2}}={\mathit{A}}_{q\mathrm{1}}{\mathit{\Phi }}_{\mathrm{1}}.\end{array}\right.（17）

式中，{\mathit{\Phi }}_{\mathrm{1}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left[\begin{array}{cccc}{q}_{\mathrm{1}}& q_{\mathrm{2}}& \cdots & q_K\end{array}\right]，{\mathit{\Phi }}_{\mathrm{1}}与{\mathit{\Psi }}_{\mathrm{1}}的特征值相同。{\mathit{\Psi }}_{\mathrm{1}}的表达式为

{\mathit{\Psi }}_{\mathrm{1}}={\left({\mathit{E}}_{q\mathrm{1}}^{\mathrm{H}}{\mathit{E}}_{q\mathrm{1}}\right)}^{-\mathrm{1}}{\mathit{E}}_{q\mathrm{1}}^{\mathrm{H}}{\mathit{E}}_{q\mathrm{2}}={\mathit{E}}_{q\mathrm{1}}^{†}{\mathit{E}}_{q\mathrm{2}}（18）

式中，{(·)}^{†}表示伪逆。因此信号的俯仰角{\theta }_k表达式为

{\theta }_k=\mathrm{c}\mathrm{o}{\mathrm{s}}^{-\mathrm{1}}\left[\frac{\lambda }{\mathrm{2}\mathrm{\pi }\delta }\mathrm{a}\mathrm{r}\mathrm{g}\left(q_k\right)\right]（19）

由于俯仰角{\theta }_k\in \left(\mathrm{0},\frac{\mathrm{\pi }}{\mathrm{2}}\right)，其不存在象限模糊问题。

为了求解{\varphi }_k，定义比值p_k和r_k为：

p_k=\frac{-\mathrm{c}\mathrm{o}\mathrm{s}{\gamma }_k\mathrm{s}\mathrm{i}\mathrm{n}{\varphi }_k+\mathrm{s}\mathrm{i}\mathrm{n}{\gamma }_k\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_k\mathrm{c}\mathrm{o}\mathrm{s}{\varphi }_k\cdot e^{\mathrm{j}{\eta }_k}}{-\mathrm{s}\mathrm{i}\mathrm{n}{\gamma }_k\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k\cdot e^{\mathrm{j}{\eta }_k}}=

\frac{\mathrm{s}\mathrm{i}\mathrm{n}{\varphi }_k}{\mathrm{t}\mathrm{a}\mathrm{n}{\gamma }_k\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k\cdot e^{\mathrm{j}{\eta }_k}}-\frac{\mathrm{c}\mathrm{o}\mathrm{s}{\varphi }_k}{\mathrm{t}\mathrm{a}\mathrm{n}{\theta }_k}（20）

r_k=\frac{\mathrm{c}\mathrm{o}\mathrm{s}{\gamma }_k\mathrm{c}\mathrm{o}\mathrm{s}{\varphi }_k+\mathrm{s}\mathrm{i}\mathrm{n}{\gamma }_k\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_k\mathrm{s}\mathrm{i}\mathrm{n}{\varphi }_k\cdot e^{\mathrm{j}{\eta }_k}}{-\mathrm{s}\mathrm{i}\mathrm{n}{\gamma }_k\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k\cdot e^{\mathrm{j}{\eta }_k}}=

-\frac{\mathrm{c}\mathrm{o}\mathrm{s}{\varphi }_k}{\mathrm{t}\mathrm{a}\mathrm{n}{\gamma }_k\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k\cdot e^{\mathrm{j}{\eta }_k}}-\frac{\mathrm{s}\mathrm{i}\mathrm{n}{\varphi }_k}{\mathrm{t}\mathrm{a}\mathrm{n}{\theta }_k}（21）

与上述求解过程类似，可以得到下式：

\left\{\begin{array}{c}{\mathit{E}}_x={\mathit{A}}_x\mathit{T};\\ {\mathit{E}}_y={\mathit{A}}_y\mathit{T};\\ {\mathit{E}}_z={\mathit{A}}_z\mathit{T}.\end{array}\right.（22）

其中：

\left\{\begin{array}{l}{\mathit{E}}_x={\mathit{E}}_z{\mathit{\Psi }}_{\mathrm{2}};\\ {\mathit{A}}_x={\mathit{A}}_z{\mathit{\Phi }}_{\mathrm{2}}.\end{array}\right.（23）

\left\{\begin{array}{l}{\mathit{E}}_y={\mathit{E}}_z{\mathit{\Psi }}_{\mathrm{3}};\\ {\mathit{A}}_y={\mathit{A}}_z{\mathit{\Phi }}_{\mathrm{3}}.\end{array}\right.（24）

同样可得表达式{\mathit{\Phi }}_{\mathrm{2}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left[\begin{array}{cccc}{p}_{\mathrm{1}}& p_{\mathrm{2}}& \cdots & p_K\end{array}\right]；{\mathit{\Phi }}_{\mathrm{3}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left[\begin{array}{cccc}{r}_{\mathrm{1}}& r_{\mathrm{2}}& \cdots & r_K\end{array}\right]。{\mathit{\Psi }}_{\mathrm{2}}和{\mathit{\Psi }}_{\mathrm{3}}的表达式为

{\mathit{\Psi }}_{\mathrm{2}}={\left({\mathit{E}}_z^{\mathrm{H}}{\mathit{E}}_z\right)}^{-\mathrm{1}}{\mathit{E}}_z^{\mathrm{H}}{\mathit{E}}_x={\mathit{E}}_z^{†}{\mathit{E}}_x（25）

{\mathit{\Psi }}_{\mathrm{3}}={\left({\mathit{E}}_z^{\mathrm{H}}{\mathit{E}}_z\right)}^{-\mathrm{1}}{\mathit{E}}_z^{\mathrm{H}}{\mathit{E}}_y={\mathit{E}}_z^{†}{\mathit{E}}_y（26）

根据p_k和r_k的特点，可以得到：

\frac{\mathrm{I}\mathrm{m}\left(p_k\right)}{\mathrm{I}\mathrm{m}\left(r_k\right)}=-\mathrm{t}\mathrm{a}\mathrm{n}{\varphi }_k（27）

因此信号的方位角为

{\varphi }_k=\mathrm{t}\mathrm{a}{\mathrm{n}}^{-\mathrm{1}}\left[-\frac{\mathrm{I}\mathrm{m}\left(p_k\right)}{\mathrm{I}\mathrm{m}\left(r_k\right)}\right]（28）

由于方位角{\varphi }_k\in \left(-\mathrm{\pi },\mathrm{\pi }\right)，其存在象限模糊的情况，为了确保正确求解，需要对{\varphi }_k进行验证。分别将p_k和r_k与\mathrm{c}\mathrm{o}\mathrm{s}{\varphi }_k和\mathrm{s}\mathrm{i}\mathrm{n}{\varphi }_k相乘，再将它们相加，得到等式：

p_k\mathrm{c}\mathrm{o}\mathrm{s}{\varphi }_k+r_k\mathrm{s}\mathrm{i}\mathrm{n}{\varphi }_k=-\frac{\mathrm{1}}{\mathrm{t}\mathrm{a}\mathrm{n}{\theta }_k}（29）

将{\theta }_k代入式（29），令\Delta =p_k\mathrm{c}\mathrm{o}\mathrm{s}{\varphi }_k+r_k\mathrm{s}\mathrm{i}\mathrm{n}{\varphi }_k+\frac{\mathrm{1}}{\mathrm{t}\mathrm{a}\mathrm{n}{\theta }_k}，根据Δ({\varphi }_k)与Δ({\varphi }_k+π)的绝对值大小比较来确定{\varphi }_k的取值具体判别方法如下：

为了分离出信号的极化参数信息，算法继续对p_k和r_k进行处理，进而得到极化参数的表征。由p_k的表达形式可得：

\mathrm{t}\mathrm{a}\mathrm{n}{\gamma }_k\cdot e^{\mathrm{j}{\eta }_k}=\frac{\mathrm{s}\mathrm{i}\mathrm{n}{\varphi }_k}{\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k\left(p_k+\frac{\mathrm{c}\mathrm{o}\mathrm{s}{\varphi }_k}{\mathrm{t}\mathrm{a}\mathrm{n}{\theta }_k}\right)}=

\frac{\mathrm{s}\mathrm{i}\mathrm{n}{\varphi }_k\mathrm{t}\mathrm{a}\mathrm{n}{\theta }_k}{\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k\left(p_k\mathrm{t}\mathrm{a}\mathrm{n}{\theta }_k+\mathrm{c}\mathrm{o}\mathrm{s}{\varphi }_k\right)}（30）

因此：

{\gamma }_{k\mathrm{1}}=\mathrm{t}\mathrm{a}{\mathrm{n}}^{-\mathrm{1}}\left[\left|\frac{\mathrm{s}\mathrm{i}\mathrm{n}{\varphi }_k\mathrm{t}\mathrm{a}\mathrm{n}{\theta }_k}{\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k\left(p_k\mathrm{t}\mathrm{a}\mathrm{n}{\theta }_k+\mathrm{c}\mathrm{o}\mathrm{s}{\varphi }_k\right)}\right|\right]（31）

{\eta }_{k\mathrm{1}}=\mathrm{a}\mathrm{r}\mathrm{g}\left[\frac{\mathrm{s}\mathrm{i}\mathrm{n}{\varphi }_k\mathrm{t}\mathrm{a}\mathrm{n}{\theta }_k}{\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k\left(p_k\mathrm{t}\mathrm{a}\mathrm{n}{\theta }_k+\mathrm{c}\mathrm{o}\mathrm{s}{\varphi }_k\right)}\right]（32）

由r_k的表达式可得：

\mathrm{t}\mathrm{a}\mathrm{n}{\gamma }_k\cdot e^{\mathrm{j}{\eta }_k}=\frac{-\mathrm{c}\mathrm{o}\mathrm{s}{\varphi }_k}{\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k\left(r_k+\frac{\mathrm{s}\mathrm{i}\mathrm{n}{\varphi }_k}{\mathrm{t}\mathrm{a}\mathrm{n}{\theta }_k}\right)}=

\frac{-\mathrm{c}\mathrm{o}\mathrm{s}{\varphi }_k\mathrm{t}\mathrm{a}\mathrm{n}{\theta }_k}{\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k\left(r_k\mathrm{t}\mathrm{a}\mathrm{n}{\theta }_k+\mathrm{s}\mathrm{i}\mathrm{n}{\varphi }_k\right)}（33）

因此：

{\gamma }_{k\mathrm{2}}=\mathrm{t}\mathrm{a}{\mathrm{n}}^{-\mathrm{1}}\left[\left|\frac{-\mathrm{c}\mathrm{o}\mathrm{s}{\varphi }_k\mathrm{t}\mathrm{a}\mathrm{n}{\theta }_k}{\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k\left(r_k\mathrm{t}\mathrm{a}\mathrm{n}{\theta }_k+\mathrm{s}\mathrm{i}\mathrm{n}{\varphi }_k\right)}\right|\right]（34）

{\eta }_{k\mathrm{2}}=\mathrm{a}\mathrm{r}\mathrm{g}\left[\frac{-\mathrm{c}\mathrm{o}\mathrm{s}{\varphi }_k\mathrm{t}\mathrm{a}\mathrm{n}{\theta }_k}{\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k\left(r_k\mathrm{t}\mathrm{a}\mathrm{n}{\theta }_k+\mathrm{s}\mathrm{i}\mathrm{n}{\varphi }_k\right)}\right]（35）

将极化参数的两种表达形式进行平均，得到最终结果如下：

{\gamma }_k=\frac{{\gamma }_{k_{\mathrm{1}}}+{\gamma }_{k_{\mathrm{2}}}}{\mathrm{2}}（36）

{\eta }_k=\frac{{\eta }_{k_{\mathrm{1}}}+{\eta }_{k_{\mathrm{2}}}}{\mathrm{2}}（37）

由于{\gamma }_k\in \left(\mathrm{0},\frac{\mathrm{\pi }}{\mathrm{2}}\right)，极化参数的求解过程不存在象限模糊问题。

对于多个入射信号，需要确定{\mathit{\Psi }}_{\mathrm{1}}、{\mathit{\Psi }}_{\mathrm{2}}和{\mathit{\Psi }}_{\mathrm{3}}的特征值分组，即需要确定{\mathit{\Psi }}_{\mathrm{1}}的哪些特征值对应于{\mathit{\Psi }}_{\mathrm{2}}、{\mathit{\Psi }}_{\mathrm{3}}的哪些特征值。由于MUSIC谱搜索法的计算量较大，采取另一种配对特征值的方法。为了将3个矩阵的特征值正确配对，可以先将{\mathit{\Psi }}_{\mathrm{1}}与{\mathit{\Psi }}_{\mathrm{2}}配对，再将{\mathit{\Psi }}_{\mathrm{2}}与{\mathit{\Psi }}_{\mathrm{3}}配对，有：

{\mathit{\Psi }}_{\mathrm{12}}={\mathit{\Psi }}_{\mathrm{1}}^{-\mathrm{1}}{\mathit{\Psi }}_{\mathrm{2}}={\mathit{T}}^{-\mathrm{1}}{\mathit{\Phi }}_{\mathrm{1}}^{-\mathrm{1}}{\mathit{\Phi }}_{\mathrm{2}}\mathit{T}（38）

式中，{\mathit{\Psi }}_{\mathrm{12}}的特征值{\overline{q}}_k对应于{\mathit{\Psi }}_{\mathrm{2}}的特征值与{\mathit{\Psi }}_{\mathrm{1}}的特征值的比率。因此，对于k=\mathrm{1,2},\cdots ,K，{\overline{q}}_k的值对应于{\mathit{\Psi }}_{\mathrm{2}}的特征值p_{k_{\mathrm{1}}}与{\mathit{\Psi }}_{\mathrm{1}}的特征值q_{k_{\mathrm{2}}}比值的所有集合中与{\overline{q}}_k相差最小的那一个配对情况：

\left\{\left|\frac{p_{k_{\mathrm{1}}}}{q_{k_{\mathrm{2}}}}-{\overline{q}}_k\right|\right\}

相似的，将{\mathit{\Psi }}_{\mathrm{2}}与{\mathit{\Psi }}_{\mathrm{3}}配对，有：

{\mathit{\Psi }}_{\mathrm{23}}={\mathit{\Psi }}_{\mathrm{2}}^{-\mathrm{1}}{\mathit{\Psi }}_{\mathrm{3}}={\mathit{T}}^{-\mathrm{1}}{\mathit{\Phi }}_{\mathrm{2}}^{-\mathrm{1}}{\mathit{\Phi }}_{\mathrm{3}}\mathit{T}（39）

得到集合：

\left\{\left|\frac{r_{k_{\mathrm{1}}}}{p_{k_{\mathrm{2}}}}-{\overline{p}}_k\right|\right\}

其中：{\overline{p}}_k为{\mathit{\Psi }}_{\mathrm{23}}的特征值；k_{\mathrm{1}},k_{\mathrm{2}}=\mathrm{1,2},\cdots ,K。只要\mathit{A}、{\mathit{A}}_{q\mathrm{1}}、{\mathit{A}}_{q\mathrm{2}}、{\mathit{A}}_x、{\mathit{A}}_y和{\mathit{A}}_z是列满秩矩阵，该方法便可以对信号的方向信息和极化信息进行联合估计。避免了谱峰搜索和迭代操作，显著降低计算量。同时，由于信号的俯仰角和方位角是根据不同的阵列流型计算得出的，信号的方向信息避免了误差传递。

3 计算步骤及复杂度分析

将二维DOA估计过程简化，算法步骤如下。

步骤1：根据阵列接收矢量\mathit{z}\left(t\right)计算快拍数为N次情况下的样本协方差矩阵\mathit{R}。

步骤2：对\mathit{R}进行特征值分解，得到信号子空间{\mathit{E}}_s。

步骤3：根据式（16）和式（22）将{\mathit{E}}_s分解为不同的子阵。

步骤4：根据式（18）、式（25）和式（26）得到矩阵{\mathit{\Psi }}_{\mathrm{1}}、{\mathit{\Psi }}_{\mathrm{2}}和{\mathit{\Psi }}_{\mathrm{3}}。

步骤5：对{\mathit{\Psi }}_{\mathrm{1}}、{\mathit{\Psi }}_{\mathrm{2}}和{\mathit{\Psi }}_{\mathrm{3}}进行特征值分解，并按照式（38）、式（39）进行特征值配对。

步骤6：按照式（19）、式（28）、式（36）和式（37）求解得到参数估计值\left({\theta }_k,{\varphi }_k,{\gamma }_k,{\eta }_k\right)。

运算量分析：在步骤1中，构造协方差矩阵的计算量为O\left\{N{\left(\mathrm{3}L\right)}^{\mathrm{2}}\right\}；步骤2中特征值分解的计算量为O\left\{{\left(\mathrm{3}L\right)}^{\mathrm{3}}\right\}；步骤4中的3次广义逆运算的计算量为O\left\{\mathrm{3}\left(L-\mathrm{1}\right)K^{\mathrm{2}}+\mathrm{2}L{K}^{\mathrm{2}}\right\}，其中的核心复乘运算量级为O\left\{\mathrm{3}\left(L-\mathrm{1}\right)K^{\mathrm{2}}+\mathrm{2}L{K}^{\mathrm{2}}\right\}；步骤5中的3次特征值分解计算量为O\left\{\mathrm{3}{K}^{\mathrm{3}}\right\}，2次参数配对的最大计算量为O\left\{\mathrm{2}{K}^{\mathrm{2}}\right\}；其他线性运算的计算量忽略不计。因此本文算法总计算量O\left\{\left(\mathrm{10}L-\mathrm{4}\right)K^{\mathrm{2}}+\mathrm{3}{K}^{\mathrm{3}}+\mathrm{9}N{L}^{\mathrm{2}}+\mathrm{27}{L}^{\mathrm{3}}\right\}。RD-MUSIC算法主要涉及构建协方差矩阵、特征值分解和2个二维搜索的过程，因此，RD-MUSIC算法计算量O\left\{\mathrm{9}N{L}^{\mathrm{2}}+\mathrm{27}{L}^{\mathrm{3}}+\left(\mathrm{9}{L}^{\mathrm{2}}+\mathrm{3}L\right)n^{\mathrm{2}}\right\}，其中n是搜索点数，由于搜索点数较多，且估计精度受搜索步长影响，因此本文算法计算复杂度远小于RD-MUSIC算法。

4 仿真实验

假设阵列排布方式与图1相同，线性极化敏感阵列沿Z轴分布，阵元数量为8，每个阵元均包含3个偶极子，阵元间距为半波长。入射信号为远场窄带信号，信号之间相互独立，噪声为高斯白噪声，信号与噪声相互独立。仿真实验硬件平台为装配Intel（R） Core（TM） i7-13650HX 处理器、主频2.60 GHz、内存24 GB的计算机，仿真软件为Matlab 2024b。

4.1 算法有效性

为验证算法有效性及参数配对性能，假设3个信号入射到极化敏感阵列，信号的俯仰角\theta分别为\mathrm{10}°,\mathrm{20}°,\mathrm{60}°；信号方位角\varphi分别为\mathrm{60}°,\mathrm{150}°,-\mathrm{30}°；信号的极化幅度角\gamma分别为\mathrm{50}°,\mathrm{20}°,\mathrm{30}°；信号的极化相位角\eta分别为\mathrm{110}°,\mathrm{60}°,\mathrm{150}°。快拍数取5 000，信噪比（SNR）取值为20 dB，蒙特卡洛试验次数为100。图2、图3为算法参数估计的星座图。可以看出，该算法可以正确估计出3个入射信号的二维DOA和极化参数，且配对正确。

4.2 计算量对比

为验证本文算法的计算效率，表2展示了3种算法的运行时间随阵元数变化的对比情况。这里假设信源数量为3，RD-MUSIC算法的搜索间隔为1°，其他实验条件不变，运行时间为100次蒙特卡洛实验的平均值。从表1可见，本文算法运行时间远低于RD-MUSIC算法，这是因为其他变量均远小于搜索点数n，同时，RD-MUSIC的估计精准度依赖n的大小。相较于Ratio-ESPRIT算法，本文算法利用了三偶极子阵元，且采用了多重旋转不变性，求逆运算的维度和次数更多，因此计算运行时间高于Ratio-ESPRIT算法。

4.3 信噪比对算法的影响

为了验证信噪比对算法的影响，假设一个入射信号，参数为\left(\theta ,\varphi ,\gamma ,\eta \right)=\left(\mathrm{10}°,\mathrm{70}°,\mathrm{20}°,\mathrm{60}°\right)，信噪比取值范围为0\sim20 dB，每间隔2 dB进行递增，快拍数取5 000，蒙特卡洛试验次数为500。考察本文方法与RD-MUSIC和Ratio-ESPRIT算法的对比，且保持每种算法的天线数目、阵元间距以及其他仿真条件均相同，RD-MUSIC算法的搜索步长设为0.1°，参数估计精度使用均方根误差（Root Mean Square Error, RMSE）。所得实验结果如图4和图5，克拉美劳下界(Cramer-Rao Lower Bound, CRLB)套样是理论最小RSME。可以看出，本文DOA估计精度优于RD-MUSIC算法和Ratio-ESPRIT算法；极化幅度角估计精度不及另外两种算法，这是因为本文算法的极化参数精度受DOA参数估计精度的影响，导致误差传递，相比于RD-MUSIC算法，牺牲了高精度，换取了较低的计算复杂度。

4.4 快拍数对算法的影响

为了验证快拍数对算法的影响，保持信号数量和信号的DOA与极化参数不变，信噪比取值为20 dB，快拍数取50\sim5 000，蒙特卡洛试验次数为500。通过比较3种算法的RMSE，从图6和图7中可以得到类似的结论。

在RD-MUSIC算法搜索步长为0.1°的条件下，本文算法的DOA估计精度优于RD-MUSIC与Ratio-ESPRIT算法，更加接近CRLB下界。极化幅度角的估计精度在快拍数为500\sim4 000时本文算法占优，其他情况时RD-MUSIC算法精度更高，极化相位角虽然均劣于RD-MUSIC算法，但3种算法的极化相位角估计精度相当，性能基本一致。

4.5 算法分辨力及配对正确率

为了验证算法的分辨率，假设信号数量为2，两个信号的方位角相同，改变俯仰角的差值，信号的参数\left(\theta ,\varphi ,\gamma ,\eta \right)分别为\left(\mathrm{10}°,\mathrm{70}°,\mathrm{20}°,\mathrm{60}°\right)、

\left(\left(\mathrm{10}+\mathrm{\Delta }i\right)°,\mathrm{70}°,\mathrm{10}°,\mathrm{150}°\right)，\mathrm{\Delta }i=\mathrm{1},\mathrm{2},\cdots ,\mathrm{10}。从图8看出，随着2个信号的俯仰角差值的增加，每个参数的估计误差逐渐降低，当角度差值增加到10°时，所有参数的均方根误差基本控制在1°以内。

其他参数不变，改变2个信号的方位角，观察不同俯仰角差值在不同的方位角时的参数配对正确率情况，如图9所示。为充分体现算法的性能，实验展示了任意方位角条件下的配对正确率情况，由于在判断方位角时存在象限模糊问题，且-90°、0°和90°为象限跳变处，因此2个信号更易受误差影响，配对成功率较低，其他角度均展现出较高正确率。

5 结束语

本文提出了一种基于ESPRIT的二维DOA参数与极化参数联合估计算法。通过子阵的不同划分方式，算法实现了DOA参数与极化参数的闭式计算，规避了搜索过程，并且仿真实验验证了该方法在二维DOA和极化参数联合估计中的高精度特性。相较于传统MUSIC和ESPRIT算法，所提方法采用均匀线性阵列，通过改进信号子空间构造方式，在节约布阵空间资源的同时，规避了全局谱搜索过程，显著降低了计算复杂度。与对比算法的实验结果证明，本文算法的DOA和极化相位角参数估计精度均优于其他对比算法，且计算效率远优于RD-MUSIC算法。

参考文献

[1] CHEN H， FANG J X， WANG W L， et al. Near-field target localization for EMVS-MIMO radar with arbitrary configuration［J］. IEEE Transactions on Aerospace and Electronic Systems， 2024，60（4）：5406-5417.

[2] WEN F Q， SHI J P， ZHANG Z J. Joint 2D-DOD， 2D-DOA， and polarization angles estimation for bistatic EMVS-MIMO radar via PARAFAC analysis［J］. IEEE Transactions on Vehicular Technology， 2020，69（2）：1626-1638.

[3] FLORES L A， LOMAS I， GUACHALÁ L， et al. Direction-of-arrival （DOA） estimation based on real field measurements and modified linear regression［J］. Engineering Proceedings， 2024，77（1）：No.11.

[4] CHEN H， LIN H G， LIU W， et al. Augmented multi-subarray dilated nested array with enhanced degrees of freedom and reduced mutual coupling［J］. IEEE Transactions on Signal Processing， 2024，72：1387-1399.

[5] FANG X Y， LI M M， HAN J Z， et al. Accurate direction-of-arrival estimation method based on space-time modulated metasurface［J］. IEEE Transactions on Antennas and Propagation， 2022，70（11）：10951-10964.

[6] SHIIBA H， KUSE R， SANO M， et al. Performance enhancement of direction of arrival estimation using a single port metasurface reflector antenna［C］//Proceedings of the 2024 IEEE International Symposium on Antennas and Propagation and INC/USNC-URSI Radio Science Meeting. Piscataway， USA： IEEE， 2024：2681-2682.

[7] DAI M H， MA X F， LIU W， et al. Monopulse based DOA and polarization estimation with polarization sensitive arrays［J］. Signal Processing， 2022，195：No.108487.

[8] YANG Y L， JIANG G J. Efficient DOA and polarization estimation for dual-polarization synthetic nested arrays［J］. IEEE Systems Journal， 2022，16（4）：6277-6288.

[9] 李金朋，赵闯，胡德秀，等.利用Lissajous图形的雷达辐射源DOA与极化参数联合估计方法［J］. 信息工程大学学报，2025，26（4）：415-422.

[10] WAGNER M， PARK Y， GERSTOFT P. Gridless DOA estimation and root-MUSIC for non-uniform linear arrays［J］. IEEE Transactions on Signal Processing， 2021，69：2144-2157.

[11] ZHU H G， FENG W K， FENG C Q， et al. Deep unfolded gridless DOA estimation networks based on atomic norm minimization［J］. Remote Sensing， 2024，15（1）：No.13.

[12] 张涛，孙昭乾，郭沐然.基于矩阵分布式重构的正交极化阵DOA估计［J］.信号处理，2024，40（10）：1822-1833.

[13] GUO R， MAO X， LI S， et al. A fast DOA estimation algorithm based on polarization MUSIC［J］. Radioengineering， 2015，24（1）：214-225.

[14] ZHANG Z Q， LAN X， WANG X P. A fast method for multidimensional joint parameter estimation of polarization-sensitive arrays［J］. Sensors， 2023，23（19）：No.8193.

[15] 马慧慧， 陶海红. 稀疏拉伸式L型极化敏感阵列的二维波达方向和极化参数联合估计［J］. 电子与信息学报， 2020，42（4）：902-909.

[16] WANG X D， ZHANG X F， LI J F， et al. Improved ESPRIT method for joint direction-of-arrival and frequency estimation using multiple-delay output［J］. International Journal of Antennas and Propagation， 2012：No.309269.

[17] ZHANG X F， CHEN C， LI J F， et al. Blind DOA and polarization estimation for polarization-sensitive array using dimension reduction MUSIC［J］. Multidimensional Systems and Signal Processing， 2014，25（1）：67-82.

[18] 李会勇，张远芳，谢菊兰.极化敏感线阵的模值约束降维Root-MUSIC算法［J］.信号处理， 2016，32（2）：173-178．

[19] 逯岩斌，陈文东，杨赟秀，等.基于互质极化敏感阵列的参数降维估计算法［J］.电子测量技术，2024，47（5）：173-180.

[20] LAN X， LIU W， NGAN H Y T. Joint 4-D DOA and polarization estimation based on linear tripole arrays［C］//Proceedings of the 2017 22nd International Conference on Digital Signal Processing. Piscataway， USA： IEEE， 2017：1-5.

[21] 刘鲁涛，王传宇. 基于极化敏感阵列均匀线阵的二维DOA估计［J］.电子与信息学报，2019，41（10）：2350-2357.