本发明属于信号处理技术领域，涉及一种信号DOA估计方法，具体涉及一种基于双层平行圆阵列天线的DOA估计方法，可用于信号DOA估计处理系统。
背景技术：
波达方向(DOA)是信号源到达阵列天线的方向，该方向可以由方位角和俯仰角来表示。DOA估计即是在某个空间阵列接收模型下对来波方向的方位角和俯仰角进行估计。空间阵列接收模型包括线阵和面阵。基于线阵的DOA估计只可估计[-π/2，π/2]的方位信息，基于面阵的DOA估计系统能估计[-π，π]的方位角和[0，π/2]的俯仰角信息，作为面阵中的一种，均匀圆面阵阵列结构具有对称性，在平面旋转时波束鲁棒性好，因此基于均匀圆阵的信号DOA估计一直是阵列信号处理领域的研究热点。在阵列信号处理中，DOA估计精度越高，对干扰的抑制能力就越强，信道容量就越大。传统的DOA估计算法包括：Bartlett波束形成法、Capon最小方差法和多重信号分类(MUSIC)法。Bartlett算法在只有一个信号时是可行的，但当存在多个信号源入射时性能估计会急剧下降，同时，该方法受到旁瓣高度和波束宽度的限制，因此空间分辨率比较低。Capon算法能提供更佳的空间分辨率，但是需要扫描整个空间，寻找空间谱的峰值，同时还需要求解接收信号的协方差矩阵的逆矩阵，计算量非常大，而且若空间中存在与期望信号相关的其他干扰信号，相关的干扰信号会恶性与期望信号合并，导致估计性能急剧下降。
MUSIC算法利用阵列模型计算信号的导向矢量矩阵，利用阵列模型接收的样本数据计算数据协方差矩阵，利用数据协方差矩阵计算噪声子空间的特征向量矩阵，利用信号的导向矢量矩阵和噪声子空间的特征向量矩阵构造空间谱的估计式，通过搜索空间谱的谱峰值，可以找到谱峰值对应的角度，该角度即为波达方向角度，MUSIC算法能够同时对多个信号进行DOA估计并且具有较高的分辨率，但它也需要进行空间谱的谱峰搜索，计算量很大，而且当信源高度相关或者入射角过于接近时，会使计算所得的包含角度信息的噪声子空间的特征向量不准确，导致空间谱估计式的分辨率低，空间谱估计式的分辨率低会使DOA估计精度低。针对MUSIC算法存在的不足，出现了一种改进型的方法，例如授权公告号为CN 103546221 B、授权公告日为2015.09.16，名称为“一种宽带相干信号波达角估计方法”的发明专利，利用变换矩阵对不同频点处的协方差矩阵进行变换，构造导向矩阵协方差矩阵；利用导向协方差矩阵估计传播算子和噪声子空间，并进行噪声子空间正交化；采用MUSIC算法计算空间谱并进行波达角估计。该方法针对传统MUSIC算法中信源高度相关时估计精度不高的问题，在不影响非相干信源DOA估计的基础上，提高了对相干信源的估计性能。该方法也存在不足之处：利用中间参数的特征向量，通过谱峰搜索来估计信号的波达方向，对空间谱谱峰搜索增加了运算量，而且因计算过程中的得到的特征向量并不能直接对应源信号的方向矢量，所以当多个信号的入射角过于接近时，该方法的分辨率较低，会使DOA估计精度不高。
技术实现要素：
本发明的目的在于克服上述现有技术存在的不足，提出了一种基于双层平行圆阵列天线的DOA估计方法，采用阵元在圆周上均匀分布的双层平行圆阵列为信号的接收模型，构造包含源信号空间基底的矩阵，对该矩阵进行特征分解可得到特征值和特征值对应的特征向量，利用该特征值对应源信号的方向元素和方向矢量，可直接求解出源信号的俯仰角和方位角，用于解决现有技术中存在的多个信号的入射角过于接近时DOA估计精度低和因谱峰搜索引起的运算量大的问题。
为实现上述目的，本发明采用的技术方案包括如下步骤：
(1)获取双层平行圆阵列天线的相关参数，包括上层子圆阵的阵元数N和半径R；获取信号源的相关参数，包括信号矢量S(t)、信号源个数M、第i个信号源的波达方向方位角和俯仰角和信号源频率f；获取噪声相关参数，包括上层子圆阵的噪声矢量N1(t)和下层子圆阵的噪声矢量N2(t)；获取双层平行圆阵列天线轴向虚拟平移操作的相关参数，包括双层平行圆阵列天线轴向虚拟平移次数p、每次平移距离d和轴向平移总次数Q；
(2)以上层子圆阵的圆心为参考点，计算上层子圆阵的导向矢量矩阵实现步骤为：
(2a)利用双层平行圆阵列天线的相关参数和信号源的相关参数，计算第i个信号源相对于参考点在上层子圆阵第k个阵元的时延τki：
其中，c表示光速；
(2b)利用第i个信号源相对于参考点在上层子圆阵第k个阵元的时延τki，计算第i个信号源的方向矢量
其中，j表示虚数单位，[·]T表示·的转置；
(2c)利用第i个信号源的方向矢量计算上层子圆阵的导向矢量矩阵
(3)利用双层平行圆阵列天线轴向虚拟平移次数p和每次平移距离d，计算第p次平移后的自协方差矩阵Rxx(p)和互协方差矩阵Ryx(p)，实现步骤为：
(3a)利用双层平行圆阵列天线虚拟平移次数p、每次平移距离d和第i个波达方向信号源的俯仰角计算对角矩阵φp：
其中，diag(a1,a2,...am)表示对角线元素为a1,a2,...am的对角矩阵；
(3b)利用上层子圆阵的导向矢量矩阵信号矢量S(t)、上层子圆阵的噪声矢量N1(t)、下层子圆阵的噪声矢量N2(t)和对角矩阵φp，计算第p次轴向虚拟平移的上层子圆阵输出矢量Xp(t)和下层子圆阵输出矢量Yp(t)：
(3c)利用信号矢量S(t)，计算信号矢量S(t)的自协方差矩阵Rss：
Rss＝E[S(t)×(S(t))H]，
其中，E[·]表示对·求数学期望，(·)H表示·的共轭转置；
(3d)利用上层子圆阵的噪声矢量N1(t)和下层子圆阵的噪声矢量N2(t)，计算上层子圆阵和下层子圆阵噪声矢量相关矩阵Rnn：
Rnn＝E[Ni1(t)×(Nj1(t))H]，
其中，i1和j1表示上层子圆阵的噪声矢量N1(t)和下层子圆阵的噪声矢量的下标N2(t)，且i1,j1＝1,2；
(3e)利用上层子圆阵的导向矢量矩阵对角矩阵φp、信号的协方差矩阵Rss和上层子圆阵和下层子圆阵噪声矢量相关矩阵Rnn，计算第p次平移后的自协方差矩阵Rxx(p)和互协方差矩阵Ryx(p)；
(4)利用轴向平移总次数Q，计算信号子空间的自相关矩阵的伪逆矩阵实现步骤为：
(4a)利用第p次平移后的自协方差矩阵Rxx(p)和互协方差矩阵Ryx(p)，计算Q次平移的协方差矩阵平均值和互协方差矩阵平均值
(4b)对协方差矩阵平均值进行特征分解，得到(N-M)个小特征值λl(l＝M+1,...N)，并计算噪声方差σ12：
(4c)利用协方差矩阵平均值和噪声方差σ12，计算信号子空间的自相关矩阵
其中，I是N×N的单位矩阵；
(4d)对信号子空间的自相关矩阵进行特征分解，得到N个特征值μk1(k1＝1,2,...N)和N个特征向量νk1(k1＝1,2,...N)，并计算信号子空间的自相关矩阵的伪逆矩阵
(5)利用信号子空间的自相关矩阵的伪逆矩阵和步骤(4a)中得到的互协方差矩阵平均值计算波达方向矩阵R；
(6)对波达方向矩阵R进行特征分解，得到M个大特征值αe(e＝1,2,...,M)和其对应的特征向量βe(e＝1,2,...,M)，并计算入射信号源的俯仰角和方位角θ1e。
本发明与现有技术相比，具有如下优点：
1、本发明利用信号子空间的自相关矩阵的伪逆矩阵和互协方差矩阵平均值得到了包含源信号空间基底的波达方向矩阵，对该矩阵进行特征分解可得到对应源信号方向元素的特征值和对应源信号方向矢量的特征向量，与现有技术中利用并不能直接对应源信号方向矢量的中间参数的特征向量求解波达方向相比，分辨率高，尤其是当多个信号的入射角过于接近时，分辨率更高，可以提高DOA估计精度。
2、本发明直接用波达方向矩阵分解得到的特征值和特征向量来求解入射信号源的方位角和俯仰角，与现有技术中需要通过谱峰搜索来估计信号的波达方向相比，更具直接性，可以减小运算量。
附图说明
图1为本发明的实现流程图。
具体实施方式
以下结合附图，对本发明作进一步详细说明。
参照图1，一种基于双层平行圆阵列天线的DOA估计方法，包括如下步骤：
步骤1，获取双层平行圆阵列天线的相关参数，包括上层子圆阵的阵元数N和半径R；获取信号源的相关参数，包括信号矢量S(t)、信号源个数M、第i个信号源的波达方向方位角和俯仰角和信号源频率f；获取噪声相关参数，包括上层子圆阵的噪声矢量N1(t)和下层子圆阵的噪声矢量N2(t)；获取双层平行圆阵列天线轴向虚拟平移操作的相关参数，包括双层平行圆阵列天线轴向虚拟平移次数p、每次平移距离d和轴向平移总次数Q，其中双层平行圆阵列天线是一种阵元均匀分布在圆周上的上下两层平行的圆阵列天线，N＞M+1，Q＞M，d＝c/(2×f)，p＝1，2，...Q，c表示光速；
步骤2，以上层子圆阵的圆心为参考点，计算上层子圆阵的导向矢量矩阵实现步骤为：
步骤2a，利用双层平行圆阵列天线的相关参数和信号源的相关参数，计算第i个信号源相对于参考点在上层子圆阵第k个阵元的时延τki：
步骤2b，利用第i个信号源相对于参考点在上层子圆阵第k个阵元的时延τki，计算第i个信号源的方向矢量
其中，j表示虚数单位，[·]T表示·的转置；
步骤2c，利用第i个信号源的方向矢量计算上层子圆阵的导向矢量矩阵
步骤3，利用双层平行圆阵列天线轴向虚拟平移次数p和每次平移距离d，计算第p次平移后的自协方差矩阵Rxx(p)和互协方差矩阵Ryx(p)，实现步骤为：
步骤3a，利用双层平行圆阵列天线虚拟平移次数p、每次平移距离d和第i个波达方向信号源的俯仰角计算对角矩阵φp：
其中，diag(a1,a2,...am)表示对角线元素为a1,a2,...am的对角矩阵；
步骤3b，利用上层子圆阵的导向矢量矩阵信号矢量S(t)、上层子圆阵的噪声矢量N1(t)、下层子圆阵的噪声矢量N2(t)和对角矩阵φp，计算第p次轴向虚拟平移的上层子圆阵输出矢量Xp(t)和下层子圆阵输出矢量Yp(t)：
步骤3c，利用信号矢量S(t)，计算信号矢量S(t)的自协方差矩阵Rss：
Rss＝E[S(t)×(S(t))H]，
其中，E[·]表示对·求数学期望，(·)H表示·的共轭转置；
步骤3d，利用上层子圆阵的噪声矢量N1(t)和下层子圆阵的噪声矢量N2(t)，计算上层子圆阵和下层子圆阵噪声矢量相关矩阵Rnn：
Rnn＝E[Ni1(t)×(Nj1(t))H]，
其中，i1和j1表示上层子圆阵的噪声矢量N1(t)和下层子圆阵的噪声矢量的下标N2(t)，且i1,j1＝1,2；
步骤3e，利用上层子圆阵的导向矢量矩阵对角矩阵φp、信号的协方差矩阵Rss和上层子圆阵和下层子圆阵噪声矢量相关矩阵Rnn，计算第p次平移后的自协方差矩阵Rxx(p)和互协方差矩阵Ryx(p)：
步骤4，利用轴向平移总次数Q，计算信号子空间的自相关矩阵的伪逆矩阵实现步骤为：
步骤4a，利用第p次平移后的自协方差矩阵Rxx(p)和互协方差矩阵Ryx(p)，计算Q次平移的协方差矩阵平均值和互协方差矩阵平均值
步骤4b，对协方差矩阵平均值进行特征分解，得到(N-M)个小特征值λl(l＝M+1,...N)，并计算噪声方差σ12：
步骤4c，利用协方差矩阵平均值和噪声方差σ12，计算信号子空间的自相关矩阵
其中，I是N×N的单位矩阵；
步骤4d，对信号子空间的自相关矩阵进行特征分解，得到N个特征值μk1(k1＝1,2,...N)和N个特征向量νk1(k1＝1,2,...N)，并计算信号子空间的自相关矩阵的伪逆矩阵
步骤5，利用信号子空间的自相关矩阵的伪逆矩阵和步骤(4a)中得到的互协方差矩阵平均值计算波达方向矩阵R：
步骤6，对波达方向矩阵R进行特征分解，得到M个大特征值αe(e＝1,2,...,M)和其对应的特征向量βe(e＝1,2,...,M)，并计算入射信号源的俯仰角和方位角θ1e：
其中，Arg(·)表示的相位,βe(·)表示βe的第·个元素。
以上描述仅是本发明的一个具体实例，显然对于本领域的专业人员来说，在了解了本发明内容和原理后，都可能在不背离本发明原理、结构的情况下，进行形式和细节上的各种修正和改变，但是这些基于本发明思想的修正和改变仍在本发明的权利要求保护范围之内。