CN109188345B - 基于去预延迟空时结构的相干信号源doa估计方法 - Google Patents

Abstract

本发明公开了一种基于去预延迟空时结构的相干信号源DOA估计方法，属于自适应阵列信号处理技术领域。本发明针对现有的宽带相干信号算法CSM只应用于相关情况，并且计算量太大。而目前的空时结构宽带DOA不能解相干，本发明公开了一种可以同时解决相干源和频域宽带测向计算量大的DOA估计方法。在基于去预延迟空时结构进行信号重组后，获取待搜索角对应的约束矩阵、约束向量，再采用迭代处理获取搜索角度的功率谱；最后，由搜索角度的功率谱的峰值位置得到信号源的DOA估计结果。本发明提出的方法在信号数据量较少和低信噪比的条件下，分辨率比前者有所提高。

Description

基于去预延迟空时结构的相干信号源DOA估计方法

技术领域

本发明涉及自适应阵列信号处理技术领域，具体是涉及基于去预延迟空时结构的一种相干信号源DOA(Direction Of Arrival)估计方法。

背景技术

现有的相关宽带信号的DOA估计算法主要是相关信号子空间算法(CSM，CoherentSignal-Subspace Method)，其基本思想就是通过傅氏变换将宽带阵列数据分解为若干个窄带数据，再通过聚焦矩阵把不同频率上的协方差矩阵聚焦到某一个参考频率上，然后利用窄带估计方法估计DOA，具体可参考文献“Focusing matrices for coherentsignal-subspaceprocessing；H.Hung,M.Kaveh；IEEE Trans.SP,June 1988,36(8),1272-1281”。此方法是基于空频结构的，基于空频结构的算法也是DOA算法中很常见的。然而在实际的应用中，求解聚焦矩阵的计算量是非常大的，并且当数据量大时，此算法测向结果误差也会增加。此外算法要求的输入信号源是基于信号之间的时延生成的，即考虑的是相关性源，对于相干信号，并不能有很好的效果。所以在完全相干并且信号采样数据较少的情况下，此方法将不再适用。具体分析如下：

对于有两个宽带相关源s₁(t)和s₂(t)，并且s₂(t)＝s₁(t-t₀)，分别来自两个不同的方向θ₁和θ₂，令

则相关函数矩阵为：

其中，R₁(τ)是s₁(t)的自相关函数，E[·]表示期望函数。对上式作傅立叶变换，得到：

当t₀≠0时，即s₂(t)不完全等于s₁(t)，对上式进行积分，得到

大多数情况下，(3)式为非奇异矩阵。可以看出，若将所有频率成分的信号功率谱密度矩阵作平均，就可以消除相关源相关矩阵的奇异性。

若是相干信号，(1)式中互相关系数不再是R₁(τ-t₀)，而是变为了只差一个系数的kR₁(τ)，此时的(3)式无法再利用上面的方式消除相关源相关矩阵的奇异性。

考虑一等间距排列的M元均匀宽带线性阵列，阵列单元间距d为信号中心频率对应的半波长。D个有相同带宽B和中心频率f₀的宽带源信号。噪声为相互独立的零均值高斯白噪声，且与信号不相关，方差为σ_n ²，并与源信号具有相同的带宽。

对持续时间T_d上的阵列输出数据向量x(t)经过快速离散傅立叶变换(FFT)分解为Q个互不重叠的窄带部分。对于K个持续时间T_d，可得到在f_q频带上得到的数据向量为：

X_k(f_q)＝A(f_q)·s_k(f_q)+N_k(f_q)q＝1,…,Q k＝1,…,K (4)

其中f_q表示第q个子带的中心频率，X_k(f_q)为M×1维子带接收信号向量，s_k(f_q)表示信号在第q个频带的谱分量，N_k(f_q)为M×1维的向量，表示噪声在f_q处的谱分量。

A(f_q)＝[a(θ₁,f_q),a(θ₂,f_q),…,a(θ_D,f_q)]表示f_q处信号的方向矩阵。θ_i表示第i个信号的入射角度，i＝1,…,D，其中

e表示自然底数，j表示虚数单位。

CSM的中心问题是构造Q个“聚焦”矩阵T(f_q)(i＝1,…,Q)，将不同频率段的方向阵聚焦到同一参考频率段f_c下的矩阵，即

T(f_q)A(f_q)＝A(f_c) (5)

聚焦后的阵列输出向量为

聚焦后的数据为

Y(f_q)＝T(f_q)X(f_q) (7)

协方差矩阵为

式中

令：

