[发明专利]一种时频DOA估计方法

权利要求书

1.一种时频DOA估计方法，其特征在于，包括以下步骤：

步骤一：建立阵列接收的时域数据模型：

在信号为远场窄带信号、阵列形式为均匀线阵且阵元各向同性的假设下，在t时刻，若N个方位角为θ＝[θ₁,θ₂,…,θ_N]的调频信号入射到含有M个阵元的阵列上，其中第n个调频信号表示为s_n(t)，n＝1,2,…,N，则阵列上第m个阵元接收到的数据y_m(t)满足：其中τ_mn为第n个信号入射到第m个阵元时相对于入射到参考阵元的空间延迟，v_m(t)为第m个阵元处在t时刻与信号独立的高斯白噪声，将M个阵元在t时刻接收到的数据排成列矢量

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

对阵列中M个阵元接收的数据进行K次快拍采样，得到快拍采样数据矩阵：y_m(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个快拍采样时间点t_k及此时信号的瞬时频率f_n,k确定时频点(t_k,f_n,k)，n＝1,2,…,N，k＝1,2,…,K，则任意两阵元的快拍采样数据使用长度为L的矩形窗的PWVD在时频点(t_k,f_n,k)处的值满足：σ为滞后变量，w＝1,2,…,M，u＝1,2,…,M，(·)^H代表共轭转置运算，其中矩形窗长度L为奇数且应足够小以保证在观察区间内阵列流型矩阵是近似非时变的；将各个阵元快拍采样数据的自PWVD项及互PWVD项排成矩阵，得到时频点(t_k,f_n,k)处快拍采样数据的PWVD矩阵：

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

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

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

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

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

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

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

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

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

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

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

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

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