基于稀疏对称阵列的混合信源定位

近年来，信源定位成为阵列信号处理中的一个热点问题并受到广泛关注。大量针对远场条件下的波达方向(Direction of Arrival, DOA)估计算法被提出[1-6]。当信源位于菲涅尔区时，信号不再以平面波传播，而是以球面波的形式传播，且信源的信息由角度和距离两个参数共同决定，此时为近场源定位。为此，国内外学者提出了大量估计方法，如2-D 多重信号分类(Multiple Signal Classification, MUSIC)算法[7]和高阶旋转不变子空间(Estimation of Signal Parameters via Rotational Invariance Techniques，ESPRIT)算法[8]都可以有效估计。

然而在实际应用中近场信源和远场信源有可能同时存在，由于信源的不匹配，上述所提到的近场或远场信源DOA估计算法将不再适用。为了解决混合信源DOA估计问题，文献[9]通过构造高阶累积量矩阵，然后使用MUSIC算法估计近场信源参数，计算量太大；为了减少运算量，文献[10]提出一种基于二阶统计量的混合信源参数估计，然而该算法在进行近场DOA估计时只利用了部分协方差矩阵的信息，造成较大的阵列孔径损失，估计精度较差；为了减少计算复杂度和阵列孔径损失，文献[11]将对称阵均匀线阵分为2分子阵，利用子阵间的关系使用ESPRIT-like算法进行参数估计，但估计结果仍难以令人满意。

以上算法都是基于均匀线阵，存在着一定的阵列孔径损失，在阵元数一定的情况下，稀疏布阵可以增加阵列孔径，从而提高信源参数估计精度。文献[12]采用互素对称阵列，通过两个阵元数互为素数的子阵构造一个四阶累积量矩阵，然后利用MUSIC算法估计信源角度，有效扩展了阵列孔径。文献[13]使用了一种特殊几何阵列结构，利用传统的MUSIC算法和构造特殊点的四阶累积量联合估计信源的角度和距离。为了进一步扩展阵列孔径，本文提出一种稀疏对称阵列模型。首先通过对不同子阵的接收数据进行四阶累积量运算剔除近场信源的距离参数，构造多个仅与信源角度有关的四阶累积量向量，通过这些累积量向量构造一个Topelize矩阵，再利用MUSIC算法估计出信源角度，最后在估计出角度的基础上进行距离搜索，根据近场远场信源位于不同区域估计出近场信源的距离参数。该算法避免了二维搜索，且稀疏布阵扩展了阵列孔径，提高了参数估计精度。

1 信号模型

本文所采用的阵列模型如图1所示，阵元总数为2(N1+N2)+1,由3个子阵列组成。其中子阵1阵元数和阵元间距分别为2N1+1和d，子阵2和子阵3的阵元数和阵元间距分别为N2和(2N1+1)d,子阵列1与子阵列2和子阵列3之间的距离分别为(N1+1)d。图中各阵元坐标为pi，pi=-N2(2N1+1),-(N2-1)(2N1+1),…,-(2N1+1),-N1,…,0,…,N1,(2N1+1),…,(N2-1)(2N1+1),N2(2N1+1),pid为第i个阵元到中心阵元的距离。当阵元数相同时，均匀线阵的阵列孔径为(2N1+2N2+1)d，本文所采用阵列的阵列孔径为(2N2+4N1N2)d，可以看出，本文的阵列具有更大的阵列孔径。

假设空间存在K个独立窄带信号入射到阵列，以中心阵元为相位参考点，则第i个阵元接收数据为

式中，sk(t)为第k个信号源，T为快拍数，ni(t)为第i个阵元接收到的噪声， μk和φk分别为[14]

式中，λ为信号的波长，rk，θk分别为第k个近场信源的距离和角度参数。

将阵元接收数据写为矢量形式为

式中,x(t)=[x-N1-N2(t),…,x-N1(t),…,xN1(t),…,xN1+N2(t)]T为阵元接收数据，s(t)=[s1(t),s2(t),…,sK(t)]T为K个信源的信号矢量，n(t)=[n-N1-N2(t),…,nN1+N2(t)]T为阵元接收的噪声矢量，A为(2(N1+N2)+1)×K维阵列流型矢量，表示为

A=[a(θ1,r1),…,a(θk,rk),…a(θK,rK)]

式中，a(θk,rk)为(2(N1+N2)+1)×1维方向矢量，当信源位于近场时，a(θk,rk)表示为

