CN110133574B - 利用多频信号二次虚拟扩展的一维doa估计方法 - Google Patents

Abstract

本发明公开了一种利用多频信号二次虚拟扩展的一维DOA估计方法，该估计方法基于均匀阵列结构、利用多频率组合信号、通过二次虚拟扩展来提高阵列自由度，该方法的基本原理是通过发射多个特定频率的窄带信号，然后在接收阵元处通过窄带滤波器分离出对应频率的时域信号，经傅里叶变换得到频域接收数据模型，接着利用频率聚焦技术和数据协方差矩阵重构，将均匀物理阵列一次扩展为虚拟非均匀阵列，最后结合现有基于非均匀阵列的DOA估计方法，进行向量化处理。该方法实现了虚拟非均匀阵列的二次虚拟扩展，获得最终的虚拟均匀阵列。本发明通过以上技术手段实现了阵列孔径的二次扩展和DOA估计精度的提高。

Description

利用多频信号二次虚拟扩展的一维DOA估计方法

技术领域

本发明涉及目标定位技术领域，具体涉及一种利用多频率组合信号实现二次虚拟阵元扩展的一维波达方向估计方法。

背景技术

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

阵列自由度是衡量一个阵列能够同时处理的最大信号数的重要指标，传统的基于均匀阵的子空间类方法如MUSIC算法等只能处理的信源目标数不超过物理阵元数。而为提高阵列处理信号的能力而增加物理阵元将导致硬件和维护成本的急剧增加。因此，如何实现采用物理阵元数少于信号数的阵列进行波达方向估计即欠定DOA估计问题，成为阵列信号处理领域的一个研究热点。经研究发现，阵列的空间结构在欠定条件下决定了在给定阵元数的情况下能够获得多少阵列自由度。随后多种结构的非均匀阵列被提出，基于非均匀阵列的DOA估计方法有效提高了阵列的自由度。非均匀阵列对阵列设计要求较高，不同结构选择很大程度影响了实际自由度。

发明内容

本发明的目的是为了降低硬件成本和进一步提升DOA估计精度，提供一种基于均匀阵列结构、利用多频率组合信号、通过二次虚拟扩展来提高阵列自由度的一维DOA估计方法，该方法的基本原理是通过发射多个特定频率的窄带信号，然后在接收阵元处通过窄带滤波器分离出对应频率的时域信号，经傅里叶变换得到频域接收数据模型，接着利用频率聚焦技术和数据协方差矩阵重构，将均匀物理阵列一次扩展为虚拟非均匀阵列，最后结合现有基于非均匀阵列的DOA估计方法，进行向量化处理。该方法实现了虚拟非均匀阵列的二次虚拟扩展，获得最终的虚拟均匀阵列。本发明通过以上技术手段实现了阵列孔径的二次扩展和DOA估计精度的提高。

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

一种利用多频信号二次虚拟扩展的一维DOA估计法，所述的估计方法包括以下步骤：

S1、根据用于一维DOA估计的一维均匀线阵结构特点构造多频率信号时域数据接收模型，同时利用傅里叶变换得到多频率信号频域数据接收模型；

S2、根据频域数据接收模型对不同频点的接收数据构造聚焦矩阵，进行频率聚焦，同时完成频域数据接收模型的重构，实现第一次虚拟阵元扩展；

S3、对重构后的频域数据协方差矩阵R_F进行向量化得到向量r，实现基于一次扩展虚拟非均匀阵列的二次虚拟阵元扩展；

S4、对向量r进行冗余和重排处理得到等效阵列信号接收向量r₁，实现以B^*⊙B为整体作为向量化方向矩阵与虚拟均匀阵列的方向矩阵的一一对应关系，同时便于下一步的空间平滑处理；

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

S6、基于协方差矩阵R_s采用一维MUSIC算法估计全部K个信源目标的来波方向角度。

进一步地，所述的步骤S1中得到时域数据接收模型和时域数据接收模型的过程如下：

