WO2021139208A1 - 一种基于特定频率组合信号的一维doa估计方法 - Google Patents

Abstract

一种基于特定频率组合信号的一维DOA估计方法，该方法通过发射多个特定频率的信号，将各阵元的接收信号分别通过多个不同频率窄带滤波器分离出对应频率信号，得到频域接收数据模型，再对源信号矩阵归一处理和重构频域接收数据矩阵，将均匀物理阵列一次扩展为虚拟非均匀阵列，最后结合现有非均匀阵列数据协方差矩阵的处理方法，对虚拟非均匀阵列进行二次虚拟扩展，最终获得二次扩展的虚拟均匀阵列，实现阵列孔径的二次扩张。该方法在无需改变传统均匀阵列结构的条件下，通过接收多个特定频率的窄带信号并依靠频率比例这一先验信息，实现了阵列孔径的虚拟扩展，使得估计精度得到提升。

Description

一种基于特定频率组合信号的一维DOA估计方法 技术领域

本发明涉及目标定位技术领域，具体涉及一种将特定频率组合信号联合处理后实现虚拟阵元扩展的一维DOA估计方法。

背景技术

波达方向估计(DOA估计)在众多领域已得到广泛应用，而DOA估计就是指在空间或放置传感器阵列利用阵列信号处理技术来对信源目标进行方向角度估计的方法。随着空间环境日趋复杂化，在目标跟踪、精密制导等各个方面确定目标信号的波达方向角并进行全面的监视和侦察，是在战争中获取主动权的首要前提。因此，研究超分辨高精度的DOA估计技术具有十分重要的意义。

现有的基于相关阵的子空间类方法如MUSIC算法和ESPRIT算法等，在低信噪比条件下性能会快速变差。为了提高阵列处理信号的能力而可以增加物理阵元数量，但这会导致硬件成本和维护成本的急剧增加。因此，如何实现在物理阵元数不变和低信噪比的条件下保证DOA估计算法的精度，成为阵列信号处理领域的一个研究难点和热点。有学者提出非均匀阵列结构，通过利用差协同阵元的更多相位差信息，使得阵列孔径得到提高，在相同物理阵元数的前提下提高了估计精度。但是非均匀阵列需要特定阵列的结构，使得传统均匀阵列失去了改造潜力。也有学者提出基于宽带信号的DOA估计算法，利用带宽内多个频点信号携带的信息，提高了DOA估计的稳健性。但是选取频点过少会导致估计精度下降，同时在相邻阵元间距保持一定的条件下，信号带宽过大会引入过多阵列孔径损失。

发明内容

本发明的目的是为了解决现有技术中的上述缺陷，提供一种基于特定频率组合信号的一维DOA估计方法。该方法在基于传统均匀阵列的前提下通过发射多个特定频率的信号，在每个接收阵元处得到的接收信号分别通过多个不同频率的窄带滤波器获取分离出对应频率信号，得到频域接收数据模型，再对源信号矩阵归一处理和重构频域接收数据矩阵，将均匀物理阵列一次扩展为虚拟非均匀阵列，最后结合现有非均匀阵列数据协方差矩阵的向量化处理，进行虚拟非均匀阵列的二次虚拟扩展，最终获得二次扩展的虚拟均匀阵列，实现了阵列孔径的二次扩张。

本发明的目的可以通过采取如下技术方案达到：

一种基于特定频率组合信号的一维DOA估计方法，所述的估计方法包括以下步骤：

S1、根据时域采样数据构造时域数据接收模型，再利用傅里叶变换得到频域数据接收模型和无噪声数据协方差矩阵

S2、根据最低频率数据进行MUSIC算法得到首次构建的角度预估值

S3、并根据角度预估值，构建归一矩阵G _p,q，利用无噪声数据协方差矩阵

和归一矩阵G _p,q重构出新的协方差矩阵

S4、对重构后的频域数据协方差矩阵

进行向量化得到接收向量r，实现基于虚拟非均匀阵列的二次阵元扩展；

