一种时频DOA估计方法与流程

本发明涉及一种时频DOA估计方法，特别是一种基于量子地雷爆炸机理的时频DOA估计方法，属于阵列信号处理领域。

背景技术：

波达方向(DOA)估计是阵列信号处理中的一个重要研究方向，在雷达、声呐、无线电天文学等领域中有广泛的应用。当前，绝大多数DoA估计方法都是针对平稳信号进行估计的，由于传统的平稳信号DOA估计方法局限于空域处理，对于非平稳信号的DOA估计则性能恶化严重，甚至失效，所以需设计专门针对非平稳信号的DOA估计方法。

对于非平稳信号的分析，分数阶傅里叶变换，循环谱估计，小波变换，时频分析等都是行之有效的方法，其中，时频分析方法的基本思想是通过设计时间和频率的联合函数，同时描述信号在时频域上不同时间及频率处的能量密度，这种方法实现简单，且具有优良的时频聚集性，将时频分析方法与阵列信号处理结合的进行时频DOA估计，是一种解决非平稳信号DoA估计问题的有效方法。

如何快速高精度获得时频DOA估计结果是一个经典的工程难题，在众多时频DOA估计方法中，使用极大似然原理的时频DOA估计是一种性能优良的方法，且可以对相干源进行有效估计，因此，尽管时频极大似然DOA估计的实现过程较为复杂，但是这种方法仍有广阔的应用前景。

为了降低时频极大似然DOA估计的复杂度，需要进行新方法的设计进而有效求解。当前，许多智能计算方法被应用于求解复杂的优化问题，其中地雷爆炸搜索是一种启发式智能计算方法，受现实世界中地雷爆炸产生弹片，引爆其他潜在地雷的现象启发，具有新颖、有效的特点，但易陷入局部收敛，因此需要改进地雷爆炸方法使其具有收敛速度快、搜索精度高的优点。

通过对现有技术文献的检索发现，Yimin Zhang等在《Journal of the Franklin Institute》(337(2000)483-497)上发表的“Time-frequency maximum likelihood methods for direction finding”中提出了一种时频极大似然DOA估计方法，可以解相干且不需要进行使阵列孔径损失的预处理，但没有解决极大似然估计实现过程复杂的难题。Siyu Chen等在《IEEE International Symposium on Microwave》(10.1109/MAPE.2015.7510425)上发表的“Direction of Arrival Estimation of Nonstationary Signal 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)满足：m＝1,2,…,M，其中τmn为第n个信号入射到第m个阵元时相对于入射到参考阵元的空间延迟，vm(t)为第m个阵元处在t时刻与信号独立的高斯白噪声，将M个阵元在t时刻接收到的数据排成列矢量

步骤二：获取快拍采样数据：

对阵列中M个阵元接收的数据进行K次快拍采样，得到快拍采样数据矩阵：ym(k)为第m个阵元在第k个快拍采样时刻的数据，m＝1,2,…,M，k＝1,2,….K，Y也可表示为Y＝A(θ)S+V，其中为阵列流型矩阵，ω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)处的值满足：σ为滞后变量，w＝1,2,…,M，u＝1,2,…,M，(·)^H代表共轭转置运算，其中矩形窗长度L为奇数且应足够小以保证在观察区间内阵列流型矩阵是近似非时变的；将各个阵元快拍采样数据的自PWVD项及互PWVD项排成矩阵，得到时频点(tk,fn,k)处快拍采样数据的PWVD矩阵：

步骤四：计算时频平均的快拍采样数据PWVD矩阵：

略去数据段两端的PWVD数据，在每个信号的全部K个时频点中，保留中间的K-L+1个时频点，对这些时频点处的PWVD数值进行算术平均，得到时频平均后的快拍采样数据PWVD矩阵

步骤五：构造时频极大似然DOA估计的极大似然方程：

极大似然方程为其中为信号方位角的最优估计值，A(θ)为解空间中一个可能的解对应的导向矢量矩阵，tr代表矩阵求迹运算；

步骤六：量子地雷量子位置的初始化：

设定量子地雷的总数量P，量子地雷量子位置的演化寻优过程分为勘探过程和探测过程，分别设定勘探过程的最大迭代次数G1，探测过程的最大迭代次数G2，当前迭代数为g∈[1,G1+G2]，第p个量子地雷的N维初代量子位置为p＝1,2,…,P，初代量子地雷量子位置每一维取值为[0,1]区间内的均匀随机数；