如果取w_j＝1，其中w_j表示与信噪比成比例的归一化权值。将(9)式与(10)式代入(8)式中，得

R_y＝A(f_c)R_SA^H(f_c)+R_N (11)

在噪声为白噪声的设置下，R_n(f_j)＝σ_nI。那么

由于T(f_q)T^H(f_q)≠I，所以在该方法聚焦后的协方差矩阵中包含了有色噪声。在这种情况下，信号子空间会泄露到噪声子空间。如果直接利用式(8)作为协方差矩阵，对其采用MUSIC类的窄带方法进行DOA估计。那么估计结果将大受影响。所以在利用窄带子空间方法进行DOA估计之前，需要对式(8)进行下述预白化处理。

将式(12)进行Cholesky分解，得到

R_N＝R_N ^H/2R_N ^1/2 (13)

将

和

分别右乘和左乘R_y，有下式

对矩阵R_y′进行特征值分解，求得噪声子空间E_N，令

E′_N＝R_N ^-1/2E_N (15)

将E′_N代替E_N代入

其中，θ代表角度扫描范围，函数P(θ)的K个峰值位置就是信号波达方向θ_D的估计。

发明内容

本发明的发明目的在于：针对现有的宽带相干信号算法CSM只应用于相关情况，并且计算量太大。而目前的空时结构宽带DOA不能解相干，本发明公开了一种可以同时解决相干源和频域宽带测向计算量大的方法。

本发明的基于去预延迟空时结构的相干信号源DOA估计方法，包括下列步骤：

步骤S1：对阵列接收信号x₁(n),x₂(n),…x_M(n)进行抽样重组，得到重组后信号的数据矩阵X：

基于预设的抽头数Q，得到阵列各阵元的接收信号x_m(n)在依次经过Q个抽头后的数据分别为：x_m0(n)＝x_m(n)，x_m1(n)＝x_m(n+1)，…，x_mQ(n)＝x_m(n+Q)，其中m＝1,…,M,M表示接收阵元数；

令x^T(n)＝[x₁₀(n),x₂₀(n)…,x_M0(n),x₁₁(n),x₂₁(n)…,x_M1(n),…,x_1Q(n),x_2Q(n)…,x_MQ(n)]，得到数据矩阵X＝[x(1)x(2)x(3)…x(N-J)]，其中x(n)中的n＝1,…,N-Q，符号“(·)^T”表示矩阵转置，N表示阵列接收信号的快拍数；

步骤S2：对角度搜索范围进行均匀离散处理，得到多个离散点，用P表示角度搜索范围的离散点数，θ_p表示各离散点对应的搜索角度，其中p＝1,2,…,P；

对接收信号的期望带宽进行均匀离散处理，得到多个离散点，用D表示期望带宽的离散点数，f_d表示各离散点对应的频率；

定义各搜索角度θ_p的空时结构的约束矩阵C_d(θ_p)＝[c(f₁,θ_p)c(f₂,θ_p)…c(f_D,θ_p)]_MQ×D，其中c(f_d,θ_p)表示搜索角度θ_p在频率f_d上的约束向量，为

符号

表示Kronecker积，并且：

其中e表示自然底数，j表示虚数单位，T_s表示采样间隔，τ_m(θ_p)表示对应搜索角度θ_p，信号到达阵元m与到达参考阵元的时间延迟；

定义约束向量

步骤S3：通过迭代处理计算每个搜索角度θ_p的功率谱P(θ_p)：

步骤S301：初始化每个搜索角度θ_p的最优化权值向量W(θ_p)为：

以及初始化接收信号的自相关矩阵R＝I；

其中I表示单位矩阵，符号“(·)^H”表示矩阵共轭；

步骤S302：计算各搜索角度θ_p的功率谱P(θ_p)为

并更新自相关矩阵R为：R＝R+P(θ_p)C_d(θ_p)C_d(θ_p)^H；

步骤S303：基于当前的功率谱对自相关矩阵R，对W(θ_p)和P(θ_p)进行迭代更新：

基于当前R，更新W(θ_p)为：

更新完所有W(θ_p)后，再令R＝I；

基于当前W(θ_p)，更新P(θ_p)为：

其中

基于当前P(θ_p)，更新R为：R＝R+P(θ_p)C_d(θ_p)C_d(θ_p)^H；

步骤S304：判断是否满足迭代收敛条件，若是，则输出当前功率谱P(θ_p)；否则继续执行步骤S303；

其中迭代收敛条件为：迭代次数达到上限，或者相邻两次迭代更新得到的功率谱的相对改变量小于或等于预设阈值；