用于一维DOA估计的一维均匀线阵结构如图3所示，假设一维均匀线阵上分布有M个阵元，阵元间距为d，接收信号为由P个不同频率的窄带信号组成的发射信号经K个信源目标反射到接收阵列的信号，窄带信号满足窄带条件，即当信号延迟远小于带宽倒数时，延迟作用相当于使基带信号产生一个相移。当有K个远场信源目标的源信号入射到一维均匀线阵，在接收端通过接收阵元和对应窄带滤波器，经H个时域快拍数采样可以得到时域接收数据矩阵X。数据矩阵X由P_,P≥2个频率信号的接收数据矩阵叠加，表示为：

其中X_p为第p,p＝0,1,…,P-1路频率信号的接收数据表达为：

X_p＝[x_p(1),x_p(2),…,x_p(H)] (2)

x_p(t)＝[x_1p(t),x_2p(t),…,x_Mp(t)]^T,t＝1,2,…,H (3)

接收端第m,m＝1,2,…,M阵元的第p路窄带滤波器得到的信号x_mp(t)如下：

其中s_pk(t)为第p路频率发射信号经第k个目标在t时刻到达阵列的信号，n_mp(t)为第m阵元第p路窄带滤波器的在t时刻的噪声信号。

a_mp(θ_k)可以表达成：

其中q_m为第m个阵元的位置，

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

那么对于角度θ_k，阵列存在方向向量a_p(θ_k)：

a_p(θ_k)＝[a_p(θ_k),a_2p(θ_k),…,a_Mp(θ_k)]^T (6)

对时域接收信号x_mp(t)的H个时域快拍数据分成L段。假设H能被L整除，则第l段时域数据有：

进行傅里叶变换，得到其频域形式

其中x_mp,l(f)为一个频域数据，是

的傅里叶变换形式。s_pk,l(f)和n_mp,l(f)分别是第l个频域快拍下的源信号和噪声信号,频域快拍数为L。

将所有M个阵元的L个频域接收信号排成一个矩阵，就得到有关频率f_p的所有数据的频域数据接收模型，即

其中方向矩阵A(f_p,θ)，信源矩阵S(f_p)和噪声矩阵N(f_p)可以表示为：

其中s_l(f_p)表示S(f_p)的第l列向量，n_l(f_p)表示N(f_p)的第l列向量。

进一步地，所述的步骤S2中对不同频点的接收数据进行频率聚焦，同时完成频域数据接收模型重构的过程如下：

将各个频点的阵列流型通过“聚焦”变换到参考频率点f₀上，最后对参考频率点的数据采用窄带处理方式进行DOA估计。对应每个频率构造一个聚焦变换矩阵T(f_p)，使得不同频率点上的阵列流型矩阵聚焦到聚焦频点f₀上，即：

T(f_p)A(f_p,θ)＝A(f₀,θ) (13)

这里可以采用现有的双边相关变换(TCT)法求出聚焦矩阵T(f_p)。接下来，使用聚焦矩阵对接收信号进行变换，将变换之后接收信号表示为Y(f_p)，则有：

Y(f_p)＝T(f_p)X(f_p) (14)

从上式可以看出，通过“聚焦”使得不同频率点信号拥有了同样的阵列流型矩阵，即拥有同一个信号子空间。接着求出对应频点聚焦变换之后的信号协方差矩阵，并对P个频率分量聚焦变换后信号协方差矩阵进行平均求和可得到最终的聚焦协方差矩阵R_Y：

其中R_s(f_p)和R_n(f_p)表示为：

为了更好地利用多个频率数据的信息，接下来进行接收协方差矩阵的扩展重构。重构出的接收协方差矩阵R_F表达为：

其中，协方差矩阵R_F的第i·M+1～(i+1)·M,i＝0,1,…,P-1行，j·M+1～(j+1)·M,j＝0,1,…,P-1列的部分表达式为：

则协方差矩阵R_F最终可以重写为：

R_F＝B·R_FS·B^H+R_FN (20)

其中B＝[A(f₀,θ),…,A(f_P-1,θ)]^T∈C^MP×K为重构后的方向矩阵，R_FS∈C^K×K为重构后的源信号矩阵，R_FN∈C^MP×MP为重构后的噪声矩阵。

重构方向矩阵B同时可以表示为：

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

其中b(θ_k)表示重构方向向量，可以表达成：

由于频率数P≥2，本发明以P＝2为例，构造类似二级嵌套阵列的虚拟非均匀阵列。设定第p路频率f_p存在比例关系：

f₁＝M·f₀ (23)

