阵列处理技术在移动通信、雷达、声呐等领域有着广泛的应用^[1-3]。实际水声环境中存在大量的相干信号，例如海底、海面反射产生的多途信号，主动探测中多个目标的回波信号，以及声干扰器、声诱饵等人为干扰信号。相干信号子空间与噪声子空间的不完全正交会导致MUSIC、ESPRIT等经典的信号子空间类算法性能严重下降^[4-5]，因此相干源的波达方向(direction of arrival, DOA)估计逐渐成为水声阵列信号处理领域的热点问题。与均匀线阵(uniform linear array, ULA)相比，均匀圆阵(uniform circular array, UCA)具有以下优点：UCA可以提供360°无模糊的方位信息，在各个方位具有近似相同的估计精度和分辨力；UCA作为平面阵列能够进行二维DOA估计^[6]。然而，UCA的阵列流形不具有ULA的Vandermonde结构^[7]，因此空间平滑算法^[8-9]、矩阵重构算法^[10]等众多ULA的解相干算法无法直接应用到UCA上。为此，Griffiths提出了基于相位模式激励的UCA方向图合成技术^[11]，可以通过模式空间变换将圆阵的阵列流形变换成具有Vandermonde结构的矩阵，给解决UCA的相干源DOA估计问题提供了新思路。在此基础上，文献[12-13]将空间平滑技术应用于UCA，提出了模式空间平滑(mode spatial smoothing, MODESS)和模式空间前后向空间平滑(mode forward-backward spatial smoothing, MODEFBSS)；但是此类算法存在计算量大、平滑次数不易确定的问题。之后，ULA的矩阵重构算法被应用到UCA上，一种无需平滑的模式空间Toeplitz矩阵重构算法(mode Toeplitz, MODETOEP)^[14-15]能够实现信号协方差矩阵的对角化，具有较低的信噪比(signal-to-noise ratio, SNR)门限；但是MODETOEP没有充分利用UCA的对称性，具有改进的空间。受上述文献启发，本文提出了一种基于UCA中心对称性的相干源DOA估计方法。

1 UCA阵列信号模型

对于一个半径为R的UCA，M个各向同性的阵元均匀分布在圆周上，N个远场窄带实信号入射到阵列上，快拍数为L。如图 1所示，用球坐标系描述信号的来波方向，以圆心O为参考点，方位角θ表示入射信号在x－y平面的投影与x轴沿顺时针方向的夹角，俯仰角φ表示入射信号与x－y平面之间的夹角，因此有$\theta \in [0, 2\pi ], \varphi \in [0, \pi /2]$信号的方向向量$\mathit{\boldsymbol{r}} = [\cos \varphi \cos \theta - \cos \varphi \sin \theta \;\;\;{\kern 1pt} \sin \varphi ]$m号阵元与x轴的夹角为θ_m=2πm/M，阵元坐标可以表示p_m=(Rcosθ_m, －Rsin θ_m, 0)，m号阵元的接收信号s_i(t), (i=1, 2, …, N)相对于参考点的时延和相移分别为：

$ {\tau _{mi}} = \frac{R}{c}\cos \left( {{\theta _i} - {\theta _m}} \right)\cos {\varphi _i} $

(1)

$ {\phi _{mi}} = \frac{{2{\rm{ \mathsf{ π} }}R}}{\lambda }\cos \left( {{\theta _i} - {\theta _m}} \right)\cos {\varphi _i} $

(2)

由此可得信号的导向矢量：

$ \mathit{\boldsymbol{a}}\left( {{\theta _i},{\varphi _i}} \right) = {\left[ {\exp \left( { - {\rm{j}}{\phi _{0i}}} \right)\quad \cdots \quad \exp \left( { - {\rm{j}}{\phi _{\left( {M - 1} \right)i}}} \right)} \right]^{\rm{T}}} $

(3)

本文只讨论所有信源都在UCA平面上的情况，即俯仰角$\varphi = 0, \mathit{\boldsymbol{a}}\left( {{\theta _i}, {\varphi _i}} \right) = \mathit{\boldsymbol{a}}\left( {{\theta _i}} \right)$。假设各阵元接收噪声n_m(t)为与信号独立的加性高斯白噪声，均值为0，方差为σ_n²。则t时刻m号阵元上的接收数据为：