a(θk,rk)=[ej[-(2N1N2+N2)μk+(-(2N1N2+N2))2φk],…,

ej[-(2N1+1)μk+(-(2N1+1))2φk],ej[-(N1)μk+(-N1)2φk],…,1,…,

ej[(2N1N2+N2)μk+(2N1N2+N2)2φk]]

当信源位于远场时，a(θk,rk)表示为

a(θk,rk)=[ej[-(2N1N2+N2)μk],…,

ej[-(2N1+1)μk],ej[-(N1)μk],…,

1,…,ej[(N1)μk],ej[(2N1+1)μk],…,

本文中，对信号作如下假设：

1) 信号为零均值、非高斯的窄带平稳随机过程，且具有非零峰度，信号之间不相关。

2) 阵元接收的噪声为零均值的高斯白噪声，并且与信号相互独立。

3) 为了确保信源参数估计的唯一性，阵列1阵元最小间距d≤λ/4。

2 基于稀疏对称阵列的信源定位

2.1 构造特殊四阶累积量矩阵

高阶累积量具有抑制高斯噪声，同时还可以扩展阵列孔径等优点，因此基于高阶累积量的空间谱估计得到广泛关注。本文所采用的四阶累积量公式为[15]

式中，

为信号的四阶累积量，*表示复共轭。为了剔除近场信源的距离参数φk，保留信源的角度参数μk，令

即i=-m,q=-n。因此，式(8)可写为

通过四阶累积量构造一个仅与信源角度有关的特殊矩阵，使其稀疏对称阵等效于阵元间距为λ/4的均匀线阵列。令m∈[-N1,…,N1]，n为0，首先，将子阵1阵元的接收数据与中心处的阵元的接收数据进行四阶累积运算得到一个(2N1+1)×1维四阶累积量向量c1，其第m个元素为

同理，将子阵1的阵元接收到的数据与子阵3的第一个阵元接收的数据进行四阶累积量的运算，构造出(2N1+1)×1维四阶累积量向量c2，其第m个元素分别为

将子阵1的阵元接收到的数据与子阵2的第一个阵元接收的数据进行四阶累积量的运算，构造出(2N1+1)×1维四阶累积量向量c3，其第m个元素分别为

令n∈[0,…,N1]，m∈[N1+2,…,N1+N2]，将子阵3阵元接收到的数据与子阵1的参考阵元右半部分阵元接收到的数据进行四阶累积量运算，构造出(N1+1)(N2-1)×1维向量c4，其第 ((m-N1-2)(N1+1)+n+1)个元素为

c4((m-N1-2)(N1+1)+n+1)=

同理，将子阵1阵元接收到的数据与子阵1的参考阵元左半部分阵元接收到的数据进行四阶累积量运算，构造出(N1+1)(N2-1)×1维向量c5，其第((m-N1-2)(N1+1)+n+1)个元素为

c5((m-N1-2)(N1+1)+n+1)=

将向量c5，c3，c1，c2，c4垒成一个(2(N1N2+2N1+N2)+1)×1维的长向量c,表示为

可以发现，向量c中元素的相位项m-n的值跟均匀线阵的相位等效，分别为：-(N1N2+2N1+N2),-(N1N2+2N1+N2-1),…,0,…,(N1N2+2N1+N2-1),(N1N2+ 2N1+N2)。通过四阶累积量运算使稀疏对称阵列的接收数据等效为一个均匀线阵，如图2所示。

矢量c1是由子阵1内各阵元的接收数据分别与子阵1中的中心参考阵元的接收数据做累积量运算的结果值组合而得;矢量c2是由子阵1内各阵元的接收数据分别与子阵3中第一个阵元的接收数据做累积量运算的结果值组合而得，矢量c3是由子阵1内各阵元的接收数据分别与子阵2中第一个阵元的接收数据做累积量运算的结果值组合而得，矢量c4是由子阵1内各阵元的接收数据分别与子阵3中除第一个阵元外的所有阵元接收数据做累积量运算的结果值组合而得，矢量c5是由子阵1内各阵元的接收数据分别与子阵2中除第一个阵元外的所有阵元接收数据做累积量运算的结果值组合而得。

在矢量的重构过程中，将矢量c1的第m个元素对应于重构矢量c的第m个元素，矢量c2的第m个元素对应于重构矢量c的第m+2N1+1个元素，矢量c3的第m个元素对应于重构矢量c的第-(m+2N1+1)个元素，矢量c4的第m个元素对应于重构矢量c的第(m-N1)(N1+1)+N1+n个元素，矢量c5的第m个元素对应于重构矢量c的第-((m-N1)(N1+1)+N1+n)个元素。

