在阵列信号处理的DOA估计方法研究过程中,色噪声的构造无统一格式,通常采用基于自回归(auto regressive, AR)模型、MA模型及自回归滑动平均(auto regressive moving average, ARMA)模型等方法,本研究选取基于MA模型来构造色噪声,具体定义如下:

假设入射信号与阵列输出数据中的色噪声不相关,不同阵元间的色噪声亦不相关,则基于L阶的MA模型的第m个阵元上色噪声的表达式^[12]为

n_m(n)=a₀w(n)+a₁w(n-1)+…+

a_Lw(n-L)(4)

其中, w(n)是零均值方差为σ²_w的高斯白噪声; m=1, 2, …, M.

由此可得n_m(n)的相关函数表达式为

r_nm(l)=E[n_m(n)n^*_m(n-l)]=

{σ²_w∑^l_i=0a_ia^*_i+1, l=0,1,…,L

0, l>L(5)

由式(5)可知,当基于MA模型色噪声相关函数的延时l大于色噪声的阶数L时,色噪声的延时相关函数值为0,即可以利用大于色噪声阶数的延时相关函数来消除阵列接收数据的色噪声.

然而,阵列的输出数据中色噪声并不是独立于入射信号而存在的,所以需要结合本小节的结论研究在阵列接收数据中去除色噪声的方法.

本研究以M元均匀线阵为例,研究阵列接收数据的延时相关函数,阵元间距为d, 第1个阵元是参考点,位于原点,设空间有P个入射信号,入射信号间不相关,则阵列的导向矢量可表示为

a(θ_i)=[1, exp(-j2πdsin θ_i/λ), …,

exp(-j2π(M-1)dsin θ_i/λ)]^T(6)

其中, θ_i为第i(i=1, 2, …, P)个信号的入射角度, λ为入射信号的波长.于是,阵列流型可以表示为

A=[a(θ₁),a(θ₂),…,a(θ_P)](7)

结合阵列输出数据的矩阵表达如式(1),可得到第m个阵元的接收数据表达式为

x_m(n)=A_mS(n)+n_m(n)(8)

结合式(6)和式(7),可得矩阵 A_m(m=1, 2, …, M)的表达式为

A_m=[a_m(θ₁), a_m(θ₂), …, a_m(θ_P)]=

[exp(-j2πdsin θ_i/λ), …,

exp(-j2π(M-1)sin θ_i/λ)]^T(9)

且有

A=[AT₁,AT₂,…,AT_M](10)

假设n_m(t)是基于L阶的MA模型的色噪声,阵元间噪声不相关,可得到第m(m=1, 2, …, M)个阵元的接收数据与第1个阵元之间的延时相关函数为

R_xmx1(l)=E[x_m(n)x^*₁(n-l)]=

∑^P_i=1exp[-j2π(m-1)dsin θ_i/λ]S_i(n)S^*_i(n-l)+

n_m(n)n^*₁(n-1)=A_mR_S(l)+r_m1(l)(11)

其中, R_S(l)=[R_S1(l),R_S2(l),…,R_SP(l)]^T为各入射信号的自延时相关函数矩阵; r_m1=n_m(n)n^*₁(n-l)为第m个阵元与第1个阵元间的色噪声延时相关函数.

根据前述色噪声延时相关函数的结论可知,若取延时l满足l >L,则可得式(11)中的r_m1(l)=0, 于是有

R_xmx1(l)=A_mR_S(l)(12)

在l >L的条件下把每个阵元与第1个阵元的延时相关函数写成矩阵形式为

R_Xx1(l)=[A₁^T, A₂^T, …, A_M^T]R_S(l)=

AR_S(l)(13)

参考上述推导,同样取l >L推导第m个阵元与第k个阵元间的延时相关函数为

R_xmxk(l)=E[x_m(n)x^*_k(n-l)]=

∑^P_i=1exp[-j2π(m-k)dsin θ_i/λ]S_i(n)S^*_i(n-l)+

n_m(n)n^*₁(n-l)=A_m-k+1R_S(l)=

A_mdiag(A^*_k)R_S(l)(14)

于是,阵列每个阵元的接收数据与第k个阵元间的延时相关函数写成矩阵形式为

R_Xxk(l)=Adiag(A^*_k)R_S(l)(15)

结合式(1)与式(15)可知,当延时l >L时,阵列的接收数据与第k个阵元间的延时相关函数可以看作阵列流型为 Adiag(A^*_k), 入射信号为 R_S(l)的一个新阵列的输出.于是可得阵列与第k个阵元的延时相关函数的导向矢量为

a_k=Adiag(A^*_k)_θi=

{exp[-j2π(1-k)sin θ_i/λ]

exp[-j2π(2-k)sin θ_i/λ] …

exp[-j2π(M-k)dsin θ_i/λ]}^T(16)

即新阵列的阵列流型为Adiag(A^*_k), 而导向矢量为其中的一列,角度是未知的θ_i. 同理,可得阵列与第M-k+1个阵元的延时相关函数的导向矢量为

a_M-k+1=Adiag(A^*_M-k+1)_θi=

[exp(j2π(M-k)sin θ_i/λ)…

exp(j2π(2-k)sin θ_i/λ)

exp(j2π(1-k)dsin θ_i/λ)]^T(17)

由式(16)和式(17)可知,阵列与第k个阵元和与其第M-k+1个阵元的延时相关函数的导向矢量满足

Adiag(A^*_k)_θi=flipup[Adiag(A^*_M-k+1)_θi]^*(18)

其中,flipup表示矩阵上下倒序翻转.

于是,若用倒序的阵列接收数据与第M-k+1个阵元求取延时相关函数,则有

R_X'xM-k+1(l)=flipup[Adiag(A^*_M-k+1)]R_S(l)=

[Adiag(A^*_k)]^*R_S(l)(19)

其中,X'表示阵列接收数据按阵元位置的反序.

于是可构造新的阵列协方差矩阵为

R_k(l)=R_Xxk(l)RH_Xxk(l)+R^*_X'xM-k+1(l)[R^*_X'xM-k+1(l)]^H=

