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

摘要

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

说明书

技术领域

本发明属于信号处理技术领域，特别涉及一种电磁信号的阵列信号波达方向角估计方法，可用于对飞机、舰船运动目标的侦察与无源定位。

背景技术

信号的波达方向角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)＝[y₁(t),…,y_i(t),…,y_M+N(t)]^T，其中，y_i(t)表示嵌套阵列的第i个阵元的输出信号，t的取值范围是1≤t≤L，L表示快拍数，i的取值范围是1≤i≤M+N，(·)^T表示矩阵转置运算；

(3)利用嵌套阵列输出信号Y(t)，计算虚拟均匀阵列协方差向量和椭圆协方差向量

(4)根据虚拟均匀阵列协方差向量得到波达角选择矩阵G，根据虚拟均匀阵列椭圆协方差向量得到波达角选择矩阵Q；

(5)分别计算波达角选择矩阵G的噪声子空间U_n和波达角选择矩阵Q的噪声子空间U_nn；

(6)根据波达角选择矩阵G的噪声子空间U_n和波达角选择矩阵Q的噪声子空间U_nn，获得不同的噪声矩阵：

(6a)提取波达角选择矩阵G的噪声子空间U_n的前L₁行(L₁+L₂-K)列的所有元素，生成L₁×(L₁+L₂-K)维的第一子矩阵，将该第一子矩阵作为第一噪声矩阵U_n1；用噪声子空间U_n的后L₂行(L₁+L₂-K)列所有元素生成L₂×(L₁+L₂-K)维的第二子矩阵，将该第二子矩阵作为第二噪声矩阵U_n2；

(6b)提取波达角选择矩阵Q的噪声子空间U_nn的前L₂行(L₁+L₂-K)列所有元素生成L₂×(L₁+L₂-K)维的第一子矩阵，将该第一子矩阵作为第三噪声矩阵U_nn1；用噪声子空间U_nn的后L₁行(L₁+L₂-K)列所有元素生成L₁×(L₁+L₂-K)维的第二子矩阵，将该第二子矩阵作为第四噪声矩阵U_nn2；

(7)根据第一噪声矩阵U_n1，第二噪声矩阵U_n2，第三噪声矩阵U_nn1和第四噪声矩阵U_nn2构造八个噪声向量c₁，c₂，c₃，c₄，c₅，c₆，c₇和c₈；

(8)根据c₁，c₂，c₃，c₄，c₅，c₆，c₇和c₈这八个噪声向量，构造第一复合向量p₁₄和第一最终向量q₂₃，得到如下多项式方程：

其中，p₁₄(j)表示第一复合向量p₁₄中第j个元素，q₂₃(j)表示第一最终向量q₂₃中第j个元素，j的取值范围是1≤j≤2(L₁+L₂)-3，z表示多项式方程的根，z＝[z₁,…,z_h,…,z_K]，z_h表示多项式方程的第h个根，该多项式最多有2L＝3L₁+L₂-4个根，且每个根都有一个与其相似的根，每对相似根中只保留其中一个，得到最多有L个根z₁，…z_n,….z_L；

(9)计算多项式方程的所有根，由多项式方程的每一个根的辐角与目标波达方向角度值的关系，得到目标波达方向角度值θ。

本发明与现有技术相比具有以下优点：

1)本发明采用了嵌套阵列模型进行波达方向角度估计，克服了现有技术中采用典型的线性均匀阵列，造成估计的信号数目低于阵元数目的缺点，提高了在阵元数目相同的条件下的阵列可识别信源数目。

2)本发明利用信号的协方差矩阵R_d，以及椭圆协方差矩阵R_s对信号进行估计，增多了可估计的非圆信号数目。

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)＝[y₁(t),…,y_i(t),…,y_M+N(t)]^T，其中，K的取值范围是K＜MN+M+N-1，y_i(t)表示嵌套阵列第i个阵元的输出信号，t的取值范围是1≤t≤L，L表示快拍数，i的取值范围是1≤i≤M+N，(·)^T表示矩阵转置运算。