$ {x_m}\left( t \right) = \sum\limits_{i = 1}^N {\exp } \left( { - {\rm{j}}\beta \cos \left( {{\theta _i} - {\theta _m}} \right)} \right){s_i}\left( t \right) + {n_m}\left( t \right) $

(4)

式中β=2πR/λ。则阵列接收信号可以表示成：

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{X}}\left( t \right) = \mathit{\boldsymbol{A}}\left( \theta \right)\mathit{\boldsymbol{S}}\left( t \right) + \mathit{\boldsymbol{N}}\left( t \right) = }\\ {{{\left[ {\begin{array}{*{20}{l}} {{x_0}(t)}&{{x_1}(t)}& \cdots &{{x_{M - 1}}(t)} \end{array}} \right]}^{\rm{T}}}} \end{array} $

(5)

式中：$\mathit{\boldsymbol{N}}(t) = \left[ {\begin{array}{*{20}{l}} {{n_0}(t)}&{{n_1}(t)}& \cdots &{{{\left. {{n_{M - 1}}(t)} \right]}^{\rm{T}}}} \end{array}} \right.$为噪声矢量；$\mathit{\boldsymbol{S}}(t) = {\left[ {\begin{array}{*{20}{l}} {{s_1}(t)}&{{s_2}(t)}& \cdots &{{s_N}(t)} \end{array}} \right]^{\rm{T}}}$为入射信号矢量；$\mathit{\boldsymbol{A}}(\mathit{{\theta }}) = \left[ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{a}}\left( {{\theta _1}} \right)}&{\mathit{\boldsymbol{a}}\left( {{\mathit{{\theta }}_\mathit{\boldsymbol{2}}}} \right)}& \cdots &{\mathit{\boldsymbol{a}}\left( {{\theta _N}} \right)} \end{array}} \right]$为UCA对信号的M×N阶阵列流形矩阵。

2 基于UCA中心对称性的相干源DOA估计算法

本文提出的解相干算法，首先采用模式空间变换算法将UCA转变为虚拟均匀线阵(virtual uniform linear array, VULA)，使其阵列流形具有易于处理的Vandermonde结构，然后利用UCA模式空间阵列流形的中心对称性对数据进行共轭平均处理，并计算相关函数；最后重构一个Hermitian Toeplitz矩阵以实现信号协方差矩阵的对角化，达到解相干的目的。

2.1 模式空间变换

图 2表示VULA的阵列模型，d表示VULA的阵元间距，K为UCA可激励的模式数。VULA的输出数据Y(t)由原阵列接收数据乘以转换矩阵 T 得到：

$ \mathit{\boldsymbol{Y}}\left( t \right) = \mathit{\boldsymbol{TX}}\left( t \right) = {\mathit{\boldsymbol{A}}_v}\left( \theta \right)\mathit{\boldsymbol{S}}\left( t \right) + \mathit{\boldsymbol{TN}}\left( t \right) $

(6)

$ {y_k}\left( t \right) = \sum\limits_{i = 1}^N {{s_i}} \left( t \right)\exp \left( {{\rm{j}}k{\theta _i}} \right) + {{\tilde n}_k}\left( t \right) $

(7)

式中k=－K, －K+1, …, K。

UCA可激励的最大模式数${K_{\max }} = \left\lfloor \beta \right\rfloor $, 且满足$2{K_{\max }} < M, \left\lfloor \cdot \right\rfloor $表示向下取整。当K取K_max时，转换后的VULA具有最大的阵列孔径。根据文献[16]，转换矩阵 T 的计算公式为：

$ \mathit{\boldsymbol{T}} = {\mathit{\boldsymbol{J}}^{ - 1}}\mathit{\boldsymbol{F}}/M $

(8)

式中：F 是空间离散傅里叶变换的子矩阵；J 是由k阶的第1类Bessel函数J_k(·)构成的对角阵：

$ \mathit{\boldsymbol{F}} = {\left[ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{w}}_{ - K}}}&{{\mathit{\boldsymbol{w}}_{ - K + 1}}}& \cdots &{{\mathit{\boldsymbol{w}}_K}} \end{array}} \right]^{\rm{H}}} $

