CN107300686A - 基于多项式求解的非圆信号波达方向角的估计方法 - Google Patents

基于多项式求解的非圆信号波达方向角的估计方法 Download PDF

Info

Publication number
CN107300686A
CN107300686A CN201710423912.3A CN201710423912A CN107300686A CN 107300686 A CN107300686 A CN 107300686A CN 201710423912 A CN201710423912 A CN 201710423912A CN 107300686 A CN107300686 A CN 107300686A
Authority
CN
China
Prior art keywords
vector
noise
matrix
array
span
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.)
Granted
Application number
CN201710423912.3A
Other languages
English (en)
Other versions
CN107300686B (zh
Inventor
蔡晶晶
武斌
张垚均
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xidian University
Original Assignee
Xidian University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Xidian University filed Critical Xidian University
Priority to CN201710423912.3A priority Critical patent/CN107300686B/zh
Publication of CN107300686A publication Critical patent/CN107300686A/zh
Application granted granted Critical
Publication of CN107300686B publication Critical patent/CN107300686B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Direction-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/02Direction-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/14Systems for determining direction or deviation from predetermined direction
    • G01S3/143Systems for determining direction or deviation from predetermined direction by vectorial combination of signals derived from differently oriented antennae
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Direction-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/78Direction-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 electromagnetic waves other than radio waves
    • G01S3/782Systems for determining direction or deviation from predetermined direction

Abstract

本发明公开了一种基于多项式求解的非圆信号波达方向角估计方法,主要解决现有技术中阵元利用率低,信号识别数量少的问题,其实现步骤是:1)获取嵌套阵列输出信号,根据该信号计算协方差矩阵和椭圆协方差矩阵,并构造等效协方差向量和等效椭圆协方差向量,计算这两个向量中所有元素的维数;2)计算虚拟阵列协方差向量和虚拟阵列椭圆协方差向量,并构造两个波达角选择矩阵,计算其噪声子空间;3)由噪声子空间获得第一噪声矩阵、第二噪声矩阵、第三噪声矩阵和第四噪声矩阵,根据这两个噪声矩阵构造多项式方程;4)计算多项式方程的根获得目标波达方向角度值。本发明在非圆信号环境下大大提高了阵列可识别的信源数目,可用于目标侦察和无源定位。

Description