S5、对向量r进行去冗余和重排处理，得到等效阵列信号接收向量r _u；

S6、对得到的接收向量r _u进行虚拟子阵划分和平均处理，得到空间平 滑后的协方差矩阵R _v，实现解相干；

S7、根据R _v采用MUSIC算法，得到第i次迭代时的目标角度估计值

若迭代次数达到最大次数或

满足收敛准则，得到最终估计值；否则将

作为新的角度预估值并跳转至步骤S3继续执行。

进一步地，所述的步骤S1过程如下：

假设用于一维DOA估计的均匀线阵上分布有M个阵元，阵元间距为d。设定P个频率的窄带发射信号，其中第p,p＝1,2,…,P路信号的频率为f _p，接收信号为发射信号经K个远场目标后返回阵列的信号，通过P个窄带滤波器，H个时域快拍数采样可以得到时域接收数据矩阵X，表示为：

其中X _p为第p路频率信号的接收数据，表示为：

X _p＝[x _p(1),x _p(2),…,x _p(H)]       (2)

x _p(t)＝[x _1,p(t),x _2,p(t),…,x _M,p(t)] ^T,t＝1,2,…,H        (3)

第m,m＝1,2,…,M个物理阵元的第p路窄带滤波器的接收快拍信号x _m,p(t)如下：

其中s _k,p(t)为第p路频率信号经第k个目标在t时刻到达阵列的源信号，n _m,p(t)为第m阵元第p路窄带滤波器的在t时刻的噪声信号，方向向量第m个元素a _m,p(θ _k)表示为：

其中q _m为第m个物理阵元位置，

为相位分量，d为阵列的阵元间距，θ _k为第k个信源目标的来波方向角度，f _p为第p路频率信号的频率，c为信号传播速度；

对于角度θ _k，阵列存在方向向量a _p(θ _k)：

a _p(θ _k)＝[a _1p(θ _k),a _2p(θ _k),…,a _Mp(θ _k)] ^T       (6)

对第m个物理阵元接收到的时域信号x _m,p(t)的H个时域快拍数据分成L段，进行傅里叶变换，得到频域快拍数据x _m,p,l(f),l＝1,2,…,L，将所有M个阵元的L个频域接收信号排成一个矩阵，得到有关频率f _p的所有数据的频域数据接收模型，即

其中

为源信号矩阵，

为噪声信号矩阵，

为阵列导向矩阵，表示为：

频率f _p对应的数据协方差矩阵R _x(f _p)表示为：

其中P _x(f _p)表示无噪声数据协方差矩阵，

表示源信号矩阵，

表示噪声矩阵；

对于数据协方差矩阵R _x(f _p)，其噪声功率

的估计值表示为：

其中λ _i(U)表示矩阵U特征分解后，第i大的特征值，因此无噪声数据协方差矩阵估计值

表示为：

其中I表示单位矩阵。

进一步地，所述的步骤S2过程如下：

将频率f ₁对应的无噪声数据协方差矩阵

进行特征分解，获得M个特征值，提取最小的M-K个特征向量构建出噪声子空间；

对于任意角度值θ，根据噪声子空间构造方向向量

构建MUSIC谱函数P _MUSIC(θ)，使θ从-90到90°变化，找到MUSIC空间谱的K个极大值对应的角度即

进一步地，所述的步骤S3过程如下：

无噪声数据协方差矩阵

对应的源信号矩阵估计值

表示为：

其中方向矩阵A _p通过引入角度预估值根据公式(8)构建，为了将频率f _p对应的源信号矩阵进行归一化处理，需要矩阵G _p,q满足下式：

因此G _p,q的解为：

根据

和G _pq重构出新的接收信号矩阵

其中

为重构后的方向矩阵，重构后的源信号矩阵为单位矩阵，

为重构后的噪声矩阵，B表示如下：

B＝[b(θ ₁),b(θ ₂)…,b(θ _K)]       (16)

其中b(θ _k)为重构方向向量：

进一步地，所述的步骤S4过程如下：