(9)

$ {\mathit{\boldsymbol{w}}_K} = {\left[ {\begin{array}{*{20}{c}} 1&{\exp \left( {\frac{{ - {\rm{j}}2{\rm{ \mathsf{ π} }}k}}{M}} \right)}& \cdots &{\exp \left( {\frac{{ - {\rm{j}}2{\rm{ \mathsf{ π} }}k\left( {M - 1} \right)}}{M}} \right)} \end{array}} \right]^{\rm{T}}} $

(10)

$ \begin{array}{*{20}{c}} {J = \left[ {\begin{array}{*{20}{c}} {{j^{ - K}}{J_{ - K}}\left( { - \beta } \right)}&{}&{}\\ {}& \ddots &{}\\ {}&{}&{{j^K}{{\rm{J}}_K}\left( { - \beta } \right)} \end{array}} \right] = }\\ {{\rm diag}\left( {{j^k}{J_k}\left( { - \beta } \right)} \right)} \end{array} $

(11)

转换后的虚拟阵列流形为：

$ {\mathit{\boldsymbol{A}}_v}\left( \theta \right) = \left[ {\begin{array}{*{20}{c}} {\exp \left( { - {\rm{j}}K{\theta _1}} \right)}& \cdots &{\exp \left( { - {\rm{j}}K{\theta _N}} \right)}\\ \vdots &{}& \vdots \\ {\exp \left( {{\rm{j}}K{\theta _1}} \right)}& \cdots &{\exp \left( {{\rm{j}}K{\theta _N}} \right)} \end{array}} \right] $

(12)

可以看出，(2K+1)×(2K+1)阶矩阵A_v(θ)具有与ULA阵列流形相同的Vandermonde结构。模式空间数据相关矩阵为：

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{R}}_Y} = {\rm{E}}\left[ {\mathit{\boldsymbol{Y}}\left( t \right){\mathit{\boldsymbol{Y}}^{\rm{H}}}\left( t \right)} \right] = {\mathit{\boldsymbol{A}}_v}\left( \theta \right){\rm{E}}\left[ {\mathit{\boldsymbol{S}}\left( t \right){\mathit{\boldsymbol{S}}^{\rm{H}}}\left( t \right)} \right]\mathit{\boldsymbol{A}}_v^{\rm{H}}\left( \theta \right) + }\\ {{\rm{E}}\left[ {\mathit{\boldsymbol{TN}}\left( t \right){\mathit{\boldsymbol{N}}^{\rm{H}}}\left( t \right){\mathit{\boldsymbol{T}}^{\rm{H}}}} \right] = {\mathit{\boldsymbol{A}}_v}\left( \theta \right){\mathit{\boldsymbol{R}}_S}\mathit{\boldsymbol{A}}_v^{\rm{H}}\left( \theta \right) + {{\mathit{\boldsymbol{\tilde R}}}_n}} \end{array} $

(13)

式中：

$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\tilde R}}}_n} = \mathit{\boldsymbol{\sigma }}_n^2\mathit{\boldsymbol{T}}{\mathit{\boldsymbol{T}}^{\rm{H}}} = }\\ {\frac{{\sigma _n^2}}{M}{\rm diag}\left( {\frac{1}{{J_K^2\left( \beta \right)}},\frac{1}{{J_{K - 1}^2\left( \beta \right)}}, \cdots ,\frac{1}{{J_K^2\left( \beta \right)}}} \right)} \end{array} $

(14)

因此，经过转换矩阵T实现了从原阵元空间中的UCA到模式空间中的ULA的转换，为在Y(t)的基础上进行解相干提供了可能。

2.2 利用中心对称性的解相干矩阵重构

经典的MODEFBSS算法受信噪比影响较大，在低信噪比下性能远不如MODETOEP；而MODETOEP算法中重构的矩阵为Toeplitz矩阵，不具有理想独立信源协方差矩阵的复共轭对称性质。本文算法利用UCA的中心对称性，对模式空间数据进行共轭平均处理，并构造与理想独立信源协方差矩阵具有相同性质的Hermitian Toeplitz矩阵，因此可以有效提高算法在低信噪比下的解相干能力和DOA估计精度。

