CN110174659B - 基于迭代近端投影的mimo雷达多测量矢量doa估计方法 - Google Patents

基于迭代近端投影的mimo雷达多测量矢量doa估计方法 Download PDF

Info

Publication number
CN110174659B
CN110174659B CN201910541040.XA CN201910541040A CN110174659B CN 110174659 B CN110174659 B CN 110174659B CN 201910541040 A CN201910541040 A CN 201910541040A CN 110174659 B CN110174659 B CN 110174659B
Authority
CN
China
Prior art keywords
matrix
mimo radar
sparse
doa estimation
function
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
CN201910541040.XA
Other languages
English (en)
Other versions
CN110174659A (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.)
Nanjing University of Information Science and Technology
Original Assignee
Nanjing University of Information Science and Technology
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 Nanjing University of Information Science and Technology filed Critical Nanjing University of Information Science and Technology
Publication of CN110174659A publication Critical patent/CN110174659A/zh
Application granted granted Critical
Publication of CN110174659B publication Critical patent/CN110174659B/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
    • G01S3/00—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received
    • G01S3/02—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received using radio waves
    • G01S3/14—Systems for determining direction or deviation from predetermined direction
    • G01S3/143—Systems for determining direction or deviation from predetermined direction by vectorial combination of signals derived from differently oriented antennae
    • 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
    • G01S7/418—Theoretical aspects
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T10/00—Road transport of goods or passengers
    • Y02T10/10—Internal combustion engine [ICE] based vehicles
    • Y02T10/40—Engine management systems

Landscapes

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

Abstract

本发明公开了一种基于迭代近端投影的MIMO雷达多测量矢量DOA估计方法，通过高维回波数据转换至低维空间以降低空域维度，并对降维后的数据进行奇异值分解，提取信号子空间以降低时域维度，利用近端函数优化模型来表示MIMO雷达多测量矢量DOA估计中的非凸非平滑稀疏优化问题，然后在迭代过程中通过外推步骤和SCAD函数获得近端算子以求解该优化问题。本发明方法在低快拍和低信噪比下相干信源的DOA估计性能优于现有算法。

Description