由于频率组合f＝[f ₁,f ₂,…,f _P] ^T，设最低频点为f _L，且所有频点都是f _L的整数倍，则有

f＝[g ₁f _L,g ₂f _L,…,g _Pf _L] ^T       (18)

其中g＝[g ₁,g ₂,…,g _P] ^T为频率比例向量，

设物理阵列位置y＝[y ₀,y ₁,…,y _M-1]，则经过矩阵重构，以f _L为参考频率的虚拟阵列位置z为

方向向量b(θ _k)重写为：

数据重构矩阵R _F为MP×MP维的协方差矩阵，第m行n列的元素表示为：

其中

是重构数据的噪声功率，δ _m,n为Kronecker delta函数，由式(21)可见，一次扩展阵列协方差矩阵R _F中的元素被视为虚拟差协同阵元的接收数据，为了进行二次扩展将协方差矩阵R _F向量化，得到：

其中，vec为向量化符号，⊙表示Khatri-Rao积，

噪声向量

B ^*⊙B表示为：

其中

表示Kronecker积。

进一步地，所述的步骤S5过程如下：

使用空间平滑技术解相干，对B ^*⊙B进行去冗余和重排处理，假设经 过去冗余和重排处理后的等效阵列信号接收向量为r _u，表示为：

其中Q为B ^*⊙B进行去冗余和重排处理后的方向矩阵，I _u是新的噪声向量，M阵元的均匀阵列经过接收数据重构一次扩展成虚拟非均匀阵列，再通过协方差矩阵矢量化，去冗余和重排处理后将二次扩展成一个新的虚拟均匀线阵，为方便阐述后续步骤，以频率比例向量

的频率组合为例，M阵元的物理阵列，得到虚拟均匀线阵位置分布范围为-(M-1)Md～(M-1)Md，即包含了2(M-1)Md+1个虚拟阵元。

进一步地，所述的步骤S6过程如下：

对接收向量r _u进行子阵划分，将二次扩展虚拟均匀阵列划分为(M-1)M+1个子阵，每个子阵包含(M-1)M+1个阵元，

第i，i＝1,2,…,(M-1)M+1个子阵的接收数据r _u,i等于阵列信号接收向量r _u的第i～(i+(M-1)M+1)列，令

计算所有R _i的值并进行平均得到空间平滑后的协方差矩阵R _v：

进一步地，所述的步骤S7过程如下：

将所获得的R _v进行特征分解，获得(M-1)M+1个特征值，把特征值依照大小实现顺序排列，提取最小的(M-1)M+1-K个特征向量构建出U _n，在搜索范围内选取角度θ构造方向向量a _v(θ)：

代入a _v(θ)得到MUSIC谱函数：

根据MUSIC谱函数，使θ从-90到90°变化，通过寻求K个极大值对 应的角度，即第i,i≥1次迭代的波达角估计值

存在收敛准则：

其中biase表示预先设定的允许偏差值，当i达到最大迭代次数或

满足上式，则

为最终角度估计值；否则将

作为新的角度预估值并跳转至步骤S3。

本发明相对于现有技术具有如下的优点及效果：

1、本发明公开的DOA估计方法无需增加物理阵元，而是通过多组不同频率信号之间的比例先验信息完成了数据协方差矩阵重构。重构矩阵的虚拟阵列等效于一组非均匀阵列，进一步对数据协方差矩阵进行向量化，实现了基于虚拟非均匀阵列的二次阵元扩展，最后的等效虚拟均匀阵列相较原物理阵列，阵列孔径增加，DOA估计精度也有较大提升。

2、本发明公开的DOA估计方法为了消除首次角度预估值带来的误差，通过迭代估计实现角度收敛，进一步提高了估计精度。

3、本发明公开的DOA估计方法在对接收向量进行去冗余过程中对位置相同的虚拟差分阵元等效接收数据进行平均处理，进一步提高了DOA估计方法的稳健性。