那么方向向量b(θ_k)可以重写为：

其中，z_w表示一次扩展后的第w,w＝1,2,…,MP虚拟阵元的位置，

由上式可见，物理阵列通过重构扩展成一个非均匀的虚拟阵列。

进一步地，所述的步骤S3中对协方差矩阵R_F进行向量化得到向量r的过程如下：

R_F为MP×MP维的协方差矩阵，第u行v列的元素可以表示为：

其中

是重构数据的第k个源信号功率，

是重构数据的噪声功率，δ_u,v为Kronecker delta函数。

一次扩展虚拟阵列的协方差矩阵中的元素可以视为是以差协同阵列作为阵元位置的阵列接收的数据，因此可以将一次扩展虚拟阵列的协方差矩阵R_F向量化，得到：

其中vec为向量化符号，⊙表示Khatri-Rao积，g为源信号功率向量，噪声向量

向量化方向矩阵B^*⊙B表示为：

其中

表示Kronecker积。

向量化后的二次虚拟阵列中阵元数远多于物理阵列中的阵元数，通过将一次扩展虚拟阵列接收数据的协方差矩阵向量化的方式，转化为二次扩展虚拟阵列下的等效单快拍接收数据，阵元数再次大大增加，故而达到提高自由度的目的。

进一步地，所述的步骤S4中对步骤S3得到的向量r进行冗余和重排处理得到等效阵列信号接收向量r₁的过程如下：

经过协方差矩阵矢量化得到r后，向量化方向矩阵B^*⊙B中只有2(M-1)M+1行是互不相同的，即互不相同行数与二次扩展虚拟阵列的阵元数相同。此时的入射信号g是K个源信号的二阶统计值，即功率值，因此新的信号向量g相当于K个相干信号入射。为了使用空间平滑技术解相干，对向量化方向矩阵B^*⊙B进行去冗余和重排处理，使处理后的方向矩阵与虚拟均匀线阵的方向矩阵对应。假设经过去冗余和重排处理后的等效阵列信号接收向量为r₁，表示为：

其中方向矩阵A₂∈C^{(2(M-1)M+1)×K}为B^*⊙B进行去冗余和重排处理后的方向矩阵，

为新的噪声向量。

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

进一步地，所述的步骤S5中对步骤S4得到的接收向量r₁进行空间平滑处理后，得到协方差矩阵R_s的过程如下：

对接收向量r₁进行子阵划分处理。这里将二次扩展虚拟均匀阵列划分为(M-1)M+1个子阵，每个子阵包含(M-1)M+1个阵元，第h子阵的接收数据r_1h等于阵列信号接收向量r₁的第h～(h+(M-1)M+1)列，其中，h＝1,2,…,(M-1)M+1，令

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

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

进一步地，所述的步骤S6中根据步骤S5得到的R_s采用一维MUSIC算法估计K个信源目标的来波方向角度的过程如下：

将所获得的R_s实现一次特征分解，获得(M-1)M+1个特征值，把特征值依照大小实现顺序排列，提取最小的(M-1)M+1-K个特征向量构建出U_n。再构造方向向量a_s(θ)，可表达为：

把a_s(θ)代入得到MUSIC谱函数：

根据MUSIC谱函数，使θ从-90到90°变化，通过寻求K个极大值，对应的角度即为信源目标的波达角估计值。

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

1、本发明公开的一维DOA估计方法在物理阵元数不变的基础上，基于多组不同频率接收信号进行数据重构实现一次虚拟阵元扩展，再利用虚拟非均匀阵列的数据协方差矩阵进行向量化，实现了二次虚拟阵元扩展，相较基于原物理阵列的普通DOA估计方法，阵列孔径大大增加，估计精度也能大幅提升。

2、本发明公开的一维DOA估计方法相较于基于非均匀物理阵列(以二级嵌套阵列为例)的一维DOA估计方法，虚拟阵元扩展的程度更高。在相同物理阵元数的条件下，本发明所提方法只要增加发射信号频率数P，就能够扩展出更多的虚拟阵元。

3、在增加窄带滤波器组之后，本发明所提方法适用于大部分现有的基于传统均匀阵列的DOA估计方法，工程可实现性比较高。

附图说明

图1是本发明装置的硬件结构模块图；

