CN113673419A - 适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法 - Google Patents

适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法 Download PDF

Info

Publication number
CN113673419A
CN113673419A CN202110953837.8A CN202110953837A CN113673419A CN 113673419 A CN113673419 A CN 113673419A CN 202110953837 A CN202110953837 A CN 202110953837A CN 113673419 A CN113673419 A CN 113673419A
Authority
CN
China
Prior art keywords
signal
power
matrix
subband
covariance 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
CN202110953837.8A
Other languages
English (en)
Other versions
CN113673419B (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.)
Northwestern Polytechnical University
Original Assignee
Northwestern Polytechnical 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 Northwestern Polytechnical University filed Critical Northwestern Polytechnical University
Priority to CN202110953837.8A priority Critical patent/CN113673419B/zh
Publication of CN113673419A publication Critical patent/CN113673419A/zh
Application granted granted Critical
Publication of CN113673419B publication Critical patent/CN113673419B/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
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/22—Source localisation; Inverse modelling
    • 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/80—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 ultrasonic, sonic or infrasonic waves
    • G01S3/802—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/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • 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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/29—Graphical models, e.g. Bayesian networks

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Biology (AREA)
  • Signal Processing (AREA)
  • Evolutionary Computation (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Computing Systems (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)

Abstract

本发明涉及一种适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法，采用MVDR‑DL对强干扰进行抑制，构建关于MVDR‑DL波束功率输出线性关系的贝叶斯概率模型，将稀疏贝叶斯算法推广至波束域。在贝叶斯框架下进行DOA估计，避免了超参数选取的问题。本发明在进行DOA估计时，每次迭代仅更新一个网格点对应的信号参数，从而避免了矩阵求逆运算，有效降低了方法的计算量。在强干扰环境下实现了对目标信号方位有效估计的同时，增强了算法的实用性。

Description

适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法
技术领域
本发明属于信号处理等领域，特别涉及一种适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法。
背景技术
基于阵列接收信号的水下目标方位(Direction of arrival，DOA)估计是被动声纳信号处理中一项主要任务。近十几年稀疏重构类DOA估计算法因其对快拍数和信噪比的要求较低受到了广泛的关注。这类算法将空间划分为离散的网格，假设信号在有限个网格点上。通过对每个网格点上的信号参数进行估计，实现DOA估计。根据估计原理不同，算法可分为基于lp范数的算法和稀疏贝叶斯算法。相比于基于lp范数的算法来说，稀疏贝叶斯算法无需选取任何超参数，因此在实际信号处理中更易实现。大多数稀疏贝叶斯算法采用期望最大算法或变分贝叶斯算法迭代估计DOA值，需要大量的矩阵求逆运算，算法运算量较大。国外学者M.E.Tipping等人提出一种快速稀疏贝叶斯方法(M.E.Tipping and A.C.Faul,“Fast marginal likelihood maximisation for sparse Bayesian models,”inProc.9th Int.Workshop Artif.Intell.Stat.,2003,vol.1.)，该方法在每次迭代时仅更新一个网格对应的信号参数，从而避免了矩阵求逆，大大提升了计算效率。不同于主动声纳自主发射信号并通过接收反射回波进行目标探测，被动声纳是通过接收舰船辐射噪声来进行目标探测的，因而具有更好的隐蔽性。随着海洋战略和经济的不断发展，水面及水下舰船越来越多，这使得被动声纳的接收信号越来越复杂。当功率较低的目标信号远离接收阵时，距离接收阵较近的大功率水面舰船的辐射噪声等可看作为强干扰信号。强干扰信号的存在将影响对目标信号的DOA估计精度，甚至掩蔽目标信号。由于强干扰的存在，导致稀疏重构类算法对目标信号的方位估计性能存在一定程度的下降。
Yang等人(Y.Yang,Y.Zhang,and L.Yang,“Wideband sparse spatial spectrumestimation using matrix filter with nulling in a strong interferenceenvironment,”J.Acoust.Soc.Am.143(6),3891–3898(2018).)提出一种基于零陷矩阵滤波器的稀疏谱估计方法，使用零陷矩阵滤波器在干扰方向上形成深凹槽抑制强干扰，并利用稀疏谱估计方法进行弱目标信号方位估计。然而，设计零陷矩阵滤波器需求解一个凸优化问题，算法计算量较大。同时，稀疏谱估计算法属于lp范数类算法，需选取一个合适的超参数保证算法性能，然而该参数选取通常较为困难，造成该算法在实际应用时存在一定局限性。
发明内容
本发明解决的技术问题是：为了高效解决强干扰环境下对目标信号的DOA估计问题，本发明给出一种基于波束功率输出的快速稀疏贝叶斯(Fast sparse Bayesianlearning based on beamformer power outputs，FSBL-BPO)方法。该方法采用基于对角加载的最小方差无失真响应(Minimum variance distortionless response with diagonalloading，MVDR-DL)波束形成器作为预处理器，在强干扰方向形成凹槽充分抑制强干扰信号，并计算波束功率输出。构建适用于MVDR-DL波束功率输出和波束响应线性关系的贝叶斯概率模型。在贝叶斯框架下对模型参数进行迭代更新，每次迭代时仅更新一个网格对应的信号参数，避免大量矩阵求逆计算，实现强干扰环境下对目标信号的快速方位估计。
本发明的技术方案是：适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法，包括以下步骤：
步骤1：假设KS个目标信号和KD个干扰信号分别从
Figure BDA0003219614900000021
和
Figure BDA0003219614900000022
方向上入射至M元的均匀线阵列，目标信号和干扰信号之间互不相关；当阵列接收到信号后，将接收信号划分为N段，对每一段进行傅里叶变换后，宽带信号划分为L个子带；第l个子带上第n段阵列接收信号对应的傅里叶变换系数记为
Figure BDA0003219614900000023
采样协方差矩阵计算为
Figure BDA0003219614900000031
上标“H”为共轭转置运算；
步骤2：建立波束域模型，包括以下步骤：
步骤2.1：通过常规波束形成算法确定目标信号所在方位区域ΘS＝[ΘSL,ΘSR]，其中ΘSL和ΘSR分别为区域ΘS的左右边界；将该区域均匀划分为KB个网格
Figure BDA0003219614900000032
Figure BDA0003219614900000033
对于第l个子带，该区域上KB个MVDR-DL波束形成器组成的波束形成矩阵Wl为
Figure BDA0003219614900000034
其中
Figure BDA0003219614900000035
为指向φk的MVDR-DL波束形成器的加权量，al(φk)为第l个子带上指向φk的阵列流形，
Figure BDA00032196149000000318
为通过求解
Figure BDA0003219614900000038
的(M-KS-KD)个小特征值的平均值得到的噪声功率估计值，IM为M维的单位矩阵，上标“-1”表示矩阵求逆；
通过该矩阵对步骤1得到的采样协方差矩阵进行滤波，得到第l个子带上波束域协方差矩阵
Figure BDA0003219614900000039
为
Figure BDA00032196149000000310
其中Pl S和Pl D表示第l个子带上的目标信号和干扰信号的协方差矩阵，σl表示第l个子带上的噪声功率，
Figure BDA00032196149000000311
和
Figure BDA00032196149000000312
表示第l个子带上目标信号和干扰信号的阵列流形矩阵，El为第l个子带上的扰动误差；
步骤2.2：对步骤2.1得到的协方差矩阵进行按列向量化运算，得到
Figure BDA00032196149000000313
其中
Figure BDA00032196149000000314
和
Figure BDA00032196149000000315
分别表示第l个子带上目标信号和干扰信号的功率向量，
Figure BDA00032196149000000316
和
Figure BDA00032196149000000317
分别表示
Figure BDA0003219614900000041
和
Figure BDA0003219614900000042
矩阵按列向量化的向量，⊙表示Khatri-Rao乘积，上标“*”为共轭运算；
步骤2.3：定义矩阵
Figure BDA0003219614900000043
对于第m行第n列的元素[J]mn，若
Figure BDA0003219614900000044
[J]mn＝1，否则[J]mn＝0；将步骤2.2得到结果与该矩阵相乘，因MVDR-DL在区域ΘS外对干扰的波束响应较低，因此第l个子带上波束功率输出的线性关系式表示为：
Figure BDA0003219614900000045
其中
Figure BDA0003219614900000046
步骤2.4：将区域ΘS均匀划分为KG个网格，网格点
Figure BDA0003219614900000047
组成的向量记为
Figure BDA0003219614900000048
基于该离散网格，步骤2.3中的公式重新表示为
Figure BDA0003219614900000049
式中
Figure BDA00032196149000000410
Figure BDA00032196149000000411
为第l个子带上该网格上的阵列流形矩阵。pl是一个稀疏向量，当
Figure BDA00032196149000000412
pl第m个元素等于
Figure BDA00032196149000000413
第n个元素，否则为0；
步骤3：在贝叶斯框架下进行迭代计算，最终输出DOA估计值，包括以下子步骤：
步骤3.1：构建贝叶斯概率模型，其中扰动概率模型为：
Figure BDA00032196149000000414
其中N(·)表示实高斯分布，
Figure BDA00032196149000000415
Figure BDA00032196149000000416
表示Hadamard积；
信号概率模型表示为：
p(pl；γ)＝N(0,Γ-1),l＝1,...,L
其中
Figure BDA00032196149000000417
为信号稀疏参数，上标“T”为转置运算，Γ＝diag(γ)为以γ中元素为对角元素的对角矩阵；
构建好后，给定信号稀疏参数γ，信号功率后验协方差矩阵Σl，信号功率后验均值向量μl，噪声功率σl，以及参数Sl,p，
Figure BDA00032196149000000418
迭代初始值；
步骤3.2：参数迭代求解，分别完成对信号稀疏参数γ，信号功率后验协方差矩阵Σl，信号功率后验均值向量μl，噪声功率σl，以及参数Sl,p，
Figure BDA0003219614900000051
的更新；
在第i次迭代中，首先对γ进行更新：
Figure BDA0003219614900000052
其中上标(i)表示第i次迭代，
Figure BDA0003219614900000053
表示γ(i)第p个元素，
Figure BDA0003219614900000054
Figure BDA0003219614900000055
p(i)表示使函数
Figure BDA0003219614900000056
最大时的下标，
Figure BDA0003219614900000057
的值为：当
Figure BDA0003219614900000058
时，
Figure BDA0003219614900000059
否则
Figure BDA00032196149000000510
为∞，∞表示无穷大；|·|表示取绝对值；记存放前(i-1)次迭代中更新下标的活跃集为Φ(i-1)，对更新分为三种情况：
(1)若
Figure BDA00032196149000000511
时，则Φ(i)＝Φ(i-1)∪{p(i)}
(2)若
Figure BDA00032196149000000512
时，则Φ(i)＝Φ(i-1)
(3)若
Figure BDA00032196149000000513
Φ(i)为将p(i)从Φ(i-1)中移除的集合；
对应三种不同情况，首先完成对信号功率后验协方差矩阵Σl,l＝1,...,L的更新，再进行参数Sl,p，
Figure BDA00032196149000000514
的更新；最终先后完成对信号功率后验均值向量μl,l＝1,...,L以及噪声功率σl,l＝1,...,L的更新；
更新后的参数值若满足迭代终止条件，则得到所估计的功率谱为
Figure BDA00032196149000000515
该功率谱中峰值对应的方位即为目标信号的DOA估计值；
若不满足迭代终止条件，则继续进行更新，直至满足迭代终止条件。
本发明进一步的技术方案是；所述步骤3.1中的各参数迭代初始值为：
信号稀疏参数的初始值γ(0)：
Figure BDA0003219614900000061
其中
Figure BDA0003219614900000062
表示γ(0)第p个元素；
Figure BDA0003219614900000063
Figure BDA0003219614900000064
表示
Figure BDA0003219614900000065
的第p列，p(0)表示使函数
Figure BDA0003219614900000066
最大时
Figure BDA0003219614900000067
对应的下标，
Figure BDA0003219614900000068
的值为：当Δ＞0时，
Figure BDA0003219614900000069
否则
Figure BDA00032196149000000610
活跃集的初始值Φ(0)：
Φ(0)＝{p(0)}
第l个子带上信号功率后验协方差矩阵的初始值
Figure BDA00032196149000000611
Figure BDA00032196149000000612
第l个子带上信号功率后验均值向量的初始值
Figure BDA00032196149000000613
Figure BDA00032196149000000614
第l个子带上噪声功率的初始值
Figure BDA00032196149000000615
Figure BDA00032196149000000616
参数Sl,p，
Figure BDA00032196149000000617
的初始值：
Figure BDA00032196149000000618
本发明进一步的技术方案是：所述步骤3.2中对应第一种情况时，第i次迭代中第l个子带上信号功率后验协方差矩阵更新为
Figure BDA0003219614900000071
其中
Figure BDA0003219614900000072
仅包含活跃集Φ(i-1)中记录的“激活”网格对应的阵列流形，
Figure BDA0003219614900000073
Figure BDA0003219614900000074
和
Figure BDA0003219614900000075
更新为
Figure BDA0003219614900000076
其中
Figure BDA0003219614900000077
本发明进一步的技术方案是：所述步骤3.2中对应第二种情况时，第i次迭代中第l个子带上信号功率后验协方差矩阵更新为
Figure BDA0003219614900000078
其中
Figure BDA0003219614900000079
为
Figure BDA00032196149000000710
的第p(i)列，
Figure BDA00032196149000000711
表示
Figure BDA00032196149000000712
在第p(i)行第p(i)列的元素，
Figure BDA00032196149000000713
为除了第p(i)个元素为1其余元素为0的向量；
Figure BDA00032196149000000714
和
Figure BDA00032196149000000715
更新为
Figure BDA00032196149000000716
本发明进一步的技术方案是：所述步骤3.2中对应第三种情况时，
Figure BDA00032196149000000717
将该公式代入第二种情况中进行计算，即可完成更新。
本发明进一步的技术方案是：所述步骤3.2中完成信号功率后验协方差矩阵Σl以及参数Sl,p，
Figure BDA0003219614900000081
的更新后，对第l个子带上信号功率后验均值向量μl的更新为
Figure BDA0003219614900000082
最后，第l个子带上噪声功率σl的更新为
Figure BDA0003219614900000083
本发明进一步的技术方案是：所述步骤3.2中的迭代终止条件为：当迭代满足||γ(i)-γ(i-1)||2/||γ(i-1)||2≤10-3，或者迭代次数大于Itermax＝1000时，其中||·||2表示l2范数，迭代终止。
发明效果
本发明的技术效果在于：采用MVDR-DL波束形成器代替零陷矩阵滤波器来充分抑制强干扰信号，避免了强干扰信号对后续目标信号方位估计影响的同时，降低了计算量；计算MVDR-DL波束功率输出，构建了关于波束功率输出与波束响应的线性关系的贝叶斯概率模型，将稀疏贝叶斯类算法推广至波束域中，提升了稀疏贝叶斯算法在强干扰环境下对目标信号方位估计的稳健性；在贝叶斯框架下进行参数迭代更新，避免了超参数选取的问题；每次迭代时仅更新一个网格点对应的信号参数，避免大量矩阵求逆计算，从而在有效解决强干扰环境下对目标信号方位估计问题的同时，实现了快速DOA估计。
附图说明
图1基于波束功率输出的快速稀疏贝叶斯方法的总体流程框图
图2 FSBL-BPO迭代流程
图3采用常规波束形成算法估计目标所在区域结果
图4快速贝叶斯方法的方位估计结果
图5基于零陷矩阵滤波器的稀疏谱估计方法方位估计结果
图6 FSBL-BPO方法的方位估计结果
具体实施方式
理解的是，术语“中心”、“纵向”、“横向”、“长度”、“宽度”、“厚度”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”、“顺时针”、“逆时针”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。
参见图1-图6，本发明解决其技术问题所采用的技术方案包括以下步骤：
步骤1：假设KS个目标信号和KD个干扰信号分别从
Figure BDA0003219614900000091
和
Figure BDA0003219614900000092
方向上入射至M元的均匀线阵列，目标信号和干扰信号之间互不相关。当阵列接收到信号后，将接收信号划分为N段，对每一段进行傅里叶变换后，宽带信号划分为L个子带。第l个子带上第n段阵列接收信号对应的傅里叶变换系数记为
Figure BDA0003219614900000093
采样协方差矩阵计算为
Figure BDA0003219614900000094
上标“H”为共轭转置运算。
步骤2：建立波束域模型，包括以下步骤：
步骤2.1：通过常规波束形成算法确定目标信号所在方位区域ΘS＝[ΘSL,ΘSR]，其中ΘSL和ΘSR分别为区域ΘS的左右边界。将该区域均匀划分为KB个网格
Figure BDA0003219614900000095
Figure BDA0003219614900000096
对于第l个子带，该区域上KB个MVDR-DL波束形成器组成的波束形成矩阵Wl为
Figure BDA0003219614900000097
其中
Figure BDA0003219614900000101
为指向φk的MVDR-DL波束形成器的加权量，al(φk)为第l个子带上指向φk的阵列流形，
Figure BDA0003219614900000102
为通过求解
Figure BDA0003219614900000103
的(M-KS-KD)个小特征值的平均值得到的噪声功率估计值，IM为M维的单位矩阵，上标“-1”表示矩阵求逆。
通过该矩阵对步骤1得到的采样协方差矩阵进行滤波，得到第l个子带上的波束域协方差矩阵
Figure BDA0003219614900000104
为
Figure BDA0003219614900000105
其中Pl S和Pl D表示第l个子带上的目标信号和干扰信号的协方差矩阵，σl表示第l个子带上的噪声功率，
Figure BDA0003219614900000106
和
Figure BDA0003219614900000107
表示第l个子带上目标信号和干扰信号的阵列流形矩阵，El为第l个子带上的扰动误差。
步骤2.2：对步骤2.1得到的协方差矩阵进行按列向量化运算，得到
Figure BDA0003219614900000108
其中
Figure BDA0003219614900000109
和
Figure BDA00032196149000001010
分别表示第l个子带上目标信号和干扰信号的功率向量，
Figure BDA00032196149000001011
和
Figure BDA00032196149000001012
分别表示Wl HW和Wl HElW矩阵按列向量化的向量，⊙表示Khatri-Rao乘积，上标“*”为共轭运算。
步骤2.3：定义矩阵
Figure BDA00032196149000001013
对于第m行第n列的元素[J]mn，若
Figure BDA00032196149000001014
[J]mn＝1，否则[J]mn＝0；将步骤2.2得到结果与该矩阵相乘，因MVDR-DL在区域ΘS外对干扰的波束响应较低，因此第l个子带上波束功率输出的线性关系式表示为：
Figure BDA00032196149000001015
其中
Figure BDA00032196149000001016
步骤2.4：将区域ΘS均匀划分为KG个网格，网格点
Figure BDA00032196149000001017
组成的向量记为
Figure BDA00032196149000001018
基于该离散网格，步骤2.3中的公式重新表示为
Figure BDA0003219614900000111
式中
Figure BDA0003219614900000112
Figure BDA0003219614900000113
为第l个子带上该网格上的阵列流形矩阵。pl是一个稀疏向量，当
Figure BDA0003219614900000114
pl第m个元素等于
Figure BDA0003219614900000115
第n个元素，否则为0。
步骤3：在贝叶斯框架下进行迭代计算，最终输出DOA估计值，包括以下子步骤：步骤3.1：构建贝叶斯概率模型，其中扰动概率模型为：
Figure BDA0003219614900000116
其中N(·)表示实高斯分布，
Figure BDA0003219614900000117
Figure BDA0003219614900000118
表示Hadamard积。信号概率模型表示为：
p(pl；γ)＝N(0,Γ-1),l＝1,...,L
其中
Figure BDA0003219614900000119
为信号稀疏参数，上标“T”为转置运算，Γ＝diag(γ)为以γ中元素为对角元素的对角矩阵。
构建好后，给定信号稀疏参数γ，信号功率后验协方差矩阵Σl，信号功率后验均值向量μl，噪声功率σl，以及参数Sl,p，
Figure BDA00032196149000001110
迭代初始值，如下信号稀疏参数的初始值γ(0)：
Figure BDA00032196149000001111
其中
Figure BDA00032196149000001112
表示γ(0)第p个元素，
Figure BDA00032196149000001113
Figure BDA00032196149000001114
表示
Figure BDA00032196149000001115
的第p列，p(0)表示使函数
Figure BDA00032196149000001116
最大时
Figure BDA00032196149000001117
对应的下标，
Figure BDA00032196149000001118
的值为：当Δ＞0时，
Figure BDA00032196149000001119
否则
Figure BDA00032196149000001120
活跃集的初始值Φ(0)：
Φ(0)＝{p(0)}
Figure BDA0003219614900000121
第l个子带上信号功率后验协方差矩阵的初始值
Figure BDA0003219614900000122
Figure BDA0003219614900000123
第l个子带上信号功率后验均值向量的初始值
Figure BDA0003219614900000124
Figure BDA0003219614900000125
第l个子带上噪声功率的初始值
Figure BDA0003219614900000126
Figure BDA0003219614900000127
参数Sl,p，
Figure BDA0003219614900000128
的初始值：
Figure BDA0003219614900000129
步骤3.2：参数迭代求解，分别完成对信号稀疏参数γ，信号功率后验协方差矩阵Σl，信号功率后验均值向量μl，噪声功率σl，以及参数Sl,p，
Figure BDA00032196149000001210
的更新。
在第i次迭代中，首先对γ进行更新：
Figure BDA00032196149000001211
其中上标(i)表示第i次迭代，
Figure BDA00032196149000001212
表示γ(i)第p个元素，
Figure BDA00032196149000001213
Figure BDA00032196149000001214
p(i)表示使函数
Figure BDA00032196149000001215
最大时的下标，
Figure BDA0003219614900000131
的值为：当
Figure BDA0003219614900000132
时，
Figure BDA0003219614900000133
否则
Figure BDA0003219614900000134
为∞，∞表示无穷大。|·|表示取绝对值。记存放前(i-1)次迭代中更新下标的活跃集为Φ(i-1)，对更新分为三种情况：
(1)若
Figure BDA0003219614900000135
时，则Φ(i)＝Φ(i-1)∪{p(i)}
(2)若
Figure BDA0003219614900000136
时，则Φ(i)＝Φ(i-1)
(3)若
Figure BDA0003219614900000137
Φ(i)为将p(i)从Φ(i-1)中移除的集合。
对应三种不同情况，首先完成对信号功率后验协方差矩阵Σl,l＝1,...,L的更新，再进行参数Sl,p，
Figure BDA0003219614900000138
的更新；最终先后完成对信号功率后验均值向量μl,l＝1,...,L以及噪声功率σl,l＝1,...,L的更新；
针对第一种情况，第i次迭代中第l个子带上信号功率后验协方差矩阵更新为
Figure BDA0003219614900000139
其中
Figure BDA00032196149000001310
仅包含活跃集Φ(i-1)中记录的“激活”网格对应的阵列流形，
Figure BDA00032196149000001311
Figure BDA00032196149000001312
和
Figure BDA00032196149000001313
更新为
Figure BDA00032196149000001314
其中
Figure BDA00032196149000001315
针对第二种情况，第i次迭代中第l个子带上信号功率后验协方差矩阵更新为
Figure BDA00032196149000001316
其中
Figure BDA00032196149000001317
为
Figure BDA00032196149000001318
的第p(i)列，
Figure BDA00032196149000001319
表示
Figure BDA00032196149000001320
在第p(i)行第p(i)列的元素，
Figure BDA0003219614900000141
为除了第p(i)个元素为1其余元素为0的向量。
Figure BDA0003219614900000142
和
Figure BDA0003219614900000143
更新为
Figure BDA0003219614900000144
针对第三种情况，
Figure BDA0003219614900000145
将该公式代入第二种情况中进行计算，即可完成更新。
完成信号功率后验协方差矩阵Σl以及参数Sl,p，
Figure BDA0003219614900000146
的更新后，对第l个子带上信号功率后验均值向量μl的更新为
Figure BDA0003219614900000147
最后，第l个子带上噪声功率σl的更新为
Figure BDA0003219614900000148
更新后的参数值若满足||γ(i)-γ(i-1)||2/||γ(i-1)||2≤10-3，其中||·||2表示l2范数，或者迭代次数大于Itermax＝1000时，迭代终止，所估计的功率谱为
Figure BDA0003219614900000149
该功率谱中峰值对应的方位即为目标信号的DOA估计值；若不满足迭代终止条件，则继续进行更新，直至满足迭代终止条件。
为了验证利用本发明给出的方法在强干扰环境下对目标信号方位估计的有效性，设计仿真实验如下：假设接收阵列为M＝20元均匀线列阵，阵元间距为4m。两个远处航船从-10°和-7°方向上入射至阵列，功率为0dB，视为目标信号，即KS＝2；一个近处水面舰从10°方向上入射至阵列，功率为20dB，大于远处航船，视为干扰信号，即KD＝1。所考虑的信号频段为[90，180]Hz。噪声为高斯白噪声，其功率在所考虑的信号频段上为0dB。
将接收信号均匀划分为N＝50段，在每一段上进行傅里叶变换，将信号划分至L＝46个子带。图3为采用常规波束形成算法对目标信号所在区域估计的结果，虚线为目标信号所在区域范围ΘS的边界ΘSL＝-14°和ΘSR＝-2°。故ΘS为[-14,-2]°。
将该区域以2°为间隔均匀划分为7个网格点，得到MVDR-DL波束指向角，即KB＝7。为了进行DOA估计，将该区域以1°为间隔均匀划分为13个网格点，即KG＝13，在该网格上进行DOA估计。图4、图5和图6分别为快速贝叶斯方法、基于零陷矩阵滤波器稀疏谱估计方法和FSBL-BPO方法的方位估计结果，图中的虚线为目标信号的真实方位。可以看出，快速贝叶斯方法(图4)由于受到了强干扰信号的影响，没有办法估计出两个信号。本发明给出的方法(图6)和基于零陷矩阵滤波器的稀疏谱估计方法(图5)避免了强干扰信号的影响，可以很好地估计出两个目标信号，说明本方法成功将稀疏贝叶斯算法推广至波束域中，增强了稀疏贝叶斯算法在强干扰环境下对目标信号方位估计的稳健性。除此之外，本发明给出的方法进行一次DOA估计的运行时间为0.2s，而基于零陷矩阵滤波器的稀疏谱估计方法的运行时间为700s。可以看出，本发明给出的方法大大提升了计算效率，实现了快速DOA估计。

Claims (7)

1.适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法，其特征在于，包括以下步骤：
步骤1：假设KS个目标信号和KD个干扰信号分别从
Figure FDA0003219614890000011
和
Figure FDA0003219614890000012
方向上入射至M元的均匀线阵列，目标信号和干扰信号之间互不相关；当阵列接收到信号后，将接收信号划分为N段，对每一段进行傅里叶变换后，宽带信号划分为L个子带；第l个子带上第n段阵列接收信号对应的傅里叶变换系数记为
Figure FDA0003219614890000013
采样协方差矩阵计算为
Figure FDA0003219614890000014
上标“H”为共轭转置运算；
步骤2：建立波束域模型，包括以下步骤：
步骤2.1：通过常规波束形成算法确定目标信号所在方位区域ΘS＝[ΘSL,ΘSR]，其中ΘSL和ΘSR分别为区域ΘS的左右边界；将该区域均匀划分为KB个网格
Figure FDA0003219614890000015
Figure FDA0003219614890000016
对于第l个子带，该区域上KB个MVDR-DL波束形成器组成的波束形成矩阵Wl为
Figure FDA0003219614890000017
其中
Figure FDA0003219614890000018
为指向φk的MVDR-DL波束形成器的加权量，al(φk)为第l个子带上指向φk的阵列流形，
Figure FDA0003219614890000019
Figure FDA00032196148900000110
为通过求解
Figure FDA00032196148900000111
的(M-KS-KD)个小特征值的平均值得到的噪声功率估计值，IM为M维的单位矩阵，上标“-1”表示矩阵求逆；
通过该矩阵对步骤1得到的采样协方差矩阵进行滤波，得到第l个子带上波束域协方差矩阵
Figure FDA00032196148900000112
为
Figure FDA00032196148900000113
其中Pl S和Pl D表示第l个子带上的目标信号和干扰信号的协方差矩阵，σl表示第l个子带上的噪声功率，
Figure FDA0003219614890000021
和
Figure FDA0003219614890000022
表示第l个子带上目标信号和干扰信号的阵列流形矩阵，El为第l个子带上的扰动误差；
步骤2.2：对步骤2.1得到的协方差矩阵进行按列向量化运算，得到
Figure FDA0003219614890000023
其中
Figure FDA0003219614890000024
和
Figure FDA0003219614890000025
分别表示第l个子带上目标信号和干扰信号的功率向量，
Figure FDA0003219614890000026
和
Figure FDA0003219614890000027
分别表示Wl HW和Wl HElW矩阵按列向量化的向量，⊙表示Khatri-Rao乘积，上标“*”为共轭运算；
步骤2.3：定义矩阵
Figure FDA0003219614890000028
对于第m行第n列的元素[J]mn，若n＝KB(m-1)+m，[J]mn＝1，否则[J]mn＝0；将步骤2.2得到结果与该矩阵相乘，因MVDR-DL在区域ΘS外对干扰的波束响应较低，因此第l个子带上波束功率输出的线性关系式表示为：
Figure FDA0003219614890000029
其中
Figure FDA00032196148900000210
步骤2.4：将区域ΘS均匀划分为KG个网格，网格点
Figure FDA00032196148900000211
组成的向量记为
Figure FDA00032196148900000212
基于该离散网格，步骤2.3中的公式重新表示为
Figure FDA00032196148900000213
式中
Figure FDA00032196148900000214
Figure FDA00032196148900000215
为第l个子带上该网格上的阵列流形矩阵。pl是一个稀疏向量，当
Figure FDA00032196148900000216
pl第m个元素等于
Figure FDA00032196148900000217
第n个元素，否则为0；
步骤3：在贝叶斯框架下进行迭代计算，最终输出DOA估计值，包括以下子步骤：
步骤3.1：构建贝叶斯概率模型，其中扰动概率模型为：
Figure FDA00032196148900000218
其中N(·)表示实高斯分布，
Figure FDA00032196148900000219
Figure FDA00032196148900000220
表示Hadamard积；
信号概率模型表示为：
p(pl；γ)＝N(0,Γ-1),l＝1,...,L
其中
Figure FDA0003219614890000031
为信号稀疏参数，上标“T”为转置运算，Γ＝diag(γ)为以γ中元素为对角元素的对角矩阵；
构建好后，给定信号稀疏参数γ，信号功率后验协方差矩阵Σl，信号功率后验均值向量μl，噪声功率σl，以及参数Sl,p，
Figure FDA0003219614890000032
Figure FDA0003219614890000033
迭代初始值；
步骤3.2：参数迭代求解，分别完成对信号稀疏参数γ，信号功率后验协方差矩阵Σl，信号功率后验均值向量μl，噪声功率σl，以及参数Sl,p，
Figure FDA0003219614890000034
Figure FDA0003219614890000035
的更新；
在第i次迭代中，首先对γ进行更新：
Figure FDA0003219614890000036
其中上标(i)表示第i次迭代，
Figure FDA0003219614890000037
表示γ(i)第p个元素，
Figure FDA0003219614890000038
Figure FDA0003219614890000039
p(i)表示使函数
Figure FDA00032196148900000310
最大时的下标，
Figure FDA00032196148900000311
的值为：当
Figure FDA00032196148900000312
时，
Figure FDA00032196148900000313
否则
Figure FDA00032196148900000314
为∞，∞表示无穷大；|·|表示取绝对值；
记存放前(i-1)次迭代中更新下标的活跃集为Φ(i-1)，对更新分为三种情况：
(1)若
Figure FDA00032196148900000315
时，则Φ(i)＝Φ(i-1)∪{p(i)}
(2)若
Figure FDA00032196148900000316
时，则Φ(i)＝Φ(i-1)
(3)若
Figure FDA00032196148900000317
Φ(i)为将p(i)从Φ(i-1)中移除的集合；
对应三种不同情况，首先完成对信号功率后验协方差矩阵Σl,l＝1,...,L的更新，再进行参数Sl,p，
Figure FDA0003219614890000041
Figure FDA0003219614890000042
的更新；最终先后完成对信号功率后验均值向量μl,l＝1,...,L以及噪声功率σl,l＝1,...,L的更新；
更新后的参数值若满足迭代终止条件，则得到所估计的功率谱为
Figure FDA0003219614890000043
该功率谱中峰值对应的方位即为目标信号的DOA估计值；
若不满足迭代终止条件，则继续进行更新，直至满足迭代终止条件。
2.如权利要求1所述的适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法，其特征在于，所述步骤3.1中的各参数迭代初始值为：
信号稀疏参数的初始值γ(0)：
Figure FDA0003219614890000044
其中
Figure FDA0003219614890000045
表示γ(0)第p个元素；
Figure FDA0003219614890000046
Figure FDA0003219614890000047
表示
Figure FDA0003219614890000048
的第p列，p(0)表示使函数
Figure FDA0003219614890000049
最大时
Figure FDA00032196148900000410
对应的下标，
Figure FDA00032196148900000411
的值为：当Δ＞0时，
Figure FDA00032196148900000412
否则
Figure FDA00032196148900000413
活跃集的初始值Φ(0)：
Φ(0)＝{p(0)}
第l个子带上信号功率后验协方差矩阵的初始值
Figure FDA00032196148900000414
Figure FDA00032196148900000415
第l个子带上信号功率后验均值向量的初始值
Figure FDA00032196148900000416
Figure FDA0003219614890000051
第l个子带上噪声功率的初始值
Figure FDA0003219614890000052
Figure FDA0003219614890000053
参数Sl,p，
Figure FDA0003219614890000054
Figure FDA0003219614890000055
的初始值：
Figure FDA0003219614890000056
3.如权利要求1所述的适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法，其特征在于，所述步骤3.2中对应第一种情况时，第i次迭代中第l个子带上信号功率后验协方差矩阵更新为
Figure FDA0003219614890000057
其中
Figure FDA0003219614890000058
仅包含活跃集Φ(i-1)中记录的“激活”网格对应的阵列流形，
Figure FDA0003219614890000059
Figure FDA00032196148900000510
和
Figure FDA00032196148900000511
更新为
Figure FDA00032196148900000512
其中
Figure FDA00032196148900000513
4.如权利要求1所述的适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法，其特征在于，所述步骤3.2中对应第二种情况时，第i次迭代中第l个子带上信号功率后验协方差矩阵更新为
Figure FDA0003219614890000061
其中
Figure FDA0003219614890000062
为
Figure FDA0003219614890000063
的第p(i)列，
Figure FDA0003219614890000064
Figure FDA0003219614890000065
表示
Figure FDA0003219614890000066
在第p(i)行第p(i)列的元素，
Figure FDA0003219614890000067
为除了第p(i)个元素为1其余元素为0的向量；
Figure FDA0003219614890000068
和
Figure FDA0003219614890000069
更新为
Figure FDA00032196148900000610
5.如权利要求1所述的适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法，其特征在于，所述步骤3.2中对应第三种情况时，
Figure FDA00032196148900000611
将该公式代入第二种情况中进行计算，即可完成更新。
6.如权利要求1所述的适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法，其特征在于，所述步骤3.2中完成信号功率后验协方差矩阵Σl以及参数Sl,p，
Figure FDA00032196148900000612
Figure FDA00032196148900000613
的更新后，对第l个子带上信号功率后验均值向量μl的更新为
Figure FDA00032196148900000614
最后，第l个子带上噪声功率σl的更新为
Figure FDA00032196148900000615
7.如权利要求1所述的适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法，其特征在于，所述步骤3.2中的迭代终止条件为：当迭代满足||γ(i)-γ(i-1)||2/||γ(i-1)||2≤10-3，或者迭代次数大于Itermax＝1000时，其中||·||2表示l2范数，迭代终止。
CN202110953837.8A 2021-08-19 2021-08-19 适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法 Active CN113673419B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110953837.8A CN113673419B (zh) 2021-08-19 2021-08-19 适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110953837.8A CN113673419B (zh) 2021-08-19 2021-08-19 适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法

Publications (2)

Publication Number Publication Date
CN113673419A true CN113673419A (zh) 2021-11-19
CN113673419B CN113673419B (zh) 2024-05-28

Family

ID=78543850

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110953837.8A Active CN113673419B (zh) 2021-08-19 2021-08-19 适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法

Country Status (1)

Country Link
CN (1) CN113673419B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114254265A (zh) * 2021-12-20 2022-03-29 军事科学院系统工程研究院网络信息研究所 基于统计流形距离的卫星通信干扰几何分析方法
CN116338574A (zh) * 2023-04-10 2023-06-27 哈尔滨工程大学 一种基于匹配波束的稀疏贝叶斯学习水下声源定位方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109116337A (zh) * 2018-07-30 2019-01-01 西北工业大学 一种基于矩阵滤波的稀疏近似最小方差方位估计方法
CN110579737A (zh) * 2019-07-17 2019-12-17 电子科技大学 一种杂波环境中基于稀疏阵列的mimo雷达宽带doa计算方法
CN112230226A (zh) * 2020-09-23 2021-01-15 浙江大学 基于贝叶斯压缩感知算法的自适应波束形成器设计方法
CN112487703A (zh) * 2020-11-09 2021-03-12 南京信息工程大学滨江学院 基于稀疏贝叶斯在未知噪声场的欠定宽带信号doa估计方法
CN112630760A (zh) * 2020-11-30 2021-04-09 海鹰企业集团有限责任公司 多目标条件下的强干扰抑制波束形成器设计方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109116337A (zh) * 2018-07-30 2019-01-01 西北工业大学 一种基于矩阵滤波的稀疏近似最小方差方位估计方法
CN110579737A (zh) * 2019-07-17 2019-12-17 电子科技大学 一种杂波环境中基于稀疏阵列的mimo雷达宽带doa计算方法
CN112230226A (zh) * 2020-09-23 2021-01-15 浙江大学 基于贝叶斯压缩感知算法的自适应波束形成器设计方法
CN112487703A (zh) * 2020-11-09 2021-03-12 南京信息工程大学滨江学院 基于稀疏贝叶斯在未知噪声场的欠定宽带信号doa估计方法
CN112630760A (zh) * 2020-11-30 2021-04-09 海鹰企业集团有限责任公司 多目标条件下的强干扰抑制波束形成器设计方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
冯杰;杨益新;孙超;: "自适应空域矩阵滤波器设计和目标方位估计", 系统仿真学报, no. 20, 20 October 2007 (2007-10-20) *
张赫;陈华伟;: "一种强干扰环境下的离格稀疏贝叶斯DOA估计方法", 数据采集与处理, no. 06, 15 November 2019 (2019-11-15) *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114254265A (zh) * 2021-12-20 2022-03-29 军事科学院系统工程研究院网络信息研究所 基于统计流形距离的卫星通信干扰几何分析方法
CN114254265B (zh) * 2021-12-20 2022-06-07 军事科学院系统工程研究院网络信息研究所 基于统计流形距离的卫星通信干扰几何分析方法
CN116338574A (zh) * 2023-04-10 2023-06-27 哈尔滨工程大学 一种基于匹配波束的稀疏贝叶斯学习水下声源定位方法
CN116338574B (zh) * 2023-04-10 2023-09-19 哈尔滨工程大学 一种基于匹配波束的稀疏贝叶斯学习水下声源定位方法

Also Published As

Publication number Publication date
CN113673419B (zh) 2024-05-28

Similar Documents

Publication Publication Date Title
CN110045323B (zh) 一种基于矩阵填充的互质阵稳健自适应波束形成算法
CN109490850B (zh) 主瓣干扰下宽带阵列自适应波束形成方法
CN109254261B (zh) 基于均匀圆阵epuma的相干信号零陷加深方法
CN109298383B (zh) 一种基于变分贝叶斯推断的互质阵波达方向角估计方法
CN113673419A (zh) 适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法
CN107576931B (zh) 一种基于协方差低维度迭代稀疏重构的相关/相干信号波达方向估计方法
CN105302936A (zh) 基于相关计算和协方差矩阵重构的自适应波束形成方法
Hossain et al. Efficient robust broadband antenna array processor in the presence of look direction errors
CN110988854A (zh) 基于交替方向乘子法的鲁棒自适应波束形成算法
CN112230226B (zh) 基于贝叶斯压缩感知算法的自适应波束形成器设计方法
CN113534151B (zh) 基于离网稀疏贝叶斯学习的双频段isar成像方法
CN107124216A (zh) 一种针对阵列误差的Capon稳健自适应波束形成方法及系统
CN113219402B (zh) 基于Modified-ALM算法的稀疏阵列DOA估计方法
CN110673119A (zh) 基于压缩感知的非正则化方位估计方法及系统
CN109541526A (zh) 一种利用矩阵变换的圆环阵方位估计方法
CN113673158B (zh) 适用于强干扰环境下的波束域变分贝叶斯方位估计方法
CN114844544B (zh) 一种基于低管秩张量分解的互质阵列波束成形方法、系统及介质
CN110161476A (zh) 基于幂迭代广义瑞利商算法的雷达波束形成方法
CN113422629B (zh) 一种协方差矩阵重构自适应波束形成方法及系统
CN110208830B (zh) 一种基于空时二维稀疏阵列的导航抗干扰方法
CN115130504A (zh) 基于稀疏贝叶斯学习的稳健波束形成方法
CN114329328B (zh) 一种基于非高斯非圆信号特性的稳健自适应波束形成方法
CN116990817B (zh) 一种雷达前视无网格重构sar成像方法和装置
CN111077493B (zh) 一种基于实值离格变分贝叶斯推理的nested阵列波达方向估计方法
CN117741559A (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