4、本发明公开的DOA估计方法只需在传统均匀阵列的基础上增加窄带滤波器组，即可适用于大部分现有的DOA估计方法的应用系统，工程可实现性比较高。

附图说明

图1是本发明实施例中公开的DOA估计方法所用的均匀线阵接收信号与窄带滤波器组模型示意图；

图2是本发明实施例中一维均匀阵元与窄带滤波器组的连接示意图；

图3是本发明实施例中公开的基于特定频率组合信号的一维DOA估计方法的流程图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例一

如附图3所示，本实施例公开了一种基于多频率发射信号的一维DOA估计方法包括以下步骤：

S1、根据时域采样数据构造时域数据接收模型，再利用傅里叶变换得到频域数据接收模型和无噪声数据协方差矩阵

假设用于一维DOA估计的均匀线阵上分布有M个阵元，阵元间距为d。设定P个频率的窄带发射信号，其中第p,p＝1,2,…,P路信号的频率为f _p。接收信号为发射信号经K个远场目标后返回阵列的信号，通过P个窄带滤波器，H个时域快拍数采样可以得到时域接收数据矩阵X，表示为：

其中X _p为第p路频率信号的接收数据，表示为：

X _p＝[x _p(1),x _p(2),…,x _p(H)]        (2)

x _p(t)＝[x _1,p(t),x _2,p(t),…,x _M,p(t)] ^T,t＝1,2,…,H     (3)

第m,m＝1,2,…,M个物理阵元的第p路窄带滤波器的接收快拍信号x _m,p(t)如下：

其中s _k,p(t)为第p路频率信号经第k个目标在t时刻到达阵列的源信号， n _m,p(t)为第m阵元第p路窄带滤波器的在t时刻的噪声信号。方向向量第m个元素a _m,p(θ _k)表示为：

其中q _m为第m个物理阵元位置，

为相位分量，d为阵列的阵元间距，θ _k为第k个信源目标的来波方向角度，f _p为第p路频率信号的频率，c为信号传播速度。

对于角度θ _k，阵列存在方向向量a _p(θ _k)：

a _p(θ _k)＝[a _1p(θ _k),a _2p(θ _k),…,a _Mp(θ _k)] ^T     (6)

对第m个物理阵元接收到的时域信号x _m,p(t)的H个时域快拍数据分成L段，进行傅里叶变换，得到频域快拍数据x _m,p,l(f),l＝1,2,…,L，将所有M个阵元的L个频域接收信号排成一个矩阵，得到有关频率f _p的频域数据接收模型，即

其中

为源信号矩阵，

为噪声信号矩阵，

为阵列导向矩阵，表示为：

频率f _p对应的数据协方差矩阵R _x(f _p)表示为：

其中P _x(f _p)表示无噪声数据协方差矩阵，

表示源信号矩阵，

表示噪声矩阵。

对于数据协方差矩阵R _x(f _p)，其噪声功率

的估计值可以表示为：

其中λ _i(U)表示矩阵U特征分解后，第i大的特征值。因此无噪声数据协方 差矩阵估计值

可以表示为：

其中I表示单位矩阵。

S2、根据最低频率数据，通过MUSIC算法得到首次构建的角度预估值

将频率f ₁对应的无噪声数据协方差矩阵

进行特征分解，获得M个特征值，提取最小的M-K个特征向量构建出噪声子空间。对于任意角度值θ，根据噪声子空间构造方向向量

进一步构建MUSIC谱函数P _MUSIC(θ)。使θ从-90到90°变化，找到MUSIC空间谱的K个极大值对应的角度即

S3、并根据角度预估值，构建归一矩阵G _p,q，利用无噪声数据协方差矩阵

和归一矩阵G _p,q重构出新的协方差矩阵

无噪声数据协方差矩阵

对应的源信号矩阵估计值

表示为：

其中方向矩阵A _p通过引入角度预估值根据公式(8)构建。为了将频率f _p对应的源信号矩阵进行归一化处理，需要矩阵G _p,q满足下式：

因此G _p,q的解为：

根据

