CN112763972B - 基于稀疏表示的双平行线阵二维doa估计方法及计算设备 - Google Patents
基于稀疏表示的双平行线阵二维doa估计方法及计算设备 Download PDFInfo
- Publication number
- CN112763972B CN112763972B CN202011615869.9A CN202011615869A CN112763972B CN 112763972 B CN112763972 B CN 112763972B CN 202011615869 A CN202011615869 A CN 202011615869A CN 112763972 B CN112763972 B CN 112763972B
- Authority
- CN
- China
- Prior art keywords
- matrix
- angle
- array
- subarray
- sparse representation
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 92
- 239000011159 matrix material Substances 0.000 claims abstract description 88
- 150000001875 compounds Chemical class 0.000 claims abstract description 25
- 239000013598 vector Substances 0.000 claims description 19
- 230000006870 function Effects 0.000 claims description 15
- 238000003491 array Methods 0.000 claims description 14
- 230000008569 process Effects 0.000 claims description 12
- 239000000654 additive Substances 0.000 claims description 6
- 230000000996 additive effect Effects 0.000 claims description 6
- 239000002131 composite material Substances 0.000 claims description 5
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 3
- 238000006243 chemical reaction Methods 0.000 claims description 2
- 238000000354 decomposition reaction Methods 0.000 abstract description 6
- 238000012545 processing Methods 0.000 abstract description 3
- 238000004422 calculation algorithm Methods 0.000 description 23
- 238000004088 simulation Methods 0.000 description 18
- 230000008859 change Effects 0.000 description 8
- 238000010586 diagram Methods 0.000 description 8
- 238000005259 measurement Methods 0.000 description 5
- 230000008901 benefit Effects 0.000 description 4
- 238000002474 experimental method Methods 0.000 description 4
- 230000009471 action Effects 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 3
- 238000013329 compounding Methods 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 2
- 238000000342 Monte Carlo simulation Methods 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 238000012733 comparative method Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000009977 dual effect Effects 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S3/00—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received
- G01S3/02—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received using radio waves
- G01S3/14—Systems for determining direction or deviation from predetermined direction
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/66—Radar-tracking systems; Analogous systems
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02D—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN INFORMATION AND COMMUNICATION TECHNOLOGIES [ICT], I.E. INFORMATION AND COMMUNICATION TECHNOLOGIES AIMING AT THE REDUCTION OF THEIR OWN ENERGY USE
- Y02D30/00—Reducing energy consumption in communication networks
- Y02D30/70—Reducing energy consumption in communication networks in wireless communication networks
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Mathematical Physics (AREA)
- Mathematical Analysis (AREA)
- Computational Mathematics (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- Computing Systems (AREA)
- Computer Networks & Wireless Communication (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)
Abstract
本发明提供一种基于稀疏表示的双平行线阵二维DOA估计方法及计算设备,涉及雷达信号处理技术领域。本发明不需要特征值分解,采用的是空间复合角解耦合思想,先利用稀疏表示方法求解其中的一个空间复合角,以此作为前提条件,另一个空间复合角就可以解耦为一维DOA估计问题,利用矩阵运算可以求解出来;最后根据已求解的两个空间复合角对方位角和俯仰角进行直接配对求解信源的DOA的数值。
Description
技术领域
本发明涉及雷达信号处理技术领域,具体涉及一种基于稀疏表示的双平行线阵二维DOA估计方法及计算设备。
背景技术
波达方向(Direction OfArrival,DOA)估计在阵列信号处理当中处于重要的地位,且得到了越来越多的研究与应用,特别是在雷达、通信、声呐领域。二维DOA估计对比一维DOA估计而言,二维DOA估计对空间的划分更加精细,对目标位置的定位更加具体,因此二维DOA估计更加具有实际意义。
现有的二维DOA估计方法大多是基于阵元间距等于半波长的简化面阵(即DOA矩阵方法),如L形阵列、双平行线阵、十字形阵列等。其中,双平行线阵由于结构简单、易于实现、具有较强的方法适用性等优点得到了广泛的关注和应用。现有的基于双平行线阵的二维DOA估计方法主要采用DOA矩阵方法,其结构示意图如图1所示,双平行线阵由阵元数为M、阵元间距为dx的两个均匀子阵Xa、Ya组成,子阵间距为dy。
目前,采用DOA矩阵方法要进行特征分解,这导致计算量大。
发明内容
(一)解决的技术问题
针对现有技术的不足,本发明提供了一种基于稀疏表示的双平行线阵二维DOA估计方法,解决了现有的DOA估计方法计算量大的技术问题。
(二)技术方案
为实现以上目的,本发明通过以下技术方案予以实现:
第一发明,本发明提供了一种基于稀疏表示的双平行线阵二维DOA估计方法,所述双平行线阵包括两个相互平行的子阵,所述二维DOA估计方法包括以下步骤:
S1、使用双平行线阵接收来自至少一个信源的信号,并定义信号的空间复合角αk、βk以及信号入射到Y轴方向的阵列流形矩阵A;
S2、基于稀疏表示方法获取βk的估计值;
S3、基于βk的估计值和阵列流形矩阵A获取阵列流形矩阵A的估计值;
S4、基于阵列流形矩阵A的估计值获取αk的估计值;
S5、基于βk的估计值和αk估计值估计所述信源的DOA的数值。
优选的,所述定义信号的空间复合角αk、βk以及信号入射到Y轴方向的阵列流形矩阵A,包括:
根据阵列的结构以及角度的定义,可得到子阵Y1、Y2对应的输出信号分别为:
Y1=As+n1 (2)
Y2=AΦs+n2 (3)其中:A=[a1,a2,…,aK]为信号入射到Y轴方向的阵列流形矩阵,s=[s1,…,sK]T为信号矢量,n1和n2分别表示子阵列Y1和子阵列Y2的接收噪声,且假设接收噪声是均值为0,方差为的高斯白噪声。
优选的,所述稀疏表示方法包括:
对所述阵列流型矩阵A进行扩展,构造一个以待估角度为参量的过完备基矩阵B。
优选的,所述对所述阵列流型矩阵进行扩展,构造一个以待估角度为参量的过完备基矩阵,包括:
设角度集合Ω={w1,w2,…,wH}包含了所有的空间角度βk,且有H>>K,H>>M,H为角度集合的元素个数,由此可以构造以wh,h=1,2,…,H为参数的过完备基矩阵B,
B=[a1(w1),a2(w2),…,ah(wh),…,aH(wH)] (4)
其中:ah(wh)被称为原子,由于H远远大于目标个数K和天线阵元个数M,且认为所有可能的空间角度βk都在角度集合Ω中,利用上式把子阵Y1对应的输出信号分别转化为稀疏表示问题:
Y1=Bq+n1 (5)
其中:q=[q1,q2,…,qh,…,qH]T,且当wh=βk时,有qh=sk,否则qh=0,通过求解矢量q,就可以根据其非零元素的位置得到复合空间角βk的估计值。
优选的,所述求解矢量q,包括:
求解矢量q等同于求解下述问题,即,
min||q||0 s.t. Y1=Bq+n1 (6)
其中:||q||0表示序列q中非零项的个数,式(6)所示的问题是非凸的,将式(6)转化为解决下述问题
min||q||0 s.t. Y1=Bq+n1 (7)
式(7)对噪声非常敏感,问题可以进一步转化为最小化目标函数
min||Y1-Bq||2+ρ||q||1 (8)
目标函数中前一项反映失配程度,后一项反映稀疏要求;由于q是复数,它的l1范数为:
采用二阶锥规划的方法,将其转化为
min||q||1 s.t. ||Y1-Bq||2≤η (10)
其中:η为正则化参数,其用于平衡B对q进行稀疏表示时的误差和矢量q的稀疏性,在正则化参数η确定后,式(10)问题转化为一个二阶锥规划问题,利用内点法可求解,得出复合空间角βk的估计值。
优选的,所述基于阵列流形矩阵A的估计值获取αk的估计值,包括:
子阵Y1的自相关矩阵为:
R11=E[Y1Y1 H]=ASAH+σ2I (11)
其中:E为数学期望算子,S为源信号的协方差矩阵,I是单位矩阵,σ2是加性噪声方差;
子阵Y1与子阵Y2的互相关矩阵为:
加性噪声是互不相关的并且是独立于源信号的零均值随机过程,因此式(12)后面三项等于零,有:
R12=E[Y1Y2 H]=AΦSAH (13)
由自相关矩阵R11和互相关矩阵R12构造DOA矩阵R为:
R=R12[R11]- (14)
其中:[R11]-表示伪逆;
如果A和S都满秩,Ф中没有相同的对角元素,那么DOA矩阵R的K个非零特征值与Ф中K个对角元素相等,并且与这些特征值对应的特征向量等于相应地信号导向矢量,即:
RA=AΦ (15)
第二方面,本发明提供了一种计算机可读存储介质,所述计算机可读存储介质用于存储程序代码,所述程序代码用于执行如上述所述的方法。
第三方面,本发明提供了一种计算设备,所述计算设备包括处理器以及存储器:
所述存储器用于存储程序代码,并将所述程序代码传输给所述处理器;
所述处理器用于根据所述程序代码中的指令执行如上述所述的方法。
(三)有益效果
本发明提供了一种基于稀疏表示的双平行线阵二维DOA估计方法。与现有技术相比,具备以下有益效果:
本发明首先使用双平行线阵接收来自至少一个信源的信号,并定义信号的空间复合角αk、βk以及信号入射到Y轴方向的阵列流形矩阵A;基于稀疏表示方法获取βk的估计值;基于βk的估计值和阵列流形矩阵A获取阵列流形矩阵A的估计值;基于阵列流形矩阵A的估计值获取αk的估计值;基于βk的估计值和αk估计值估计所述信源的DOA的数值。本发明不需要特征值分解,采用的是空间复合角解耦合思想,先利用稀疏表示方法求解其中的一个空间复合角,以此作为前提条件,另一个空间复合角就可以解耦为一维DOA估计问题,利用矩阵运算可以求解出来;最后根据已求解的两个空间复合角对方位角和俯仰角进行直接配对求解信源的DOA的数值。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其他的附图。
图1为现有的DOA矩阵方法的双平行线阵的结构示意图;
图2为本发明实施例的基于稀疏表示的双平行线阵二维DOA估计方法的框图;
图3为本发明实施例的双线性平行阵列的结构示意图;
图4为仿真1的实验结果角度估计图;
图5为仿真2的各算法测角RMSE随SNR的变化的示意图;
图6为仿真3的各算法测角RMSE随快拍数的变化的示意图;
图7为仿真4的各算法测角RMSE随方位角的变化的示意图;
图8为仿真4的各算法测角RMSE随俯仰角的变化的示意图;
图9为仿真4的各算法测角RMSE随方位和俯仰角的变化的示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
本申请实施例通过提供一种基于稀疏表示的双平行线阵二维DOA估计方法,解决了现有的DOA估计方法计算量大的技术问题,实现降低计算量。
本申请实施例中的技术方案为解决上述技术问题,总体思路如下:
目前稀疏表示技术在DOA估计领域得到了广泛的研究,并取得了较好的结果。现在大部分DOA算法都是针对一维来展开研究的,更多时候往往需要对二维DOA进行估计。为了解决这个问题,一个可行的方法就是将现有的一维DOA方法直接拓展到二维DOA估计中,虽然也能够估计出角度信息,但是这种方法在进行二维DOA估计时需要将阵元由一维结构拓展到二维平面,如面阵、圆阵等,这将会导致阵元个数大大增加,相应地算法复杂度也会增加。另一个可行的方法是将二维DOA估计问题转换成两个一维DOA估计,这样就能够利用已经成熟的一维DOA估计算法。本发明实施例基于第二种思路,采用双线性平行阵列,能够在较少阵元的前提下,利用稀疏表示技术,不需要进行需谱峰搜索以及特征值分解,且估计角度参数自动配对,确保了算法的估计性能,降低计算复杂度,提高了方法的适用性。
为了更好的理解上述技术方案,下面将结合说明书附图以及具体的实施方式对上述技术方案进行详细的说明。
本发明实施例提供了一种基于稀疏表示的双平行线阵二维DOA估计方法,双平行线阵包括两个相互平行的子阵。如图2所示,该方法包括:
S1、使用双平行线阵接收来自至少一个信源的信号,并定义信号的空间复合角αk、βk以及信号入射到Y轴方向的阵列流形矩阵A;
S2、基于稀疏表示方法获取βk的估计值;
S3、基于βk的估计值和阵列流形矩阵A获取阵列流形矩阵A的估计值;
S4、基于阵列流形矩阵A的估计值获取αk的估计值;
S5、基于βk的估计值和αk估计值估计所述信源的DOA的数值。
本发明实施例不需要特征值分解,采用的是空间复合角解耦合思想,先利用稀疏表示方法求解其中的一个空间复合角,以此作为前提条件,另一个空间复合角就可以解耦为一维DOA估计问题,利用矩阵运算可以求解出来;最后根据已求解的两个空间复合角对方位角和俯仰角进行直接配对求解信源的DOA的数值。
下面对各个步骤进行详细说明:
在步骤S1中,使用双平行线阵接收来自至少一个信源的信号,并定义信号的空间复合角αk、βk以及信号入射到Y轴方向的阵列流形矩阵A。具体实施过程如下:
本发明实施例采用的双平行线阵如图3所示,由两个子阵Y1和Y2组成,子阵的阵元数为M、子阵的阵元间距为dY,子阵Y1与子阵Y2之间的间距为dX。假设有K个非相干窄带信号入射到该双平行线阵,第k个信号的入射方向可描述为(αk,βk),αk和βk叫做空间复合角,分别定义为入射方向同X轴和Y轴正方向的夹角。同时定义第k个信号的方位角和俯仰角分别为那么由图3中的几何位置关系可知:
根据阵列的结构以及角度的定义,可得到子阵列Y1、Y2对应的输出信号分别为
Y1=As+n1 (2)
Y2=AФs+n2 (3)其中:A=[a1,a2,…,aK]为信号入射到Y轴方向的阵列流形矩阵,s=[s1,…,sK]T为信号矢量,n1和n2分别表示子阵列Y1和子阵列Y2的接收噪声,且假设接收噪声是均值为0,方差为的高斯白噪声。
在步骤S2中,基于稀疏表示方法获取βk的估计值。具体实施过程如下:
为了用稀疏表示方法对空间角度βk进行估计,通常需要把阵列流型矩阵A进行扩展,形成过完备基矩阵B,因而需要构造一个以待估角度为参量的过完备基矩阵。假设角度集合Ω={w1,w2,…,wH}包含了所有的空间角度βk,且有H>>K,H>>M。H为角度集合的元素个数,由此可以构造以wh,h=1,2,…,H,为参数的过完备基矩阵B。
B=[a1(w1),a2(w2),…,ah(wh),…,aH(wH)] (4)
其中:ah(wh)被称为原子。由于H远远大于目标个数K和天线阵元个数M,且认为所有可能的空间角度βk都在角度集合Ω中,利用式(4)可以把式(2)转化为稀疏表示问题:
Y1=Bq+n1 (5)
其中:q=[q1,q2,…,qh,…,qH]T,且当ωh=βk时,有qh=sk,否则qh=0,通过求解矢量q,就可以根据其非零元素的位置得到复合空间角βk的估计值。
求解矢量q的过程如下:
由于H远远大于目标个数K和天线阵元个数M,对式(5)中的q进行求解为一欠定问题,即存在无穷多组解。但是根据前面的分析可以知道矢量q是稀疏的,则矢量q等同于求解下述问题,即,
min||q||0 s.t. Y1=Bq+n1 (6)
其中:||q||0表示序列q中非零项的个数,由稀疏表示基础理论可知,式(6)所示的问题是非凸的,很难进行直接求解。为了解决这个难点,将其转化为解决下述问题:
min||q||1 s.t. Y1=Bq+n1 (7)
式(7)对噪声非常敏感,微量的噪声就可能导致问题最优解发生很大改变,问题可以进一步转化为最小化目标函数:
min||Y1-Bq||2+ρ||q||1 (8)
目标函数中前一项反映失配程度,后一项反映稀疏要求。由于q是复数,它的l1范数为:
由于上式的平方仍不能消除平方根项,直接导致不能使用二次规划的方式最小化目标函数。为了解决这个问题,采用二阶锥规划的方法,将其转化为:
min||q||1 s.t.||Y1-Bq||2≤η (10)
其中:η为正则化参数,其用于平衡B对q进行稀疏表示时的误差和矢量q的稀疏性。在正则化参数η确定后,式(10)问题转化为一个二阶锥规划问题,利用内点法可求解,求解出βk的估计值。
在步骤S3中,基于βk的估计值和阵列流形矩阵A获取阵列流形矩阵A的估计值。具体实施过程如下:
在步骤S4中,基于阵列流形矩阵A的估计值获取αk的估计值。具体实施过程如下:
子阵列Y1的自相关矩阵为:
R11=E[Y1Y1 H]=ASAH+σ2I (11)
其中:E为数学期望算子,S为源信号的协方差矩阵,I是单位矩阵,σ2是加性噪声方差。
子阵列Y1与子阵列Y2的互相关矩阵为:
由前面假设,加性噪声是互不相关的并且是独立于源信号的零均值随机过程,因此式(12)后面三项等于零。有:
R12=E[Y1Y2 H]=AΦSAH (13)
由自相关矩阵R11和互相关矩阵R12构造DOA矩阵R为
R=R12[R11]- (14)
其中:[R11]-表示伪逆。
定理1如果A和S都满秩,Φ中没有相同的对角元素,那么DOA矩阵R的K个非零特征值与Φ中K个对角元素相等,并且与这些特征值对应的特征向量等于相应地信号导向矢量,即:
RA=AΦ (15)
在步骤S5中,基于βk的估计值和αk估计值估计所述信源的DOA的数值。具体实施过程如下:
利用分维的方法进行方位角和俯仰角估计时,一个难题就是方位角和俯仰角的配对问题。例如算法要进行特征值分解的时候,特征矢量的排序有可能不同,导致估算出来的方位角和俯仰角不能直接配对。本发明实施例依次求解两个空间复合角αk和βk,然后根据式(1)可以反推出方位角θ和俯仰角φ:
为了验证本发明实施例的有效性,下面通过仿真实验进行对比说明。
在相同实验条件下同时给出了DOA矩阵算法,ESPRIT算法和MEMP算法的实验结果。
实验中假定子阵Y1与子阵Y2的阵元个数M为10,两子阵列的阵元沿Y轴方向的间距dY等于半波长,子阵列与子阵列的间距dX等于半波长。设各目标入射信号的功率均相等,各算法的估计均方根误差(RMSE)定义为
仿真1:下面首先对所提方法(即本发明实施例的方法)的估计角度情况进行分析。假设有两个非相干远场窄带信号入射到双线性平行阵列上面,并且有两入射信号的入射方位角和俯仰角分别为和接收的快拍数设定为500,信噪比SNR=15dB。仿真实验结果如图4所示,由图可知所提方法能较精确的分辨这两个信号的方位角和俯仰角。
仿真2:下面对各算法估计角度RMSE随信噪比的变化情况进行分析。在该仿真实验中设定两入射信号的入射方位角和俯仰角同仿真实验1相同,接收的快拍数为500,单个阵元接收信号的信噪比由0dB变化到20dB,则由此给出各算法的估计角度RMSE随信噪比的变化情况如图5所示。
由图5所示的仿真结果可知,当信噪比较高时,所提方法的估计角度RMSE具有较好的结果,而当信噪比较低时所提方法的估计角度RMSE变化较大,性能下降。其原因在于信噪比较低时,正则化参数的值有偏差,导致算法性能下降。
仿真3:下面对各算法估计角度RMSE随快拍数的变化情况进行分析。在该仿真实验中设定两入射信号的入射方位角和俯仰角同仿真实验1相同,且入射信号的信噪比SNR=15dB,改变接收的快拍数,由50变化到1000,则由此给出各算法的估计角度RMSE随快拍数的变化情况如图6所示。
由图6可知,所提方法的性能较对比方法受快拍数的影响较小,因而在快拍数较小时,其优势更加明显。另外,DOA矩阵法和MEMP算法在快拍数较大时能够对目标角度进行高精度的估计,快拍数较小时估计的偏差较大。所提方法具有较好的稳健型,这得益于两子阵列接收信号互相关信息的利用及稀疏表示技术的应用。
仿真4:本实验对入射信号的空域角度间隔对各算法的性能影响进行分析。在实验中仿真以下三种情况:(1)两个入射信号的俯仰角间隔大并且固定不变,考查各算法性能随方位角间隔变化情况,此时设定(2)两个入射信号的方位角间隔较大且固定不变,考查各算法性能随俯仰角间隔变化情况,此时设定(3)两个入射信号的方位角和俯仰角间隔同时由小变大,此时设定其中Δθ、Δ均由2°变化到20°。此外入射信号的信噪比SNR=15dB,接收的快拍数设为500。在上述条件下,得到各算法的估计角度RMSE随方位角间隔Δθ、俯仰角间隔以及方位和俯仰角间隔Δ的变化情况分别如图7、图8、图9所示。
图7~图9的仿真结果可以看出各算法性能受到俯仰角间隔影响较大,而受方位角间隔影响较小。其原因在于各算法进行角度测量时,首先进行目标俯仰角的测量,然后再根据俯仰角测量结果对方位角进行计算。因此,当俯仰角间隔足够大时,目标俯仰角间隔可以准确估计,使得方位角测量性能几乎不受方位角间隔的影响。
基于同一发明构思,本申请实施例还提供了一种计算机可读存储介质,其特征在于,计算机可读存储介质用于存储程序代码,程序代码用于执行权利要求任一项所述的方法。
基于同一发明构思,本申请实施例还提供了一种计算设备,其特征在于,计算设备包括处理器以及存储器:
存储器用于存储程序代码,并将所述程序代码传输给所述处理器;
处理器用于根据所述程序代码中的指令执行权利要求任一项所述的方法。
综上所述,与现有技术相比,具备以下有益效果:
1、本发明实施例的方法不需要特征值分解,采用的是空间复合角解耦合思想,先利用稀疏表示方法求解其中的一个空间复合角,以此作为前提条件,另一个空间复合角就可以解耦为一维DOA估计问题,利用矩阵运算可以求解出来;最后根据已求解的两个空间复合角对方位角和俯仰角进行直接配对求解信源的DOA的数值。
2、和现有算法相比较,本发明实施例的方法受快拍数的影响较小,在信噪比较高、角度间隔较大的情况下,具有良好的性能。
3、本发明实施例的方法适用于其它阵列,如T型阵列、L型阵列上,适用性较强。
需要说明的是,在本文中,诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来,而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且,术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含,从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素,而且还包括没有明确列出的其他要素,或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下,由语句“包括一个……”限定的要素,并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。
以上实施例仅用以说明本发明的技术方案,而非对其限制;尽管参照前述实施例对本发明进行了详细的说明,本领域的普通技术人员应当理解:其依然可以对前述各实施例所记载的技术方案进行修改,或者对其中部分技术特征进行等同替换;而这些修改或者替换,并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。
Claims (6)
1.一种基于稀疏表示的双平行线阵二维DOA估计方法,其特征在于,所述双平行线阵包括两个相互平行的子阵,子阵Y1和Y2,子阵的阵元数为M、子阵的阵元间距为dY,子阵Y1与子阵Y2之间的间距为dX;所述二维DOA估计方法包括以下步骤:
S1、使用双平行线阵接收来自至少一个信源的信号,并定义信号的空间复合角αk、βk以及信号入射到Y轴方向的阵列流形矩阵A,包括:
根据阵列的结构以及角度的定义,可得到子阵Y1、Y2对应的输出信号分别为:
Y1=As+n1 (2)
Y2=AΦs+n2 (3)
其中:A=[a1,a2,…,aK]为信号入射到Y轴方向的阵列流形矩阵,
S2、基于稀疏表示方法获取βk的估计值;
S3、基于βk的估计值和阵列流形矩阵A获取阵列流形矩阵A的估计值;
S4、基于阵列流形矩阵A的估计值获取αk的估计值,包括:
子阵Y1的自相关矩阵为:
R11=E[Y1Y1 H]=ASAH+σ2I (11)
其中:E为数学期望算子,S为源信号的协方差矩阵,I是单位矩阵,σ2是加性噪声方差;
子阵Y1与子阵Y2的互相关矩阵为:
加性噪声是互不相关的并且是独立于源信号的零均值随机过程,因此式(12)后面三项等于零,有:
R12=E[Y1Y2 H]=AΦSAH (13)
由自相关矩阵R11和互相关矩阵R12构造DOA矩阵R为:
R=R12[R11]- (14)
其中:[R11]-表示伪逆;
如果A和S都满秩,Φ中没有相同的对角元素,那么DOA矩阵R的K个非零特征值与Φ中K个对角元素相等,并且与这些特征值对应的特征向量等于相应地信号导向矢量,即:
RA=AΦ (15)
S5、基于βk的估计值和αk估计值估计所述信源的DOA的数值。
2.如权利要求1所述的基于稀疏表示的双平行线阵二维DOA估计方法,其特征在于,所述稀疏表示方法包括:
对所述阵列流型矩阵A进行扩展,构造一个以待估角度为参量的过完备基矩阵B。
3.如权利要求2所述的基于稀疏表示的双平行线阵二维DOA估计方法,其特征在于,所述对所述阵列流型矩阵进行扩展,构造一个以待估角度为参量的过完备基矩阵,包括:
设角度集合Ω={w1,w2,…,wH}包含了所有的空间角度βk,且有H>>K,H>>M,H为角度集合的元素个数,由此构造以wh,h=1,2,…,H为参数的过完备基矩阵B,
B=[a1(w1),a2(w2),…,ah(wh),…,aH(wH)] (4)
其中:ah(wh)被称为原子,由于H远远大于目标个数K和天线阵元个数M,且认为所有可能的空间角度βk都在角度集合Ω中,利用上式把子阵Y1对应的输出信号分别转化为稀疏表示问题:
Y1=Bq+n1 (5)
其中:q=[q1,q2,…,qh,…,qH]T,且当wh=βk时,有qh=sk,否则qh=0,通过求解矢量q,就根据其非零元素的位置得到复合空间角βk的估计值。
4.如权利要求3所述的基于稀疏表示的双平行线阵二维DOA估计方法,其特征在于,所述求解矢量q,包括:
求解矢量q等同于求解下述问题,即,
min||q||0 s.t.Y1=Bq+n1 (6)
其中:||q||0表示序列q中非零项的个数,式(6)所示的问题是非凸的,将式(6)转化为解决下述问题
min||q||1 s.t.Y1=Bq+n1 (7)
式(7)对噪声非常敏感,问题进一步转化为最小化目标函数
min||Y1-Bq||2+ρ||q||1 (8)
目标函数中前一项反映失配程度,后一项反映稀疏要求;由于q是复数,它的l1范数为:
采用二阶锥规划的方法,将其转化为
min||q||1 s.t.||Y1-Bq||2≤η (10)
其中:η为正则化参数,其用于平衡B对q进行稀疏表示时的误差和矢量q的稀疏性,在正则化参数η确定后,式(10)问题转化为一个二阶锥规划问题,利用内点法可求解,得出复合空间角βk的估计值。
5.一种计算机可读存储介质,其特征在于,所述计算机可读存储介质用于存储程序代码,所述程序代码用于执行权利要求1~4任一项所述的方法。
6.一种计算设备,其特征在于,所述计算设备包括处理器以及存储器:
所述存储器用于存储程序代码,并将所述程序代码传输给所述处理器;
所述处理器用于根据所述程序代码中的指令执行权利要求1~4任一项所述的方法。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202011615869.9A CN112763972B (zh) | 2020-12-30 | 2020-12-30 | 基于稀疏表示的双平行线阵二维doa估计方法及计算设备 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202011615869.9A CN112763972B (zh) | 2020-12-30 | 2020-12-30 | 基于稀疏表示的双平行线阵二维doa估计方法及计算设备 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN112763972A CN112763972A (zh) | 2021-05-07 |
CN112763972B true CN112763972B (zh) | 2023-05-09 |
Family
ID=75697830
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202011615869.9A Active CN112763972B (zh) | 2020-12-30 | 2020-12-30 | 基于稀疏表示的双平行线阵二维doa估计方法及计算设备 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN112763972B (zh) |
Families Citing this family (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113341371B (zh) * | 2021-05-31 | 2022-03-08 | 电子科技大学 | 一种基于l阵和二维esprit算法的doa估计方法 |
CN113885004A (zh) * | 2021-09-29 | 2022-01-04 | 哈尔滨工程大学 | 一种稀疏面阵二维方位估计方法 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108594164A (zh) * | 2017-11-30 | 2018-09-28 | 山东农业大学 | 一种平面阵列doa估计方法及设备 |
JP2019090749A (ja) * | 2017-11-16 | 2019-06-13 | 富士通株式会社 | 到来方向推定装置および到来方向推定方法 |
CN110297209A (zh) * | 2019-04-08 | 2019-10-01 | 华南理工大学 | 一种基于平行互质阵列时空扩展的二维波达方向估计方法 |
Family Cites Families (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106483493B (zh) * | 2016-09-13 | 2018-12-18 | 电子科技大学 | 一种稀疏双平行线阵及二维波达方向估计方法 |
CN111707985A (zh) * | 2020-06-15 | 2020-09-25 | 浙江理工大学 | 基于协方差矩阵重构的off-grid DOA估计方法 |
-
2020
- 2020-12-30 CN CN202011615869.9A patent/CN112763972B/zh active Active
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP2019090749A (ja) * | 2017-11-16 | 2019-06-13 | 富士通株式会社 | 到来方向推定装置および到来方向推定方法 |
CN108594164A (zh) * | 2017-11-30 | 2018-09-28 | 山东农业大学 | 一种平面阵列doa估计方法及设备 |
CN110297209A (zh) * | 2019-04-08 | 2019-10-01 | 华南理工大学 | 一种基于平行互质阵列时空扩展的二维波达方向估计方法 |
Non-Patent Citations (1)
Title |
---|
Jianfeng Li 等.Sparse representation based two-dimensional direction of arrival estimation using co-prime array.《Multidimensional Systems and Signal Processing》.2016,(第29期),35-47. * |
Also Published As
Publication number | Publication date |
---|---|
CN112763972A (zh) | 2021-05-07 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN110197112B (zh) | 一种基于协方差修正的波束域Root-MUSIC方法 | |
CN109633522B (zh) | 基于改进的music算法的波达方向估计方法 | |
CN107576940B (zh) | 一种低复杂度单基地mimo雷达非圆信号角度估计方法 | |
CN105445709B (zh) | 一种稀布阵列近场无源定位幅相误差校正方法 | |
CN107870314B (zh) | 基于极化敏感阵列的完备电磁分量加权融合测向优化方法 | |
CN104407335B (zh) | 一种3轴交叉阵列的doa估计方法 | |
CN110161489B (zh) | 一种基于伪框架的强弱信号测向方法 | |
CN112763972B (zh) | 基于稀疏表示的双平行线阵二维doa估计方法及计算设备 | |
CN111046591B (zh) | 传感器幅相误差与目标到达角度的联合估计方法 | |
CN109116293A (zh) | 一种基于离格稀疏贝叶斯的波达方向估计方法 | |
CN113835063B (zh) | 一种无人机阵列幅相误差与信号doa联合估计方法 | |
CN103353588B (zh) | 基于天线均匀平面阵的二维波达方向角估计方法 | |
CN112379327A (zh) | 一种基于秩损估计的二维doa估计与互耦校正方法 | |
CN109490820A (zh) | 一种基于平行嵌套阵的二维doa估计方法 | |
CN109254272B (zh) | 一种共点式极化mimo雷达的两维角度估计方法 | |
CN104020440B (zh) | 基于l型干涉式线性阵列的二维波达角估计方法 | |
CN106970348B (zh) | 电磁矢量传感器阵列解相干二维music参数估计方法 | |
CN106249196A (zh) | 三分量声矢量传感器稀疏阵列四元数解模糊方法 | |
CN104793177B (zh) | 基于最小二乘法的麦克风阵列测向方法 | |
CN110286350A (zh) | 一种l型稀疏阵doa估计的精确配对方法及装置 | |
CN108983145A (zh) | 电磁矢量传感器阵列宽带相干源定位方法 | |
CN110579737B (zh) | 一种杂波环境中基于稀疏阵列的mimo雷达宽带doa计算方法 | |
CN109696651B (zh) | 一种基于m估计的低快拍数下波达方向估计方法 | |
CN111368256B (zh) | 一种基于均匀圆阵的单快拍测向方法 | |
CN114460531A (zh) | 一种均匀线阵music空间谱估计方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |