CN110046326B - 一种时频doa估计方法 - Google Patents

一种时频doa估计方法 Download PDF

Info

Publication number
CN110046326B
CN110046326B CN201910349676.4A CN201910349676A CN110046326B CN 110046326 B CN110046326 B CN 110046326B CN 201910349676 A CN201910349676 A CN 201910349676A CN 110046326 B CN110046326 B CN 110046326B
Authority
CN
China
Prior art keywords
quantum
time
landmine
frequency
signal
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
CN201910349676.4A
Other languages
English (en)
Other versions
CN110046326A (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.)
Jiaxing Nuoaidi Communication Technology Co ltd
Original Assignee
Harbin Engineering 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 Harbin Engineering University filed Critical Harbin Engineering University
Priority to CN201910349676.4A priority Critical patent/CN110046326B/zh
Publication of CN110046326A publication Critical patent/CN110046326A/zh
Application granted granted Critical
Publication of CN110046326B publication Critical patent/CN110046326B/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
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Theoretical Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Evolutionary Biology (AREA)
  • Computing Systems (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)

Abstract

本发明公开一种时频DOA估计方法，包括：建立阵列接收的时域数据模型；对时域数据进行快拍采样；对快拍采样数据进行时频分析得到PWVD矩阵；计算时频平均的快拍采样数据PWVD矩阵；构造极大似然方程；初始化量子地雷量子位置；由极大似然方程构造适应度函数；模拟量子地雷爆炸过程获得量子弹片的量子位置；计算量子弹片量子位置映射态的适应度函数值，选择适应度大的优秀量子位置作为放置量子地雷的量子位置，用于引爆下一代的量子地雷，根据所有量子位置的适应度更新全局最优量子位置；达到最大迭代次数后，输出信号方位角最优估计值，本发明能在较短时间内得到较准确的非平稳信号DOA估计结果，并且在信号源为相干源的条件下仍有效。

Description