和G _pq重构出新的接收信号矩阵

其中

为重构后的方向矩阵，重构后的源信号矩阵为单位矩阵，

为重构后的噪声矩阵。B表示如下：

B＝[b(θ ₁),b(θ ₂)…,b(θ _K)]       (16)

其中b(θ _k)为重构方向向量：

S4、对重构后的频域数据协方差矩阵

进行向量化得到接收向量r，实现基于虚拟非均匀阵列的二次阵元扩展：

由于频率组合f＝[f ₁,f ₂,…,f _P] ^T，设最低频点为f _L，且所有频点都是f _L的整数倍，则有

f＝[g ₁f _L,g ₂f _L,…,g _Pf _L] ^T        (18)

其中g＝[g ₁,g ₂,…,g _P] ^T为频率比例向量。

设物理阵列位置y＝[y ₀,y ₁,…,y _M-1]，则经过矩阵重构，以f _L为参考频率的虚拟阵列位置z为

方向向量b(θ _k)重写为：

数据重构矩阵R _F为MP×MP维的协方差矩阵，第m行n列的元素表示为：

其中

是重构数据的噪声功率，δ _m,n为Kronecker delta函数。由式(21)可见，一次扩展阵列协方差矩阵R _F中的元素被视为虚拟差协同阵元的接收数据。为了进行二次扩展将协方差矩阵R _F向量化，得到：

其中，vec为向量化符号，⊙表示Khatri-Rao积，

噪声向量

B ^*⊙B表示为：

其中

表示Kronecker积。

通过将一次扩展虚拟阵列接收数据的协方差矩阵向量化，转化为二次扩展虚拟阵列下的等效单快拍接收数据，等效阵元数增加。

S5、对向量r进行去冗余和重排处理，得到等效阵列信号接收向量r _u：以频率比例向量

的频率组合为例，则经过协方差矩阵向量化得到r后，B ^*⊙B中只有2(M-1)M+1行是互不相同的，即二次扩展虚拟阵列的阵元数为2(M-1)M+1。此时的入射信号相当于K个相干信号入射。为了使用空间平滑技术解相干，对B ^*⊙B进行去冗余和重排处理，使处理后的方向矩阵与虚拟均匀线阵的方向矩阵对应。为了提高算法的稳健性，本方法在去冗余过程中需要对位置相同的虚拟阵元等效接收数据做平均求和处理。假设经过去冗余和重排处理后的等效阵列信号接收向量为r _u，表示为：

其中

为B ^*⊙B进行去冗余和重排处理后的方向矩阵，

是新的噪声向量。

M阵元的均匀阵列经过接收数据重构一次扩展成一个虚拟非均匀阵列，再通过协方差矩阵矢量化，去冗余和重排处理后二次扩展成新的虚拟均匀线阵。虚拟均匀线阵的位置分布范围为-(M-1)Md～(M-1)Md，即包含了2(M-1)Md+1个虚拟阵元。

S6、对得到的接收向量r _u进行虚拟子阵划分和平均处理，得到空间平滑后的协方差矩阵R _v，实现解相干：

对接收向量r _u进行子阵划分，这里将二次扩展虚拟均匀阵列划分为(M-1)M+1个子阵，每个子阵包含(M-1)M+1个阵元，

第i，i＝1,2,…,(M-1)M+1个子阵的接收数据r _u,i等于阵列信号接收向量r _u的第i～(i+(M-1)M+1)列。令

计算所有R _i的值并进行平均就可以得到空间平滑后的协方差矩阵R _v：

协方差矩阵R _v为((M-1)M+1)×((M-1)M+1)维矩阵，当信源目标个数K满足K≤(M-1)M+1时，基于空间平滑的估计方法仍然能分辨出全部K个信源目标。

S7、基于R _v采用MUSIC算法得到目标角估计值

若迭代次数达到最大次数或

满足收敛准则，得到最终估计值；否则将