图2是本发明中阵列的接收阵元与处理器连接示意图；

图3是本发明所用的均匀线阵接收信号与窄带滤波器组模型示意图；

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

图5是本发明公开的领域多频信号二次虚拟扩展的一维DOA估计方法的流程图。

具体实施方式

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

实施例一

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

S1、根据用于一维DOA估计的一维均匀线阵结构特点构造多频率信号时域数据接收模型，同时利用傅里叶变换得到多频率信号频域数据接收模型。

用于一维DOA估计的一维均匀线阵结构如图3所示，线阵上分布有M个阵元，阵元间距为d。设接收信号为由P个不同频率的窄带信号组成的发射信号经K个信源目标反射到接收阵列的信号,窄带信号满足窄带条件，即当信号延迟远小于带宽倒数时，延迟作用相当于使基带信号产生一个相移。当有K个远场信源目标的源信号入射到均匀阵列，在接收端通过接收阵元和对应窄带滤波器，经H个时域快拍数采样可以得到时域接收数据矩阵X。数据矩阵X由P,P≥2个频率信号的接收数据矩阵叠加，表示为：

其中X_p为第p,p＝0,1,…,P-1路频率信号的接收数据表达为：

X_p＝[x_p(1),x_p(2),…,x_p(H)] (2)

x_p(t)＝[x_1p(t),x_2p(t),…,x_Mp(t)]^T,t＝1,2,…,H (3)

接收端第m,m＝1,2,…,M阵元的第p路窄带滤波器得到的信号x_mp(t)如下：

其中s_pk(t)为第p路频率发射信号经第k个目标在t时刻到达一维均匀线阵的信号，n_mp(t)为第m阵元第p路窄带滤波器的在t时刻的噪声信号。

a_mp(θ_k)可以表达成：

其中q_m为第m个阵元的位置，

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

那么对于角度θ_k，阵列存在方向向量a_p(θ_k)：

a_p(θ_k)＝[a_p(θ_k),a_2p(θ_k),…,a_Mp(θ_k)]^T (6)

对时域接收信号x_mp(t)的H个时域快拍数据分成L段。假设H能被L整除，则第l段时域数据有：

进行傅里叶变换，得到其频域形式

其中x_mp,l(f)为一个频域数据，是

的傅里叶变换形式。s_pk,l(f)和n_mp,l(f)分别是第l个频域快拍下的源信号和噪声信号,频域快拍数为L。

将所有M个阵元的L个频域接收信号排成一个矩阵，得到有关频率f_p的所有数据的阵列频域接收模型，即

其中方向矩阵A(f_p,θ)，信源矩阵S(f_p)和噪声矩阵N(f_p)可以表示为：

其中s_l(f_p)表示S(f_p)的第l列向量，n_l(f_p)表示N(f_p)的第l列向量。

S2、根据频域数据接收模型对不同频点的接收数据构造聚焦矩阵，进行频率聚焦，同时完成接收模型的重构，实现第一次虚拟阵元扩展。

通过将各个频点的阵列流型通过“聚焦”变换到参考频率点f₀上，最后对参考频率点的数据采用窄带处理方式进行DOA估计。对应每个频率构造一个聚焦变换矩阵T(f_p)，使得不同频率点上的阵列流型矩阵聚焦到聚焦频点f₀上，即：

T(f_p)A(f_p,θ)＝A(f₀,θ) (13)

这里可以采用现有的双边相关变换(TCT)法求出聚焦矩阵T(f_p)。接下来，使用聚焦矩阵对接收信号进行变换，将变换之后接收信号表示为Y(f_p)，则有：

Y(f_p)＝T(f_p)X(f_p) (14)

从上式可以看出，通过“聚焦”使得不同频率点信号拥有了同样的阵列流型矩阵，即拥有同一个信号子空间。接着求出对应频点聚焦变换之后的信号协方差矩阵，并对P个频率分量聚焦变换后信号协方差矩阵进行平均求和可得到最终的聚焦协方差矩阵R_Y

其中，

为了更好地利用多个频率数据的信息，接下来进行接收协方差矩阵的扩展重构。重构出的接收协方差矩阵R_F表达为：

其中，协方差矩阵R_F的第i·M+1～(i+1)·M行，j·M+1～(j+1)·M列(i，j＝0,1,…,P-1)的部分表达式为：