步骤3：计算协方差矩阵R_d和椭圆协方差矩阵R_s。

利用嵌套阵列输出信号Y(t)，计算协方差矩阵R_d和椭圆协方差矩阵R_s：

其中，(·)^H表示矩阵共轭转置运算。

步骤4：构造等效协方差向量r_d和等效椭圆协方差向量r_s。

根据协方差矩阵R_d和椭圆协方差矩阵R_s中的元素，分别构造等效协方差向量r_d和等效椭圆协方差向量r_s：

r_d＝[R_d(1,1),R_d(2,1),…,R_d(M+N,1),R_d(1,2),…,R_d(M+N,2),

...,R_d(i,j),…,R_d(1,M+N),…,R_d(M+N,M+N)]^T

r_s＝[R_s(1,1),R_s(2,1),...,R_s(M+N,1),R_s(1,2),…,R_s(M+N,2),

…,R_s(i,j),…,R_s(1,M+N),…,R_s(M+N,M+N)]^T

其中，R_d(i,j)表示协方差矩阵R_d中位于第i行、第j列的元素，i的取值范围为1≤i≤M+N，j的取值范围为1≤j≤M+N，R_s(i,j)表示椭圆协方差矩阵R_s中位于第i行，第j列的元素。

步骤5：计算等效协方差向量和等效椭圆协方差向量中所有元素的维数。

根据等效协方差向量r_d和等效椭圆协方差向量r_s中每一个元素所在的行和列在嵌套阵列中对应的阵元位置，计算等效协方差向量r_d中所有元素的维数E_i,j和等效椭圆协方差向量r_s中所有元素的维数F_i,j：

E_i,j＝D(j)-D(i)，

F_i,j＝D(j)+D(i)，

其中，D(i)表示嵌套阵列中第i个阵元的位置，D(j)表示嵌套阵列中第j个阵元的位置。

步骤6：获得虚拟均匀阵列协方差向量和虚拟均匀阵列椭圆协方差向量

根据等效协方差向量r_d中所有元素的维数，删除等效协方差向量r_d中维数相同的元素和维数不连续的元素，并将剩余元素按维数从小到大排列，得到虚拟均匀阵列协方差向量

根据等效椭圆协方差向量r_s中所有元素的维数，删除等效椭圆协方差向量r_s中维数相同的元素和维数不连续的元素，并将剩余元素按维数从小到大排列，得到虚拟均匀阵列椭圆协方差向量

步骤7：构造波达角选择矩阵G和波达角选择矩阵Q。

根据虚拟均匀阵列协方差向量和虚拟均匀阵列椭圆协方差向量中所有元素进行行列排列，得到波达角选择矩阵G和Q：

其中，中间变量L₁＝(C_d+1)/2,L₂＝C_s+1-(C_d+1)/2，且有L₁＞L₂，C_d表示虚拟均匀阵列协方差向量中元素的个数，C_d的取值为2MN+2N-1，C_s表示虚拟均匀阵列椭圆协方差向量中元素的个数，C_s的取值为MN+M+N，(·)^*表示向量的共轭运算。

步骤8：计算波达角选择矩阵G的噪声子空间U_n，波达角选择矩阵Q的噪声子空间U_nn。

(8a)将波达角选择矩阵G进行如下特征分解：

G＝U·Λ·U^H，

其中，Λ为波达角选择矩阵G的特征值矩阵，U为矩阵G的特征值所对应的特征向量矩阵，(·)^H表示矩阵的共轭转置运算；

(8b)将特征值矩阵Λ中的特征值按从大到小排序，取后(L₁+L₂-K)个较小特征值对应的特征向量矩阵作为噪声子空间U_n；

(8c)将波达角选择矩阵Q进行如下特征分解：

Q＝U·Λ·U^H，

其中，Λ为波达角选择矩阵Q的特征值矩阵，U为矩阵Q的特征值所对应的特征向量矩阵，(·)^H表示矩阵的共轭转置运算；