基于多项式求解的非圆信号波达方向角的估计方法
技术领域
本发明属于信号处理技术领域,特别涉及一种电磁信号的阵列信号波达方向角估计方法,可用于对飞机、舰船运动目标的侦察与无源定位。
背景技术
信号的波达方向角DOA估计是阵列信号处理领域的一个重要分支,它是指利用天线阵列对空间声学信号、电磁信号进行感应接收,再运用现代信号处理方法快速准确的估计出信号源的方向,在雷达、声纳、无线通信等领域具有重要应用价值。
在现代通信中,二相相移键控以及M进制幅移键控等非圆信号的应用越来越多,因此有关非圆信号的DOA估计受到了越来越多的关注。目前已有关于利用阵列天线处理非圆信号的一些方法被提出,比较有代表性的是P Charge等人发表的论文“A non-circularsources direction finding method using polynomial rooting”(《SignalProcessing》,VOL 81,pp.1765-1770 2001)中公开了一种利用多项式求解进行非圆信号DOA估计的方法,该方法是基于均匀阵列的。
另一方面,为了在较少的阵元条件下得到尽量大的角度自由度,检测更多的信源,一些新的阵列结构被提出,比较有代表性的是嵌套阵列以及互质阵列。P Piya等人在其发表的论文“Nested Arrays:A Novel Approach to Array Processing With EnhancedDegrees of Freedom”(《IEEE transactions on signal processing》,VOL58,NO.8,August 2010)中公开了一种基于嵌套阵列的DOA估计方法,该方法能够使用M+N个阵元,生成2MN+2N-1个虚拟阵元,可检测MN+N-1个信号。该方法具有估计多于阵元数目的信号数的能力,但是,该阵列的讨论都集中在接收信号为圆信号的条件下,对于如何利用该阵列进行非圆信号的处理目前还没有研究。
在实际应用中,具有非圆信号环境,给定一定数量的阵元,如果不能合理利用这些阵元以及信号的非圆特性,就不能估计足够多的信号,造成侦察和定位资源的浪费。
发明内容
本发明的目的在于针对上述现有技术存在的不足,提出一种基于多项式求解的非圆信号波达方向角估计方法,以在非圆信号环境下,合理利用阵元及信号的非圆特性对足够多的信号进行估计,避免资源浪费。
为实现上述目的,本发明技术方案包括如下:
(1)用M+N个天线接收机形成嵌套阵列,其中M、N分别表示两个天线接收阵列的阵元数,其取值范围为M≥1,N≥1;
(2)假设空间中有K个非圆目标信号,通过采样和滤波得到嵌套阵列输出信号:Y(t)=[y1(t),…,yi(t),…,yM+N(t)]T,其中,yi(t)表示嵌套阵列的第i个阵元的输出信号,t的取值范围是1≤t≤L,L表示快拍数,i的取值范围是1≤i≤M+N,(·)T表示矩阵转置运算;
(3)利用嵌套阵列输出信号Y(t),计算虚拟均匀阵列协方差向量和椭圆协方差向量
(4)根据虚拟均匀阵列协方差向量得到波达角选择矩阵G,根据虚拟均匀阵列椭圆协方差向量得到波达角选择矩阵Q;
(5)分别计算波达角选择矩阵G的噪声子空间Un和波达角选择矩阵Q的噪声子空间Unn
(6)根据波达角选择矩阵G的噪声子空间Un和波达角选择矩阵Q的噪声子空间Unn,获得不同的噪声矩阵:
(6a)提取波达角选择矩阵G的噪声子空间Un的前L1行(L1+L2-K)列的所有元素,生成L1×(L1+L2-K)维的第一子矩阵,将该第一子矩阵作为第一噪声矩阵Un1;用噪声子空间Un的后L2行(L1+L2-K)列所有元素生成L2×(L1+L2-K)维的第二子矩阵,将该第二子矩阵作为第二噪声矩阵Un2
(6b)提取波达角选择矩阵Q的噪声子空间Unn的前L2行(L1+L2-K)列所有元素生成L2×(L1+L2-K)维的第一子矩阵,将该第一子矩阵作为第三噪声矩阵Unn1;用噪声子空间Unn的后L1行(L1+L2-K)列所有元素生成L1×(L1+L2-K)维的第二子矩阵,将该第二子矩阵作为第四噪声矩阵Unn2
(7)根据第一噪声矩阵Un1,第二噪声矩阵Un2,第三噪声矩阵Unn1和第四噪声矩阵Unn2构造八个噪声向量c1,c2,c3,c4,c5,c6,c7和c8
(8)根据c1,c2,c3,c4,c5,c6,c7和c8这八个噪声向量,构造第一复合向量p14和第一最终向量q23,得到如下多项式方程:
其中,p14(j)表示第一复合向量p14中第j个元素,q23(j)表示第一最终向量q23中第j个元素,j的取值范围是1≤j≤2(L1+L2)-3,z表示多项式方程的根,z=[z1,…,zh,…,zK],zh表示多项式方程的第h个根,该多项式最多有2L=3L1+L2-4个根,且每个根都有一个与其相似的根,每对相似根中只保留其中一个,得到最多有L个根z1,…zn,….zL
(9)计算多项式方程的所有根,由多项式方程的每一个根的辐角与目标波达方向角度值的关系,得到目标波达方向角度值θ。
本发明与现有技术相比具有以下优点:
1)本发明采用了嵌套阵列模型进行波达方向角度估计,克服了现有技术中采用典型的线性均匀阵列,造成估计的信号数目低于阵元数目的缺点,提高了在阵元数目相同的条件下的阵列可识别信源数目。
2)本发明利用信号的协方差矩阵Rd,以及椭圆协方差矩阵Rs对信号进行估计,增多了可估计的非圆信号数目。
3)本发明利用了嵌套阵列的特点和非圆信号的特性,将非圆信号在嵌套阵列上进行信号处理,在M+N个阵元上能估计出(MN+M+N-1)/2+MN+N-1个信号,大大提高了阵列的利用率,增加了阵列可识别信源的数目。
附图说明
图1是本发明的实现流程图;
图2是本发明中嵌套阵列的结构示意图。
具体实施方式
参照图1,本实例的实现步骤如下:
步骤1:用M+N个天线接收机形成嵌套阵列。
(1a)将每个天线接收机称为一个阵元,用M个天线接收机形成第一均匀线性阵列a,其阵元间距为d,定义第一均匀线性阵列a的第一个阵元为起始阵元,定义起始阵元位置D(1)=1,第一均匀线性阵列a的其它阵元位置依次为D(2)=2,D(3)=3,D(4)=4,…,D(M)=M;其中,M的取值范围为M≥1,d的取值范围为0<d≤λ/2,λ为入射到阵列的窄带信号波长;
(1b)用N个天线接收机形成第二均匀线性阵列b,其阵元间距为(M+1)d,第二均匀线性阵列b的阵元位置依次设置为D(M+1)=M+1,D(M+2)=2(M+1),D(M+2)=3(M+1),…,D(M+N)=N(M+1),其中,N的取值范围为N≥1;
(1c)将第二均匀线性阵列b的第一个阵元放置于与起始阵元相距为Md的位置,将第二均匀线性阵列b的所有阵元依次放置于第一均匀线性阵列a之后,形成嵌套阵列。
步骤2:获得嵌套阵列输出信号Y(t)。
假设空间中有K个非圆目标信号,使用嵌套阵列天线接收机,对空间目标信号进行快拍采样和匹配滤波操作,得到嵌套阵列输出信号:Y(t)=[y1(t),…,yi(t),…,yM+N(t)]T,其中,K的取值范围是K<MN+M+N-1,yi(t)表示嵌套阵列第i个阵元的输出信号,t的取值范围是1≤t≤L,L表示快拍数,i的取值范围是1≤i≤M+N,(·)T表示矩阵转置运算。
步骤3:计算协方差矩阵Rd和椭圆协方差矩阵Rs
利用嵌套阵列输出信号Y(t),计算协方差矩阵Rd和椭圆协方差矩阵Rs
其中,(·)H表示矩阵共轭转置运算。
步骤4:构造等效协方差向量rd和等效椭圆协方差向量rs
根据协方差矩阵Rd和椭圆协方差矩阵Rs中的元素,分别构造等效协方差向量rd和等效椭圆协方差向量rs
rd=[Rd(1,1),Rd(2,1),…,Rd(M+N,1),Rd(1,2),…,Rd(M+N,2),
...,Rd(i,j),…,Rd(1,M+N),…,Rd(M+N,M+N)]T
rs=[Rs(1,1),Rs(2,1),...,Rs(M+N,1),Rs(1,2),…,Rs(M+N,2),
…,Rs(i,j),…,Rs(1,M+N),…,Rs(M+N,M+N)]T
其中,Rd(i,j)表示协方差矩阵Rd中位于第i行、第j列的元素,i的取值范围为1≤i≤M+N,j的取值范围为1≤j≤M+N,Rs(i,j)表示椭圆协方差矩阵Rs中位于第i行,第j列的元素。
步骤5:计算等效协方差向量和等效椭圆协方差向量中所有元素的维数。
根据等效协方差向量rd和等效椭圆协方差向量rs中每一个元素所在的行和列在嵌套阵列中对应的阵元位置,计算等效协方差向量rd中所有元素的维数Ei,j和等效椭圆协方差向量rs中所有元素的维数Fi,j
Ei,j=D(j)-D(i),
Fi,j=D(j)+D(i),
其中,D(i)表示嵌套阵列中第i个阵元的位置,D(j)表示嵌套阵列中第j个阵元的位置。
步骤6:获得虚拟均匀阵列协方差向量和虚拟均匀阵列椭圆协方差向量
根据等效协方差向量rd中所有元素的维数,删除等效协方差向量rd中维数相同的元素和维数不连续的元素,并将剩余元素按维数从小到大排列,得到虚拟均匀阵列协方差向量
根据等效椭圆协方差向量rs中所有元素的维数,删除等效椭圆协方差向量rs中维数相同的元素和维数不连续的元素,并将剩余元素按维数从小到大排列,得到虚拟均匀阵列椭圆协方差向量
步骤7:构造波达角选择矩阵G和波达角选择矩阵Q。
根据虚拟均匀阵列协方差向量和虚拟均匀阵列椭圆协方差向量中所有元素进行行列排列,得到波达角选择矩阵G和Q:
其中,中间变量L1=(Cd+1)/2,L2=Cs+1-(Cd+1)/2,且有L1>L2,Cd表示虚拟均匀阵列协方差向量中元素的个数,Cd的取值为2MN+2N-1,Cs表示虚拟均匀阵列椭圆协方差向量中元素的个数,Cs的取值为MN+M+N,(·)*表示向量的共轭运算。
步骤8:计算波达角选择矩阵G的噪声子空间Un,波达角选择矩阵Q的噪声子空间Unn
(8a)将波达角选择矩阵G进行如下特征分解:
G=U·Λ·UH
其中,Λ为波达角选择矩阵G的特征值矩阵,U为矩阵G的特征值所对应的特征向量矩阵,(·)H表示矩阵的共轭转置运算;
(8b)将特征值矩阵Λ中的特征值按从大到小排序,取后(L1+L2-K)个较小特征值对应的特征向量矩阵作为噪声子空间Un
(8c)将波达角选择矩阵Q进行如下特征分解:
Q=U·Λ·UH
其中,Λ为波达角选择矩阵Q的特征值矩阵,U为矩阵Q的特征值所对应的特征向量矩阵,(·)H表示矩阵的共轭转置运算;
(8d)将特征值矩阵Λ中的特征值按从大到小排序,取后(L1+L2-K)个较小特征值对应的特征向量矩阵作为噪声子空间Unn
步骤9:根据G的噪声子空间Un和波达角选择矩阵Q的噪声子空间Unn,获得不同的噪声矩阵。
(9a)提取波达角选择矩阵G的噪声子空间Un的前L1行(L1+L2-K)列的所有元素,生成L1×(L1+L2-K)维的第一子矩阵,将该第一子矩阵作为第一噪声矩阵Un1;用噪声子空间Un的后L2行(L1+L2-K)列所有元素生成L2×(L1+L2-K)维的第二子矩阵,将该第二子矩阵作为第二噪声矩阵Un2
(9b)提取波达角选择矩阵Q的噪声子空间Unn的前L2行(L1+L2-K)列所有元素生成L2×(L1+L2-K)维的第一子矩阵,将该第一子矩阵作为第三噪声矩阵Unn1;用噪声子空间Unn的后L1行(L1+L2-K)列所有元素生成L1×(L1+L2-K)维的第二子矩阵,将该第二子矩阵作为第四噪声矩阵Unn2
步骤10:根据第一噪声矩阵Un1,第二噪声矩阵Un2,第三噪声矩阵Unn1和第四噪声矩阵Unn2构造八个噪声向量c1,c2,c3,c4,c5,c6,c7和c8
(10a)根据第一噪声矩阵Un1和第三噪声矩阵Unn1,分别计算第一噪声向量c1和第二噪声向量c2
c1=[c1(1),c1(2),…,c1(u),...,c1(2L1-1)],
c2=[c2(1),c2(2),...,c2(r),…,c2(2L2-1)],
其中,c1(u)表示第一噪声向量c1中的第u个元素,u的取值范围是1≤u≤2L1-1,的取值范围是:c2(r)表示第二噪声向量c2中的第r个元素,r的取值范围是1≤r≤2L2-1,的取值范围是:
(10b)根据第一噪声矩阵Un1和第三噪声矩阵Unn1,以及第二噪声矩阵Un2和第四噪声矩阵Unn2,分别计算第三噪声向量c3和第四噪声向量c4
c3=[c3(1),c3(2),…,c3(v),…,c3(L1+L2-1)],
c4=[c4(1),c4(2),…,c4(m),…,c4(L1+L2-1)],
其中,c3(v)表示第三噪声向量c3中第v个元素,v的取值范围是1≤v≤L1+L2-1,的取值范围是c4(m)表示第四噪声向量c4中第m个元素,m的取值范围是1≤m≤L1+L2-1,的取值范围是
(10c)根据第一噪声矩阵Un1和第三噪声矩阵Unn1,以及第二噪声矩阵Un2和第四噪声矩阵Unn2,分别计算第五噪声向量c5和第六噪声向量c6
c5=[c5(1),c5(2),…,c5(w),…,c5(L1+L2-1)],
c6=[c6(1),c6(2),…,c6(f),…,c6(L1+L2-1)],
其中,c5(w)表示第五噪声向量c5中第w个元素,w的取值范围是1≤w≤L1+L2-1,的取值范围是c6(f)表示第六噪声向量c6中第f个元素,f的取值范围是1≤f≤L1+L2-1,的取值范围是
(10d)根据第一噪声矩阵Un1和第三噪声矩阵Unn1,分别计算第七噪声向量c7和第八噪声向量c8
c7=[c7(1),c7(2),...,c7(z),...,c7(2L2-1)],
c8=[c8(1),c8(2),…,c8(h),…,c8(2L2-1)],
其中,c7(z)表示第七噪声向量c7中的第z个元素,z的取值范围是1≤z≤2L2-1,的取值范围是c8(h)表示第八噪声向量c8中的第h个元素,h的取值范围是1≤h≤2L2-1,的取值范围是:
步骤11:根据c1,c2,c3,c4,c5,c6,c7和c8这八个噪声向量,构造第一复合向量p14和第一最终向量q23,构造多项式方程。
(11a)根据第二噪声向量c2和第七噪声向量c7构造两个暂用向量e1和e2
在第二噪声向量c2的前和后分别补(L1-L2)个0,成为第一暂用向量e1,其长度为2L1-1;
在第七噪声向量c7的前和后分别补(L1-L2)个0,成为第二暂用向量e2,其长度为2L1-1;
(11b)根据六个噪声向量c1,c3,c4,c5,c6,c8和两个暂用向量e1和e2,构造四个过渡向量m1,m2,m3和m4
根据第一噪声向量c1和第一暂用向量e1,计算第一过渡向量m1=c1+e1
根据第三噪声向量c3和第四噪声向量c4,计算第二过渡向量m2=c3+c4
根据第五噪声向量c5和第六噪声向量c6,计算第三过渡向量m3=c5+c6
根据第八噪声向量c8和第二暂用向量e2,计算第四过渡向量m4=c8+e2
(11c)根据第一过渡向量m1和第四过渡向量m4,计算第一复合向量p14
p14=[p14(1),p14(2),…,p14(x),…,p14(4L1-3)],
其中,p14(x)表示第一复合向量p14中的第x个元素,x的取值范围是1≤x≤4L1-3,的取值范围是
(11d)根据第二过渡向量m2和第三过渡向量m3,计算第二复合向量p23
p23=[p23(1),p23(2),...,p23(y),...,p23(2L1+2L2-3)],
其中,p23(y)表示第二复合向量p23中第y个元素,y的取值范围是1≤y≤2L1+2L2-3,的取值范围是
(11e)在第二复合向量p23的前后各补(L1-L2)个0,成为第一最终向量q23,将第一复合向量p14和第一最终向量q23中的元素做为多项式方程的系数,得到多项式方程:
步骤12:获得目标波达方向角度值θ。
(11a)根据多项式方程,计算多项式方程的所有根z:
该多项式最多有2L=3L1+L2-4个根,其中每个根都有一个与其相似的根,每对相似根中只保留其中一个,就得到了最多有L个根z1,...zn,....zL,如果信号数K<L,则此处得到根的数量应该是K个,分别为z1,…,zh,…,zK,将其表示为:
z=[z1,…,zh,…,zK],
其中,z表示多项式方程的根,zh表示多项式方程的第h个根,h的取值范围是1≤h≤K。
(11b)由多项式方程的每一个根的辐角与相应的目标波达方向角度值的关系,得到相应的目标波达方向角度值:
θh=arcsin(λ/(2πd)arg(zh)),
其中,θh表示第h个目标信号波达方向角度值;
(11c)由每一个的目标波达方向角度值,得到目标波达方向角度值θ:
θ=[θ12,…,θh,…,θK],
完成对目标波达方向角度值θ的估计。
综上,本发明解决了现有技术无法利用嵌套阵列求解非圆信号DOA角度的问题,均匀阵列识别信源数目少,或没有利用非圆信号特性等问题,降低了对阵元数目的要求,保证了阵元数目使用的高效性,提高了一定阵元数情况下阵列可识别的信源数目以及低信噪比下对非圆信号方向角的估计性能。