作为新的角度预估值并跳转至S3：将所获得的R _v进行特征分解，获得(M-1)M+1个特征值，把特征值依照大小实现顺序排列，提取最小的(M-1)M+1-K个特征向量构建出U _n。在搜索范围内选取角度θ构造方向向量a _v(θ)：

代入a _v(θ)得到MUSIC谱函数：

根据MUSIC谱函数，使θ从-90到90°变化，通过寻求K个极大值对应的角度，即第i,i≥1次迭代的波达角估计值

存在收敛准则：

其中biase表示预先设定的允许偏差值。当i达到最大迭代次数或

满足上式，则

为最终角度估计值；否则将

作为新的角度预估值并跳转至步骤S3。

实施例二

本实施例公开了一种基于多频率发射信号的一维DOA估计方法，具体工作步骤如下：

T1、假设均匀线阵中的阵元个数M统一定为6。

发射信号为

频率数P＝2，频率组合取f＝[f ₁,f ₂] ^T＝[1,6] ^TkHz；声速取c＝1500m/s，组合信号数据经重构后频率统一为f ₁＝1kHz，因此最小半波长为0.75m。任意两相邻线阵之间的距离必须小于等于0.75m，在满足此限制条件下可以任意选取阵元间距，所以设置两个均匀线阵的平均间距取0.75m，即第一个阵元和最后一个阵元相隔3.75m。在水下放置K＝2个目标声源，信源目标信号入射的波达方向角分别为(20°，60°)。

T2、进行多次采样；每次采样均匀线阵接收到的信号为x(t)＝[x _0,1(t),x _1,1(t),…,x _4,1(t),x _5,1(t),x _0,2(t),x _1,2(t),…,x _5,2(t)]，共时域采样接收1024次，频域快拍数设置为256，并将接收到的信号进行分析处理。

T3、设定MUSIC算法最大迭代次数48，收敛偏差0.02°，MUSIC算法处理步骤具体如下：

T31、将时域采样数据分为L段，对每个子段数据进行离散傅里叶变换得到频域接收数据

根据X _p构造频域数据协方差矩阵R _x(f _p)，根据公式(11)得到无噪声数据协方差矩阵

T32、利用一维MUSIC算法得到角度预估值

T33、根据公式(12)求出源信号矩阵估计值

根据公式(14)求出频率数据两两对应的归一矩阵G _p,q，根据公式(15)求得重构协方差矩阵

T34、对

进行向量化得到接收向量r，即虚拟非均匀阵列二次扩展得到2(M-1)M+1个虚拟均匀阵元；

T35、根据公式(24)得到接收向量r _u，完成B ^*⊙B虚拟阵列方向矩阵去冗余和重排序；

T36、对r _u进行子阵划分和平均处理，按式(25)得到空间平滑后的最终协方差矩阵R _v；

T37、按照式(27)进行谱峰搜索，得到第i次迭代的角度估计值

若迭代次数达到最大次数或

满足收敛准则(28)，则得到最终估计值；否则将这组

作为新的角度预估值跳转至步骤T33。

步骤T4、根据本实施例算法，经12次迭代估计出的二维波达方向角(20.15°,60.03°)，对目标估计达到了预期精度，说明估计结果正确，本发明方法可行。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (8)

1. 一种基于特定频率组合信号的一维DOA估计方法，其特征在于，所述的估计方法包括以下步骤：

   S1、根据时域采样数据构造时域数据接收模型，再利用傅里叶变换得到频域数据接收模型和无噪声数据协方差矩阵

   

   S2、根据最低频率数据进行MUSIC算法得到首次构建的角度预估值

   

   S3、并根据角度预估值，构建归一矩阵G _p,q，利用无噪声数据协方差矩阵

   

   和归一矩阵G _p,q重构出新的协方差矩阵

   

   S4、对重构后的频域数据协方差矩阵

   

   进行向量化得到接收向量r，实现基于虚拟非均匀阵列的二次阵元扩展；

   S5、对向量r进行去冗余和重排处理，得到等效阵列信号接收向量r _u；

   S6、对得到的接收向量r _u进行虚拟子阵划分和平均处理，得到空间平滑后的协方差矩阵R _v，实现解相干；

   S7、根据R _v采用MUSIC算法，得到第i次迭代时的目标角度估计值

   

   若迭代次数达到最大次数或

   

   满足收敛准则，得到最终估计值；否则将

   

   作为新的角度预估值并跳转至步骤S3继续执行。