用[R_Y]_{p, q}表示R_Y中第p行第q列的元素，σ_si²表示第i个信源的功率，σ_yn²为模式空间的噪声功率。当入射信号为理想情况下的独立信号时，有式(15)成立，其中r_Y(·)表示不同阵元接收数据间的相关函数。

$ \begin{array}{*{20}{c}} {{{\left[ {{\mathit{\boldsymbol{R}}_Y}} \right]}_{p,q}} = {\rm{E}}\left[ {{y_p}\left( t \right)y_q^{\rm{H}}\left( t \right)} \right] = {r_Y}\left( {p - q} \right) = }\\ {\sum\limits_{i = 1}^N {\sigma _{si}^2} \exp \left( {{\rm{j}}(p - q){\theta _i}} \right) + \sigma _n^2{\delta _{p,q}}/M = }\\ {r_Y^*(q - p) = \left[ {{\mathit{\boldsymbol{R}}_Y}} \right]_{q,p}^*} \end{array} $

(15)

$ {\delta _{p,q}} = \left\{ {\begin{array}{*{20}{l}} {1,}&{p = q}\\ {0,}&{p \ne q} \end{array}} \right. $

(16)

用r_Y(·)来表示模式空间数据的相关矩阵 R_Y：

$ {\mathit{\boldsymbol{R}}_Y} = \left[ {\begin{array}{*{20}{c}} {{r_Y}(0)}&{{r_Y}( - 1)}& \cdots &{{r_Y}( - 2K)}\\ {{r_Y}(1)}&{{r_Y}(0)}& \cdots &{{r_Y}(1 - 2K)}\\ \vdots & \vdots &{}& \vdots \\ {{r_Y}(2K)}&{{r_Y}(2K - 1)}& \cdots &{{r_Y}(0)} \end{array}} \right] $

(17)

可见，独立信号源的协方差矩阵任意一个对角线上的元素取相同值，且矩阵中的元素满足复共轭对称关系，即为Hermitian Toeplitz矩阵。在实际阵列信号处理中，受到系统误差以及信源相干性的影响，R_Y一般只是对角占优的矩阵，而不是Hermitian Toeplitz矩阵。

由UCA的中心对称性可知：

$ {\mathit{\boldsymbol{A}}_v}(\theta ) = {\mathit{\boldsymbol{B}}_{2K + 1}}\mathit{\boldsymbol{A}}_v^*(\theta ) $

(18)

式中：B为副对角线上元素全为1，其他元素全为0的置换矩阵。UCA在模式空间的阵列流形具有对称性，因此可以将中心阵元(0号阵元)一侧的阵列数据进行共轭重排后与另一侧阵元的数据进行平均。中心阵元右侧的接收数据${\mathit{\boldsymbol{Y}}_r}\left( t \right) = {\left[ {\begin{array}{*{20}{l}} {{y_0}(t)}& \cdots &{{y_K}(t)} \end{array}} \right]^{\rm{T}}}$, 即Y_r(t)为k=0, 1，…, K的模式空间接收数据组合。对于中心阵元左侧的模式空间阵列接收数据$\mathit{\boldsymbol{Y}}(t) = {\left[ {\begin{array}{*{20}{l}} {{y_{ - K}}(t)}& \cdots &{{y_0}(t)} \end{array}} \right]^ \top }$，按照式(19)进行倒序排列：

$ {\mathit{\boldsymbol{Y}}_l}(t) = {\mathit{\boldsymbol{B}}_{K + 1}}{\mathit{\boldsymbol{Y}}^*}(t) = {\left[ {\begin{array}{*{20}{l}} {y_0^*(t)}& \cdots &{y_{ - K}^*(t)} \end{array}} \right]^{\rm{T}}} $

(19)

对处理后的2部分接收数据进行平均，由于实际中的发射信号和噪声均为实信号，单个阵元上的接收数据可以表示为：

$ \begin{array}{*{20}{c}} {{z_k}(t) = \left( {{y_k}(t) + y_{ - k}^*(t)} \right)/2 = }\\ {\left( {\sum\limits_{i = 1}^N {{s_i}} (t)\exp \left( {{\rm{j}}k{\theta _i}} \right) + {{\tilde n}_k}(t) + } \right.}\\ {\left. {\sum\limits_{i = 1}^N {s_i^*} (t)\exp \left( {{\rm{j}}k{\theta _i}} \right) + \tilde n_{ - k}^*(t)} \right)/2 = }\\ {\sum\limits_{i = 1}^N {{s_i}} (t)\exp \left( {{\rm{j}}k{\theta _i}} \right) + \left( {{{\tilde n}_k}(t) + {{\tilde n}_{ - k}}(t)} \right)/2} \end{array} $

(20)

比较式(7)和式(20)，可以看出信号部分没有变化，随机噪声被平均。当k≠0时，因为不同阵元上接收的噪声是不相关的，所以噪声功率减半。将处理后的数据表示为式(21)中的矢量形式，式(22)为新的阵列流形矩阵，可见共轭平均处理损失了一半的阵列孔径，阵元个数减少到K+1个。

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{Z}}(t) = \left( {{\mathit{\boldsymbol{Y}}_r}(t) + {\mathit{\boldsymbol{Y}}_l}(t)} \right)/2 = }\\ {{{\left[ {\begin{array}{*{20}{l}} {{z_0}(t)}&{{z_1}(t)}& \cdots &{{z_K}(t)} \end{array}} \right]}^{\rm{T}}} = }\\ {{{\mathit{\boldsymbol{\bar A}}}_v}(\theta )\mathit{\boldsymbol{S}}(t) + {\mathit{\boldsymbol{N}}_z}(t)} \end{array} $

