CN114167347A - 冲击噪声环境下互质阵列的幅相误差校正和测向方法 - Google Patents

冲击噪声环境下互质阵列的幅相误差校正和测向方法 Download PDF

Info

Publication number
CN114167347A
CN114167347A CN202111421629.XA CN202111421629A CN114167347A CN 114167347 A CN114167347 A CN 114167347A CN 202111421629 A CN202111421629 A CN 202111421629A CN 114167347 A CN114167347 A CN 114167347A
Authority
CN
China
Prior art keywords
quantum
array
eagle
harris
matrix
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
CN202111421629.XA
Other languages
English (en)
Other versions
CN114167347B (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.)
Harbin Engineering University
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 CN202111421629.XA priority Critical patent/CN114167347B/zh
Publication of CN114167347A publication Critical patent/CN114167347A/zh
Application granted granted Critical
Publication of CN114167347B publication Critical patent/CN114167347B/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
    • 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
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • 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

Abstract

本发明提供一种冲击噪声环境下互质阵列的幅相误差校正和测向方法，首先采用外加辅助源的校正算法得出幅相误差的粗估计值，再用量子哈里斯鹰算法在粗估计相位误差周围进行搜索，可以实现在极低信噪比下对幅相误差进行更精确的估计。同时，在冲击噪声下互质阵列的波达方向估计问题上，本发明所设计的基于量子哈里斯鹰机制的分数低阶协方差结合虚拟矩阵的极大似然测向方法，可在相同信噪比下取得比其他传统算法更低的均方根误差，其中引入的虚拟阵列和空间平滑算法，可有效提高互质阵列的空间自由度。

Description

