技术领域
[0001] 本
发明属于阵列
信号处理技术领域，尤其涉及一种非均匀阵列设计和波达方向估计方法，通过灵活布阵，对多个
空域目标源进行波达方向估计。
背景技术
[0002] 在现代阵列
信号处理领域中，波达方向估计(以下简称DOA估计)占有一个重要地位，传统DOA估计主要通过MUSIC方法、ESPRIT方法及其衍生方法进行。对于一个总阵元数为N的均匀线性阵列，采用传统DOA估计方法所能识别的最大空域目标个数为N-1个。在实际中，空域目标个数大于阵元数情况时有发生，采用少于目标个数的均匀线性阵列的传统DOA估计方法将失效。因此，空域目标个数大于阵元数的DOA估计被许多学者广泛研究，最近提出的嵌套阵列和互质阵列通过合理布阵构造一等效的虚拟阵列，提高阵列
自由度(Degree of freedom，DOF)，通过利用虚拟阵列而不是原始阵列来进行DOA 估计。嵌套阵列是由两个或多个具有不同阵元间隔的线性子阵组成，由它构造的虚拟阵列是一个完全填充的均匀线性阵列，但其第一个子阵阵元间距通常较小(阵元排布
密度较大)，会引起阵元间的互耦问题。互质阵列是由两个阵元数分别为M、N的均匀线性子阵组成，其中M、N为互质整数，虽然互质阵列能获得多余阵元数的自由度DOF，但其构造的虚拟阵列不是一个完全填充的均匀线性阵列，导致在DOA估计时，只能利用互质阵列构造的虚拟阵列的一部分，这就带来了自由度和阵列孔径的降低。
[0003] 由于嵌套阵列、互质阵列每个子阵均由均匀线性子阵组成，这就导致在机载雷达或物理空间受限的情况下，时常出现无法为其子阵阵元安装找到合理
位置的问题，同时关于非均匀阵列提出的新
算法，
发明名称为基于嵌套式子阵阵列的波达方向估计方法，公开(公告)号：CN105824002A，提出的对数据协方差矩阵特征分解后利用MUSIC算法谱峰搜索进行DOA估计；以及发明名称为嵌套阵列基于K-R子空间的波达方向估计方法，公开(公告)号：CN107544051A，提出的
稀疏信号重构法利用了信号谱的稀疏性来进行DOA估计；这两种方法在进行DOA估计时，通常需要巨大的计算量。本发明的目的在于克服上述已有非均匀阵列以及DOA估计方法不足，设计一种新的非均匀阵列和波达方向估计方法，该非均匀阵列无需由均匀子阵列组成，且其生成的虚拟阵列是一个完全填充的均匀线性阵列，同时为降低计算量，通过构造一线性算子得到信号子空间进行DOA估计。本发明在保证可识别的空域目标个数大于阵元数前提下，提高了布阵的灵活性，降低DOA估计运算量。
发明内容
[0004] 为实现上述目的，本发明的技术思路是：根据嵌套阵列经过最优布阵获取阵列最大自由度(最多虚拟阵元个数)思想，来确定非均匀阵列第一个阵元和最后一个阵元位置系数，通过计算非均匀阵列其余阵元位置系数，构造非均匀阵列；对非均匀阵列构造的虚拟差分合成阵列去冗余、排序得到虚拟阵列；构造线性算子，根据虚拟阵列进行DOA估计。其实现方案包括如下：
[0005] 1)构造非均匀阵列：
[0006] 1a)根据嵌套阵列最大自由度DOF确定非均匀阵列位置系数x1和xN；
[0007] 假定非均匀阵列阵元位置为：
[0008] D＝[d1，d2，...，dn，...，dN]＝d×[x1，x2，...，xn，...，xN]
[0009] 根据嵌套阵列构造的虚拟阵列：
[0010] DULA＝[-(N2(N1+1)-1)d，...，0，...，(N2(N1+1)-1)d]
[0011] 获取最大自由度DOF＝2N2(N1+1)-1，得到非均匀阵列位置系数则本发明阵列阵元位置为
[0012] 其中d为阵元间隔，取值为入射信号最小半
波长；N1与N2为嵌套阵列子阵阵元个数，N1＝N2＝N/2；xn表示非均匀阵列第n个阵元的位置系数，n＝1，2，...，N，N为总阵元数；元素为{1～(xN-1)}的随机递增互异整数；
[0013] 1b)计算满足非均匀阵列差合矩阵 中包含所有虚拟阵元位置的位置系数[0014] 为保证本发明阵列在差合处理之后包含所有的虚拟阵元位置，且与嵌套阵列自由度一致，构造一向量 P＝[D，D，...，D]，向量P由N个D 按行排列组成，再构造另一矩阵 即矩阵 为非均匀阵列产生的差合矩阵：
[0015]
[0016] 矩阵 为反对称矩阵，只研究其上三
角形元素即可，将其按行排布得到一个行向量公式如下：
[0017]
[0018] 为使非均匀阵列位置 所产生的虚拟差合阵列是完备的，即包含所有虚拟阵元位置，构造的差合矩阵 所对应的 中应包含 个不同元素；
[0019] 1c)根据位置系数 获取非均匀阵列阵元位置；
[0020] 在1a)中，随机从{1～(xN-1)}中选取N-2个互异递增数构成 得到一组均匀阵列阵元位置D，进行
迭代，当D中阵元位置满足1b)中的差合矩阵 所对应 中包含 个不同值即可停止迭代，否则继续从1a)选取一组阵元位置，进行迭代，直至满足迭代停止条件，退出迭代，最终得到阵元位置
[0021] 2)根据所设计非均匀阵列得到接收数据X，进而得到数据协方差矩阵 Rx，向量化该协方差矩阵，得到虚拟差分合成阵列接收数据r；
[0022] 3)根据虚拟差分合成阵列接收数据r，对r进行去冗余、排序操作，最终得到虚拟阵列接收数据
[0023] 4)根据虚拟阵列接收数据 构造选择矩阵Jz对 进行秩恢复操作得到满秩矩阵V；
[0024] 5)根据虚拟阵列的满秩矩阵V，构造一线性算子Q，通过估计线性算子Q 得到信号子空间ES；
[0025] 6)基于旋转因子不变法的思想定义两个选择矩阵Jg1和Jg2，获得ES1和ES2，通过ES1和ES2的
相位关系获得旋转矩阵Φ，最后通过旋转矩阵估计波达方向
[0026] 本发明与现有阵列结构相比较有如下优点：
[0027] 1.本发明阵列虽与嵌套阵列自由度相同，但比嵌套布阵、互质布阵方式更灵活，不需要由两个或多个均匀线性子阵构成，且可产生完全填充型均匀线性虚拟阵列进行DOA估计，更适于工程实现；
[0028] 2.本发明可满足在机载或者物理空间受限的环境下，从本发明中选取一种性能相同但却适合现场位置的非均匀阵列，进行灵活稀疏布阵，得到多余阵元数的自由度，扩大了阵列孔径，使得可估计的空域目标数远超过实际阵元数，具有更好的测角性能；
[0029] 3.本发明无需对数据协方差矩阵进行特征分解，也不需通过对整个空域角度进行谱峰搜索，通过构造线性算子Q得到信号子空间ES，利用旋转因子不变法思想便可得到信号的DOA估计值，解决了在
采样点数增多、空域网格步长划分较细，使得特征分解以及普峰搜索运算量大、耗时长、实时性差的问题；
[0030] 4.本发明阵列为非均匀阵列，与传统均匀阵列相比，在阵列面积不变的情况下，能够有效减少阵元数目，降低成本，且阵元完全相同，更易于生产和降低成本。
附图说明
[0031] 图1为本发明的三种六阵元非均匀阵列及虚拟阵列排布方式图；
[0032] 图2为六阵元二阶嵌套阵列及产生虚拟阵元阵列排布方式图；
[0034] 图4为本发明方法空域目标DOA估计散点图；
[0035] 图5为本发明和二阶嵌套阵、均匀线阵对DOA估计均方根误差与
信噪比关系比较图；
具体实施方式
[0036] 为了让本发明的上述和其它目的、特征及优点能更明显，下文特举本发明
实施例，并配合所附图示，做详细说明如下：
[0037] 参照图1为本发明三种六阵元非均匀阵列及产生虚拟阵列排布方式图，从图1可看到存在多种满足迭代条件的阵元位置系数 即存在多种阵元排布方式，要说明的是，嵌套阵列排布只是本发明阵列的一个子集；
[0038] 参照图2为六阵元二阶嵌套阵列及产生虚拟阵元阵列排布方式图，本发明的非均匀阵列与嵌套阵列构造的虚拟阵列是一个完全填充的均匀线性阵列，即全部位置上虚拟阵元都是连续的，二者虚拟阵列的最大自由度DOF相同。因此，本发明阵列具有嵌套阵列的优点，不仅可识别多于阵元数的目标个数而且提高了目标方向识别的空间分辨
力，同时具有布阵灵活的优点。
[0039] 参照图3，由于嵌套阵列构造的虚拟阵列是一个完全填充的均匀线性阵列，可提高DOA估计的自由度DOF，本发明根据嵌套阵列经过最优布阵获取阵列最大自由度即最多虚拟阵元个数思想来设计一种新的非均匀阵列，构造一线性算子，利用连续的虚拟阵列进行DOA估计。为便于叙述，假定本发明阵列与二阶嵌套阵列总阵元数均为N，N为偶数，N为奇数的情况可类比得到，本发明方法步骤如下：
[0040] 步骤一、根据嵌套阵列构造的虚拟阵列获取最大自由度DOF，计算非均匀阵列阵元位置系数，得到非均匀阵列；
[0041] 1a)根据嵌套阵列最大自由度DOF确定非均匀阵列位置系数x1和xN；
[0042] 二阶嵌套阵为获取最大的阵列自由度，则每阶阵元数为N1＝N2＝N/2，其最大自由度为DOF＝2N2(N1+1)-1，由其构造的虚拟阵列可等效为一个阵元位置为：
[0043] DULA＝[-(N2(N1+1)-1)d，...，0，...，(N2(N1+1)-1)d] 的均匀线阵，[0044] 在总阵元数相同情况下为得到与二阶嵌套阵相同的阵列自由度，非均匀阵列位置矢量可表示为：
[0045] D＝[d1，d2，...，dn，...，dN]＝d×[x1，x2，...，xn，...，xN]
[0046] 第一个阵元为参考阵元其位置系数x1＝0，最后一个阵元位置系数 则[0047] 其中d为阵元间隔，取值为入射信号的最小半波长；xn表示非均匀阵列第 n个阵元的位置系数，n＝1，2，...，N，N为总阵元数； 元素为{1～(xN-1)}的随机递增互异整数；
[0048] 1b)计算满足非均匀阵列差合矩阵 中包含所有虚拟阵元位置的位置系数[0049] 为保证本发明阵列与嵌套阵列构造的虚拟阵列一致，则本发明阵列必须满足在差合之后包含所有的虚拟阵元位置，即是完备的，其具体实现如下：
[0050] 首先构造一向量 P＝[D，D，...，D]，构造另一个向量即矩阵 为非均匀阵
列产生的差合矩阵：
[0051]
[0052] 由于矩阵 为反对称矩阵，故只研究其上三角形元素即可，将其按行排布得到一个行向量 为：
[0053]
[0054] 为使非均匀阵列位置 所产生的虚拟差合阵列是完备的，即包含所有虚拟阵元位置，构造的差合矩阵 所对应的 中应包含 个不同元素；
[0055] 1c)根据位置系数 获取非均匀阵列阵元位置
[0056] 在1a)中，随机从{1～(xN-1)}中选取N-2个互异递增数构成 得到一组均匀阵列阵元位置D，进行迭代，当D中阵元位置满足1b)中的差合矩阵 所对应 中包含 个不同值即可停止迭代，否则继续从1a)选取一组阵元位置，进行迭代，直至满足迭代停止条件，退出迭代，最终得到阵元位置
[0057] 步骤二、根据所设计的非均匀阵列，得到阵列接收数据X，进而得到数据协方差矩阵Rx，向量化该协方差矩阵，得到虚拟差分合成阵列接收数据r；
[0058] 2a)本发明阵列接收数据信号为远场窄带信号，通过L次同步采样获取数据X；
[0059] 2b)根据非均匀阵列接收数据X，估计得到数据协方差矩阵Rx：
[0060]
[0061] 其中L表示快拍数，X∈CN×L，Rx∈CN×N，(·)H表示共轭转置；
[0062] 2c)根据阵列协方差矩阵Rx，计算虚拟差分合成阵列接收数据r：
[0063] r＝vec(Rx)
[0064] 其中， vec(·)表示对矩阵进行向量化操作；
[0065] 步骤三、根据虚拟差分合成阵列接收数据r，对r进行去冗余、排序操作，最终得到虚拟阵列接收数据
[0066] 由于接收数据r中包含噪声，故无法直接对r进行去冗余、排序操作，只有通过合理设计才能得到虚拟阵列接收数据 其原理是由于接收数据r元素位置与阵列产生的差合阵元位置一一对应，通过对差合阵元位置进行处理得到索引集，由索引集对r内部元素进行选取便可达到对r去冗余、排序目的，进而得到 具体实现如下：
[0067] 3a)根据1)得到的非均匀阵列，构造一列向量 列向量 由1)中差合矩阵 的元素按行
叠加组成一个列向量：
[0068]
[0069] 再构造另一向量U，其内部 个元素均匀连续变化，与非均匀阵列产生的虚拟阵元排布相同：
[0070] 对列向量 中每一个元素 与向量U的Uj进行比较，在 时，记录 在 中的索引值i，并令Uj＝N2，要说明的是，在进行比较时，只要满足 便更新索引集Γ＝Γ∪{i}，最终得到索引集Γ，其内部包含 中 个虚拟阵元的索引值，所对应的元素构成一个新的向量 其中
[0071] i＝1，2，...，N2，
[0072]
[0073] 式中N为总阵元数；
[0074] 3b)根据对 排序得到新的索引集 获取虚拟阵列接收数据
[0075] 由于索引集Γ所对应 中元素并不是按从小到大顺序排列的，因此通过对 排序即 向量 与 所对应F中的 个元素一一对应；
[0076] 至此，得到排序后的索引值 通过接收数据r元素位置与阵列产生的差合阵元位置一一对应关系，用索引值 选取接收数据r便可达到去冗余、排序目的，最终得到虚拟阵列接收数据
[0077] 其中，sort(·)表示排序运算，F表示内部元素按顺序排列后的向量，表示F中每个元素对应索引值的集合；
[0078] 步骤四、根据虚拟阵列接收数据 构造选择矩阵Jz对 进行秩恢复操作得到满秩矩阵V，其具体过程如下：
[0079] 4a)根据虚拟阵列接收数据 构造选择矩阵Jz将 分割为 个子矩阵：
[0080]
[0081] 每个子矩阵用 表示，由接收数据 的 到 行组成，其中Jz表示为：
[0082] 4b)根据 叠加得到满秩矩阵V，V可表示为：
[0083]
[0084] 其中，
[0085] 步骤五、非均匀阵列接收数据快拍数L、维度增加时，为了降低由特征分解获得信号子空间的运算量，根据满秩矩阵V，构造一线性算子，通过估计线性算子Q得到信号子空间ES，其具体过程如下：
[0086] 5a)列出线性算子Q与信号子空间ES关系式：
[0087]
[0088] 其中，IK表示K×K的单位矩阵，K为目标个数，(·)H表示共轭转置，因此只需估计出线性算子即可获得信号子空间ES；
[0089] 5b)根据满秩矩阵V，估计线性算子Q，通过线性算子Q与信号子空间 ES关系式获得信号子空间ES；
[0090] 为得到线性算子如下：
[0091]
[0092] 分割满秩矩阵
[0093] 其中 由V的前K行组成， 由V的 行组成，则线性算子能被估计如下：
[0094]
[0095] 获得线性算子Q后，通过 便可得到信号子空间；
[0096] 步骤六、基于旋转因子不变法的思想定义两个选择矩阵Jg1和Jg2，获得 ES1和ES2，通过ES1和ES2的相位关系获得旋转矩阵Φ，最后利用旋转矩阵估计波达方向 其具体实现如下：
[0097] 6a)定义两个选择矩阵：
[0098]
[0099]
[0100] 其中 表示的单位矩阵， 表示 的全1向量， 和 分别表示
和 零矩阵；
[0101] 6b)通过选择矩阵Jg1和Jg2获取信号子空间ES1和 ES2，ES1＝Jg1ES，ES2＝Jg2ES，基于旋转因子不变法的思想，利用信号子空间ES1和ES2关系获得 对Ψ∈CK×K进行特征分解 [v，q]＝eig(Ψ)得到旋转矩阵Φ∈CK×K对角线元素：
[0102]
[0103] 其中eig(·)表示特征分解， 表示对矩阵求伪逆，diag(·)表示对角阵，v和q＝[q1，q2，...，qk，...，qK]分别表示矩阵Ψ的
特征向量矩阵和特征值矩阵，θk表示第k个信号的波达方向，λ为入射信号波长；
[0104] 6c)通过旋转矩阵Φ估计波达方向
[0105]
[0106] 其中 表示k个波达方向估计值，Φk，k表示矩阵Φ的第k行第k列的元素，d为入射信号半波长，arg(·)表示取相位运算。
[0107] 本发明的效果可以通过以下的仿真结果进一步说明：
[0108] 仿真1：利用本发明非均匀阵列能估计数量大于接收阵列阵元数的目标数，图4仿真为6阵元非均匀阵列，阵元位置矢量D＝d×[0，1，3，4，9，11]，阵元个数为N＝6，d/λ＝1/2，对K＝11个远场窄带目标进行方向识别，这11个目标的波达方向为：
[0109] (-50°，-40°，-30°，-20°，-10°，0°，10°，20°，30°40°，50°)
[0110] 在信噪比SNR＝-10dB，快拍数为L＝1024，进行50次Monte-Carlo实验；
[0111] 从图4中可以看到，在低信噪比下，本发明非均匀阵列在阵元数为6的情况下能正确无误的识别出11个空域目标源，而传统六阵元均匀线阵在此情况下将不能进行目标方向识别；
[0112] 仿真2：本发明非均匀阵列与二阶嵌套阵、均匀线阵在不同信噪比条件下的性能比较如图5所示，本发明非均匀阵列阵元个数为N＝6，阵元位置矢量D＝d×[0，1，3，4，9，11]；二阶嵌套阵总阵元数为N1+N2＝6，N1和N2分别为每阶子阵阵元个数，N1＝N2＝3，其二阶嵌套阵阵元位置矢量为 Z＝d×[0，1，2，3，7，11]；均匀线阵总阵元数为6，，对空间K＝2个远场窄带目标进行方向识别，两个目标的波达方向为(θ1，θ2)＝(20°，40°)，将目标 DOA估计的结果和实际值比较，最后得到误差统计。实验信噪比SNR从0dB 到40dB，快拍数为L＝1024，进行
500次Monte-Carlo实验；
[0113] 从图5中可以看到，当信噪比在0dB下，本发明的非均匀阵列和二阶嵌套阵DOA估计
精度在0.05°附近，均匀线阵的估计精度在0.34°附近，可以得到本发明方法优于均匀线阵；当信噪比在[0，20]dB时，本发明方法的估计精度优于二阶嵌套阵方法，随着信噪比的升高，均匀线阵估计性能上升，曲线下降加快，在40dB时三者性能趋近相同；
[0114] 以上所述，仅是本发明的较佳实施例而已，并非对本发明做任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容做出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单
修改、等同变化与修饰，均仍属于本发明技术方案的范围内。