向量c构造一个(N1N2+2N1+N2+1)×(N1N2+2N1+N2+1)维的Topelize矩阵C，C的第m列可以表示为

C(:,m)=c(N1N2+2N1+N2+2-

m:2(N1N2+2N1+N2)+2-m)

因此，我们构造了一个仅包含信源角度信息的特殊矩阵C，类似于远场协方差矩阵，可以表示为

C4,S=diag[c4,s1,…,c4,sP]

A(θ)=[a(θ1),…,a(θp)]

a(θp)=[1,ej2μp,…,ej2(N1N2+2N1+N2)μp]T

针对接收端对回波数据的接收方式，与已有的将稀疏对称阵列划分为4个子阵来分别接收数据的方式不同[13]，本文中将稀疏对称阵列划分为3个子阵来分别接收数据。

2.2 信源DOA估计

由式(17)～式(20)可以看出，矩阵A(θ)类似于K个远场信号入射到阵元数为(N1N2+2N1+N2+1)的均匀线阵所产生的阵列流型矩阵， a(θk)为第k个信号的类远场阵列流型矢量。由于信源sk的四阶累积量非零，且A(θ)和 C4,S均为满秩，因此可以运用MUSIC算法来估计信源角度，对构造的四阶累积量矩阵C进行特征值分解为

Vs=R(N1N2+2N1+N2+1)×K为较大K个特征值构成的信号子空间，Vn=R(N1N2+2N1+N2+1)×((N1N2+2N1+N2+1)-K)为较小(N1N2+2N1+N2+1-K)个特征值构成的噪声子空间。由式(22)所示的MUSIC谱可估计出信源的方位角：

2.3 近场信源距离估计

阵列接收数据的协方差矩阵为

对R进行特征值分解

式中， Un∈R((2(N1+N2)+1)×(2(N1+N2)+1)-K)为较小(2(N1+N2)+1-K)个特征值构成的噪声子空间。将估计出的角度带入阵列流型向量，由距离MUSIC谱可以估计出近场信源的距离参数：

如果rk∈[0.62(D3/λ),2D2/λ]，其中D为阵列孔径，则与rk对应的信源位于菲涅耳区，属于近场信源。若rk>2D2/λ，则与rp对应的信源为远场信源。由此可以估计出近场信源参数，无需二维搜索，且近场信源的角度和距离参数自动配对。

由于采用稀疏对称阵列接收数据的协方差矩阵来估计近场信源的距离参数，直接用稀疏阵列估计信源参数是需要考虑空间模糊问题。信源θk方向上产生距离模糊的条件为

式中，pk为第k个信源的位置，l为整数，|l|≥1。由式(23)可得

要使式(26)右边取得最小值，令cos2θk=1， rk=0.62(D3/λ)1/2, r′k=∞，化简得

假设阵列孔径D=λ,最小阵元间距d=λ/4代入式(25)可得pk≥4.45。因此pk≥5或pk≤-5时阵列流型向量的第k个元素会产生模糊，然而当D=λ，d=λ/4时，-4≤pk≤4，所以整个阵列流型向量都不会产生模糊现象，估计出的每一个近场信源的距离参数都是唯一的。

2.4 计算复杂度分析

对本文、文献[11]、文献[12]及文献[13]四种算法的计算复杂度进行对比。由于运算量主要来自于构造累积量矩阵、特征值分解以及角度和距离的谱峰搜索这几部分，因此我们主要比较这几个方面的运算量。假设角度搜索的步长为Δθ，距离搜索的步长为Δr，快拍数为T。文献[12]需要构造一个(2N2+2N-1)×1的四阶累积向量以及一个(4N-1)×(4N-1)的协方差矩阵，对它们进行特征值分解，然后用MUSIC进行角度和距离的搜索，因此文献[12]的总计算量为O{9(2N2+2N-1)T+4/3(N2+N)3+Δθ(N2+N)2+(4N-1)2T+4/3(4N-1)3+KΔr(4N-1)2}。文献[11]需要构造一个(2N+1)×(2N+1) 的协方差矩阵，并对其进行特征值分解，最后使用求根算法进行DOA估计，文献[11]的计算量为

其中

