极化敏感阵列方位依赖误差校正算法

刘鲁涛; 赵梓君; 李利; LIU Lutao; ZHAO Zijun; LI Li

极化敏感阵列方位依赖误差校正算法

doi: 10.11887/j.cn.202406019

刘鲁涛¹ ，赵梓君¹ ，李利²

1. 哈尔滨工程大学信息与通信工程学院，黑龙江哈尔滨 150001

2. 海军大连舰艇学院信息系统系，辽宁大连 116018

基金项目: 航空科学基金资助项目（201901012005）

详细信息

作者简介

刘鲁涛（1977—），男，黑龙江哈尔滨人，副教授，博士，博士生导师，E-mail: liulutao@hrbeu.edu.cn

中图分类号: TN911.7

文献标识码: A

文章编号: 1001-2486(2024)06-174-10

Orientation-dependent error calibration algorithm for polarization sensitive array

LIU Lutao¹ ， ZHAO Zijun¹ ， LI Li²

1. College of Information and Communication Engineering, Harbin Engineering University, Harbin 150001 , China

2. Department of Information System, Dalian Naval Academy, Dalian 116018 , China

摘要

在实际应用中多种类型阵列误差同时存在，针对这种情况下阵列误差方位依赖的特点，提出了一种基于流形分离技术(manifold separation technique, MST)的改进多重信号分类(multiple signal classification, MUSIC)算法，可以有效解决多种阵列误差影响下的波达方向估计问题。利用MST获得包含阵列非理想特性的采样矩阵，从而进行精准测向；通过二维傅里叶变换求解二维空间谱，与现有MUSIC校正算法相比，减少了谱峰搜索的运算量。理论分析和仿真验证了该算法的有效性，可为实际问题的解决提供参考。

关键词

极化敏感阵列 / 流形分离技术 / 阵列误差校正 / 波达方向估计

Abstract

In practical applications, multiple types of array errors exist simultaneously. In view of the orientation dependence of array errors in this case, an improved MUSIC(multiple signal classification) algorithm based on MST(manifold separation technique) was proposed, which can effectively solve DOA(direction of arrival) estimation problem under the influence of multiple array errors. The sampling matrix, which contains all the non-ideal characteristics of the array was obtained by using MST, so as to achieve accurate direction finding. By using two-dimensional Fourier transform, the 2D spatial spectrum was solved. Compared with existing MUSIC algorithm, the computation amount of spectrum peak searching was reduced. Theoretical analysis and simulation verify the effectiveness of the algorithm, which can provide reference for the solution of practical problems.

Keywords

polarization sensitive array / manifold separation technique / array error calibration / direction of arrival estimation

1 模型的建立 1.1 极化敏感阵列模型的建立 1.2 误差模型的建立 1.2.1 通道幅相误差 1.2.2 阵元位置误差 1.2.3 阵元互耦误差 1.2.4 阵元方向图误差 1.3 MST模型的建立 2 基于MST的DOA估计方法 2.1 2D-PES-MUSIC算法 2.2 改进的2D-FFT-MUSIC算法 3 仿真实验及分析 3.1 MST重构精度分析 3.2 空域-极化域联合估计精度分析 3.2.1 单信源估计精度分析 3.2.2 双信源估计精度分析 3.3 运算复杂度对比 4 结论

极化敏感阵列（polarization sensitive array，PSA）是利用内部阵元极化选择特性来获取信号的，在雷达、声呐、通信、导航、地震勘探、射电天文、工业探伤、医学成像等众多领域引起了广泛关注^[1]。实际的工程应用中，阵列不能满足理想的条件，不可避免地会存在阵元幅度相位误差、阵元位置误差等影响，这就导致传统测向算法的性能严重恶化，因而研究阵列存在误差条件下的波达方向（direction of arrival，DOA）估计技术是十分必要的。误差的校正技术主要分为两类：有源校正^[2-4]和自校正^[5-7]。有源校正是指通过引入方位精确已知的辅助信源来求解误差参数，在已知误差参数的前提下进行精准测向，这种算法的优点是计算的复杂度较低，但在引入辅助源信息时会产生新的误差；自校正则是利用某些原理构造代价函数，通过迭代的方式对DOA信息和阵列误差进行联合校正，这种算法虽然不会引入新的误差，但运算复杂度较高。上述方法利用误差参数对阵列扰动进行建模，将阵列误差校正转化为参数估计问题，然而，在实际应用中，极化天线单元并不是全向天线，接收机对各个方向来波信号的增益不同且未知，故实际的综合误差矩阵无法根据误差参数写出明确的解析表达式。

流形分离技术（manifold separation technique，MST）是将波形域建模的方式运用到阵列接收数据建模的一种新方法，即将阵列的接收数据中与信号有关的部分表示为阵列对入射信号波长的采样^[8]，故可以利用多重信号分类（multiple signal classification，MUSIC）算法对存在误差的阵列进行精准测向。文献^[9]将信号的导向矢量分解成一个采样矩阵和一个具有Vandermonde结构矢量的乘积形式，进而利用root-MUSIC算法对空域信息进行估计。文献^[10]在此基础上进一步对方位角和俯仰角联合估计的问题进行了研究，并提出了一种基于二维傅里叶变换的MUSIC算法，避免了谱峰搜索的过程，大大减少了运算量。文献^[8]对基于MST的DOA和极化参数的联合估计算法进行了研究，并针对相干信源的测向提出了极化修正变量投影信号子空间拟合（modified variable projection-signal subspace fitting，MVP-SSF）算法。文献^[11]重点研究了基于流形分离技术的相干源的DOA-极化联合估计，提出基于流形分离的广义子空间拟合算法对空间角和极化信息进行估计，并针对广义子空间拟合算法运算量大的问题，研究了斜投影算法对极化信息进行估计。