(21)

$ {{\mathit{\boldsymbol{\bar A}}}_v}(\theta ) = \left[ {\begin{array}{*{20}{c}} 1& \cdots &1\\ \vdots &{}& \vdots \\ {\exp \left( {{\rm{j}}K{\theta _1}} \right)}& \cdots &{\exp \left( {{\rm{j}}K{\theta _N}} \right)} \end{array}} \right] $

(22)

定义以下相关函数表示各阵元接收数据与中心阵元的相关性：

$ \begin{array}{*{20}{c}} {r(k) = {\rm{E}}\left[ {{\mathit{\boldsymbol{z}}_0}(t)\mathit{\boldsymbol{z}}_k^{\rm{H}}(t)} \right] = \sum\limits_{i = 1}^N {{d_i}} \exp \left( { - {\rm{j}}k{\theta _i}} \right) + }\\ {\sigma _n^2{\delta _{0,k}}/M} \end{array} $

(23)

式中：d_i表示第i个信号与所有信号的互相关之和：

$ {d_i} = \sum\limits_{l = 1}^N {\rm{E}} \left[ {{\mathit{\boldsymbol{s}}_i}(t)\mathit{\boldsymbol{s}}_l^{\rm{H}}(t)} \right]\;\;\;\;i = 1,2, \cdots ,N $

(24)

比较式(15)和式(23)可以看出，相关函数r(k)具有与r_Y(k)相似的结构，同样包含了所有信号的入射角信息。

利用式(23)可以得到K+1个相关值，构造如下Hermitian Toeplitz矩阵，式中I_K+1表示单位阵：

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{R}}_{{\rm{HT}}}} = \left[ {\begin{array}{*{20}{c}} {r(0)}&{r(1)}& \cdots &{r(K)}\\ {{r^*}(1)}&{r(0)}& \cdots &{r(K - 1)}\\ \vdots & \vdots &{}& \vdots \\ {{r^*}(K)}&{{r^*}(K - 1)}& \cdots &{r(0)} \end{array}} \right] = }\\ {\overline {{\mathit{\boldsymbol{A}}_v}} \mathit{\boldsymbol{D}}\overline {\mathit{\boldsymbol{A}}_v^{\rm{H}}} + \frac{{\sigma _n^2}}{M}{\mathit{\boldsymbol{I}}_{K + 1}}} \end{array} $

(25)

式中

$ \mathit{\boldsymbol{D}} = diag\left( {{d_1},{d_2}, \cdots ,{d_N}} \right) $