为近场信源个数。文献[13]需要构造一个 (2(N2+1)(N1+1)+1)×1四阶累积向量以及(2N1+2N2+1)×(2N1+2N2+1)的协方差矩阵，分别对它们进行特征值分解，最后进行谱峰搜索，文献[13]算法的运算量为O{9(2N1N2+2N1+2N2+3)T+4/3(N1N2+N1+N2+1)3+ Δθ(N1N2+N1+N2+1)2+(2N1+2N2+1)2T+(2N1+2N2+1)2T+4/3(2N1+2N2+1)3+KΔr(2N1+2N2+1)2}。本文算法需要构造一个(2N1N2+2N1+N2+1)×1的四阶累积向量以及(2N1+2N2+1)×(2N1+2N2+1)的协方差矩阵，分别对它们进行特征值分解，最后进行谱峰搜索，因此本文算法的运算量为O{9(2N1N2+4N1+2N2+1)T+4/3(N1N2+2N1+N2+1)3+ Δθ(N1N2+2N1+N2+1)2+(2N1+2N2+1)2T +(2N1+2N2+1)2T+4/3(2N1+2N2+1)3+KΔr(2N1+2N2+1)2}。在实际的应用中，需要较大的快拍数，从以上分析可以看出当快拍数较大时，本文算法的计算复杂度要略大于文献[12]和文献[13]的计算复杂度，由于文献[11]使用二阶统计量来估计信源参数，其计算复杂度最低，但是相对本文算法它的阵列孔径更小，其估计精度更低。

3 计算机仿真

为了验证所提算法的优良性能，将本文所提算法与文献[11]中基于二阶统计量的ESPRIT算法、文献[12]中基于互素对称阵的近场源定位和文献[13]中基于特殊的几何阵列结构的远近场混合定位进行对比。假设阵元总数7(N1=1，N2=2)，最小阵元间距为λ/4，采用角度和距离的均方根误差(RMSE)作为衡量标准，分别定义为

式中：H为蒙特卡洛次数；

分别为角度和距离第h次实验的估计值；θk，rk分别为角度和距离第k个信源参数的实际值。

实验1：信号源为近场远场混合信源。入射信号为2个4QAM窄带信号，近场信源与远场信源位置分别为(θ1=-13°,r1=3.5λ)和(θ2=21°, r2=∞)。

在本实验中，验证本文算法DOA估计精度随信噪比和快拍数变化的情况。第一，信噪比从-5 dB到10 dB不断变化，步长为2 dB，快拍数为700，蒙特卡洛次数为500；第二，信噪比为7 dB，快拍数从100到1 100不断变化，步长为200。图3为信源角度均方根误差随信噪比的变化情况，图4为近场信源距离的均方根误差随信噪比的变化。图5和图6分别是方位角和近场信源距离随快拍数的变化曲线图。从图3可以看出，信源角度的估计精度随信噪比的增加而提高。由于本文算法采用稀疏对称阵列，其阵列孔径比均匀线阵[11]、互质对称阵列[12]以及特殊的几何阵列结构[13]大，因此本文算法角度估计精度要高于对比的算法。从图5和图6可以看出，4种算法在快拍数较小时参数估计精度较差，均随着快拍数的增加而提高。从图4和图6可以看出，本文算法中的近场信源距离估计精度在随着信噪比和快拍数的变化中都高于对比算法。本文算法和对比算法都是采用MUSIC算法来估计近场信源的距离，由于距离估计是基于角度估计，所以所提算法也具有更高的距离估计精度。

实验2：考虑信源均为近场时的情况。入射信号为2个4QAM近场窄带信号，其位置分别为(θ1=-13°,r1=3.5λ)和(θ1=21°,r1=4.5λ)。

在本实验中，验证本文算法DOA估计精度随信噪比和快拍数变化的情况。第一，信噪比从-5 dB到10 dB不断变化，步长为2 dB，快拍数为700，蒙特卡洛次数为500；第二，信噪比为7 dB，快拍数从100到1 100不断变化，步长为200。图7为DOA估计的角度均方根误差随信噪比的变化情况，图8为近场信源距离的均方根误差随信噪比的变化。图9和图10分别是方位角和距离随快拍数的变化曲线图。从图7可以看出，本文算法的信源角度估计精度高于对比算法。由于本文算法所采用稀疏对称阵列，其阵列孔径要略大于互质对称阵列[11]和特殊的几何阵列结构[13]，且远大于均匀线阵[12]，因此本文算法具有更高的角度估计精度。从图9和图10可以看出，4种算法在快拍数较小时估计精度较差，随着快拍数的增大性能有所改善并趋于稳定。从图8和图10可以看出，本文算法的近场信源距离估计精度在随着信噪比和快拍数变化下都高于对比算法。由于距离估计是基于角度估计，所提算法角度估计性能最好，所以也具有更高的距离估计精度。