Claims (8)

1.一种基于多项式求解的非圆信号波达方向角估计方法,其特征在于,包括:
(1)用M+N个天线接收机形成嵌套阵列,其中M、N分别表示两个天线接收阵列的阵元数,其取值范围为M≥1,N≥1;
(2)假设空间中有K个非圆目标信号,通过采样和滤波得到嵌套阵列输出信号:Y(t)=[y1(t),…,yi(t),…,yM+N(t)]T,其中,yi(t)表示嵌套阵列的第i个阵元的输出信号,t的取值范围是1≤t≤L,L表示快拍数,i的取值范围是1≤i≤M+N,(·)T表示矩阵转置运算;
(3)利用嵌套阵列输出信号Y(t),计算虚拟均匀阵列协方差向量和椭圆协方差向量
(4)根据虚拟均匀阵列协方差向量得到波达角选择矩阵G,根据虚拟均匀阵列椭圆协方差向量得到波达角选择矩阵Q;
(5)利用矩阵特征值分解的方法,分别计算波达角选择矩阵G的噪声子空间Un和波达角选择矩阵Q的噪声子空间Unn
(6)根据波达角选择矩阵G的噪声子空间Un和波达角选择矩阵Q的噪声子空间Unn,获得不同的噪声矩阵:
(6a)提取波达角选择矩阵G的噪声子空间Un的前L1行(L1+L2-K)列的所有元素,生成L1×(L1+L2-K)维的第一子矩阵,将该第一子矩阵作为第一噪声矩阵Un1;用噪声子空间Un的后L2行(L1+L2-K)列所有元素生成L2×(L1+L2-K)维的第二子矩阵,将该第二子矩阵作为第二噪声矩阵Un2
(6b)提取波达角选择矩阵Q的噪声子空间Unn的前L2行(L1+L2-K)列所有元素生成L2×(L1+L2-K)维的第一子矩阵,将该第一子矩阵作为第三噪声矩阵Unn1;用噪声子空间Unn的后L1行(L1+L2-K)列所有元素生成L1×(L1+L2-K)维的第二子矩阵,将该第二子矩阵作为第四噪声矩阵Unn2
(7)根据第一噪声矩阵Un1,第二噪声矩阵Un2,第三噪声矩阵Unn1和第四噪声矩阵Unn2构造八个噪声向量c1,c2,c3,c4,c5,c6,c7和c8
(8)根据c1,c2,c3,c4,c5,c6,c7和c8这八个噪声向量,构造第一复合向量p14和第一最终向量q23,得到如下多项式方程:
其中,p14(j)表示第一复合向量p14中第j个元素,q23(j)表示第一最终向量q23中第j个元素,j的取值范围是1≤j≤2(L1+L2)-3,z表示多项式方程的根,z=[z1,…,zh,…,zK],zh表示多项式方程的第h个根,该多项式最多有2L=3L1+L2-4个根,且每个根都有一个与其相似的根,每对相似根中只保留其中一个,得到最多有L个根z1,...zn,....zL
(9)计算多项式方程的所有根,由多项式方程的每一个根的辐角与目标波达方向角度值的关系,得到目标波达方向角度值θ。
2.根据权利要求1所述的方法,其中步骤(1)中用M+N个天线接收机形成嵌套阵列,按如下步骤进行:
(1a)将每个天线接收机称为一个阵元,用M个天线接收机形成第一均匀线性阵列a,其阵元间距为d,定义第一均匀线性阵列a的第一个阵元为起始阵元,定义起始阵元位置D(1)=1,第一均匀线性阵列a的其它阵元位置依次为D(2)=2,D(3)=3,D(4)=4,…,D(M)=M;
(1b)用N个天线接收机形成第二均匀线性阵列b,其阵元间距为(M+1)d,第二均匀线性阵列b的阵元位置依次设置为D(M+1)=M+1,D(M+2)=2(M+1),D(M+2)=3(M+1),…,D(M+N)=N(M+1),其中,M≥1,N≥1,0<d≤λ/2,λ为入射到阵列的窄带信号波长;
(1c)将第二均匀线性阵列b的第一个阵元放置于与起始阵元相距为Md的位置,将第二均匀线性阵列b的所有阵元依次放置在第一均匀线性阵列a之后,形成嵌套阵列。
3.根据权利要求1所述的方法,其中步骤(3)中计算虚拟均匀阵列协方差向量和椭圆协方差向量按如下步骤进行:
(3a)根据输出信号Y(t)计算协方差矩阵Rd和椭圆协方差矩阵Rs
其中,(·)H表示共轭转置运算;
(3b)将等效协方差矩阵Rd和椭圆协方差矩阵Rs中的元素分别进行排列,得到等效协方差向量rd和等效椭圆协方差向量rs
rd=[Rd(1,1),Rd(2,1),…,Rd(M+N,1),Rd(1,2),…,Rd(M+N,2),…,Rd(1,M+N),…,Rd(M+N,M+N)]T
rs=[Rs(1,1),Rs(2,1),…,Rs(M+N,1),Rs(1,2),…,Rs(M+N,2),…,Rs(1,M+N),…,Rs(M+N,M+N)]T
其中,Rd(i,j)表示协方差矩阵Rd中位于第i行,第j列的元素,i的取值范围为1≤i≤M+N,j的取值范围为1≤j≤M+N;Rs(i,j)表示椭圆协方差矩阵Rs中位于第i行,第j列的元素;
(3c)计算等效协方差向量rd中所有元素的维数Ei,j和等效椭圆协方差向量rs中所有元素的维数Fi,j
Ei,j=D(j)-D(i)
Fi,j=D(j)+D(i)
其中,D(i)表示嵌套阵列中第i个阵元的位置,D(j)表示嵌套阵列中第j个阵元的位置;
(3d)删除等效协方差向量rd中维数相同的元素和维数不连续的元素,并将剩余元素按维数从小到大排列,得到虚拟均匀阵列协方差向量删除等效椭圆协方差向量rs中维数相同的元素和维数不连续的元素,并将剩余元素按维数从小到大排列,得到虚拟均匀阵列椭圆协方差向量其中中元素个数为Cd=2MN+2N-1,中元素个数为Cs=MN+M+N。
4.根据权利要求1所述的方法,其中步骤(4)中的波达角选择矩阵G和波达角选择矩阵Q的构造方法如下:
其中,中间变量L1=(Cd+1)/2,L2=Cs+1-(Cd+1)/2,且有L1>L2,Cd表示虚拟均匀阵列协方差向量中元素的个数,Cd的取值为2MN+2N-1,Cs表示虚拟均匀阵列椭圆协方差向量中元素的个数,Cs的取值为MN+M+N,(·)*表示向量的共轭运算。
5.根据权利要求1所述的方法,其中步骤(5)中利用矩阵特征值分解的方法,计算波达角选择矩阵G的噪声子空间Un,Q的噪声子空间Unn,按如下步骤进行:
(5a)将波达角选择矩阵G进行如下特征分解:
G=U·Λ·UH
其中,Λ为波达角选择矩阵G的特征值矩阵,U为矩阵G的特征值所对应的特征向量矩阵,(·)H表示矩阵的共轭转置运算;
(5b)将特征值矩阵Λ中的特征值按从大到小排序,取后(L1+L2-K)个较小特征值对应的特征向量矩阵作为噪声子空间Un
(5c)将波达角选择矩阵Q进行如下特征分解:
Q=U·Λ·UH
其中,Λ为波达角选择矩阵Q的特征值矩阵,U为矩阵Q的特征值所对应的特征向量矩阵,(·)H表示矩阵的共轭转置运算;
(5d)将特征值矩阵Λ中的特征值按从大到小排序,取后(L1+L2-K)个较小特征值对应的特征向量矩阵作为噪声子空间Unn
6.根据权利要求1所述的方法,其中步骤(7)中根据第一噪声矩阵Un1,第二噪声矩阵Un2,第三噪声矩阵Unn1和第四噪声矩阵Unn2构造八个噪声向量c1,c2,c3,c4,c5,c6,c7和c8,按如下步骤进行:
(7a)根据第一噪声矩阵Un1和第三噪声矩阵Unn1,分别计算第一噪声向量c1和第二噪声向量c2
c1=[c1(1),c1(2),…,c1(u),…,c1(2L1-1)]
c2=[c2(1),c2(2),…,c2(r),…,c2(2L2-1)]
其中,c1(u)表示第一噪声向量c1中的第u个元素,u的取值范围是1≤u≤2L1-1,的取值范围是:c2(r)表示第二噪声向量c2中的第r个元素,r的取值范围是1≤r≤2L2-1,的取值范围是:
(7b)根据第一噪声矩阵Un1和第三噪声矩阵Unn1,以及第二噪声矩阵Un2和第四噪声矩阵Unn2,分别计算第三噪声向量c3和第四噪声向量c4
c3=[c3(1),c3(2),…,c3(v),…,c3(L1+L2-1)],
c4=[c4(1),c4(2),…,c4(m),…,c4(L1+L2-1)],
其中,c3(v)表示第三噪声向量c3中第v个元素,v的取值范围是1≤v≤L1+L2-1,的取值范围是c4(m)表示第四噪声向量c4中第m个元素,m的取值范围是1≤m≤L1+L2-1,的取值范围是
(7c)根据第一噪声矩阵Un1和第三噪声矩阵Unn1,以及第二噪声矩阵Un2和第四噪声矩阵Unn2,分别计算第五噪声向量c5和第六噪声向量c6
c5=[c5(1),c5(2),…,c5(w),…,c5(L1+L2-1)]
c6=[c6(1),c6(2),…,c6(f),…,c6(L1+L2-1)]
其中,c5(w)表示第五噪声向量c5中第w个元素,w的取值范围是1≤w≤L1+L2-1,的取值范围是c6(f)表示第六噪声向量c6中第f个元素,f的取值范围是1≤f≤L1+L2-1,的取值范围是
(7d)根据第一噪声矩阵Un1和第三噪声矩阵Unn1,分别计算第七噪声向量c7和第八噪声向量c8
c7=[c7(1),c7(2),…,c7(z),…,c7(2L2-1)]
c8=[c8(1),c8(2),…,c8(h),…,c8(2L2-1)]
其中,c7(z)表示第七噪声向量c7中的第z个元素,z的取值范围是1≤z≤2L2-1,的取值范围是c8(h)表示第八噪声向量c8中的第h个元素,h的取值范围是1≤h≤2L2-1,的取值范围是:
7.根据权利要求1所述的方法,其中步骤(8)中根据八个噪声向量c1,c2,c3,c4,c5,c6,c7和c8构造第一复合向量p14和第一最终向量q23,得到多项式方程,按如下步骤进行:
(8a)根据第二噪声向量c2和第七噪声向量c7构造两个暂用向量e1和e2
在第二噪声向量c2的前和后分别补(L1-L2)个0,成为第一暂用向量e1,其长度为2L1-1;
在第七噪声向量c7的前和后分别补(L1-L2)个0,成为第二暂用向量e2,其长度为2L1-1;
(8b)根据六个噪声向量c1,c3,c4,c5,c6,c8和两个暂用向量e1和e2,构造四个过渡向量m1,m2,m3和m4
根据第一噪声向量c1和第一暂用向量e1,计算第一过渡向量m1=c1+e1
根据第三噪声向量c3和第四噪声向量c4,计算第二过渡向量m2=c3+c4
根据第五噪声向量c5和第六噪声向量c6,计算第三过渡向量m3=c5+c6
根据第八噪声向量c8和第二暂用向量e2,计算第四过渡向量m4=c8+e2
(8c)根据第一过渡向量m1和第四过渡向量m4,计算第一复合向量p14
p14=[p14(1),p14(2),…,p14(x),…,p14(4L1-3)],
其中,p14(x)表示第一复合向量p14中的第x个元素,x的取值范围是1≤x≤4L1-3,的取值范围是
(8d)根据第二过渡向量m2和第三过渡向量m3,计算第二复合向量p23
p23=[p23(1),p23(2),…,p23(y),…,p23(2L1+2L2-3)],
其中,p23(y)表示第二复合向量p23中第y个元素,y的取值范围是1≤y≤2L1+2L2-3,的取值范围是
(8e)在第二复合向量p23前后各补(L1-L2)个0,成为第一最终向量q23,将第一复合向量p14和第一最终向量q23中的元素做为多项式方程的系数,得到多项式方程:
8.根据权利要求1所述的方法,其中步骤(9)中由多项式方程的每一个根的辐角与目标波达方向角度值的关系,得到目标波达方向角度值θ,按如下步骤进行:
(9a)计算多项式方程的所有根z:
如果非圆目标信号数K<L,则此处得到根的数量为K个,分别为z1,...zn,....zk,即:
z=[z1,…,zh,…,zK]
(9b)由多项式方程的每一个根的辐角与相应的目标波达方向角度值的关系,得到相应的目标波达方向角度值:
θh=arcsin(λ/(2πd)arg(zh)),
其中,θh表示第h个目标信号波达方向角度值,λ表示入射到阵列的窄带信号波长,d表示第一均匀线性阵列a的阵元间距;
(9c)由每一个的目标波达方向角度值,得到目标波达方向角度值θ:
θ=[θ12,…,θh,…,θK]。
CN201710423912.3A 2017-06-07 2017-06-07 基于多项式求解的非圆信号波达方向角的估计方法 Active CN107300686B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710423912.3A CN107300686B (zh) 2017-06-07 2017-06-07 基于多项式求解的非圆信号波达方向角的估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710423912.3A CN107300686B (zh) 2017-06-07 2017-06-07 基于多项式求解的非圆信号波达方向角的估计方法