本文利用流形分离技术将入射信号的导向矢量进行分解，在未知阵列结构实际参数的情况下，通过实测数据获得采样矩阵；在计算二维空间谱时，提出了一种改进的二维快速傅里叶变换的多重信号分类（2 dimensional fast Fourier transformation-multiple signal classification，2D-FFT-MUSIC）算法，无须计算系数矩阵，大大减少了运算量。该算法适用于任意配置的阵列，可实现存在阵列误差时的精准测向，具有重要的实际应用价值。

1 模型的建立

1.1 极化敏感阵列模型的建立

分布式极化敏感阵列的摆放示意如图1所示，阵列中各阵元均是单偶极子，假设各阵元的波束指向一致。

假设有M个远场窄带信号入射至该N元阵列，其中第m个信源的空域信息和极化域信息为（θ_m，φ_m，γ_m，η_m），θ和φ分别为电磁波的方位角和俯仰角，γ和η分别为信号的极化辅助角和极化相位差，则极化敏感阵列接收模型可表示为：

\begin{matrix} X (t) = A s (t) + n (t) \\ = \sum_{m = 1}^{M} a (θ_{m}, ϕ_{m}, γ_{m}, η_{m}) s_{m} (t) + n (t) \end{matrix}

(1)

图1极化敏感阵列空间摆放示意图

Fig.1Spatial layout diagram of polarization sensitive array

式中：s_m（t）为第m个远场窄带信号；n（t）为高斯白噪声矩阵；A=[a（θ₁，φ₁，γ₁，η₁）··· a（θ_M，φ_M，γ_M，η_M）]为空域-极化域联合导向矢量矩阵，其中a（θ_m，φ_m，γ_m，η_m）为第m个信号的导向矢量，构造如式（2）所示。

\begin{matrix} a (θ_{m}, ϕ_{m}, γ_{m}, η_{m}) \\ = u (θ_{m}, ϕ_{m}) \cdot B e t a \cdot Ψ (θ_{m}, ϕ_{m}, γ_{m}, η_{m}) \\ = [\begin{matrix} a_{h} (θ_{m}, ϕ_{m}) & a_{v} (θ_{m}, ϕ_{m}) \end{matrix}] h (γ_{m}, η_{m}) \end{matrix}

(2)

式中：

u (θ_{m} ， ϕ_{m}) = d i a g (u_{1} (θ_{m} ， ϕ_{m}) ， \dots ， u_{N} (θ_{m} ， ϕ_{m}))

为第m个信号的空间相位因子，

u_{n} (θ_{m} ， ϕ_{m}) = e x p [- j ω_{m} (x_{n} c o s θ_{m} s i n ϕ_{m} + y_{n} s i n θ_{m} s i n ϕ_{m} + z_{n} c o s ϕ_{m})]

为第m个信号在n个阵元处形成的空间延迟，

ω_{m}

为第m个信号的频率，

(x_{n} ， y_{n} ， z_{n})

为第n个阵元坐标。

Beta = {[\begin{matrix} G_{1}^{sd} \dots G_{N}^{sd} \end{matrix}]}^{T} \oplus {[\begin{matrix} b_{1} \dots b_{N} \end{matrix}]}^{T}

为极化敏感矩阵，其中

\oplus

Hadamard积；

G_{n}^{s d}

为信号和第n个阵元极化完全匹配时的增益，若所有天线单元均为全向天线，且对任意方位来波信号的增益一致，同时接收机不存在失配现象时，

G_{n}^{sd} = 1; b_{n} = [c o s α_{n}^{（ 1 ）} s i n α_{n}^{（ 2 ）} ， s i n α_{n}^{（ 1 ）} s i n α_{n}^{（ 2 ）} ， c o s α_{n}^{（ 2 ）} ， 0_{1 \times 3}]

为第n个阵元在三个方向上电场和磁场的分量，

(α_{n}^{（ 1 ）} ， α_{n}^{（ 2 ）})

为第n个阵元的摆放角，示意图如图2所示；对于图1所示阵列，各阵元摆方均在XOY平面内，满足

z_{n} = 0 （ n = 1 ， \dots ， N ）

。

a_{h} (θ_{m} ， ϕ_{m}) 和 a_{v} (θ_{m} ， ϕ_{m})

分别表示水平极化分量和垂直极化分量对应的导向矢量。

图2阵元摆放角示意图

Fig.2Diagram of array elements′ placement angle

Ψ (θ_{m} ， ϕ_{m} ， γ_{m} ， η_{m})

为第m个信号的极化-角度域导向矢量，表示为：

\begin{matrix} Ψ (θ_{m}, ϕ_{m}, γ_{m}, η_{m}) \\ = [\begin{matrix} - s i n θ_{m} & c o s θ_{m} c o s ϕ_{m} \\ c o s θ_{m} & s i n θ_{m} c o s ϕ_{m} \end{matrix}] [\begin{matrix} c o s γ_{m} \\ s i n γ_{m} e^{j η_{m}} \end{matrix}] \\ = Ξ (θ_{m}, ϕ_{m}) \cdot h (γ_{m}, η_{m}) \end{matrix}

(3)

矢量

Ψ (θ_{m} ， ϕ_{m} ， γ_{m} ， η_{m})

描述了信号的角度域和极化域的相干结构，与矢量阵元的空间位置无关；式中的h（γ_m，η_m）为第m个信号的极化矢量，与信号的极化状态有关，但与阵列的空间结构与响应特性无关。

1.2 误差模型的建立

实际应用中的误差形式并不单一，阵列接收数据包含了多种误差的影响，包括通道幅相误差、阵元位置误差、阵元互耦误差和阵元方向图误差等因素。

1.2.1 通道幅相误差

阵列的通道幅相误差是一种与信号空间方位无关的复增益误差，但是与系统的工作频率有关。对于窄带信号而言，由于窄带信号随时间是缓变的，可以假设工作频带内的幅相误差是稳定的，则幅相误差矩阵可表示为：

\begin{matrix} Γ = d i a g (Γ_{1}, Γ_{2}, \dots, Γ_{N}) \\ = d i a g (1, ρ_{2} e x p (j β_{2}), \dots, ρ_{N} e x p (j β_{N})) \end{matrix}