一种时频DOA估计方法
技术领域
本发明涉及一种时频DOA估计方法，特别是一种基于量子地雷爆炸机理的时频DOA估计方法，属于阵列信号处理领域。
背景技术
波达方向(DOA)估计是阵列信号处理中的一个重要研究方向，在雷达、声呐、无线电天文学等领域中有广泛的应用。当前，绝大多数DoA估计方法都是针对平稳信号进行估计的，由于传统的平稳信号DOA估计方法局限于空域处理，对于非平稳信号的DOA估计则性能恶化严重，甚至失效，所以需设计专门针对非平稳信号的DOA估计方法。
对于非平稳信号的分析，分数阶傅里叶变换，循环谱估计，小波变换，时频分析等都是行之有效的方法，其中，时频分析方法的基本思想是通过设计时间和频率的联合函数，同时描述信号在时频域上不同时间及频率处的能量密度，这种方法实现简单，且具有优良的时频聚集性，将时频分析方法与阵列信号处理结合的进行时频DOA估计，是一种解决非平稳信号DoA估计问题的有效方法。
如何快速高精度获得时频DOA估计结果是一个经典的工程难题，在众多时频DOA估计方法中，使用极大似然原理的时频DOA估计是一种性能优良的方法，且可以对相干源进行有效估计，因此，尽管时频极大似然DOA估计的实现过程较为复杂，但是这种方法仍有广阔的应用前景。
为了降低时频极大似然DOA估计的复杂度，需要进行新方法的设计进而有效求解。当前，许多智能计算方法被应用于求解复杂的优化问题，其中地雷爆炸搜索是一种启发式智能计算方法，受现实世界中地雷爆炸产生弹片，引爆其他潜在地雷的现象启发，具有新颖、有效的特点，但易陷入局部收敛，因此需要改进地雷爆炸方法使其具有收敛速度快、搜索精度高的优点。
通过对现有技术文献的检索发现，Yimin Zhang等在《Journal of the FranklinInstitute》(337(2000)483-497)上发表的“Time-frequency maximum likelihoodmethods for direction finding”中提出了一种时频极大似然DOA估计方法，可以解相干且不需要进行使阵列孔径损失的预处理，但没有解决极大似然估计实现过程复杂的难题。Siyu Chen等在《IEEE International Symposium on Microwave》(10.1109/MAPE.2015.7510425)上发表的“Direction of Arrival Estimation of NonstationarySignal Based on Time-frequency Music”中提出了一种时频MUSIC估计方法，此方法若使用较小的计算量，则需大的扫描间隔，不能保证对非平稳信号的精确估计，且不能对相干源进行有效估计。已有文献检索结果表明，对于非平稳信号的DOA估计问题，缺少一种能在较短时间内得到较准确的估计结果，并且在信号源为相干源的条件下仍有效的时频DOA估计方法。
发明内容
针对上述现有技术，本发明要解决的技术问题是提供一种能够对非平稳信号进行快速、稳健的DOA估计且可以对多个相干的非平稳信号进行DOA估计的时频DOA估计方法。
为解决上述技术问题，本发明提供一种时频DOA估计方法，包括以下步骤：
步骤一：建立阵列接收的时域数据模型：
在信号为远场窄带信号、阵列形式为均匀线阵且阵元各向同性的假设下，在t时刻，若N个方位角为θ＝[θ1,θ2,…,θN]的调频信号入射到含有M个阵元的阵列上，其中第n个调频信号表示为sn(t)，n＝1,2,…,N，则阵列上第m个阵元接收到的数据ym(t)满足：
Figure BDA0002043512500000021
m＝1,2,…,M，其中τmn为第n个信号入射到第m个阵元时相对于入射到参考阵元的空间延迟，vm(t)为第m个阵元处在t时刻与信号独立的高斯白噪声，将M个阵元在t时刻接收到的数据排成列矢量
Figure BDA0002043512500000022
步骤二：获取快拍采样数据：
对阵列中M个阵元接收的数据进行K次快拍采样，得到快拍采样数据矩阵：
Figure BDA0002043512500000023
ym(k)为第m个阵元在第k个快拍采样时刻的数据，m＝1,2,…,M，k＝1,2,…K，Y也可表示为Y＝A(θ)S+V，其中
Figure BDA0002043512500000024
为阵列流型矩阵，ωn，n＝1,2,…,N，为快拍采样时间点处第n个调频信号的瞬时角频率，S是大小为N×K的信号矩阵，V是大小为M×K的噪声矩阵；
步骤三：对快拍采样数据进行时频分析：
使用离散时间形式的Cohen类时频分布中的PWVD对快拍采样数据进行时频分析，在信号的时频域平面上，第n个信号的第k个快拍采样时间点tk及此时信号的瞬时频率fn,k确定时频点(tk,fn,k)，n＝1,2,…,N，k＝1,2,…,K，则任意两阵元的快拍采样数据使用长度为L的矩形窗的PWVD在时频点(tk,fn,k)处的值满足：
Figure BDA0002043512500000031
σ为滞后变量，w＝1,2,…,M，u＝1,2,…,M，(·)H代表共轭转置运算，其中矩形窗长度L为奇数且应足够小以保证在观察区间内阵列流型矩阵是近似非时变的；将各个阵元快拍采样数据的自PWVD项及互PWVD项排成矩阵，得到时频点(tk,fn,k)处快拍采样数据的PWVD矩阵：
Figure BDA0002043512500000032
步骤四：计算时频平均的快拍采样数据PWVD矩阵：
略去数据段两端的PWVD数据，在每个信号的全部K个时频点中，保留中间的K-L+1个时频点，对这些时频点处的PWVD数值进行算术平均，得到时频平均后的快拍采样数据PWVD矩阵
Figure BDA0002043512500000033
步骤五：构造时频极大似然DOA估计的极大似然方程：
极大似然方程为
Figure BDA0002043512500000034
其中
Figure BDA0002043512500000035
为信号方位角的最优估计值，A(θ)为解空间中一个可能的解对应的导向矢量矩阵，tr代表矩阵求迹运算；
步骤六：量子地雷量子位置的初始化：
设定量子地雷的总数量P，量子地雷量子位置的演化寻优过程分为勘探过程和探测过程，分别设定勘探过程的最大迭代次数G1，探测过程的最大迭代次数G2，当前迭代数为g∈[1,G1+G2]，第p个量子地雷的N维初代量子位置为
Figure BDA0002043512500000036
p＝1,2,…,P，初代量子地雷量子位置每一维取值为[0,1]区间内的均匀随机数；
步骤七：由极大似然方程构造适应度函数，计算每个量子地雷量子位置适应度，根据适应度值获得全局最优量子位置：
在第g次迭代中，将量子地雷量子位置
Figure BDA0002043512500000041
映射到信号方位角的解空间
Figure BDA0002043512500000042
范围内，得到量子地雷量子位置的映射态
Figure BDA0002043512500000043
p＝1,2,…,P，构造其适应度函数为
Figure BDA0002043512500000044
p＝1,2,…,P，其中
Figure BDA0002043512500000045
为第g次迭代中第p个量子地雷量子位置的映射态对应的阵列流型矩阵，计算所有量子地雷量子位置适应度，适应度值最大的量子地雷量子位置为全局最优量子位置，记作
Figure BDA0002043512500000046
b∈[1,P]；
步骤八：量子地雷爆炸产生量子弹片，使用模拟量子旋转门模拟量子地雷爆炸过程获得量子弹片的量子位置：
每个量子地雷爆炸产生的量子弹片个数为Q，使用模拟量子旋转门计算第p个量子地雷爆炸产生的第q个量子弹片的第n维量子位置为
Figure BDA0002043512500000047
p＝1,2,…,P，q＝1,2,…,Q，n＝1,2,…,N，在勘探过程中g∈[1,G1]，第p个量子地雷的第n维量子位置对应的量子旋转角为
Figure BDA0002043512500000048
其中
Figure BDA0002043512500000049
为[0,1]区间内的均匀随机数；在探测过程中g∈[G1+1,G1+G2]，第p个量子地雷的第n维量子位置对应的量子旋转角为
Figure BDA00020435125000000410
p≠b；另外规定当p＝b时，
Figure BDA00020435125000000411
其值取正或负的概率相等，其中
Figure BDA00020435125000000412
为[0,1]区间内的均匀随机数；
步骤九：对于每个量子地雷爆炸产生的量子弹片量子位置，计算量子弹片量子位置映射态的适应度函数值，应用贪婪选择机制，选择适应度大的优秀量子位置作为放置量子地雷的量子位置，并在此量子位置上引爆下一代的量子地雷，并根据所有量子位置的适应度更新全局最优量子位置；
计算量子弹片量子位置的映射态
Figure BDA00020435125000000413
p＝1,2,…,P，q＝1,2,…,Q，量子弹片量子位置的映射态和量子地雷量子位置的映射态的适应度函数定义相同，计算第p个量子地雷爆炸产生的第q个量子弹片量子位置映射态的适应度函数值
Figure BDA0002043512500000051
p＝1,2,…,P，q＝1,2,…,Q，应用贪婪选择，选择最优量子位置用于引爆下一代量子地雷，逐个判断量子弹片量子位置的适应度值是否大于当前量子地雷量子位置的适应度值
Figure BDA0002043512500000052
若
Figure BDA0002043512500000053
则有
Figure BDA0002043512500000054
q＝1,2,…,Q；否则
Figure BDA0002043512500000055
p＝1,2,…,P；将全局最优量子位置更新为适应度最大的量子地雷量子位置
Figure BDA0002043512500000056
步骤十：判断是否达到最大迭代次数，若g＜G1+G2则令g＝g+1，返回步骤八；若g＝G1+G2，则将最后一代中最优的量子地雷量子位置的映射态作为信号方位角的最优估计值输出。
本发明有益效果：
本发明所设计的基于量子地雷爆炸机理的时频DOA估计方法，涉及智能计算和阵列信号处理两个领域。本发明针对现有的时频子空间分解类DOA估计方法无法对相干信号准确测向，而时频极大似然DOA估计方法运算量大、实现复杂的难题，设计了一种能快速获得高精度时频DoA估计结果，且能对相干源进行有效测向的DOA估计方法，本发明基于量子地雷爆炸机理，实现了对非平稳信号快速、稳健的DOA估计，且可以对多个相干的非平稳信号进行DOA估计
本发明针对现有的时频子空间分解类的DOA估计方法无法对相干信号准确测向，而时频极大似然DOA估计方法运算量大的难题，设计了一种基于量子地雷爆炸机理的时频DOA估计方法，实现了对非平稳信号快速、稳健的DOA估计，且可以对多个相干的非平稳信号进行DOA估计，避免了对时变方向向量进行汇聚、插值等复杂的变换。利用地雷爆炸机理收敛速度快的优点，并引入量子原理，不仅使用量子编码，而且设计了量子演化方程，使用模拟的量子旋转门演化量子地雷位置，在减少计算量的同时获得了性能的提高，实现了在较短时间内得到较准确的DOA估计结果。
多数DOA估计方法是基于数据协方差矩阵的估计，本发明所设计的基于量子地雷爆炸机理的时频DOA估计方法是利用快拍采样数据的时频平均的PWVD矩阵进行估计，利用了二次型时频分布的时频聚集特性，构造极大似然时频测向方法，解决了高斯白噪声环境下的非平稳信号DOA估计问题。
本发明利用了时频极大似然DOA估计的优势，在相干信源情况下仍能有效估计来波方向，同时，为了减少极大似然估计所需的计算量，快速准确地获得时频DOA估计结果，本发明使用量子地雷爆炸机理进行寻优搜索，相比于传统的粒子群等优化方法具有更高的估计精度。
附图说明
图1为本发明所设计的基于量子地雷爆炸机理的时频DOA估计方法示意图；
图2为两个线性调频信号经阵列流型矩阵合成后的PWVD；
图3为信号估计角度的均方根误差与信噪比的关系曲线；
图4为两个全相干的调频信号经阵列流型矩阵合成后的PWVD；
图5为使用基于量子地雷爆炸机理的时频DOA估计方法得到的两个信号来波方向的角度估计值；
图6为使用粒子群极大似然时频DOA估计方法得到的两个信号来波方向的角度估计值。
具体实施方式
下面结合附图对本发明具体实施方式做进一步说明。
如图1所示，本发明包括以下步骤：
步骤一，建立阵列接收的时域数据模型。
在信号为远场窄带信号、阵列形式为均匀线阵且阵元各向同性的假设下，在t时刻，若N个方位角为θ＝[θ1,θ2,…,θN]的调频信号入射到含有M个阵元的阵列上，其中第n个调频信号表示为sn(t)，n＝1,2,…,N，则阵列上第m个阵元接收到的数据
Figure BDA0002043512500000061
m＝1,2,…,M，其中τmn为第n个信号入射到第m个阵元时，相对于入射到参考阵元的空间延迟，vm(t)为第m个阵元处在t时刻与信号独立的高斯白噪声。将M个阵元在t时刻接收到的数据排成列矢量
Figure BDA0002043512500000062
步骤二，对时域数据进行快拍采样。
对阵列中M个阵元接收的数据进行K次快拍采样，得到快拍采样数据矩阵，
Figure BDA0002043512500000063
ym(k)为第m个阵元在第k个快拍采样时刻的数据，m＝1,2,…,M，k＝1,2,…K，也可表示为Y＝A(θ)S+V，其中
Figure BDA0002043512500000071
为阵列流型矩阵，ωn，n＝1,2,…,N，为快拍采样时间点处第n个调频信号的瞬时角频率，在足够短的时间区间内，可以认为ωn近似不变，S是大小为N×K的信号矩阵，V是大小为M×K的噪声矩阵。
步骤三，对快拍采样数据进行时频分析得到快拍采样数据的PWVD矩阵。
使用离散时间形式的Cohen类时频分布中的PWVD(pseudo-Wigner-Villedistribution)对快拍采样数据进行时频分析。在信号的时频域平面上，第n个信号的第k个快拍采样时间点tk及此时信号的瞬时频率fn,k确定时频点(tk,fn,k)，n＝1,2,…,N，k＝1,2,…,K，则任意两阵元的快拍采样数据使用长度为L的矩形窗的PWVD在时频点(tk,fn,k)处的值
Figure BDA0002043512500000072
σ为滞后变量，w＝1,2,…,M，u＝1,2,…,M，(·)H代表共轭转置运算，其中矩形窗长度L为奇数且应足够小以保证在观察区间内阵列流型矩阵是近似非时变的。将各个阵元快拍采样数据的自PWVD项及互PWVD项排成矩阵，得到时频点(tk,fn,k)处快拍采样数据的PWVD矩阵
Figure BDA0002043512500000073
步骤四，计算时频平均的快拍采样数据PWVD矩阵。
选取多个时频点进行时频平均有利于更充分地利用时频信息，进而提高有效信噪比。由于快拍采样数据的PWVD在时频点处的数值在数据段的中间趋于常数，在观察时间的两端则性能下降，为减小边缘效应的影响，略去数据段两端的PWVD数据，在每个信号的全部K个时频点中，保留中间的K-L+1个时频点，对这些时频点处的PWVD数值进行算术平均，得到时频平均后的快拍采样数据PWVD矩阵
Figure BDA0002043512500000074
步骤五，构造时频极大似然DoA估计的极大似然方程。
极大似然方程为
Figure BDA0002043512500000075
其中
Figure BDA0002043512500000076
为信号方位角的最优估计值，A(θ)为解空间中一个可能的解对应的导向矢量矩阵，tr代表矩阵求迹运算。
步骤六，初始化量子地雷量子位置。
设定量子地雷的总数量P，量子地雷量子位置的演化寻优过程分为勘探过程和探测过程，分别设定勘探过程的最大迭代次数G1，探测过程的最大迭代次数G2，当前迭代数为g∈[1,G1+G2]，第p个量子地雷的N维初代量子位置为
Figure BDA0002043512500000081
p＝1,2,…,P，初代量子地雷量子位置每一维取值为[0,1]区间内的均匀随机数。
步骤七，构造适应度函数，计算每个量子地雷量子位置适应度，根据适应度值获得全局最优量子位置。
在第g次迭代中，将量子地雷量子位置
Figure BDA0002043512500000082
映射到信号方位角的解空间
Figure BDA0002043512500000083
范围内，得到量子地雷量子位置的映射态
Figure BDA0002043512500000084
p＝1,2,…,P，构造其适应度函数为
Figure BDA0002043512500000085
p＝1,2,…,P，其中
Figure BDA0002043512500000086
为第g次迭代中第p个量子地雷量子位置的映射态对应的阵列流型矩阵。计算所有量子地雷量子位置适应度，适应度值最大的量子地雷量子位置为全局最优量子位置记作
Figure BDA0002043512500000087
b∈[1,P]。
步骤八，量子地雷爆炸产生量子弹片，根据模拟量子旋转门模拟量子地雷爆炸过程获得量子弹片的量子位置。
每个量子地雷爆炸产生的量子弹片个数为Q，使用模拟量子旋转门计算第p个量子地雷爆炸产生的第q个量子弹片的第n维量子位置为
Figure BDA0002043512500000088
p＝1,2,…,P，q＝1,2,…,Q，n＝1,2,…,N，在勘探过程中g∈[1,G1]，第p个量子地雷的第n维量子位置对应的量子旋转角为
Figure BDA0002043512500000089
其中
Figure BDA00020435125000000810
为[0,1]区间内的均匀随机数；在探测过程中g∈[G1+1,G1+G2]，第p个量子地雷的第n维量子位置对应的量子旋转角为
Figure BDA00020435125000000811
p≠b；另外规定当p＝b时，
Figure BDA00020435125000000812
其值取正或负的概率相等，其中
Figure BDA00020435125000000813
为[0,1]区间内的均匀随机数。
步骤九，对于每个量子地雷爆炸产生的量子弹片量子位置，计算量子弹片量子位置映射态的适应度函数值，应用贪婪选择，选择适应度大的优秀量子位置作为放置量子地雷的量子位置，并在此量子位置上引爆下一代的量子地雷，并根据所有量子位置的适应度更新全局最优量子位置。
计算量子弹片量子位置的映射态
Figure BDA0002043512500000091
p＝1,2,…,P，q＝1,2,…,Q，量子弹片量子位置的映射态和量子地雷量子位置的映射态的适应度函数定义相同，计算第p个量子地雷爆炸产生的第q个量子弹片量子位置映射态的适应度函数值
Figure BDA0002043512500000092
p＝1,2,…,P，q＝1,2,…,Q，为保证迭代寻优过程中量子地雷量子位置的适应度值随迭代次数的增大单调递增，应用贪婪选择，选择最优量子位置用于引爆下一代量子地雷，逐个判断量子弹片量子位置的适应度值是否大于当前量子地雷量子位置的适应度值
Figure BDA0002043512500000093
若
Figure BDA0002043512500000094
则有
Figure BDA0002043512500000095
q＝1,2,…,Q；否则
Figure BDA0002043512500000096
p＝1,2,…,P。将全局最优量子位置更新为适应度最大的量子地雷量子位置
Figure BDA0002043512500000097
步骤十，判断是否达到最大迭代次数，若g＜G1+G2则令g＝g+1，返回步骤八。若g＝G1+G2，则将最后一代中最优的量子地雷量子位置的映射态作为信号方位角的最优估计值输出。
基于量子地雷爆炸机理的时频DOA估计方法的参数设置如下：N＝2，M＝8，K＝1024，Q＝10，L＝103，使用阵元间距为二分之一信号波长的均匀线阵。
图2、图3中，两个线性调频信号的归一化瞬时频率分别为
Figure BDA0002043512500000098
和
Figure BDA0002043512500000099
k＝1,2,…,K，来波方向分别为θ1＝-10°和θ2＝20°，图2为无噪声情况下参考阵元输出的PWVD。图3为基于量子地雷爆炸机理的时频DOA估计方法与粒子群极大似然时频DOA估计方法性能对比，基于量子地雷爆炸机理的时频DOA估计方法的参数设置如下：P＝50，G1＝10，G2＝50，粒子群极大似然时频DOA估计的个体数为50，最大迭代次数为60，惯性权重为0.8，Monte Carlo实验次数为100次。由实验结果显示出量子地雷爆炸机理比粒子群搜索具有更高的估计精度。
在图4、图5和图6中，两个全相干的调频信号的来波方向分别为-20°和20°，全相干调频信号的归一化瞬时频率为
Figure BDA0002043512500000101
k＝1,2,…,K，无噪声情况下参考阵元输出的PWVD如图4所示，信噪比均为10dB的情况下，使用基于量子地雷爆炸机理的时频DOA估计方法得到的两个信号来波方向的角度估计值如图5所示，参数设置如下：P＝30，G1＝10，G2＝50，使用粒子群极大似然时频DOA估计方法得到的两个信号来波方向的角度估计值如图6所示，粒子群极大似然DOA估计的个体数为30，最大迭代次数为60，惯性权重为0.8，Monte Carlo实验次数均为20次。对比图5、6可以看出，本发明设计的方法能对相干源的来波方向进行有效估计，且能更好地得到近似全局最优解。

