CN114755628A - 非均匀噪声下声矢量传感器阵列波达方向估计方法 - Google Patents

非均匀噪声下声矢量传感器阵列波达方向估计方法 Download PDF

Info

Publication number
CN114755628A
CN114755628A CN202210353033.9A CN202210353033A CN114755628A CN 114755628 A CN114755628 A CN 114755628A CN 202210353033 A CN202210353033 A CN 202210353033A CN 114755628 A CN114755628 A CN 114755628A
Authority
CN
China
Prior art keywords
covariance matrix
signal
sparse
noise
vector
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.)
Pending
Application number
CN202210353033.9A
Other languages
English (en)
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.)
Henan University of Technology
Original Assignee
Henan University of Technology
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 Henan University of Technology filed Critical Henan University of Technology
Priority to CN202210353033.9A priority Critical patent/CN114755628A/zh
Publication of CN114755628A publication Critical patent/CN114755628A/zh
Pending legal-status Critical Current

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/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/8003—Diversity systems specially adapted for direction finding
    • 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

Landscapes

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

Abstract

本发明涉及阵列信号处理技术领域，公开了一种非均匀噪声下声矢量传感器阵列波达方向估计方法，所述方法包括以下步骤：步骤1，接收声矢量传感器阵列输出的输出矢量X；步骤2，基于稀疏信号模型，根据所述输出矢量X计算信号协方差矩阵P和噪声协方差矩阵Q，利用所得的信号协方差矩阵P和噪声协方差矩阵Q，构建稀疏协方差矩阵R；步骤3，基于所述稀疏协方差矩阵R构造代价函数，使用代价函数估计稀疏信号功率，并得到信号功率矢量
Figure DDA0003581551110000012
步骤4，对信号功率矢量
Figure DDA0003581551110000011
进行谱峰搜索，所得谱峰对应的声源位置即为估计的目标方位角。该方法能够提高非均匀噪声下声矢量传感器阵列的方位估计的性能。

Description