Publications (2)

Publication Number Publication Date
CN107300686A true CN107300686A (zh) 2017-10-27
CN107300686B CN107300686B (zh) 2019-08-06

Family

ID=60135405

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710423912.3A Active CN107300686B (zh) 2017-06-07 2017-06-07 基于多项式求解的非圆信号波达方向角的估计方法

Country Status (1)

Country Link
CN (1) CN107300686B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107870315A (zh) * 2017-11-06 2018-04-03 重庆邮电大学 一种利用迭代相位补偿技术估计任意阵列波达方向方法
CN108490383A (zh) * 2018-03-07 2018-09-04 大连理工大学 一种基于有界非线性协方差的非圆信号波达方向估计方法
CN109490820A (zh) * 2018-11-13 2019-03-19 电子科技大学 一种基于平行嵌套阵的二维doa估计方法
CN113740802A (zh) * 2021-07-30 2021-12-03 香港中文大学(深圳) 以自适应噪声估计进行矩阵补全的信号源定位方法和系统

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103353588A (zh) * 2013-06-13 2013-10-16 西安电子科技大学 基于天线均匀平面阵的二维波达方向角估计方法
CN104020439A (zh) * 2014-06-20 2014-09-03 西安电子科技大学 基于空间平滑协方差矩阵稀疏表示的波达方向角估计方法
CN106443574A (zh) * 2016-11-08 2017-02-22 西安电子科技大学 基于双层嵌套阵列的波达方向角估计方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103353588A (zh) * 2013-06-13 2013-10-16 西安电子科技大学 基于天线均匀平面阵的二维波达方向角估计方法
CN104020439A (zh) * 2014-06-20 2014-09-03 西安电子科技大学 基于空间平滑协方差矩阵稀疏表示的波达方向角估计方法
CN106443574A (zh) * 2016-11-08 2017-02-22 西安电子科技大学 基于双层嵌套阵列的波达方向角估计方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
PASCAL CHARGE ET AL.: "A non-circular sources direction finding method using polynomial rooting", 《SIGNAL PROCESSING》 *
蔡晶晶 等: "基于空域平滑稀疏重构的DOA估计算法", 《电子与信息学报》 *

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107870315A (zh) * 2017-11-06 2018-04-03 重庆邮电大学 一种利用迭代相位补偿技术估计任意阵列波达方向方法
CN107870315B (zh) * 2017-11-06 2021-07-30 重庆邮电大学 一种利用迭代相位补偿技术估计任意阵列波达方向方法
CN108490383A (zh) * 2018-03-07 2018-09-04 大连理工大学 一种基于有界非线性协方差的非圆信号波达方向估计方法
CN109490820A (zh) * 2018-11-13 2019-03-19 电子科技大学 一种基于平行嵌套阵的二维doa估计方法
CN109490820B (zh) * 2018-11-13 2021-04-27 电子科技大学 一种基于平行嵌套阵的二维doa估计方法
CN113740802A (zh) * 2021-07-30 2021-12-03 香港中文大学(深圳) 以自适应噪声估计进行矩阵补全的信号源定位方法和系统
CN113740802B (zh) * 2021-07-30 2022-07-22 香港中文大学(深圳) 以自适应噪声估计进行矩阵补全的信号源定位方法和系统