步骤七：由极大似然方程构造适应度函数，计算每个量子地雷量子位置适应度，根据适应度值获得全局最优量子位置：

在第g次迭代中，将量子地雷量子位置映射到信号方位角的解空间范围内，得到量子地雷量子位置的映射态p＝1,2,…,P，构造其适应度函数为p＝1,2,…,P，其中为第g次迭代中第p个量子地雷量子位置的映射态对应的阵列流型矩阵，计算所有量子地雷量子位置适应度，适应度值最大的量子地雷量子位置为全局最优量子位置，记作b∈[1,P]；

步骤八：量子地雷爆炸产生量子弹片，使用模拟量子旋转门模拟量子地雷爆炸过程获得量子弹片的量子位置：

每个量子地雷爆炸产生的量子弹片个数为Q，使用模拟量子旋转门计算第p个量子地雷爆炸产生的第q个量子弹片的第n维量子位置为p＝1,2,…,P，q＝1,2,…,Q，n＝1,2,…,N，在勘探过程中g∈[1,G1]，第p个量子地雷的第n维量子位置对应的量子旋转角为其中为[0,1]区间内的均匀随机数；在探测过程中g∈[G1+1,G1+G2]，第p个量子地雷的第n维量子位置对应的量子旋转角为p≠b；另外规定当p＝b时，其值取正或负的概率相等，其中为[0,1]区间内的均匀随机数；

步骤九：对于每个量子地雷爆炸产生的量子弹片量子位置，计算量子弹片量子位置映射态的适应度函数值，应用贪婪选择机制，选择适应度大的优秀量子位置作为放置量子地雷的量子位置，并在此量子位置上引爆下一代的量子地雷，并根据所有量子位置的适应度更新全局最优量子位置；

计算量子弹片量子位置的映射态p＝1,2,…,P，q＝1,2,…,Q，量子弹片量子位置的映射态和量子地雷量子位置的映射态的适应度函数定义相同，计算第p个量子地雷爆炸产生的第q个量子弹片量子位置映射态的适应度函数值p＝1,2,…,P，q＝1,2,…,Q，应用贪婪选择，选择最优量子位置用于引爆下一代量子地雷，逐个判断量子弹片量子位置的适应度值是否大于当前量子地雷量子位置的适应度值若则有q＝1,2,…,Q；否则p＝1,2,…,P；将全局最优量子位置更新为适应度最大的量子地雷量子位置

步骤十：判断是否达到最大迭代次数，若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个阵元接收到的数据m＝1,2,…,M，其中τmn为第n个信号入射到第m个阵元时，相对于入射到参考阵元的空间延迟，vm(t)为第m个阵元处在t时刻与信号独立的高斯白噪声。将M个阵元在t时刻接收到的数据排成列矢量

步骤二，对时域数据进行快拍采样。

对阵列中M个阵元接收的数据进行K次快拍采样，得到快拍采样数据矩阵，ym(k)为第m个阵元在第k个快拍采样时刻的数据，m＝1,2,…,M，k＝1,2,….K，也可表示为Y＝A(θ)S+V，其中为阵列流型矩阵，ωn，n＝1,2,…,N，为快拍采样时间点处第n个调频信号的瞬时角频率，在足够短的时间区间内，可以认为ωn近似不变，S是大小为N×K的信号矩阵，V是大小为M×K的噪声矩阵。

步骤三，对快拍采样数据进行时频分析得到快拍采样数据的PWVD矩阵。

使用离散时间形式的Cohen类时频分布中的PWVD(pseudo-Wigner-Ville distribution)对快拍采样数据进行时频分析。在信号的时频域平面上，第n个信号的第k个快拍采样时间点tk及此时信号的瞬时频率fn,k确定时频点(tk,fn,k)，n＝1,2,…,N，k＝1,2,…,K，则任意两阵元的快拍采样数据使用长度为L的矩形窗的PWVD在时频点(tk,fn,k)处的值σ为滞后变量，w＝1,2,…,M，u＝1,2,…,M，(·)^H代表共轭转置运算，其中矩形窗长度L为奇数且应足够小以保证在观察区间内阵列流型矩阵是近似非时变的。将各个阵元快拍采样数据的自PWVD项及互PWVD项排成矩阵，得到时频点(tk,fn,k)处快拍采样数据的PWVD矩阵