非均匀噪声下声矢量传感器阵列波达方向估计方法
技术领域
本发明涉及阵列信号处理技术领域，具体涉及一种非均匀噪声下声矢量传感器阵列波达方向估计方法。
背景技术
声矢量传感器是一种新型传感器，可空间共点、同步地测量声场中的声压和振速信息。与传统声压传感器相比，声矢量传感器不仅可获取声场中的声压信息，还可获取声振速信息，通过利用额外的声振速信息，声矢量传感器具有更好的信号检测和方位估计性能。相较于单声矢量传感器，由其构成的声矢量传感器阵列具有高信号增益、强干扰抑制能力以及高空间分辨率等优点。现已广泛应用于声纳、雷达等探测领域。
估计水下声源的波达方向(Direction Of Arrival,DOA)是被动探测领域的一个热门研究方向。波束形成方法是一种比较经典的DOA估计方法，具有结构简单、计算量小的优点，但受到物理孔径的限制，分辨性能较差。子空间类算法，如多重信号分类(MultipleSignal Classification,MUSIC)算法及其改进算法，虽然可以实现高分辨估计，但在低信噪比、少快拍及相干源场景下，该类方法的DOA估计性能会进一步恶化。稀疏重构类方法利用入射信号相对于整个空间具有稀疏分布的特性，对信号所在角度空间进行网格划分，把方位估计问题转化为稀疏信号重构问题。该类算法弥补了传统方位估计方法在相干源、低信噪比及少快拍情况下方位估计性能恶化的缺陷。l1范数优化类方法是一种比较经典的稀疏信号重构类方法，在低信噪比和少快拍下具有比传统方法更高的估计精度和鲁棒性，然而l1范数优化类方法在复杂的水声环境中面临着正则化参数难以恰当选择的难题。迭代自适应方法(Iteration Adaptive Approach,IAA)是一种无需设置用户参数的稀疏方法，不仅具有稀疏重构类算法的优点，而且具有较小的计算量，已广泛应用于雷达***中。然而，上述方法是基于声压、声振速各分量的噪声功率相等的理想假设下，对声矢量传感器阵列的方位估计性能进行的研究。由于声矢量传感器接收到的噪声主要是环境噪声而不是传感器内部的噪声，加之声矢量传感器的振速轴具有指向性，使得振速轴方向测量的部分环境噪声被滤掉，导致振速分量的噪声功率不等于声压分量的噪声功率。在此情况下，上述方法的方位估计性能将会严重恶化。
为解决非均匀噪声情况下，声压传感器阵列方位估计性能恶化的问题，Bin Liao等人在论文“Iterative Methods for Subspace and DOA Estimation in NonuniformNoise”(IEEE Transactions on Signal Processing Volume 64,Issue 12.2016.PP3008-3020)中提出一种迭代子空间估计方法，采用最大似然(Maximum-likelihood,ML)和最小二乘(Least-squares,LS)的方法估计信号子空间和噪声协方差矩阵。尽管该方法能够应用于声矢量传感器阵列，解决非均匀噪声情况下声矢量传感器阵列方位估计性能恶化的难题，然而，这两种方法均采用迭代方式估计噪声协方差矩阵，导致该算法的计算复杂度较高，难以满足***对实时性的需求。为了解决此难题，Aifei Liu等人在论文“Augmentedsubspace MUSIC method for DOA estimation using acoustic vector sensor array”(IET Radar,Sonar&Navigation Volume 13,Issue 6.2019.PP 969-975)中，提出了增广子空间MUSIC(Augmented subspace MUSIC，ASMUSIC)方法，通过将虚拟源的数目扩展到信号子空间来提高声矢量传感器阵列在非均匀噪声情况下的方位估计精度。但是，该算法仅适用于大孔径阵列及声矢量传感器中声压和振速之间噪声功率不等的情况。在传感器之间相应不同时，声矢量传感器之间的声压和振速及振速和振速通道之间接收到的噪声也会各不相同，在此情况下，ASMUSIC方法的方位估计性能也会严重恶化。
发明内容
针对上述情况，本发明的目的是提供一种正则化加权稀疏协方差矩阵拟合(Weighted Sparse Covariance Matrix Fitting,WSCMF)方位估计方法，以提高非均匀噪声下声矢量传感器阵列的方位估计的精度。
为了实现上述目的，本发明通过以下技术方案实现的：
非均匀噪声下声矢量传感器阵列波达方向估计方法，所述方法包括以下步骤：
步骤1，接收声矢量传感器阵列输出的输出矢量X；
步骤2，基于稀疏信号模型，根据所述输出矢量X计算信号协方差矩阵P和噪声协方差矩阵Q，利用所得的信号协方差矩阵P和噪声协方差矩阵Q，构建稀疏协方差矩阵R；
步骤3，基于所述稀疏协方差矩阵R构造代价函数，使用代价函数估计稀疏信号功率，并得到信号功率矢量
Figure BDA0003581551090000041
步骤4，对信号功率矢量
Figure BDA0003581551090000042
进行谱峰搜索，所得谱峰对应的声源位置即为估计的目标方位角。
进一步的，上述非均匀噪声下声矢量传感器阵列波达方向估计方法中，步骤2包括以下子步骤：
步骤21，构建非均匀噪声情况下声矢量传感器阵列的稀疏信号模型：
假设K个波长为λ的远场等功率窄带信号，入射到由M个声矢量传感器组成的均匀线阵上，阵元间距为d，窄带信号沿着x轴和y轴方向的振速分量分别为ux和uy，且x轴和y轴具有正交性，在快拍数为L的情况下，声矢量传感器阵列的输出矢量X可表示为
X＝A(θ)S+W
式中，X＝[x(1)，x(2)，…，x(L)]表示声矢量传感器阵列的输出矢量，A(θ)＝[a(θ1)，a(θ2)，…，a(θK)]表示理想声矢量传感器阵列的阵列流形矩阵，θ＝[θ1，θ2，…，θK]表示目标的方位角，
Figure BDA0003581551090000043
Figure BDA0003581551090000044
表示第k个目标的流形矢量，
Figure BDA0003581551090000045
表示笛卡尔乘积，h(θk)＝[1 cosθk sinθk]T表示第k个目标入射到声矢量传感器上的阵列流形矢量，
Figure BDA0003581551090000046
Figure BDA0003581551090000047
表示第k个目标入射在第m个声矢量传感器上的声压响应系数，W＝[w(1)，w(2)，…，w(L)]表示噪声矢量，S＝[s(1)，s(2)，…，s(L)]表示输入信号波形矢量；
假设角度空间(-180°～180°)被均匀地划分为N个离散的网格，且满足离散网格数N远大于声源数K，则所有声源的来波方向可表示为
Figure BDA0003581551090000051
声矢量传感器阵列的稀疏信号模型可表示为
Figure BDA0003581551090000052
式中，
Figure BDA0003581551090000053
表示稀疏信号模型下声矢量传感器阵列的阵列流形矩阵，
Figure BDA0003581551090000054
是一个行稀疏矩阵，并且它的非零行表征着真实信号的位置；
步骤22，计算信号协方差矩阵P和噪声协方差矩阵Q：
信号协方差矩阵P可表示为
Figure BDA0003581551090000055
式中，E{·}表示期望运算，diag{·}表示对角矩阵操作，
Figure BDA0003581551090000056
表示第n个网格对应信号的功率；
噪声协方差矩阵Q可表示为
Q＝diag{q1，…，qM}
式中，
Figure BDA0003581551090000057
假设各入射信号之间互不相关，则稀疏信号模型下声矢量传感器阵列的协方差矩阵一般可表示为
Figure BDA0003581551090000058
一般的，若背景噪声为高斯白噪声，则噪声协方差矩阵Q＝σ2I3M，σ2表示均匀噪声环境下，声矢量传感器阵列每个通道输出的噪声功率，I3M表示3M×3M的单位矩阵。但在实际应用中，由于声矢量传感器接收的噪声主要是环境噪声而不是传感器内部的噪声，加之声矢量传感器的振速轴具有指向性，使得振速轴方向测量的部分环境噪声被滤掉，导致声振速分量的噪声功率并不等于声压分量的噪声功率。此外，考虑到传感器之间响应的不同，使得声矢量传感器之间声压通道和声振速通道之间接收到的噪声功率各不相同，即
Figure BDA0003581551090000061
为了解决上述难题，本发明定义了一个新的声矢量传感器阵列的阵列流形矩阵，重新构建了稀疏协方差矩阵R，使其能够同时估计每个网格上的信号功率和声矢量传感器阵列中每个通道输出的噪声功率。需要说明的是，所述声矢量传感器阵列的每个通道指的是每个声矢量传感器阵元的声压通道、以及每个声矢量传感器阵元分别沿着x轴和y轴方向的声振速分量对应的声振速通道。
步骤23，构建稀疏协方差矩阵R：
为了估计稀疏信号功率和声矢量传感器阵列中每个通道输出的噪声功率，定义一个新的声矢量传感器阵列的阵列流形矩阵
Figure BDA0003581551090000062
其中I3M表示3M×3M的单位矩阵，设定各入射信号之间互不相关，则稀疏协方差矩阵R可表示为
R＝DΞDH
其中，
Figure BDA0003581551090000063
Figure BDA0003581551090000071
式中，Ξ对角线元素的前N项表示信号功率，第N+1至N+3M表示声矢量传感器阵列中每个通道输出的噪声功率。
于是，噪声协方差矩阵Q可以进一步表示为
Figure BDA0003581551090000072
式中，Dn表示D的第n列数据，DN+l＝ul，l＝1，...，3M，ul表示I3M的第l列。
进一步的，上述非均匀噪声下声矢量传感器阵列波达方向估计方法中，步骤3包括以下子步骤：
步骤31，根据稀疏协方差矩阵R构建干扰加噪声协方差矩阵Bn：
根据步骤23构建的稀疏协方差矩阵R，第n个网格对应的信号协方差矩阵可表示为
Figure BDA0003581551090000073
则相应于第n个网格上信号的干扰加噪声协方差矩阵为Bn＝R-Rn＝DΞDH-Rn
步骤32，基于稀疏协方差矩阵拟合准则，构造如下代价函数：
Figure BDA0003581551090000081
式中，第一项
Figure BDA0003581551090000082
表示协方差矩阵拟合项，第二项
Figure BDA0003581551090000083
表示稀疏信号补偿项，||·||F表示取F范数运算，
Figure BDA0003581551090000084
表示样本协方差矩阵，
Figure BDA0003581551090000085
q表示用于补偿信号的用户参数，λs表示平衡信号稀疏性与拟合误差关系的正则化参数，且λs＞0；
步骤33，使用迭代法求解稀疏信号功率，并得到信号功率矢量
Figure BDA0003581551090000086
进一步的，上述非均匀噪声下声矢量传感器阵列波达方向估计方法中，步骤33包括以下子步骤：
步骤331，为了降低计算复杂度和噪声对上述代价函数中协方差矩阵拟合项效果的影响，对样本协方差矩阵
Figure BDA0003581551090000087
进行奇异值分解，得
Figure BDA0003581551090000088
式中，
Figure BDA0003581551090000089
表示信号子空间，
Figure BDA00035815510900000810
表示噪声子空间；
则代价函数可以表示为
Figure BDA00035815510900000811
步骤332，使用迭代法求解稀疏信号功率，设定第j+1次迭代得到唯一解，可得稀疏信号功率
Figure BDA00035815510900000812
式中，an为
Figure BDA00035815510900000813
的简写形式，
Figure BDA00035815510900000814
Figure BDA00035815510900000815
Figure BDA00035815510900000816
Figure BDA0003581551090000091
具体的，使用迭代法求解稀疏信号功率的过程如下：
从代价函数的表达式可以看出，用户参数q的引入，使得代价函数是关于待求变量
Figure BDA0003581551090000092
的非线性函数。所以，从代价函数中直接求解
Figure BDA0003581551090000093
是十分困难的。为解决此问题，假设f(x)＝xq是关于变量x的凹函数。当变量x＞0时，对于x0＞0，可得
f(x)≤f(x0)+f′(x0)(x-x0) (1)
令
Figure BDA0003581551090000094
上标j表示迭代次数，代入式(1)可得
Figure BDA0003581551090000095
对式(2)两边分别乘以λs，再加上协方差矩阵拟合项
Figure BDA0003581551090000096
可进一步转化为
Figure BDA0003581551090000097
式中，
Figure BDA0003581551090000098
因此，代价函数关于
Figure BDA0003581551090000099
最小化的问题最终等价于
Figure BDA00035815510900000910
在式(4)中，
Figure BDA00035815510900000911
是关于
Figure BDA00035815510900000912
的函数，由
Figure BDA00035815510900000913
可知，
Figure BDA00035815510900000914
也是关于
Figure BDA00035815510900000915
的函数。因此，可以认为在式(4)中，
Figure BDA00035815510900000916
是未知变量。假设在式(4)中，
Figure BDA0003581551090000101
是在第j次迭代中已估计，则式(5)最小化
Figure BDA0003581551090000102
可以转化为最小化
Figure BDA0003581551090000103
式中，
Figure BDA0003581551090000104
由式(7)易知，
Figure BDA0003581551090000105
是关于
Figure BDA0003581551090000106
的线性函数。因此，可以直接从
Figure BDA0003581551090000107
中求解稀疏信号功率和噪声功率。
根据Frobenius范数的性质，对式(7)进一步化简得
Figure BDA0003581551090000108
式中，Tr{·}表示矩阵求迹运算。对
Figure BDA0003581551090000109
关于
Figure BDA00035815510900001010
求一阶导数，可得
Figure BDA00035815510900001011
对
Figure BDA00035815510900001012
关于
Figure BDA00035815510900001013
求二阶导数，可得
Figure BDA00035815510900001014
由式(10)可得，
Figure BDA00035815510900001015
恒成立。因此，
Figure BDA00035815510900001016
有唯一解。令
Figure BDA0003581551090000111
可得函数
Figure BDA0003581551090000112
在第j+1次迭代中的唯一解为
Figure BDA0003581551090000113
从式(11)中可以看出，在计算N个离散网格上对应的稀疏信号功率时，需要计算N次干扰加噪声协方差矩阵Bn的逆，极大地增加了计算量。为了减少计算量，根据Woodbury公式，可得
Figure BDA0003581551090000114
将式(12)代入式(11)，可得
Figure BDA0003581551090000115
式中，
Figure BDA0003581551090000116
Figure BDA0003581551090000117
Figure BDA0003581551090000118
根据D的定义可得稀疏信号功率
Figure BDA0003581551090000119
式中，an为
Figure BDA00035815510900001110
的简写形式。
则噪声功率
Figure BDA00035815510900001111
式中，
Figure BDA0003581551090000121
Figure BDA0003581551090000122
因此，声压通道和声振速通道的噪声功率为
Figure BDA0003581551090000123
Figure BDA0003581551090000124
Figure BDA0003581551090000125
步骤333，根据步骤332求得的稀疏信号功率，得到N个网格对应信号的信号功率矢量为
Figure BDA0003581551090000126
与现有技术相比，本发明的有益效果为：
1、本发明定义了新的声矢量传感器阵列的阵列流形矩阵，根据新定义的声矢量传感器阵列的阵列流形矩阵，重构了包含稀疏信号功率和噪声功率的稀疏协方差矩阵，并基于稀疏协方差矩阵拟合准则和稀疏信号补偿机理，构建了关于稀疏信号功率和噪声功率的代价函数，然后，采用正则化加权协方差矩阵拟合方法估计噪声功率与信号功率。在每次迭代中，通过本发明提供的方法估计的稀疏信号功率和噪声功率对稀疏协方差矩阵进行重构，使其在下一次迭代中估计的稀疏信号功率和噪声功率更精确，待迭代终止时对稀疏信号功率进行谱峰搜索可实现更加精确的方位估计。
相较于现有的波达方向估计方法中只对信号功率进行估计的方法，本发明通过重构包含稀疏信号功率和噪声功率的稀疏协方差矩阵，同时对信号功率和噪声功率进行估计，使得利用精确的噪声功率所估计的信号功率也比现有估计方法更加精确。。
2、相较于现有的波达方向估计方法中对声矢量传感器阵列的输出矢量的白化处理，本发明提供的波达方向估计方法具有更少的计算量和更高的估计精度。在实际应用中，本发明所提供的方法可以提高非均匀噪声下声矢量传感器阵列的方位估计精度。
附图说明
图1为本发明提出的WSCMF方法的流程图；
图2为声矢量传感器阵列的均匀线阵模型；
图3为本发明提出的WSCMF方法和MUSIC方法、IAA方法、IMLSE方法及ASMUSIC方法的归一化空间谱图；
图4为本发明提出的WSCMF方法和MUSIC方法、IAA方法、IMLSE方法及ASMUSIC方法的均方根误差和角度间隔的关系曲线图；
图5为本发明提出的WSCMF方法和MUSIC方法、IAA方法、IMLSE方法及ASMUSIC方法的均方根误差和信噪比的关系曲线图；
图6为本发明提出的WSCMF方法和MUSIC方法、IAA方法、IMLSE方法及ASMUSIC方法的均方根误差和最差噪声功率比的关系曲线图。
具体实施方式
下面将结合本实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，均属于本发明的保护范围。
实施例
如图1-6所示，非均匀噪声下声矢量传感器阵列波达方向估计方法，所述方法包括以下步骤：
步骤1，接收声矢量传感器阵列输出的输出矢量X；
步骤2，基于稀疏信号模型，根据所述输出矢量X计算信号协方差矩阵P和噪声协方差矩阵Q，利用所得的信号协方差矩阵P和噪声协方差矩阵Q，构建稀疏协方差矩阵R；
具体的，步骤2包括以下子步骤：
步骤21，构建非均匀噪声情况下声矢量传感器阵列的稀疏信号模型：
假设K个波长为λ的远场等功率窄带信号，入射到由M个声矢量传感器组成的均匀线阵上，声矢量传感器阵列的配置如图2所示，其中，阵元间距为d，窄带信号沿着x轴和y轴方向的振速分量分别为ux和uy，且x轴和y轴具有正交性，在快拍数为L的情况下，声矢量传感器阵列的输出矢量X可表示为
X＝A(θ)S+W
式中，X＝[x(1)，x(2)，…，x(L)]表示声矢量传感器阵列的输出矢量，A(θ)＝[a(θ1)，a(θ2)，…，a(θK)]表示理想声矢量传感器阵列的阵列流形矩阵，θ＝[θ1，θ2，…，θK]表示目标的方位角，
Figure BDA0003581551090000141
Figure BDA0003581551090000142
表示第k个目标的流形矢量，
Figure BDA0003581551090000143
表示笛卡尔乘积，h(θk)＝[1 cosθk sinθk]T表示第k个目标入射到声矢量传感器上的阵列流形矢量，
Figure BDA0003581551090000151
Figure BDA0003581551090000152
表示第k个目标入射在第m个声矢量传感器上的声压响应系数，W＝[w(1)，w(2)，…，w(L)]表示噪声矢量，S＝[s(1)，s(2)，…，s(L)]表示输入信号波形矢量；
假设角度空间(-180°～180°)被均匀地划分为N个离散的网格，且满足离散网格数N远大于声源数K，则所有声源的来波方向可表示为
Figure BDA0003581551090000153
声矢量传感器阵列的稀疏信号模型可表示为
Figure BDA0003581551090000154
式中，
Figure BDA0003581551090000155
表示稀疏信号模型下声矢量传感器阵列的阵列流形矩阵，
Figure BDA0003581551090000156
是一个行稀疏矩阵，并且它的非零行表征着真实信号的位置；
步骤22，计算信号协方差矩阵P和噪声协方差矩阵Q：
信号协方差矩阵P可表示为
Figure BDA0003581551090000157
式中，E{·}表示期望运算，diag{·}表示对角矩阵操作，
Figure BDA0003581551090000158
表示第n个网格对应信号的功率；
噪声协方差矩阵Q可表示为
Q＝diag{q1，…，qM}
式中，
Figure BDA0003581551090000159
步骤23，构建稀疏协方差矩阵R：
定义一个新的声矢量传感器阵列的阵列流形矩阵
Figure BDA00035815510900001510
其中I3M表示3M×3M的单位矩阵，则稀疏协方差矩阵R可表示为
R＝DΞDH
其中，
Figure BDA0003581551090000161
式中，Ξ的对角线元素的前N项表示信号功率，第N+1至N+3M表示声矢量传感器阵列中每个通道输出的噪声功率。于是，噪声协方差矩阵Q可以进一步表示为
Figure BDA0003581551090000162
式中，Dn表示D的第n列数据，DN+l＝ul，l＝1，...，3M，ul表示I3M的第l列。
步骤3，基于所述稀疏协方差矩阵R构造代价函数，使用代价函数估计稀疏信号功率，并得到信号功率矢量
Figure BDA0003581551090000163
具体的，步骤3包括以下子步骤：
步骤31，根据稀疏协方差矩阵R构建干扰加噪声协方差矩阵Bn：
根据步骤23构建的稀疏协方差矩阵R，第n个网格对应的信号协方差矩阵可表示为
Figure BDA0003581551090000171
则相应于第n个网格上信号的干扰加噪声协方差矩阵为Bn＝R-Rn＝DΞDH-Rn
步骤32，基于稀疏协方差矩阵拟合准则，构造如下代价函数：
Figure BDA0003581551090000172
式中，||·||F表示取F范数运算，
Figure BDA0003581551090000173
表示样本协方差矩阵，q表示用于补偿信号的用户参数，λs表示平衡信号稀疏性与拟合误差关系的正则化参数，且λs＞0；
步骤33，使用迭代法求解稀疏信号功率，并得到信号功率矢量
Figure BDA0003581551090000174
具体的，步骤33包括以下子步骤：
步骤331，对样本协方差矩阵进行奇异值分解，得
Figure BDA0003581551090000175
式中，
Figure BDA0003581551090000176
表示信号子空间，
Figure BDA0003581551090000177
表示噪声子空间；
则代价函数可以表示为
Figure BDA0003581551090000178
步骤332，使用迭代法求解稀疏信号功率，设定第j+1次迭代得到唯一解，可得稀疏信号功率
Figure BDA0003581551090000179
式中，an为
Figure BDA00035815510900001710
的简写形式，
Figure BDA0003581551090000181
Figure BDA0003581551090000182
Figure BDA0003581551090000183
步骤333，根据步骤332求得的稀疏信号功率，得到N个网格对应信号的信号功率矢量为
Figure BDA0003581551090000184
步骤4，对信号功率矢量
Figure BDA0003581551090000185
进行谱峰搜索，所得谱峰对应的声源位置即为估计的目标方位角。
与现有非均匀噪声功率估计方法相比，本发明所设计的方法具有更少的计算量和更高的估计精度。同时，利用精确的噪声功率所估计的信号功率也比现有估计方法更加精确。
本发明的效果可以通过仿真实验与现有的MUSIC方法、IAA方法、迭代最大似然子空间估计(IterativeMLSubspaceEstimation，IMLSE)方法及ASMUSIC方法作比较和说明：
对比例1
从归一化空间谱上比较本发明提出的WSCMF方法与现有的MUSIC方法、IAA方法、IMLSE方法及ASMUSIC方法的分辨性能。
在本对比例中，设定声矢量传感器阵列为由M＝4个声矢量传感器组成的均匀线列阵，声矢量传感器阵列的背景噪声为非均匀噪声，噪声协方差矩阵为
Q＝diag{3，1.5，2，5，2，1.5，4，2，5，2，30，2，1.5}
因此，最差噪声功率比(Worst Noise Power Ratio，WNPR)为
Figure BDA0003581551090000191
式中，
Figure BDA0003581551090000192
表示Q中的最大值，
Figure BDA0003581551090000193
表示Q中的最小值。信噪比(SNR)定义为
Figure BDA0003581551090000194
式中，
Figure BDA0003581551090000195
表示信号功率。
在本对比例中，设定三个信噪比为10dB的远场窄带信号分别从方位角θ1＝-15°，θ2＝38°，θ3＝50°入射在声矢量传感器阵列上，声速为1500m/s，***采样频率为5kHz，信号的中心频率为700Hz，快拍数为300。在本发明提出的WSCMF方法中，迭代次数为30，收敛门限为10-3。
从归一化空间谱上比较本发明提出的WSCMF方法和现有的MUSIC方法、IAA方法、IMLSE方法及ASMUSIC方法的分辨性能，仿真实验结果如图3所示。
从图3所示的归一化空间谱上可以看出，MUSIC方法、IAA方法及IMLSE方法都无法分辨出位于θ2＝38°的目标方位。尽管ASMUSIC方法可以分辨出方位角位于θ2＝38°出的目标，但是估计的结果明显偏离真实的目标方位角度。而本发明提出的WSCMF方法能够分辨出三个目标，并且估计结果与真实目标的方位角度最接近，展示出本发明方法具有较好的分辨性能。
对比例2
以均方根误差(Root Mean Square Error，RMSE)作为方位估计性能的评判标准，比较本发明提出的WSCMF方法与现有的MUSIC方法、IAA方法、IMLSE方法及ASMUSIC方法的估计精度。
通常的，均方根误差越小表示估计越精确，均方根误差越大表示估计精度越低。均方根误差可以表示为
Figure BDA0003581551090000201
式中，θk为第k个目标的真实水平方位角度，
Figure BDA0003581551090000202
第u次蒙特卡洛实验中估计的第k个目标的水平方位角，U＝3000实验中进行的独立蒙特卡洛试验次数。
在本对比例中，分别从均方根误差和角度间隔的关系曲线、均方根误差和信噪比的关系曲线及均方根误差和最差噪声功率比的关系曲线研究几种方位估计方法的估计性能，具体设定仿真实验条件如下：
(1)信噪比为5dB的两个等功率远场窄带信号入射在声矢量传感器阵列上，两个入射信号的角度分别为θ1＝-14°，θ2＝-14°+Δ，Δ表示两个入射信号的方位角度间隔，从15°到55°以5°为增量进行变化，比较应用本发明提出的WSCMF方法与现有的MUSIC方法、IAA方法、IMLSE方法及ASMUSIC方法所得估计结果的均方根误差与角度间隔的关系，结果如图4所示；
(2)两个等功率的远场窄带信号分别从θ1＝-14°，θ2＝38°入射在声矢量传感器阵列上，信噪比从-15dB到15dB变化，比较应用本发明提出的WSCMF方法与现有的MUSIC方法、IAA方法、IMLSE方法及ASMUSIC方法所得估计结果的均方根误差与信噪比的关系，结果如图5所示；
(3)两个等功率的远场窄带信号分别从方位角θ1＝-14°，θ2＝38°入射在声矢量传感器阵列上，信噪比为5dB，最差噪声功率比从10到80以10为增量进行变化，比较应用本发明提出的WSCMF方法与现有的MUSIC方法、IAA方法、IMLSE方法及ASMUSIC方法所得估计结果的均方根误差与最差噪声功率比的关系，结果如图6所示。
从图4可以看出，IMLSE方法和ASMUSIC方法无法分辨出方位角度间隔分别小于45°和35°的两个目标。IAA方法无法分辨出角度间隔小于30°的两个目标。MUSIC方法无法分辨出角度间隔小于25°的两个目标。本发明提出的WSCMF方法可以分辨出角度间隔不小于20°的两个目标。尽管本发明方法在方位角度间隔小于25°时，均方根误差有所增加，但是相较于IMLSE方法、ASMUSIC方法、IAA方法和MUSIC方法，本发明方法的估计精度仍是最高的。此外，易看出，方位角度间隔大于50°时，本发明提出的WSCMF方法与IMLSE方法、ASMUSIC方法、IAA方法和MUSIC方法的方位估计精度都在可接受范围之内，但是本发明提出的WSCMF方法的方位估计精度是最高的。
从图5可以发现，随着信噪比的增加，MUSIC方法的估计性能提高得的较小。尽管IAA方法、IMLSE方法、ASMUSIC方法的性能有所提高，但在信噪比较低时，它们的方位估计精度都不够高。本发明提出的WSCMF方法采用正则化加权稀疏协方差矩阵拟合方法估计噪声功率和信号功率，并通过用户参数较好的平衡了信号在空域的稀疏性与拟合误差之间的关系，使得本发明提出的WSCMF方法在低信噪比时声矢量传感器阵列的方位估计精度得以提高。
从图6可以看出，未考虑非均匀噪声的MUSIC方法和IAA方法的方位估计估计性能随着最差噪声功率比的增加严重恶化。尽管随着最差噪声功率比的增大，IMLSE方法、ASMUSIC方法及本发明所提方法的均方根误差受影响较小，但是，相较于IMLSE方法和ASMUSIC方法，本发明所提方法的均方根误差最小，这表明在非均匀噪声场景下，本发明所提方法具较好的方位估计性能。
以上已经描述了本发明的实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的实施例。在不偏离所说明实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。