Claims (1)

1.一种时频DOA估计方法，其特征在于，包括以下步骤：
步骤一：建立阵列接收的时域数据模型：
在信号为远场窄带信号、阵列形式为均匀线阵且阵元各向同性的假设下，在t时刻，若N个方位角为θ＝[θ1,θ2,…,θN]的调频信号入射到含有M个阵元的阵列上，其中第n个调频信号表示为sn(t)，n＝1,2,…,N，则阵列上第m个阵元接收到的数据ym(t)满足：
Figure FDA0002043512490000011
其中τmn为第n个信号入射到第m个阵元时相对于入射到参考阵元的空间延迟，vm(t)为第m个阵元处在t时刻与信号独立的高斯白噪声，将M个阵元在t时刻接收到的数据排成列矢量
Figure FDA0002043512490000012
步骤二：获取快拍采样数据：
对阵列中M个阵元接收的数据进行K次快拍采样，得到快拍采样数据矩阵：
Figure FDA0002043512490000013
ym(k)为第m个阵元在第k个快拍采样时刻的数据，m＝1,2,…,M，k＝1,2,…K，Y也可表示为Y＝A(θ)S+V，其中
Figure FDA0002043512490000014
为阵列流型矩阵，ωn，n＝1,2,…,N，为快拍采样时间点处第n个调频信号的瞬时角频率，S是大小为N×K的信号矩阵，V是大小为M×K的噪声矩阵；
步骤三：对快拍采样数据进行时频分析：
使用离散时间形式的Cohen类时频分布中的PWVD对快拍采样数据进行时频分析，在信号的时频域平面上，第n个信号的第k个快拍采样时间点tk及此时信号的瞬时频率fn,k确定时频点(tk,fn,k)，n＝1,2,…,N，k＝1,2,…,K，则任意两阵元的快拍采样数据使用长度为L的矩形窗的PWVD在时频点(tk,fn,k)处的值满足：
Figure FDA0002043512490000015
σ为滞后变量，w＝1,2,…,M，u＝1,2,…,M，(·)H代表共轭转置运算，其中矩形窗长度L为奇数且应足够小以保证在观察区间内阵列流型矩阵是近似非时变的；将各个阵元快拍采样数据的自PWVD项及互PWVD项排成矩阵，得到时频点(tk,fn,k)处快拍采样数据的PWVD矩阵：
Figure FDA0002043512490000021
步骤四：计算时频平均的快拍采样数据PWVD矩阵：
略去数据段两端的PWVD数据，在每个信号的全部K个时频点中，保留中间的K-L+1个时频点，对这些时频点处的PWVD数值进行算术平均，得到时频平均后的快拍采样数据PWVD矩阵
Figure FDA0002043512490000022
步骤五：构造时频极大似然DOA估计的极大似然方程：
极大似然方程为
Figure FDA0002043512490000023
其中
Figure FDA0002043512490000024
为信号方位角的最优估计值，A(θ)为解空间中一个可能的解对应的导向矢量矩阵，tr代表矩阵求迹运算；
步骤六：量子地雷量子位置的初始化：
设定量子地雷的总数量P，量子地雷量子位置的演化寻优过程分为勘探过程和探测过程，分别设定勘探过程的最大迭代次数G1，探测过程的最大迭代次数G2，当前迭代数为g∈[1,G1+G2]，第p个量子地雷的N维初代量子位置为
Figure FDA0002043512490000025
初代量子地雷量子位置每一维取值为[0,1]区间内的均匀随机数；
步骤七：由极大似然方程构造适应度函数，计算每个量子地雷量子位置适应度，根据适应度值获得全局最优量子位置：
在第g次迭代中，将量子地雷量子位置
Figure FDA0002043512490000026
映射到信号方位角的解空间
Figure FDA0002043512490000027
范围内，得到量子地雷量子位置的映射态
Figure FDA0002043512490000028
构造其适应度函数为
Figure FDA0002043512490000029
其中
Figure FDA00020435124900000210
为第g次迭代中第p个量子地雷量子位置的映射态对应的阵列流型矩阵，计算所有量子地雷量子位置适应度，适应度值最大的量子地雷量子位置为全局最优量子位置，记作
Figure FDA0002043512490000031
步骤八：量子地雷爆炸产生量子弹片，使用模拟量子旋转门模拟量子地雷爆炸过程获得量子弹片的量子位置：
每个量子地雷爆炸产生的量子弹片个数为Q，使用模拟量子旋转门计算第p个量子地雷爆炸产生的第q个量子弹片的第n维量子位置为
Figure FDA0002043512490000032
在勘探过程中g∈[1,G1]，第p个量子地雷的第n维量子位置对应的量子旋转角为
Figure FDA0002043512490000033
其中
Figure FDA0002043512490000034
为[0,1]区间内的均匀随机数；在探测过程中g∈[G1+1,G1+G2]，第p个量子地雷的第n维量子位置对应的量子旋转角为
Figure FDA0002043512490000035
p≠b；另外规定当p＝b时，
Figure FDA0002043512490000036
其值取正或负的概率相等，其中
Figure FDA0002043512490000037
为[0,1]区间内的均匀随机数；
步骤九：对于每个量子地雷爆炸产生的量子弹片量子位置，计算量子弹片量子位置映射态的适应度函数值，应用贪婪选择机制，选择适应度大的优秀量子位置作为放置量子地雷的量子位置，并在此量子位置上引爆下一代的量子地雷，并根据所有量子位置的适应度更新全局最优量子位置；
计算量子弹片量子位置的映射态
Figure FDA0002043512490000038
Figure FDA0002043512490000039
量子弹片量子位置的映射态和量子地雷量子位置的映射态的适应度函数定义相同，计算第p个量子地雷爆炸产生的第q个量子弹片量子位置映射态的适应度函数值
Figure FDA00020435124900000310
应用贪婪选择，选择最优量子位置用于引爆下一代量子地雷，逐个判断量子弹片量子位置的适应度值是否大于当前量子地雷量子位置的适应度值
Figure FDA00020435124900000311
若
Figure FDA00020435124900000312
则有
Figure FDA00020435124900000313
否则
Figure FDA00020435124900000314
Figure FDA00020435124900000315
将全局最优量子位置更新为适应度最大的量子地雷量子位置
Figure FDA00020435124900000316
步骤十：判断是否达到最大迭代次数，若g＜G1+G2则令g＝g+1，返回步骤八；若g＝G1+G2，则将最后一代中最优的量子地雷量子位置的映射态作为信号方位角的最优估计值输出。
CN201910349676.4A 2019-04-28 2019-04-28 一种时频doa估计方法 Active CN110046326B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910349676.4A CN110046326B (zh) 2019-04-28 2019-04-28 一种时频doa估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910349676.4A CN110046326B (zh) 2019-04-28 2019-04-28 一种时频doa估计方法