(8d)将特征值矩阵Λ中的特征值按从大到小排序，取后(L₁+L₂-K)个较小特征值对应的特征向量矩阵作为噪声子空间U_nn。

步骤9：根据G的噪声子空间U_n和波达角选择矩阵Q的噪声子空间U_nn，获得不同的噪声矩阵。

(9a)提取波达角选择矩阵G的噪声子空间U_n的前L₁行(L₁+L₂-K)列的所有元素，生成L₁×(L₁+L₂-K)维的第一子矩阵，将该第一子矩阵作为第一噪声矩阵U_n1；用噪声子空间U_n的后L₂行(L₁+L₂-K)列所有元素生成L₂×(L₁+L₂-K)维的第二子矩阵，将该第二子矩阵作为第二噪声矩阵U_n2；

(9b)提取波达角选择矩阵Q的噪声子空间U_nn的前L₂行(L₁+L₂-K)列所有元素生成L₂×(L₁+L₂-K)维的第一子矩阵，将该第一子矩阵作为第三噪声矩阵U_nn1；用噪声子空间U_nn的后L₁行(L₁+L₂-K)列所有元素生成L₁×(L₁+L₂-K)维的第二子矩阵，将该第二子矩阵作为第四噪声矩阵U_nn2。

步骤10：根据第一噪声矩阵U_n1，第二噪声矩阵U_n2，第三噪声矩阵U_nn1和第四噪声矩阵U_nn2构造八个噪声向量c₁，c₂，c₃，c₄，c₅，c₆，c₇和c₈。

(10a)根据第一噪声矩阵U_n1和第三噪声矩阵U_nn1，分别计算第一噪声向量c₁和第二噪声向量c₂：

c₁＝[c₁(1),c₁(2),…,c₁(u),...,c₁(2L₁-1)]，

c₂＝[c₂(1),c₂(2),...,c₂(r),…,c₂(2L₂-1)]，

其中，c₁(u)表示第一噪声向量c₁中的第u个元素，u的取值范围是1≤u≤2L₁-1，的取值范围是：c₂(r)表示第二噪声向量c₂中的第r个元素，r的取值范围是1≤r≤2L₂-1，的取值范围是：

(10b)根据第一噪声矩阵U_n1和第三噪声矩阵U_nn1，以及第二噪声矩阵U_n2和第四噪声矩阵U_nn2，分别计算第三噪声向量c₃和第四噪声向量c₄：

c₃＝[c₃(1),c₃(2),…,c₃(v),…,c₃(L₁+L₂-1)]，

c₄＝[c₄(1),c₄(2),…,c₄(m),…,c₄(L₁+L₂-1)]，

其中，c₃(v)表示第三噪声向量c₃中第v个元素，v的取值范围是1≤v≤L₁+L₂-1，的取值范围是c₄(m)表示第四噪声向量c₄中第m个元素，m的取值范围是1≤m≤L₁+L₂-1，的取值范围是

(10c)根据第一噪声矩阵U_n1和第三噪声矩阵U_nn1，以及第二噪声矩阵U_n2和第四噪声矩阵U_nn2，分别计算第五噪声向量c₅和第六噪声向量c₆：

c₅＝[c₅(1),c₅(2),…,c₅(w),…,c₅(L₁+L₂-1)]，

c₆＝[c₆(1),c₆(2),…,c₆(f),…,c₆(L₁+L₂-1)]，

其中，c₅(w)表示第五噪声向量c₅中第w个元素，w的取值范围是1≤w≤L₁+L₂-1，的取值范围是c₆(f)表示第六噪声向量c₆中第f个元素，f的取值范围是1≤f≤L₁+L₂-1，的取值范围是

(10d)根据第一噪声矩阵U_n1和第三噪声矩阵U_nn1，分别计算第七噪声向量c₇和第八噪声向量c₈：

c₇＝[c₇(1),c₇(2),...,c₇(z),...,c₇(2L₂-1)]，

c₈＝[c₈(1),c₈(2),…,c₈(h),…,c₈(2L₂-1)]，