冲击噪声环境下互质阵列的幅相误差校正和测向方法
技术领域
本发明涉及的是一种针对互质阵列在冲击噪声下的幅相误差校正和测向方法，属于阵列信号处理领域。
背景技术
阵列信号的测向是阵列信号处理中一项研究重点，可广泛应用于军事和民用通信中。随着现代民用通信和军事通信的电磁环境日益复杂，干扰手段越来越多，对接收阵列的信号处理能力的要求也越来越高。在传统均匀线阵中DOA估计算法最多可估计出小于信源个数的波达角，其空间自由度较小。因此近几年提出一种拓展互质阵列，在估计波达角时，可通过构建其虚拟均匀阵列的方法，估计出大于原互质阵元个数的信源波达角数，从而达到增加其空间自由度的目的。在实际应用当中，接收阵列可能存在各种误差，其中最常见的幅相误差对于测向算法的精确度影响较大，然而在大部分互质阵列的幅相误差校正和测向方法中，均假定互质阵列处于高斯白噪声环境之下，真实通信场景却可能存在大量冲击噪声，会使此类方法的准确性降低，甚至完全失效。因此研究在冲击噪声下对互质阵列幅相误差的校正方法和测向方法具有重要意义和价值。
单辅助源幅相误差校正算法是最经典的幅相误差校正的方法，可以将其推广至互质阵列的幅相误差校正上。但单辅助源的幅相误差校正方法要求较大的信噪比，否则无法对幅相误差进行精确估计，且通常情况下仅适用于高斯噪声环境下；同样地，基于虚拟阵列和极大似然原理的互质阵列测向算法，虽然能有效增大阵列孔径，但仍然不适用于冲击噪声环境下的测向。
根据查阅已有文献发现，盘敏容等在《雷达科学与技术》(2020,Vol.18,NO.1,j.issn/1672-2337)上发表的“基于协方差矩阵重构的互质阵列DOA估计”中，所提出的波达估计方法在冲击噪声环境下性能较差。孙兵等在《系统工程与电子技术》(2020，ISSN 1001-506X,CN 11-2422/TN)上发表的“幅度相位误差条件下的互质阵列DOA估计方法”中，所提出的幅相误差校正方法不适用于冲击噪声环境，且在极低信噪比时性能较差。因此本发明提出一种在极低信噪比的冲击环境下有良好性能的幅相误差校正及测向方法，以解决传统的幅相误差矫正和测向方法在该环境下性能恶化的技术难题。
发明内容
针对现有冲击噪声环境下互质阵列的幅相误差校正和波达方向估计方法的缺点和不足，本发明设计了一种冲击噪声环境下基于量子哈里斯鹰算法的幅相误差校正和波达方向估计的方法，通过仿真可知，此方法在极低信噪比环境下相对于特殊阵列的幅相误差矫正和波达方向估计的传统方法有更好的鲁棒性。
本发明的目的是这样实现的：步骤如下：
步骤一：建立冲击噪声下存在幅相误差时互质阵列接收信号的数学模型，并采取有源校正的方法估计其幅相误差；
步骤二：初始化量子哈里斯鹰种群并构造目标函数和适应度函数；计算初始化后量子哈里斯鹰的适应度值并更新最优适应度和最优量子哈里斯鹰的位置；
步骤三：执行量子哈里斯鹰算法；
步骤四：检测循环次数t是否等于Ta，若是则输出最佳适应度
Figure BDA0003377598380000021
和最佳适应度位置
Figure BDA0003377598380000022
并按照公式
Figure BDA0003377598380000023
输出估计的第m个阵元的相位误差；若循环次数t小于Ta，则令t＝t+1，回到步骤三；
步骤五：用互质阵列接收信号，计算其分数低阶协方差矩阵，并进行幅相误差补偿；
步骤六：初始化量子哈里斯鹰种群并给出适应度函数；
步骤七：执行步骤三，在每一次迭代结束后，将第i只量子哈里斯鹰的位置向量
Figure BDA0003377598380000024
代入适应度函数中计算其适应度值；并更新最佳适应度
Figure BDA0003377598380000025
和最佳适应度的鹰所在的量子位置
Figure BDA0003377598380000026
检查迭代次数t是否等于设定最大迭代次数Tb，若是，则将
Figure BDA0003377598380000027
映射为估计的波达角
Figure BDA0003377598380000028
其映射关系为
Figure BDA0003377598380000029
若不是，则令t＝t+1，继续执行步骤三直至迭代次数达到要求为止。
进一步地，步骤一具体包括：计算互质阵列接收数据x(k)的分数低阶协方差矩阵Cfloc，该矩阵第g行第h列的元素为
Figure BDA00033775983800000210
其中(·)*是求共轭操作，E(·)是求数学期望操作，xg(t)、xh(t)分别是第g个阵元和第h个阵元接收数据，
Figure BDA00033775983800000211
p是任选计算常数，且0＜p≤1；分数低阶协方差矩阵的元素可由互质阵列接收到的有限快拍数据估计为
Figure BDA00033775983800000212
其中xg(k)、xh(k)是第g个、第h个阵元接收的第k次快拍数据；
对分数低阶协方差矩阵进行特征值分解，并将其特征值从大至小排序，其特征值对应的特征向量分别为ζ1至
Figure BDA00033775983800000213
由特征分解性质可知Γa(θ0)＝εζ1，其中
Figure BDA00033775983800000214
是分数低阶协方差矩阵最大的特征值对应的特征向量；该向量等式展开后可写为
Figure BDA0003377598380000031
等式中，
Figure BDA0003377598380000032
Figure BDA0003377598380000033
是第m个阵元的幅相误差值；等式左侧μm是待估计的值，右侧ε是未知参数，其余都是已知量；根据对应项相等的原则，可由εξ1,1＝1计算出参数ε的值，进一步即可估计出幅相误差矩阵记为Γb。
进一步地，步骤二具体包括：
首先，提取估计出的幅相误差矩阵Γb中阵元的相位值，记为
Figure BDA0003377598380000034
其中
Figure BDA0003377598380000035
是第m个阵元相位误差的估计值，
Figure BDA0003377598380000036
此时相位误差估计值的单位是弧度制，将其转化为角度制，转化关系为
Figure BDA0003377598380000037
随后初始化量子哈里斯鹰的量子位置并设定算法的参数；设定种群中，量子哈里斯鹰的数目为Na，最大迭代次数为Ta，迭代标号为t，t∈[1,Ta]；第t次迭代时，第i只量子哈里斯鹰的位置为
Figure BDA0003377598380000038
其中，
Figure BDA0003377598380000039
为阵元个数，且
Figure BDA00033775983800000310
初代量子哈里斯鹰的量子位置是[0,1]之间的均匀随机数；对于每一只量子哈里斯鹰，均设定其第m维的物理位置上限
Figure BDA00033775983800000311
其物理位置下限为
Figure BDA00033775983800000312
其中，δ是设定的量子哈里斯鹰搜索区间长度的一半；第t次迭代时，第i只量子哈里斯鹰第m维的量子位置与物理位置关系为
Figure BDA00033775983800000313
Figure BDA00033775983800000314
是第i只哈里斯鹰第t次迭代时，估计的第m个阵元的相位误差值；则第t次迭代时，由第i只量子哈里斯鹰的量子位置计算得出的相位误差估计对角度阵应当为
Figure BDA00033775983800000315
其中，
Figure BDA00033775983800000316
用矩阵Γg对步骤二中计算出的分数低阶协方差矩阵Cfloc进行误差补偿，记
Figure BDA00033775983800000317
上标H代表求矩阵的共轭转置，(·)-1是矩阵求逆的操作；然后生成原互质阵列的虚拟阵列，再由误差补偿后的分数低阶协方差矩阵
Figure BDA00033775983800000318
生成该虚拟阵列阵元接收的虚拟快拍矢量，并从中截取连续部分的虚拟阵元所接收的虚拟快拍矢量；
然后对误差补偿后的矩阵
Figure BDA00033775983800000319
进行向量化，记向量化后得到的
Figure BDA00033775983800000320
维的列向量为
Figure BDA00033775983800000321
同时对
Figure BDA00033775983800000322
η两个向量去冗余、排序，生成两个新的列向量
Figure BDA00033775983800000323
此时的向量
Figure BDA00033775983800000324
中含有连续的
Figure BDA0003377598380000041
行元素，该部分元素即是均匀虚拟阵列的位置集合，而
Figure BDA0003377598380000042
中对应位置的元素就是均匀虚拟阵列所接收的虚拟快拍数据，截取
Figure BDA0003377598380000043
中相应位置的元素，生成一个新的列向量
Figure BDA0003377598380000044
即是虚拟均匀阵元的单快拍数据矢量；
随后采用空间平滑算法，将虚拟阵列分为
Figure BDA0003377598380000045
个互相重叠的子阵列，每个子阵列有
Figure BDA0003377598380000046
个均匀阵元，将所有子阵列的虚拟接收数据矢量拼接成一个
Figure BDA0003377598380000047
维的接收矩阵
Figure BDA0003377598380000048
其中矩阵Z的第l列矢量Zl的元素对应的是虚拟均匀阵列接收数据矢量
Figure BDA0003377598380000049
中的第
Figure BDA00033775983800000410
行；
将Z视作一个
Figure BDA00033775983800000411
次快拍的接收信号矩阵，则其协方差矩阵RS可由矩阵Z的列向量估算为
Figure BDA00033775983800000412
矩阵RS会根据由量子哈里斯鹰算法所估计出的相位误差变化而变化；故可将RS记为
Figure BDA00033775983800000413
则其极大似然方程可写为
Figure BDA00033775983800000414
其中，tr(·)是矩阵求迹，向量
Figure BDA00033775983800000415
是第1个至第
Figure BDA00033775983800000416
个虚拟均匀线阵的导向矢量，
Figure BDA00033775983800000417
是第
Figure BDA00033775983800000418
个参与计算的虚拟均匀阵元的位置；根据极大似然方程写出适应度函数：
Figure BDA00033775983800000419
量子哈里斯鹰算法中定义变量
Figure BDA00033775983800000420
用于储存历史最佳的适应度值，定义向量
Figure BDA00033775983800000421
记录拥有最佳适应度值的鹰的量子位置；将所有初始化的量子哈里斯鹰的位置代入适应度函数中计算其适应度值，并将拥有最大适应度值的量子哈里斯鹰的适应度存入变量
Figure BDA00033775983800000422
中，将其量子位置存入
Figure BDA00033775983800000423
中。
进一步地，步骤三具体包括：量子哈里斯鹰对猎物的逃逸能量E进行评估，针对不同的逃逸能量，将选择不同的方法进行狩猎；猎物逃逸能量
Figure BDA00033775983800000424
其中E0是一个[-1,1]间的随机数；
(1)若|E|＞1,此时猎物的逃逸能量较大，第i只量子哈里斯鹰执行探索阶段，该阶段中第i只量子哈里斯鹰鹰第m维的量子旋转角的公式定义为：
Figure BDA0003377598380000051
其中，
Figure BDA0003377598380000052
是随机选择的一只量子哈里斯鹰的第m维量子坐标，
Figure BDA0003377598380000053
是当前所有哈里斯鹰第m维坐标的平均值，
Figure BDA0003377598380000054
分别代表量子哈里斯鹰第m维坐标的最大值和最小值；
Figure BDA0003377598380000055
和
Figure BDA0003377598380000056
都是取值在[0,1]之间的随机数，
(2)若|E|＜1，此时猎物逃逸能量小，则量子鹰哈里斯鹰i执行开发阶段，该阶段中算法生成随机数
Figure BDA0003377598380000057
并和猎物逃逸能量共同进行判断以选择四种不同的开发策略：
策略一：软围攻，当|E|≥0.5且
Figure BDA0003377598380000058
时第i只量子哈里斯鹰执行此策略；
执行软围攻时，第i只量子哈里斯鹰第m维的量子旋转角公式应为
Figure BDA0003377598380000059
其中，
Figure BDA00033775983800000510
Figure BDA00033775983800000511
为[0,1]间的随机数；
策略二：渐进式软围攻，当|E|≥0.5且
Figure BDA00033775983800000512
时第i只量子哈里斯鹰执行此策略；
此时，第i只量子鹰有两种不同的方法更新其量子旋转角，首先采取方法A1，此时的量子旋转角为
Figure BDA00033775983800000513
根据此量子旋转角改变该量子鹰的量子位置并计算其适应度，如果该值小于迭代次数为t时该量子鹰的适应度值，则转而执行方法B1，方法B1的量子旋转角公式为
Figure BDA00033775983800000514
其中
Figure BDA00033775983800000515
是[0,1]之间的随机数，F为归一化后的levy飞行函数；
策略三：硬包围，当|E|＜0.5且
Figure BDA00033775983800000516
时，量子哈里斯鹰i将执行硬包围策略；
此时的量子旋转角为
Figure BDA00033775983800000517
策略四：渐进式硬包围，当|E|＜0.5且
Figure BDA00033775983800000518
时，量子鹰i将执行渐进式硬包围；
此时，该量子鹰会采取两种不同的方法进行量子旋转角的更新；首先采用方法A2，方法A2的量子旋转角为
Figure BDA00033775983800000519
用此量子旋转角更新该量子鹰的量子位置并代入适应度函数计算其适应度，如果该值小于迭代次数为t时该量子鹰的适应度值，则转而执行方法B2；方法B2中第i只量子哈里斯鹰第m维的量子旋转角计算公式为
Figure BDA00033775983800000520
用以上策略计算得出的量子旋转角更新量子哈里斯鹰的量子位置，其量子位置更新公式为
Figure BDA0003377598380000061
根据更新后的量子位置
Figure BDA0003377598380000062
带入适应度函数中计算
Figure BDA0003377598380000063
和更新
Figure BDA0003377598380000064
进一步地，步骤五具体包括：将步骤四中最终输出的最佳量子鹰的位置映射为相位误差估计值
Figure BDA0003377598380000065
至
Figure BDA0003377598380000066
从而得出相位误差估计矩阵
Figure BDA0003377598380000067
同时提取步骤二中粗估计的幅相误差阵的各元素的模值，记作abs(Γb)，其中abs(·)表示求括号内矩阵元素的模值的操作；生成最终的幅相误差估计矩阵为Cg＝Γg·abs(Γb)；
互质阵列的接收信号矢量可表示为
Figure BDA0003377598380000068
其中
Figure BDA0003377598380000069
是Q×1维的信号的第o次快拍数据矢量，o∈[1,O]，Q为信号源的个数；
Figure BDA00033775983800000610
是互质阵列在第o次快拍下所接收到的冲击噪声矢量；
Figure BDA00033775983800000611
是Q个信号源的导向矢量阵，其中
Figure BDA00033775983800000612
是第q个信源对应的导向矢量，且q∈[1,Q]；设
Figure BDA00033775983800000613
是接收信号的分数低阶协方差矩阵，其中各个位置的元素可根据互质阵列的接收数据进行估算，具体的估算公式为
Figure BDA00033775983800000614
其中
Figure BDA00033775983800000615
分别是第g个、第h个阵元接收到的第o次快拍数据；随后对该共变矩阵进行幅相误差补偿，幅相误差补偿后的分数低阶协方差矩阵应为
Figure BDA00033775983800000616
进一步地，步骤六具体包括：依据步骤二中的方法，将幅相误差补偿后的分数低阶协方差矩阵C0进行向量化操作，截取其中虚拟均匀阵元的虚拟接收数据，组成一个
Figure BDA00033775983800000617
维的虚拟接收向量
Figure BDA00033775983800000618
再对其按照步骤二中的方法进行空间平滑算法，可以得到一个
Figure BDA00033775983800000619
维的矩阵
Figure BDA00033775983800000620
该虚拟快拍数据矩阵的协方差矩阵可估算为
Figure BDA00033775983800000621
这里的
Figure BDA00033775983800000622
是矩阵
Figure BDA00033775983800000623
的第l维列向量；
令量子哈里斯鹰的种群规模为Nb，最大迭代次数为Tb，第t次迭代时第i只哈里斯鹰的量子位置设定为
Figure BDA00033775983800000624
初次迭代时
Figure BDA00033775983800000625
是一个[0,1]间的随机数，Q为信源个数，量子鹰的每一维估计一个信源的到达角，将第i只鹰的第q维转化为波达角为
Figure BDA0003377598380000071
其中，
Figure BDA0003377598380000072
是量子哈里斯鹰算法的搜索下限和上限，由第i只量子哈里斯鹰算法计算出的信源导向矢量为
Figure BDA0003377598380000073
其中，
Figure BDA0003377598380000074
根据求出的虚拟快拍数据的协方差矩阵Rss，可写出极大似然方程为
Figure BDA0003377598380000075
因此定义第t次迭代时第i只量子哈里斯鹰适应度函数为
Figure BDA0003377598380000076
定义第t次迭代时，量子哈里斯鹰种群的最佳适应度为
Figure BDA0003377598380000077
最佳适应度的鹰所在的量子位置为
Figure BDA0003377598380000078
与现有技术相比，本发明的有益效果是：与现有的互质阵列在冲击噪声下的幅相误差校正技术相比，本发明针对传统的有源幅相误差估计方法在极低信噪比下性能恶化的问题，设计出更具鲁棒性的幅相误差校正方法。首先采用外加辅助源的校正算法得出幅相误差的粗估计值，再用量子哈里斯鹰算法在粗估计相位误差周围进行搜索，可以实现在极低信噪比下对幅相误差进行更精确的估计。同时，在冲击噪声下互质阵列的波达方向估计问题上，本发明所设计的基于量子哈里斯鹰机制的分数低阶协方差结合虚拟矩阵的极大似然测向方法，可在相同信噪比下取得比其他传统算法更低的均方根误差，其中引入的虚拟阵列和空间平滑算法，可有效提高互质阵列的空间自由度。
附图说明
图1是本发明所设计的基于量子哈里斯鹰机制的幅相误差校正和测向方法示意图。
图2是量子哈里斯鹰算法的具体步骤示意图。
图3是互质阵列示意图。
图4是在α＝1.25时，三种不同的方法下所估计相位误差的均方根误差与辅助源广义信噪比的关系曲线。
图5是采取本方法校正幅相误差前后，对其进行基于极大似然方程的谱峰搜索所得到的谱峰图。
图6是在α＝1.25下，分别用三种不同方法进行幅相误差校正和信号源测向的均方根误差与辅助源广义信噪比的关系曲线。
图7是在α＝1.25下，采取本文所设计的幅相误差校正方法后，用三种不同方法进行测向的均方根误差与信源广义信噪比的关系曲线。
具体实施方式
下面结合附图与具体实施方式对本发明作进一步详细描述。
结合图1-图7，本发明的步骤如下：
步骤一，建立冲击噪声下存在幅相误差时互质阵列接收信号的数学模型，设置辅助源，采取有源校正的方法估计互质阵列的幅相误差。
构建的互质阵列由A,B两个均匀子阵列组成。假定用于校正的信号源信号波长为λ，设两个阵元的单位间距
Figure BDA0003377598380000081
为信号波长的一半，阵列A含有
Figure BDA0003377598380000082
个阵元，阵元间距为
Figure BDA0003377598380000083
阵列A的阵元位置是
Figure BDA0003377598380000084
阵列B含
Figure BDA0003377598380000085
有个阵元，阵元间距为
Figure BDA0003377598380000086
则其阵元位置为
Figure BDA0003377598380000087
其中，
Figure BDA0003377598380000088
且二数互质。
将A、B两个均匀阵列组合在一起，由于
Figure BDA0003377598380000089
互质，除了第一个阵元外其它阵元均不重合。最后组成的互质阵列共由
Figure BDA00033775983800000810
个阵元组成。
假定用于校正的外部辅助源信号来波方向与阵列法线夹角为θ0，信号和噪声均服从同向复对称的SαS分布且不相关。互质阵列接受的第k次快拍数据矢量可以表示为x(k)＝Γa0(θ0)s(k)+N(k)。式中，
Figure BDA00033775983800000811
是互质阵列接收到的第k次快拍数据矢量，上标T是矩阵求转置的操作。s(k)是辅助信号源发出的第k次快拍数据，N(k)是
Figure BDA00033775983800000812
维服从复冲击噪声分布矢量，设互质阵列一共接收到L次快拍，则有k∈[1,L]。
Figure BDA00033775983800000813
是互质阵列的幅相误差对角矩阵，ρm、
Figure BDA00033775983800000814
为第m个阵元的幅度误差和相位误差，
Figure BDA00033775983800000815
以首阵元为参考阵元，并对其幅相误差做归一化处理即令ρ1＝1,
Figure BDA00033775983800000816
当互质阵列不存在幅相误差时，
Figure BDA00033775983800000817
维的导向矢量阵为
Figure BDA00033775983800000818
其中dm是第m个互质阵元相对于首个阵元的坐标位置。
计算互质阵列接收数据x(k)的分数低阶协方差矩阵Cfloc，该矩阵第g行第h列的元素为
Figure BDA00033775983800000819
其中(·)*是求共轭操作，E(·)是求数学期望操作，xg(t)、xh(t)分别是第g个阵元和第h个阵元接收数据，
Figure BDA00033775983800000820
p是任选计算常数，且0＜p≤1。实际计算中，分数低阶协方差矩阵的元素可由互质阵列接收到的有限快拍数据估计为
Figure BDA0003377598380000091
其中xg(k)、xh(k)是第g个、第h个阵元接收的第k次快拍数据。
对分数低阶协方差矩阵进行特征值分解，并将其特征值从大至小排序，其特征值对应的特征向量分别为ζ1至
Figure BDA0003377598380000092
由特征分解性质可知Γa(θ0)＝εζ1，其中
Figure BDA0003377598380000093
是分数低阶协方差矩阵最大的特征值对应的特征向量。该向量等式展开后可写为
Figure BDA0003377598380000094
等式中，
Figure BDA0003377598380000095
Figure BDA0003377598380000096
是第m个阵元的幅相误差值。等式左侧μm是待估计的值，右侧ε是未知参数，其余都是已知量。根据对应项相等的原则，可由εξ1,1＝1计算出参数ε的值，进一步即可估计出幅相误差矩阵记为Γb。由于冲击噪声的影响，这一步中估算出的幅相误差值会存在误差。
步骤二，通过极大似然方程给出目标函数和适应度函数，初始化量子哈里斯鹰量子位置并设定参数。
由于相位误差在低信噪比的冲击噪声环境下估计精度较低，因此采取量子哈里斯鹰算法对相位误差进行进一步估计。首先，提取估计出的幅相误差矩阵Γb中阵元的相位值，记为
Figure BDA0003377598380000097
其中
Figure BDA0003377598380000098
是第m个阵元相位误差的估计值，
Figure BDA0003377598380000099
此时相位误差估计值的单位是弧度制，将其转化为角度制，转化关系为
Figure BDA00033775983800000910
随后初始化量子哈里斯鹰的量子位置并设定算法的参数。设定种群中，量子哈里斯鹰的数目为Na，最大迭代次数为Ta，迭代标号为t，t∈[1,Ta]。第t次迭代时，第i只量子哈里斯鹰的位置为
Figure BDA00033775983800000911
其中，
Figure BDA00033775983800000912
为阵元个数，且
Figure BDA00033775983800000913
初代量子哈里斯鹰的量子位置是[0,1]之间的均匀随机数。对于每一只量子哈里斯鹰，均设定其第m维的物理位置上限
Figure BDA00033775983800000914
其物理位置下限为
Figure BDA00033775983800000915
其中，δ是设定的量子哈里斯鹰搜索区间长度的一半。第t次迭代时，第i只量子哈里斯鹰第m维的量子位置与物理位置关系为
Figure BDA00033775983800000916
Figure BDA00033775983800000917
是第i只哈里斯鹰第t次迭代时，估计的第m个阵元的相位误差值。则第t次迭代时，由第i只量子哈里斯鹰的量子位置计算得出的相位误差估计对角度阵应当为
Figure BDA00033775983800000918
其中，
Figure BDA0003377598380000101
用矩阵Γg对步骤二中计算出的分数低阶协方差矩阵Cfloc进行误差补偿，记
Figure BDA0003377598380000102
这里上标H代表求矩阵的共轭转置，(·)-1是矩阵求逆的操作。然后生成原互质阵列的虚拟阵列，再由误差补偿后的分数低阶协方差矩阵
Figure BDA0003377598380000103
生成该虚拟阵列阵元接收的虚拟快拍矢量，并从中截取连续部分的虚拟阵元所接收的虚拟快拍矢量。具体方法是，定义一个
Figure BDA0003377598380000104
维向量η，
Figure BDA0003377598380000105
其中，ηm＝dm-ν，是一个
Figure BDA0003377598380000106
维的行向量，其中dm是互质阵列中第m个阵元的位置，
Figure BDA0003377598380000107
ν是
Figure BDA0003377598380000108
维的原互质阵列所有阵元的位置矢量，向量η即是虚拟阵列阵元的位置矢量。然后对误差补偿后的矩阵
Figure BDA0003377598380000109
进行向量化，记向量化后得到的
Figure BDA00033775983800001010
维的列向量为
Figure BDA00033775983800001011
同时对
Figure BDA00033775983800001012
η两个向量去冗余、排序，生成两个新的列向量
Figure BDA00033775983800001013
此时的向量
Figure BDA00033775983800001014
中含有连续的
Figure BDA00033775983800001015
行元素，该部分元素即是均匀虚拟阵列的位置集合，而
Figure BDA00033775983800001016
中对应位置的元素就是均匀虚拟阵列所接收的虚拟快拍数据，截取
Figure BDA00033775983800001017
中相应位置的元素，生成一个新的列向量
Figure BDA00033775983800001018
即是虚拟均匀阵元的单快拍数据矢量。
随后采用空间平滑算法，将虚拟阵列分为
Figure BDA00033775983800001019
个互相重叠的子阵列，每个子阵列有
Figure BDA00033775983800001020
个均匀阵元，将所有子阵列的虚拟接收数据矢量拼接成一个
Figure BDA00033775983800001021
维的接收矩阵
Figure BDA00033775983800001022
其中矩阵Z的第l列矢量Zl的元素对应的是虚拟均匀阵列接收数据矢量
Figure BDA00033775983800001023
中的第
Figure BDA00033775983800001024
行。
将Z视作一个
Figure BDA00033775983800001025
次快拍的接收信号矩阵，则其协方差矩阵RS可由矩阵Z的列向量估算为
Figure BDA00033775983800001026
矩阵RS会根据由量子哈里斯鹰算法所估计出的相位误差变化而变化。故可将RS记为
Figure BDA00033775983800001027
则其极大似然方程可写为
Figure BDA00033775983800001028
其中，tr(·)是矩阵求迹，向量
Figure BDA00033775983800001029
是第1个至第
Figure BDA00033775983800001030
个虚拟均匀线阵的导向矢量，
Figure BDA00033775983800001031
是第
Figure BDA00033775983800001032
个参与计算的虚拟均匀阵元的位置。可根据极大似然方程写出适应度函数：
Figure BDA0003377598380000111
量子哈里斯鹰算法中定义变量
Figure BDA0003377598380000112
用于储存历史最佳的适应度值，定义向量
Figure BDA0003377598380000113
记录拥有最佳适应度值的鹰的量子位置。将所有初始化的量子哈里斯鹰的位置代入适应度函数中计算其适应度值，并将拥有最大适应度值的量子哈里斯鹰的适应度存入变量
Figure BDA0003377598380000114
中，将其量子位置存入
Figure BDA0003377598380000115
中。
步骤三，执行量子哈里斯鹰算法。
量子哈里斯鹰对猎物的逃逸能量E进行评估，针对不同的逃逸能量，将选择几种不同的方法进行狩猎。猎物逃逸能量
Figure BDA0003377598380000116
其中E0是一个[-1,1]间的随机数。
(1)若|E|＞1,此时猎物的逃逸能量较大，第i只量子哈里斯鹰执行探索阶段
该阶段中第i只量子哈里斯鹰鹰第m维的量子旋转角的公式定义为
Figure BDA0003377598380000117
其中，
Figure BDA0003377598380000118
是随机选择的一只量子哈里斯鹰的第m维量子坐标。
Figure BDA0003377598380000119
是当前所有哈里斯鹰第m维坐标的平均值。
Figure BDA00033775983800001110
分别代表量子哈里斯鹰第m维坐标的最大值和最小值。
Figure BDA00033775983800001111
和
Figure BDA00033775983800001112
都是取值在[0,1]之间的随机数。
(2)若|E|＜1，此时猎物逃逸能量小，则量子鹰哈里斯鹰i执行开发阶段
该阶段中算法生成随机数
Figure BDA00033775983800001113
并和猎物逃逸能量共同进行判断以选择四种不同的开发策略。
策略一：软围攻，当|E|≥0.5且
Figure BDA00033775983800001114
时第i只量子哈里斯鹰执行此策略。
执行软围攻时，第i只量子哈里斯鹰第m维的量子旋转角公式应为
Figure BDA00033775983800001115
其中，
Figure BDA00033775983800001116
Figure BDA00033775983800001117
为[0,1]间的随机数。
策略二：渐进式软围攻，当|E|≥0.5且
Figure BDA00033775983800001118
时第i只量子哈里斯鹰执行此策略。
此时，第i只量子鹰有两种不同的方法更新其量子旋转角。首先采取方法A1，此时的量子旋转角为
Figure BDA00033775983800001119
根据此量子旋转角改变该量子鹰的量子位置并计算其适应度，如果该值小于迭代次数为t时该量子鹰的适应度值，则转而执行方法B1。方法B1的量子旋转角公式为
Figure BDA0003377598380000121
其中
Figure BDA0003377598380000122
是[0,1]之间的随机数，F为归一化后的levy飞行函数。
策略三：硬包围，当|E|＜0.5且
Figure BDA0003377598380000123
时，量子哈里斯鹰i将执行硬包围策略。
此时的量子旋转角为
Figure BDA0003377598380000124
策略四：渐进式硬包围，当|E|＜0.5且
Figure BDA0003377598380000125
时，量子鹰i将执行渐进式硬包围。
此时，该量子鹰会采取两种不同的方法进行量子旋转角的更新。首先采用方法A2，方法A2的量子旋转角为
Figure BDA0003377598380000126
用此量子旋转角更新该量子鹰的量子位置并代入适应度函数计算其适应度，如果该值小于迭代次数为t时该量子鹰的适应度值，则转而执行方法B2。方法B2中第i只量子哈里斯鹰第m维的量子旋转角计算公式为
Figure BDA0003377598380000127
用以上策略计算得出的量子旋转角更新量子哈里斯鹰的量子位置，其量子位置更新公式为
Figure BDA0003377598380000128
根据更新后的量子位置
Figure BDA0003377598380000129
带入适应度函数中计算
Figure BDA00033775983800001210
和更新
Figure BDA00033775983800001211
步骤四，检测循环次数t是否等于Ta，若是则输出最佳适应度
Figure BDA00033775983800001212
和最佳适应度位置
Figure BDA00033775983800001213
并按照公式
Figure BDA00033775983800001214
输出估计的第m个阵元的相位误差。若循环次数t小于Ta，则令t＝t+1，回到步骤三。
步骤五，建立冲击噪声下互质阵列接收信号的模型，计算其分数低阶协方差并对其进行幅相误差矫正
将步骤四中最终输出的最佳量子鹰的位置映射为相位误差估计值
Figure BDA00033775983800001215
至
Figure BDA00033775983800001216
从而得出相位误差估计矩阵即
Figure BDA00033775983800001217
同时提取步骤二中粗估计的幅相误差阵的各元素的模值，记作abs(Γb)，其中abs(·)表示求括号内矩阵元素的模值的操作。生成最终的幅相误差估计矩阵为Cg＝Γg·abs(Γb)。
互质阵列的接收信号矢量可表示为
Figure BDA00033775983800001218
其中
Figure BDA00033775983800001219
是Q×1维的信号的第o次快拍数据矢量，o∈[1,O]，Q为信号源的个数。
Figure BDA00033775983800001220
是互质阵列在第o次快拍下所接收到的冲击噪声矢量。
Figure BDA0003377598380000131
是Q个信号源的导向矢量阵，其中
Figure BDA0003377598380000132
是第q个信源对应的导向矢量，且q∈[1,Q]。设
Figure BDA0003377598380000133
是接收信号的分数低阶协方差矩阵，其中各个位置的元素可根据互质阵列的接收数据进行估算，具体的估算公式为
Figure BDA0003377598380000134
其中
Figure BDA0003377598380000135
分别是第g个、第h个阵元接收到的第o次快拍数据。随后对该共变矩阵进行幅相误差补偿，幅相误差补偿后的分数低阶协方差矩阵应为
Figure BDA0003377598380000136
步骤六，由极大似然方程给出目标函数和适应度函数，并初始化量子哈里斯鹰种群，设定算法参数。
依据步骤二中的方法，将幅相误差补偿后的分数低阶协方差矩阵C0进行向量化操作，截取其中虚拟均匀阵元的虚拟接收数据，组成一个
Figure BDA0003377598380000137
维的虚拟接收向量
Figure BDA0003377598380000138
再对其按照步骤二中的方法进行空间平滑算法，可以得到一个
Figure BDA0003377598380000139
维的矩阵
Figure BDA00033775983800001310
该虚拟快拍数据矩阵的协方差矩阵可估算为
Figure BDA00033775983800001311
这里的
Figure BDA00033775983800001312
是矩阵
Figure BDA00033775983800001313
的第l维列向量。
令量子哈里斯鹰的种群规模为Nb，最大迭代次数为Tb。第t次迭代时第i只哈里斯鹰的量子位置设定为
Figure BDA00033775983800001314
初次迭代时
Figure BDA00033775983800001315
是一个[0,1]间的随机数，Q为信源个数，量子鹰的每一维估计一个信源的到达角，将第i只鹰的第q维转化为波达角为
Figure BDA00033775983800001316
其中，
Figure BDA00033775983800001317
是量子哈里斯鹰算法的搜索下限和上限，由第i只量子哈里斯鹰算法计算出的信源导向矢量为
Figure BDA00033775983800001318
其中，
Figure BDA00033775983800001319
根据求出的虚拟快拍数据的协方差矩阵Rss，可写出极大似然方程为
Figure BDA00033775983800001320
因此定义第t次迭代时第i只量子哈里斯鹰适应度函数为
Figure BDA00033775983800001321
定义第t次迭代时，量子哈里斯鹰种群的最佳适应度为
Figure BDA00033775983800001322
最佳适应度的鹰所在的量子位置为
Figure BDA00033775983800001323
步骤七，执行步骤三，在每一次迭代结束后，将第i只量子哈里斯鹰的位置向量
Figure BDA0003377598380000141
代入适应度函数中计算其适应度值。并更新最佳适应度
Figure BDA0003377598380000142
和最佳适应度的鹰所在的量子位置
Figure BDA0003377598380000143
检查迭代次数t是否等于设定最大迭代次数Tb，若是，则将
Figure BDA0003377598380000144
映射为估计的波达角
Figure BDA0003377598380000145
其映射关系为
Figure BDA0003377598380000146
若不是，则令t＝t+1，继续执行步骤三直至迭代次数达到要求为止。
图4中，本发明所设计的基于量子哈里斯鹰机制的分数低阶协方差单辅助源误差校正方法记作基于FLOC-QHHO的误差校正；基于分数低阶协方差的单辅助源误差校正方法记作基于FLOC的误差校正；基于共变矩阵的单辅助源误差校正方法记作基于ROC的误差校正。本例中，量子哈里斯鹰算法的个体数Na为100，迭代次数Ta为200次，搜索区间半径设为δ＝10。设定组成互质阵列的阵元A有3个均匀阵元，组成互质阵列的阵元B有5个均匀阵元。互质阵列共有10个阵元组成。阵元幅度误差为0.85至1.15之间的随机数，相位误差为-60°至60°间的随机数，分数低阶协方差的计算参数p＝0.6，共变矩阵的计算参数p1＝1.2，蒙特卡洛实验次数为300。可见，相比另外两种传统的单辅助源幅相误差校正方法，本发明在低信噪比下所估计的相位误差值均方根误差更小，更具鲁棒性。
图5中辅助源的广义信噪比设为1dB，信号源的广义信噪比均设置为0dB，两个信号源的角度分别为-10°和50°，可见采用本方法后，原本畸变的谱峰图得到校正，可使后续的测向算法能达到更高的精确度。
在图6中，本发明中设计的幅相误差校正方法记作FLOC-QHHO误差校正，设计的基于量子哈里斯鹰机制的分数低阶协方差结合虚拟阵列的极大似然测向方法记作FLOC-QHHO-SS-ML；第二种对比方法中，采取FLOC-QHHO误差校正方法进行幅相误差校正，并采取基于量子哈里斯鹰机制的共变矩阵结合虚拟阵列的极大似然测向方法进行测向，该测向方法记作ROC-QHHO-SS-ML；第三种对比方法中，采取基于FLOC校正方法进行误差校正并采用FLOC-QHHO-SS-ML方法进行测向。本例中，首次执行量子哈里斯鹰算法的个体数Na为100，迭代次数Ta为200次，搜索区间半径设为δ＝10。第二次执行量子哈里斯鹰算法时，个体数Nb为100，迭代次数Tb为200次。设定组成互质阵列的阵元A有3个均匀阵元，组成互质阵列的阵元B有5个均匀阵元。互质阵列共有10个阵元组成。阵元幅度误差为0.85至1.15之间的随机数，相位误差为-60°至60°间的随机数，分数低阶协方差的计算参数p＝0.6，共变矩阵的计算参数p1＝1.2，两个信号源的角度分别是-10°和10°，其广义信噪比为1dB。蒙特卡洛实验次数为300，从-8dB至5dB改变校正源的广义信噪比。从以上三种方法的均方根误差变化曲线可以看出，本发明中设计的幅相误差矫正和测向方法在恶劣信噪比下表现更好，相比其他方法更具鲁棒性。理论上，FLOC-QHHO-SS-ML测向方法通常比ROC-QHHO-SS-ML测向方法性能更好，然而采取本发明所设计的FLOC-QHHO误差校正后，FLOC-QHHO-SS-ML在校正源的信噪比较低的情况下其均方根误差低于基于FLOC的传统误差校正方法进行误差校正后再进行FLOC-QHHO-SS-ML测向的均方根误差，充分体现了本发明设计的幅相误差校正方法的鲁棒性。
在图7中，本发明所设计的方法记作FLOC-SS-QHHO-ML，基于量子哈里斯鹰机制的分数低阶矩结合虚拟阵列的极大似然测向方法记作FLOM-SS-QHHO-ML，基于量子哈里斯鹰机制的共变矩阵结合虚拟阵列的极大似然测向方法记作ROC-SS-QHHO-ML，本例中，两次量子哈里斯鹰算法中的参数、幅相误差的范围、互质阵列的设置、共变矩阵和分数低阶协方差的计算参数均与图6中的实验相同，并设定用于校正的辅助源广义信噪比为0dB。两个待估计的信号源的角度角度分别是-10°和10°，从-15dB到10dB依次改变信源的广义信噪比，进行150次蒙特卡洛实验后，从最终绘制的曲线中可以看出，本发明中设计的测向方法在信源的广义信噪比较小时，仍能保持较低的均方根误差，相比其它两种测向方法更具鲁棒性。