Publications (2)

Publication Number Publication Date
CN110046326A CN110046326A (zh) 2019-07-23
CN110046326B true CN110046326B (zh) 2022-09-27

Family

ID=67280038

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910349676.4A Active CN110046326B (zh) 2019-04-28 2019-04-28 一种时频doa估计方法

Country Status (1)

Country Link
CN (1) CN110046326B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113111304B (zh) * 2021-04-01 2022-09-27 哈尔滨工程大学 强冲击噪声下基于量子射线机理的相干分布源测向方法
CN113759303B (zh) * 2021-08-04 2024-05-24 中山大学 一种基于粒子群算法的无网格波达角估计方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5610612A (en) * 1994-03-07 1997-03-11 Piper; John E. Method for maximum likelihood estimations of bearings
CN104616059A (zh) * 2015-01-29 2015-05-13 无锡职业技术学院 一种基于量子粒子群的波达方向估计方法
CN109633558A (zh) * 2018-10-25 2019-04-16 上海无线电设备研究所 一种基于极化时频分布的波达方向估计算法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5610612A (en) * 1994-03-07 1997-03-11 Piper; John E. Method for maximum likelihood estimations of bearings
CN104616059A (zh) * 2015-01-29 2015-05-13 无锡职业技术学院 一种基于量子粒子群的波达方向估计方法
CN109633558A (zh) * 2018-10-25 2019-04-16 上海无线电设备研究所 一种基于极化时频分布的波达方向估计算法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Direction of Arrival Estimation of Nonstationary;Siyu Chen等;《IEEE International Symposium on Microwave》;20160714;749-752 *
Time-frequency maximum likelihood methods;Yimin Zhang,Weifeng Mu,Moeness G.Amin*;《Journal of The Franklin Institute》;20000731(第337期);483-497 *
基于极化时频分布的改进DOA估计算法;李芬等;《上海航天》;20180430;第35卷(第4期);108-113 *