2. 根据权利要求1所述的一种基于特定频率组合信号的一维DOA估计方法，其特征在于，所述的步骤S1过程如下：

   假设用于一维DOA估计的均匀线阵上分布有M个阵元，阵元间距为d。设定P个频率的窄带发射信号，其中第p,p＝1,2,…,P路信号的频率为f _p，接收信号为发射信号经K个远场目标后返回阵列的信号，通过P个窄带滤波器，H个时域快拍数采样可以得到时域接收数据矩阵X，表示为：

   

   其中X _p为第p路频率信号的接收数据，表示为：

   X _p＝[x _p(1),x _p(2),…,x _p(H)]  (2)

   x _p(t)＝[x _1,p(t),x _2,p(t),…,x _M,p(t)] ^T,t＝1,2,…,H  (3)

   第m,m＝1,2,…,M个物理阵元的第p路窄带滤波器的接收快拍信号x _m,p(t)如下：

   

   其中s _k,p(t)为第p路频率信号经第k个目标在t时刻到达阵列的源信号，n _m,p(t)为第m阵元第p路窄带滤波器的在t时刻的噪声信号，方向向量第m个元素a _m,p(θ _k)表示为：

   

   其中q _m为第m个物理阵元位置，

   

   为相位分量，d为阵列的阵元间距，θ _k为第k个信源目标的来波方向角度，f _p为第p路频率信号的频率，c为信号传播速度；

   对于角度θ _k，阵列存在方向向量a _p(θ _k)：

   a _p(θ _k)＝[a _1p(θ _k),a _2p(θ _k),…,a _Mp(θ _k)] ^T  (6)

   对第m个物理阵元接收到的时域信号x _m,p(t)的H个时域快拍数据分成L段，进行傅里叶变换，得到频域快拍数据x _m,p,l(f),l＝1,2,…,L，将所有M个阵元的L个频域接收信号排成一个矩阵，得到有关频率f _p的所有数据的频域数据接收模型，即

   

   其中

   

   为源信号矩阵，

   

   为噪声信号矩阵，

   

   为阵列导向矩阵，表示为：

   

   频率f _p对应的数据协方差矩阵R _x(f _p)表示为：

   

   其中P _x(f _p)表示无噪声数据协方差矩阵，

   

   表示源信号矩阵，

   

   表示噪声矩阵；

   对于数据协方差矩阵R _x(f _p)，其噪声功率

   

   的估计值表示为：

   

   其中λ _i(U)表示矩阵U特征分解后，第i大的特征值，因此无噪声数据协方差矩阵估计值

   

   表示为：

   

   其中I表示单位矩阵。
3. 根据权利要求2所述的一种基于特定频率组合信号的一维DOA估计方法，其特征在于，所述的步骤S2过程如下：

   将频率f ₁对应的无噪声数据协方差矩阵

   

   进行特征分解，获得M个特征值，提取最小的M-K个特征向量构建出噪声子空间；

   对于任意角度值θ，根据噪声子空间构造方向向量

   

   构建MUSIC谱函数P _MUSIC(θ)，使θ从-90到90°变化，找到MUSIC空间谱的K个极大值对应的角度即

   
4. 根据权利要求3所述的一种基于特定频率组合信号的一维DOA估计方法，其特征在于，所述的步骤S3过程如下：

   无噪声数据协方差矩阵

   

   对应的源信号矩阵估计值

   

   表示为：

   

   其中方向矩阵A _p通过引入角度预估值根据公式(8)构建，为了将频率f _p对应 的源信号矩阵进行归一化处理，需要矩阵G _p,q满足下式：

   

   因此G _p,q的解为：

   

   根据

   

   和G _pq重构出新的接收信号矩阵

   

   

   其中

   

   为重构后的方向矩阵，重构后的源信号矩阵为单位矩阵，

   

   为重构后的噪声矩阵，B表示如下：

   B＝[b(θ ₁),b(θ ₂)…,b(θ _K)]  (16)

   其中b(θ _k)为重构方向向量：

   