(26)

R_HT与理论上独立信源协方差矩阵具有相同的复共轭对称性质，可以用其代替数据协方差矩阵进行特征分解。R_HT由对角阵D和Vandermonde矩阵相乘得到，因此秩数rank(R_HT)=min(N, K+1)。只要K+1>N，且信号方位角各不相同，重构的协方差矩阵的秩数就等于信号源数，可对信号子空间与噪声子空间做出正确的估计，解决了高分辨算法在信号相干时性能严重下降的问题。

对R_HT进行广义特征值分解，得到信号子空间U_S和噪声子空间U_N后，按照式(27)所示的MUSIC谱峰搜索公式便可以得到相干源DOA估计结果，搜索矢量${\mathit{\boldsymbol{\bar a}}_v}(\theta ) = \left[ {\begin{array}{*{20}{l}} 1&{\exp ({\rm{j}}\theta )}& \cdots &{\exp ({\rm{j}}K\theta )} \end{array}} \right]$。

$ {P_{{\rm{MUSIC}}}}(\theta ) = 1/\left( {\mathit{\boldsymbol{\bar a}}_v^{\rm{H}}(\theta ){\mathit{\boldsymbol{U}}_N}\mathit{\boldsymbol{U}}_N^{\rm{H}}{{\mathit{\boldsymbol{\bar a}}}_v}(\theta )} \right) $

(27)

2.3 算法流程及计算量分析

基于以上分析，本文算法的具体流程可以归纳为4步：1)利用式(8)~(11)求模式空间转换矩阵T，由式(6)将UCA接收数据X (t)变换为模式空间数据Y(t)；2)利用式(22)求共轭平均后的数据Z(t)；3)根据式(23)计算Z(t)的相关函数r(k)，利用式(25)重构Hermitian Toeplitz矩阵R_HT；4)对R_HT进行特征分解，由式(27)得到DOA估计结果。

UCA解相干算法的主要计算量集中在模式空间变换、解相干矩阵构造和对矩阵进行特征分解上。根据上文中的算法流程，由于模式激励法的变换矩阵可以离线计算，本文算法在模式变换时的计算复杂度为O((2K+1)²L)；忽略求和与平均运算，求相关函数时的计算复杂度为O((K+1)L)；特征分解的计算复杂度为O((K+1)³)。所以本文算法的计算复杂度为O((2K+1)²L+(K+1)L+(K+1)³)。当MODEFBSS算法的平滑次数为(K+1)时，解相干处理矩阵的阶数与本文算法相同，计算复杂度为O((2K+1)²L+(K+1)³L+(K+1)³)，MODETOEP的计算量为O((2K+1)²L+(2K+1)L+ (K+1)³)。因此，本文算法的计算量略小于以上2种算法。

3 仿真验证

为验证本文算法的有效性，并与MODETOEP和MODEFBSS进行性能比较，本文设计了仿真实验。实验中UCA的阵元间隔为半波长，阵元数为21，2个远场窄带信号的相关系数为1，中心频率为30 kHz，默认的快拍数为600，SNR为0 dB。

3.1 分辨率

图 3为本文算法、MODETOEP、MODEFBSS以及MUSIC算法在3种情况下的空间谱估计结果。由图 3(a)可知，信号源相干时MUSIC算法已经失效，无法得到可靠的DOA估计结果；本文算法可以成功地估计相干源的方向，仿真结果验证了算法的有效性；其他2种算法也可以成功分辨相干源，但本文算法得到的谱峰最尖锐明显，对噪声的抑制能力最强，分辨率最高。在图 3(b)中，随着SNR下降到-10 dB，3种算法的分辨率均有所降低，MODEFBSS性能下降最为严重。在图 3(c)中，信号入射角度间隔从25°减小到15°，MODEFBSS无法分辨2个相干信源，本文算法和MODETOEP仍然有效。综上所述，在相同的仿真条件下，本文算法的分辨率高于MODETOEP和MODEFBSS。

3.2 解相干能力

