CN113673317A - 基于原子范数最小化可降维的二维离格doa估计方法 - Google Patents

基于原子范数最小化可降维的二维离格doa估计方法 Download PDF

Info

Publication number
CN113673317A
CN113673317A CN202110783040.8A CN202110783040A CN113673317A CN 113673317 A CN113673317 A CN 113673317A CN 202110783040 A CN202110783040 A CN 202110783040A CN 113673317 A CN113673317 A CN 113673317A
Authority
CN
China
Prior art keywords
matrix
signal
array
snapshot
follows
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.)
Granted
Application number
CN202110783040.8A
Other languages
English (en)
Other versions
CN113673317B (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.)
University of Electronic Science and Technology of China
Original Assignee
University of Electronic Science and Technology of China
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 University of Electronic Science and Technology of China filed Critical University of Electronic Science and Technology of China
Priority to CN202110783040.8A priority Critical patent/CN113673317B/zh
Publication of CN113673317A publication Critical patent/CN113673317A/zh
Application granted granted Critical
Publication of CN113673317B publication Critical patent/CN113673317B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/02—Preprocessing
    • G06F2218/04—Denoising
    • 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
    • Y02D—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN INFORMATION AND COMMUNICATION TECHNOLOGIES [ICT], I.E. INFORMATION AND COMMUNICATION TECHNOLOGIES AIMING AT THE REDUCTION OF THEIR OWN ENERGY USE
    • Y02D30/00—Reducing energy consumption in communication networks
    • Y02D30/70—Reducing energy consumption in communication networks in wireless communication networks

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)

Abstract

本发明公开了一种基于原子范数最小化可降维的二维离格DOA估计方法，该方法利用Kronecker积的性质，对阵列接收数据模型进行变形，从而将二维联合角度估计通过降维分为两个一维DOA估计问题，目的是为了降低计算复杂度，降维后分别引用原子范数最小化ANM理论，建立多快拍下的原子范数最小化问题，并将非凸问题转为半正定规划问题，使用CVX工具箱求解，最后通过esprit算法实现方位角和俯仰角的估计。本发明可实现角度的自动配对的二维平面阵列的DOA估计，可解决网格失配问题，可解相干性，在低信噪比下也能得到很好的估计效果，以及使得计算量降低。

Description