则协方差矩阵R_F最终可以重写为：

R_F＝B·R_FS·B^H+R_FN (20)

其中B＝[A(f₀,θ),…,A(f_P-1,θ)]^T∈C^MP×K为重构后的方向矩阵，R_FS∈C^K×K为重构后的源信号矩阵，R_FN∈C^MP×MP为重构后的噪声矩阵。

重构方向矩阵B同时可以表示为：

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

其中b(θ_k)表示重构方向向量，可以表达成：

由于频率数P≥2，本实施例以P＝2为例，构造类似二级嵌套阵列的虚拟非均匀阵列。设定第p路频率f_p存在比例关系：

f₁＝M·f₀ (23)

那么方向向量b(θ_k)可以重写为：

其中z_w表示一次扩展后的第w,w＝1,2,…,MP虚拟阵元的位置，

由上式可见，物理阵列通过重构扩展成一个非均匀的虚拟阵列。

S3、重构后的频域数据协方差矩阵R_F为MP×MP维的协方差矩阵，第u行v列的元素可以表示为：

其中

是重构数据的第k个源信号功率，

是重构数据的噪声功率，δ_u,v为Kronecker delta函数。

由式(28)可见，一次扩展虚拟阵列的协方差矩阵中的元素可以视为是以差协同阵列作为阵元位置的阵列接收的数据，因此可以将一次扩展虚拟阵列的协方差矩阵R_F向量化，得到：

其中vec为向量化符号，⊙表示Khatri-Rao积，g为源信号功率向量，噪声向量

向量化方向矩阵B^*⊙B表示为：

其中

表示Kronecker积。

向量化后的二次虚拟阵列中阵元数远多于物理阵列中的阵元数，通过将一次扩展虚拟阵列接收数据的协方差矩阵向量化的方式，转化为二次扩展虚拟阵列下的等效单快拍接收数据，阵元数再次大大增加，故而达到提高自由度的目的。

S4、经过协方差矩阵矢量化得到r后，向量化方向矩阵B^*⊙B中只有2(M-1)M+1行是互不相同的，即互不相同行数与二次扩展虚拟阵列的阵元数相同。此时的入射信号g是K个源信号的二阶统计值，即功率值，因此新的信号向量g相当于K个相干信号入射。为了使用空间平滑技术解相干，对B^*⊙B进行去冗余和重排处理，使处理后的方向矩阵与虚拟均匀线阵的方向矩阵对应。假设经过去冗余和重排处理后的等效阵列信号接收向量为r₁，表示为：

其中方向矩阵A₂∈C^{(2(M-1)M+1)×K}为B^*⊙B进行去冗余和重排处理后的方向矩阵，

为新的噪声向量。

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

S5、对得到的接收向量r₁进行子阵划分处理。这里将二次扩展虚拟均匀阵列划分为(M-1)M+1个子阵，每个子阵包含(M-1)M+1个阵元，第h,h＝1,2,…,(M-1)M+1子阵的接收数据r_1h等于阵列信号接收向量r₁的第h～(h+(M-1)M+1)列。

令

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

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

S6、基于R_s采用一维MUSIC算法估计全部K个信源目标的来波方向角度。将所获得的R_s实现一次特征分解，获得(M-1)M+1个特征值，把特征值依照大小实现顺序排列，提取最小的(M-1)M+1-K个特征向量构建出U_n。再构造方向向量a_s(θ)，可表达为：

把a_s(θ)代入得到MUSIC谱函数：

根据MUSIC谱函数，使θ从-90到90°变化，通过寻求K个极大值，对应的角度即为信源目标的波达角估计值。

实施例二

本实例公开了一种基于多频率组合信号的高精度一维波达方向估计装置，用于作为上述实施例中一维DOA估计方法的实施装置，包括数据处理与控制模块、发射模块、接收模块、输出模块和电源模块，如图1和图2所示。

数据处理与控制模块由一对A/D、D/A转换器和一个处理器组成，是整个装置的核心部分，其它所有模块都与它直接相连。它可以控制发射模块，使发射模块发射我们指定的信号；同时能够对接收模块传过来的信号进行处理，通过本发明的算法计算出波达方向角，然后将结果传输至输出模块。