基于迭代近端投影的MIMO雷达多测量矢量DOA估计方法
技术领域
本发明涉及多输入多输出雷达目标参数估计方法，特别是涉及一种基于迭代近端投影的MIMO雷达多测量矢量DOA估计方法。
背景技术
多输入多输出(Multiple Input and Multiple Output,MIMO)雷达是基于MIMO通信技术发展起来的一种新体制雷达。相比于传统相控阵雷达，MIMO雷达在目标探测、抗干扰、目标参数估计和目标识别等方面显示出诸多优势，因此已成为学术界研究的热点。
波达方向(Direction of Arrival,DOA)估计是MIMO雷达参数估计的一个重要内容。如今关于MIMO雷达DOA估计的算法已经多不胜数，诸如Capon算法，多重信号分类(Multiple Signal Classification,MUSIC)算法，基于旋转不变技术的信号参数估计(Estimation of Signal Parameters via Rotational Invariance Technique,ESPRIT)算法等。
通常，目标相对整个观测空间高度稀疏，因此许多学者将压缩感知理论应用到MIMO雷达DOA估计中。然而，求解l0范数最小化问题属于NP-hard问题，需要组合搜索、当维度增加时难以实现。Malioutov等人在论文“A sparse signal reconstructionperspective for source localization with sensor arrays”(IEEE Transactions onSignal Processing,2005,53(8):3010-3022)中针对传统阵列提出一种基于l1-SVD(l1norm-Singular Value Decomposition)的DOA估计方法，该方法通过对阵列接收数据进行奇异值分解，并建立信号子空间构建l2,1范数联合稀疏模型，然后采用二阶锥规划(Second Order Cone Programming,SOCP)求解该模型。Liu等人在论文“Direction ofarrival estimation via reweighted l1norm penalty algorithm for monostaticMIMO radar”(Multidimensional Systems and Signal Processing,2018,29(2):733-744)中提出一种基于加权l1范数算法的单基地MIMO雷达DOA估计方法，该算法利用降维Capon(Reduced-Dimensional Capon,RD-Capon)空间谱的系数构造加权矩阵，对小稀疏向量加以较大权重惩罚，对大稀疏向量加以较小权重惩罚以促进解的稀疏性，但由于采用SOCP求解方法导致其计算复杂度较高。Mohimani等人在论文“A Fast Approach forOvercomplete Sparse Decomposition Based on Smoothed L0Norm”(IEEE Transactionson Signal Processing,2009,57(1):289-301)中提出一种平滑l0范数(Smoothed l0norm，SL0)算法，该方法利用梯度投影原理以及最速下降法，采用一个平滑的高斯函数近似l0范数，从而将离散函数的优化问题转化为连续函数的优化问题。该算法的运算效率较高，在保证相同精度的条件下，能够比基追踪算法的重构速度快2～3倍。为了进一步提高目标的重构性能，Jing等人在论文“Joint Smoothed-Norm DOA Estimation Algorithm forMultiple Measurement Vectors in MIMO Radar”(Sensors,2017,17(5):1068)中提出一种基于平滑l0范数的MIMO雷达多测量矢量(Multiple Measurement Vector，MMV)DOA估计方法，利用接收数据的四阶累积量矩阵构建稀疏表示模型以有效抑制噪声，并将MMV问题转化为一个联合平滑函数的求解问题，且该方法相比l1-SVD方法快约两个数量级。Sadeghi等人在论文“Iterative Sparsification-Projection:Fast and Robust Sparse SignalApproximation”(IEEE Transactions on Signal Processing,2016,64(21):5536-5548)中基于近端方法提出一种相对SL0算法更为广泛的稀疏重构算法以解决非凸平滑问题，称为迭代稀疏投影(Iterative Sparsification-Projection,ISP)算法。针对非凸非平滑优化问题，Ghayem等人在论文“Sparse Signal Recovery Using Iterative ProximalProjection”(IEEE Transactions on Signal Processing,2018,66(4):879-894)中提出一种迭代近端投影(Iterative Proximal Projection,IPP)算法，该算法通过构造近端投影优化模型，并利用SCAD(Smoothly Clipped Absolute Deviation Penalty)阈值函数获得近端算子来求解该模型以促进解的稀疏度，同时引用外推步骤改善该算法的收敛性能。传统基于压缩感知的MIMO雷达DOA估计算法将非凸非平滑稀疏表示问题近似成凸或平滑函数问题进行求解，因此在稀疏表示模型上存在一定程度上的误差，从而导致传统基于压缩感知的DOA估计算法的性能不能达到最优。
发明内容
发明目的：本发明解决的技术问题是针对基于压缩感知的MIMO雷达DOA估计问题，提供一种能有效解决非凸非平滑稀疏表示问题的方法以提高DOA估计性能。
技术方案：本发明提供了一种基于迭代近端投影的MIMO雷达多测量矢量DOA估计方法，包括如下步骤：
步骤(1)：对MIMO雷达接收阵列信号进行匹配滤波，取多个快拍下虚拟阵列的输出信号矩阵；
步骤(2)：对MIMO雷达虚拟阵列输出信号矩阵进行降维变换，得到降维后的接收数据矩阵；
步骤(3)：对MIMO雷达降维后的输出信号进行奇异值分解，得到由多测量矢量构成的矩阵；
步骤(4)：根据稀疏重构理论，将搜索空域按等角度间隔划分，将步骤(3)得到的矩阵转换成稀疏表示模型；
步骤(5)：利用近端函数模型建立MIMO雷达多测量矢量DOA估计的稀疏优化问题，在迭代过程中通过外推步骤和SCAD函数获得近端算子以求解该优化问题；
步骤(6)：步骤(5)获得稀疏解之后，通过搜索其谱峰所在位置得到真实目标DOA估计值。
进一步地，步骤(1)具体为：对具有M个发射阵元和N个接收阵元的MIMO雷达接收阵列信号进行匹配滤波，取J个快拍下MIMO雷达虚拟阵列的输出信号X＝AS+N，其中
Figure BDA0002102569020000031
为接收信号矩阵；/>
Figure BDA0002102569020000032
为信号矩阵，其中/>
Figure BDA0002102569020000033
表示复数域；/>
Figure BDA0002102569020000034
为高斯噪声矩阵；
Figure BDA0002102569020000035
为发射接收联合导向矩阵，其中
Figure BDA0002102569020000036
为对应第p(p＝1,2,...,P)个目标的发射阵列的导向向量，/>
Figure BDA0002102569020000037
为对应第p(p＝1,2,...,P)个目标的接收阵列的导向向量，(·)T表示矩阵转置，θp为第p(p＝1,2,...,P)个目标的方位角，/>
Figure BDA0002102569020000038
表示Kronecker积，P是相干目标的数目。
步骤(2)具体为：对MIMO雷达虚拟阵列输出信号矩阵X进行降维变换，得到降维后的接收数据矩阵
Figure BDA0002102569020000039
其中，/>
Figure BDA00021025690200000310
为降维矩阵，/>
Figure BDA00021025690200000311
Figure BDA00021025690200000312
为转换矩阵，/>
Figure BDA00021025690200000313
Figure BDA00021025690200000314
0N×M为N×M维的零矩阵；(·)H表示共轭转置运算
步骤(3)具体为：为了进一步降低计算复杂度和对噪声的敏感性，对MIMO雷达降维后的输出信号矩阵
Figure BDA00021025690200000319
进行奇异值分解，得到/>
Figure BDA00021025690200000315
其中USV由P个大特征值对应的左奇异值矢量组成的信号子空间矩阵，UNV由其余M+N-1-P个小特征值对应的左奇异特征值矢量组成的噪声子空间矩阵，V为右奇异特征值矢量组成的矩阵，Λ为
Figure BDA00021025690200000316
的特征值构成的对角矩阵。令/>
Figure BDA00021025690200000317
则/>
Figure BDA00021025690200000318
为由多测量矢量构成的矩阵，可表示为
Figure BDA0002102569020000041
其中/>
Figure BDA0002102569020000042
为降维后的阵列流形矩阵，
Figure BDA0002102569020000043
为虚拟均匀线阵导向矩阵，SSV＝SVDP，
Figure BDA0002102569020000044
ΛP×P由P个大特征值组成的对角矩阵，0P×(J-P)为P×(J-P)维的零矩阵。
步骤(4)具体为：根据稀疏重构理论，将搜索空域[-90°,90°]按等角度间隔划分为L个单元，且L＞＞P，
Figure BDA0002102569020000045
表示空域内所有可能的入射方向，定义冗余字典
Figure BDA0002102569020000046
其中/>
Figure BDA0002102569020000047
则/>
Figure BDA0002102569020000048
又可转换成稀疏表示模型：/>
Figure BDA0002102569020000049
其中,/>
Figure BDA00021025690200000410
与SSV具有相同的行支撑，即Sθ是P行稀疏矩阵，Sθ中的非零行元素对应冗余字典中目标的DOA，因此DOA估计问题可以转化为求解稀疏矩阵Sθ中非零行元素位置的问题。
步骤(5)具体为：利用近端函数优化模型建立MIMO雷达多测量矢量DOA估计的稀疏优化问题：
Figure BDA00021025690200000411
其中，/>
Figure BDA00021025690200000412
是非凸非平滑函数，z为辅助变量，/>
Figure BDA00021025690200000413
是由矩阵z的每一行向量的l2范数构成的列向量，/>
Figure BDA00021025690200000414
是由矩阵Sθ的每一行向量的l2范数构成的列向量，/>
Figure BDA00021025690200000415
定义为可行集/>
Figure BDA00021025690200000416
的指示函数。该稀疏优化问题需要多次迭代来求解，其中在第k次迭代中的稀疏解可表示为/>
Figure BDA00021025690200000417
其中
Figure BDA00021025690200000418
为非凸非平滑函数/>
Figure BDA00021025690200000419
的近端算子，/>
Figure BDA00021025690200000420
代表可行集/>
Figure BDA00021025690200000421
的投影。
步骤(6)：在获得稀疏解
Figure BDA00021025690200000422
之后，通过搜索其谱峰所在位置得到真实目标DOA估计值。
进一步地：步骤(5)中的通过多次迭代来求解稀疏优化问题的具体步骤为：
步骤(a)：为了避免迭代近端投影算法在求解上述非凸非平滑稀疏优化问题时易陷入局部最小值的问题，利用外推步骤改善迭代近端投影算法性能，则在第k次迭代中的稀疏解表示为
Figure BDA00021025690200000423
其中w≥0为权重常数。该模型在第k次迭代中的稀疏解可进一步表示为/>
Figure BDA0002102569020000051
其中
Figure BDA0002102569020000052
为非凸非平滑函数/>
Figure BDA0002102569020000053
的近端算子。
步骤(b)：由于非平滑函数的近端算子
Figure BDA0002102569020000054
可通过SCAD(Smoothly Clipped Absolute Deviation Penalty)惩罚函数
Figure BDA0002102569020000055
产生相应的SCAD阈值函数
Figure BDA0002102569020000056
来计算，其中，/>
Figure BDA0002102569020000057
为矢量/>
Figure BDA0002102569020000058
中第l(l＝1,2,…,L)个元素，λ为调整参数，a为常量，一般取值为a＞2，sign(·)为符号函数，(α)+＝max(α，0)。
发明原理：由于传统基于压缩感知的MIMO雷达DOA估计算法通常将非凸非平滑稀疏优化问题近似成凸或平滑函数优化问题进行求解，由于稀疏优化问题模型误差的存在而导致DOA估计性能不理想。因此，本发明提出一种基于迭代近端投影的MIMO雷达多测量矢量DOA估计方法。该方法首先将高维回波数据转换至低维空间以降低空域维度，并对降维后的数据进行奇异值分解，提取信号子空间以降低时域维度，利用近端函数优化模型来表示MIMO雷达多测量矢量DOA估计中的非凸非平滑稀疏优化问题，然后在迭代过程中通过外推步骤和SCAD函数获得近端算子以求解该模型。本发明方法在低快拍和低信噪比下相干信源的DOA估计性能优于现有算法。
有益效果：与现有技术相比
(1)基于压缩感知的MIMO雷达DOA估计优化问题是一种非凸非平滑稀疏优化问题，由于该问题求解难度较大，因此传统方法通常将此问题近似成凸或平滑函数优化问题进行求解，而本发明利用了近端函数优化模型较好地表示了该非凸非平滑稀疏优化问题，尽可能地降低稀疏优化问题模型的误差，从而提高了MIMO雷达DOA估计性能；
(2)本发明采用SCAD函数获得近端算子来求解以近端函数模型所表示的非凸非平滑稀疏优化问题，不仅克服了硬阈值收缩函数对数据中微小波动的敏感性，而且避免了软阈值收缩函数带来的偏差，因此能够进一步促进解的稀疏性，有效提高了本发明算法对MIMO雷达DOA估计的精度；
(3)本发明方法将高维回波数据转换至低维空间以降低空域维度，并对降维后的数据进行奇异值分解(SVD)，提取信号子空间以降低时域维度，这样不仅提高了算法的实时性，而且降低了算法对噪声的敏感性。
附图说明
图1是本发明的流程图；
图2是不同算法的DOA估计均方根误差随信噪比变化的曲线图；
图3是不同算法的DOA估计均方根误差随快拍数的变化关系；
图4是不同算法的DOA估计均方根误差随角度间隔变化关系。
具体实施方式
下面结合实施例对本发明作更进一步的说明。
如图1所示，本发明的基于迭代近端投影的MIMO雷达多测量矢量DOA估计方法包括如下步骤：
步骤1：对具有M个发射阵元和N个接收阵元的MIMO雷达接收阵列信号进行匹配滤波，取J个快拍下MIMO雷达虚拟阵列的输出信号X＝AS+N，其中
Figure BDA0002102569020000061
为接收信号矩阵；/>
Figure BDA0002102569020000062
为信号矩阵，其中/>
Figure BDA0002102569020000063
表示复数域；
Figure BDA0002102569020000064
为高斯噪声矩阵；/>
Figure BDA0002102569020000065
为发射接收联合导向矩阵，其中/>
Figure BDA0002102569020000066
为对应第p(p＝1,2,...,P)个目标的发射阵列的导向向量，/>
Figure BDA0002102569020000067
为对应第p(p＝1,2,...,P)个目标的接收阵列的导向向量，(·)T表示矩阵转置，θp为第p(p＝1,2,...,P)个目标的方位角，/>
Figure BDA0002102569020000068
表示Kronecker积，P是相干目标的数目。
假设窄带单基地MIMO雷达系统具有M个发射阵元和N个接收阵元，发射和接收阵列均为均匀线阵，其阵元间隔分别为dt＝λ/2和dr＝λ/2，λ为接收信号波长。假设存在P个远场窄带相干目标，其入射角度分别为θ1,θ2,...,θP，MIMO雷达接收阵列信号经匹配滤波后可表示为：
x(t)＝As(t)+n(t) (1)
式中，
Figure BDA0002102569020000071
为发射接收联合导向矩阵；
Figure BDA0002102569020000072
为对应第p(p＝1,2,...,P)个目标的发射阵列的导向向量；/>
Figure BDA0002102569020000073
为对应第p(p＝1,2,...,P)个目标的接收阵列的导向向量；(·)T表示矩阵转置；/>
Figure BDA0002102569020000074
表示复数域；/>
Figure BDA0002102569020000075
表示Kronecker积；
Figure BDA0002102569020000076
为相干信源信号矢量，其中，sp(t)＝αps1(t)，复常数αp表示sp(t)相对于s1(t)的相干系数；/>
Figure BDA0002102569020000077
为接收阵列的噪声矩阵，服从零均值，方差为σ2的高斯分布，即n(t)～Nc(0,σ2IMN)，IMN表示MN×MN维的单位矩阵。取J个快拍下MIMO雷达虚拟阵列的输出信号，即
X＝AS+N (2)
式中，
Figure BDA0002102569020000078
为接收信号矩阵；/>
Figure BDA0002102569020000079
为信号矩阵；/>
Figure BDA00021025690200000710
为高斯噪声矩阵。
步骤2：对MIMO雷达虚拟阵列输出信号矩阵X进行降维变换，得到降维后的接收数据矩阵
Figure BDA00021025690200000711
其中，/>
Figure BDA00021025690200000712
为降维矩阵，/>
Figure BDA00021025690200000713
Figure BDA00021025690200000714
为转换矩阵，/>
Figure BDA00021025690200000715
Figure BDA00021025690200000716
0N×M为N×M维的零矩阵；(·)H表示共轭转置运算。
对于发射阵列和接收阵列均为阵元间距等于半波长的均匀线阵的单基地MIMO雷达来说，其有效虚拟阵元个数为M+N-1，因此，MN×1维的目标均匀线阵导向矢量可由(M+N-1)×1维的虚拟均匀线阵导向矢量通过线性变换来表示，即
Figure BDA00021025690200000717
式中，
Figure BDA0002102569020000081
表示(M+N-1)×1维的虚拟均匀线阵导向矢量，/>
Figure BDA0002102569020000082
为转换矩阵，其中，
Figure BDA0002102569020000083
Figure BDA0002102569020000084
0N×M为N×M维的零矩阵。根据式(3)，阵列流形矩阵A可进一步表示为
A＝GB (4)
式中，B＝[b(θ1),b(θ2),...,b(θP)]为(M+N-1)×P维的虚拟均匀线阵导向矩阵。为了降低算法的计算复杂度，可定义一个大小为(M+N-1)×MN的降维矩阵对接收数据进行降维预处理。为了使得降维后的噪声是服从N～(0,σ2IM+N-1)的高斯白噪声，降维矩阵T需满足TTH＝IM+N-1，因此降维矩阵可选取为
Figure BDA0002102569020000085
其中
Figure BDA0002102569020000086
式中，min(·)表示取最小的元素；diag(·)表示对角化操作；(·)H表示共轭转置运算。经过降维变换后，式(2)的接收数据矩阵可表示为：
Figure BDA0002102569020000087
/>
式中，
Figure BDA0002102569020000088
为降维后的接收数据矩阵；/>
Figure BDA0002102569020000089
为降维后的阵列流形矩阵；/>
Figure BDA00021025690200000810
为降维后的高斯白噪声矩阵。
步骤3：为了进一步降低计算复杂度和对噪声的敏感性，对MIMO雷达降维后的输出信号矩阵
Figure BDA00021025690200000819
进行奇异值分解，得到/>
Figure BDA00021025690200000811
其中USV由P个大特征值对应的左奇异值矢量组成的信号子空间矩阵，UNV由其余M+N-1-P个小特征值对应的左奇异特征值矢量组成的噪声子空间矩阵，V为右奇异特征值矢量组成的矩阵，Λ为/>
Figure BDA00021025690200000812
的特征值构成的对角矩阵。令/>
Figure BDA00021025690200000813
则/>
Figure BDA00021025690200000814
为由多测量矢量构成的矩阵，可表示为
Figure BDA00021025690200000815
其中/>
Figure BDA00021025690200000816
为降维后的阵列流形矩阵，
Figure BDA00021025690200000817
为虚拟均匀线阵导向矩阵，SSV＝SVDP，/>
Figure BDA00021025690200000818
Figure BDA0002102569020000091
ΛP×P由P个大特征值组成的对角矩阵，0P×(J-P)为P×(J-P)维的零矩阵。
为了进一步降低计算复杂度和对噪声的敏感性，对接收数据矩阵
Figure BDA00021025690200000917
进行奇异值分解(SVD)，
Figure BDA0002102569020000092
其中，USV由P个大特征值对应的左奇异值矢量组成的信号子空间矩阵，UNV由其余M+N-1-P个小特征值对应的左奇异特征值矢量组成的噪声子空间矩阵，V为右奇异特征值矢量组成的矩阵，Λ为X的特征值组成的对角矩阵。令
Figure BDA0002102569020000093
进一步推导式(7)得：
Figure BDA0002102569020000094
其中，
Figure BDA0002102569020000095
ΛP×P由P个大特征值组成的对角矩阵，0P×(J-P)为P×(J-P)维的零矩阵。令SSV＝SVDP，/>
Figure BDA0002102569020000096
则降维后的接收数据模型(8)可简化为：
Figure BDA0002102569020000097
对比式(6)和式(9)可知，阵列接收数据矩阵的维数从(M+N-1)×J降到(M+N-1)×P，而实际应用中，目标个数P远远小于快拍数J，即P＜＜J，因此，式(9)的信号矩阵维度得到显著降低，有利于目标DOA估计的快速实现。
步骤4：根据稀疏重构理论，将搜索空域[-90°,90°]按等角度间隔划分为L个单元，且L＞＞P，
Figure BDA0002102569020000098
表示空域内所有可能的入射方向，定义冗余字典/>
Figure BDA0002102569020000099
其中/>
Figure BDA00021025690200000910
则/>
Figure BDA00021025690200000911
又可转换成稀疏表示模型：/>
Figure BDA00021025690200000912
其中,/>
Figure BDA00021025690200000913
与SSV具有相同的行支撑，即Sθ是P行稀疏矩阵，Sθ中的非零行元素对应冗余字典中目标的DOA，因此DOA估计问题可以转化为求解稀疏矩阵Sθ中非零行元素位置的问题。
根据稀疏重构理论，将搜索空域[-90°,90°]按等角度间隔划分为L个单元，且L＞＞P，
Figure BDA00021025690200000914
表示空域内所有可能的入射方向，定义冗余字典/>
Figure BDA00021025690200000915
其中
Figure BDA00021025690200000916
则式(9)可转换成稀疏表示模型：/>
Figure BDA0002102569020000101
其中,
Figure BDA0002102569020000102
与SSV具有相同的行支撑，即Sθ是P行稀疏矩阵，Sθ中的非零行元素对应冗余字典中目标的DOA，因此DOA估计问题可以转化为求解稀疏矩阵Sθ中非零行元素位置的问题。从多个测量矢量/>
Figure BDA0002102569020000103
中恢复多个未知的稀疏信号源Sθ，该重建模型(10)可以称为多测量矢量(MMV)模型。该模型采用多快拍下的阵列接收数据，提高了在低信噪比下DOA估计性能。
步骤5：利用近端函数优化模型建立MIMO雷达多测量矢量DOA估计的稀疏优化问题：
Figure BDA0002102569020000104
其中，/>
Figure BDA0002102569020000105
是非凸非平滑函数，z为辅助变量，/>
Figure BDA0002102569020000106
是由矩阵z的每一行向量的l2范数构成的列向量，/>
Figure BDA0002102569020000107
是由矩阵Sθ的每一行向量的l2范数构成的列向量，/>
Figure BDA0002102569020000108
定义为可行集/>
Figure BDA0002102569020000109
的指示函数。该稀疏优化问题需要多次迭代来求解，其中在第k次迭代中的稀疏解可表示为/>
Figure BDA00021025690200001010
其中
Figure BDA00021025690200001011
为非凸非平滑函数/>
Figure BDA00021025690200001012
的近端算子，/>
Figure BDA00021025690200001013
代表可行集/>
Figure BDA00021025690200001014
的投影。
为获得式(10)的稀疏解，可将MIMO雷达的多测量矢量DOA估计问题转化为近端函数优化问题，其稀疏优化问题为：
Figure BDA00021025690200001015
其中，
Figure BDA00021025690200001016
是非凸非平滑函数，z为辅助变量，/>
Figure BDA00021025690200001017
是由矩阵z的每一行向量的l2范数构成的列向量，/>
Figure BDA00021025690200001018
是由矩阵Sθ的每一行向量的l2范数构成的列向量。/>
Figure BDA00021025690200001019
定义为可行集
Figure BDA00021025690200001020
的指示函数，即：
Figure BDA00021025690200001021
在式(11)的稀疏优化模型中引入惩罚函数方法，即：
Figure BDA00021025690200001022
式中，λ＞0为惩罚因子，||·||2表示l2范数。通过交替最小化方法将式(13)转化为关于
Figure BDA00021025690200001023
和/>
Figure BDA00021025690200001024
的两个子问题的迭代求解，即：
Figure BDA0002102569020000111
根据近端投影的定义，式(14)进一步化简为：
Figure BDA0002102569020000112
其中，
Figure BDA0002102569020000113
为非平滑函数/>
Figure BDA0002102569020000114
的近端算子。该模型第k次迭代中的稀疏解可表示为：
Figure BDA0002102569020000115
其中，
Figure BDA0002102569020000116
代表可行集/>
Figure BDA0002102569020000117
的投影。
为了避免在求解非凸非平滑问题时易陷入局部最小值的问题，利用外推步骤改善迭代近端投影(IPP)算法的收敛性能，则在第k次迭代中的稀疏解表示为：
Figure BDA0002102569020000118
其中，w≥0为权重常数。则式(13)在第k次迭代中的的稀疏解可进一步表示为：
Figure BDA0002102569020000119
其中，
Figure BDA00021025690200001110
为非凸非平滑函数/>
Figure BDA00021025690200001111
的近端算子。
非平滑函数的近端算子
Figure BDA00021025690200001112
可通过不同的惩罚函数/>
Figure BDA00021025690200001113
产生相应的阈值函数获得。采用l0惩罚函数来表示/>
Figure BDA00021025690200001114
即/>
Figure BDA00021025690200001115
则非平滑函数的近端算子可表示为
Figure BDA00021025690200001116
式中，硬阈值函数
Figure BDA00021025690200001117
可表示为
Figure BDA00021025690200001118
式中，λ为调整参数。采用l1惩罚函数来表示
Figure BDA00021025690200001119
即/>
Figure BDA00021025690200001120
则非平滑函数的近端算子可表示为：
Figure BDA0002102569020000121
式中，软阈值函数
Figure BDA0002102569020000122
可表示为：
Figure BDA0002102569020000123
/>
当采用SCAD惩罚函数来表示
Figure BDA0002102569020000124
即
Figure BDA0002102569020000125
式中，a为常量，一般取值为a＞2。则近端算子
Figure BDA0002102569020000126
可采用SCAD阈值函数
Figure BDA0002102569020000127
来表示，即
Figure BDA0002102569020000128
式中，SCAD阈值函数
Figure BDA0002102569020000129
可表示为
Figure BDA00021025690200001210
式中，sign(·)为符号函数，(α)+＝max(α，0)。该SCAD阈值函数
Figure BDA00021025690200001211
不仅克服了硬阈值收缩函数/>
Figure BDA00021025690200001212
对数据中微小波动的敏感性，而且避免了软阈值收缩函数/>
Figure BDA00021025690200001213
带来的偏差，因此能够进一步促进解的稀疏性。
在步骤5中采用外推步骤和SCAD阈值函数
Figure BDA00021025690200001214
来求解稀疏优化问题的具体步骤如下：
初始化：
(a)设初始值
Figure BDA00021025690200001215
(b)选取一组合适的序列[λ1,λ2,...,λK]，且λk+1＝cλk，0＜c＜1，λ1＝max{|u0|}，a＝2000，w＝0.9
算法迭代：
Fork＝1,2,...,K,其中K为外循环次数
(a)令λ＝λk
(b)进行Q次迭代求解全局最小值，并将该最小值投影到可行集上
(1)令
Figure BDA0002102569020000131
(2)Forq＝1,2,...,Q,其中Q为内循环次数
a)计算
Figure BDA0002102569020000132
/>
b)计算
Figure BDA0002102569020000133
c)将
Figure BDA0002102569020000134
投影到可行集，/>
Figure BDA0002102569020000135
Figure BDA0002102569020000136
代表可行集/>
Figure BDA0002102569020000137
的投影
(3)令
Figure BDA0002102569020000138
最终解
Figure BDA0002102569020000139
步骤6：在获得稀疏解
Figure BDA00021025690200001310
之后，通过搜索其谱峰所在位置得到真实目标DOA估计值。
本发明的技术效果可通过以下仿真进一步说明，为了验证本发明方法在MIMO雷达DOA估计方面的优势，选取MUSIC算法、前后向空域平滑MUSIC(Forward and BackwardSpatial Smoothing_Music,FBSS_MUSIC)算法、SL0(Smoothedl0norm,SL0)算法与l1-SVD(l1-norm Singular Value Decomposition,l1-SVD)算法进行对比，将本发明称为IPP_SCAD_SVD算法。假设均匀线阵单基地MIMO雷达的发射和接收阵元数均为5，收发阵元间隔为dt＝dr＝λ/2，在空间角度范围[-90°,90°]以角度间隔0.05°等分。假设存在3个远场窄带相干目标，设置各目标的DOA分别为θ1＝-19.8°，θ2＝0°，θ3＝20.6°。DOA估计的均方根误差定义为
Figure BDA0002102569020000141
其中，/>
Figure BDA0002102569020000142
表示第p个目标在第mt次蒙特卡罗实验中的目标DOA估计值，MT为蒙特卡罗实验次数。信噪比定义为/>
Figure BDA0002102569020000143
在本发明方法中，设置SCAD阈值函数/>
Figure BDA0002102569020000144
中的参数a＝2000，权重常数w＝0.9，外循环次数K＝3，衰减因子c＝0.8。
仿真实验1：图1为各种算法的DOA估计均方根误差随信噪比变化的曲线图。设置信噪比在-10～15dB之间变化，快拍数J＝100，进行100次蒙特卡罗实验。从图1中可以看出，由于信源相干导致协方差矩阵的秩亏缺，信号特征向量发散到噪声子空间，造成MUSIC算法的DOA估计方法失效；FBSS_MUSIC算法利用前后向空间平滑技术实现信号的解相干，能够有效估计DOA，但在高信噪比下该方法的DOA估计性能较差。SL0_SVD算法、L1_SVD算法以及IPP_SCAD_SVD算法均属于基于压缩感知的DOA估计方法，这些方法能够对相干信源的DOA进行有效估计，且DOA估计的均方根误差随信噪比的增加而减小，而本发明算法的DOA估计精度明显优于其他算法。
仿真实验2：图2为各种算法的DOA估计均方根误差随快拍数的变化关系曲线。设置信噪比为-5dB，进行100次蒙特卡罗实验，快拍数J在50～350之间变化。从图2中可以看出，MUSIC算法无法处理相干信源，其他各种算法的角度估计精度均随着快拍数的增加而提高，而本发明算法相比其他算法具有更高的DOA估计精度。
仿真实验3：图3为各种算法的DOA估计均方根误差随两个相干目标的入射角度间隔Δθ的变化关系。两个相干目标的入射角度分别为θ1＝6°，θ2＝6°+Δθ，其中Δθ∈[14°22°]，信噪比为5dB，快拍数J＝100，进行100次蒙特卡罗实验。从图3中可以看出，MUSIC算法无法分辨相干目标，FBSS_MUSIC算法、SL0_SVD算法、L1_SVD算法和本发明算法的DOA估计精度均随着目标角度间隔的增大而提高，而本发明算法的DOA估计精度始终高于其他算法，表明本发明算法相比其他算法能获得更高的空间角度分辨率。