步骤S4：根据功率谱P(θ_p)的峰值位置得到信号源的DOA估计结果，即对角度搜索范围按固定步长进行遍历，若当前搜索角度的功率谱出现峰值，则将该搜索角度作为DOA估计结果。

综上所述，由于采用了上述技术方案，本发明的有益效果是：本发明在低信噪比和信号数据较少的条件下也能成功进行DOA估计并且误差比现有方式小，即本发明在信号数据量较少和低信噪比的条件下，有效的提升了DOA估计的分辨率。

附图说明

图1为现有的Frost空时处理结构的结构示意图；

图2为同等条件下，信号源为全相干情况下的本发明方法的DOA估计(基于空时去预延迟结构IAA方法)结果与传统基于空频结构的CSM算法DOA估计对比图；

图3为同等条件下，信号源为相关情况下的本发明方法的DOA估计结果与传统基于空频结构的CSM算法DOA估计对比图；

图4为同等条件下，信号源为相关情况下的本发明方法的DOA估计与基于空频结构的CSM算法DOA估计的均方根误差(RMSE)随输入SNR的变化对比图；

图5为同等条件下，信号源为相关情况下的本发明方法的DOA估计与基于空频结构的CSM算法DOA估计的均方根误差(RMSE)随快拍数N的变化对比图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面结合实施方式和附图，对本发明作进一步地详细描述。

本发明基于空时Frost处理结构处理宽带信号，然后利用迭代自适应思想(IAA,IterativeAdaptive Approach)选择初始值对其迭代从而到达解相干的目的。

本发明对宽带信号处理部分运用的是Frost空时处理结构，如图1所示。其核心思想是通过在每个阵元后面增加Q个延迟器或时域抽头，这些延迟期或抽头组成了有限长单位冲激响应(FIR)滤波器，形成了一组与频率有关的响应，用来补偿因为频率的变化形成的相位差。当阵元接收到的信号为宽带信号时，根据孔径渡越效应，会发生不同频率对应的相移值不同的现象，因而传统的移相器方法无法形成正确的波束指向；而若是采用Frost结构，即用时间延迟器来替代移相单元，则可通过补偿回波信号到达不同阵元之间的时间差，从而将主瓣对准期望方向。因此，经预延迟处理的空时处理结构能很好的波束形成，却不能用于信号源DOA估计。

当去掉空时处理结构中的预延迟模块后，用x_k,m(n)表示入射角度为θ_k的接收信号，则可以得到其输出为：

其中，M表示接收阵元数，ω_mq表示自适应权值，m＝1,…,M,q＝0,…,Q，Q表示去预延迟模块的空时结构的时域抽头数；τ_m(θ_k)表示第k个信号到达阵元m时，相对于该信号到达参考阵元的时间延迟，T_s表示采样间隔。

式(17)的频域形式如下：

可以得到信号x_k,m(n)的空-频域转换函数H_d(f,θ_k)：

其中

f_q＝s(f,θ_k)^TW_q。与传统空时处理结构相似，f_q可以看作为集成所有分支的TDL(Tapped Delay Line，抽头延迟)滤波器的集合TDL滤波器的权重。为了输出信号信号不失真，应该有H(f)＝e^{-j2πf(Q-1)/2}，频率f必须包含在信号的期望带宽内，即f_d∈[f_l f_h],d＝1,2,…,D，其中f_l和f_h分别代表了期望信号的最低频率和最高频率，D表示选择频率的数量。

为了求解权值W_q，采用线性约束最小方差(LCMV)准则，线性约束最小方差准则是通过选择权值W，在一定的线性约束条件下，使得输出的方差或者功率最小，输出功率或者方差为

所以线性约束最小方差(LCMV)准则就是求解如下最优化问题：

其中，R表示信号源经过时域抽头后数据的自相关矩阵，C_d为约束矩阵，f为约束向量。利用拉格朗日乘子法，可以得到基于LCMV(基于线性约束最小方差)准则的最优化权向量为：

空时处理结构中的时域抽头或延迟器可以通过对信号源数据进行抽样重组来达到同等效果，已知阵列输入信号为x₁(n),x₂(n),…x_M(n)，则阵元1的接收信号经过各抽头后的数据为x₁₀(n)＝x₁(n),x₁₁(n)＝x₁(n+1),…,x_1Q(n)＝x₁(n+Q)，相应的权为w_1,0,w_1,2,…,w_1,Q。阵元m后接收信号经过各抽头后的数据为x_m0(n)＝x_m(n),x_m1(n)＝x_m(n+1),…,x_mQ(n)＝x_m(n+Q)，相应的权为w_m,1,w_m,2,…,w_m,Q。令