接收模块包括以M个物理阵元均匀间距摆放的超声波探头阵列和P组不同频率的窄带滤波器，其中阵列的各个阵元均匀排布于坐标系x轴上保持固定，每个阵元的输出再作为输入分别经过P个窄带滤波器，得到M·P路输出；

发射模块由一个阻抗匹配电路和一个超声波发射探头组成，通过D/A转换器与处理器相连，能够根据处理器发出的指令发射指定的携带P组不同频率的混合信号。

输出模块由一个USB接口和一个显示器组成，并且与数据处理与控制模块和电源模块相连。它能够提供人机交互，将数据处理与控制模块中处理好的数据通过USB接口输出到外部装置或者在显示器上显示出来。

电源模块由一个电源组成，并且与数据处理与控制模块、发射模块、接收模块和输出模块相连。它能够为这些模块供电。

本装置的主要工作流程如下：在实测过程中根据我们想要发射的信号参数，通过数据处理与控制模块输入对应的参数，使处理器产生相应的数字信号，然后通过D/A转换后传给发射模块，超声波发射探头就能产生我们需要的信号并进行发射。均匀线阵1和均匀线阵2之间的夹角值δ可以通过数据处理与控制模块进行设定，处理器发送特定的脉冲信号到步进电机驱动电路，然后驱动步进电机转动至我们需要的角度。接收模块中的接收阵列收到从目标声源反射回来的信号后将其通过A/D转换成数字信号后发送给处理器，然后处理器根据本发明提供的算法计算出结果。最后数据处理与控制模块将计算结果传给输出模块，输出模块将结果通过USB接口传给外部设备或者通过显示器显示出来。电源模块为所有其它模块供电。

实施例三

本实施例公开了一种基于多频率组合信号的高精度一维波达方向估计装置，所述的估计装置包括数据处理与控制模块、发射模块、接收模块、输出模块和电源模块，如图1和图2所示。

数据处理与控制模块可以用DSP芯片实现(如：TI公司TMS320VC5509A型号的DSP芯片)，此DSP芯片可实现A/D转换和D/A转换的功能，并能够实现三维均匀线阵的旋转算子和最终波达方向的计算；

此外接收模块使用1个固定均匀直线阵列和P组中心频率分别为{f₁,f₂,…,f_P}的窄带滤波器，其中每个阵列包括多个超声接收探头，并且数量相同，按图4所示组装；发射模块使用一个超声波发射探头；输出模块使用一个USB接口和一个LCD显示屏。图1即为本实施例估计装置的硬件结构模块图。

本估计装置主要工作步骤具体如下：

步骤T1、按图2连接好具体装置，其中接收模块中的每个均匀线阵中的阵元个数M统一定为5。利用数据处理与控制模块发送指令，控制超声发射探头发射超声信号

频率数P＝2，频率分别取f₀＝3kHz,f₁＝15kHz；声速取c＝1500m/s，可以求出最大半波长为25cm。在本发明中，需满足任意两相邻线阵之间的距离必须小于等于最大半波长的条件，所以设置两个均匀线阵的平均间距取25cm，即第一个阵元和最后一个阵元相隔100cm。在水下放置K＝2个目标声源，信源目标信号入射的波达方向角分别为(40，°70°)。

步骤T2、对超声接收探头阵接收到的源信号进行采样；阵列接收到的信号为x₁₀(t),x₂₀(t),x₃₀(t),x₄₀(t),x₅₀(t),x₁₁(t),x₂₁(t),x₃₁(t),x₄₁(t),x₅₁(t)，共时域采样接收1024次，频域快拍数设置为16，并将接收到的信号传递给数据处理与控制模块进行分析处理。

步骤T3、信号在数据处理与控制模块中的分析处理步骤具体如下：

T31、根据阵列特点，得到多频率信号时域接收数据模型和频域接收数据模型。

T32、根据频域的数据接收模型对不同频点的接收数据构造聚焦矩阵，进行频率聚焦，同时完成接收模型的重构，实现第一次虚拟阵元扩展。物理阵列进行一次扩展得到10阵元的虚拟非均匀阵列。

T33、对重构后的频域接收数据协方差矩阵R_F进行向量化得到接收向量r，虚拟非均匀阵列二次扩展得到41个虚拟均匀阵元。