Claims (7)

1.一种基于迭代近端投影的MIMO雷达多测量矢量DOA估计方法，其特征在于，包括如下步骤：
步骤(1)：对MIMO雷达接收阵列信号进行匹配滤波，取多个快拍下虚拟阵列的输出信号矩阵；
步骤(2)：对MIMO雷达虚拟阵列输出信号矩阵进行降维变换，得到降维后的接收数据矩阵；
步骤(3)：对MIMO雷达降维后的输出信号进行奇异值分解，得到由多测量矢量构成的矩阵；
步骤(4)：根据稀疏重构理论，将搜索空域按等角度间隔划分，将步骤(3)得到的矩阵转换成稀疏表示模型；
步骤(5)：利用近端函数模型建立MIMO雷达多测量矢量DOA估计的稀疏优化问题，在迭代过程中通过外推步骤和SCAD函数获得近端算子以求解该优化问题；
步骤(6)：步骤(5)获得稀疏解之后，通过搜索其谱峰所在位置得到真实目标DOA估计值。
2.根据权利要求1所述的基于迭代近端投影的MIMO雷达多测量矢量DOA估计方法，其特征在于：步骤(1)具体为：
对具有M个发射阵元和N个接收阵元的MIMO雷达接收阵列信号进行匹配滤波，取J个快拍下虚拟阵列的输出信号矩阵X＝AS+N；其中
Figure FDA0004168421160000011
为接收信号矩阵；/>
Figure FDA0004168421160000012
为信号矩阵，其中/>
Figure FDA0004168421160000013
表示复数域；
Figure FDA0004168421160000014
为高斯噪声矩阵；/>
Figure FDA0004168421160000015
为发射接收联合导向矩阵，其中/>
Figure FDA0004168421160000016
为对应第p个目标的发射阵列的导向向量，/>
Figure FDA0004168421160000017
为对应第p个目标的接收阵列的导向向量，(·)T表示矩阵转置，θp为第p个目标的方位角，/>
Figure FDA0004168421160000018
表示Kronecker积，p＝1,2,...,P，P是相干目标的数目。
3.根据权利要求2所述的基于迭代近端投影的MIMO雷达多测量矢量DOA估计方法，其特征在于，步骤(2)具体为：
对MIMO雷达虚拟阵列输出信号矩阵X进行降维变换，得到降维后的接收数据矩阵
Figure FDA0004168421160000021
其中，/>
Figure FDA0004168421160000022
为降维矩阵，/>
Figure FDA0004168421160000023
Figure FDA0004168421160000024
为转换矩阵；
Figure FDA0004168421160000025
0N×M为N×M维的零矩阵；(·)H表示共轭转置运算。
4.根据权利要求3所述的基于迭代近端投影的MIMO雷达多测量矢量DOA估计方法，其特征在于，步骤(3)具体为：
对MIMO雷达降维后的输出信号
Figure FDA0004168421160000026
进行奇异值分解，得到/>
Figure FDA0004168421160000027
/>
Figure FDA0004168421160000028
其中USV由P个大特征值对应的左奇异值矢量组成的信号子空间矩阵，UNV由其余M+N-1-P个小特征值对应的左奇异特征值矢量组成的噪声子空间矩阵，V为右奇异特征值矢量组成的矩阵，Λ为/>
Figure FDA0004168421160000029
的特征值构成的对角矩阵；令/>
Figure FDA00041684211600000210
则/>
Figure FDA00041684211600000211
为由多测量矢量构成的矩阵，可表示为/>
Figure FDA00041684211600000212
其中/>
Figure FDA00041684211600000213
为降维后的阵列流形矩阵，
Figure FDA00041684211600000214
为虚拟均匀线阵导向矩阵，SSV＝SVDP，
Figure FDA00041684211600000215
ΛP×P由P个大特征值组成的对角矩阵，0P×(J-P)为P×(J-P)维的零矩阵。
5.根据权利要求4所述的基于迭代近端投影的MIMO雷达多测量矢量DOA估计方法，其特征在于，步骤(4)具体为：
根据稀疏重构理论，将搜索空域[-90°,90°]按等角度间隔划分为L个单元，且L＞＞P，
Figure FDA00041684211600000216
表示空域内所有可能的入射方向，定义冗余字典/>
Figure FDA00041684211600000217
其中
Figure FDA00041684211600000218
则/>
Figure FDA00041684211600000219
又可转换成稀疏表示模型：/>
Figure FDA00041684211600000220
其中,/>
Figure FDA00041684211600000221
与SSV具有相同的行支撑，即Sθ是P行稀疏矩阵，Sθ中的非零行元素对应冗余字典中目标的DOA。
6.根据权利要求5所述的基于迭代近端投影的MIMO雷达多测量矢量DOA估计方法，其特征在于，步骤(5)具体为：
利用近端函数优化模型建立MIMO雷达多测量矢量DOA估计的稀疏优化问题：
Figure FDA00041684211600000222
其中，/>
Figure FDA00041684211600000223
是非凸非平滑函数，z为辅助变量，/>
Figure FDA00041684211600000224
是由矩阵z的每一行向量的l2范数构成的列向量，/>
Figure FDA0004168421160000031
是由矩阵Sθ的每一行向量的l2范数构成的列向量，
Figure FDA0004168421160000032
定义为可行集/>
Figure FDA0004168421160000033
的指示函数；该稀疏优化问题需要多次迭代求解，其中在第k次迭代中的稀疏解可表示为/>
Figure FDA0004168421160000034
其中
Figure FDA0004168421160000035
为非凸非平滑函数/>
Figure FDA0004168421160000036
的近端算子，/>
Figure FDA0004168421160000037
代表可行集/>
Figure FDA0004168421160000038
的投影。
7.根据权利要求6所述的基于迭代近端投影的MIMO雷达多测量矢量DOA估计方法，其特征在于，步骤(5)中的通过多次迭代求解稀疏优化问题的具体步骤为：
步骤(a)：利用外推步骤改善算法性能，在第k次迭代中的稀疏解表示为
Figure FDA0004168421160000039
其中w≥0为权重常数；该模型在第k次迭代中的稀疏解可进一步表示为/>
Figure FDA00041684211600000310
其中/>
Figure FDA00041684211600000311
为非凸非平滑函数/>
Figure FDA00041684211600000312
的近端算子；
步骤(b)：非平滑函数的近端算子
Figure FDA00041684211600000313
可通过SCAD惩罚函数/>
Figure FDA00041684211600000314
产生相应的SCAD阈值函数/>
Figure FDA00041684211600000315
来计算，其中，/>
Figure FDA00041684211600000316
为矢量/>
Figure FDA00041684211600000317
中第l个元素，l＝1,2,…,L，λ为调整参数，a为常量，取值为a＞2，sign(·)为符号函数，(α)+＝max(α，0)。/>
CN201910541040.XA 2019-05-08 2019-06-21 基于迭代近端投影的mimo雷达多测量矢量doa估计方法 Active CN110174659B (zh)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN201910382095 2019-05-08
CN2019103820950 2019-05-08