则

X＝[x(1) x(2) x(3) … x(N-Q)] (24)

所得数据矩阵X即可等效于信号源数据通过时域抽头后的结果，自相关矩阵

用数据矩阵X的统计平均进行估计，即

根据式(22)，(20)可以得到权向量，进而得到最终的输出功率谱，而输出功率谱对应着角度估计谱。

此时，如果让式(26)中的自相关函数R直接等于信号自相关的估计量

能对非相干信号有很好的效果，但若是宽带信号相干，则该方法会失效，若是用该方法来进行DOA估计则无法得到准确的结果。

针对宽带相干信号的DOA估计问题，本发明采用迭代处理去相干的空时去预研迟处理的结构的DOA估计，具体为：

定义Θ表示角度搜索范围，对Θ等间隔的取出P个点，得到P个搜索角度θ_p，p＝1,2,…,P，定义空时结构的约束矩阵为：

C_d(θ_p)＝[c(f₁,θ_p)c(f₂,θ_p)…c(f_D,θ_p)]_MQ×D (28)

其中c(f_d,θ_p)为选择频率f_d上的约束向量。定义c(f_d,θ_p)为：

其中

表示Kronecker积，并且：

约束向量为：

(1)初始化每个搜索角度θ_p的最优化权值向量W(θ_p)：

以及初始化接收信号的自相关矩阵R＝I，其中I表示单位矩阵；

(2)计算每个搜索角度θ_p的功率谱P(θ_p)：

并更新自相关矩阵R为：

R＝R+P(θ_p)C_d(θ_p)C_d(θ_p)^H,p＝1,...,P (34)

(3)基于当前的功率谱对自相关矩阵R，对最优化权值向量W、功率谱P进行迭代更新：

基于当前R，更新W(θ_p)为：

更新完P个W(θ_p)后，再令R＝I；

基于当前W(θ_p)，更新P(θ_p)为：

基于当前P(θ_p)，更新R为：

R＝R+P(θ_p)C_d(θ_p)C_d(θ_p)^H,p＝1,...,P (37)

(4)判断是否收敛，即当相邻两次迭代更新得到的功率谱的相对改变量小于或等于预设阈值时，则满足迭代收敛；或者当迭代次数达到预设上限时满足迭代收敛；当满足迭代收敛时，输出当前功率谱P(θ_p)；否则继续执行步骤(3)。

根据上述方式得到最终收敛后的每一个搜索角度的P(θ_p)，则可以得到峰值对应的角度则是信源入射的角度，即由P(θ_p)的峰值位置得到信号源的DOA估计结果。

为了进一步验证本发明方法与现有方法的性能，对基于空频结构的CSM算法DOA估计方法以及本发明提出的DOA估计方法(基于空时去预延迟结构IAA方法)进行了仿真，仿真结果证明了本发明提出的方法在信号数据量较少和低信噪比的条件下，分辨率比前者有所提高。

仿真试验1：10元均匀线阵，信号源包含两个完全相干的宽带信号，信号带宽均为400MHz，入射方向分别为10°、-12°，中心频率均为2200MHz，对信号进行带通采样，采样频率为800MHz，阵元间距为二分之一的最小波长；对于去预延迟处理的空时结构，设置时域抽头数为16，频率选取点数为14；对于空频结构，设置FFT点数为400；其中信号源方向可能存在的区间为Θ＝[-20°,20°]，仿真结果如图2所示。

仿真试验2：10元均匀线阵，信号源包含两个相关性不太强的相关的宽带信号，信号带宽均为400MHz，入射方向分别为10°、-12°，中心频率均为2200MHz，对信号进行带通采样，采样频率为800MHz，阵元间距为二分之一的最小波长；对于去预延迟处理的空时结构，设置时域抽头数为16，频率选取点数为14；对于空频结构，设置FFT点数为400；其中信号源方向可能存在的区间为Θ＝[-20°,20°]，仿真结果如图3所示。

仿真试验3：10元均匀线阵，信号源包含两个相关性不太强的相关的宽带信号，信号带宽均为400MHz，入射方向分别为15°、-10°，中心频率均为2200MHz，对信号进行带通采样，采样频率为800MHz，阵元间距为二分之一的最小波长；对于去预延迟处理的空时结构，设置时域抽头数为16，频率选取点数为14；对于空频结构，设置FFT点数为400；其中信号源方向可能存在的区间为Θ＝[-20°,20°]，信号快拍数为3000，输入信噪比变化范围是-10dB到30dB，进行100次蒙特卡罗实验得到两种算法的信号源DOA估计结果的均方根误差随输入信噪比的变化而变化的对比图，仿真结果如图4所示。