T34、对向量r进行冗余和重排处理得到等效阵列信号接收向量r₁，实现向量化方向矩阵B^*⊙B与虚拟均匀阵列的方向矩阵的一一对应。

T35、对得到的接收向量r₁进行虚拟子阵划分和平均处理，按式(29)得到空间平滑后的协方差矩阵R_s。

T36、基于R_s，按照一维MUSIC谱函数即式(31)进行谱峰搜索，得到全部K个信源目标的波达方向角度的估计值。

步骤T4、根据一维DOA估计方法，估计出的二维波达方向角(40.15°,70.23°)，对目标估计达到了预期精度，说明估计结果正确，本发明方法及装置可行。

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

Claims (5)

1.一种利用多频信号二次虚拟扩展的一维DOA估计方法，其特征在于，所述的一维DOA估计方法包括以下步骤：

S1、根据用于一维DOA估计的一维均匀线阵结构特点构造多频率信号时域数据接收模型，同时利用傅里叶变换得到多频率信号频域数据接收模型；

S2、根据频域数据接收模型对不同频点的接收数据构造聚焦矩阵，进行频率聚焦，同时完成频域数据接收模型的重构，实现第一次虚拟阵元扩展；

S3、对重构后的频域数据协方差矩阵R_F进行向量化得到向量r，实现基于一次扩展虚拟非均匀阵列的二次虚拟阵元扩展；

S4、对向量r进行冗余和重排处理得到等效阵列信号接收向量r₁，实现以B^*⊙B为整体作为向量化方向矩阵与虚拟均匀阵列的方向矩阵的一一对应关系，同时便于下一步的空间平滑处理，其中，⊙表示Khatri-Rao积；

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

S6、基于协方差矩阵R_s采用一维MUSIC算法估计全部K个信源目标的来波方向角度；

其中，所述的步骤S1中得到时域数据接收模型和频域数据接收模型的过程如下：

假设一维均匀线阵上分布有M个阵元，阵元间距为d，接收信号为由P个不同频率的窄带信号组成的发射信号经K个信源目标反射到接收阵列的信号，窄带信号满足窄带条件，即当信号延迟远小于带宽倒数时，延迟作用相当于使基带信号产生一个相移；

当有K个远场信源目标的源信号入射到一维均匀线阵，在接收端通过接收阵元和对应窄带滤波器，经H个时域快拍数采样得到时域接收数据矩阵X，该数据矩阵X由P,P≥2个频率信号的接收数据矩阵叠加，表示为：

其中X_p为第p,p＝0,1,…,P-1路频率信号的接收数据表达为：

X_p＝[x_p(1),x_p(2),…,x_p(H)] (2)

x_p(t)＝[x_1p(t),x_2p(t),…,x_Mp(t)]^T,t＝1,2,…,H (3)

接收端第m,m＝1,2,…,M阵元的第p路窄带滤波器得到的信号x_mp(t)如下：

其中s_pk(t)为第p路频率发射信号经第k个目标在t时刻到达阵列的信号，n_mp(t)为第m阵元第p路窄带滤波器的在t时刻的噪声信号，其中，

a_mp(θ_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)

对时域接收信号x_mp(t)的H个时域快拍数据分成L段，假设H能被L整除，则第l段时域数据有：

进行傅里叶变换，得到其频域形式

其中x_mp,l(f)为一个频域数据，是

的傅里叶变换形式，s_pk,l(f)和n_mp,l(f)分别是第l个频域快拍下的源信号和噪声信号,频域快拍数为L；

将所有M个阵元的L个频域接收信号排成一个矩阵，得到有关频率f_p的所有数据的频域数据接收模型X(f_p)，即

其中方向矩阵A(f_p,θ)，信源矩阵S(f_p)和噪声矩阵N(f_p)可以表示为：

其中s_l(f_p)表示S(f_p)的第l列向量，n_l(f_p)表示N(f_p)的第l列向量；

其中，所述的步骤S2中对不同频点的接收数据进行频率聚焦，同时完成频域数据接收模型重构的过程如下：

将各个频点的阵列流型通过“聚焦”变换到参考频率点f₀上，对参考频率点的数据采用窄带处理方式进行DOA估计，对应每个频率构造一个聚焦变换矩阵T(f_p)，使得不同频率点上的阵列流型矩阵聚焦到聚焦频点f₀上，即：