Publications (2)

Publication Number Publication Date
CN110174659A CN110174659A (zh) 2019-08-27
CN110174659B true CN110174659B (zh) 2023-05-23

Family

ID=67698839

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910541040.XA Active CN110174659B (zh) 2019-05-08 2019-06-21 基于迭代近端投影的mimo雷达多测量矢量doa估计方法

Country Status (1)

Country Link
CN (1) CN110174659B (zh)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111107023B (zh) * 2019-09-23 2024-02-02 南京邮电大学 在大规模mimo中基于光滑范数的压缩感知信道估计方法
CN111060885B (zh) * 2019-12-13 2023-06-20 河海大学 一种mimo雷达的参数估计方法
CN111142063B (zh) * 2020-01-06 2023-04-07 西安邮电大学 一种基于降维优化的快速压缩感知低空目标测角方法
CN111680262B (zh) * 2020-05-07 2023-07-14 北京控制工程研究所 一种自适应增益梯度投影辨识方法
CN112462352B (zh) * 2020-10-30 2022-10-18 哈尔滨工程大学 一种适用于低信噪比条件下的线谱增强方法
CN113835085B (zh) * 2021-09-30 2023-07-25 南京信息工程大学 一种基于复杂地形补偿的雷达快速测高方法
CN116362316B (zh) * 2023-05-29 2023-12-12 成都阿加犀智能科技有限公司 一种模型转换方法、装置、存储介质及电子设备

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101272168B (zh) * 2007-03-23 2012-08-15 中国科学院声学研究所 一种信源数估计方法及其波达方向估计方法
US20120259590A1 (en) * 2011-04-11 2012-10-11 Jong Chul Ye Method and apparatus for compressed sensing with joint sparsity
CN103954950B (zh) * 2014-04-25 2016-09-07 西安电子科技大学 一种基于样本协方差矩阵稀疏性的波达方向估计方法
CN103983958A (zh) * 2014-05-16 2014-08-13 哈尔滨工程大学 基于多测量矢量稀疏表示的mimo雷达连续目标角度估计方法
JP6395677B2 (ja) * 2015-08-10 2018-09-26 三菱電機株式会社 到来方向推定装置
CN106772225B (zh) * 2017-01-20 2019-03-26 大连大学 基于压缩感知的波束域doa估计
CN108957388B (zh) * 2018-05-21 2022-03-25 南京信息工程大学 一种基于协方差匹配sl0算法的mimo雷达相干信源doa估计方法

