CN113673158B - 适用于强干扰环境下的波束域变分贝叶斯方位估计方法 - Google Patents

适用于强干扰环境下的波束域变分贝叶斯方位估计方法 Download PDF

Info

Publication number
CN113673158B
CN113673158B CN202110954450.4A CN202110954450A CN113673158B CN 113673158 B CN113673158 B CN 113673158B CN 202110954450 A CN202110954450 A CN 202110954450A CN 113673158 B CN113673158 B CN 113673158B
Authority
CN
China
Prior art keywords
signal
matrix
power
noise power
updating
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
CN202110954450.4A
Other languages
English (en)
Other versions
CN113673158A (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 CN202110954450.4A priority Critical patent/CN113673158B/zh
Publication of CN113673158A publication Critical patent/CN113673158A/zh
Application granted granted Critical
Publication of CN113673158B publication Critical patent/CN113673158B/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
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
    • 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/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N7/00—Computing arrangements based on specific mathematical models
    • G06N7/01—Probabilistic graphical models, e.g. probabilistic networks
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/08—Probabilistic or stochastic CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • 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

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Computation (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Software Systems (AREA)
  • Computing Systems (AREA)
  • Artificial Intelligence (AREA)
  • Algebra (AREA)
  • Probability & Statistics with Applications (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Medical Informatics (AREA)
  • Computer Hardware Design (AREA)
  • Geometry (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Databases & Information Systems (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)

Abstract

本发明涉及一种适用于强干扰环境下的波束域变分贝叶斯方位估计方法，有效解决强干扰环境下对目标信号的方位估计问题。方法采用基于对角加载的最小方差无失真响应波束形成器抑制强干扰信号，并计算波束功率输出。基于波束功率输出的线性关系构建后验概率模型，将变分贝叶斯推断推广至波束域。在贝叶斯框架下对目标方位进行估计，避免了选取合适超参数的难题，在强干扰环境下实现了对目标信号方位估计的同时，增强了算法的实用性。

Description

适用于强干扰环境下的波束域变分贝叶斯方位估计方法
技术领域
本发明属于信号处理等领域，特别涉及一种适用于强干扰环境下的波束域变分贝叶斯方位估计方法。
背景技术
基于阵列接收信号的水下目标方位(Direction of arrival，DOA)估计是被动声纳信号处理中一项主要任务。近十几年稀疏重构类DOA估计算法因其对快拍数和信噪比的要求较低受到了广泛的关注。根据原理不同，算法可分为基于lp范数的算法和稀疏贝叶斯算法。相比于基于lp范数的算法来说，稀疏贝叶斯算法无需选取任何超参数，因此在实际信号处理中更易实现。
在实际水下信号处理系统中，进行DOA估计前采用波束形成器对接收信号进行预处理是一个必要过程，通过波束形成可有效提升信噪比、减小后续计算量以及降低对系统存储量的要求。为在小快拍和低信噪比下保持高分辨性，许多基于波束功率输出的稀疏重构类算法被提出。这些算法采用常规波束形成器(Conventional beamforming，CBF)作为预处理器使感兴趣区域的信号通过，基于此构建关于波束功率输出和由波束响应构成的稀疏字典的线性关系，并采用基于lp范数的算法完成DOA估计。
然而上述算法均需选取一个合适的超参数进行DOA估计，该参数的选取通常较为困难，造成这些算法在实际声纳系统中难以使用。同时这些算法仅考虑窄带DOA估计问题，在现代声纳系统中，宽带信号处理被广泛使用。若将这些算法用于宽带DOA估计，可将接收信号通过傅里叶变换划分至各个子带，在各子带中分别进行窄带DOA估计。然而，这将忽略各子带信号的稀疏性是一致的这一先验信息，从而在一定程度上影响了DOA估计的精度。最重要的是，被动声纳主要通过接收舰船辐射噪声来进行目标探测，当目标信号功率较低时，其周围环境中存在的强干扰信号如拖船噪声将影响对目标信号的DOA估计，甚至掩蔽目标信号。若将基于CBF波束功率输出的lp范数算法用于上述强干扰环境下，CBF由于旁瓣较高，无法完全抑制强干扰信号，干扰剩余量严重影响后续DOA估计精度，甚至导致算法完全失效。
发明内容
本发明解决的技术问题是：为了在强干扰环境下对目标信号的方位进行有效估计，同时避免超参数的选取问题，本发明给出一种基于波束功率输出的变分贝叶斯(Variational Bayesian inference based on beamformer power outputs，VBI-BPO)方法。方法采用基于对角加载的最小方差无失真响应(Minimum variance distortionlessresponse with diagonal loading，MVDR-DL)波束形成器作为预处理器，在强干扰方向形成凹槽充分抑制强干扰信号，并计算MVDR-DL波束功率输出。构建适用于波束功率输出和MVDR-DL波束响应的线性关系的后验概率模型，采用变分贝叶斯推断(VariationalBayesian inference，VBI)进行DOA估计，将VBI推广至波束域中，避免了超参数的选取。
本发明的技术方案是：适用于强干扰环境下的波束域变分贝叶斯方位估计方法，其特征在于，包括以下步骤：
步骤1：假设KS个目标信号和KD个干扰信号分别从
Figure GDA0004089042410000021
和
Figure GDA0004089042410000022
方向上入射至M元的均匀线阵列，其中
Figure GDA0004089042410000023
表示第
Figure GDA0004089042410000024
个目标信号方位，
Figure GDA0004089042410000025
表示第
Figure GDA0004089042410000026
个干扰信号方位，目标信号和干扰信号之间互不相关；当阵列接收到信号后，将接收信号划分为N段，对每一段进行傅里叶变换后，宽带信号划分为L个子带；第l个子带上第n段阵列接收信号对应的傅里叶变换系数记为
Figure GDA0004089042410000027
采样协方差矩阵计算为
Figure GDA0004089042410000028
上标“H”为共轭转置运算；
步骤2：建立波束域模型，包括以下步骤：
步骤2.1：通过常规波束形成算法确定目标信号所在方位区域ΘS＝[ΘSL,ΘSR]，其中ΘSL和ΘSR分别为区域ΘS的左右边界；将该区域均匀划分为KB个网格
Figure GDA0004089042410000031
Figure GDA00040890424100000317
对于第l个子带，该区域上KB个MVDR-DL波束形成器组成的波束形成矩阵Wl为
Figure GDA0004089042410000032
其中
Figure GDA0004089042410000033
为指向φk的MVDR-DL波束形成器的加权量，al(φk)为第l个子带上指向φk的阵列流形，
Figure GDA0004089042410000034
Figure GDA0004089042410000035
为通过求解
Figure GDA0004089042410000036
的(M-KS-KD)个小特征值的平均值得到的噪声功率估计值，IM为M维的单位矩阵，上标“-1”表示矩阵求逆；
通过该矩阵对步骤1得到的采样协方差矩阵进行滤波，得到波束域上协方差矩阵
Figure GDA0004089042410000037
为
Figure GDA0004089042410000038
其中Pl S和Pl D表示第l个子带上的目标信号和干扰信号的协方差矩阵，σl表示第l个子带上的噪声功率，
Figure GDA0004089042410000039
和
Figure GDA00040890424100000310
表示第l个子带上目标信号和干扰信号的阵列流形矩阵，El为第l个子带上的扰动误差；
步骤2.2：对步骤2.1得到的协方差矩阵进行按列向量化运算，得到
Figure GDA00040890424100000311
其中
Figure GDA00040890424100000312
和
Figure GDA00040890424100000313
分别表示目标信号和干扰信号的功率向量，
Figure GDA00040890424100000314
和
Figure GDA00040890424100000315
分别表示Wl HW和Wl HElW矩阵按列向量化的向量，⊙表示Khatri-Rao乘积，上标“*”为共轭运算；
步骤2.3：定义矩阵
Figure GDA00040890424100000316
对于第m行第n列的元素[J]mn，若n＝KB(m-1)+m，[J]mn＝1，否则[J]mn＝0；将步骤2.2得到结果与该矩阵相乘，因MVDR-DL在区域ΘS外对干扰的波束响应较低，因此波束功率输出的线性关系式表示为：
Figure GDA0004089042410000041
其中
Figure GDA0004089042410000042
步骤2.4：将区域ΘS均匀划分为KG个网格，网格点
Figure GDA0004089042410000043
组成的向量记为
Figure GDA0004089042410000044
基于该离散网格，步骤2.3中的公式重新表示为
Figure GDA0004089042410000045
式中
Figure GDA0004089042410000046
Figure GDA0004089042410000047
为该网格上的阵列流形矩阵，
Figure GDA0004089042410000048
表示
Figure GDA0004089042410000049
的第h列；pl是一个稀疏向量，当
Figure GDA00040890424100000410
pl第m个元素等于
Figure GDA00040890424100000411
第n个元素，否则为0；
步骤3：在贝叶斯框架下进行迭代计算，最终输出DOA估计值，包括以下子步骤：
步骤3.1：构建后验概率密度函数模型，表示为：
Figure GDA00040890424100000412
式中
Figure GDA00040890424100000413
为联合概率密度函数，
Figure GDA00040890424100000414
N(·)表示实高斯分布，
Figure GDA00040890424100000415
表示Hadamard积；p(pl；γ)＝N(0,Γ-1)，
Figure GDA00040890424100000416
为信号稀疏参数，上标“T”为转置运算，Γ＝diag(γ)为以γ中元素为对角元素的对角矩阵；
Figure GDA00040890424100000417
Figure GDA00040890424100000418
为噪声功率的方差；
Figure GDA00040890424100000419
为边缘概率密度函数；
构建好后，给定参数迭代初始值：
信号稀疏参数初始值
Figure GDA00040890424100000420
其中
Figure GDA00040890424100000421
表示元素为1的KG×1维向量；
噪声功率后验均值初始值
Figure GDA00040890424100000422
噪声功率方差
Figure GDA00040890424100000423
步骤3.2：进行迭代更新，分别完成对信号功率后验协方差矩阵Σl，噪声功率后验方差
Figure GDA00040890424100000424
信号功率后验均值向量μl，信号稀疏参数γ，噪声功率后验均值
Figure GDA00040890424100000425
以及噪声功率方差
Figure GDA0004089042410000051
更新；
在第i次迭代中，第l个子带上信号功率后验协方差矩阵更新：
Figure GDA0004089042410000052
其中上标(i)表示第i次迭代，Γ(i-1)＝diag(γ(i-1))；
第l个子带上噪声功率后验方差更新：
Figure GDA0004089042410000053
第l个子带上信号功率后验均值向量更新：
Figure GDA0004089042410000054
信号稀疏参数更新：
Figure GDA0004089042410000055
其中
Figure GDA0004089042410000056
表示γ(i)第p个元素，
Figure GDA0004089042410000057
为
Figure GDA0004089042410000058
第p个元素，
Figure GDA0004089042410000059
为
Figure GDA00040890424100000510
第p行第p列的元素；
第l个子带上噪声功率后验均值更新：
Figure GDA00040890424100000511
第l个子带上噪声功率方差更新：
Figure GDA00040890424100000512
更新后的参数值若满足迭代终止条件，所估计的功率谱为
Figure GDA00040890424100000513
该功率谱中峰值对应的方位即为目标信号的DOA估计值；
若不满足迭代终止条件，则继续进行更新，直至满足迭代终止条件。
本发明进一步的技术方案是：所述步骤3.2中的迭代终止条件为：当迭代满足||γ(i)-γ(i-1)||2/||γ(i-1)||2≤10-3，其中||·||2表示l2范数，或者迭代次数大于Itermax＝1000时，迭代终止。
发明效果
本发明的技术效果在于：本发明采用MVDR-DL代替CBF对阵列接收信号进行滤波，充分抑制强干扰信号，避免干扰剩余量对后续DOA估计的影响；计算波束功率输出，构建了适用于MVDR-DL波束功率输出与波束响应的线性关系的后验概率模型，采用VBI自动迭代估计各个参数，将VBI推广至波束域中，避免了选取合适超参数的难题，同时也增强了VBI在强干扰环境下对目标信号方位估计的稳健性。
附图说明
图1基于波束功率输出的变分贝叶斯推断的总体流程框图
图2构建的贝叶斯概率图模型
图3VBI-BPO方法迭代流程
图4采用常规波束形成算法估计目标所在区域结果
图5基于CBF功率输出的lp范数算法的方位估计结果
图6VBI-BPO方法的方位估计结果
具体实施方式
在本发明的描述中，需要理解的是，术语“中心”、“纵向”、“横向”、“长度”、“宽度”、“厚度”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”、“顺时针”、“逆时针”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。
参见图1-图6，基于波束功率输出的变分贝叶斯推断的总体流程框图总结如图1，本发明解决其技术问题所采用的技术方案包括以下步骤：
步骤1：计算阵列接收信号采样协方差矩阵
假设KS个目标信号和KD个干扰信号分别从
Figure GDA0004089042410000061
和
Figure GDA0004089042410000062
方向上入射至M元的均匀线阵列，其中
Figure GDA0004089042410000063
表示第
Figure GDA0004089042410000064
个目标信号方位，
Figure GDA0004089042410000065
表示第
Figure GDA0004089042410000066
个干扰信号方位，目标信号和干扰信号之间互不相关。当阵列接收到信号后，将接收信号划分为N段，对每一段信号进行傅里叶变换后，宽带信号划分为L个子带。第l个子带上第n段阵列接收信号对应的傅里叶变换系数记为
Figure GDA0004089042410000067
采样协方差矩阵计算为
Figure GDA0004089042410000071
上标“H”为共轭转置运算。
步骤2：建立波束域模型，包括以下步骤：
步骤2.1：通过常规波束形成算法确定目标信号所在方位区域ΘS＝[ΘSL,ΘSR]，其中ΘSL和ΘSR分别为区域ΘS的左右边界。将该区域均匀划分为KB个网格
Figure GDA0004089042410000072
Figure GDA0004089042410000073
对于第l个子带，该区域上KB个MVDR-DL波束形成器组成的波束形成矩阵Wl为
Figure GDA0004089042410000074
其中
Figure GDA0004089042410000075
为指向φk的MVDR-DL波束形成器的加权量，al(φk)为第l个子带上指向φk的阵列流形，
Figure GDA0004089042410000076
Figure GDA0004089042410000077
为通过求解
Figure GDA0004089042410000078
的(M-KS-KD)个小特征值的平均值得到的噪声功率估计值，IM为M维的单位矩阵，上标“-1”表示矩阵求逆。
通过该矩阵对步骤1得到的采样协方差矩阵进行滤波，得到波束域上协方差矩阵
Figure GDA0004089042410000079
为
Figure GDA00040890424100000710
其中Pl S和Pl D表示第l个子带上的目标信号和干扰信号的协方差矩阵，σl表示第l个子带上的噪声功率，
Figure GDA00040890424100000711
和
Figure GDA00040890424100000712
表示第l个子带上目标信号和干扰信号的阵列流形矩阵，El为第l个子带上的扰动误差。
步骤2.2：对步骤2.1得到的协方差矩阵进行按列向量化运算，得到
Figure GDA00040890424100000713
其中
Figure GDA00040890424100000714
和
Figure GDA00040890424100000715
分别表示目标信号和干扰信号的功率向量，
Figure GDA00040890424100000716
和
Figure GDA00040890424100000717
分别表示Wl HW和Wl HElW矩阵按列向量化的向量，⊙表示Khatri-Rao乘积，上标“*”为共轭运算。
步骤2.3：定义一个矩阵
Figure GDA0004089042410000081
其第m行第n列的元素[J]mn为
Figure GDA0004089042410000082
将步骤2.2得到结果与该矩阵相乘，因MVDR-DL在区域ΘS外对干扰的波束响应较低，因此波束功率输出的线性关系式表示为：
Figure GDA0004089042410000083
其中
Figure GDA0004089042410000084
步骤2.4：将区域ΘS均匀划分为KG个网格，网格点
Figure GDA0004089042410000085
组成的向量记为
Figure GDA0004089042410000086
基于该离散网格，步骤2.3中的公式重新表示为
Figure GDA0004089042410000087
式中
Figure GDA0004089042410000088
Figure GDA00040890424100000820
为该网格上的阵列流形矩阵，
Figure GDA0004089042410000089
表示
Figure GDA00040890424100000810
的第h列。pl是一个稀疏向量，当
Figure GDA00040890424100000811
pl第m个元素等于
Figure GDA00040890424100000812
第n个元素，否则为0。
步骤3：在贝叶斯框架下进行迭代计算，最终输出DOA估计值，包括以下子步骤：
步骤3.1：如图2所示的构建的贝叶斯概率图模型，所构建后验概率密度函数模型表示为：
Figure GDA00040890424100000813
式中
Figure GDA00040890424100000814
为联合概率密度函数，
Figure GDA00040890424100000815
N(·)表示实高斯分布，○表示Hadamard积；p(pl；γ)＝N(0,Γ-1)，
Figure GDA00040890424100000816
为信号稀疏参数，上标“T”为转置运算，Γ＝diag(γ)为以γ中元素为对角元素的对角矩阵；
Figure GDA00040890424100000817
Figure GDA00040890424100000818
为噪声功率的方差；
Figure GDA00040890424100000819
为边缘概率密度函数。
构建好后，给定参数迭代初始值：
信号稀疏参数初始值
Figure GDA0004089042410000091
其中
Figure GDA0004089042410000092
表示元素为1的KG×1维向量；
噪声功率后验均值初始值
Figure GDA0004089042410000093
噪声功率方差
Figure GDA0004089042410000094
步骤3.2：进行迭代更新。分别完成对信号功率后验协方差矩阵Σl，噪声功率后验方差
Figure GDA0004089042410000095
信号功率后验均值向量μl，信号稀疏参数γ，噪声功率后验均值
Figure GDA0004089042410000096
以及噪声功率方差
Figure GDA0004089042410000097
的更新。图3总结了VBI-BPO方法迭代流程。
在第i次迭代中，第l个子带上信号功率后验协方差矩阵更新：
Figure GDA0004089042410000098
其中上标(i)表示第i次迭代，Γ(i-1)＝diag(γ(i-1))；
第l个子带上噪声功率后验方差更新：
Figure GDA0004089042410000099
第l个子带上信号功率后验均值向量更新：
Figure GDA00040890424100000910
信号稀疏参数更新：
Figure GDA00040890424100000911
其中
Figure GDA00040890424100000912
表示γ(i)第p个元素，
Figure GDA00040890424100000913
为
Figure GDA00040890424100000919
第p个元素，
Figure GDA00040890424100000914
为
Figure GDA00040890424100000915
第p行第p列的元素；
第l个子带上噪声功率后验均值更新：
Figure GDA00040890424100000916
第l个子带上噪声功率方差更新：
Figure GDA00040890424100000917
更新后的参数值若满足||γ(i)-γ(i-1)||2/||γ(i-1)||2≤10-3，其中||·||2表示l2范数，或者迭代次数大于Itermax＝1000时，迭代终止，所估计的功率谱为
Figure GDA00040890424100000918
该功率谱中峰值对应的方位即为目标信号的DOA估计值；若不满足迭代终止条件，则继续进行更新，直至满足迭代终止条件。
为了验证利用本发明给出的方法在强干扰环境下对目标信号方位估计的有效性，设计仿真实验如下：假设接收阵列为M＝32元均匀线列阵，阵元间距为4m。两个远处航船从-10°和-7°方向上入射至阵列，功率为0dB，视为目标信号，即KS＝2；一个近处水面舰从10°方向上入射至阵列，功率为20dB，大于远处航船，视为干扰信号，即KD＝1。所考虑的信号频段为[90，180]Hz。噪声为高斯白噪声，其功率在所考虑的信号频段上为0dB。
将接收信号均匀划分为N＝50段，在每一段上进行傅里叶变换，将信号划分至L＝46个子带。图4为采用常规波束形成算法对目标信号所在区域估计的结果，虚线为目标信号所在区域范围ΘS的边界ΘSL＝-14°和ΘSR＝-2°。故ΘS为[-14,-2]°。
将该区域以2°为间隔均匀划分为7个网格点，得到MVDR-DL波束指向角，即KB＝7。为了进行DOA估计，将该区域以1°为间隔均匀划分为13个网格点，即KG＝13，在该网格上进行DOA估计。图5和图6分别为基于CBF功率输出的lp范数方法和VBI-BPO方法的方位估计结果，图中的虚线为目标信号的真实方位。可以看出，本发明给出的方法(图6)可以很好地估计出两个目标信号，而基于CBF功率输出的lp范数方法在所考虑的仿真条件下完全失效，由此证明了本发明给出的方法能够有效解决强干扰环境下对目标信号的方位估计问题。

Claims (2)

1.适用于强干扰环境下的波束域变分贝叶斯方位估计方法，其特征在于，包括以下步骤：
步骤1：假设KS个目标信号和KD个干扰信号分别从
Figure FDA0004089042400000011
和
Figure FDA0004089042400000012
方向上入射至M元的均匀线阵列，其中
Figure FDA0004089042400000013
表示第
Figure FDA0004089042400000014
个目标信号方位，
Figure FDA0004089042400000015
表示第
Figure FDA0004089042400000016
个干扰信号方位，目标信号和干扰信号之间互不相关；当阵列接收到信号后，将接收信号划分为N段，对每一段进行傅里叶变换后，宽带信号划分为L个子带；第l个子带上第n段阵列接收信号对应的傅里叶变换系数记为
Figure FDA0004089042400000017
采样协方差矩阵计算为
Figure FDA0004089042400000018
上标“H”为共轭转置运算；
步骤2：建立波束域模型，包括以下步骤：
步骤2.1：通过常规波束形成算法确定目标信号所在方位区域ΘS＝[ΘSL,ΘSR]，其中ΘSL和ΘSR分别为区域ΘS的左右边界；将该区域均匀划分为KB个网格
Figure FDA00040890424000000114
Figure FDA00040890424000000115
对于第l个子带，该区域上KB个MVDR-DL波束形成器组成的波束形成矩阵Wl为
Figure FDA00040890424000000116
其中
Figure FDA0004089042400000019
为指向φk的MVDR-DL波束形成器的加权量，al(φk)为第l个子带上指向φk的阵列流形，
Figure FDA00040890424000000110
为通过求解
Figure FDA00040890424000000111
的(M-KS-KD)个小特征值的平均值得到的噪声功率估计值，IM为M维的单位矩阵，上标“-1”表示矩阵求逆；
通过该矩阵对步骤1得到的采样协方差矩阵进行滤波，得到波束域上协方差矩阵
Figure FDA00040890424000000112
为
Figure FDA00040890424000000113
其中Pl S和Pl D表示第l个子带上的目标信号和干扰信号的协方差矩阵，σl表示第l个子带上的噪声功率，
Figure FDA0004089042400000021
和
Figure FDA0004089042400000022
表示第l个子带上目标信号和干扰信号的阵列流形矩阵，El为第l个子带上的扰动误差；
步骤2.2：对步骤2.1得到的协方差矩阵进行按列向量化运算，得到
Figure FDA0004089042400000023
其中
Figure FDA0004089042400000024
和
Figure FDA0004089042400000025
分别表示目标信号和干扰信号的功率向量，
Figure FDA0004089042400000026
和
Figure FDA0004089042400000027
分别表示Wl HW和Wl HElW矩阵按列向量化的向量，⊙表示Khatri-Rao乘积，上标“*”为共轭运算；
步骤2.3：定义矩阵
Figure FDA0004089042400000028
对于第m行第n列的元素[J]mn，若n＝KB(m-1)+m，[J]mn＝1，否则[J]mn＝0；将步骤2.2得到结果与该矩阵相乘，因MVDR-DL在区域ΘS外对干扰的波束响应较低，因此波束功率输出的线性关系式表示为：
Figure FDA0004089042400000029
其中
Figure FDA00040890424000000210
步骤2.4：将区域ΘS均匀划分为KG个网格，网格点
Figure FDA00040890424000000211
组成的向量记为
Figure FDA00040890424000000212
基于该离散网格，步骤2.3中的公式重新表示为
Figure FDA00040890424000000213
式中
Figure FDA00040890424000000214
Figure FDA00040890424000000215
为该网格上的阵列流形矩阵，
Figure FDA00040890424000000216
h＝1,...,KG表示
Figure FDA00040890424000000217
的第h列；pl是一个稀疏向量，当
Figure FDA00040890424000000218
pl第m个元素等于
Figure FDA00040890424000000219
第n个元素，否则为0；
步骤3：在贝叶斯框架下进行迭代计算，最终输出DOA估计值，包括以下子步骤：
步骤3.1：构建后验概率密度函数模型，表示为：
Figure FDA00040890424000000220
式中
Figure FDA00040890424000000221
为联合概率密度函数，
Figure FDA0004089042400000031
N(·)表示实高斯分布，表示Hadamard积；p(pl；γ)＝N(0,Γ-1)，
Figure FDA0004089042400000032
为信号稀疏参数，上标“T”为转置运算，Γdiag(γ)为以γ中元素为对角元素的对角矩阵；
Figure FDA0004089042400000033
为噪声功率的方差；
Figure FDA0004089042400000034
为边缘概率密度函数；
构建好后，给定参数迭代初始值：
信号稀疏参数初始值
Figure FDA0004089042400000035
其中
Figure FDA0004089042400000036
表示元素为1的KG×1维向量；
噪声功率后验均值初始值
Figure FDA0004089042400000037
噪声功率方差
Figure FDA0004089042400000038
步骤3.2：进行迭代更新，分别完成对信号功率后验协方差矩阵Σl，噪声功率后验方差
Figure FDA0004089042400000039
信号功率后验均值向量μl，信号稀疏参数γ，噪声功率后验均值
Figure FDA00040890424000000310
以及噪声功率方差
Figure FDA00040890424000000311
更新；
在第i次迭代中，第l个子带上信号功率后验协方差矩阵更新：
Figure FDA00040890424000000312
其中上标(i)表示第i次迭代，Γ(i-1)＝diag(γ(i-1))；
第l个子带上噪声功率后验方差更新：
Figure FDA00040890424000000313
第l个子带上信号功率后验均值向量更新：
Figure FDA00040890424000000314
信号稀疏参数更新：
Figure FDA00040890424000000315
其中
Figure FDA00040890424000000316
表示γ(i)第p个元素，
Figure FDA00040890424000000317
为
Figure FDA00040890424000000318
第p个元素，
Figure FDA00040890424000000319
为
Figure FDA00040890424000000320
第p行第p列的元素；第l个子带上噪声功率后验均值更新：
Figure FDA00040890424000000321
第l个子带上噪声功率方差更新：
Figure FDA00040890424000000322
更新后的参数值若满足迭代终止条件，所估计的功率谱为
Figure FDA00040890424000000323
该功率谱中峰值对应的方位即为目标信号的DOA估计值；
若不满足迭代终止条件，则继续进行更新，直至满足迭代终止条件。
2.如权利要求1所述的适用于强干扰环境下的波束域变分贝叶斯方位估计方法，其特征在于，所述步骤3.2中的迭代终止条件为：当迭代满足||γ(i)-γ(i-1)||2/||γ(i-1)||2≤10-3，其中||·||2表示l2范数，或者迭代次数大于Itermax＝1000时，迭代终止。
CN202110954450.4A 2021-08-19 2021-08-19 适用于强干扰环境下的波束域变分贝叶斯方位估计方法 Active CN113673158B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110954450.4A CN113673158B (zh) 2021-08-19 2021-08-19 适用于强干扰环境下的波束域变分贝叶斯方位估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110954450.4A CN113673158B (zh) 2021-08-19 2021-08-19 适用于强干扰环境下的波束域变分贝叶斯方位估计方法

Publications (2)

Publication Number Publication Date
CN113673158A CN113673158A (zh) 2021-11-19
CN113673158B true CN113673158B (zh) 2023-05-26

Family

ID=78543974

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110954450.4A Active CN113673158B (zh) 2021-08-19 2021-08-19 适用于强干扰环境下的波束域变分贝叶斯方位估计方法

Country Status (1)

Country Link
CN (1) CN113673158B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114401168B (zh) * 2021-12-17 2023-11-03 郑州中科集成电路与系统应用研究院 适用复杂强噪声环境下短波莫尔斯信号的语音增强方法
CN115130504B (zh) * 2022-06-22 2024-09-13 西北工业大学 基于稀疏贝叶斯学习的稳健波束形成方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109116337A (zh) * 2018-07-30 2019-01-01 西北工业大学 一种基于矩阵滤波的稀疏近似最小方差方位估计方法

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
AU2003903826A0 (en) * 2003-07-24 2003-08-07 University Of South Australia An ofdm receiver structure
KR101209908B1 (ko) * 2011-08-04 2012-12-11 광주과학기술원 희소 신호 전송 방법 및 장치, 그리고 희소 신호 복구 방법 및 장치
WO2014066360A1 (en) * 2012-10-22 2014-05-01 Saab-Sensis Corporation Sensor system and method for determining target location using sparsity-based processing
CN104749553B (zh) * 2015-04-10 2017-03-08 西安电子科技大学 基于快速稀疏贝叶斯学习的波达方向角估计方法
CN109407046A (zh) * 2018-09-10 2019-03-01 西北工业大学 一种基于变分贝叶斯推断的嵌套阵列波达方向角估计方法
CN110133578B (zh) * 2019-05-08 2023-02-28 西北工业大学 一种基于半圆柱体积阵的海底反射声线入射角度估计方法
CN110208735B (zh) * 2019-06-12 2022-11-11 西北工业大学 一种基于稀疏贝叶斯学习的相干信号doa估计方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109116337A (zh) * 2018-07-30 2019-01-01 西北工业大学 一种基于矩阵滤波的稀疏近似最小方差方位估计方法

Also Published As

Publication number Publication date
CN113673158A (zh) 2021-11-19

Similar Documents

Publication Publication Date Title
CN110275166B (zh) 基于admm的快速稀疏孔径isar自聚焦与成像方法
CN109490850B (zh) 主瓣干扰下宽带阵列自适应波束形成方法
CN110113085B (zh) 一种基于协方差矩阵重构的波束形成方法及系统
CN108375763B (zh) 一种应用于多声源环境的分频定位方法
CN113673158B (zh) 适用于强干扰环境下的波束域变分贝叶斯方位估计方法
CN109116337B (zh) 一种基于矩阵滤波的稀疏近似最小方差方位估计方法
CN113673419B (zh) 适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法
CN110045323B (zh) 一种基于矩阵填充的互质阵稳健自适应波束形成算法
CN113422630B (zh) 一种自适应聚焦宽带波束形成方法及系统
CN109298383B (zh) 一种基于变分贝叶斯推断的互质阵波达方向角估计方法
CN105445709B (zh) 一种稀布阵列近场无源定位幅相误差校正方法
CN105137399A (zh) 基于斜投影滤波的雷达自适应波束形成方法
CN112612006B (zh) 基于深度学习的机载雷达非均匀杂波抑制方法
CN106707250A (zh) 基于互耦校正的雷达阵列自适应波束形成方法
CN106842135B (zh) 基于干扰加噪声协方差矩阵重构的自适应波束形成方法
CN110535519A (zh) 一种基于空间平滑的稳健自适应波束形成方法
CN109541572B (zh) 一种基于线性环境噪声模型的子空间方位估计方法
CN114487985B (zh) 一种基于差-和信号的波束锐化方法及系统
CN114371441B (zh) 虚拟阵列波达方向估计方法、装置、产品及存储介质
CN105572631A (zh) 基于多波位联合处理的最大似然目标doa估计方法
CN114152918A (zh) 基于压缩感知的抗间歇式主瓣干扰方法
CN111880167A (zh) 一种基于先随机后优化的波达方向估计方法
CN113933779B (zh) 一种基于s变换的未知声源个数doa估计方法
WO2024182916A1 (en) Adaptating a microphone array to a target beamformer
CN117590363A (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