CN110174657B - 基于秩一降维模型和块矩阵恢复的波达方向估计方法 - Google Patents

基于秩一降维模型和块矩阵恢复的波达方向估计方法 Download PDF

Info

Publication number
CN110174657B
CN110174657B CN201910516081.3A CN201910516081A CN110174657B CN 110174657 B CN110174657 B CN 110174657B CN 201910516081 A CN201910516081 A CN 201910516081A CN 110174657 B CN110174657 B CN 110174657B
Authority
CN
China
Prior art keywords
matrix
vector
representing
signal
shift
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.)
Active
Application number
CN201910516081.3A
Other languages
English (en)
Other versions
CN110174657A (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 CN201910516081.3A priority Critical patent/CN110174657B/zh
Publication of CN110174657A publication Critical patent/CN110174657A/zh
Application granted granted Critical
Publication of CN110174657B publication Critical patent/CN110174657B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO 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
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明公开了一种基于秩一降维模型和块矩阵恢复的波达方向估计方法，该方法的具体步骤为：建立雷达的接收信号模型，确定接收信号的测量矩阵；根据接收信号的测量矩阵，构造基于秩一降维的信号协方差矩阵模型；对基于秩一去噪的信号协方差矩阵模型进行稀疏重构，得到稀疏信号矢量；根据稀疏重构信号矢量，采用交替网格优化算法估计目标信源的波达方向。本发明能够重建无噪声信号协方差矩阵，避免了信号协方差矩阵和自由度的信息丢失，提高了对密集信源和多目标的DOA估计性能。

Description

基于秩一降维模型和块矩阵恢复的波达方向估计方法
技术领域
本发明属于雷达信号处理技术领域，尤其涉及基于秩一降维模型和块矩阵恢复的波达方向估计方法。
背景技术
众所周知，在非均匀噪声和有限快拍数的条件下，经典的波达方向(DOA)估计方法性能较差。
研究表明，用于波达方向估计的传统信号子空间方法，在具有足够快拍的高斯白噪声的条件下，对于波达方向估计的性能表现良好。利用稀疏重构算法的波达方向估计方法可以通过找到目标函数的稀疏解来实现信号的重建。顾名思义，稀疏重构算法计算量很大，对于快拍数的要求同样也很高。基于压缩感知理论的L1奇异值分解(L1-SVD)算法，该算法利用1范数约束和SVD可以在低信噪比环境下提高DOA估计性能。然而，上述稀疏重建算法计算复杂度大，需要的快拍数也相当多。基于具有反馈的零空间调整算法(NST+HT+FB)的交替网格优化(AGO)算法的提出很好的解决了快拍数和耗时问题。但是由于该方法会导致阵列孔径减小，因此造成DOA估计性能某种程度的损失。
基于子空间和稀疏重构的算法的优越性仅体现于均匀白噪声条件下，该方法在非均匀高斯白噪声条件下不能体现出其优势。可以通过相对于信号和噪声参数处理对数似然函数来迭代地确定输入信号的DOA估计。该方法包含的迭代过程和高度非线性优化问题导致耗时量大。
以迭代方式确定信号/噪声子空间的DOA估计方法的原理是通过求解对数似然函数或最小二乘问题来估计噪声协方差矩阵和信号/噪声子空间，然后可以从估计的信号/噪声子空间或噪声协方差矩阵完成DOA估计。但是该方法由于包含迭代过程导致耗时较大，而通过对阵列协方差矩阵进行矢量化并充分利用虚拟阵列结构能够改善耗时问题。但是该方法通过消除包括噪声方差的元素从而获得用于DOA估计的无噪声协方差稀疏向量，其中的移除操作可能导致信号协方差矩阵中包含的部分信息丢失。因此，现有传统的稀疏重构算法对于密集信源和多目标的DOA估计性能不是很理想。
发明内容
针对上述问题，本发明的目的是提出一种基于秩一降维模型和块矩阵恢复的波达方向估计方法。本发明基于秩一降维模型，得到降维的信号协方差矩阵，以消除非均匀噪声，再通过直接计算块矩阵，从降维的信号协方差矩阵中估计噪声协方差矩阵，从而重建无噪声信号协方差矩阵，避免了信号协方差矩阵和自由度的信息丢失；再采用稀疏重构和交替网格优化(AGO)算法，显著提高了对密集信源和多目标的DOA估计性能。
为了达到上述目的，本发明采用以下技术方案予以解决。
基于秩一降维模型和块矩阵恢复的波达方向估计方法，包括以下步骤：
步骤1，建立雷达的接收信号模型x(t)，确定接收信号的测量矩阵X。
(1.1)设置信号接收雷达为包含M个阵元的均匀直线阵列，则t时刻接收到的目标信源的回波信号x(t)，其表达式为：
x(t)＝As(t)+n(t)；
其中，A为导向矩阵，s(t)是信号波形矢量，n(t)＝[n1(t),n2(t),…,nM(t)]T为零均值加性非均匀复数高斯白噪声向量，且n(t)～CN(0,Q)，Q是与n(t)相关的噪声协方差矩阵。
(1.2)导向矩阵A的表达式为：
A＝[a(θ1),a(θ2),a(θl),…,a(θL)]；
其中，L是目标信源数，l＝1,2,...,L；a(θl)是M×1的导向矢量，M代表雷达阵元数，θl代表第l个入射信号角度，θl∈Θ，Θ表示角度搜索范围，其表达式为：
a(θl)＝[1,e-jα,…,e-j(M-1)α]T；
其中，α表示对应第l个DOA的阵列几何形状，α＝2πdsin(θl)/λ，d代表阵元间距，λ代表波长；
s(t)是信号波形矢量，其表达式为：
s(t)＝[s1(t),s2(t),…,sL(t)]T；
(1.3)根据接收信号模型x(t)，得到接收信号的测量矩阵X：
Figure BDA0002095080900000031
其中，J是快拍数，tz表示第z个快拍时刻，x(tz)表示tz快拍时刻天线阵列所接收的M×1维回波信号数据，x1(tz)表示tz快拍时刻天线阵列中第一个天线阵元所接收到的回波信号数据，*表示共轭操作。
步骤2，根据接收信号的测量矩阵X，构造基于秩一降维模型的信号协方差矩阵，得到降维的信号协方差矩阵。
(2.1)根据接收信号的测量矩阵X，构建基于秩一降维模型的信号协方差矩阵的第m列
Figure BDA0002095080900000032
其表达式为：
Figure BDA0002095080900000041
其中，Em是统计期望；
Figure BDA0002095080900000042
是一个秩一相关矢量；X((1:M)m,tz)表示除去第m个阵元的测量矩阵；
Figure BDA0002095080900000043
是第m个阵元的测量矩阵的共轭转置；tz表示第z个快拍时刻。
(2.2)由M个秩一相关矢量
Figure BDA0002095080900000044
得到基于秩一降维模型的信号协方差矩阵的组合表达式：
Figure BDA0002095080900000045
其中，A((1:M)1,:)表示导向矩阵中第1到M行中的第1行的所有列；
Figure BDA0002095080900000046
表示信号功率矢量；Pl表示第l个信号功率。
(2.3)根据基于秩一降维模型的信号协方差矩阵的组合表达式，当矩阵中的行序号与列序号不相等时，即x≠y，基于秩一降维模型的信号协方差矩阵
Figure BDA0002095080900000047
的第x行第y列元素
Figure BDA0002095080900000048
的表达式为：
Figure BDA0002095080900000049
其中，d代表阵元间距，λ代表波长；
(2.4)根据基于秩一降维模型的信号协方差矩阵的第x行第y列元素
Figure BDA00020950809000000410
的表达式，得到降维的信号协方差矩阵模型
Figure BDA00020950809000000411
为：
Figure BDA00020950809000000412
其中，
Figure BDA0002095080900000051
表示L个信号功率的总和，
Figure BDA0002095080900000052
Figure BDA0002095080900000053
j为虚数单位，d代表阵元间距，λ代表波长。
步骤3，根据降维的信号协方差矩阵模型
Figure BDA0002095080900000054
得到全信号协方差矩阵
Figure BDA0002095080900000055
通过块矩阵恢复出无噪声信号协方差矩阵R0。
(3.1)根据降维的信号协方差矩阵
Figure BDA0002095080900000056
构造全信号协方差矩阵
Figure BDA0002095080900000057
的表达式为：
Figure BDA0002095080900000058
其中，G表示选择矩阵，且
Figure BDA0002095080900000059
Figure BDA00020950809000000510
表示第M行为0的M×(M-1)维交换矩阵，删除交换矩阵
Figure BDA00020950809000000511
的第M行元素即形成(M-1)×(M-1)维的单位矩阵，上标T表示矩阵的转置。
则全信号协方差矩阵
Figure BDA00020950809000000512
的元素形式的表达式为：
Figure BDA00020950809000000513
其中，r12表示
Figure BDA00020950809000000514
的第1行第2列元素，当行序号与列序号相同时，对应矩阵
Figure BDA00020950809000000515
中的元素为0；当行序号与列序号不相同，矩阵
Figure BDA00020950809000000516
中的所有项可表示为
Figure BDA00020950809000000517
上标H表示矩阵的共轭转置。
(3.2)通过块矩阵将全信号协方差矩阵
Figure BDA00020950809000000518
表示为块矩阵形式：
Figure BDA00020950809000000519
其中，块矩阵
Figure BDA00020950809000000520
的维数分别为L×L，块矩阵
Figure BDA00020950809000000521
的维数分别为L×(M-2L)，块矩阵
Figure BDA0002095080900000061
的维数分别为(M-2L)×L，块矩阵
Figure BDA0002095080900000062
的维数为(M-2L)×(M-2L)。
(3.3)将导向矩阵A表示为块矩阵形式：
Figure BDA0002095080900000063
其中，块矩阵A11、A21的维数分别为L×L，块矩阵A31的维数分别为(M-2L)×L。
(3.4)根据行序号与列序号不相同时的全信号协方差矩阵的表达式：
Figure BDA0002095080900000064
将行序号与列序号不相同时的全信号协方差矩阵
Figure BDA0002095080900000065
表示为：
Figure BDA0002095080900000066
其中，上标H表示矩阵的共轭转置。
(3.5)设定导向矩阵A是列满秩矩阵，将全信号协方差矩阵
Figure BDA0002095080900000067
中的块矩阵
Figure BDA0002095080900000068
表示为：
Figure BDA0002095080900000069
其中，
Figure BDA00020950809000000610
表示全信号协方差矩阵
Figure BDA00020950809000000611
的第3行3列。
(3.6)依据全信号协方差矩阵
Figure BDA00020950809000000612
主对角线上的元素相等，由
Figure BDA00020950809000000613
得到全信号协方差矩阵
Figure BDA00020950809000000614
主对角线上第
Figure BDA00020950809000000615
项的信号功率
Figure BDA00020950809000000616
Figure BDA0002095080900000071
其中，
Figure BDA0002095080900000072
表示块矩阵
Figure BDA0002095080900000073
的第
Figure BDA0002095080900000074
行第
Figure BDA0002095080900000075
列元素；
(3.7)根据全信号协方差矩阵
Figure BDA0002095080900000076
和全信号协方差矩阵
Figure BDA0002095080900000077
主对角线上第
Figure BDA0002095080900000078
项的信号功率
Figure BDA0002095080900000079
估计无噪声信号协方差矩阵R0：
Figure BDA00020950809000000710
其中，
Figure BDA00020950809000000711
是一个维数为M×M的单位矩阵。
步骤4，对无噪声信号协方差矩阵R0进行稀疏重构，得到稀疏重构信号矢量r。
(4.1)根据无噪声信号协方差矩阵R0，构造一个维数为2(M-1)×1的信号矢量r'，信号矢量r'的第
Figure BDA00020950809000000712
项的表达式为：
Figure BDA00020950809000000713
其中，
Figure BDA00020950809000000714
表示R0的第
Figure BDA00020950809000000715
行第
Figure BDA00020950809000000716
列的元素。
(4.2)将信号矢量r'进行线性表示为：
r'＝B(θ)P；
其中，P表示信号功率矢量，P＝[P1,P2,…,PL]T；B(θ)表示虚拟流形矩阵，B(θ)＝[b(θ1),…,b(θl),…,b(θL)]，第l列导向矢量b(θl)为
Figure BDA0002095080900000081
(4.3)根据实际回波信号在空域中的稀疏性，采用重构算法对信号矢量r'进行稀疏重构，得到稀疏重构信号矢量r：
Figure BDA0002095080900000082
其中，
Figure BDA0002095080900000083
表示划分的一组字典角度，一般
Figure BDA0002095080900000084
Figure BDA0002095080900000085
是具有2(M-1)＜＜Nθ属性的过完备字典，且满足受限等距属性(RIP)；
Figure BDA0002095080900000086
表示字典角度为
Figure BDA0002095080900000087
时的稀疏导向矢量，即第1列稀疏导向矢量，Nθ是字典数；
Figure BDA0002095080900000088
表示稀疏信号功率矢量，其仅有L列非零系数，Nθ＞＞L；
Figure BDA0002095080900000089
表示第1个稀疏信号功率。
步骤5，根据稀疏重构信号矢量r，采用交替网格优化算法，估计目标信源的波达方向。
与现有技术相比，本发明的有益效果为：
(1)本发明采用基于秩一降维模型获得降维的信号协方差矩阵，通过对该矩阵进行块矩阵恢复方法得到无噪声信号协方差矩阵，可以避免信号协方差矩阵和自由度的信息丢失。
(2)本发明通过对无噪声信号协方差矩阵进行矢量化稀疏重构，可以显著提高DOA估计性能和抗均匀噪声的鲁棒性。
附图说明
下面结合附图和具体实施例对本发明做进一步详细说明。
图1是本发明的实现流程图；
图2是本发明实施例中分别采用交替网格优化算法(AGO)、基于压缩感知理论的L1奇异值分解算法(L1-SVD)和He算法，在信噪比为0时的目标信源波达方向的估计结果图；
图3是本发明实施例中采用本发明方法在信噪比为0时的目标信源波达方向的估计结果图；
图4是本发明实施例中对网格优化算法(AGO)、基于压缩感知理论的L1奇异值分解算法(L1-SVD)、He算法和本发明方法，在不同信噪比条件下，目标信源波达方向的估计均方根误差与信噪比的关系图。
具体实施方式
下面结合附图对本发明的实施例及效果作进一步详细描述。
参照图1，本发明的实现步骤如下：
步骤1，建立雷达的接收信号模型x(t)，确定接收信号的测量矩阵X。
(1.1)设置信号接收雷达为包含M个阵元的均匀直线阵列，则t时刻接收到的目标信源的回波信号x(t)，其表达式为：
x(t)＝As(t)+n(t)
其中，A为导向矩阵，s(t)是信号波形矢量，n(t)＝[n1(t),n2(t),…,nM(t)]T为零均值加性非均匀复数高斯白噪声向量，且n(t)～CN(0,Q)，Q是与n(t)相关的噪声协方差矩阵。
(1.2)导向矩阵A的表达式为：
A＝[a(θ1),a(θ2),a(θl),…,a(θL)]；
其中，L是目标信源数，l＝1,2,...,L；a(θl)是M×1的导向矢量，M代表雷达阵元数，θl代表第l个入射信号角度，θl∈Θ，Θ表示角度搜索范围，其表达式为：
a(θl)＝[1,e-jα,…,e-j(M-1)α]T；
其中，α＝2πdsin(θl)/λ对应第l个DOA的阵列几何形状，d代表阵元间距，λ代表波长；
信号波形矢量s(t)的表达式为：
s(t)＝[s1(t),s2(t),…,sL(t)]T；
(1.3)根据接收信号模型x(t)，得到接收信号的测量矩阵X：
Figure BDA0002095080900000101
其中，J是快拍数，tz表示第z个快拍时刻，x(tz)表示tz快拍时刻天线阵列所接收的M×1维回波信号数据，x1(tz)表示tz快拍时刻天线阵列中第一个天线阵元所接收到的回波信号数据，*表示共轭操作。
步骤2，根据接收信号的测量矩阵X，构造基于秩一降维模型的信号协方差矩阵，得到降维的信号协方差矩阵
Figure BDA0002095080900000102
(2.1)根据接收信号的测量矩阵X，构建基于秩一降维模型的信号协方差矩阵的第m列
Figure BDA0002095080900000103
其表达式为：
Figure BDA0002095080900000104
其中，Em是统计期望；
Figure BDA0002095080900000105
是一个秩一相关矢量；X((1:M)m,tz)表示除去第m个阵元的测量矩阵；
Figure BDA0002095080900000106
是第m个阵元的测量矩阵的共轭转置；tz表示第z个快拍时刻。
(2.2)由M个秩一相关矢量
Figure BDA0002095080900000107
得到基于秩一降维模型的信号协方差矩阵的组合表达式：
Figure BDA0002095080900000111
其中，A((1:M)1,:)表示导向矩阵中第1到M行中的第1行的所有列；
Figure BDA0002095080900000112
表示信号功率矢量；Pl表示第l个信号功率。
(2.3)根据基于秩一降维模型的信号协方差矩阵的组合表达式，当矩阵中的行序号与列序号不相等时，即x≠y，基于秩一降维模型的信号协方差矩阵的第x行第y列元素
Figure BDA0002095080900000113
的表达式为：
Figure BDA0002095080900000114
(2.4)根据基于秩一降维模型的信号协方差矩阵的第x行第y列元素
Figure BDA0002095080900000115
的表达式，得到降维的信号协方差矩阵
Figure BDA0002095080900000116
为：
Figure BDA0002095080900000117
其中，
Figure BDA0002095080900000118
表示L个信号功率的总和，
Figure BDA0002095080900000119
步骤3，根据降维的信号协方差矩阵
Figure BDA00020950809000001110
得到全信号协方差矩阵
Figure BDA00020950809000001111
通过块矩阵恢复出无噪声信号协方差矩阵R0：
(3.1)根据降维的信号协方差矩阵
Figure BDA00020950809000001112
构造全信号协方差矩阵
Figure BDA00020950809000001113
的表达式为：
Figure BDA0002095080900000121
其中，G表示选择矩阵，且
Figure BDA0002095080900000122
Figure BDA0002095080900000123
表示第M行为0的M×(M-1)维交换矩阵，删除交换矩阵
Figure BDA0002095080900000124
的第M行元素即形成(M-1)×(M-1)维的单位矩阵，上标T表示矩阵的转置。
则全信号协方差矩阵
Figure BDA0002095080900000125
的元素形式的表达式为：
Figure BDA0002095080900000126
其中，当行序号与列序号相同时，对应矩阵
Figure BDA0002095080900000127
中的元素为0；当行序号与列序号不相同，矩阵
Figure BDA0002095080900000128
中的所有项可表示为
Figure BDA0002095080900000129
上标H表示矩阵的共轭转置。
(3.2)通过块矩阵将全信号协方差矩阵
Figure BDA00020950809000001210
表示为块矩阵形式：
Figure BDA00020950809000001211
其中，块矩阵
Figure BDA00020950809000001212
的维数分别为L×L，块矩阵
Figure BDA00020950809000001213
的维数分别为L×(M-2L)，块矩阵
Figure BDA00020950809000001214
的维数分别为(M-2L)×L，块矩阵
Figure BDA00020950809000001215
的维数为(M-2L)×(M-2L)。
(3.3)将导向矩阵A表示为块矩阵形式：
Figure BDA00020950809000001216
其中，块矩阵A11、A21的维数分别为L×L，块矩阵A31的维数分别为(M-2L)×L。
(3.4)根据行序号与列序号不相同时的全信号协方差矩阵的表达式：
Figure BDA0002095080900000131
将行序号与列序号不相同时的全信号协方差矩阵
Figure BDA0002095080900000132
表示为：
Figure BDA0002095080900000133
其中，上标H表示矩阵的共轭转置。
(3.5)设定导向矩阵A是列满秩矩阵，将全信号协方差矩阵
Figure BDA0002095080900000134
中的块矩阵
Figure BDA0002095080900000135
表示为：
Figure BDA0002095080900000136
其中，
Figure BDA0002095080900000137
表示全信号协方差矩阵
Figure BDA0002095080900000138
的第3行3列。
(3.6)依据全信号协方差矩阵
Figure BDA0002095080900000139
主对角线上的元素相等，由
Figure BDA00020950809000001310
得到全信号协方差矩阵
Figure BDA00020950809000001311
主对角线上第
Figure BDA00020950809000001312
项的信号功率
Figure BDA00020950809000001313
Figure BDA00020950809000001314
其中，
Figure BDA00020950809000001315
表示块矩阵
Figure BDA00020950809000001316
的第
Figure BDA00020950809000001317
行第
Figure BDA00020950809000001318
列元素；
(3.7)根据全信号协方差矩阵
Figure BDA00020950809000001319
和全信号协方差矩阵
Figure BDA00020950809000001320
主对角线上第
Figure BDA00020950809000001321
项的信号功率
Figure BDA00020950809000001322
估计无噪声信号协方差矩阵R0：
Figure BDA0002095080900000141
其中，
Figure BDA0002095080900000142
是一个M×M的单位矩阵。
步骤4，对无噪声信号协方差矩阵R0进行稀疏重构，得到稀疏重构信号矢量r。
(4.1)根据无噪声信号协方差矩阵R0，构造一个维数为2(M-1)×1的信号矢量r'，信号矢量r'的第
Figure BDA0002095080900000143
项的表达式为：
Figure BDA0002095080900000144
其中，
Figure BDA0002095080900000145
表示无噪声信号协方差矩阵R0的第
Figure BDA0002095080900000146
行第
Figure BDA0002095080900000147
列的元素。
(4.2)将信号矢量r'进行线性表示为：
r'＝B(θ)P；
其中，P表示信号功率矢量，P＝[P1,P2,…,PL]T；B(θ)表示虚拟流形矩阵，B(θ)＝[b(θ1),…,b(θl),…,b(θL)]，第l列导向矢量b(θl)为
Figure BDA0002095080900000151
(4.3)根据实际回波信号在空域中的稀疏性，采用重构算法对信号矢量r'进行稀疏重构，得到稀疏重构信号矢量r:
Figure BDA0002095080900000152
其中，
Figure BDA0002095080900000153
Figure BDA0002095080900000154
表示划分的一组字典角度，一般
Figure BDA0002095080900000155
Figure BDA0002095080900000156
是具有2(M-1)＜＜Nθ属性的过完备字典，且满足受限等距属性(RIP)，Nθ是字典数；
Figure BDA0002095080900000157
表示稀疏信号功率矢量，其仅有L列非零系数，Nθ＞＞L。
步骤5，根据稀疏重构信号矢量r，采用交替网格优化算法估计目标信源的波达方向。
其具体实施步骤为：
(5.1)按照下式，计算导向矩阵A的右逆矩阵：
Φa＝AH(A·AH)-1；
其中，Φa表示导向矩阵A的右逆矩阵，A表示导向矩阵，H表示共轭转置操作，上标-1表示求逆操作；
(5.2)按照下式，计算导向矩阵A的正交投影矩阵：
Φ＝I-ΦaA；
其中，Φ表示导向矩阵A的正交投影矩阵，I表示单位矩阵，Φa表示导向矩阵A的右逆矩阵，A表示导向矩阵；
(5.3)按照下式，计算初始恢复矢量:
Figure BDA0002095080900000158
其中，
Figure BDA0002095080900000161
表示初始恢复矢量，Φa表示导向矩阵A的右逆矩阵，r表示稀疏重构信号；
(5.4)估计信源数：
(5.4a)将动态信源数初始化为1；
(5.4b)按照下式，对第k次内循环中的恢复矢量进行排序操作：
Figure BDA0002095080900000162
其中，
Figure BDA0002095080900000163
表示恢复矢量
Figure BDA0002095080900000164
取模值后按降序重排的矢量，l表示外循环次数，k表示内循环次数，T表示排序操作后记录
Figure BDA0002095080900000165
中每一个元素在恢复矢量
Figure BDA0002095080900000166
中对应元素的下标组成的索引指标集，|·|表示取模值操作，sort(|·|,'descend')表示降序排列操作；
(5.4c)按照下式，计算第k+1次内循环中的恢复矢量：
Figure BDA0002095080900000167
其中，
Figure BDA0002095080900000168
表示第k+1次内循环中的恢复矢量，
Figure BDA0002095080900000169
表示第k次内循环中的恢复矢量，uk表示第k次内循环中的中间辅助矢量，Φ表示导向矩阵A的正交投影矩阵；
(5.4d)按照下式，计算内循环相对误差值：
Figure BDA00020950809000001610
其中，H2表示内循环相对误差值，
Figure BDA00020950809000001611
和
Figure BDA00020950809000001612
分别表示内循环次数为k+1和k时的恢复矢量，||·||2表示取2范数操作；
(5.4e)判断内循环相对误差值H2是否大于10-3，若是，则执行步骤(5.4b)，否则，执行步骤(5.4f)；
(5.4f)将第l次外循环中的动态信源数加1，用加1后的动态信源数作为下一次外循环中的动态信源数；
(5.4g)按照下式，计算失配相对误差：
Figure BDA0002095080900000171
其中，γl+1表示第l+1次外循环中的失配相对误差，l表示外循环次数，A表示导向矩阵，uk表示第k次内循环中的中间辅助矢量，X表示测量矩阵，||·||2表示取2范数操作；
(5.4h)按照下式，计算外循环相对误差值：
H1＝γl+1-γl|；
其中，H1表示外循环相对误差值，γl+1和γl分别表示外循环次数为l+1和l时的失配相对误差，|·|表示取模值操作；
(5.4i)判断外循环相对误差值H1是否大于0.05，若是，则执行步骤(4.4b)，否则，执行步骤(4.4j)；
(5.4j)将外循环结束时的动态信源数值作为信源数的估计值；
(5.5)第一次估计目标到达角：
(5.5a)按照下式，寻找峰值位置矢量：
Figure BDA0002095080900000172
其中，pV表示恢复矢量
Figure BDA0002095080900000173
取模值后按降序重排的峰值矢量，pI表示恢复矢量元素取模值降序重排后，恢复矢量原下标值经过重新排列得到的峰值位置矢量，findpeaks(·,'descend')表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作；
(5.5b)将峰值位置矢量中的第一个元素值，放入信源位置矢量的第一个位置，将峰值位置矢量中的下一个元素值依次放入信源位置矢量的第二个位置，直至放入信源位置矢量的元素个数与估计信源数的值相同时停止取值，得到最终信源位置矢量；
(5.5c)取出角度搜索范围Θ中与最终信源位置矢量中元素值相同的下标值对应的角度值，将取出的角度值放入第一次估计的目标到达角矢量
Figure BDA00020950809000001814
中；
(5.6)按照下式，寻找峰值位置矢量：
Figure BDA0002095080900000181
其中，pV表示恢复矢量
Figure BDA0002095080900000182
取模值后按降序重排的峰值矢量，pI表示恢复矢量元素取模值降序重排后，恢复矢量原下标值经过重新排列得到的峰值位置矢量，findpeaks(·,'descend')表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作；
(5.7)按照下式，计算当前代价函数：
Figure BDA0002095080900000183
其中，F(M)表示当前代价函数，pV(s)表示峰值矢量pV中的第s个元素，s表示动态信源数，∑(·)表示求和，
Figure BDA0002095080900000184
表示恢复矢量
Figure BDA0002095080900000185
中的第i个元素，
Figure BDA0002095080900000186
Figure BDA0002095080900000187
表示不属于符号，
Figure BDA0002095080900000188
表示峰值位置矢量pI的前
Figure BDA0002095080900000189
个元素，
Figure BDA00020950809000001810
表示信源数的估计值，|·|表示取模值操作；
(5.8)按照下式，计算当前峰值：
Figure BDA00020950809000001811
其中，
Figure BDA00020950809000001812
表示当前峰值，
Figure BDA00020950809000001813
表示当前峰值对应的位置下标，pV(s)表示峰值矢量pV中的第s个元素，s表示动态信源数，pI(s)表示峰值位置矢量pI中的第s个元素；
(5.9)计算左移代价函数：
(5.9a)用角度搜索范围中当前位置下标对应的角度减去可调栅格步长后得到的搜索角度范围作为左移角度搜索范围；
(5.9b)按照下式，计算左移导向矩阵：
A(L)＝[a(θ1),a(θ2),…,a(θi),…]
其中，A(L)表示左移导向矩阵，M表示天线阵列中的阵元个数，θi∈Θ(L)，∈表示属于符号，Θ(L)表示左移角度搜索范围，i表示搜索角度θi在左移角度搜索范围Θ(L)中的序号；
(5.9c)按照下式，计算左移导向矩阵的正交投影矩阵：
Φ(L)＝I-(A(L))H(A(L)(A(L))H)-1A(L)
其中，Φ(L)表示左移导向矩阵A(L)的正交投影矩阵，I表示单位矩阵，A(L)表示左移导向矩阵，H表示共轭转置操作，-1表示求逆操作；
(5.9d)按照下式，计算左移恢复矢量：
Figure BDA0002095080900000191
其中，
Figure BDA0002095080900000192
表示左移恢复矢量，u(L)表示左移辅助矢量，Φ(L)表示左移导向矩阵A(L)的正交投影矩阵；
(5.9e)按照下式，寻找左移峰值位置矢量：
Figure BDA0002095080900000193
其中，
Figure BDA0002095080900000194
表示左移恢复矢量
Figure BDA0002095080900000195
取模值后按降序重排的左移峰值矢量，
Figure BDA0002095080900000196
表示左移恢复矢量元素取模值降序重排后，左移恢复矢量原下标值经过重新排列得到的左移峰值位置矢量，findpeaks(·,'descend')表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作；
(5.9f)按照下式，计算左移代价函数：
Figure BDA0002095080900000201
其中，F(L)表示左移代价函数，
Figure BDA0002095080900000202
表示左移峰值矢量
Figure BDA0002095080900000203
中的第s个元素，s表示动态信源数，∑(·)表示求和，
Figure BDA0002095080900000204
表示左移恢复矢量
Figure BDA0002095080900000205
中的第i个元素，
Figure BDA0002095080900000206
Figure BDA0002095080900000207
表示不属于符号，
Figure BDA0002095080900000208
表示左移峰值位置矢量
Figure BDA0002095080900000209
的前
Figure BDA00020950809000002010
个元素，
Figure BDA00020950809000002011
表示信源数的估计值，|·|表示取模值操作；
(5.10)计算右移代价函数：
(5.10a)用角度搜索范围中当前位置下标对应的角度加上可调栅格步长后得到的搜索角度范围作为右移角度搜索范围；
(5.10b)按照下式，计算右移导向矩阵：
A(R)＝[a(θ1),a(θ2),…,a(θi),…]
其中，A(R)为右移导向矩阵，M表示天线阵列中的阵元个数，θi∈Θ(R)，∈表示属于符号，Θ(R)表示右移角度搜索范围，i表示搜索角度θi在右移角度搜索范围Θ(R)中的序号；
(5.10c)按照下式，计算右移导向矩阵的正交投影矩阵：
Φ(R)＝I-(A(R))H(A(R)(A(R))H)-1A(R)
其中，Φ(R)表示右移导向矩阵A(R)的正交投影矩阵，I表示单位矩阵，A(R)表示右移导向矩阵，H表示共轭转置操作，-1表示求逆操作；
(5.10d)按照下式，计算右移恢复矢量：
Figure BDA0002095080900000211
其中，
Figure BDA0002095080900000212
表示右移恢复矢量，u(R)表示右移辅助矢量，Φ(R)表示右移导向矩阵A(R)的正交投影矩阵；
(5.10e)按照下式，寻找右移峰值位置矢量：
Figure BDA0002095080900000213
其中，
Figure BDA0002095080900000214
表示右移恢复矢量
Figure BDA0002095080900000215
取模值后按降序重排的右移峰值矢量，
Figure BDA0002095080900000216
表示右移峰值矢量元素取模值降序重排后，右移恢复矢量原下标值经过重新排列得到的右移峰值位置矢量，findpeaks(·,'descend')表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作；
(5.10f)按照下式，计算右移代价函数：
Figure BDA0002095080900000217
其中，F(R)表示右移代价函数，
Figure BDA0002095080900000218
表示右移峰值矢量
Figure BDA0002095080900000219
中的第s个元素，s表示动态信源数，∑(·)表示求和，
Figure BDA00020950809000002110
表示右移恢复矢量
Figure BDA00020950809000002111
中的第i个元素，
Figure BDA00020950809000002112
Figure BDA00020950809000002113
表示不属于符号，
Figure BDA00020950809000002114
表示右移峰值位置矢量
Figure BDA00020950809000002115
的前
Figure BDA00020950809000002116
个元素，
Figure BDA00020950809000002117
表示信源数的估计值，|·|表示取模值操作；
(5.11)更新当前栅格参数：
(5.11a)判断左移代价函数是否同时大于当前代价函数和右移代价函数，若是，则执行步骤(5.11b)，否则，执行步骤(5.11e)；
(5.11b)将步骤(5.7)的当前代价函数的值更新为左移代价函数的值；
(5.11c)将恢复矢量的元素值更新为左移恢复矢量的元素值；
(5.11d)按照下式，更新目标到达角矢量：
Figure BDA0002095080900000221
其中，
Figure BDA0002095080900000222
表示更新后的目标到达角矢量，Θ(L)表示左移角度搜索范围，
Figure BDA0002095080900000223
表示左移峰值位置矢量
Figure BDA0002095080900000224
的第s个元素，s表示动态信源数；
(5.11e)判断右移代价函数是否同时大于当前代价函数和左移代价函数，若是，则执行步骤(5.11f)，否则，执行步骤(5.12)；
(5.11f)将当前代价函数的值更新为右移代价函数的值；
(5.11g)将恢复矢量的元素值更新为右移恢复矢量的元素值；
(5.11h)按照下式，更新目标到达角估计矢量：
Figure BDA0002095080900000225
其中，
Figure BDA0002095080900000226
表示更新后的目标到达角估计矢量，Θ(R)表示右移角度搜索范围，
Figure BDA0002095080900000227
表示右移峰值位置矢量
Figure BDA0002095080900000228
的第s个元素，s表示动态信源数；
(5.12)判断可调栅格步长△是否大于栅格优化精度ξθ，若是，则执行步骤(5.9)，否则，执行步骤(5.13)；
(5.13)按照下式，更新目标到达角估计矢量：
Figure BDA0002095080900000229
其中，
Figure BDA00020950809000002210
表示更新后的目标到达角估计矢量的第s个元素值，s表示动态信源数，
Figure BDA00020950809000002211
表示更新后的目标到达角估计矢量，
Figure BDA00020950809000002212
表示当前峰值对应的位置下标；
(5.14)判断动态信源数是否小于信源数，若是，则执行步骤(5.6)，否则，执行步骤(5.15)；
(5.15)获得第二次目标到达角估计矢量：
将步骤(5.13)中更新后的目标到达角矢量作为第二次的目标到达角估计矢量。
仿真实验
本发明的效果可通过以下仿真实验进一步说明。
(1)仿真参数：
本发明采用M＝24个阵元，且阵元间距为半波长(d＝λ2)的均匀线性阵列。采用在-90°至90°范围内具有1°间隔的一般网格。阈值τ和误差常数ξ分别为5M和5×10-4。假定非均匀噪声协方差矩阵为Q＝diag{Qn,3*Qn,1.5*Qn}，其中Qn＝[2.0,10,2.5,5.0,0.5,1.5,3.0,5.0]。采用三个不相关的密集间隔信号源，其中DOA为{-5°,1°,3°}。
参数设置如表1：
表1***仿真参数
Figure BDA0002095080900000231
(2)仿真内容：
仿真1，在上述仿真参数下，分别采用AGO算法、L1-SVD算法和He算法在信噪比为0的条件下，对目标的回波信号进行波达方向估计，结果如图2所示。
由图2可以看出，当信号源密集放置时，AGO和He方法无法区分信号；此外，经典的L1-SVD方法性能显著恶化。
仿真2，在上述仿真参数下，采用本发明方法在信噪比为0的条件下，对回波信号进行波达方向估计，结果如图3所示。
由图3可以看出，本发明方法能够确定三个密集的不相关信号；比较图2和图3可以看出，本发明方法在旁瓣较低的情况下，由于获得了无噪声信号协方差矩阵，估计出了三个密集不相关信号的波达方向，说明本发明方法的波达方向估计性能优于AGO、L1-SVD和He方法。
在上述仿真参数下，分别采用AGO算法、L1-SVD算法、He算法和本发明方法在信噪比为分别为-6dB、-4dB、-2dB、0dB、2dB、4dB、6dB的条件下，对回波信号进行波达方向估计，并采用均方根误差RMSE对估计结果进行评价：
Figure BDA0002095080900000241
其中，K是蒙特卡洛实验次数，
Figure BDA0002095080900000242
表示第k次蒙特卡洛实验中第l个目标的估计角度，θl表示第l个目标的真实角度。
以上试验结果如图4所示，由图4可以看出，随着信噪比SNR的增大，四种方法的估计均方根误差RMSE均呈现下降趋势，且本发明方法的估计均方根误差在各信噪比时，较其他算法的估计均方根误差最小，说明本发明方法在不同信噪比的情况下，均具有较优的波达方向估计性能。
本领域普通技术人员可以理解：实现上述方法实施例的全部或部分步骤可以通过程序指令相关的硬件来完成，前述的程序可以存储于一计算机可读取存储介质中，该程序在执行时，执行包括上述方法实施例的步骤；而前述的存储介质包括：ROM、RAM、磁碟或者光盘等各种可以存储程序代码的介质。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。

Claims (8)

1.基于秩一降维模型和块矩阵恢复的波达方向估计方法，其特征在于，包括以下步骤：
步骤1，建立雷达的接收信号模型x(t)，确定接收信号的测量矩阵X；
步骤2，根据接收信号的测量矩阵X，构造基于秩一降维模型的信号协方差矩阵，得到降维的信号协方差矩阵模型
Figure FDA0002095080890000011
步骤3，根据降维的信号协方差矩阵模型
Figure FDA0002095080890000012
得到全信号协方差矩阵
Figure FDA0002095080890000013
通过块矩阵恢复出无噪声信号协方差矩阵R0；
步骤4，对无噪声信号协方差矩阵R0进行稀疏重构，得到稀疏重构信号矢量r；
步骤5，根据稀疏重构信号矢量r，采用交替网格优化算法，估计目标信源的波达方向。
2.根据权利要求1所述的基于秩一降维模型和块矩阵恢复的波达方向估计方法，其特征在于，步骤1的具体步骤为：
(1.1)设置信号接收雷达为包含M个阵元的均匀直线阵列，则t时刻接收到的目标信源的回波信号，即t时刻的接收信号x(t)的表达式为：
x(t)＝As(t)+n(t)；
其中，A为导向矩阵，s(t)是信号波形矢量，n(t)＝[n1(t)，n2(t)，…，nM(t)]T为零均值加性非均匀复数高斯白噪声向量，且n(t)～CN(0，Q)，Q是与n(t)相关的噪声协方差矩阵；
(1.2)根据接收信号模型x(t)，得到接收信号的测量矩阵X：
Figure FDA0002095080890000014
其中，J是快拍数，tz表示第z个快拍时刻，x(tz)表示tz快拍时刻天线阵列所接收的M×1维回波信号数据，x1(tz)表示tz快拍时刻天线阵列中第一个天线阵元所接收到的回波信号数据，*表示共轭操作。
3.根据权利要求2所述的基于秩一降维模型和块矩阵恢复的波达方向估计方法，其特征在于，所述导向矩阵A的表达式为：
A＝[a(θ1)，a(θ2)，a(θl)，…，a(θL)]；
其中，L是目标信源数，l＝1，2，...，L；M代表雷达阵元数，θl代表第l个入射信号角度，θl∈Θ，Θ表示角度搜索范围；a(θl)是M×1的导向矢量，其表达式为：
a(θl)＝[1，e-jα，…，e-j(M-1)α]T；
其中，α表示对应第l个DOA的阵列几何形状，α＝2πdsin(θl)/λ，d代表阵元间距，λ代表波长；
所述信号波形矢量s(t)的表达式为：
s(t)＝[s1(t)，s2(t)，…，sL(t)]T；
其中，sL(t)表示第L个信号波形。
4.根据权利要求1所述的基于秩一降维模型和块矩阵恢复的波达方向估计方法，其特征在于，步骤2按照以下步骤实施：
(2.1)根据接收信号的测量矩阵X，构建基于秩一降维模型的信号协方差矩阵的第m列
Figure FDA0002095080890000021
其表达式为：
Figure FDA0002095080890000022
其中，Em是统计期望；
Figure FDA0002095080890000023
是一个秩一相关矢量；X((1：M)m，tz)表示除去第m个阵元的测量矩阵；
Figure FDA0002095080890000024
是第m个阵元的测量矩阵的共轭转置；tz表示第z个快拍时刻；
(2.2)由M个秩一相关矢量
Figure FDA0002095080890000031
得到基于秩一降维模型的信号协方差矩阵的组合表达式：
Figure FDA0002095080890000032
其中，A((1：M)1，：)表示导向矩阵中第1到M行中的第1行的所有列；
Figure FDA0002095080890000033
表示信号功率矢量；Pl表示第l个信号功率；
(2.3)根据基于秩一降维模型的信号协方差矩阵的组合表达式，当矩阵中的行序号与列序号不相等时，即x≠y，基于秩一降维模型的信号协方差矩阵的第x行第y列元素
Figure FDA0002095080890000034
的表达式为：
Figure FDA0002095080890000035
其中，d代表阵元间距，λ代表波长，θl代表第l个入射信号角度；
(2.4)根据基于秩一降维模型的信号协方差矩阵的第x行第y列元素
Figure FDA0002095080890000036
的表达式，得到降维的信号协方差矩阵
Figure FDA0002095080890000037
为：
Figure FDA0002095080890000038
其中，
Figure FDA0002095080890000039
表示L个信号功率的总和，
Figure FDA00020950808900000310
j为虚数单位，M代表雷达阵元数。
5.根据权利要求1所述的基于秩一降维模型和块矩阵恢复的波达方向估计方法，其特征在于，所述根据降维的信号协方差矩阵模型
Figure FDA0002095080890000041
得到全信号协方差矩阵
Figure FDA0002095080890000042
其具体步骤为：
根据降维的信号协方差矩阵
Figure FDA0002095080890000043
构造全信号协方差矩阵
Figure FDA0002095080890000044
的表达式为：
Figure FDA0002095080890000045
其中，G表示选择矩阵，且
Figure FDA0002095080890000046
Figure FDA0002095080890000047
表示第M行元素都为0的M×(M-1)维交换矩阵，删除交换矩阵
Figure FDA0002095080890000048
的第M行元素即形成(M-1)×(M-1)维的单位矩阵，上标T表示矩阵的转置；
则全信号协方差矩阵
Figure FDA0002095080890000049
的元素形式的表达式为：
Figure FDA00020950808900000410
其中，r12表示全信号协方差矩阵
Figure FDA00020950808900000411
中的第1行第2列元素，当行序号与列序号相同时，对应矩阵
Figure FDA00020950808900000412
中的元素为0；当行序号与列序号不相同，
Figure FDA00020950808900000413
上标H表示矩阵的共轭转置。
6.根据权利要求5所述的基于秩一降维模型和块矩阵恢复的波达方向估计方法，其特征在于，所述通过块矩阵恢复出无噪声信号协方差矩阵R0，其具体步骤为：
(3.1)通过块矩阵将全信号协方差矩阵
Figure FDA00020950808900000414
表示为块矩阵形式：
Figure FDA00020950808900000415
其中，块矩阵
Figure FDA0002095080890000051
的维数分别为L×L，块矩阵
Figure FDA0002095080890000052
的维数分别为L×(M-2L)，块矩阵
Figure FDA0002095080890000053
的维数分别为(M-2L)×L，块矩阵
Figure FDA0002095080890000054
的维数为(M-2L)×(M-2L)；
(3.2)将导向矩阵A表示为块矩阵形式：
Figure FDA0002095080890000055
其中，块矩阵A11、A21的维数分别为L×L，块矩阵A31的维数为(M-2L)×L；
(3.3)根据行序号与列序号不相同时的全信号协方差矩阵的表达式：
Figure FDA0002095080890000056
将行序号与列序号不相同时的全信号协方差矩阵
Figure FDA0002095080890000057
表示为：
Figure FDA0002095080890000058
其中，上标H表示矩阵的共轭转置；
(3.4)设定导向矩阵A是列满秩矩阵，将全信号协方差矩阵
Figure FDA0002095080890000059
中的块矩阵
Figure FDA00020950808900000510
表示为：
Figure FDA00020950808900000511
(3.5)依据全信号协方差矩阵
Figure FDA00020950808900000512
主对角线上的元素相等，由
Figure FDA00020950808900000513
得到全信号协方差矩阵
Figure FDA00020950808900000514
主对角线上第
Figure FDA00020950808900000515
项的信号功率
Figure FDA00020950808900000516
Figure FDA0002095080890000061
其中，
Figure FDA0002095080890000062
表示块矩阵
Figure FDA0002095080890000063
的第
Figure FDA0002095080890000064
行第
Figure FDA0002095080890000065
列元素；
(3.6)根据全信号协方差矩阵
Figure FDA0002095080890000066
和全信号协方差矩阵
Figure FDA0002095080890000067
主对角线上第
Figure FDA0002095080890000068
项的信号功率
Figure FDA0002095080890000069
估计无噪声信号协方差矩阵R0：
Figure FDA00020950808900000610
其中，
Figure FDA00020950808900000611
是一个维数为M×M的单位矩阵，
Figure FDA00020950808900000612
j为虚数单位，M代表雷达阵元数，θl代表第l个入射信号角度，d代表阵元间距，λ代表波长。
7.根据权利要求1所述的基于秩一降维模型和块矩阵恢复的波达方向估计方法，其特征在于，步骤4按以下步骤实施：
(4.1)根据无噪声信号协方差矩阵R0，构造一个维数为2(M-1)×1的信号矢量r′，信号矢量r′的第
Figure FDA00020950808900000613
项的表达式为：
Figure FDA00020950808900000614
其中，
Figure FDA0002095080890000071
表示无噪声信号协方差矩阵R0的第
Figure FDA0002095080890000072
行第
Figure FDA0002095080890000073
列的元素，M代表雷达阵元数；
(4.2)将信号矢量r′进行线性表示为：
r′＝B(θ)P；
其中，P表示信号功率矢量，P＝[P1，P2，…，PL]T；B(θ)表示虚拟流形矩阵，B(θ)＝[b(θ1)，…，b(θl)，…，b(θL)]，第l列导向矢量b(θl)为
Figure FDA0002095080890000074
θl代表第l个入射信号角度，d代表阵元间距，λ代表波长；
(4.3)根据实际回波信号在空域中的稀疏性，采用重构算法对信号矢量r′进行稀疏重构，得到稀疏重构信号矢量r：
Figure FDA0002095080890000075
其中，
Figure FDA0002095080890000076
表示划分的一组字典角度；
Figure FDA0002095080890000077
是具有2(M-1)＜＜Nθ属性的过完备字典，且满足受限等距属性；
Figure FDA0002095080890000078
表示字典角度为
Figure FDA0002095080890000079
时的稀疏导向矢量，即第1列稀疏导向矢量，Nθ是字典数；
Figure FDA00020950808900000710
表示稀疏信号功率矢量，
Figure FDA00020950808900000711
表示第1个稀疏信号功率。
8.根据权利要求1所述的基于秩一降维模型和块矩阵恢复的波达方向估计方法，其特征在于，步骤5的按照以下步骤实施：
(5.1)按照下式，计算导向矩阵A的右逆矩阵：
Φa＝AH(A·AH)-1，
其中，Φa表示导向矩阵A的右逆矩阵，A表示导向矩阵，H表示共轭转置操作，上标-1表示求逆操作；
(5.2)按照下式，计算导向矩阵A的正交投影矩阵：
Φ＝I-ΦaA；
其中，Φ表示导向矩阵A的正交投影矩阵，I表示单位矩阵，Φa表示导向矩阵A的右逆矩阵，A表示导向矩阵；
(5.3)按照下式，计算初始恢复矢量：
Figure FDA0002095080890000081
其中，
Figure FDA0002095080890000082
表示初始恢复矢量，Φa表示导向矩阵A的右逆矩阵，r表示稀疏重构信号；
(5.4)估计信源数：
(5.4a)将动态信源数初始化为1；
(5.4b)按照下式，对第k次内循环中的恢复矢量进行排序操作：
Figure FDA0002095080890000083
其中，
Figure FDA0002095080890000084
表示恢复矢量
Figure FDA0002095080890000085
取模值后按降序重排的矢量，l表示外循环次数，k表示内循环次数，T表示排序操作后记录
Figure FDA0002095080890000086
中每一个元素在恢复矢量
Figure FDA0002095080890000087
中对应元素的下标组成的索引指标集，|·|表示取模值操作，sort(|·|，′descend′)表示降序排列操作；
(5.4c)按照下式，计算第k+1次内循环中的恢复矢量：
Figure FDA0002095080890000088
其中，
Figure FDA0002095080890000089
表示第k+1次内循环中的恢复矢量，
Figure FDA00020950808900000810
表示第k次内循环中的恢复矢量，uk表示第k次内循环中的中间辅助矢量，Φ表示导向矩阵A的正交投影矩阵；
(5.4d)按照下式，计算内循环相对误差值：
Figure FDA0002095080890000091
其中，H2表示内循环相对误差值，
Figure FDA0002095080890000092
和
Figure FDA0002095080890000093
分别表示内循环次数为k+1和k时的恢复矢量，||·||2表示取2范数操作；
(5.4e)判断内循环相对误差值H2是否大于10-3，若是，则执行步骤(5.4b)，否则，执行步骤(5.4f)；
(5.4f)将第l次外循环中的动态信源数加1，用加1后的动态信源数作为下一次外循环中的动态信源数；
(5.4g)按照下式，计算失配相对误差：
Figure FDA0002095080890000094
其中，γl+1表示第l+1次外循环中的失配相对误差，l表示外循环次数，A表示导向矩阵，uk表示第k次内循环中的中间辅助矢量，X表示测量矩阵，||·||2表示取2范数操作；
(5.4h)按照下式，计算外循环相对误差值：
H1＝|γl+1-γl|；
其中，H1表示外循环相对误差值，γl+1和γl分别表示外循环次数为l+1和l时的失配相对误差，|·|表示取模值操作；
(5.4i)判断外循环相对误差值H1是否大于0.05，若是，则执行步骤(4.4b)，否则，执行步骤(4.4j)；
(5.4j)将外循环结束时的动态信源数值作为信源数的估计值；
(5.5)第一次估计目标到达角：
(5.5a)按照下式，寻找峰值位置矢量：
Figure FDA0002095080890000101
其中，pV表示恢复矢量
Figure FDA0002095080890000102
取模值后按降序重排的峰值矢量，pI表示恢复矢量元素取模值降序重排后，恢复矢量原下标值经过重新排列得到的峰值位置矢量，findpeaks(·，′descend′)表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作；
(5.5b)将峰值位置矢量中的第一个元素值，放入信源位置矢量的第一个位置，将峰值位置矢量中的下一个元素值依次放入信源位置矢量的第二个位置，直至放入信源位置矢量的元素个数与估计信源数的值相同时停止取值，得到最终信源位置矢量；
(5.5c)取出角度搜索范围Θ中与最终信源位置矢量中元素值相同的下标值对应的角度值，将取出的角度值放入第一次估计的目标到达角矢量
Figure FDA0002095080890000103
中；
(5.6)按照下式，寻找峰值位置矢量：
Figure FDA0002095080890000104
其中，pV表示恢复矢量
Figure FDA0002095080890000105
取模值后按降序重排的峰值矢量，pI表示恢复矢量元素取模值降序重排后，恢复矢量原下标值经过重新排列得到的峰值位置矢量，findpeaks(·，′descend′)表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作；
(5.7)按照下式，计算当前代价函数：
Figure FDA0002095080890000106
其中，F(M)表示当前代价函数，pV(s)表示峰值矢量pV中的第s个元素，s表示动态信源数，∑(·)表示求和，
Figure FDA0002095080890000111
表示恢复矢量
Figure FDA0002095080890000112
中的第i个元素，
Figure FDA0002095080890000113
Figure FDA0002095080890000114
表示不属于符号，
Figure FDA0002095080890000115
表示峰值位置矢量pI的前
Figure FDA0002095080890000116
个元素，
Figure FDA0002095080890000117
表示信源数的估计值，|·|表示取模值操作；
(5.8)按照下式，计算当前峰值：
Figure FDA0002095080890000118
其中，
Figure FDA0002095080890000119
表示当前峰值，
Figure FDA00020950808900001110
表示当前峰值对应的位置下标，pV(s)表示峰值矢量pV中的第s个元素，s表示动态信源数，pI(s)表示峰值位置矢量pI中的第s个元素；
(5.9)计算左移代价函数：
(5.9a)用角度搜索范围中当前位置下标对应的角度减去可调栅格步长后得到的搜索角度范围作为左移角度搜索范围；
(5.9b)按照下式，计算左移导向矩阵：
A(L)＝[a(θ1)，a(θ2)，…，a(θi)，…]
其中，A(L)表示左移导向矩阵，M表示天线阵列中的阵元个数，θi∈Θ(L)，∈表示属于符号，Θ(L)表示左移角度搜索范围，i表示搜索角度θi在左移角度搜索范围Θ(L)中的序号；
(5.9c)按照下式，计算左移导向矩阵的正交投影矩阵：
Φ(L)＝I-(A(L))H(A(L)(A(L))H)-1A(L)
其中，Φ(L)表示左移导向矩阵A(L)的正交投影矩阵，I表示单位矩阵，A(L)表示左移导向矩阵，H表示共轭转置操作，-1表示求逆操作；
(5.9d)按照下式，计算左移恢复矢量：
Figure FDA0002095080890000121
其中，
Figure FDA0002095080890000122
表示左移恢复矢量，u(L)表示左移辅助矢量，Φ(L)表示左移导向矩阵A(L)的正交投影矩阵；
(5.9e)按照下式，寻找左移峰值位置矢量：
Figure FDA0002095080890000123
其中，
Figure FDA0002095080890000124
表示左移恢复矢量
Figure FDA0002095080890000125
取模值后按降序重排的左移峰值矢量，
Figure FDA0002095080890000126
表示左移恢复矢量元素取模值降序重排后，左移恢复矢量原下标值经过重新排列得到的左移峰值位置矢量，findpeaks(·，′descend′)表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作；
(5.9f)按照下式，计算左移代价函数：
Figure FDA0002095080890000127
其中，F(L)表示左移代价函数，
Figure FDA0002095080890000128
表示左移峰值矢量
Figure FDA0002095080890000129
中的第s个元素，s表示动态信源数，∑(·)表示求和，
Figure FDA00020950808900001210
表示左移恢复矢量
Figure FDA00020950808900001211
中的第i个元素，
Figure FDA00020950808900001212
Figure FDA00020950808900001213
表示不属于符号，
Figure FDA00020950808900001214
表示左移峰值位置矢量
Figure FDA00020950808900001215
的前
Figure FDA00020950808900001216
个元素，
Figure FDA00020950808900001217
表示信源数的估计值，|·|表示取模值操作；
(5.10)计算右移代价函数：
(5.10a)用角度搜索范围中当前位置下标对应的角度加上可调栅格步长后得到的搜索角度范围作为右移角度搜索范围；
(5.10b)按照下式，计算右移导向矩阵：
A(R)＝[a(θ1)，a(θ2)，…，a(θi)，…]
其中，A(R)为右移导向矩阵，M表示天线阵列中的阵元个数，θi∈Θ(R)，∈表示属于符号，Θ(R)表示右移角度搜索范围，i表示搜索角度θi在右移角度搜索范围Θ(R)中的序号；
(5.10c)按照下式，计算右移导向矩阵的正交投影矩阵：
Φ(R)＝I-(A(R))H(A(R)(A(R))H)-1A(R)
其中，Φ(R)表示右移导向矩阵A(R)的正交投影矩阵，I表示单位矩阵，A(R)表示右移导向矩阵，H表示共轭转置操作，-1表示求逆操作；
(5.10d)按照下式，计算右移恢复矢量：
Figure FDA0002095080890000131
其中，
Figure FDA0002095080890000132
表示右移恢复矢量，u(R)表示右移辅助矢量，Φ(R)表示右移导向矩阵A(R)的正交投影矩阵；
(5.10e)按照下式，寻找右移峰值位置矢量：
Figure FDA0002095080890000133
其中，
Figure FDA0002095080890000134
表示右移恢复矢量
Figure FDA0002095080890000135
取模值后按降序重排的右移峰值矢量，
Figure FDA0002095080890000136
表示右移峰值矢量元素取模值降序重排后，右移恢复矢量原下标值经过重新排列得到的右移峰值位置矢量，findpeaks(·，′descend′)表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作；
(5.10f)按照下式，计算右移代价函数：
Figure FDA0002095080890000137
其中，F(R)表示右移代价函数，
Figure FDA0002095080890000141
表示右移峰值矢量
Figure FDA0002095080890000142
中的第s个元素，s表示动态信源数，∑(·)表示求和，
Figure FDA0002095080890000143
表示右移恢复矢量
Figure FDA0002095080890000144
中的第i个元素，
Figure FDA0002095080890000145
Figure FDA0002095080890000146
表示不属于符号，
Figure FDA0002095080890000147
表示右移峰值位置矢量
Figure FDA0002095080890000148
的前
Figure FDA0002095080890000149
个元素，
Figure FDA00020950808900001410
表示信源数的估计值，|·|表示取模值操作；
(5.11)更新当前栅格参数：
(5.11a)判断左移代价函数是否同时大于当前代价函数和右移代价函数，若是，则执行步骤(5.11b)，否则，执行步骤(5.11e)；
(5.11b)将步骤(5.7)的当前代价函数的值更新为左移代价函数的值；
(5.11c)将恢复矢量的元素值更新为左移恢复矢量的元素值；
(5.11d)按照下式，更新目标到达角矢量：
Figure FDA00020950808900001411
其中，
Figure FDA00020950808900001412
表示更新后的目标到达角矢量，Θ(L)表示左移角度搜索范围，
Figure FDA00020950808900001413
表示左移峰值位置矢量
Figure FDA00020950808900001414
的第s个元素，s表示动态信源数；
(5.11e)判断右移代价函数是否同时大于当前代价函数和左移代价函数，若是，则执行步骤(5.11f)，否则，执行步骤(5.12)；
(5.11f)将当前代价函数的值更新为右移代价函数的值；
(5.11g)将恢复矢量的元素值更新为右移恢复矢量的元素值；
(5.11h)按照下式，更新目标到达角估计矢量：
Figure FDA00020950808900001415
其中，
Figure FDA00020950808900001416
表示更新后的目标到达角估计矢量，Θ(R)表示右移角度搜索范围，
Figure FDA00020950808900001417
表示右移峰值位置矢量
Figure FDA00020950808900001418
的第s个元素，s表示动态信源数；
(5.12)判断可调栅格步长Δ是否大于栅格优化精度ξθ，若是，则执行步骤(5.9)，否则，执行步骤(5.13)；
(5.13)按照下式，更新目标到达角估计矢量：
Figure FDA0002095080890000151
其中，
Figure FDA0002095080890000152
表示更新后的目标到达角估计矢量的第s个元素值，s表示动态信源数，
Figure FDA0002095080890000153
表示更新后的目标到达角估计矢量，
Figure FDA0002095080890000154
表示当前峰值对应的位置下标；
(5.14)判断动态信源数是否小于信源数，若是，则执行步骤(5.6)，否则，执行步骤(5.15)；
(5.15)获得第二次目标到达角估计矢量：
将步骤(5.13)中更新后的目标到达角矢量作为第二次的目标到达角估计矢量。
CN201910516081.3A 2019-06-14 2019-06-14 基于秩一降维模型和块矩阵恢复的波达方向估计方法 Active CN110174657B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910516081.3A CN110174657B (zh) 2019-06-14 2019-06-14 基于秩一降维模型和块矩阵恢复的波达方向估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910516081.3A CN110174657B (zh) 2019-06-14 2019-06-14 基于秩一降维模型和块矩阵恢复的波达方向估计方法

Publications (2)

Publication Number Publication Date
CN110174657A CN110174657A (zh) 2019-08-27
CN110174657B true CN110174657B (zh) 2022-12-27

Family

ID=67698455

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910516081.3A Active CN110174657B (zh) 2019-06-14 2019-06-14 基于秩一降维模型和块矩阵恢复的波达方向估计方法

Country Status (1)

Country Link
CN (1) CN110174657B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111814096B (zh) * 2020-06-28 2023-10-20 海南大学 基于子空间拟合的加权块稀疏恢复的mimo雷达定位方法
CN113255579B (zh) * 2021-06-18 2021-09-24 上海建工集团股份有限公司 一种施工监测异常采集数据自动识别与处理的方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103399291A (zh) * 2013-07-22 2013-11-20 西安电子科技大学 基于快速稀疏恢复的超分辨波达方向估计方法
CN103954950A (zh) * 2014-04-25 2014-07-30 西安电子科技大学 一种基于样本协方差矩阵稀疏性的波达方向估计方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103399291A (zh) * 2013-07-22 2013-11-20 西安电子科技大学 基于快速稀疏恢复的超分辨波达方向估计方法
CN103954950A (zh) * 2014-04-25 2014-07-30 西安电子科技大学 一种基于样本协方差矩阵稀疏性的波达方向估计方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于协方差矩阵降维稀疏表示的二维波达方向估计;李文杰等;《计算机应用》;20160810(第08期);第2197-2201页 *

Also Published As

Publication number Publication date
CN110174657A (zh) 2019-08-27

Similar Documents

Publication Publication Date Title
CN110244272B (zh) 基于秩一去噪模型的波达方向估计方法
Zhou et al. Direction-of-arrival estimation for coprime array via virtual array interpolation
CN110174658B (zh) 基于秩一降维模型和矩阵补全的波达方向估计方法
CN112363110B (zh) 一种基于嵌套交叉偶极子阵列的无网格单比特doa估计方法
CN109116293B (zh) 一种基于离格稀疏贝叶斯的波达方向估计方法
CN110109050B (zh) 嵌套阵列下基于稀疏贝叶斯的未知互耦的doa估计方法
CN111337893A (zh) 一种基于实值稀疏贝叶斯学习的离格doa估计方法
CN107544051A (zh) 嵌套阵列基于k‑r子空间的波达方向估计方法
Wu et al. Direction-of-arrival estimation based on Toeplitz covariance matrix reconstruction
CN105335615B (zh) 一种低复杂度的二维角度和极化参数联合估计方法
Bohme Source-parameter estimation by approximate maximum likelihood and nonlinear regression
CN104375133B (zh) 一种空间二维doa的估算方法
CN111239678A (zh) 一种基于l型阵列的二维doa估计方法
CN110174657B (zh) 基于秩一降维模型和块矩阵恢复的波达方向估计方法
KR101958337B1 (ko) 신호의 도래각을 추정하는 방법 및 장치
CN114720938A (zh) 基于深度展开的大规模天线阵列单比特采样doa估计方法
CN110941980B (zh) 一种密集环境中基于压缩感知的多径时延估计方法及装置
Feng et al. An off-grid iterative reweighted approach to one-bit direction of arrival estimation
Sahnoun et al. Tensor polyadic decomposition for antenna array processing
CN107544050A (zh) 白噪声背景下一种构造自适应阈值估计信号源数目的方法
CN112731280A (zh) 互质阵列混合噪声环境下的esprit-doa估计方法
Yang et al. A correlation-aware sparse Bayesian perspective for DOA estimation with off-grid sources
CN107135026B (zh) 未知互耦存在时基于矩阵重构的稳健波束形成方法
CN115421119A (zh) 一种基于收发翻转互质mimo雷达结构的doa估计方法
CN114844544A (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