Also Published As

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

Similar Documents

Publication Publication Date Title
CN110174659B (zh) 基于迭代近端投影的mimo雷达多测量矢量doa估计方法
CN110261841B (zh) 基于迭代加权近端投影的mimo雷达单测量矢量doa估计方法
CN106980106B (zh) 阵元互耦下的稀疏doa估计方法
Wang et al. A sparse representation scheme for angle estimation in monostatic MIMO radar
CN109375154B (zh) 一种冲击噪声环境下基于均匀圆阵的相干信号参数估计方法
CN108303683B (zh) 单基地mimo雷达实值esprit非圆信号角度估计方法
CN104407335B (zh) 一种3轴交叉阵列的doa估计方法
CN110197112B (zh) 一种基于协方差修正的波束域Root-MUSIC方法
CN107870315B (zh) 一种利用迭代相位补偿技术估计任意阵列波达方向方法
CN112881972B (zh) 一种阵列模型误差下基于神经网络的波达方向估计方法
CN107991659B (zh) 基于字典学习的米波雷达低仰角目标测高方法
CN112379327A (zh) 一种基于秩损估计的二维doa估计与互耦校正方法
CN113391260B (zh) 一种基于低秩和稀疏先验的mimo雷达doa估计方法
Liu et al. Reweighted smoothed l0-norm based DOA estimation for MIMO radar
CN113567913A (zh) 基于迭代重加权可降维的二维平面doa估计方法
CN110095750B (zh) 基于准平稳信号稀疏重构的快速二维欠定测角方法
Tang et al. Gridless DOD and DOA estimation in bistatic MIMO radar using 2D-ANM and its low complexity algorithms
Li et al. DOD and DOA estimation for MIMO radar based on combined MUSIC and sparse Bayesian learning
Zhang et al. Estimation of fading coefficients in the presence of multipath propagation
CN109507634B (zh) 一种任意传感器阵列下的基于传播算子的盲远场信号波达方向估计方法
Ning et al. A velocity independent MUSIC algorithm for DOA estimation
CN110412535B (zh) 一种序贯的空时自适应处理参数估计方法
CN115421098A (zh) 嵌套面阵下降维求根music的二维doa估计方法
Xiaozhi et al. An effective DOA estimation method of coherent signals based on reconstruct weighted noise subspace
CN113093098A (zh) 基于lp范数补偿的轴向不一致矢量水听器阵列测向方法

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