Also Published As

Publication number Publication date
CN107300686B (zh) 2019-08-06

Similar Documents

Publication Publication Date Title
CN107037393B (zh) 基于嵌套阵列的非圆信号波达方向角估计方法
CN105182293B (zh) 基于互质阵列mimo雷达doa与dod估计方法
Wen et al. Angle estimation and mutual coupling self-calibration for ULA-based bistatic MIMO radar
CN107300686B (zh) 基于多项式求解的非圆信号波达方向角的估计方法
CN104931931B (zh) 互耦条件下基于张量实值子空间的双基地mimo雷达角度估计方法
CN104991236B (zh) 一种单基地mimo雷达非圆信号相干源波达方向估计方法
CN105259550B (zh) 基于压缩感知的多输入多输出雷达二维角度估计方法
CN110463147A (zh) 解码符号的方法以及接收并解码符号的接收器
CN103245956B (zh) 一种基于稳健波束形成算法的gps抗多径方法
CN105403874B (zh) 非均匀阵列欠定波达方向估计方法
CN106569171B (zh) 基于双层混合阵列的波达方向角估计方法
CN103901395B (zh) 一种冲击噪声环境下相干信号波达方向动态跟踪方法
CN106443574A (zh) 基于双层嵌套阵列的波达方向角估计方法
CN106646344A (zh) 一种利用互质阵的波达方向估计方法
CN106772226A (zh) 基于压缩感知时间调制阵列的doa估计方法
CN106785486B (zh) 一种广义互质面阵天线结构
CN103344940B (zh) 低复杂度的doa估计方法及系统
CN107589399A (zh) 基于多采样虚拟信号奇异值分解的互质阵列波达方向估计方法
CN108020812B (zh) 基于特殊三平行线阵结构的二维doa估计方法
CN106019215A (zh) 基于四阶累量的嵌套阵列波达方向角估计方法
CN104020439A (zh) 基于空间平滑协方差矩阵稀疏表示的波达方向角估计方法
CN106972882A (zh) 基于虚拟域空间功率谱估计的互质阵列自适应波束成形方法
CN103983952A (zh) 一种非圆信号双基地mimo雷达低复杂度收发角度联合估计方法
CN108896954A (zh) 互质阵中一种基于联合实值子空间的波达角估计方法
CN107315161A (zh) 基于压缩感知的非圆信号波达方向角估计方法

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