5. 根据权利要求4所述的一种基于特定频率组合信号的一维DOA估计方法，其特征在于，所述的步骤S4过程如下：

   根据频率组合f＝[f ₁,f ₂,…,f _P] ^T，设最低频点为f _L，且所有频点都是f _L的整数倍，则有

   f＝[g ₁f _L,g ₂f _L,…,g _Pf _L] ^T  (18)

   其中g＝[g ₁,g ₂,…,g _P] ^T为频率比例向量，

   设物理阵列位置y＝[y ₀,y ₁,…,y _M-1]，则经过矩阵重构，以f _L为参考频率的虚拟阵列位置z为

   

   方向向量b(θ _k)重写为：

   

   数据重构矩阵R _F为MP×MP维的协方差矩阵，第m行n列的元素表示为：

   

   其中

   

   是重构数据的噪声功率，δ _m,n为Kronecker delta函数，由式(21)可见，一次扩展阵列协方差矩阵R _F中的元素被视为虚拟差协同阵元的接收数据，为了进行二次扩展将协方差矩阵R _F向量化，得到：

   

   其中，vec为向量化符号，⊙表示Khatri-Rao积，

   噪声向量

   

   B ^*⊙B表示为：

   

   其中

   

   表示Kronecker积。
6. 根据权利要求5所述的一种基于特定频率组合信号的一维DOA估计方法，其特征在于，所述的步骤S5过程如下：

   使用空间平滑技术解相干，对B ^*⊙B进行去冗余和重排处理，假设经过去冗余和重排处理后的等效阵列信号接收向量为r _u，表示为：

   

   其中Q为B ^*⊙B进行去冗余和重排处理后的方向矩阵，I _u是新的噪声向量，M阵元的均匀阵列经过接收数据重构一次扩展成虚拟非均匀阵列，再通过协方差矩阵矢量化，去冗余和重排处理后将二次扩展成一个新的虚拟均匀 线阵，以频率比例向量

   

   的频率组合为例，M阵元的物理阵列，得到虚拟均匀线阵位置分布范围为-(M-1)Md～(M-1)Md，即包含2(M-1)Md+1个虚拟阵元。
7. 根据权利要求6所述的一种基于特定频率组合信号的一维DOA估计方法，其特征在于，所述的步骤S6过程如下：

   对接收向量r _u进行子阵划分，将二次扩展虚拟均匀阵列划分为(M-1)M+1个子阵，每个子阵包含(M-1)M+1个阵元，

   第i，i＝1,2,…,(M-1)M+1个子阵的接收数据r _u,i等于阵列信号接收向量r _u的第i～(i+(M-1)M+1)列，令

   

   计算所有R _i的值并进行平均得到空间平滑后的协方差矩阵R _v：

   
8. 根据权利要求7所述的一种基于特定频率组合信号的一维DOA估计方法，其特征在于，所述的步骤S7过程如下：

   将所获得的R _v进行特征分解，获得(M-1)M+1个特征值，把特征值依照大小实现顺序排列，提取最小的(M-1)M+1-K个特征向量构建出U _n，在搜索范围内选取角度θ构造方向向量a _v(θ)：

   

   代入a _v(θ)得到MUSIC谱函数：

   

   根据MUSIC谱函数，使θ从-90到90°变化，通过寻求K个极大值对应的角度，即第i,i≥1次迭代的波达角估计值

   

   存在收敛准则：

   

   其中biase表示预先设定的允许偏差值，当i达到最大迭代次数或

   

   满足上式，则

   

   为最终角度估计值；否则将

   

   作为新的角度预估值并跳转至步骤S3。

Images