仿真试验4：10元均匀线阵，信号源包含两个相关性不太强的相关的宽带信号，信号带宽均为400MHz，入射方向分别为15°、-10°，中心频率均为2200MHz，对信号进行带通采样，采样频率为800MHz，阵元间距为二分之一的最小波长；对于去预延迟处理的空时结构，设置时域抽头数为16，频率选取点数为14；对于空频结构，设置FFT点数为400；其中信号源方向可能存在的区间为Θ＝[-20°,20°]，输入信噪比是5dB，信号快拍数变化范围为500到4000，进行100次蒙特卡罗实验得到两种算法的信号源DOA估计结果的均方根误差随信号快拍数的变化而变化的对比图，仿真结果如图5所示。

由图4和图5可以看出，本发明在低信噪比和信号数据较少的条件下也能成功进行DOA估计并且误差比传统方法小。

以上所述，仅为本发明的具体实施方式，本说明书中所公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的替代特征加以替换；所公开的所有特征、或所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以任何方式组合。

Claims (1)

1.基于去预延迟空时结构的相干信号源DOA估计方法，其特征在于，包括下列步骤：

步骤S1：对阵列接收信号x₁(n),x₂(n),…x_M(n)进行抽样重组，得到重组后信号的数据矩阵X：

基于预设的抽头数Q，得到阵列各阵元的接收信号x_m(n)在依次经过Q个抽头后的数据分别为：x_m0(n)＝x_m(n)，x_m1(n)＝x_m(n+1)，…，x_mQ(n)＝x_m(n+Q)，其中m＝1,…,M,M表示接收阵元数；

令x^T(n)＝[x₁₀(n),x₂₀(n)…,x_M0(n),x₁₁(n),x₂₁(n)…,x_M1(n),…,x_1Q(n),x_2Q(n)…,x_MQ(n)]，得到数据矩阵X＝[x(1) x(2) x(3) … x(N-J)]，其中x(n)中的n＝1,…,N-Q，符号“(·)^T”表示矩阵转置，N表示阵列接收信号的快拍数；

步骤S2：对角度搜索范围进行均匀离散处理，得到多个离散点，用P表示角度搜索范围的离散点数，θ_p表示各离散点对应的搜索角度，其中p＝1,2,…,P；

对接收信号的期望带宽进行均匀离散处理，得到多个离散点，用D表示期望带宽的离散点数，f_d表示各离散点对应的频率；

定义各搜索角度θ_p的空时结构的约束矩阵C_d(θ_p)＝[c(f₁,θ_p) c(f₂,θ_p) … c(f_D,θ_p)]_MQ×D，其中c(f_d,θ_p)表示搜索角度θ_p在频率f_d上的约束向量，为

符号

表示Kronecker积，并且：

其中e表示自然底数，j表示虚数单位，T_s表示采样间隔，τ_m(θ_p)表示对应搜索角度θ_p，信号到达阵元m与到达参考阵元的时间延迟；

定义约束向量

步骤S3：通过迭代处理计算每个搜索角度θ_p的功率谱P(θ_p)：

步骤S301：初始化每个搜索角度θ_p的最优化权值向量W(θ_p)为：

以及初始化接收信号的自相关矩阵R＝I；

其中I表示单位矩阵，符号“(·)^H”表示矩阵共轭；

步骤S302：计算各搜索角度θ_p的功率谱P(θ_p)为

并更新自相关矩阵R为：R＝R+P(θ_p)C_d(θ_p)C_d(θ_p)^H；

步骤S303：基于当前的功率谱对自相关矩阵R，对W(θ_p)和P(θ_p)进行迭代更新：

基于当前R，更新W(θ_p)为：

更新完所有W(θ_p)后，再令R＝I；

基于当前W(θ_p)，更新P(θ_p)为：

其中

基于当前P(θ_p)，更新R为：R＝R+P(θ_p)C_d(θ_p)C_d(θ_p)^H；

步骤S304：判断是否满足迭代收敛条件，若是，则输出当前功率谱P(θ_p)；否则继续执行步骤S303；

其中迭代收敛条件为：迭代次数达到上限，或者相邻两次迭代更新得到的功率谱的相对改变量小于或等于预设阈值；

步骤S4：根据功率谱P(θ_p)的峰值位置得到信号源的DOA估计结果。

Images