Adiag(A^*_k)R_S(l)RH_S(l)[Adiag(A^*_k)]^H+Adiag(A^*_k)R^*_S(l)RT_S(l)[Adiag(A^*_k)]^H=

Adiag(A^*_k)[R_S(l)RH_S(l)+R^*_S(l)RT_S(l)][Adiag(A^*_k)]^H=

Adiag(A^*_k)R_SS(l)[Adiag(A^*_k)]^H(20)

其中,R_SS(l)=R_S(l)RH_S(l)+R^*_S(l)RT_S(l)=2Re[R_S(l)RH_S(l)].

由式(20)可得到新的协方差矩阵.对式(20)进行特征分解,仿照式(3),即可通过谱峰搜索得到信号的入射角度.由式(20)的构造又可知,有M个阵元的阵列,并可得到k组结果. 其中,当M为偶数时,k=1, 2, …, M/2; 当M为奇数时, k=1, 2, …,(M+1)/2. 将k组的搜索结果求取统计平均可以提高搜索精度.

由式(20)得到新的协方差矩阵后,对矩阵进行特征分解可得对应的信号子空间U'_S、噪声子空间 U'_N和特征值λ_m. 其中,当m=1,2, …, P时,λ_m为信号特征值; 当m=P+1, P+2, …, M时, λ_m为信号特征值. 对应第k(M-k+1)个阵元的导向矢量表达如式(16),仿照MUSIC算法可得到在色噪声情况下的入射信号角度搜索的谱函数为

P_CLR-MUSIC=1/(a_kU'_N(U'_N)^HaH_k)(21)

由式(21)可见,特征分解得到的信号子空间及特征值等结果并未在谱函数中有充分应用,而且信号子空间处理具有稳健性高的特点,文献[13]也证实了信号子空间的加权可以提高MUSIC算法谱函数的分辨能力,并给出了加权矩阵的选取方法.本研究仿照子空间加权的方法将式(21)谱函数改写为

P_CLR2nd-MUSIC=(a_kU'_SW_S(U'_S)^HaH_k)/(a_kU'_N(U'_N)^HaH_k)(22)

其中, W_S=diag((λ_m)/(λ_m-σ²_rms)),m=1, 2, …, P, 这里, σ²_rms是噪声的平均功率,实际估计过程中σ_rms是未知量.根据本研究的噪声背景通过噪声特征值对其进行估计,则

σ_rms=(1/(M-P)∑^M_m=P+1λ²_m)^1/2(23)

于是,在阵元间延时相关阶数大于色噪声阶数L的情况下,通过对接收数据的二阶预处理及新的协方差矩阵的构造,可得到去除色噪声的延时二阶预处理MUSIC方法(CLR2nd-MUSIC算法).

实验1 对3种方法的谱函数图进行比较.信号入射角度分别设为1.5°和-1.5°,实验中信噪比设置为20 dB,快拍数为500,1次实验的谱函数图如图1.其中,纵坐标P表示3种算法以dB为单位的谱值.

图1 三种方法的谱函数图
Fig.1 Spectrum functions of three methods

实验2 信噪比、快拍数对分辨能力的影响.当信号入射角度为1.5°和-1.5°,快拍数为500时,不同信噪比情况的分辨成功概率如图2.当信噪比为15 dB时,不同快拍数情况下的分辨成功概率如图3.

图2 不同信噪比条件下的分辨概率
Fig.2 Resolution probability with different SNR

图3 不同快拍数条件下的分辨概率
Fig.3 Resolution probability with different number of snaps

实验3 不同信号入射间隔对成功概率的影响.两个入射信号对称分布在视轴两侧,快拍数设置为500,信噪比为15 dB,仿真结果见图4.

图4 不同角度间隔的分辨概率
Fig.4 Resolution probability with different angle interval

综合实验1到实验3的仿真结果,显然在色噪声背景下,MUSIC算法难以实现对入射信号的分辨,而累积量方法及CLR2nd-MUSIC算法均能在没有孔径损失、不预先估计噪声空间的情况下,实现对入射信号的分辨且估计性能相近.

由算法的分辨力仿真结论可以看出,累积量方法及CLR2nd-MUSIC算法对色噪声均能体现出较好的抑制作用.为此,对算法的运算量复杂度进行分析,并通过仿真实验进行验证.

对实验中3种算法的运算复杂度进行协方差的计算、特征值分解及谱峰搜索3个主要步骤分析,见表1.

表1 三种算法的运算复杂度1)
Table 1 Computational complexity of three methods

由表1可见,与MUSIC算法相比,CLR2nd-MUSIC算法仅在构造协方差矩阵时有复数乘加运算次数的增加,其余步骤均无计算量的增加,而累积量算法在角度估计的各个步骤中均增加了大量的运算复杂度.

实验4 运算时间的比较.在计算机相同的运行条件下,对3种方法估计2个入射信号所消耗的时间进行统计,结果如表2.

表2 相同运行条件下3种方法的计算时间
Table 2 Computing times for three methods under same operation conditions