4 结束语

本文在稀疏对称阵列模型的基础上，提出了一种混合信源DOA估计算法。通过不同子阵列接收数据的四阶累积量运算，构造一个仅与信源角度有关的特殊四阶累积量矩阵，然后使用MUSIC算法得到所有信源的角度估计；在每个角度估计的基础上进行距离维的搜索，进而估计出近场信源的距离参数。本文算法避免了二维搜索和参数配对，稀疏对称阵列扩展了阵列孔径，仿真实验结果表明，在相同的阵元情况下，本文算法具有更高的估计精度。

[1] MA Yan, CHEN Baixiao, YANG Minglei, et al. A Novel ESPRIT-Based Algorithm for DOA Estimation with Distributed Subarray Antenna[J]. Circuits Systems & Signal Processing, 2015 , 34(9):2951-2972.

[2] MA Xiurong, DONG Xuhao, XIE Yufeng. An Improved Spatial Differencing Method for DOA Estimation with the Coexistence of Uncorrelated and Coherent Signals[J]. IEEE Sensors Journal, 2016, 16(10):3719-3723.

[3]曾文浩, 朱晓华, 李洪涛,等. 一种稀疏阵列下的二维DOA估计方法[J]. 航空学报, 2016，37(7):2269-2275.

[4] ZHENG Zhi, SUN Jie, WANG Wenqin,et al. Classification and Localization of Mixed Near-Field and Far-Field Sources Using Mixed-Order Statistics[J]. Signal Processing,2018,143:134-139.

[5] PAIK J W, LEE J.Performance Analysis of ML-Based DOA Estimation Algorithm in Bistatic MIMO Radar System[C]∥2017 7th IEEE International Symposium on Microwave, Antenna, Propagation, and EMC Technologies (MAPE), Xi’an:IEEE,2017:543-548.

[6] SUN Fenggang, GAO Bin, CHEN Lizhen, et al. A Low-Complexity ESPRIT-Based DOA Estimation Method for Co-Prime Linear Arrays[J]. Sensors, 2016, 16(9):1367-1376.

[7] HUANG Y D, BARKAT M. Near-Field Multiple Source Localization by Passive Sensor Array[J]. IEEE Trans on Antennas & Propagation, 1991, 39(7):968-975.

[8]YUEN N, FRIEDLANDER B. Performance Analysis of Higher Order ESPRIT for Localization of Near-Field Sources[J]. IEEE Trans on Signal Processing, 1998, 46(3):709-719.

[9] LIANG J L, LIU D.Passive Localization of Mixed Near-Field and Far-Field Sources Using Two-Stage MUSIC Algorithm[J]. IEEE Trans on Signal Processing, 2010, 58(1):108-120.

[10] JIN H, SWAMY M N S, AHMAD M O. Efficient Application of MUSIC Algorithm Under the Coexistence of Far-Field and Near-Field Sources[J]. IEEE Trans on Signal Processing, 2012, 60(4):2066-2070.

[11] JIANG Jiajia, DUAN Fajie, CHEN Jin, et al. Mixed Near-Field and Far-Field Sources Localization Using the Uniform Linear Sensor Array[J]. IEEE Sensors Journal, 2013, 13(8):3136-3143.

[12] 梁国龙，韩博. 基于互素对称阵的近场源定位[J]. 电子与信息学报, 2014, 36(1):135-139.

[13] WANG Bo, ZHAO Yanping, LIU Juanjuan. Mixed-Order MUSIC Algorithm for Localization of Far-Field and Near-Field Sources[J]. IEEE Signal Processing Letters, 2013, 20(4):311-314.

[14] TANAKA K, KIKUMA N, SAKAKIBARA K. Influence of Mutual Coupling Between Array Elements

in Location Estimation of Radio Sources Using Near-Field DOA-Matrix Method[C]∥2016 International Symposium on Antennas and Propagation (ISAP), Okinawa,Japan:IEEE, 2016:1024-1025.

[15] 刘庆华,周秀清,晋良念.相关色噪声下无冗余累积量稀疏表示DOA估计[J].航空学报,2017,38(4):212-221.