步骤四，计算时频平均的快拍采样数据PWVD矩阵。

选取多个时频点进行时频平均有利于更充分地利用时频信息，进而提高有效信噪比。由于快拍采样数据的PWVD在时频点处的数值在数据段的中间趋于常数，在观察时间的两端则性能下降，为减小边缘效应的影响，略去数据段两端的PWVD数据，在每个信号的全部K个时频点中，保留中间的K-L+1个时频点，对这些时频点处的PWVD数值进行算术平均，得到时频平均后的快拍采样数据PWVD矩阵

步骤五，构造时频极大似然DoA估计的极大似然方程。

极大似然方程为其中为信号方位角的最优估计值，A(θ)为解空间中一个可能的解对应的导向矢量矩阵，tr代表矩阵求迹运算。

步骤六，初始化量子地雷量子位置。

设定量子地雷的总数量P，量子地雷量子位置的演化寻优过程分为勘探过程和探测过程，分别设定勘探过程的最大迭代次数G1，探测过程的最大迭代次数G2，当前迭代数为g∈[1,G1+G2]，第p个量子地雷的N维初代量子位置为p＝1,2,…,P，初代量子地雷量子位置每一维取值为[0,1]区间内的均匀随机数。

步骤七，构造适应度函数，计算每个量子地雷量子位置适应度，根据适应度值获得全局最优量子位置。

在第g次迭代中，将量子地雷量子位置映射到信号方位角的解空间范围内，得到量子地雷量子位置的映射态p＝1,2,…,P，构造其适应度函数为p＝1,2,…,P，其中为第g次迭代中第p个量子地雷量子位置的映射态对应的阵列流型矩阵。计算所有量子地雷量子位置适应度，适应度值最大的量子地雷量子位置为全局最优量子位置记作b∈[1,P]。

步骤八，量子地雷爆炸产生量子弹片，根据模拟量子旋转门模拟量子地雷爆炸过程获得量子弹片的量子位置。

每个量子地雷爆炸产生的量子弹片个数为Q，使用模拟量子旋转门计算第p个量子地雷爆炸产生的第q个量子弹片的第n维量子位置为p＝1,2,…,P，q＝1,2,…,Q，n＝1,2,…,N，在勘探过程中g∈[1,G1]，第p个量子地雷的第n维量子位置对应的量子旋转角为其中为[0,1]区间内的均匀随机数；在探测过程中g∈[G1+1,G1+G2]，第p个量子地雷的第n维量子位置对应的量子旋转角为p≠b；另外规定当p＝b时，其值取正或负的概率相等，其中为[0,1]区间内的均匀随机数。

步骤九，对于每个量子地雷爆炸产生的量子弹片量子位置，计算量子弹片量子位置映射态的适应度函数值，应用贪婪选择，选择适应度大的优秀量子位置作为放置量子地雷的量子位置，并在此量子位置上引爆下一代的量子地雷，并根据所有量子位置的适应度更新全局最优量子位置。

计算量子弹片量子位置的映射态p＝1,2,…,P，q＝1,2,…,Q，量子弹片量子位置的映射态和量子地雷量子位置的映射态的适应度函数定义相同，计算第p个量子地雷爆炸产生的第q个量子弹片量子位置映射态的适应度函数值p＝1,2,…,P，q＝1,2,…,Q，为保证迭代寻优过程中量子地雷量子位置的适应度值随迭代次数的增大单调递增，应用贪婪选择，选择最优量子位置用于引爆下一代量子地雷，逐个判断量子弹片量子位置的适应度值是否大于当前量子地雷量子位置的适应度值若则有q＝1,2,…,Q；否则p＝1,2,…,P。将全局最优量子位置更新为适应度最大的量子地雷量子位置

步骤十，判断是否达到最大迭代次数，若g＜G1+G2则令g＝g+1，返回步骤八。若g＝G1+G2，则将最后一代中最优的量子地雷量子位置的映射态作为信号方位角的最优估计值输出。

基于量子地雷爆炸机理的时频DOA估计方法的参数设置如下：N＝2，M＝8，K＝1024，Q＝10，L＝103，使用阵元间距为二分之一信号波长的均匀线阵。

图2、图3中，两个线性调频信号的归一化瞬时频率分别为和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°，全相干调频信号的归一化瞬时频率为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可以看出，本发明设计的方法能对相干源的来波方向进行有效估计，且能更好地得到近似全局最优解。