Claims (4)

1.非均匀噪声下声矢量传感器阵列波达方向估计方法，其特征在于，所述方法包括以下步骤：
步骤1，接收声矢量传感器阵列输出的输出矢量X；
步骤2，基于稀疏信号模型，根据所述输出矢量X计算信号协方差矩阵P和噪声协方差矩阵Q，利用所得的信号协方差矩阵P和噪声协方差矩阵Q，构建稀疏协方差矩阵R；
步骤3，基于所述稀疏协方差矩阵R构造代价函数，使用代价函数估计稀疏信号功率，并得到信号功率矢量
Figure FDA0003581551080000011
步骤4，对信号功率矢量
Figure FDA0003581551080000012
进行谱峰搜索，所得谱峰对应的声源位置即为估计的目标方位角。
2.根据权利要求1所述的非均匀噪声下声矢量传感器阵列波达方向估计方法，其特征在于，步骤2包括以下子步骤：
步骤21，构建非均匀噪声情况下声矢量传感器阵列的稀疏信号模型：
假设K个波长为λ的远场等功率窄带信号，入射到由M个声矢量传感器组成的均匀线阵上，阵元间距为d，窄带信号沿着x轴和y轴方向的振速分量分别为ux和uy，且x轴和y轴具有正交性，在快拍数为L的情况下，声矢量传感器阵列的输出矢量X可表示为
X＝A(θ)S+W
式中，X＝[x(1)，x(2)，…，x(L)]表示声矢量传感器阵列的输出矢量，A(θ)＝[a(θ1)，a(θ2)，…，a(θK)]表示理想声矢量传感器阵列的阵列流形矩阵，θ＝[θ1，θ2，…，θK]表示目标的方位角，
Figure FDA0003581551080000021
Figure FDA0003581551080000022
表示第k个目标的流形矢量，
Figure FDA0003581551080000023
表示笛卡尔乘积，h(θk)＝[1 cosθk sinθk]T表示第k个目标入射到声矢量传感器上的阵列流形矢量，
Figure FDA0003581551080000024
Figure FDA0003581551080000025
表示第k个目标入射在第m个声矢量传感器上的声压响应系数，W＝[w(1)，w(2)，…，w(L)]表示噪声矢量，S＝[s(1)，s(2)，…，s(L)]表示输入信号波形矢量；
假设角度空间(-180°～180°)被均匀地划分为N个离散的网格，且满足离散网格数N远大于声源数K，则所有声源的来波方向可表示为
Figure FDA0003581551080000026
声矢量传感器阵列的稀疏信号模型可表示为
Figure FDA0003581551080000027
式中，
Figure FDA0003581551080000028
表示稀疏信号模型下声矢量传感器阵列的阵列流形矩阵，
Figure FDA0003581551080000029
是一个行稀疏矩阵，并且它的非零行表征着真实信号的位置；
步骤22，计算信号协方差矩阵P和噪声协方差矩阵Q：
信号协方差矩阵P可表示为
Figure FDA00035815510800000210
式中，E{·}表示期望运算，diag{·}表示对角矩阵操作，
Figure FDA00035815510800000211
表示第n个网格对应信号的功率；
噪声协方差矩阵Q可表示为
Q＝diag{q1，…，qM}
式中，
Figure FDA00035815510800000212
步骤23，构建稀疏协方差矩阵R：
定义一个新的声矢量传感器阵列的阵列流形矩阵
Figure FDA0003581551080000031
其中I3M表示3M×3M的单位矩阵，则稀疏协方差矩阵R可表示为
R＝DΞDH
其中，
Figure FDA0003581551080000032
式中，Ξ的对角线元素的前N项表示信号功率，第N+1至N+3M表示声矢量传感器阵列中每个通道输出的噪声功率。于是，噪声协方差矩阵Q可以进一步表示为
Figure FDA0003581551080000033
式中，Dn表示D的第n列数据，DN+l＝ul，l＝1，...，3M，ul表示I3M的第l列。
3.根据权利要求2所述的非均匀噪声下声矢量传感器阵列波达方向估计方法，其特征在于，步骤3包括以下子步骤：
步骤31，根据稀疏协方差矩阵R构建干扰加噪声协方差矩阵Bn：
根据步骤23构建的稀疏协方差矩阵R，第n个网格对应的信号协方差矩阵可表示为
Figure FDA0003581551080000041
则相应于第n个网格上信号的干扰加噪声协方差矩阵为
Bn＝R-Rn＝DΞDH-Rn
步骤32，基于稀疏协方差矩阵拟合准则，构造如下代价函数：
Figure FDA0003581551080000042
式中，||·||F表示取F范数运算，
Figure FDA0003581551080000043
表示样本协方差矩阵，q表示用于补偿信号的用户参数，λs表示平衡信号稀疏性与拟合误差关系的正则化参数，且λs＞0；
步骤33，使用迭代法求解稀疏信号功率，并得到信号功率矢量
Figure FDA0003581551080000044
4.根据权利要求3所述的非均匀噪声下声矢量传感器阵列波达方向估计方法，其特征在于，步骤33包括以下子步骤：
步骤331，对样本协方差矩阵进行奇异值分解，得
Figure FDA0003581551080000045
式中，
Figure FDA0003581551080000046
表示信号子空间，
Figure FDA0003581551080000047
表示噪声子空间；
则代价函数可以表示为
Figure FDA0003581551080000048
步骤332，使用迭代法求解稀疏信号功率，设定第j+1次迭代得到唯一解，可得稀疏信号功率
Figure FDA0003581551080000051
式中，an为
Figure FDA0003581551080000052
的简写形式，
Figure FDA0003581551080000053
Figure FDA0003581551080000054
Figure FDA0003581551080000055
步骤333，根据步骤332求得的稀疏信号功率，得到N个网格对应信号的信号功率矢量为
Figure FDA0003581551080000056
CN202210353033.9A 2022-04-06 2022-04-06 非均匀噪声下声矢量传感器阵列波达方向估计方法 Pending CN114755628A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210353033.9A CN114755628A (zh) 2022-04-06 2022-04-06 非均匀噪声下声矢量传感器阵列波达方向估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210353033.9A CN114755628A (zh) 2022-04-06 2022-04-06 非均匀噪声下声矢量传感器阵列波达方向估计方法