(4)

式中，ρ_n和β_n分别为幅度误差和相位误差。存在幅相误差的阵列流形表示为：

\tilde{A} (θ_{m}, ϕ_{m}, γ_{m}, η_{m}) = Γ \cdot A (θ_{m}, ϕ_{m}, γ_{m}, η_{m})

(5)

式中，A（θ_m，φ_m，γ_m，η_m）为理想的阵列流形矩阵。

1.2.2 阵元位置误差

在实际环境中，阵列会受到环境温度、制造和安装误差等因素的影响，导致实际阵元位置与标称值有一定偏差，这种偏差就是阵元位置误差。假设第n个阵元的位置标称值为（x_n，y_n，z_n），阵元位置偏差为（Δx_n，Δy_n，Δz_n），存在位置误差的阵列流形表示为：

\tilde{A} (θ_{m}, ϕ_{m}, γ_{m}, η_{m}) = Δ A_{p} (θ_{m}, ϕ_{m}) \cdot A (θ_{m}, ϕ_{m}, γ_{m}, η_{m})

(6)

式中，

Δ A_{p}

=diag（exp（-j2πfΔτ₁），···，exp（-j2πfΔτ_N））为阵元位置误差矩阵，其中Δτ_n=（Δx_ncosθ_msinφ_m+Δy_nsinθ_msinφ_m+Δz_ncosφ_m）是由位置误差引起的波程差。

1.2.3 阵元互耦误差

实际天线工作时，阵元间的相互影响即互耦效应是不可避免的。存在互耦误差的阵列流形表示为：

\tilde{A} (θ_{m}, ϕ_{m}, γ_{m}, η_{m}) = Z \cdot A (θ_{m}, ϕ_{m}, γ_{m}, η_{m})

(7)

式中，Ζ为反映阵元互耦效应的互耦矩阵。根据互易原理，互耦矩阵Ζ通常为一对称矩阵。

1.2.4 阵元方向图误差

实际的极化天线单元并不是全向天线，对各个方向来波信号的增益并不相同，且存在失配现象，因此G^sd无法给出确定值，极化敏感矩阵Beta也无法给出具体计算公式。方向图的不一致性可通过引入方位依赖的幅相误差来描述，此时阵列流形表示为：

\tilde{A} (θ_{m}, ϕ_{m}, γ_{m}, η_{m}) = W (θ_{m}, ϕ_{m}) \cdot A (θ_{m}, ϕ_{m}, γ_{m}, η_{m})

(8)

式中，

W (θ_{m} ， ϕ_{m})

=diag（W₁

(θ_{m} ， ϕ_{m})

，···，W_N（θ_m，φ_m）），W_n（θ_m，φ_m）为第n个阵元在方位（θ_m，φ_m）处的天线响应。

通过对上述各种典型的阵元误差的讨论和建模可知，综合多种误差的影响，可以将误差矩阵等效为一个与来波方向有关的幅相误差模型，且该综合误差矩阵无法根据误差参数进行明确的解析表示^[12]。此时真实的阵列流形矩阵可以表示为：

\tilde{A} (θ_{m}, ϕ_{m}, γ_{m}, η_{m}) = C (θ_{m}, ϕ_{m}) \cdot A (θ_{m}, ϕ_{m}, γ_{m}, η_{m})

(9)

式中，

C (θ_{m} ， ϕ_{m})

为综合误差矩阵，与信号的来波方向有关，可表示为

C (θ_{m}, ϕ_{m}) = Γ \cdot Δ A_{p} (θ_{m}, ϕ_{m}) \cdot W (θ_{m}, ϕ_{m}) \cdot Z

(10)

当系统中不存在上述的某种误差时，可将其对应的误差矩阵视为单位阵。

1.3 MST模型的建立

利用MST可以将实际的导向矢量分解为采样矩阵和两个Vandermonde结构矢量的乘积，即：

\tilde{a} (θ, ϕ) = G d (θ, ϕ) + ε

(11)

式中：ε表示模型误差；G表示仅与阵列配置有关的采样矩阵；d（θ，φ）表示仅与来波方向有关的矢量，由两个Vandermonde结构矢量构成，在二维角度搜索的情况下，d（θ，φ）可表示为

d (θ, ϕ) = d (θ) \otimes d (ϕ) \in C^{(M a \times M e) \times 1}

式中，

\otimes

表示Kronecker积运算，Ma和Me分别表示在方位角维度和俯仰角维度截断的模数，d（θ）和d（φ）具体表达式为