Claims (6)

1.冲击噪声环境下互质阵列的幅相误差校正和测向方法，其特征在于，步骤如下：
步骤一：建立冲击噪声下存在幅相误差时互质阵列接收信号的数学模型，并采取有源校正的方法估计其幅相误差；
步骤二：初始化量子哈里斯鹰种群并构造目标函数和适应度函数；计算初始化后量子哈里斯鹰的适应度值并更新最优适应度和最优量子哈里斯鹰的位置；
步骤三：执行量子哈里斯鹰算法；
步骤四：检测循环次数t是否等于Ta，若是则输出最佳适应度
Figure FDA0003377598370000011
和最佳适应度位置
Figure FDA0003377598370000012
并按照公式
Figure FDA0003377598370000013
输出估计的第m个阵元的相位误差；若循环次数t小于Ta，则令t＝t+1，回到步骤三；
步骤五：用互质阵列接收信号，计算其分数低阶协方差矩阵，并进行幅相误差补偿；
步骤六：初始化量子哈里斯鹰种群并给出适应度函数；
步骤七：执行步骤三，在每一次迭代结束后，将第i只量子哈里斯鹰的位置向量
Figure FDA0003377598370000014
代入适应度函数中计算其适应度值；并更新最佳适应度
Figure FDA0003377598370000015
和最佳适应度的鹰所在的量子位置
Figure FDA0003377598370000016
检查迭代次数t是否等于设定最大迭代次数Tb，若是，则将
Figure FDA0003377598370000017
映射为估计的波达角
Figure FDA0003377598370000018
其映射关系为
Figure FDA0003377598370000019
若不是，则令t＝t+1，继续执行步骤三直至迭代次数达到要求为止。
2.根据权利要求1所述的冲击噪声环境下互质阵列的幅相误差校正和测向方法，其特征在于，步骤一具体包括：计算互质阵列接收数据x(k)的分数低阶协方差矩阵Cfloc，该矩阵第g行第h列的元素为
Figure FDA00033775983700000110
其中(·)*是求共轭操作，E(·)是求数学期望操作，xg(t)、xh(t)分别是第g个阵元和第h个阵元接收数据，
Figure FDA00033775983700000111
p是任选计算常数，且0＜p≤1；分数低阶协方差矩阵的元素可由互质阵列接收到的有限快拍数据估计为
Figure FDA00033775983700000112
其中xg(k)、xh(k)是第g个、第h个阵元接收的第k次快拍数据；
对分数低阶协方差矩阵进行特征值分解，并将其特征值从大至小排序，其特征值对应的特征向量分别为ζ1至
Figure FDA00033775983700000113
由特征分解性质可知Γa(θ0)＝εζ1，其中
Figure FDA00033775983700000114
是分数低阶协方差矩阵最大的特征值对应的特征向量；该向量等式展开后可写为
Figure FDA0003377598370000021
等式中，
Figure FDA0003377598370000022
Figure FDA0003377598370000023
是第m个阵元的幅相误差值；等式左侧μm是待估计的值，右侧ε是未知参数，其余都是已知量；根据对应项相等的原则，可由εξ1,1＝1计算出参数ε的值，进一步即可估计出幅相误差矩阵记为Γb。
3.根据权利要求1所述的冲击噪声环境下互质阵列的幅相误差校正和测向方法，其特征在于，步骤二具体包括：
首先，提取估计出的幅相误差矩阵Γb中阵元的相位值，记为
Figure FDA0003377598370000024
其中
Figure FDA0003377598370000025
是第m个阵元相位误差的估计值，
Figure FDA0003377598370000026
此时相位误差估计值的单位是弧度制，将其转化为角度制，转化关系为
Figure FDA0003377598370000027
随后初始化量子哈里斯鹰的量子位置并设定算法的参数；设定种群中，量子哈里斯鹰的数目为Na，最大迭代次数为Ta，迭代标号为t，t∈[1,Ta]；第t次迭代时，第i只量子哈里斯鹰的位置为
Figure FDA0003377598370000028
其中，
Figure FDA0003377598370000029
为阵元个数，且
Figure FDA00033775983700000210
初代量子哈里斯鹰的量子位置是[0,1]之间的均匀随机数；对于每一只量子哈里斯鹰，均设定其第m维的物理位置上限
Figure FDA00033775983700000211
其物理位置下限为
Figure FDA00033775983700000212
其中，δ是设定的量子哈里斯鹰搜索区间长度的一半；第t次迭代时，第i只量子哈里斯鹰第m维的量子位置与物理位置关系为
Figure FDA00033775983700000213
Figure FDA00033775983700000214
是第i只哈里斯鹰第t次迭代时，估计的第m个阵元的相位误差值；则第t次迭代时，由第i只量子哈里斯鹰的量子位置计算得出的相位误差估计对角度阵应当为
Figure FDA00033775983700000215
其中，
Figure FDA00033775983700000216
用矩阵Γg对步骤二中计算出的分数低阶协方差矩阵Cfloc进行误差补偿，记
Figure FDA00033775983700000217
上标H代表求矩阵的共轭转置，(·)-1是矩阵求逆的操作；然后生成原互质阵列的虚拟阵列，再由误差补偿后的分数低阶协方差矩阵
Figure FDA00033775983700000218
生成该虚拟阵列阵元接收的虚拟快拍矢量，并从中截取连续部分的虚拟阵元所接收的虚拟快拍矢量；
然后对误差补偿后的矩阵
Figure FDA00033775983700000219
进行向量化，记向量化后得到的
Figure FDA00033775983700000220
维的列向量为
Figure FDA0003377598370000031
同时对
Figure FDA0003377598370000032
η两个向量去冗余、排序，生成两个新的列向量
Figure FDA0003377598370000033
此时的向量
Figure FDA0003377598370000034
中含有连续的
Figure FDA0003377598370000035
行元素，该部分元素即是均匀虚拟阵列的位置集合，而
Figure FDA0003377598370000036
中对应位置的元素就是均匀虚拟阵列所接收的虚拟快拍数据，截取
Figure FDA0003377598370000037
中相应位置的元素，生成一个新的列向量
Figure FDA0003377598370000038
即是虚拟均匀阵元的单快拍数据矢量；
随后采用空间平滑算法，将虚拟阵列分为
Figure FDA0003377598370000039
个互相重叠的子阵列，每个子阵列有
Figure FDA00033775983700000310
个均匀阵元，将所有子阵列的虚拟接收数据矢量拼接成一个
Figure FDA00033775983700000311
维的接收矩阵
Figure FDA00033775983700000312
其中矩阵Z的第l列矢量Zl的元素对应的是虚拟均匀阵列接收数据矢量
Figure FDA00033775983700000313
中的第
Figure FDA00033775983700000314
行；
将Z视作一个
Figure FDA00033775983700000315
次快拍的接收信号矩阵，则其协方差矩阵RS可由矩阵Z的列向量估算为
Figure FDA00033775983700000316
矩阵RS会根据由量子哈里斯鹰算法所估计出的相位误差变化而变化；故可将RS记为
Figure FDA00033775983700000317
则其极大似然方程可写为
Figure FDA00033775983700000318
其中，tr(·)是矩阵求迹，向量
Figure FDA00033775983700000319
是第1个至第
Figure FDA00033775983700000320
个虚拟均匀线阵的导向矢量，
Figure FDA00033775983700000321
是第
Figure FDA00033775983700000322
个参与计算的虚拟均匀阵元的位置；根据极大似然方程写出适应度函数：
Figure FDA00033775983700000323
量子哈里斯鹰算法中定义变量
Figure FDA00033775983700000328
用于储存历史最佳的适应度值，定义向量
Figure FDA00033775983700000324
记录拥有最佳适应度值的鹰的量子位置；将所有初始化的量子哈里斯鹰的位置代入适应度函数中计算其适应度值，并将拥有最大适应度值的量子哈里斯鹰的适应度存入变量
Figure FDA00033775983700000325
中，将其量子位置存入
Figure FDA00033775983700000326
中。
4.根据权利要求1所述的冲击噪声环境下互质阵列的幅相误差校正和测向方法，其特征在于，步骤三具体包括：量子哈里斯鹰对猎物的逃逸能量E进行评估，针对不同的逃逸能量，将选择不同的方法进行狩猎；猎物逃逸能量
Figure FDA00033775983700000327
其中E0是一个[-1,1]间的随机数；
(1)若|E|＞1,此时猎物的逃逸能量较大，第i只量子哈里斯鹰执行探索阶段，该阶段中第i只量子哈里斯鹰鹰第m维的量子旋转角的公式定义为：
Figure FDA0003377598370000041
其中，
Figure FDA0003377598370000042
是随机选择的一只量子哈里斯鹰的第m维量子坐标，
Figure FDA0003377598370000043
是当前所有哈里斯鹰第m维坐标的平均值，
Figure FDA0003377598370000044
分别代表量子哈里斯鹰第m维坐标的最大值和最小值；
Figure FDA0003377598370000045
和
Figure FDA0003377598370000046
都是取值在[0,1]之间的随机数，
(2)若|E|＜1，此时猎物逃逸能量小，则量子鹰哈里斯鹰i执行开发阶段，该阶段中算法生成随机数
Figure FDA0003377598370000047
并和猎物逃逸能量共同进行判断以选择四种不同的开发策略：
策略一：软围攻，当|E|≥0.5且
Figure FDA0003377598370000048
时第i只量子哈里斯鹰执行此策略；
执行软围攻时，第i只量子哈里斯鹰第m维的量子旋转角公式应为
Figure FDA0003377598370000049
其中，
Figure FDA00033775983700000410
Figure FDA00033775983700000411
为[0,1]间的随机数；
策略二：渐进式软围攻，当|E|≥0.5且
Figure FDA00033775983700000412
时第i只量子哈里斯鹰执行此策略；
此时，第i只量子鹰有两种不同的方法更新其量子旋转角，首先采取方法A1，此时的量子旋转角为
Figure FDA00033775983700000413
根据此量子旋转角改变该量子鹰的量子位置并计算其适应度，如果该值小于迭代次数为t时该量子鹰的适应度值，则转而执行方法B1，方法B1的量子旋转角公式为
Figure FDA00033775983700000414
其中
Figure FDA00033775983700000415
是[0,1]之间的随机数，F为归一化后的levy飞行函数；
策略三：硬包围，当|E|＜0.5且
Figure FDA00033775983700000416
时，量子哈里斯鹰i将执行硬包围策略；
此时的量子旋转角为
Figure FDA00033775983700000417
策略四：渐进式硬包围，当|E|＜0.5且
Figure FDA00033775983700000418
时，量子鹰i将执行渐进式硬包围；
此时，该量子鹰会采取两种不同的方法进行量子旋转角的更新；首先采用方法A2，方法A2的量子旋转角为
Figure FDA00033775983700000419
用此量子旋转角更新该量子鹰的量子位置并代入适应度函数计算其适应度，如果该值小于迭代次数为t时该量子鹰的适应度值，则转而执行方法B2；方法B2中第i只量子哈里斯鹰第m维的量子旋转角计算公式为
Figure FDA0003377598370000051
用以上策略计算得出的量子旋转角更新量子哈里斯鹰的量子位置，其量子位置更新公式为
Figure FDA0003377598370000052
根据更新后的量子位置
Figure FDA0003377598370000053
带入适应度函数中计算
Figure FDA0003377598370000054
和更新
Figure FDA0003377598370000055
5.根据权利要求1所述的冲击噪声环境下互质阵列的幅相误差校正和测向方法，其特征在于，步骤五具体包括：将步骤四中最终输出的最佳量子鹰的位置映射为相位误差估计值
Figure FDA0003377598370000056
至
Figure FDA0003377598370000057
从而得出相位误差估计矩阵
Figure FDA0003377598370000058
同时提取步骤二中粗估计的幅相误差阵的各元素的模值，记作abs(Γb)，其中abs(·)表示求括号内矩阵元素的模值的操作；生成最终的幅相误差估计矩阵为Cg＝Γg·abs(Γb)；
互质阵列的接收信号矢量可表示为
Figure FDA0003377598370000059
其中
Figure FDA00033775983700000510
是Q×1维的信号的第o次快拍数据矢量，o∈[1,O]，Q为信号源的个数；
Figure FDA00033775983700000511
是互质阵列在第o次快拍下所接收到的冲击噪声矢量；
Figure FDA00033775983700000512
是Q个信号源的导向矢量阵，其中
Figure FDA00033775983700000513
是第q个信源对应的导向矢量，且q∈[1,Q]；设
Figure FDA00033775983700000514
是接收信号的分数低阶协方差矩阵，其中各个位置的元素可根据互质阵列的接收数据进行估算，具体的估算公式为
Figure FDA00033775983700000515
其中
Figure FDA00033775983700000516
Figure FDA00033775983700000517
分别是第g个、第h个阵元接收到的第o次快拍数据；随后对该共变矩阵进行幅相误差补偿，幅相误差补偿后的分数低阶协方差矩阵应为
Figure FDA00033775983700000518
6.根据权利要求1所述的冲击噪声环境下互质阵列的幅相误差校正和测向方法，其特征在于，步骤六具体包括：依据步骤二中的方法，将幅相误差补偿后的分数低阶协方差矩阵C0进行向量化操作，截取其中虚拟均匀阵元的虚拟接收数据，组成一个
Figure FDA00033775983700000519
维的虚拟接收向量
Figure FDA00033775983700000520
再对其按照步骤二中的方法进行空间平滑算法，可以得到一个
Figure FDA00033775983700000521
维的矩阵
Figure FDA00033775983700000522
该虚拟快拍数据矩阵的协方差矩阵可估算为
Figure FDA0003377598370000061
这里的
Figure FDA0003377598370000062
是矩阵
Figure FDA0003377598370000063
的第l维列向量；
令量子哈里斯鹰的种群规模为Nb，最大迭代次数为Tb，第t次迭代时第i只哈里斯鹰的量子位置设定为
Figure FDA0003377598370000064
初次迭代时
Figure FDA0003377598370000065
是一个[0,1]间的随机数，Q为信源个数，量子鹰的每一维估计一个信源的到达角，将第i只鹰的第q维转化为波达角为
Figure FDA0003377598370000066
其中，
Figure FDA0003377598370000067
是量子哈里斯鹰算法的搜索下限和上限，由第i只量子哈里斯鹰算法计算出的信源导向矢量为
Figure FDA0003377598370000068
其中，
Figure FDA0003377598370000069
根据求出的虚拟快拍数据的协方差矩阵Rss，可写出极大似然方程为
Figure FDA00033775983700000610
因此定义第t次迭代时第i只量子哈里斯鹰适应度函数为
Figure FDA00033775983700000611
定义第t次迭代时，量子哈里斯鹰种群的最佳适应度为
Figure FDA00033775983700000612
最佳适应度的鹰所在的量子位置为
Figure FDA00033775983700000613
CN202111421629.XA 2021-11-26 2021-11-26 冲击噪声环境下互质阵列的幅相误差校正和测向方法 Active CN114167347B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111421629.XA CN114167347B (zh) 2021-11-26 2021-11-26 冲击噪声环境下互质阵列的幅相误差校正和测向方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111421629.XA CN114167347B (zh) 2021-11-26 2021-11-26 冲击噪声环境下互质阵列的幅相误差校正和测向方法