T(f_p)A(f_p,θ)＝A(f₀,θ)(13)

采用双边相关变换法求出聚焦矩阵T(f_p)，然后使用聚焦矩阵对接收信号进行变换，将变换之后接收信号表示为Y(f_p)，则有：

Y(f_p)＝T(f_p)X(f_p)(14)

通过“聚焦”使得不同频率点信号拥有同样的阵列流型矩阵，即拥有同一个信号子空间，接着求出对应频点聚焦变换之后的信号协方差矩阵，并对P个频率分量聚焦变换后信号协方差矩阵进行平均求和得到最终的聚焦协方差矩阵R_Y：

其中R_s(f_p)和R_n(f_p)表示为：

对接收协方差矩阵的扩展重构，重构出的接收协方差矩阵R_F表达为：

其中，协方差矩阵R_F的第i·M+1～(i+1)·M,i＝0,1,…,P-1行，j·M+1～(j+1)·M,j＝0,1,…,P-1列的部分表达式为：

则协方差矩阵R_F重写为：

R_F＝B·R_FS·B^H+R_FN(20)

其中B＝[A(f₀,θ),…,A(f_P-1,θ)]^T∈C^MP×K为重构后的方向矩阵，R_FS∈C^K×K为重构后的源信号矩阵，R_FN∈C^MP×MP为重构后的噪声矩阵；

重构方向矩阵B同时表示为：

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

其中b(θ_k)表示重构方向向量，表达如下：

设定第p路频率f_p存在比例关系：

f₁＝M·f₀(23)

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

其中，z_w表示一次扩展后的第w,w＝1,2,…,MP虚拟阵元的位置，

2.根据权利要求1所述的利用多频信号二次虚拟扩展的一维DOA估计方法，其特征在于，所述的步骤S3中对协方差矩阵R_F进行向量化得到向量r的过程如下：

协方差矩阵R_F为MP×MP维，第u行v列的元素表示为：

其中

是重构数据的第k个源信号功率，

是重构数据的噪声功率，δ_u,v为Kroneckerdelta函数；

将一次扩展虚拟阵列的协方差矩阵R_F向量化，得到：

其中，vec为向量化符号，⊙表示Khatri-Rao积，g为源信号功率向量，噪声向量

向量化方向矩阵B^*⊙B表示为：

其中

表示Kronecker积。

3.根据权利要求2所述的利用多频信号二次虚拟扩展的一维DOA估计方法，其特征在于，所述的步骤S4中对向量r进行冗余和重排处理得到等效阵列信号接收向量r₁的过程如下：

对向量化方向矩阵B^*⊙B进行去冗余和重排处理，使处理后的方向矩阵与虚拟均匀线阵的方向矩阵对应；假设经过去冗余和重排处理后的等效阵列信号接收向量为r₁，表示为：

其中方向矩阵

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

为新的噪声向量；

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

4.根据权利要求3所述的利用多频信号二次虚拟扩展的一维DOA估计方法，其特征在于，所述的步骤S5中对接收向量r₁进行空间平滑处理后，得到协方差矩阵R_s的过程如下：

对接收向量r₁进行子阵划分处理； 这里将二次扩展虚拟均匀阵列划分为(M-1)M+1个子阵，每个子阵包含(M-1)M+1个阵元，第h子阵的接收数据r_1h等于阵列信号接收向量r₁的第h～(h+(M-1)M+1)列，其中，h＝1,2,…,(M-1)M+1，令

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

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

5.根据权利要求4所述的利用多频信号二次虚拟扩展的一维DOA估计方法，其特征在于，所述的步骤S6中根据协方差矩阵R_s采用一维MUSIC算法估计K个信源目标的来波方向角度的过程如下：

将所获得的协方差矩阵R_s实现一次特征分解，获得(M-1)M+1个特征值，把特征值依照大小实现顺序排列，提取最小的(M-1)M+1-K个特征向量构建出U_n，再构造方向向量a_s(θ)，表达为：

把a_s(θ)代入得到MUSIC谱函数：

根据MUSIC谱函数，使θ从-90到90°变化，通过寻求K个极大值，对应的角度即为信源目标的波达角估计值。

Images