\{\begin{matrix} d (θ) = [e^{j \frac{M_{a} - 1}{2} θ}, e^{j \frac{M_{a} - 1}{2} (θ - Δ θ)}, e^{j \frac{M_{a} - 1}{2} (θ - 2 Δ θ)}, \dots, e^{j \frac{M_{a} - 1}{2} \cdot 0}, \\ {e^{- j \frac{M_{a} - 1}{2} \cdot Δ θ}, \dots, e^{- j \frac{M_{a} - 1}{2} (θ - Δ θ)}, e^{- j \frac{M_{a} - 1}{2} θ}]}^{T} \in C^{M_{a} \times 1} \\ d (ϕ) = [e^{\frac{M_{e - 1}}{2} ϕ}, e^{\frac{M_{e - 1}}{2} (ϕ - Δ ϕ)}, e^{\frac{M_{e - 1}}{2} (ϕ - 2 Δ ϕ)}, \dots, e^{j \frac{M_{e - 1}}{2} \cdot 0}, \\ {e^{- j \frac{M_{e - 1}}{2} \cdot Δ ϕ}, \dots, e^{- j \frac{M_{e - 1}}{2} (ϕ - Δ ϕ)}, e^{- j \frac{M_{e - 1}}{2} ϕ}]}^{T} \in C^{M_{e} \times 1} \end{matrix}

(13)

关于采样矩阵的计算，可以采用有效孔径分布函数（effective aperture distribution function，EADF）。如果阵列测量数据采用等宽采样网格，则可以通过二维傅里叶变换（2-dimensional Fourier transformation，2DFT）得到采样矩阵。

当满足θ∈[0，2π），φ∈[0，π）时，可以描述完整的球面波束图，但由于二维傅里叶变换的基函数d（θ，φ）在圆环体上而非单位圆上满足正交关系（如图3所示，圆环上的每个点表示二维傅里叶变换系数），故要将俯仰角维度的测量数据扩展至

\tilde{ϕ}

∈[0，2π），即测量到的波束图被反射在南极点的切向平面上，并在方位角上发生π角度偏移。若从矢量场来定义，补充数据的波束图是原数据波束图乘-1得到的^[13]。通过上述操作可以得到定义在圆环面上的周期波束图，如图4所示。

图32DFT拓扑结构圆环体

Fig.32DFT topological torus

图4测量数据单位圆示意图

Fig.4Diagram of unit circle of measurement data

根据图中可以看出，A（θ′，φ′）=-A（θ，φ），图中所示角度关系为：

\{\begin{matrix} θ^{'} = θ + π \\ ϕ^{'} = π - ϕ \end{matrix}

(14)

\tilde{ϕ}

∈[π，2π）时导向矢量可以通过

\tilde{ϕ}

∈[0，π）时的导向矢量在方位角进行角度π的偏移，在俯仰角进行翻转得到。即完整周期的数据导向矢量可以通过下式确定：

a_{c a l} (θ, \tilde{ϕ}) = \{\begin{matrix} a (θ, \tilde{ϕ}) \tilde{ϕ} \in [0, π), \forall θ \\ a (θ + π, 2 π - \tilde{ϕ}) \tilde{ϕ} \in [π, 2 π), \forall θ \end{matrix}

(15)

分别在θ∈[0，2π）和φ∈[0，π）取Qa和Qe个测量点，Acal_n（θ_c，φ_c）表示第n个天线的校准测量矩阵。为了使测量数据在俯仰角上以[0，2π）为周期，构造一个矩阵

A_{n}^{r}

（θ_c，φ_c）∈C^（^Qe^-2）×^Qa，通过对

A_{n}

（θ_c，φ_c）∈C^Qe^×^Qa在方位角进行角度π的偏移，在俯仰角进行翻转，再去掉第一行和最后一行以避免φ=0、φ=π处的冗余，然后组合这两个矩阵得到

A {c a l}_{n} (θ_{c}, ϕ_{c}) = [\begin{matrix} A_{n} (θ_{c}, ϕ_{c}) \\ A_{n}^{r} (θ_{c}, ϕ_{c}) \end{matrix}] \in C^{(2 Q_{e - 2) \times Q_{a}}}

(16)

第n个阵元的2D有效孔径分布函数G_n∈C^（2^Qe^-2）×^Qa对应Acal_n（θ_c，φ_c）的二维离散傅里叶反变换系数。

二维傅里叶系数具有超指数衰减的特征，导致高阶的EADF模很小，当Qa和Qe趋于无穷大时，ε可以忽略。现定义Gt_n∈C^Me^×^Ma为EADF矩阵G_n∈C^（2^Qe^-2）×^Qa的截断矩阵，由此可求得采样矩阵G为：

G = [\begin{matrix} v e c {\{G t_{1}\}}^{T} \\ ⋮ \\ v e c {\{G t_{N}\}}^{T} \end{matrix}] \in C^{N \times (M_{a} \times M_{e})}

(17)

此时阵列接收到的数据矩阵X可改写为：

X = G d (θ, ϕ) S + N

(18)

式中，S和N为信号和噪声的矩阵形式。根据接收数据获得的采样矩阵G包含了阵列误差的有关信息，可解决存在误差的阵列测向问题。

对于极化敏感阵列，利用MST将导向矢量进行正交展开，即：

\{\begin{matrix} {\tilde{a}}_{h} (θ, ϕ) = G_{h} d (θ, ϕ) \\ {\tilde{a}}_{v} (θ, ϕ) = G_{v} d (θ, ϕ) \end{matrix}

(19)

式中，G_h，G_v∈C^M^×（^Ma^×^Me^）为前述的采样矩阵，其完全由阵列结构决定，与信源DOA信息无关，故所有阵列误差的影响均包含在采样矩阵内。此时极化敏感阵列接收数据可改写为：

\begin{matrix} X = \tilde{A} (θ, ϕ, γ, η) S + N \\ = [G_{h} d (θ, ϕ) G_{v} d (θ, ϕ)] h (γ, η) S + N \\ = G [I_{2} \otimes d (θ, ϕ)] h (γ, η) S + N \end{matrix}

(20)

式中，I₂是二阶单位阵。对上式进行特征值分解，可以得到：

R = \frac{1}{K} X X^{H} = E_{s} Λ_{s} E_{s}^{H} + E_{n} Λ_{n} E_{n}^{H}

(21)

式中，

Λ_{s}

和

Λ_{n}

分别为信号子空间和噪声子空间对应的特征值，

E_{s}

和

E_{n}

为对应特征值的特征矢量，K为信号的快拍数。

2 基于MST的DOA估计方法

2.1 2D-PES-MUSIC算法

文献^[13]提出了一种基于MST的空域-极化域角度联合估计方法，首先利用MUSIC算法求出信号的DOA信息，表达式如下：

\{{\hat{θ}}_{m}, {\hat{ϕ}}_{m}\} = a r g \underset{θ_{m}, ϕ_{m}}{m a x} {d e t}^{- 1} \{[I_{2} \otimes d^{H} (θ, ϕ)] G^{H} E_{n} E_{n}^{H} .

G [I_{2} \otimes d (θ, ϕ)]\}

(22)

获得空域角度信息后，将求得的（

{\hat{θ}}_{m} ， {\hat{ϕ}}_{m}

）代入下式可得信号极化域信息：

\{{\hat{γ}}_{m}, {\hat{η}}_{m}\} = a r g \underset{θ_{m}, ϕ_{m}}{m a x} {d e t}^{- 1} \{h^{H} (γ, η) M ({\hat{θ}}_{m}, {\hat{ϕ}}_{m}) h (γ, η)

(23)

式中，矩阵

M ({\hat{θ}}_{m} ， {\hat{ϕ}}_{m}) = [I_{2} \otimes d^{H} ({\hat{θ}}_{m} ， {\hat{ϕ}}_{m})] . G^{H} E_{n} E_{n}^{H} G [I_{2} \otimes d ({\hat{θ}}_{m} ， {\hat{ϕ}}_{m})]

。

由于MUSIC算法存在的二维搜索过程将导致该算法存在大量的运算复杂度，故下节将对此进行改进。

2.2 改进的2D-FFT-MUSIC算法

对上述2D-PES-MUSIC算法进行改进，将空间谱按照傅里叶变换表达式的形式进行改写。对空域角度进行估计的表达式为：

\{{\hat{θ}}_{m}, {\hat{ϕ}}_{m}\} = a r g \underset{θ_{m}, ϕ_{m}}{m a x} {d e t}^{- 1} (M (θ, ϕ))

式中，

M （ θ ， ϕ ）

写成矩阵形式为

\begin{matrix} M (θ, ϕ) \\ = [\begin{matrix} d^{H} (θ, ϕ) & 0_{1 \times MaMe} \\ 0_{1 \times MaMe} & d^{H} (θ, ϕ) \end{matrix}] G^{H} E_{n} E_{n}^{H} G [\begin{matrix} d (θ, ϕ) & 0_{MaMe \times 1} \\ 0_{MaMex1} & d (θ, ϕ) \end{matrix}] \end{matrix}

\begin{matrix} = D^{H} (θ, ϕ) G^{H} E_{n} E_{n}^{H} G D (θ, ϕ) \\ = Q^{H} Q \end{matrix}

(25)

式中，0为零矩阵，D（θ，φ）=

[\begin{matrix} d （ θ ， ϕ ） & 0_{MaMe \times 1} \\ 0_{MaMex1} & d （ θ ， ϕ ） \end{matrix}]

，

Q = E_{n}^{H} G D （ θ ， ϕ ）

。由式（19）知采样矩阵可以写成如下的分块矩阵形式：

G = [\begin{matrix} G_{h} G_{v} \end{matrix}]

(26)

定义矩阵

V = E_{n}^{H} G

，将矩阵

V

写成分块矩阵的形式如下：

V = E_{n}^{H} \cdot [\begin{matrix} G_{h} G_{v} \end{matrix}] = [\begin{matrix} V_{1} V_{2} \end{matrix}]

(27)

式中，

V_{1}

E_{n}^{H} G_{h}

∈C^（^N-M^）×^MaMe，

V_{2}

E_{n}^{H} G_{v}

∈C^（^N-M^）×^MaMe。将矩阵

V_{1}

和

V_{2}

改写为矢量的形式，即

V_{i} = {[\begin{matrix} v_{i ， 1}^{T} v_{i ， 2}^{T} \dots v_{i ， N - M}^{T} \end{matrix}]}^{T}

，其中

v_{i ， j} \in C^{1 \times MaMe}; i = 1，2; j = 1，2 ， \dots ， N - M

。则将Q改写为如下形式：

Q = [\begin{matrix} V_{1} V_{2} \end{matrix}] \cdot [\begin{matrix} d (θ, ϕ) & 0 \\ 0 & d (θ, ϕ) \end{matrix}]

= [\begin{matrix} d^{T} (ϕ) unvec (v_{1,1}) d (θ) & d^{T} (ϕ) unvec (v_{2,1}) d (θ) \\ d^{T} (ϕ) unvec (v_{1,2}) d (θ) & d^{T} (ϕ) unvec (v_{2,2}) d (θ) \\ ⋮ & ⋮ \\ d^{T} (ϕ) unvec (v_{1, N - M}) d (θ) & d^{T} (ϕ) unvec (v_{2, N - M}) d (θ) \end{matrix}]

(28)

式中，unvec（·）为矢量矩阵化操作。注意到式（28）满足2DFT的一般表达形式，故可改写为：

Q = [\begin{matrix} 2 D F T (unvec (v_{1,1})) & 2 D F T (unvec (v_{2,1})) \\ 2 D F T (unvec (v_{1,2})) & 2 D F T (unvec (v_{2,2})) \\ ⋮ & ⋮ \\ 2 D F T (unvec (v_{1, N - M})) & 2 D F T (unvec (v_{2, N - M})) \end{matrix}]

(29)

故空间谱的计算可以通过对矩阵V₁和矩阵V₂的行矢量做2DFT获得，式（25）可改写为：

\begin{matrix} M (θ, ϕ) = Q^{H} Q \\ = [\begin{matrix} {(2 D F T (v_{1, j}))}^{H} \\ {(2 D F T (v_{2, j}))}^{H} \end{matrix}] [2 D F T (v_{1, j}) 2 D F T (v_{2, j})] \end{matrix}

(30)

改进的2D-FFT-MUSIC算法利用2DFT来计算空间谱，避免了2D-PES-MUSIC算法复杂的二维搜索过程，大大减少了计算量，提高了算法的运行速度和应用时的实效性。

3 仿真实验及分析

3.1 MST重构精度分析

根据理论分析可知，2DFT的系数随着阶数的增加具有超指数衰减的特征，即EADF矩阵的模值具有超指数衰减的特征，如图5所示。由于高阶的EADF矩阵模值很小，故可以将得到的EADF矩阵进行截断，从而减少运算量，分别在方位角维度和俯仰角维度截取Ma和Me阶数据作为采样矩阵，关于Ma和Me的确定可以选择不同的值。

图5EADF矩阵幅度与截断模数的关系

Fig.5Relation of EADF matrix amplitude and truncation mode

利用MST重构的导向矢量精度与截断模数Ma、Me的选取有关，现改变截断模数值，分析重构精度与模数的关系。定义归一化误差来描述MST重构精度：

n e = \frac{∥ \tilde{a} (θ, ϕ, γ, η) - \hat{a} (θ, ϕ, γ, η) ∥}{∥ \tilde{a} (θ, ϕ, γ, η) ∥}

(31)

式中，

\tilde{a} （ θ ， ϕ ， γ ， η ）

表示实际的导向矢量，

\hat{a} （ θ ， ϕ ， γ ， η ）

表示由MST生成的导向矢量。控制截断模数从2到180变化，取值间隔为1；方位角取θ∈[0，2π），俯仰角取φ∈[0，π），步进均为1°，求取

\hat{a} （ θ ， ϕ ， γ ， η ）

相对于

\tilde{a} （ θ ， ϕ ， γ ， η ）

的归一化误差，再求和取平均值，则重构精度与截断模数的关系如图6所示。

图6重构精度与截断模数的关系

Fig.6Relation of reconstruction precision and truncation mode

由图6可以看出，随着截断模数的增加，归一化误差整体呈下降趋势，但运算量随之增大；当截断模数增大到一定程度后，归一化误差可以忽略不计。综合重构精度和运算复杂度两个因素，选取截断模数Ma=Me=60进行仿真。

3.2 空域-极化域联合估计精度分析

3.2.1 单信源估计精度分析

首先考虑单个信源入射至非理想极化敏感阵列的情况，设置信源空域信息和极化域信息为（θ，φ，γ，η）=（150°，23°，46°，120°），快拍数为100，信噪比为20 dB，截断模数Ma=Me=60，八元共形阵的阵元均由单偶极子构成，阵列半径r为0.02 m，信号频率为3 GHz。设置各通道幅度误差在0~10 dBV范围内随机生成，相位误差在-10°~10°范围内随机生成，其中幅相误差依赖角度变化。利用传统MUSIC算法对空域和极化域信息进行联合估计，得到谱峰图如图7所示。

图7未校正算法二维谱图（单信源）

Fig.72D spectrum of uncorrected algorithm (single information source)

使用改进的2D-FFT-MUSIC校正算法，得到的单信源测向结果如图8所示。

图8改进的2D-FFT-MUSIC算法二维谱图（单信源）

Fig.82D spectrum of revised 2D-FFT-MUSIC algorithm (single information source)

通过图7和图8的对比可以看出，改进的2D-FFT-MUSIC算法可以在阵列误差存在的情况下得到尖锐谱峰，说明该校正算法可以对存在幅相误差和阵元位置误差的阵列进行精准测向。

接下来对2D-PES-MUSIC和改进的2D-FFT-MUSIC两种算法的性能进行分析。设置信噪比分别为0 dB、5 dB、10 dB、15 dB、20 dB、25 dB、30 dB，信源信息随机生成，在不同信噪比条件下做500次蒙特卡罗实验，分别得到两种算法的均方根误差随信噪比变化关系，如图9所示。假设每次估计结果与事先设置的角度相差在3°内视为估计成功，两种算法的测向成功概率随信噪比变化关系如图10所示。

根据图9和图10可以看出，与2D-PES-MUSIC算法相比，改进的2D-FFT-MUSIC算法可获得良好的估计性能和较高的测向成功率；随着信噪比的提升，估计结果的误差减小，测向成功的概率增大，估计性能改善。改进的2D-FFT-MUSIC算法在10 dB低信噪比条件下依旧能够获得较好的测向结果和较高的成功概率。

图9单信源均方根误差变化情况

Fig.9Change of root mean square error of single information source

图10单信源成功率变化情况

Fig.10Change of success rate of single information source

3.2.2 双信源估计精度分析

考虑多个信源入射至阵列的情况，设置两个来波方向分别为（θ₁，φ₁）=（150°，23°）、（θ₂，φ₂）=（220°，45°），极化信息分别为（γ₁，η₁）=（46°，20°）、（γ₂，η₂）=（65°，120°），其余条件不做改动，分别利用传统MUSIC算法和改进的2D-FFT-MUSIC算法对空域和极化域信息进行联合估计，得到谱峰图如图11所示。使用改进的2D-FFT-MUSIC校正算法，得到的多信源测向结果如图12所示。

通过图11和图12的对比可以看出，对于存在阵列误差的多信源测向问题，不进行误差校正直接测向会出现错误；改进的2D-FFT-MUSIC算法可以在阵列误差存在的情况下得到尖锐谱峰，说明该校正算法可以对存在幅相误差和阵元位置误差的阵列进行精准测向。

图11未校正算法二维谱图（双信源）

Fig.112D spectrum of uncorrected algorithm (dual information sources)

接下来考虑双信元入射时2D-PES-MUSIC和改进的2D-FFT-MUSIC两种算法的测向性能，两信源信息随机生成，其余条件不做改动，在不同信噪比条件下做500次蒙特卡罗实验，分别得到两种算法的均方根误差随信噪比的变化关系，如图13所示。假设每次估计结果与事先设置的角度相差在3°内视为估计成功，两种算法的测向成功率随信噪比变化关系如图14所示。

图12改进的2D-FFT-MUSIC算法二维谱图（双信源）

Fig.122D spectrum of revised 2D-FFT-MUSIC algorithm (dual information sources)

图13双信源均方根误差变化情况

Fig.13Change of root mean square error of dual information sources

根据图13和图14可以看出，对于多信源测向的情况，与2D-PES-MUSIC算法相比，改进的2D-FFT-MUSIC算法也可获得良好的估计性能和较高的测向成功率；随着信噪比的提升，估计结果的误差减小，测向成功的概率增大，估计性能改善。当信噪比提升至20 dB时，改进的2D-FFT-MUSIC算法能够获得较好的测向结果和较高的成功概率。

3.3 运算复杂度对比

现对基于MST的2D-PES-MUSIC算法和改进的2D-FFT-MUSIC算法的复杂度进行分析。由于两种算法在进行极化域信息估计时步骤相同，故不再考虑这一步的运算复杂度，只考虑两种算法在进行空域估计时的运算复杂度。

图14双信源成功率变化情况

Fig.14Change of success rate of dual information sources

利用快速子空间分解算法^[14]计算M×N维信号子空间E_s的复杂度为O（MN²），则基于MST的2D-PES-MUSIC算法的复杂度为O（MN²+L_aL_eMN），其中L_a和L_e为沿方位角和俯仰角搜索的数量，通常L_a和L_e远远大于M和N的值，故最终运算复杂度是由搜索点数L_a和L_e的值主导的。

计算二维FFT的运算复杂度为O（L²log₂L），其中L代表FFT点数，MATLAB中运行fft（n）函数点数取的是最接近数据长度的2的整数次方，则改进的2D-FFT-MUSIC算法的复杂度为O（MN²+2（N-M）MaMe+2（N-M）L²log₂L+QaQe）。

截取模数Ma=Me=60，则实际运算时的点数L=64，搜索点数取L_a=360、L_e=180，信源数M分别取1、2、3、4，不同阵元数量下运算复杂度对比如图15所示。

图15两种算法的运算复杂度对比

Fig.15Comparison of computational complexity of two algorithms

通过图15的对比可以看出，与2D-PES-MUSIC算法相比，改进的2D-FFT-MUSIC算法在运算复杂度方面有明显优势；信源数越多，改进算法的优势越明显，这是因为改进算法是直接对

V = E_{n}^{H} G

这一矩阵进行操作，在其他条件不变时，信源数的增加会导致矩阵V维度的减小，从而大幅度减少运算复杂度。

4 结论

本文针对存在多种阵列误差的极化敏感阵列，提出了一种基于MST的改进2D-FFT-MUSIC算法，该算法首先利用MST得到包含阵列误差的采样矩阵进行测向，并通过二维傅里叶变换求解二维空间谱，减少了谱峰搜索时的计算量。通过仿真讨论该算法的性能，并将其与已有的算法性能进行对比，证实了该算法的有效性。对于实际测向问题中只有部分数据已知的情况，该算法同样适用，且可以减少事先存储的数据量。该算法适用于存在任何阵列误差的精准测向，且适用于任意形式的阵列，对于存在阵列误差的测向方法研究具有普遍性意义。

图1极化敏感阵列空间摆放示意图

Fig.1Spatial layout diagram of polarization sensitive array

下载: 全尺寸图片

图2阵元摆放角示意图

Fig.2Diagram of array elements′ placement angle

下载: 全尺寸图片

图32DFT拓扑结构圆环体

Fig.32DFT topological torus

下载: 全尺寸图片

图4测量数据单位圆示意图

Fig.4Diagram of unit circle of measurement data

下载: 全尺寸图片

图5EADF矩阵幅度与截断模数的关系

Fig.5Relation of EADF matrix amplitude and truncation mode

下载: 全尺寸图片

图6重构精度与截断模数的关系

Fig.6Relation of reconstruction precision and truncation mode

下载: 全尺寸图片

图7未校正算法二维谱图（单信源）

Fig.72D spectrum of uncorrected algorithm (single information source)

下载: 全尺寸图片

图8改进的2D-FFT-MUSIC算法二维谱图（单信源）

Fig.82D spectrum of revised 2D-FFT-MUSIC algorithm (single information source)

下载: 全尺寸图片

图9单信源均方根误差变化情况

Fig.9Change of root mean square error of single information source

下载: 全尺寸图片

图10单信源成功率变化情况

Fig.10Change of success rate of single information source

下载: 全尺寸图片

图11未校正算法二维谱图（双信源）

Fig.112D spectrum of uncorrected algorithm (dual information sources)

下载: 全尺寸图片

图12改进的2D-FFT-MUSIC算法二维谱图（双信源）

Fig.122D spectrum of revised 2D-FFT-MUSIC algorithm (dual information sources)

下载: 全尺寸图片

图13双信源均方根误差变化情况

Fig.13Change of root mean square error of dual information sources

下载: 全尺寸图片

图14双信源成功率变化情况

Fig.14Change of success rate of dual information sources

下载: 全尺寸图片

图15两种算法的运算复杂度对比

Fig.15Comparison of computational complexity of two algorithms

下载: 全尺寸图片

图1极化敏感阵列空间摆放示意图

Fig.1Spatial layout diagram of polarization sensitive array

图2阵元摆放角示意图

Fig.2Diagram of array elements′ placement angle

图32DFT拓扑结构圆环体

Fig.32DFT topological torus

图4测量数据单位圆示意图

Fig.4Diagram of unit circle of measurement data

图5EADF矩阵幅度与截断模数的关系

Fig.5Relation of EADF matrix amplitude and truncation mode

图6重构精度与截断模数的关系

Fig.6Relation of reconstruction precision and truncation mode

图7未校正算法二维谱图（单信源）

Fig.72D spectrum of uncorrected algorithm (single information source)

图8改进的2D-FFT-MUSIC算法二维谱图（单信源）

Fig.82D spectrum of revised 2D-FFT-MUSIC algorithm (single information source)

图9单信源均方根误差变化情况

Fig.9Change of root mean square error of single information source

图10单信源成功率变化情况

Fig.10Change of success rate of single information source

图11未校正算法二维谱图（双信源）

Fig.112D spectrum of uncorrected algorithm (dual information sources)

图12改进的2D-FFT-MUSIC算法二维谱图（双信源）

Fig.122D spectrum of revised 2D-FFT-MUSIC algorithm (dual information sources)

图13双信源均方根误差变化情况

Fig.13Change of root mean square error of dual information sources

图14双信源成功率变化情况

Fig.14Change of success rate of dual information sources

图15两种算法的运算复杂度对比

Fig.15Comparison of computational complexity of two algorithms

图1极化敏感阵列空间摆放示意图

Fig.1Spatial layout diagram of polarization sensitive array

图2阵元摆放角示意图

Fig.2Diagram of array elements′ placement angle

图32DFT拓扑结构圆环体

Fig.32DFT topological torus

图4测量数据单位圆示意图

Fig.4Diagram of unit circle of measurement data

图5EADF矩阵幅度与截断模数的关系

Fig.5Relation of EADF matrix amplitude and truncation mode

图6重构精度与截断模数的关系

Fig.6Relation of reconstruction precision and truncation mode

图7未校正算法二维谱图（单信源）

Fig.72D spectrum of uncorrected algorithm (single information source)

图8改进的2D-FFT-MUSIC算法二维谱图（单信源）

Fig.82D spectrum of revised 2D-FFT-MUSIC algorithm (single information source)

图9单信源均方根误差变化情况

Fig.9Change of root mean square error of single information source

图10单信源成功率变化情况

Fig.10Change of success rate of single information source

图11未校正算法二维谱图（双信源）

Fig.112D spectrum of uncorrected algorithm (dual information sources)

图12改进的2D-FFT-MUSIC算法二维谱图（双信源）

Fig.122D spectrum of revised 2D-FFT-MUSIC algorithm (dual information sources)

图13双信源均方根误差变化情况

Fig.13Change of root mean square error of dual information sources

图14双信源成功率变化情况

Fig.14Change of success rate of dual information sources

图15两种算法的运算复杂度对比

Fig.15Comparison of computational complexity of two algorithms

LIU C, YIN S N. An efficient 2-D DOA estimation for a cylindrical conformal array with unknown mutual coupling[J]. International Journal of Antennas and Propagation,2018,2018:4015980.

汪争志, 杨志伟, 范志豪. 基于响应估计频域拟合的通道均衡方法[J]. 系统工程与电子技术,2022,44(3):713-721. WANG Z Z, YANG Z W, FAN Z H. Channel equalization method based on response estimation and frequency domain fitting[J]. Systems Engineering and Electronics,2022,44(3):713-721.(in Chinese)

杨志伟, 张攀, 陈颖, 等. 导向矢量和协方差矩阵联合迭代估计的稳健波束形成算法[J]. 电子与信息学报,2018,40(12):2874-2880. YANG Z W, ZHANG P, CHEN Y,et al. Steering vector and covariance matrix joint iterative estimations for robust beamforming[J]. Journal of Electronics & Information Technology,2018,40(12):2874-2880.(in Chinese)

张东兴. 基于稀疏贝叶斯理论的波达方向估计及阵列误差校正[D]. 哈尔滨: 哈尔滨工业大学,2020. ZHANG D X. Sparse off-grid DOA estimation method with array error calibration[D]. Harbin: Harbin Institute of Technology,2020.(in Chinese)

林潇, 薛敬宏, 乔晓林. 基于降维的极化敏感阵列幅相误差自校正算法[J]. 现代雷达,2020,42(10):51-56. LIN X, XUE J H, QIAO X L. Reduced-dimension self-correction algorithm for amplitude-phase error in polarization-sensitive array[J]. Modern Radar,2020,42(10):51-56.(in Chinese)

窦慧晶, 郭彩环, 张少飞, 等. 极化敏感阵列取向误差校正[J]. 吉林大学学报(工学版),2016,46(6):2087-2093. DOU H J, GUO C H, ZHANG S F,et al. Study on misorientation calibration of polarization sensitive arrays[J]. Journal of Jilin University(Engineering and Technology Edition),2016,46(6):2087-2093.(in Chinese)

洪贝. 阵列DOA估计及误差校准[D]. 成都: 电子科技大学,2020. HONG B. Array DOA estimation and error calibration[D]. Chengdu: University of Electronic Science and Technology of China,2020.(in Chinese)

胡显智, 王宁, 戴旭初. 基于流形分离技术的DOA和极化参数联合估计方法[J]. 测试技术学报,2018,32(5):435-442. HU X Z, WANG N, DAI X C. Joint estimation algorithm of DOA and polarization based on the manifold separation technique[J]. Journal of Test and Measurement Technology,2018,32(5):435-442.(in Chinese)

BELLONI F, RICHTER A, KOIVUNEN V. DoA estimation via manifold separation for arbitrary array structures[J]. IEEE Transactions on Signal Processing,2007,55(10):4800-4810.

COSTA M, KOIVUNEN V, RICHTER A. Azimuth,elevation,and polarization estimation for arbitrary polarimetric array configurations[C]//Proceedings of the IEEE/SP 15th Workshop on Statistical Signal Processing,2009.

杜毓杰. 非理想阵列条件下DOA估计相关问题研究[D]. 哈尔滨: 哈尔滨工业大学,2020. DU Y J. Research on problems regarding DOA estimation under non-ideal array conditions[D]. Harbin: Harbin Institute of Technology,2020.(in Chinese)

王永良, 陈辉, 彭应宁, 等. 空间谱估计理论与算法[M]. 北京: 清华大学出版社,2004. WANG Y L, CHEN H, PENG Y N,et al. Theory and algorithm of spatial spectrum estimation[M]. Beijing: Tsinghua University Press,2004.(in Chinese)

LANDMANN M, RICHTER A, THOMA R S. DoA resolution limits in MIMO channel sounding[C]//Proceedings of IEEE Antennas and Propagation Society Symposium,2004.

ZHUANG J, DUAN C H, WANG W,et al. Joint estimation of azimuth and elevation via manifold separation for arbitrary array structures[J]. IEEE Transactions on Vehicular Technology,2018,60(7):5585-5596.

首页

期刊介绍

投稿指南

编委会

出版声明

开放获取声明

联系我们

期刊订阅

Rss

English

1 模型的建立

2 基于MST的DOA估计方法

3 仿真实验及分析

4 结论