算法的解相干能力可以用分辨概率来衡量，分辨概率表示算法成功分辨出2个相干信号的概率，当所有DOA估计结果的偏差都在允许范围内，认为是一次成功的解相干。一个信号从-10°方向入射，改变另一个信号的入射角，偏差门限取角度间隔的一半，进行300次蒙特卡洛独立实验，比较3种算法的分辨概率，结果如图 4(a)所示。随着角度间隔的增大，MODETOEP、MODEFBSS和本文算法分辨相干信源的成功概率均逐渐增加至1；当分辨概率大于0.9时，本文算法对应的角度间隔为12°，MODETOEP为14°，MODEFBSS为18°，因此本文算法对角度间隔的要求最低。

相干信号分别从-10°和15°入射，偏差门限取2°，改变SNR进行300次蒙特卡洛实验，结果如图 4(b)所示。随着SNR的增加，3种算法的分辨概率逐渐增加；当SNR为-9 dB时，本文算法的分辨概率为1，MODETOEP的分辨概率为0.98，MODEFBSS分辨概率为0.18；在-5 dB信噪比以下，本文算法的解相干能力显著优于MODEFBSS，略优于MODETOEP。

当SNR取0 dB，偏差门限为1°时，快拍数从10变化到560，分辨概率随快拍数变化情况如图 4(c)所示。此种仿真条件下，相同快拍数时MODEFBSS的分辨概率最小，本文算法的分辨概率大于另外2种算法。

3.3 DOA估计精度

本文利用均方根误差(root mean square error, RMSE)衡量算法的DOA估计精度，具体计算过程：

$ {\rm RMSE} = \sqrt {\frac{1}{{{N_{{\rm{sig}}}}{N_{{\rm{det}}}}}}\sum\limits_{i = 1}^{{N_{{\rm{sig}}}}} {\sum\limits_{j = 1}^{{N_{{\rm{det}}}}} {{{\left( {{{\bar \theta }_{ij}} - {\theta _i}} \right)}^2}} } } $

(28)

式中：N_sig表示实验中已知的信号个数；N_det表示成功分辨出2个相干信号的次数。在分辨概率都达到0.8以上时，比较3.2节中3种算法DOA估计的RMSE。

从图 5 DOA估计结果可以看出，随着角度间隔、SNR和快拍数的增加，3种算法的RMSE均逐渐减小，DOA估计精度逐渐提高；在相同仿真条件下时，本文方法DOA估计的RMSE比MODEFBSS和MODETOEP小，即估计精度更高。

3.4 算法稳健性

以上仿真结果都是在理想阵列模型下得到的，而在实际中阵列幅相误差是普遍存在的，这就对算法的稳健性提出了考验。幅相误差对阵列流形矩阵中元素的幅度和相位产生扰动，用u_k和ϑ_k分别表示阵元k的幅度误差和相位误差，仿真实验中依式(29)获得^[17]：

$ \left\{ {\begin{array}{*{20}{l}} {{u_k} = 1 + \sqrt {12} {a_k}{\sigma _u}}\\ {{\vartheta _k} = \sqrt {12} {b_k}{\sigma _\vartheta }} \end{array}} \right. $

(29)

式中：a_k和b_k是在[-0.5, 0.5]内均匀分布的随机数，σ_u和σ_ϑ分别表示u_m和ϑ_m的标准差。信号从-10°和15°入射，改变σ_u，保持σ_ϑ等于0；改变σ_ϑ，保持σ_u等于0。分别进行500次随机实验，结果如图 6所示。可见，随着幅度误差和相位误差的增大，3种算法的相干源分辨概率均逐渐减小；本文算法的稳健性不如MODEFBSS和MODETOEP，且受相位误差的影响较大。

4 结论

1) 与空间平滑类算法相比，本文算法无需多次平滑、无需计算接收数据协方差矩阵，计算量更小；

2) 与MODETOEP相比，本文算法利用了UCA的中心对称性，解相干能力更强；

3) 本文算法能够估计相干信源的DOA，其分辨率和DOA估计精度高于MODEFBSS和MODETOEP，但幅相误差存在时的稳健性较差。

本文算法对中心对称阵列的信号处理有借鉴意义，但是本文算法依赖于预处理后的中间阵元，因此如何降低对中间阵元的依赖、充分利用其它阵元间的互相关信息，值得今后进一步研究。