Publications (2)

Publication Number Publication Date
CN114167347A true CN114167347A (zh) 2022-03-11
CN114167347B CN114167347B (zh) 2023-07-21

Family

ID=80481420

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111421629.XA Active CN114167347B (zh) 2021-11-26 2021-11-26 冲击噪声环境下互质阵列的幅相误差校正和测向方法

Country Status (1)

Country Link
CN (1) CN114167347B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116016085A (zh) * 2023-01-04 2023-04-25 哈尔滨工程大学 一种基于快速哈里斯鹰优化的otfs单音干扰抑制方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110007266A (zh) * 2019-04-22 2019-07-12 哈尔滨工程大学 一种冲击噪声下的任意阵列相干源测向方法
CN110940949A (zh) * 2019-12-11 2020-03-31 哈尔滨工程大学 强冲击噪声环境下基于量子企鹅搜索机制的互质阵列doa估计方法
CN113240069A (zh) * 2021-05-14 2021-08-10 江苏科技大学 一种基于改进哈里斯鹰算法的rbf神经网络优化方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110007266A (zh) * 2019-04-22 2019-07-12 哈尔滨工程大学 一种冲击噪声下的任意阵列相干源测向方法
CN110940949A (zh) * 2019-12-11 2020-03-31 哈尔滨工程大学 强冲击噪声环境下基于量子企鹅搜索机制的互质阵列doa估计方法
CN113240069A (zh) * 2021-05-14 2021-08-10 江苏科技大学 一种基于改进哈里斯鹰算法的rbf神经网络优化方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
DAYONG ZHAO 等: "Direction Finding of Maximum Likelihood Algorithm Using Artificial Bee Colony in the Impulsive Noise" *
MING DIAO 等: "Direction Finding of Propagator Method Using Adaptive Cultural Algorithm in the Impulsive Noise" *
高洪元 等: "重构分数低阶协方差的子空间拟合测向算法" *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116016085A (zh) * 2023-01-04 2023-04-25 哈尔滨工程大学 一种基于快速哈里斯鹰优化的otfs单音干扰抑制方法
CN116016085B (zh) * 2023-01-04 2023-10-10 哈尔滨工程大学 一种基于快速哈里斯鹰优化的otfs单音干扰抑制方法