Publications (1)

Publication Number Publication Date
CN114755628A true CN114755628A (zh) 2022-07-15

Family

ID=82328220

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210353033.9A Pending CN114755628A (zh) 2022-04-06 2022-04-06 非均匀噪声下声矢量传感器阵列波达方向估计方法

Country Status (1)

Country Link
CN (1) CN114755628A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114966640A (zh) * 2022-07-29 2022-08-30 宁波博海深衡科技有限公司 基于阵列背景噪声统计协方差估计的方位估计方法及***

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114966640A (zh) * 2022-07-29 2022-08-30 宁波博海深衡科技有限公司 基于阵列背景噪声统计协方差估计的方位估计方法及***
CN114966640B (zh) * 2022-07-29 2022-12-09 宁波博海深衡科技有限公司 基于阵列背景噪声统计协方差估计的方位估计方法及***

Similar Documents

Publication Publication Date Title
CN108375763B (zh) 一种应用于多声源环境的分频定位方法
CN109188344B (zh) 脉冲噪声环境下基于互循环相关music算法信源个数与来波方向角估计方法
CN110113085B (zh) 一种基于协方差矩阵重构的波束形成方法及***
CN107315162B (zh) 基于内插变换和波束形成的远场相干信号doa估计方法
CN110045321B (zh) 基于稀疏和低秩恢复的稳健doa估计方法
CN110045323B (zh) 一种基于矩阵填充的互质阵稳健自适应波束形成算法
CN109557504B (zh) 一种近场窄带信号源的定位方法
CN108761380B (zh) 一种用于提高精度的目标波达方向估计方法
CN112180339A (zh) 一种基于稀疏处理的雷达回波信号精确测向方法
CN114814830B (zh) 一种基于鲁棒主成分分析降噪的米波雷达低仰角测高方法
CN108828586B (zh) 一种基于波束域的双基地mimo雷达测角优化方法
CN114755628A (zh) 非均匀噪声下声矢量传感器阵列波达方向估计方法
CN110535519A (zh) 一种基于空间平滑的稳健自适应波束形成方法
Zuo et al. Improved capon estimator for high-resolution doa estimation and its statistical analysis
CN109541573A (zh) 一种弯曲水听器阵列的阵元位置校准方法
CN113625220A (zh) 一种多径信号波达方向和扩散角快速估计新方法
CN109541572B (zh) 一种基于线性环境噪声模型的子空间方位估计方法
CN108594165B (zh) 一种基于期望最大化算法的窄带信号波达方向估计方法
CN113093098B (zh) 基于lp范数补偿的轴向不一致矢量水听器阵列测向方法
CN115980721A (zh) 一种无误差协方差矩阵分离的阵列自校正方法
CN114167346B (zh) 基于协方差矩阵拟合阵元扩展的doa估计方法及***
CN114035157B (zh) 一种基于期望最大化算法的分频带时延估计方法及其***
CN111323750B (zh) 一种基于声矢量阵列网络的直接定位方法
CN113589223B (zh) 基于互耦情况下嵌套阵列的测向方法
CN114184999B (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