Also Published As

Publication number Publication date
CN110046326A (zh) 2019-07-23

Similar Documents

Publication Publication Date Title
Melvin et al. Improving practical space-time adaptive radar
CN103399292B (zh) 一种基于软稀疏表示的doa估计方法
Shahbazpanahi et al. A covariance fitting approach to parametric localization of multiple incoherently distributed sources
CN103744076B (zh) 基于非凸优化的mimo雷达动目标检测方法
CN109061554A (zh) 一种基于空间离散网格动态更新的目标到达角度估计方法
CN108375751A (zh) 多信源波达方向估计方法
CN106021637A (zh) 互质阵列中基于迭代稀疏重构的doa估计方法
CN105699950B (zh) 基于自适应迭代前后向平滑共轭梯度的雷达杂波抑制方法
CN109239646B (zh) 一种冲击噪声环境下连续量子水蒸发的二维动态测向方法
CN105403874A (zh) 非均匀阵列欠定波达方向估计方法
CN110046326B (zh) 一种时频doa估计方法
Tao et al. A knowledge aided SPICE space time adaptive processing method for airborne radar with conformal array
CN112346030A (zh) 无人机群的超分辨波达方向估计方法
Qi et al. Time-frequency DOA estimation of chirp signals based on multi-subarray
CN109212466B (zh) 一种基于量子蜻蜓演化机制的宽带测向方法
CN113671485B (zh) 基于admm的米波面阵雷达二维doa估计方法
CN112904298B (zh) 一种基于局部网格分裂的网格偏离空时自适应处理方法
CN113759303A (zh) 一种基于粒子群算法的无网格波达角估计方法
Sun et al. Airborne radar STAP using sparse recovery of clutter spectrum
Reaz et al. A comprehensive analysis and performance evaluation of different direction of arrival estimation algorithms
Hua et al. Target detection in sea clutter via weighted averaging filter on the Riemannian manifold
CN112800596A (zh) 强冲击噪声下基于嵌套阵列的鲁棒动态测向方法
TWI677696B (zh) 一種雷達估測目標方位之方法
Li et al. Space-Time Adaptive Processing Using Convolutional Neural Network
CN112363106B (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
TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20230727

Address after: Room 203, Building 1, Science and Technology Entrepreneurship Service Center, Nanhu District, Jiaxing City, Zhejiang Province, 314000

Patentee after: JIAXING NUOAIDI COMMUNICATION TECHNOLOGY CO.,LTD.

Address before: 150001 Intellectual Property Office, Harbin Engineering University science and technology office, 145 Nantong Avenue, Nangang District, Harbin, Heilongjiang

Patentee before: HARBIN ENGINEERING University