Also Published As

Publication number Publication date
CN114167347B (zh) 2023-07-21

Similar Documents

Publication Publication Date Title
CN108872946B (zh) 导向矢量和协方差矩阵联合迭代的稳健波束形成方法
CN110208735B (zh) 一种基于稀疏贝叶斯学习的相干信号doa估计方法
CN108680891B (zh) 非均匀噪声条件下考虑互耦效应的doa估计方法
CN110113085B (zh) 一种基于协方差矩阵重构的波束形成方法及系统
CN107907852B (zh) 基于空间平滑的协方差矩阵秩最小化doa估计方法
CN107450047B (zh) 嵌套阵下基于未知互耦信息的压缩感知doa估计方法
CN110007266B (zh) 一种冲击噪声下的任意阵列相干源测向方法
CN107315162B (zh) 基于内插变换和波束形成的远场相干信号doa估计方法
CN110196410B (zh) 一种阵列天线主瓣干扰抑制方法及系统
CN110244272B (zh) 基于秩一去噪模型的波达方向估计方法
CN108872930B (zh) 扩展孔径二维联合对角化doa估计方法
CN111337873A (zh) 一种基于稀疏阵的doa估计方法
CN112255629A (zh) 基于联合uca阵列的序贯esprit二维不相干分布源参数估计方法
CN113835063B (zh) 一种无人机阵列幅相误差与信号doa联合估计方法
CN110196417B (zh) 基于发射能量集中的双基地mimo雷达角度估计方法
CN109696651B (zh) 一种基于m估计的低快拍数下波达方向估计方法
CN114167347B (zh) 冲击噪声环境下互质阵列的幅相误差校正和测向方法
Luo et al. Mainlobe anti-jamming via eigen-projection processing and covariance matrix reconstruction
CN109212466B (zh) 一种基于量子蜻蜓演化机制的宽带测向方法
CN112668155B (zh) 一种基于二次重构的稳健波束形成方法及系统
CN116774162A (zh) 用于均匀线性阵列的抗主副瓣干扰自适应单脉冲测角方法
CN109683128B (zh) 冲击噪声环境下的单快拍测向方法
CN115453487A (zh) 一种相控阵雷达鲁棒波束形成方法
CN115079090A (zh) 基于esprit与加权降维搜索的多阵列非圆源快速定位方法
CN115421098A (zh) 嵌套面阵下降维求根music的二维doa估计方法

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