本发明公开了一种基于压缩感知的嵌套多输入多输出雷达波达方向DOA估计方法，其实现过程是：（1）利用嵌套阵的阵元排布形式确定嵌套多输入多输出雷达的发射阵形和接收阵形；（2）在雷达的每个脉冲重复周期内匹配滤波回波信号；（3）对匹配滤波输出数据的协方差矩阵进行向量化；（4）采用矩阵重构法和奇异值分解得到噪声子空间；（5）构造稀疏恢复问题中的L1范数约束项和L2范数约束项；（6）通过求解稀疏恢复问题确定目标波达方向。本发明同时利用了嵌套阵的自由度扩展特性和压缩感知的高分辨特性，可用于雷达对个数多于虚拟阵元总数的目标的精确定位。
1.一种基于压缩感知的嵌套多输入多输出雷达DOA估计方法，包括如下步骤：(1)确定收发阵形：
(1a)按照嵌套阵元位置坐标集合的计算公式，计算嵌套阵元的位置坐标集合，得到嵌套阵形；
其中，xi表示第i个嵌套阵元的位置坐标，d表示雷达发射信号的半个波长，Q表示嵌套阵元总数，{·}表示集合操作；
(1b)按照虚拟阵元位置坐标集合的计算公式，对多组不同收发阵形的多输入多输出雷达，分别计算各组雷达的虚拟阵元位置坐标集合，得到多组虚拟阵元位置坐标集合；
{xm+xn},m＝1,…,M,n＝1,…,N其中，xm表示第m个雷达发射阵元的位置坐标，xn表示第n个雷达接收阵元的位置坐标，M表示雷达发射阵元总数，N表示雷达接收阵元总数，{·}表示集合操作；
(1c)从多组虚拟阵元位置坐标集合中，选取包含所有嵌套阵形阵元位置坐标的虚拟阵元位置坐标集合，将所选取的多个虚拟阵对应的收发阵形组成多组嵌套收发阵集合；
(1d)从多组嵌套收发阵集合中，选取收发阵元总数最少的阵形，将该阵形作为嵌套多输入多输出雷达的发射阵形和接收阵形；
在雷达的每个发射脉冲重复周期内，匹配滤波嵌套多输入多输出雷达所有接收阵元上的回波信号，得到匹配滤波输出向量；
(3a)将每个雷达发射脉冲重复周期的匹配滤波输出向量与其共轭转置相乘，得到每个发射脉冲重复周期的自相关矩阵；
(3b)取多个雷达发射脉冲重复周期的自相关矩阵的平均值，得到协方差矩阵；
(3c)利用Khatri-Rao积对向量化的协方差矩阵变形，得到Khatri-Rao向量；
其中，A*表示对嵌套多输入多输出雷达的阵列流形矩阵做复共轭操作，得到的共轭阵列流形矩阵，⊙表示求Khatri-Rao积操作，g表示目标反射系数向量， 表示雷达接收阵元噪声方差，表示向量化的单位矩阵；
(3d)去除Khatri-Rao向量中的冗余行，得到标准Khatri-Rao向量；
(4a)采用矩阵重构法，求标准Khatri-Rao向量的平滑相关矩阵；
第一步，将标准Khatri-Rao向量以中心元素作为对称轴，划分成两个子向量；
第二步，以标准Khatri-Rao向量的中心元素作为主对角线元素，以与中心元素相邻的两个子向量中的元素分别作为+1和-1对角线上的元素，依此类推，得到(J+1)×(J+1)维的平滑相关矩阵，其中，J表示每个子向量中的元素个数；
(4b)对平滑相关矩阵进行奇异值分解，得到由所有小特征值对应的特征矢量组成的噪声子空间；
(5a)对Khatri-Rao向量的估计误差进行约束，构造稀疏恢复问题中的L2范数约束项；
其中，W-1表示对误差矩阵做求逆操作得到的加权矩阵，y表示Khatri-Rao向量， 表示雷达接收阵元噪声方差， 表示向量化的单位矩阵， 表示嵌套多输入多输出雷达的完备扩展阵列流形矩阵，Θ表示覆盖所有可能目标方位的空域离散角度集，gΘ表示稀疏信号向量，||·||2表示取2范数操作，γ表示L2范数约束项的上限阈值；
(5b)利用噪声子空间对稀疏信号向量加权，构造稀疏恢复问题中的L1范数约束项；
其中，diag(·)表示求对角矩阵操作，w表示对噪声子空间加权矩阵，逐行求2范数得到的加权向量，gΘ表示稀疏信号向量，Θ表示覆盖所有可能目标方位的空域离散角度集，||·||1表示取1范数操作；
(5c)在满足L2范数约束项的条件下，最小化L1范数约束项，得到稀疏恢复问题；
求解构造出的稀疏恢复问题，确定目标波达方向DOA。
2.根据权利要求1所述的基于压缩感知的嵌套多输入多输出雷达DOA估计方法，其特征在于：步骤(2)所述匹配滤波是按照下式进行的：
其中，yn(t)表示第n个雷达接收阵元在时刻t的回波信号，n＝1，2，…，N，N表示雷达接收阵元总数， 表示对第m个雷达发射阵元在时刻t的发射波形复包络，做复共轭操作后得到的共轭发射信号，m＝1，2，…，M，M表示雷达发射阵元总数。
基于压缩感知的嵌套多输入多输出雷达DOA估计方法
技术领域
[0001] 本发明属于通信技术领域，更进一步涉及雷达技术领域中的一种基于压缩感知的嵌套多输入多输出雷达波达方向(Direction-of-Arrival DOA)估计方法。本发明可用于实现对个数多于多输入多输出雷达虚拟阵元总数的目标进行精确定位。
背景技术
[0002] 多输入多输出雷达由于采用正交发射信号形成空间分集，因此相较相控阵雷达，能够有效扩展阵列孔径，提高自由度，有利于探测多目标。现有关于多输入多输出雷达阵形设计的技术主要集中在均匀布阵，通过大间距发射阵和小间距接收阵来虚拟多个等距阵元，从而提高系统的参数可辨识度。在多输入多输出雷达的波达方向DOA估计问题中，由于压缩感知作为一种高效的稀疏信号重构方法，可有效解决传统子空间类方法(如MUSIC方法、ESPRIT方法)在少快拍数、低信噪比情形下存在的性能局限，因此基于压缩感知的波达方向DOA估计方法(如JLZA方法、SPICE方法)得到了广泛的应用。
[0003] Piya Pal和P P Vaidyanathan在论文“Nested arrays:a novel approach to array processing with enhanced degrees of freedom”(IEEE Transactions on Signal Processing,2010,58(8):4167-4181)中提出一种可有效扩展阵列系统自由度的阵形设计方法—嵌套阵。该方法设计出一组具有精确闭式解的嵌套阵元位置坐标，同时利用KR积形式的接收信号协方差矩阵解决待分辨信源数多于阵元数目的参数估计问题，在工程应用中便于实现。该方法存在的不足是，仅涉及普通阵列系统的阵形优化，没有将嵌套阵形应用于多输入多输出雷达系统中，从而无法进一步扩展多输入多输出雷达系统的虚拟孔径，进而不能分辨多于多输入多输出雷达虚拟阵元数的目标。
[0004] 北京邮电大学在其申请的专利“一种压缩感知信号的重构方法”(专利申请号201210090845.5，公开号CN102624399A)中公开了一种压缩感知信号的重构方法。该方法设计出一种运算复杂度较低的迭代算法来精确恢复稀疏信号向量。该方法的具体步骤为：首先，计算经过压缩感知处理的信号与传感矩阵每一列的内积，作为第一次迭代的内积结果，然后根据本次迭代的内积结果判断是否停止迭代，若不停止，则从本次迭代的内积结果中找到绝对值最大的元素对应的索引值，将索引值加入索引集合，计算下次迭代的内积结果，并进入下次迭代过程；若停止，则将当前索引集合中各个互不相同的索引值对应的传感矩阵的列构成矩阵，并根据构成的矩阵和所述信号重构压缩感知信号。该方法存在的不足是，没有充分利用稀疏阵(如最小冗余阵，嵌套阵等)的潜在自由度扩展特性，因此不能应用于欠定波达方向DOA估计问题。
[0005] Nan Hu，Zhongfu Ye，Xu Xu和Ming Bao在论文“DOA estimation for sparse array via sparse signal reconstruction”(IEEE Transactions on Aerospace and Electronic Systems,2013,49(2):760-773)中提出一种适用于嵌套阵的压缩感知波达方向DOA估计方法。该方法利用Khatri-Rao积形式的接收信号协方差矩阵构造出运算复杂度较低的单测量矢量稀疏问题，并通过联合L1、L2范数求解该问题。该方法(以下简称L1-L2方法)适用于欠定波达方向DOA估计问题，存在的不足是，求解结果易出现伪峰，给确定目标真实方位造成困难。
发明内容
[0006] 本发明的目的在于克服上述已有技术的不足，提出一种基于压缩感知的嵌套多输入多输出雷达波达方向DOA估计方法。
[0007] 实现本发明目的的基本思路是：利用Khatri-Rao积表示向量化的协方差矩阵，得到嵌套阵的自由度扩展特性；同时使用压缩感知方法的高分辨特性，通过合理构造稀疏恢复问题，来精确估计个数多于多输入多输出雷达虚拟阵元数的目标的波达方向DOA。
[0008] 为实现上述目的，本发明具体实现步骤包括如下：
[0009] (1)确定收发阵形：
[0010] (1a)按照嵌套阵元位置坐标集合的计算公式，计算嵌套阵元的位置坐标集合，得到嵌套阵形；
[0011] (1b)按照虚拟阵元位置坐标集合的计算公式，对多组不同收发阵形的多输入多输出雷达，分别计算各组雷达的虚拟阵元位置坐标集合，得到多组虚拟阵元位置坐标集合；
[0012] (1c)从多组虚拟阵元位置坐标集合中，选取包含所有嵌套阵形阵元位置坐标的虚拟阵元位置坐标集合，将所选取的多个虚拟阵对应的收发阵形组成多组嵌套收发阵集合；
[0013] (1d)从多组嵌套收发阵集合中，选取收发阵元总数最少的阵形，将该阵形作为嵌套多输入多输出雷达的发射阵形和接收阵形。
[0014] (2)匹配滤波回波信号：
[0015] 在雷达的每个发射脉冲重复周期内，匹配滤波嵌套多输入多输出雷达所有接收阵元上的回波信号，得到匹配滤波输出向量。
[0016] (3)协方差矩阵向量化：
[0017] (3a)将每个雷达发射脉冲重复周期的匹配滤波输出向量与其共轭转置相乘，得到每个发射脉冲重复周期的自相关矩阵；
[0018] (3b)取多个雷达发射脉冲重复周期的自相关矩阵的平均值，得到协方差矩阵；
[0019] (3c)利用Khatri-Rao积对向量化的协方差矩阵变形，得到Khatri-Rao向量；
[0020] (3d)去除Khatri-Rao向量中的冗余行，得到标准Khatri-Rao向量。
[0021] (4)分离噪声子空间：
[0022] (4a)采用矩阵重构法，求标准Khatri-Rao向量的平滑相关矩阵；
[0023] (4b)对平滑相关矩阵进行奇异值分解，得到由所有小特征值对应的特征矢量组成的噪声子空间。
[0024] (5)构造稀疏恢复问题：
[0025] (5a)对Khatri-Rao向量的估计误差进行约束，构造稀疏恢复问题中的L2范数约束项；
[0026] (5b)利用噪声子空间对稀疏信号向量加权，构造稀疏恢复问题中的L1范数约束项；
[0027] (5c)在满足L2范数约束项的条件下，最小化L1范数约束项，得到稀疏恢复问题。
[0028] (6)求解目标波达方向DOA：
[0029] 利用凸优化软件包求解构造出的稀疏恢复问题，确定目标波达方向DOA。
[0030] 本发明与现有技术相比具有以下优点：
[0031] 第一，本发明在对嵌套多输入多输出雷达收发阵形的设计过程中，使用了嵌套阵的自由度扩展特性，克服了现有技术无法分辨个数多于虚拟阵元数的目标的缺点，使得本发明能够分辨比传统多输入多输出雷达更多的目标。
[0032] 第二，本发明在对目标的波达方向DOA的估计过程中采用了压缩感知方法，能够充分利用稀疏信号重构方法的高分辨特性，克服了现有技术在少快拍数、低信噪比情形下存在的性能局限，使得本发明能够提高目标波达方向DOA的估计精度。
附图说明
[0033] 图1为本发明的流程图；
[0034] 图2为本发明对15个目标波达方向DOA估计的空间谱图；
[0035] 图3为本发明与现有技术对3个目标波达方向DOA估计的均方根误差随信噪比变化曲线图。
具体实施方式
[0036] 下面结合附图对本发明做进一步的描述。
[0037] 参照附图1，本发明的具体步骤如下。
[0038] 步骤1，确定收发阵形。
[0039] 按照如下嵌套阵元位置坐标集合的计算公式，计算嵌套阵元的位置坐标集合，得到嵌套阵形：
[0040]
[0041] 其中，xi表示第i个嵌套阵元的位置坐标，d表示雷达发射信号的半个波长，Q表示嵌套阵元总数，{·}表示集合操作。
[0042] 按照如下虚拟阵元位置坐标集合的计算公式，对多组不同收发阵形的多输入多输出雷达，分别计算各组雷达的虚拟阵元位置坐标，得到多组虚拟阵元位置坐标集合：
[0043] {xm+xn},m＝1,…,M,n＝1,…,N
[0044] 其中，xm表示第m个雷达发射阵元的位置坐标，xn表示第n个雷达接收阵元的位置坐标，M表示雷达发射阵元总数，N表示雷达接收阵元总数，{·}表示集合操作。
[0045] 从多组虚拟阵元位置坐标集合中，选取包含所有嵌套阵形阵元位置坐标的虚拟阵元位置坐标集合，将所选取的多个虚拟阵对应的收发阵形组成多组嵌套收发阵集合。
[0046] 从多组嵌套收发阵集合中，选取收发阵元总数最少的阵形，将该阵形作为嵌套多输入多输出雷达的发射阵形和接收阵形。
[0047] 步骤2，匹配滤波回波信号。
[0048] 在雷达的每个发射脉冲重复周期内，按照下式匹配滤波嵌套多输入多输出雷达所有接收阵元上的回波信号，得到匹配滤波输出向量：
[0049]
[0050] 其中，yn(t)表示第n个雷达接收阵元在时刻t的回波信号，n＝1,2,…,N，N表示雷达接收阵元总数， 表示对第m个雷达发射阵元在时刻t的发射波形复包络，做复共轭操作后得到的共轭发射信号，m＝1,2,…,M，M表示雷达发射阵元总数。
[0051] 步骤3，协方差矩阵向量化。
[0052] 将每个雷达发射脉冲重复周期的匹配滤波输出向量与其共轭转置相乘，得到每个发射脉冲重复周期的自相关矩阵。
[0053] 取多个雷达发射脉冲重复周期的自相关矩阵的平均值，得到协方差矩阵。
[0054] 利用Khatri-Rao积对向量化的协方差矩阵变形，按照下式得到Khatri-Rao向量：
[0055]
[0056] 其中，A*表示对嵌套多输入多输出雷达的阵列流形矩阵做复共轭操作，得到的共轭阵列流形矩阵，⊙表示求Khatri-Rao积操作，g表示目标反射系数向量， 表示雷达接收阵元噪声方差， 表示向量化的单位矩阵。
[0057] 去除Khatri-Rao向量中的冗余行，得到标准Khatri-Rao向量。
[0058] 步骤4，分离噪声子空间。
[0059] 采用矩阵重构法，求标准Khatri-Rao向量的平滑相关矩阵，其具体步骤为：
[0060] 第一步，将标准Khatri-Rao向量以中心元素作为对称轴，划分成两个子向量；
[0061] 第二步，以标准Khatri-Rao向量的中心元素作为主对角线元素，以与中心元素相邻的两个子向量中的元素分别作为+1和-1对角线上的元素，依此类推，得到(J+1)×(J+1)维的平滑相关矩阵，其中，J表示每个子向量中的元素个数。
[0062] 对平滑相关矩阵进行奇异值分解，得到由所有小特征值对应的特征矢量组成的噪声子空间。
[0063] 步骤5，构造稀疏恢复问题。
[0064] 按照下式计算误差矩阵：
[0065]
[0066] 其中，P表示雷达发射脉冲重复周期总个数，L*表示对由协方差矩阵进行Cholesky分解得到的下三角阵做复共轭操作，得到的共轭下三角阵，表示求Kronecker积操作。
[0067] 按照下式计算完备扩展阵列流形矩阵：
[0068]
[0069] 其中， 表示对嵌套多输入多输出雷达的扩展阵列流形矩阵做复共轭操作，得到的共轭扩展阵列流形矩阵，Θ表示覆盖所有可能目标方位的空域离散角度集，⊙表示求Khatri-Rao积操作。
[0070] 对Khatri-Rao向量的估计误差进行约束，按照下式构造稀疏恢复问题中的L2范数约束项：
[0071]
[0072] 其中，W-1表示对误差矩阵做求逆操作得到的加权矩阵，y表示Khatri-Rao向量，表示雷达接收阵元噪声方差，表示向量化的单位矩阵， 表示嵌套多输入多输出雷达的完备扩展阵列流形矩阵，Θ表示覆盖所有可能目标方位的空域离散角度集，gΘ表示稀疏信号向量，||·||2表示取2范数操作，γ表示L2范数约束项的上限阈值。
[0073] 利用噪声子空间对稀疏信号向量加权，按照下式构造稀疏恢复问题中的L1范数约束项：
[0074] ||diag(w)gΘ||1
[0075] 其中，diag(·)表示求对角矩阵操作，w表示对噪声子空间加权矩阵，即由完备扩展阵列流形矩阵与噪声子空间作内积运算得到的矩阵，逐行求2范数得到的加权向量，gΘ表示稀疏信号向量，Θ表示覆盖所有可能目标方位的空域离散角度集，||·||1表示取1范数操作。
[0076] 在满足L2范数约束项的条件下，最小化L1范数约束项，得到稀疏恢复问题。
[0077] 步骤6，求解目标波达方向DOA。
[0078] 利用凸优化软件包求解构造出的稀疏恢复问题，确定目标波达方向DOA。
[0079] 下面通过仿真对本发明的效果做进一步说明。
[0080] 1.仿真条件：
[0081] 本发明的仿真是在MATLAB R2010a的软件环境下进行的。#
[0082] 2.仿真内容：
[0083] 嵌套多输入多输出雷达的3个发射阵元的位置坐标分别为：0，1.5d，6.5d，4个接收阵元的位置坐标分别为：0，0.5d，d，3d。d为雷达发射信号的半个波长。雷达脉冲重复周期总个数为300，Monte Carlo实验次数为300。假设3个独立目标的波达方向DOA分别为：-10.3°，3.2°，5.4°。波达方向DOA估计性能用均方根误差来表示，定义均方根误差
为开平方操作，E为求期望值操作，θ和 分别为波达方向DOA真实值和估计值。
[0084] 3.仿真效果分析：
[0085] 图2给出了15个独立目标波达方向DOA估计的空间谱曲线。图2中的横坐标表示波达方向，纵坐标表示空间谱幅度。图2中实线表示本发明所述波达方向DOA估计方法获得的空间谱曲线。图2中标示点划线的曲线表示所有目标的真实波达方向DOA。从图2中空间谱曲线的15个谱峰可以看出，嵌套MIMO雷达可以精确估计出个数多于虚拟阵元数的目标的波达方向DOA。在相同的发射阵元数和接收阵元数情况下，传统多输入多输出雷达最多可分辨11个独立目标，而本发明设计的嵌套多输入多输出雷达可提高多输入多输出雷达的自由度，改善目标的参数可辨识性。
[0086] 图3给出了3个独立目标波达方向DOA均方根误差随信噪比的变化曲线。图3中的横坐标表示信噪比，纵坐标表示均方根误差。图3中带有圆圈的曲线表示本发明对目标波达方向DOA的估计均方根误差获得的曲线。图3中带有方形的曲线表示MUSIC方法对目标波达方向DOA的估计均方根误差获得的曲线。图3中带有加号的曲线表示JLZA方法对目标波达方向DOA的估计均方根误差获得的曲线。图3中带有三角形的曲线表示SPICE方法对目标波达方向DOA的估计均方根误差获得的曲线。图3中带有星形的曲线表示克拉美罗下界随信噪比变化的曲线。比较图3中的五条曲线可以看出，在相同信噪比情况下，本发明所述方法对目标波达方向DOA估计的均方根误差要小于现有方法，故本发明的方法对目标波达方向DOA的估计性能要优于现有方法。