基于原子范数最小化可降维的二维离格DOA估计方法
技术领域
本发明属于雷达通信和阵列信号处理技术领域，尤其涉及一种基于原子范数最小化可降维的二维离格DOA估计方法。
背景技术
波达方向(Direction of Arrival,DOA)估计广泛应用于雷达、声呐、无线通信等阵列信号处理领域。其中用于波达方向估计的经典方法有Schmidt.R.O等人提出的多重信号分类(Multiple Signal Classification,MUSIC)算法，Roy.R和Kailath.T提出的旋转不变子空间(Estimation of Signal Parameter by Rotational Invariant Techniques,ESPRIT)算法，但这些算法不适用于信号相干时的角度估计。2006年，Donoho正式提出了压缩感知理论，Malio utov等人在压缩感知理论上提出了L1-SVD算法。Cotter等人在MP(Matching Pursuits)算法基础上提出基于正交匹配追踪(Orthogonal matchingpursuit,OMP)算法，美国学者Tang等人在2013年首次将原子范数理论引入线谱估计当中即单快拍的DOA估计问题，并提出一种基于半正定规划的估计方法。
基于压缩感知类的算法不仅适用于独立信号的估计，而且可以解相干信号，在低信噪比下使用多快拍数据能得到很好的估计效果。相比一维空间估计，研究稀疏重构类二维平面DOA估计具有重要意义，如果直接应用一维的稀疏重构算法解决二维平面的DOA估计问题，则需要构造一个二维的过完备基矩阵，导致计算量过大，不便于工程实现。因此，赵光辉等人提出一种基于阵列流形可分离的DOA估计算法，通过对来波方向的方位角和俯仰角进行新的定义，然后对方位角和俯仰角进行解耦操作，得到两个降维后的字典矩阵，利用压缩感知框架建立稀疏信号重构优化问题，应用交替迭代思想进行求解，降维后计算量相对得到降低，但是迭代交替过程又会增加计算量，而且该算法只适用于单快拍数据，在低信噪比下估计效果不佳，(参见文献：A Sparse Representation-Based DOA EstimationAlgorithm With Separable Observation Model[J],Zhao G,Shi G,Shen F,et al,IEEEAntennas and Wireless Propagation Letters,2015,14:1586-1589)；同时赵光辉等人基于该分离模型提出一种基于稀疏表示的稳健二维波达方向估计方法，该方法在阵列流形可分离的DOA估计算法的基础上，应用一阶泰勒展开近似方法解决二维DOA估计中的基失配问题，建立稀疏信号优化问题，采用交替迭代的思想校正俯仰维和方位维的偏差，解决基失配问题，但只适用于单快拍数据模型，虽然可实现无网格DOA估计，但依旧存在低信噪比下估计效果不佳以及迭代交替过程在一定程度上会增加计算量的问题(参见文献：赵光辉等，基于稀疏表示的稳健二维波达方向估计方法:中国，201510494789.5[P].2015.12.23)。目前，许多谱峰搜索的算法都是假设入射信号的角度恰好落在划分的网格精度上，但是在实际应用中，信号的角度并不总是精确的落在已划分好的网格上，此时用谱峰搜索类的算法进行角度估计会存在一定的偏差，如果通过将网格角度划分的足够细来减少偏差，必然会带来计算量的问题，特别是对二维平面来说，势必会增加字典矩阵的长度，字典矩阵长度的增加势必会对稀疏重构算法带来巨大的计算量，增加工程负担，因此有许多学者针对无网格划分的DOA估计算法进行了研究，其中有学者提出一种二维多快拍无网格压缩波束形成方法，该方法在多快拍无网格框架下定义了接收信号的原子范数，构造一个半正定规划来求解原子范数最小化问题，该过程需要重构一个双重Toep litz矩阵，结构复杂，计算量大，最后引入矩阵束和配对方法(MaPP)来处理这个双重Toep litz矩阵并重构信源分布(参考文献：Two-dimensional multiple snapshot grid-free compressi ve beamforming[J],Y.Yang,Z.Chu,G.Ping,Mechanical Systems and Signal Processing,Vol.124,pp.524-540,2019)。为解决在尽可能降低计算量的基础上实现二维平面联合角度的无网格估计问题，有必要对已有的方法进行改进和完善。
发明内容
目前，现有的一些基于二维平面DOA估计的稀疏重构算法需要构造二维的超完备基字典矩阵，此字典矩阵使得稀疏信号重构时计算复杂度高，另外基于原子范数最小化的二维算法在重构信号时需要用CVX工具箱重构一个双重托普利兹(Toeplitz)矩阵，该矩阵结构复杂，同样存在计算量大的问题。针对以上问题，以及针对很多算法很难在多快拍接收数据模型下同时达到降低计算量、解决相干信号估计问题、解决低信噪比估计问题和解决格点失配问题的目的，提出了一种基于原子范数最小化的可降维的二维平面无网格DOA估计方法，该方法是一种二维平面阵列可降维的无网格波达方向估计方法。
本发明公开了一种基于原子范数最小化可降维的二维离格DOA估计方法，该方法无需重新对方位角和俯仰角进行定义，基于降维的思想，对阵列接收数据模型变形，将二维联合角度估计转为两个一维角度估计，从而避免双重Toeplitz矩阵的重构，以相对较低的计算量以及高分辨率实现二维平面无网格DOA估计。具体包括如下步骤：
步骤1、首先通过均匀矩形平面阵列获得信号接收数据矩阵Y＝AS+E，其中
Figure BDA0003157909420000021
为方向矩阵，
Figure BDA0003157909420000022
为多快拍下的空间信号矩阵，
Figure BDA0003157909420000023
是白噪声矩阵；
步骤2、利用克罗内克(Kronecker)积的性质，对二维均匀矩形平面阵列接收信号模型即对模型Y＝AS+E变形，分离出只含有俯仰角信息的方向矩阵
Figure BDA0003157909420000031
同时引入辅助变量X，得到对二维均匀矩形平面阵列接收信号模型变形后的模型，从而实现降维操作。将该变形后的模型分为两个一维阵列接收信号模型求解，从而将二维联合角度估计转变为两个一维角度估计。因此对二维平面阵列接收信号模型变形后的模型如下：
Figure BDA0003157909420000032
其中
Figure BDA0003157909420000033
为变形后的阵列信号接收数据矩阵，
Figure BDA0003157909420000034
是与俯仰角有关的方向矩阵，
Figure BDA0003157909420000035
是与方位角有关的方向矩阵，
Figure BDA0003157909420000036
为与信号有关的矩阵，其由S矩阵中的数据构成，辅助变量
Figure BDA0003157909420000037
是变形后的白噪声矩阵，该矩阵实际是噪声矩阵
Figure BDA0003157909420000038
的重排，(·)T表示转置，M为均匀矩形平面阵列的行数，N为均匀矩形平面阵列的列数，T为快拍数，K为入射信号个数。
步骤3、已知
Figure BDA0003157909420000039
包含T个快拍的数据信息，若T＞K，直接使用Y进行计算，会增加计算复杂度。因此，当快拍数T大于入射信号个数K，为降低计算量，可对信号接收数据矩阵Y进行奇异值分解即Y＝U′ΛVH，其中U′、V是奇异值分解的左右特征向量，Λ是奇异值分解的特征值，(·)H表示共轭转置，利用特征向量得到信号子空间的接收数据矩阵Ys＝U′ΛDK，其中DK＝[IK O]T，IK为K×K维单位矩阵，O为K×(T-K)维零矩阵，另外此时的噪声矩阵为E′s＝EV(VHV)-1DK。对Ys进行重排后变成矩阵Yss，则降低计算量后且经过变形后的阵列接收数据模型和上式类似，可表示为：
Figure BDA00031579094200000310
其中θ〞是经过降维后与方位角有关的方向矩阵，S〞是降维后与信号有关的理论矩阵，Ess是由
Figure BDA00031579094200000311
重排得到的与噪声有关的矩阵；
对Y做奇异值分解的推导过程相对较简单，对Y′的奇异值分解的表达式推导过程还没实现，相较上面的分解过程会比较复杂，且对Y′奇异值分解后，得不到理论上信号子空间的计算公式，因此不采取对Y′分解，而是选择对Y分解。Y和Y′矩阵的元素值是相同的，只是排列的方式不同。奇异值分解后的YS近似等价于Y，所以Y变形得到Y′，那么类似的YS变形为YSS。本发明中的奇异值分解都是针对原始接收信号矩阵Y，本发明相当于做了两次降维处理，对Y奇异值分解是降低列的维度(这是降低快拍数太大带来的高计算量)，Y′是降低行的维度(这是为了分别求方位角和俯仰角)。
步骤4、针对Yss＝(ψXss+Ess)的理论模型估计俯仰角，使用原子范数最小化理论，建立关于俯仰角的一维半正定规划问题，使用CVX工具箱求解，重构出接收信号矩阵
Figure BDA0003157909420000041
和具备Toeplitz性质的协方差矩阵
Figure BDA0003157909420000042
是由向量u构成的Toeplitz矩阵，使用esprit算法对
Figure BDA0003157909420000043
进行特征值分解，得到信号子空间
Figure BDA0003157909420000044
分别取Es的前M-1行向量和后M-1行向量得到矩阵
Figure BDA0003157909420000045
和
Figure BDA0003157909420000046
构造矩阵
Figure BDA0003157909420000047
并计算
Figure BDA0003157909420000048
的特征值分解
Figure BDA0003157909420000049
Λ′是矩阵
Figure BDA00031579094200000410
特征分解后的特征值，U是矩阵
Figure BDA00031579094200000411
的特征向量构成的矩阵，将U均匀分块，每一块矩阵的维度是K×K，即
Figure BDA00031579094200000412
定义矩阵
Figure BDA00031579094200000413
对ΨTLS特征值分解得到K个特征值
Figure BDA00031579094200000414
且根据公式
Figure BDA00031579094200000415
计算得到俯仰角的K个估计值
Figure BDA00031579094200000416
λ为波长，
Figure BDA00031579094200000417
为Z轴上的阵元间距；
步骤5、利用步骤4求解的接收信号矩阵
Figure BDA00031579094200000418
和俯仰角估计值
Figure BDA00031579094200000419
根据
Figure BDA00031579094200000420
计算得到包含入射信号和方位角信息的估计值
Figure BDA00031579094200000421
是由步骤4求得的
Figure BDA00031579094200000422
构成的方向矩阵；
步骤6：取出
Figure BDA00031579094200000423
矩阵的第k(k＝1,...,K)行向量，记为
Figure BDA00031579094200000424
将zk重新排列(即对zk向量按顺序抽取，每次依次取N个值作为Zk的一个列向量，第k(k＝1,...,K)次取出的N个值作为Zk的第k个列向量，连续取K次，从而得到一个矩阵
Figure BDA00031579094200000425
然后根据第k个入射信号的与方位角有关的接收信号矩阵理论模型
Figure BDA00031579094200000426
为导向矢量，Sk为与信号有关的矩阵，再次引用原子范数最小化理论，建立第k个入射信号关于方位角的一维半正定规划问题，重构出具备Toeplitz性质的协方差矩阵Tk(u′)，Tk(u′)是通过向量u′构成的Toeplitz矩阵，使用esprit算法求解方位角度估计值；首先对Tk(u′)进行特征值分解得到信号子空间
Figure BDA0003157909420000051
分别取Ek,s的前N-1行向量和后N-1行向量得到矩阵
Figure BDA0003157909420000052
和
Figure BDA0003157909420000053
构造矩阵
Figure BDA0003157909420000054
并计算
Figure BDA0003157909420000055
的特征值分解
Figure BDA0003157909420000056
Λk表示矩阵
Figure BDA0003157909420000057
分解后的特征值，Uk是矩阵
Figure BDA0003157909420000058
的特征向量构成的矩阵，将Uk均匀分块，每一块矩阵的维度是1×1，即
Figure BDA0003157909420000059
定义矩阵
Figure BDA00031579094200000510
对Ψk,TLS特征值分解得到特征值
Figure BDA00031579094200000511
且根据公式
Figure BDA00031579094200000512
计算得到与俯仰角
Figure BDA00031579094200000513
对应的方位角的估计值即
Figure BDA00031579094200000514
此时的方位角和俯仰角一一对应，可自动实现角度的配对，重复执行K次步骤6，可得K个俯仰角一一对应的方位角的估计值。
本发明公开了一种基于原子范数最小化可降维的二维离格DOA估计方法，该方法巧妙地利用了Kronecker积的性质对阵列接收数据模型进行变形，从而将二维联合角度估计转为两个一维的DOA估计，分别使用一维的原子范数最小化(ANM)算法实现二维的无网格DOA估计。该方法能够实现解相干操作，并尽可能地降低计算量，克服低信噪比下单快拍数据估计效果差的问题。
附图说明
图1为本发明的方法流程图
图2为本发明的均匀矩形平面阵列的几何结构图
图3为本发明各个算法下俯仰角均方根误差随信噪比变化的对比图
图4为本发明各个算法下方位角均方根误差随信噪比变化的对比图
图5为本发明各个算法下俯仰角和方位角均方根误差随信噪比变化的对比图
图6为本发明各个算法下俯仰角均方根误差随快拍数变化的对比图
图7为本发明各个算法下方位角均方根误差随快拍数变化的对比图
图8为本发明各个算法下俯仰角和方位角均方根误差随快拍数变化的对比图
图9为本发明各个算法下俯仰角均方根误差随阵元数变化的对比图
图10为本发明各个算法下方位角均方根误差随阵元数变化的对比图
图11为本发明各个算法下俯仰角和方位角均方根误差随阵元数变化的对比图
具体实施方式
下面结合说明书附图对本发明的具体实施方式以及工作原理作进一步详细说明。
如图2所示，X轴与YOZ平面垂直，为了更好地描述，首先进行如下定义：
方位角θ(-90°,90°)：射线在XOY面投影与X轴(法线)的夹角；
俯仰角ψ(-90°,90°)：射线与XOY面投影的夹角；
均匀矩形平面阵列：阵元均匀分布在YOZ平面上，形成均匀矩形平面阵列，以坐标原点作为参考点，Y轴上的阵元间距dy和Z轴上阵元间距dz都是半波长，并且平行于Y轴方向上均匀排列N个阵元，平行于Z轴方向上均匀排列M个阵元。
假设信源个数已知，摆放在YOZ面上的均匀矩形平面阵列的阵元个数为M×N，M为矩形平面阵列的行数，N为均匀矩形平面阵列的列数，K个窄带远场信号入射到该均匀矩形平面阵列上，快拍数为T；θ是窄带远场信号的入射方向在XOY面投影与X轴的夹角，俯仰角
Figure BDA0003157909420000061
是窄带远场信号的入射方向与XOY面投影的夹角；阵列的空域导向向量为
Figure BDA0003157909420000062
Y轴上的均匀线阵导向矢量为
Figure BDA0003157909420000063
Z轴上的均匀线阵导向矢量为
Figure BDA0003157909420000064
如图1所示，一种基于原子范数最小化可降维的二维离格DOA估计方法具体包含以下步骤：
步骤1、假设K个窄带远场信号分别从
Figure BDA0003157909420000065
方向同时入射到阵元个数为M×N的均匀矩形平面阵列上，则单快拍阵列接收信号模型定义如下：
Figure BDA0003157909420000066
上式中，
Figure BDA0003157909420000067
为单快拍下阵列接收信号向量，
Figure BDA0003157909420000068
为第k个入射信号对应的空域导向向量，
Figure BDA0003157909420000069
表示复数域，
Figure BDA00031579094200000610
为方向矩阵，
Figure BDA0003157909420000071
为第t个快拍时刻入射到矩形平面阵列的所有空间信号向量，sk(t)是第k个信号的第t个快拍数据，
Figure BDA0003157909420000072
是第t个快拍数据的白噪声向量，其中A定义如下：
Figure BDA0003157909420000073
其中
Figure BDA0003157909420000074
是Kronecker积，(·)T表示转置，
Figure BDA0003157909420000075
是第k个入射信号在Y轴上均匀线阵导向矢量，
Figure BDA0003157909420000076
是第k个入射信号在Z轴上均匀线阵导向矢量，
Figure BDA0003157909420000077
是Y轴上阵元间距，
Figure BDA0003157909420000078
是Z轴上阵元间距，λ是入射信号波长，m＝1,2,...,M，n＝1,2,...,N，其中
Figure BDA0003157909420000079
和
Figure BDA00031579094200000710
定义如下：
Figure BDA00031579094200000711
Figure BDA00031579094200000712
当t＝1,…,T，快拍数为T时，得到每个快拍的阵列接收数据矢量如下：
Figure BDA00031579094200000713
根据式(5)得到多快拍阵列接收信号模型定义如下：
Y＝AS+E (6)
上式中，
Figure BDA00031579094200000714
为均匀矩形平面阵列的多快拍接收信号矩阵，
Figure BDA00031579094200000715
为多快拍下的空间信号矩阵，
Figure BDA00031579094200000716
为方向矩阵，其定义为式(2)，
Figure BDA00031579094200000717
是多快拍下的白噪声矩阵；
步骤2、为了降低计算量，需要将二维平面DOA估计降维成两个一维的角度估计，可通过使用Kronecker积的性质，即使用性质
Figure BDA0003157909420000081
对阵列接收信号模型进行变形，vec(·)表示矩阵向量化，则首先对单快拍下的阵列接收信号模型变形，即对式(1)的第三个等式
Figure BDA0003157909420000082
变形如下：
Figure BDA0003157909420000083
其中，e′(t)表示变形后的第t个快拍数据的白噪声向量。
将式(1)的第三个等式变形为式(7)，实际是将常规的单快拍接收信号向量
Figure BDA0003157909420000084
重新排列成一个
Figure BDA0003157909420000085
矩阵，其接收信号向量重排之后如下所示：
Figure BDA0003157909420000086
其中，ym,n(t)表示矩阵yss(t)中第m行第n列的元素，m＝1,2,...,M，n＝1,2,...,N。
步骤3、以步骤2的单快拍阵列接收信号模型的变形为基础，利用Kronecker积的性质可对式(6)变形，得到变形后的多快拍下的阵列接收信号模型如下：
Y′＝ψS′θ+E′ (9)
其中，
Figure BDA0003157909420000087
Figure BDA0003157909420000088
Figure BDA0003157909420000091
其中
Figure BDA0003157909420000092
为变形后的多快拍阵列接收信号矩阵，
Figure BDA0003157909420000093
是与俯仰角有关的方向矩阵，
Figure BDA0003157909420000094
是与方位角有关的方向矩阵，
Figure BDA0003157909420000095
为变形后的与信号有关的矩阵，其由S矩阵中的数据元素构成，
Figure BDA0003157909420000096
是变形后的多快拍下的白噪声矩阵，该矩阵实际是对白噪声矩阵
Figure BDA0003157909420000097
的重排。
引入辅助变量
Figure BDA0003157909420000098
则公式(9)可以简化成如下的模型：
Figure BDA0003157909420000099
步骤4、已知
Figure BDA00031579094200000910
包含着T个快拍数据的信息，若快拍数T＞K，则计算复杂度越大。因此当快拍数T大于入射信号个数K时，为减少运算复杂度以及避免随机噪声对算法的影响，对均匀矩形平面阵列的多快拍接收信号矩阵Y进行奇异值分解，即Y＝U′ΛVH，其中U′、V是奇异值分解的左右特征向量，Λ是奇异值分解的特征值，(·)H表示共轭转置，从而采样快拍数目从T快拍数减少为K，利用信号子空间得到降维后包含入射信号的接收信号矩阵
Figure BDA00031579094200000911
且降维后的接收信号矩阵为Ys＝U′ΛDK，噪声矩阵为E′s＝EV(VHV)-1DK。其中DK＝[IK O]T，(·)T表示转置，IK为K×K维单位矩阵，O为K×(T-K)维零矩阵，那么经过奇异值分解之后的多个快拍下的阵列接收信号矩阵
Figure BDA00031579094200000912
重排如下：
Figure BDA0003157909420000101
其中，y′m,n(k)表示矩阵Yss中第m行第nk列的元素，m＝1,2,...,M，n＝1,2,...,N，k＝1,2,...,K；
Y′是快拍数为T的阵列接收信号矩阵，Yss是由奇异值分解后的信号子空间计算加重排后得到的阵列接收信号矩阵，其和矩阵Y′一样包含着入射信号的信息以及其方向信息，但Yss矩阵的维度比Y′小，因此为减少计算复杂度，可用奇异值分解后的Yss求解信号的方位角和俯仰角，而奇异值分解后的阵列接收信号理论模型与式(10)类似，可表示如下：
Figure BDA0003157909420000102
其中θ″是经过降维后与方位角有关的方向矩阵，S″是降维后与信号有关的理论矩阵，Ess是由E′s重排得到的与噪声有关的矩阵；
因此，为进一步降低计算量，避免重构双重Toeplitz矩阵，针对模型(12)，可将二维联合角度估计转为两个一维的空间角度估计，先针对Yss＝(ψXss+Ess)求解俯仰角，再针对
Figure BDA0003157909420000103
求解方位角。
步骤5、针对Yss＝(ψXss+Ess)模型估计俯仰角，此时可利用一维原子范数最小化(Atom ic Norm Minimization，ANM)算法进行估计，原子范数理论如下：
在多测量(MMV)模型下，不考虑噪声Ess时的接收信号可以写为：
Figure BDA0003157909420000104
其中
Figure BDA0003157909420000105
是
Figure BDA0003157909420000106
的第k行向量，
Figure BDA0003157909420000107
和
Figure BDA0003157909420000108
且||Φk||2＝1，其中||·||2表示求向量的2范数，第k个入射信号的俯仰角
Figure BDA0003157909420000109
的范围是(-90°,90°)，
Figure BDA00031579094200001010
是与俯仰角有关的方向矩阵，等价于ψ，
Figure BDA00031579094200001011
是包含第k个入射信号方位角和信号源信息的矩阵，等价于Xss。
空域中任意一个俯仰角
Figure BDA00031579094200001012
的范围是(-90°,90°)，定义原子集合如下：
Figure BDA0003157909420000111
上式中的
Figure BDA0003157909420000112
和Φ属于泛指，带下标k的表示接收到的第k个入射信号的参数。
Figure BDA0003157909420000113
是集合
Figure BDA00031579094200001116
的线性组合，如果rk≥0且每个入射信号的俯仰角不相同，可认为
Figure BDA0003157909420000114
的分解是K阶原子分解，通过压缩感知理论，在MMV模型的l0范数的原子方法可对空间入射信号恢复和估计出入射信号的俯仰角，则
Figure BDA0003157909420000115
的原子范数定义为：
Figure BDA0003157909420000116
其中inf表示求下界，式(15)的下界以原子系数rk的和作为优化目标，类似于l1范数，式(15)也被称为原子l1范数。
步骤6、为解决网格有限离散化所导致的网格失配问题，引入了基于连续时间信号的原子范数理论，针对式(15)的原子范数，建立原子范数最小化问题：
Figure BDA0003157909420000117
其中
Figure BDA0003157909420000118
为待求的中间变量，
Figure BDA0003157909420000119
为待求的中间变量,
Figure BDA00031579094200001110
是关于u的Toeplitz矩阵，
Figure BDA00031579094200001111
为由奇异值分解后的信号子空间计算加重排后得到的阵列接收信号矩阵，M为均匀矩形平面阵列的行数，tr(·)求矩阵的迹。
步骤7、使用凸松弛方法求解式(16)的非凸问题，将该非凸问题转化为MMV模型下的半正定规划(SDP)问题求解，如下所示：
Figure BDA00031579094200001112
Figure BDA00031579094200001113
为待求的变量，表示接收信号矩阵的估计值，τ为正则化参数，
Figure BDA00031579094200001114
表示求矩阵的2范数的平方，
Figure BDA00031579094200001115
Yss、M的含义与步骤6相同。
步骤8、使用CVX工具箱求解(17)问题，求得
Figure BDA0003157909420000121
和
Figure BDA0003157909420000122
然后使用esprit算法求解得到俯仰角的估计。
即对
Figure BDA0003157909420000123
进行特征值分解得到信号子空间
Figure BDA0003157909420000124
分别取Es的前M-1行向量和后M-1行向量得到矩阵
Figure BDA0003157909420000125
和
Figure BDA0003157909420000126
构造矩阵
Figure BDA0003157909420000127
并计算
Figure BDA0003157909420000128
的特征值分解
Figure BDA0003157909420000129
Λ′是特征值矩阵，U是特征向量构成的矩阵，将U均匀分块，每一块矩阵的维度是K×K，即
Figure BDA00031579094200001210
定义矩阵
Figure BDA00031579094200001211
对ΨTLS特征值分解得到K个特征值
Figure BDA00031579094200001212
且根据公式
Figure BDA00031579094200001213
计算得到俯仰角的估计值如下：
Figure BDA00031579094200001214
步骤9、针对模型
Figure BDA00031579094200001215
求解方位角，首先利用步骤8求解得到的信号矩阵
Figure BDA00031579094200001216
和俯仰角的估计值
Figure BDA00031579094200001217
来得到Xss的估计值
Figure BDA00031579094200001218
由
Figure BDA00031579094200001219
求解得到估计值
Figure BDA00031579094200001220
其中：
Figure BDA00031579094200001221
Figure BDA00031579094200001222
因为
Figure BDA00031579094200001223
为行满秩，
Figure BDA00031579094200001224
是列满秩，则求得
Figure BDA00031579094200001225
如下：
Figure BDA00031579094200001226
其中(·)+是广义逆矩阵。
步骤10、从
Figure BDA00031579094200001227
矩阵中取第k,k＝1,…,K行向量，记为
Figure BDA00031579094200001228
将zk重新排列(即对zk向量按顺序抽取，每次依次取N个值作为Zk的一个列向量，第k(k＝1,...,K)次取出的N个值作为Zk的第k个列向量，连续取K次，从而得到一个矩阵
Figure BDA00031579094200001229
Zk表示如下：
Figure BDA0003157909420000131
其中，zk,N(k-1)+1,zk,N(k-1)+2,...,zk,Nk表示从zk中第k次依次取出的N个值，k＝1,...,K；
步骤11、推导出第k个信号与方位角有关的接收信号矩阵理论模型为
Figure BDA0003157909420000132
Sk由S″中的元素构成，k＝1,…,K，与步骤6类似，可利用原子范数方法得到方位角的估计值，原子范数理论如下：
Figure BDA0003157909420000133
其中rk′＝||Sk||2，和
Figure BDA00031579094200001313
且||Θk||2＝1，其中||·||2表示求向量的2范数。
空域中任意一个方位角θ的范围是(-90°,90°)，定义原子集合如下：
Figure BDA0003157909420000134
因此，Zk的原子范数为:
Figure BDA0003157909420000135
步骤12、针对(25)问题，建立原子范数最小化问题，同时转为半正定规化(SDP)问题，如下所示：
针对每一个俯仰角
Figure BDA0003157909420000136
分别求出与其对应的方位角
Figure BDA0003157909420000137
则可分别建立优化问题如下：
Figure BDA0003157909420000138
其中
Figure BDA0003157909420000139
为待求的中间变量，
Figure BDA00031579094200001310
为待求的中间变量，
Figure BDA00031579094200001311
是关于u′的Toeplitz矩阵，
Figure BDA00031579094200001312
为待求的变量，在优化函数里，表示对Zk的近似估计值；N为均匀矩形平面阵列的列数，τ为正则化参数。
对于式(26)的求解，可使用CVX工具箱来求解得到Tk(u′)矩阵，再使用esprit算法求解出方位角的估计值。首先对Tk(u′)进行特征值分解得到信号子空间
Figure BDA0003157909420000141
分别取Ek,s的前N-1行向量和后N-1行向量得到矩阵
Figure BDA0003157909420000142
和
Figure BDA0003157909420000143
构造矩阵
Figure BDA0003157909420000144
并计算
Figure BDA0003157909420000145
的特征值分解
Figure BDA0003157909420000146
Uk是特征向量，Λk是特征值矩阵，将Uk均匀分块，每一块矩阵的维度是1×1，即
Figure BDA0003157909420000147
定义矩阵
Figure BDA0003157909420000148
对Ψk,TLS特征值分解得到特征值
Figure BDA0003157909420000149
且根据公式
Figure BDA00031579094200001410
可计算得到俯仰角
Figure BDA00031579094200001411
对应的方位角的估计值如下：
Figure BDA00031579094200001412
因此重复执行步骤10到步骤12共K次，可得到K个俯仰角一一对应的方位角估计值。
求解方位角时，利用每一个俯仰角对应矩阵
Figure BDA00031579094200001413
的一行重构数据，分别对应求解出方位角，实现角度的自动配对。另一种处理方式是将矩阵
Figure BDA00031579094200001414
的每一行重新排列成一个矩阵
Figure BDA00031579094200001415
将这K个矩阵
Figure BDA00031579094200001416
重新拼接成一个新的矩阵
Figure BDA00031579094200001417
然后建立原子范数最小化问题，转为半正定规划问题求解，利用esprit算法求解出K个特征值
Figure BDA00031579094200001418
然后利用信号方向向量与空间噪声子空间的正交性实现角度配对。
均匀矩形平面阵列位于坐标轴的YOZ平面，以坐标原点作为参考点，信号的入射方向的单位向量是
Figure BDA00031579094200001419
若阵元摆放在YOZ面，则每个阵元位置坐标向量为pm×n＝[0 ndy mdz]，若阵元摆放在XOY面，则每个阵元位置坐标向量为pm×n＝[mdx ndy 0]。本发明针对阵列摆放在YOZ面，若考虑阵元摆放在XOY面时，令
Figure BDA00031579094200001420
利用Kronecker积的性质分离f1,f2，建立原子范数最小化问题，使用凸松弛方法求解出f1,f2，利用信号方向向量与空间相关矩阵噪声子空间的正交性进行f1,f2的匹配，从而计算信号方位角和俯仰角的估计值。
为使本发明的目的、技术方案和技术效果更加清楚，通过仿真实验对本发明作进一步地详细描述。
本次实验针对发明基于原子范数最小化可降维的二维离格DOA估计方法进行了仿真实验，以下仿真实验中，阵列均为位于YOZ面上的均匀矩形面阵，如图2所示，入射信号均为窄带信号，均匀矩形面阵的阵元个数为16×16，Y轴和Z轴对应的阵元间距都为半波长，二维MUSIC算法、二维OMP算法、分两步估计的二维L1-SVD算法的网格搜索间隔为0.5°，蒙特卡洛实验次数为50。
用于计算各算法计算复杂度的字母定义如下：M是矩形平面阵列的行数；N是矩形平面阵列的列数，T是快拍数；K是信源数；Q是空间俯仰角从-90°到90°划分的网格数；P是空间方位角从-90°到90°划分的网格数，且Q＞＞K,P＞＞K,Q＞＞M,P＞＞N；L是蒙特卡洛实验次数，ε是期望恢复精度。
用于对比的方法有二维的MUSIC算法、二维ESPRIST算法、二维OMP算法、分两步估计的二维L1-SVD算法，最小均方误差计算公式如下：
Figure BDA0003157909420000151
仿真实验条件一：4个远场信号分别以
Figure BDA0003157909420000152
Figure BDA0003157909420000153
入射到16×16均匀矩形面阵上，假定信噪比为10dB，快拍数T＝200，正则化参数τ＝0.25，采用本发明的方法得到多次实验的角度均值估计为表1所示。
表1信噪比10dB，入射信号相互独立和信号相干的角度估计值
Figure BDA0003157909420000154
仿真实验条件二：4个远场信号分别以
Figure BDA0003157909420000155
Figure BDA0003157909420000161
入射到16×16均匀矩形面阵上，假定信噪比为-10dB，快拍数T＝200，正则化参数τ＝0.25，采用本发明的方法得到多次实验的角度均值估计为表2所示。
表2信噪比-10dB，入射信号相互独立和信号相干的角度估计值
Figure BDA0003157909420000162
仿真实验条件三：2个远场独立信号分别以
Figure BDA0003157909420000163
入射到16×16均匀矩形面阵上，信噪比为10dB，快拍数T＝200，正则化参数τ＝0.25。为验证本发明算法的性能，得到各个算法的单次运行时间比较，如表3所示，以及各个算法的计算复杂度分析如表4所示。
表3信噪比10dB时各个算法单次执行所需时间对比
Figure BDA0003157909420000164
表4各个算法的计算复杂度分析
算法 计算复杂度
2D_MUSIC O{TM<sup>2</sup>N<sup>2</sup>+M<sup>3</sup>N<sup>3</sup>+QP(MN-K)(MN+1)}
2D_ESPRIT O{TM<sup>2</sup>N<sup>2</sup>+M<sup>3</sup>N<sup>3</sup>+2K<sup>2</sup>(M-1)N+2K<sup>2</sup>(N-1)M+6K<sup>3</sup>}
2D_OMP O{MNTPQK}
2D_L1_SVD O{K<sup>3</sup>M<sup>3</sup>Q<sup>3</sup>+K<sup>3</sup>P<sup>3</sup>}
2D_ANM(未降维) O{(MN+L)<sup>3</sup>log(1/ε)+M<sup>3</sup>N<sup>3</sup>+2K<sup>2</sup>(M-1)N+2K<sup>2</sup>(N-1)M+6K<sup>3</sup>}
本发明方法 O{(M+NK)<sup>3</sup>log(1/ε)+K(N+K)<sup>3</sup>log(1/ε)+M<sup>3</sup>+2K<sup>2</sup>(M-1)+2K<sup>2</sup>(N-1)+6K<sup>3</sup>}
仿真实验条件四：2个远场独立信号分别以
Figure BDA0003157909420000171
入射到16×16均匀矩形面阵上，快拍数T＝200，正则化参数τ＝0.25，信噪比从-5:5:20变化。为验证本发明算法的性能，得到各个算法方位角、俯仰角以及两个角度差和的最小均方误差(RMSE，Root Mean Squared Error)与信噪比的关系图，如图3-图5所示。
仿真实验条件五：2个远场独立信号分别以
Figure BDA0003157909420000172
入射到16×16均匀矩形面阵上，信噪比为10dB，正则化参数τ＝0.25，快拍数从20:20:200变化。为验证本发明算法的性能，得到各个算法方位角、俯仰角以及两个角度差和的最小均方误差(RMSE，Root Mean Squared Error)与快拍数的关系图，如图6-图8所示。
仿真实验条件六：2个远场独立信号分别以
Figure BDA0003157909420000173
入射到均匀矩形面阵上，正则化参数τ＝0.25，阵元数为[4×4,6×6,8×8,10×10,12×12,14×14,16×16]。为验证本发明算法的性能，信噪比为10dB，快拍数T＝200，得到各个算法方位角、俯仰角以及两个角度差和的最小均方误差(RMSE，Ro ot Mean Squared Error)与阵元数的关系图，如图9-图11所示。
从上述这些仿真实验中可以看出，本发明在信号独立或相干时都可以得到精确的角度估计值，说明该方法具备解相干的能力，在低信噪比的环境下，也可以实现二维角度的估计。由表3可知，二维ESPRIT算法的执行时间最短，本发明所提出的方法次之，但本发明方法的角度估计误差更小；由表4可知，相比于常规未降维的二维原子范数(ANM)算法的计算复杂度，本发明方法的计算复杂度更低。由最小均方误差分别与信噪比、快拍数、阵元数的关系图可知，本发明方法的最小均方误差的值比其他四种方法都要小，说明本方法的角度估计效果优于其他方法，同时本发明方法不需要进行谱峰搜索，而且能估计出入射信号角度没有落在网格上时的场景，二维的MUSIC、OMP、L1-SVD算法的RMSE比较大，是因为这些方法不能解决网格失配带来的误差，同时这些方法中的谱峰搜索使得计算复杂度高。综上可知，本发明方法具备解相干的能力，同时低信噪比下具有很好的估计效果，也可以解决网格失配带来的误差问题，以较低的计算量实现了二维平面阵列的波达方向估计。
以上所述，仅为本发明的具体实施方式，本说明书中所公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的替代特征加以替换；所公开的所有特征、或所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以任何方式组合；本领域的技术人员根据本发明技术方案的技术特征所做出的任何非本质的添加、替换，均属于本发明的保护范围。

Claims (3)

1.一种基于原子范数最小化可降维的二维离格DOA估计方法，其特征在于，该方法具体包括如下步骤：
步骤1、假设K个窄带远场信号分别从
Figure FDA0003157909410000011
方向同时入射到阵元个数为M×N的均匀矩形平面阵列上，阵元均匀分布在YOZ平面上，X轴与YOZ平面垂直，则单快拍阵列接收信号模型定义如下：
Figure FDA0003157909410000012
上式中，方位角θk是第k个窄带远场信号的入射方向在XOY面投影与X轴的夹角，俯仰角
Figure FDA0003157909410000013
是第k个窄带远场信号的入射方向与XOY面投影的夹角；
Figure FDA0003157909410000014
为单快拍下阵列接收信号向量，
Figure FDA0003157909410000015
为第k个入射信号对应的空域导向向量，
Figure FDA0003157909410000016
表示复数域，
Figure FDA0003157909410000017
为方向矩阵，
Figure FDA0003157909410000018
为第t个快拍时刻入射到均匀矩形平面阵列的所有空间信号向量，sk(t)是第k个入射信号的第t个快拍数据，
Figure FDA0003157909410000019
是第t个快拍数据的白噪声向量，其中A定义如下：
Figure FDA00031579094100000110
其中
Figure FDA00031579094100000111
是Kronecker积，(·)T表示转置，
Figure FDA00031579094100000112
是第k个入射信号在Y轴上均匀线阵导向矢量，
Figure FDA00031579094100000113
是第k个入射信号在Z轴上均匀线阵导向矢量，
Figure FDA00031579094100000114
是Y轴上阵元间距，
Figure FDA00031579094100000115
是Z轴上阵元间距，λ是入射信号波长，m＝1,2,...,M，n＝1,2,...,N，其中
Figure FDA00031579094100000116
和
Figure FDA00031579094100000117
定义如下：
Figure FDA00031579094100000118
Figure FDA0003157909410000021
当t＝1,…,T，快拍数为T时，得到每个快拍的阵列接收数据矢量如下：
Figure FDA00031579094100000211
根据式(5)得到多快拍阵列接收信号模型定义如下：
Y＝AS+E (6)
上式中，
Figure FDA0003157909410000022
为均匀矩形平面阵列的多快拍接收信号矩阵，
Figure FDA0003157909410000023
为多快拍下的空间信号矩阵，
Figure FDA0003157909410000024
为方向矩阵，其定义为式(2)，
Figure FDA0003157909410000025
是多快拍下的白噪声矩阵；
步骤2、通过使用Kronecker积的性质对单快拍阵列接收信号模型变形，即对式(1)的第三个等式
Figure FDA0003157909410000026
变形如下：
Figure FDA0003157909410000027
其中，e′(t)表示变形后的第t个快拍数据的白噪声向量；
将式(1)的第三个等式变形为式(7)，实际是将常规的单快拍下阵列接收信号向量
Figure FDA0003157909410000028
重新排列成一个
Figure FDA0003157909410000029
矩阵，其接收信号向量重排之后如下所示：
Figure FDA00031579094100000210
其中，ym,n(t)表示矩阵yss(t)中第m行第n列的元素，m＝1,2,...,M，n＝1,2,...,N；
步骤3、以步骤2的单快拍阵列接收信号模型的变形为基础，利用Kronecker积的性质对式(6)变形，得到变形后的多快拍阵列接收信号模型如下：
Y′＝ψS′θ+E′ (9)
其中，
Figure FDA0003157909410000031
Figure FDA0003157909410000032
Figure FDA0003157909410000033
其中
Figure FDA0003157909410000034
为变形后的多快拍阵列接收信号矩阵，
Figure FDA0003157909410000035
是与俯仰角有关的方向矩阵，
Figure FDA0003157909410000036
是与方位角有关的方向矩阵，
Figure FDA0003157909410000037
为变形后的与信号有关的矩阵，其由S矩阵中的数据元素构成，
Figure FDA0003157909410000038
是变形后的多快拍下的白噪声矩阵，该矩阵实际是对白噪声矩阵
Figure FDA0003157909410000039
的重排；
引入辅助变量
Figure FDA00031579094100000310
则公式(9)简化成如下的模型：
Figure FDA00031579094100000311
步骤4、已知
Figure FDA00031579094100000312
包含着T个快拍数据的信息，当快拍数T大于入射信号个数K时，为减少运算复杂度以及避免随机噪声对算法的影响，对均匀矩形平面阵列的多快拍接收信号矩阵Y进行奇异值分解，即Y＝U′ΛVH，其中U′、V是奇异值分解的左右特征向量，Λ是奇异值分解的特征值，(·)H表示共轭转置，从而采样快拍数目从T快拍数减少为K，利用信号子空间得到降维后包含入射信号的接收信号矩阵
Figure FDA00031579094100000313
且降维后的接收信号矩阵为Ys＝U′ΛDK，降维后的白噪声矩阵为E′s＝EV(VHV)-1DK，其中DK＝[IK O]T，(·)T表示转置，IK为K×K维单位矩阵，O为K×(T-K)维零矩阵，那么经过奇异值分解之后的多快拍接收信号矩阵
Figure FDA0003157909410000041
重排如下：
Figure FDA0003157909410000042
其中，y′m,n(k)表示矩阵Yss中第m行第nk列的元素，m＝1,2,...,M，n＝1,2,...,N，k＝1,2,...,K；
Y′是快拍数为T的阵列接收信号矩阵，Yss是由奇异值分解后的信号子空间计算加重排后得到的阵列接收信号矩阵，其和矩阵Y′一样包含着入射信号的信息以及其方向信息，但Yss矩阵的维度比Y′小，因此为减少计算复杂度，采用奇异值分解后的Yss求解信号的方位角和俯仰角，而奇异值分解后的阵列接收信号理论模型表示如下：
Figure FDA0003157909410000043
其中θ″是经过降维后与方位角有关的方向矩阵，S″是降维后与信号有关的理论矩阵，Ess是由E′s重排得到的与噪声有关的矩阵；
为进一步降低计算量，避免重构双重Toeplitz矩阵，针对模型(12)，将二维联合角度估计转为两个一维的空间角度估计，先针对Yss＝(ψXss+Ess)求解俯仰角，再针对
Figure FDA0003157909410000044
求解方位角；
步骤5、针对Yss＝(ψXss+Ess)模型估计俯仰角，此时利用一维原子范数最小化ANM算法进行估计，原子范数理论如下：
在多测量MMV模型下，不考虑噪声Ess时的接收信号表示为：
Figure FDA0003157909410000045
其中
Figure FDA0003157909410000046
是
Figure FDA0003157909410000047
的第k行向量，
Figure FDA0003157909410000048
且||Φk||2＝1，其中||·||2表示求向量的2范数，第k个入射信号的俯仰角
Figure FDA0003157909410000051
的范围是(-90°,90°)，
Figure FDA0003157909410000052
是与俯仰角有关的方向矩阵，等价于ψ，
Figure FDA0003157909410000053
是包含第k个入射信号方位角和信号源信息的矩阵，等价于Xss；
空域中任意一个俯仰角
Figure FDA0003157909410000054
的范围是(-90°,90°)，定义原子集合如下：
Figure FDA0003157909410000055
上式中的
Figure FDA0003157909410000056
和Φ属于泛指；
Figure FDA0003157909410000057
是集合
Figure FDA0003157909410000058
的线性组合，如果rk≥0且每个入射信号的俯仰角不相同，则认为
Figure FDA0003157909410000059
的分解是K阶原子分解，通过压缩感知理论，在MMV模型的l0范数的原子方法对空间入射信号恢复和估计出入射信号的俯仰角，则
Figure FDA00031579094100000510
的原子范数定义为：
Figure FDA00031579094100000511
其中inf表示求下界，式(15)的下界以原子系数rk的和作为优化目标，式(15)也被称为原子l1范数；
步骤6、为解决网格有限离散化所导致的网格失配问题，引入基于连续时间信号的原子范数理论，针对式(15)的原子范数，建立原子范数最小化问题：
Figure FDA00031579094100000512
其中
Figure FDA00031579094100000513
为待求的中间变量，
Figure FDA00031579094100000514
为待求的中间变量,
Figure FDA00031579094100000515
是关于u的Toeplitz矩阵，
Figure FDA00031579094100000516
为由奇异值分解后的信号子空间计算加重排后得到的阵列接收信号矩阵，M为均匀矩形平面阵列的行数，tr(·)求矩阵的迹；
步骤7、使用凸松弛方法求解式(16)的非凸问题，将该非凸问题转化为MMV模型下的半正定规划SDP问题求解，如下所示：
Figure FDA0003157909410000061
Figure FDA0003157909410000062
为待求的变量，表示接收信号矩阵的估计值，τ为正则化参数，
Figure FDA0003157909410000063
表示求矩阵的2范数的平方；
步骤8、使用CVX工具箱求解(17)问题，求得
Figure FDA0003157909410000064
和
Figure FDA0003157909410000065
然后使用esprit算法求解得到俯仰角的估计，即对
Figure FDA0003157909410000066
进行特征值分解得到信号子空间
Figure FDA0003157909410000067
分别取Es的前M-1行向量和后M-1行向量得到矩阵
Figure FDA0003157909410000068
和
Figure FDA0003157909410000069
构造矩阵
Figure FDA00031579094100000610
并计算
Figure FDA00031579094100000611
的特征值分解
Figure FDA00031579094100000612
Λ′是特征值矩阵，U是特征向量构成的矩阵，将U均匀分块，每一块矩阵的维度是K×K，即
Figure FDA00031579094100000613
定义矩阵
Figure FDA00031579094100000614
对ΨTLS特征值分解得到K个特征值e-jφk,k＝1,…,K，且根据公式
Figure FDA00031579094100000615
计算得到俯仰角的估计值如下：
Figure FDA00031579094100000616
步骤9、针对模型
Figure FDA00031579094100000617
求解方位角，首先利用步骤8求解得到的信号矩阵
Figure FDA00031579094100000618
和俯仰角的估计值
Figure FDA00031579094100000619
来得到Xss的估计值
Figure FDA00031579094100000620
由
Figure FDA00031579094100000621
求解得到估计值
Figure FDA00031579094100000622
其中：
Figure FDA00031579094100000623
Figure FDA00031579094100000624
因为
Figure FDA00031579094100000625
为行满秩，
Figure FDA00031579094100000626
是列满秩，则求得
Figure FDA00031579094100000627
如下：
Figure FDA00031579094100000628
其中(·)+是广义逆矩阵；
步骤10、从
Figure FDA00031579094100000629
矩阵中取第k,k＝1,…,K行向量，记为
Figure FDA00031579094100000630
将zk重新排列，即对zk向量按顺序抽取，每次依次取N个值作为Zk的一个列向量，第k次取出的N个值作为Zk的第k个列向量，连续取K次，从而得到一个矩阵
Figure FDA0003157909410000071
Zk表示如下：
Figure FDA0003157909410000072
其中，zk,N(k-1)+1,zk,N(k-1)+2,...,zk,Nk表示从zk中第k次依次取出的N个值，k＝1,...,K；
步骤11、推导出第k个信号与方位角有关的接收信号矩阵理论模型为
Figure FDA0003157909410000073
Sk由S〞中的元素构成，k＝1,…,K，利用原子范数方法得到方位角的估计值，原子范数理论如下：
Figure FDA0003157909410000074
其中r′k＝||Sk||2，Θk＝r′k -1Sk，且||Θk||2＝1，其中||·||2表示求向量的2范数；
空域中任意一个方位角θ的范围是(-90°,90°)，定义原子集合如下：
Figure FDA0003157909410000075
上式中的θ和Θ属于泛指；
因此，Zk的原子范数为:
Figure FDA0003157909410000076
步骤12、针对式(25)，建立原子范数最小化问题，同时转为半正定规化SDP问题，如下所示：
针对每一个俯仰角
Figure FDA0003157909410000077
分别求出与其对应的方位角
Figure FDA0003157909410000078
分别建立优化问题如下：
Figure FDA0003157909410000079
其中
Figure FDA00031579094100000710
为待求的中间变量，
Figure FDA00031579094100000711
为待求的中间变量，
Figure FDA00031579094100000712
是关于u′的Toeplitz矩阵，
Figure FDA0003157909410000081
为待求的变量，在优化函数里，表示对Zk的近似估计值；N为均匀矩形平面阵列的列数，τ为正则化参数；
对于式(26)的求解，使用CVX工具箱求解得到Tk(u′)矩阵，再使用esprit算法求解出方位角的估计值，即首先对Tk(u′)进行特征值分解得到信号子空间
Figure FDA0003157909410000082
分别取Ek,s的前N-1行向量和后N-1行向量得到矩阵
Figure FDA0003157909410000083
和
Figure FDA0003157909410000084
构造矩阵
Figure FDA0003157909410000085
并计算
Figure FDA0003157909410000086
的特征值分解
Figure FDA0003157909410000087
Uk是特征向量，Λk是特征值矩阵，将Uk均匀分块，每一块矩阵的维度是1×1，即
Figure FDA0003157909410000088
定义矩阵
Figure FDA0003157909410000089
对Ψk,TLS特征值分解得到特征值
Figure FDA00031579094100000813
且根据公式
Figure FDA00031579094100000810
计算得到俯仰角
Figure FDA00031579094100000811
对应的方位角的估计值如下：
Figure FDA00031579094100000812
重复执行步骤10到步骤12共K次，得到与K个俯仰角一一对应的方位角估计值。
2.根据权利要求1所述的基于原子范数最小化可降维的二维离格DOA估计方法，其特征在于，所述均匀矩形平面阵列中M＝16，N＝16。
3.根据权利要求2所述的基于原子范数最小化可降维的二维离格DOA估计方法，其特征在于，所述正则化参数τ＝0.25。
CN202110783040.8A 2021-07-12 2021-07-12 基于原子范数最小化可降维的二维离格doa估计方法 Active CN113673317B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110783040.8A CN113673317B (zh) 2021-07-12 2021-07-12 基于原子范数最小化可降维的二维离格doa估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110783040.8A CN113673317B (zh) 2021-07-12 2021-07-12 基于原子范数最小化可降维的二维离格doa估计方法

Publications (2)

Publication Number Publication Date
CN113673317A true CN113673317A (zh) 2021-11-19
CN113673317B CN113673317B (zh) 2023-04-07

Family

ID=78538870

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110783040.8A Active CN113673317B (zh) 2021-07-12 2021-07-12 基于原子范数最小化可降维的二维离格doa估计方法

Country Status (1)

Country Link
CN (1) CN113673317B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114280544A (zh) * 2021-12-02 2022-04-05 电子科技大学 一种基于松弛优化的最小过渡带宽度方向图赋形方法
CN116299150A (zh) * 2022-12-27 2023-06-23 南京航空航天大学 一种均匀面阵中降维传播算子的二维doa估计方法
CN116430347A (zh) * 2023-06-13 2023-07-14 成都实时技术股份有限公司 一种雷达数据采集与存储方法

Citations (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20160057556A1 (en) * 2013-03-22 2016-02-25 Thomson Licensing Method and apparatus for enhancing directivity of a 1st order ambisonics signal
CN106569172A (zh) * 2016-10-13 2017-04-19 北京邮电大学 二维doa估计方法
WO2017161874A1 (zh) * 2016-03-23 2017-09-28 中兴通讯股份有限公司 一种mimo雷达波达方向估计方法和装置
CN107703477A (zh) * 2017-09-11 2018-02-16 电子科技大学 基于块稀疏贝叶斯学习的准平稳宽带阵列信号波达方向估计方法
WO2018094565A1 (zh) * 2016-11-22 2018-05-31 深圳大学 脉冲噪声下的波束成形方法及装置
CN109917328A (zh) * 2019-04-15 2019-06-21 南京邮电大学 一种基于原子范数最小化的l型阵列波达方向估计方法
CN111337893A (zh) * 2019-12-19 2020-06-26 江苏大学 一种基于实值稀疏贝叶斯学习的离格doa估计方法
CN111398890A (zh) * 2020-03-24 2020-07-10 南京信息工程大学 一种基于协方差矩阵重构的布谷鸟搜索doa估计方法
CN111812581A (zh) * 2020-06-16 2020-10-23 重庆大学 基于原子范数的球面阵列声源波达方向估计方法
CN111983553A (zh) * 2020-08-20 2020-11-24 上海无线电设备研究所 一种基于互质多载频稀疏阵列的无网格doa估计方法
CN112162280A (zh) * 2020-08-25 2021-01-01 中国人民解放军空军预警学院雷达士官学校 基于原子范数最小化的sf isar一维高分辨距离成像方法
CN112305495A (zh) * 2020-10-22 2021-02-02 南昌工程学院 一种基于原子范数最小的互质阵列协方差矩阵重构方法
CN112363110A (zh) * 2020-11-30 2021-02-12 海南大学 一种基于嵌套交叉偶极子阵列的无网格单比特doa估计方法
WO2021068496A1 (zh) * 2020-05-03 2021-04-15 浙江大学 基于结构化虚拟域张量信号处理的互质面阵二维波达方向估计方法
CN112800599A (zh) * 2021-01-15 2021-05-14 吉林大学 一种阵元失配情况下基于admm的无网格doa估计方法

Patent Citations (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20160057556A1 (en) * 2013-03-22 2016-02-25 Thomson Licensing Method and apparatus for enhancing directivity of a 1st order ambisonics signal
WO2017161874A1 (zh) * 2016-03-23 2017-09-28 中兴通讯股份有限公司 一种mimo雷达波达方向估计方法和装置
CN106569172A (zh) * 2016-10-13 2017-04-19 北京邮电大学 二维doa估计方法
WO2018094565A1 (zh) * 2016-11-22 2018-05-31 深圳大学 脉冲噪声下的波束成形方法及装置
CN107703477A (zh) * 2017-09-11 2018-02-16 电子科技大学 基于块稀疏贝叶斯学习的准平稳宽带阵列信号波达方向估计方法
CN109917328A (zh) * 2019-04-15 2019-06-21 南京邮电大学 一种基于原子范数最小化的l型阵列波达方向估计方法
CN111337893A (zh) * 2019-12-19 2020-06-26 江苏大学 一种基于实值稀疏贝叶斯学习的离格doa估计方法
CN111398890A (zh) * 2020-03-24 2020-07-10 南京信息工程大学 一种基于协方差矩阵重构的布谷鸟搜索doa估计方法
WO2021068496A1 (zh) * 2020-05-03 2021-04-15 浙江大学 基于结构化虚拟域张量信号处理的互质面阵二维波达方向估计方法
CN111812581A (zh) * 2020-06-16 2020-10-23 重庆大学 基于原子范数的球面阵列声源波达方向估计方法
CN111983553A (zh) * 2020-08-20 2020-11-24 上海无线电设备研究所 一种基于互质多载频稀疏阵列的无网格doa估计方法
CN112162280A (zh) * 2020-08-25 2021-01-01 中国人民解放军空军预警学院雷达士官学校 基于原子范数最小化的sf isar一维高分辨距离成像方法
CN112305495A (zh) * 2020-10-22 2021-02-02 南昌工程学院 一种基于原子范数最小的互质阵列协方差矩阵重构方法
CN112363110A (zh) * 2020-11-30 2021-02-12 海南大学 一种基于嵌套交叉偶极子阵列的无网格单比特doa估计方法
CN112800599A (zh) * 2021-01-15 2021-05-14 吉林大学 一种阵元失配情况下基于admm的无网格doa估计方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
TANG WENGEN等: ""Grid-Free DOA and DOA Estimation for MIMO radar via duality-based 2D atomic norm minimization"" *
WEN-GEN TANG等: ""Gridless DOD and DOA estimation in bistatic MIMO radar using 2D-ANM and its low complexity algorithms"" *
ZIYAN CHENG等: ""Jiont Design for MIMO Radar and Down Link Communication Systems Coexisece"" *
卢爱红等: ""基于原子范数最小化的二维稀疏阵列波达角估计算法"" *
徐峰等: ""基于快速张量分解的波束空间MIMO雷达二维DOA估计算法"" *
陈海华等: ""基于宽带均匀同心球阵列的低复杂度二维波达方向估计算法"" *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114280544A (zh) * 2021-12-02 2022-04-05 电子科技大学 一种基于松弛优化的最小过渡带宽度方向图赋形方法
CN116299150A (zh) * 2022-12-27 2023-06-23 南京航空航天大学 一种均匀面阵中降维传播算子的二维doa估计方法
CN116299150B (zh) * 2022-12-27 2023-12-01 南京航空航天大学 一种均匀面阵中降维传播算子的二维doa估计方法
CN116430347A (zh) * 2023-06-13 2023-07-14 成都实时技术股份有限公司 一种雷达数据采集与存储方法
CN116430347B (zh) * 2023-06-13 2023-08-22 成都实时技术股份有限公司 一种雷达数据采集与存储方法

Also Published As

Publication number Publication date
CN113673317B (zh) 2023-04-07

Similar Documents

Publication Publication Date Title
CN113673317B (zh) 基于原子范数最小化可降维的二维离格doa估计方法
CN109444810B (zh) 一种非负稀疏贝叶斯学习框架下的互质阵列非网格doa估计方法
CN107037392B (zh) 一种基于压缩感知的自由度增加型互质阵列波达方向估计方法
CN109655799B (zh) 基于iaa的协方差矩阵向量化的非均匀稀疏阵列测向方法
Yan et al. Real-valued root-MUSIC for DOA estimation with reduced-dimension EVD/SVD computation
CN110082708A (zh) 非均匀阵列设计和波达方向估计方法
US11422177B2 (en) Spatial spectrum estimation method with enhanced degree-of-freedom based on block sampling tensor construction for coprime planar array
CN111610486A (zh) 基于平面互质阵列虚拟域张量空间谱搜索的高分辨精确二维波达方向估计方法
CN112904272B (zh) 基于互相关张量的三维互质立方阵列波达方向估计方法
CN109143151B (zh) 部分阵元损坏的均匀面阵张量重构方法及信源定位方法
CN111610485A (zh) 基于平面互质阵列块采样张量信号构造的自由度增强型空间谱估计方法
CN112269172A (zh) 一种基于张量结构的嵌套mimo雷达角度估计方法和装置
CN111965591A (zh) 一种基于四阶累积量矢量化dft的测向估计方法
CN112731275A (zh) 一种基于零化插值的互质阵部分极化信号参数估计方法
CN111983554A (zh) 非均匀l阵下的高精度二维doa估计
CN113567913B (zh) 基于迭代重加权可降维的二维平面doa估计方法
Yan et al. Two‐Step Root‐MUSIC for Direction of Arrival Estimation without EVD/SVD Computation
US11300648B2 (en) High-resolution, accurate, two-dimensional direction-of-arrival estimation method based on coarray tensor spatial spectrum searching with co-prime planar array
Liu et al. Real-valued sparse Bayesian learning algorithm for off-grid DOA estimation in the beamspace
CN114648041A (zh) 一种基于平行稀疏阵列的二维欠定doa估计算法
CN113281698A (zh) 一种嵌套阵中基于级联的非高斯信源测向方法
CN112327244A (zh) 一种基于l型阵列的二维非相干分布式目标参数估计方法
Li et al. DOA Estimation Based on Sparse Reconstruction via Acoustic Vector Sensor Array under Non-uniform Noise
Zhan et al. Increasing utilization of redundant virtual array for DOA estimation based on coprime array
Zhang et al. Gridless sparse methods based on fourth-order cumulant for DOA estimation

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