其中，c₇(z)表示第七噪声向量c₇中的第z个元素，z的取值范围是1≤z≤2L₂-1，的取值范围是c₈(h)表示第八噪声向量c₈中的第h个元素，h的取值范围是1≤h≤2L₂-1，的取值范围是：

步骤11：根据c₁，c₂，c₃，c₄，c₅，c₆，c₇和c₈这八个噪声向量，构造第一复合向量p₁₄和第一最终向量q₂₃，构造多项式方程。

(11a)根据第二噪声向量c₂和第七噪声向量c₇构造两个暂用向量e₁和e₂：

在第二噪声向量c₂的前和后分别补(L₁-L₂)个0，成为第一暂用向量e₁，其长度为2L₁-1；

在第七噪声向量c₇的前和后分别补(L₁-L₂)个0，成为第二暂用向量e₂，其长度为2L₁-1；

(11b)根据六个噪声向量c₁，c₃，c₄，c₅，c₆，c₈和两个暂用向量e₁和e₂，构造四个过渡向量m₁，m₂，m₃和m₄：

根据第一噪声向量c₁和第一暂用向量e₁，计算第一过渡向量m₁＝c₁+e₁；

根据第三噪声向量c₃和第四噪声向量c₄，计算第二过渡向量m₂＝c₃+c₄；

根据第五噪声向量c₅和第六噪声向量c₆，计算第三过渡向量m₃＝c₅+c₆；

根据第八噪声向量c₈和第二暂用向量e₂，计算第四过渡向量m₄＝c₈+e₂；

(11c)根据第一过渡向量m₁和第四过渡向量m₄，计算第一复合向量p₁₄：

p₁₄＝[p₁₄(1),p₁₄(2),…,p₁₄(x),…,p₁₄(4L₁-3)]，

其中，p₁₄(x)表示第一复合向量p₁₄中的第x个元素，x的取值范围是1≤x≤4L₁-3，的取值范围是

(11d)根据第二过渡向量m₂和第三过渡向量m₃，计算第二复合向量p₂₃：

p₂₃＝[p₂₃(1),p₂₃(2),...,p₂₃(y),...,p₂₃(2L₁+2L₂-3)]，

其中，p₂₃(y)表示第二复合向量p₂₃中第y个元素，y的取值范围是1≤y≤2L₁+2L₂-3，的取值范围是

(11e)在第二复合向量p₂₃的前后各补(L₁-L₂)个0，成为第一最终向量q₂₃，将第一复合向量p₁₄和第一最终向量q₂₃中的元素做为多项式方程的系数，得到多项式方程：

步骤12：获得目标波达方向角度值θ。

(11a)根据多项式方程，计算多项式方程的所有根z：

该多项式最多有2L＝3L₁+L₂-4个根，其中每个根都有一个与其相似的根，每对相似根中只保留其中一个，就得到了最多有L个根z₁，...z_n,....z_L，如果信号数K＜L，则此处得到根的数量应该是K个，分别为z₁,…,z_h,…,z_K，将其表示为：

z＝[z₁,…,z_h,…,z_K]，

其中，z表示多项式方程的根，z_h表示多项式方程的第h个根，h的取值范围是1≤h≤K。

(11b)由多项式方程的每一个根的辐角与相应的目标波达方向角度值的关系，得到相应的目标波达方向角度值：

θ_h＝arcsin(λ/(2πd)arg(z_h))，

其中，θ_h表示第h个目标信号波达方向角度值；

(11c)由每一个的目标波达方向角度值，得到目标波达方向角度值θ：

θ＝[θ₁,θ₂,…,θ_h,…,θ_K]，

完成对目标波达方向角度值θ的估计。

综上，本发明解决了现有技术无法利用嵌套阵列求解非圆信号DOA角度的问题，均匀阵列识别信源数目少，或没有利用非圆信号特性等问题，降低了对阵元数目的要求，保证了阵元数目使用的高效性，提高了一定阵元数情况下阵列可识别的信源数目以及低信噪比下对非圆信号方向角的估计性能。