CN113504505A - 一种适用于低信噪比环境下的一维doa估计方法 - Google Patents

一种适用于低信噪比环境下的一维doa估计方法 Download PDF

Info

Publication number
CN113504505A
CN113504505A CN202110614512.7A CN202110614512A CN113504505A CN 113504505 A CN113504505 A CN 113504505A CN 202110614512 A CN202110614512 A CN 202110614512A CN 113504505 A CN113504505 A CN 113504505A
Authority
CN
China
Prior art keywords
matrix
signal
convex optimization
vector
iteration
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
CN202110614512.7A
Other languages
English (en)
Other versions
CN113504505B (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.)
South China University of Technology SCUT
Original Assignee
South China University of Technology SCUT
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 South China University of Technology SCUT filed Critical South China University of Technology SCUT
Priority to CN202110614512.7A priority Critical patent/CN113504505B/zh
Publication of CN113504505A publication Critical patent/CN113504505A/zh
Application granted granted Critical
Publication of CN113504505B publication Critical patent/CN113504505B/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/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/86—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 with means for eliminating undesired waves, e.g. disturbing noises
    • 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

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)

Abstract

本发明公开了一种适用于低信噪比环境下的一维DOA估计方法，该方法通过对均匀线阵的接收信号进行处理，在对接收信号的协方差矩阵做降维处理的同时减去噪声信号的特征值达到降噪的效果，并在稀疏向量解的凸优化过程中迭代地更新加权矩阵的系数来逐步修正角度位置的偏差，减小了角度的估计误差，提高压缩感知类算法在低信噪比环境下的DOA估计精度。本发明基于压缩感知估计理论对接收信号进行采样，使得采样频率不受奈奎斯特定理的限制，通过对接收数据进行降噪预处理和迭代更新权重系数，提高了水下DOA的估计精度，在低信噪比环境下具有较好的降噪效果。

Description

一种适用于低信噪比环境下的一维DOA估计方法
技术领域
本发明涉及水下目标定位技术领域，具体涉及一种适用于低信噪比环境下的一维DOA估计方法。
背景技术
波达方向估计(DOA估计)是一种根据阵列接收信号获取目标信号的空域位置信息的技术，在水下目标定位等领域有着重要的应用和发展。传统的DOA估计算法在信号采样过程中采样频率需要高于奈奎斯特频率，这对于采样后信号的存储空间和传输过程带来一定的负担和资源浪费。近年来提出的基于压缩感知的DOA估计算法可以突破传统采样机制对于信号频率的限制，大大降低了信号采样的最低频率，从而弥补传统算法的不足。现有的平滑L0压缩感知算法利用平滑的L0范数优化稀疏重构过程，提高了信号稀疏重构的精度的同时具有较低的运算复杂度；还有一种加权的压缩感知算法利用噪声子空间与信号特征向量间的正交性构建加权矩阵，代入到凸优化过程中以提高角度的估计精度。但以上算法在快拍数较少的情况下对于多信源波达方向角的分辨率较低。
此外在水声通信领域中，由于环境因素的复杂性，水声通信往往处于信噪比较低的环境，压缩感知类算法在低信噪比环境下存在估计性能下降的问题，尤其对于色噪声信号的抑制效果较差。目前亟待提出一种DOA估计方法解决上述算法的缺陷以及在低信噪比环境下的估计精度问题。
发明内容
本发明的目的是为了解决现有技术中的上述缺陷，提供一种适用于低信噪比环境下的一维DOA估计方法。该方法通过对均匀线阵的接收信号进行处理，在对接收信号的协方差矩阵做降维处理的同时减去噪声信号的特征值达到降噪的效果，并在稀疏向量解的凸优化过程中迭代地更新加权矩阵的系数来逐步修正角度位置的偏差，最终减小角度的估计误差，提高压缩感知类算法在低信噪比环境下的DOA估计精度。
本发明的目的可以通过采取如下技术方案达到：
一种适用于低信噪比环境下的一维DOA估计方法，所述估计方法包括以下步骤：
S1、建立一维的均匀线阵的阵列信号模型，经过L次采样得到接收信号矩阵Y，并根据阵列流型矩阵A(θ)构建凸优化过程所需的观测矩阵
Figure BDA0003096843170000021
S2、对接收信号矩阵Y进行奇异值分解，得到从大到小排序的M个奇异值，对后M-K个奇异值取平均值作为噪声信号的平均功率；
S3、降低接收信号和稀疏信号的矩阵维数，同时对接收信号矩阵Y做降噪处理，并初始化加权矩阵为W(0)＝IN，令i＝1，i表示迭代轮数，开始迭代更新权重系数，IN为N阶单位矩阵，N为空域划分的角度数；
S4、对凸优化方程求解得到最优解
Figure BDA0003096843170000022
并根据
Figure BDA0003096843170000023
构造虚拟的接收信号矩阵
Figure BDA0003096843170000024
对
Figure BDA0003096843170000025
做奇异值分解，得到第i轮的噪声子空间
Figure BDA0003096843170000026
S5、利用噪声子空间与观测矩阵中目标信号对应矢量的正交关系更新第i轮迭代的加权矩阵W(i)，令i＝i+1，重复步骤S4到S5，直到收敛条件终止迭代；
S6、输出最后一次的凸优化方程解
Figure BDA0003096843170000027
通过谱峰搜索获取前K个峰值相应的位置即为角度估计值。
进一步地，所述步骤S1的过程如下：假设均匀线阵上排布M个接收阵元，相邻阵元的间距为d，设有K个远场窄带声波信号入射到该均匀线阵上，阵元个数M和信源个数K满足关系K<M，每个信号对应的波长为λk，k＝1,2,…,K，真实信号的波达方向角分别θ＝{θ1,θ2,…,θK}，定义为信号入射方向与线阵所在直线的法线的夹角，接收阵元的采样信号y以列向量形式表示为：
Figure BDA0003096843170000031
其中，
Figure BDA0003096843170000032
表示阵列流型矩阵，
Figure BDA0003096843170000033
表示K个信源对应的目标信号矢量，
Figure BDA0003096843170000034
表示M个阵元的噪声信号矢量，
Figure BDA0003096843170000035
表示复数域，经过L次采样后得到的接收信号矩阵Y表达式如公式(2)所示：
Y＝A(θ)S+G 公式(2)
其中，S表示包含K个信源信息的K×L维目标信号矩阵，G表示M×L维的噪声矩阵，A(θ)表示M×K维的阵列流型矩阵，其表达式为：
Figure BDA0003096843170000036
通过压缩感知类算法求解凸优化方程获得的N×K维稀疏解，该稀疏解中包含K个真实信源的角度信息，其中，凸优化方程的构建需要根据阵列流型矩阵中的元素构建M×N维的观测矩阵
Figure BDA0003096843170000037
该观测矩阵
Figure BDA0003096843170000038
的表示形式为：
Figure BDA0003096843170000039
其中，N表示接收阵元上空域划分的角度个数，满足M<<N，λ表示信号的中心频率对应的波长值；
根据压缩感知理论，通过对接收信号矩阵Y进行凸优化处理得到稀疏表示信号的最优解，该最优解中包含目标信号的入射角度，凸优化方程如公式(5)所示：
Figure BDA0003096843170000041
在凸优化方程中的观测矩阵
Figure BDA0003096843170000042
基于阵列流型矩阵A(θ)中原子的参数构建；
针对接收信号模型中第一个阵元的接收信号，将目标信号矢量
Figure BDA0003096843170000043
扩展成包含K个真实信号的N×1维虚拟原始信号矢量
Figure BDA0003096843170000044
N表示空域划分的角度数，角度集合Θ所对应的虚拟原始信号矢量
Figure BDA0003096843170000045
是一个K-稀疏信号，表明
Figure BDA0003096843170000046
是可压缩的，在原始信号的稀疏解中包含来波方向角信息的标量位于虚拟原始信号矢量
Figure BDA0003096843170000047
的非零元素位置。
进一步地，所述步骤S2过程如下：
首先对接收信号矩阵Y进行奇异值分解构建降噪的奇异值矩阵Σ′：
Y＝UΣVH 公式(6)
Σ′＝Σ-ΣN 公式(7)
Figure BDA0003096843170000048
其中，
Figure BDA0003096843170000049
表示主对角线元素为奇异值构成的矩阵，
Figure BDA00030968431700000410
表示实数域，
Figure BDA00030968431700000411
表示Σ的主对角线元素中第K+1个到第M个奇异值之和的平均值作为噪声信号的平均功率。
进一步地，所述步骤S3过程如下：
对接收信号矩阵Y降维处理的同时，减去噪声平均功率
Figure BDA0003096843170000051
得到处理后的接收信号矩阵YSVW：
YSVW＝UΣ′DK 公式(9)
Figure BDA0003096843170000052
其中，
Figure BDA0003096843170000053
凸优化方程中待估计信号的L1范数替代公式(5)中的L0范数作为衡量最优结果的参数，将约束条件转换为误差项并加入加权矩阵，凸优化方程变为：
Figure BDA0003096843170000054
其中，
Figure BDA0003096843170000055
表示使凸优化方程解为最小值的XSV取值，
Figure BDA0003096843170000056
和
Figure BDA0003096843170000057
分别表示稀疏表示信号的初始值和凸优化方程解，
Figure BDA0003096843170000058
表示矩阵的Frobenius范数的平方，||·||2,1表示矩阵中各列向量L2范数组成向量的L1范数值，
Figure BDA0003096843170000059
表示加权矩阵，公式(11)对应的二阶锥规划问题表示为公式(12)：
min p+ηq 公式(12)
min表示求解p+ηq的最小值，其约束条件为存在K个N×1维的矢量[Z1,Z2,…,Zk,…,ZK]和N×1维的矢量[γ1,γ2,…,γj,…,γN]T满足关系
Figure BDA00030968431700000510
且||γ1,γ2,…,γj,…,γN||1≤q，Zk和γn的表达式为：
Figure BDA0003096843170000061
||WXSV(j,:)||2≤γj,j＝1,2,…,N 公式(14)
其中，p、q为约束条件中两个范数和的临时变量，η表示正则化参数因子，用来平衡误差和稀疏性之间的关系，
Figure BDA0003096843170000062
表示经过降噪和降维处理后的接收信号矩阵，YSVW(:,k)为YSVW的第k列列向量，XSV(:,k)表示XSV的第k列列向量，
Figure BDA0003096843170000063
表示
Figure BDA0003096843170000064
的第j行行向量，
Figure BDA0003096843170000065
表示加权矩阵，
Figure BDA0003096843170000066
表示的Zk共轭向量，||·||1、||·||2和
Figure BDA0003096843170000067
分别表示向量的L1范数和L2范数以及Frobenius范数的平方，同时开始迭代前的加权矩阵W(0)初始化为N×N维的单位矩阵IN，重新令i＝1，开始进入第一轮迭代。
进一步地，所述步骤S4过程如下：
将第i-1轮迭代的加权矩阵W(i-1)代入公式(11)中进行凸优化方程的求解，得到第i轮凸优化方程解
Figure BDA0003096843170000068
根据其估计值构建出虚拟的接收信号矩阵
Figure BDA0003096843170000069
并对
Figure BDA00030968431700000610
进行奇异值分解：
Figure BDA00030968431700000611
Figure BDA00030968431700000612
其中，U(i)、Σ(i)、V(i)分别为第i轮迭代中
Figure BDA00030968431700000613
的左奇异矩阵、奇异值矩阵和右奇异矩阵，U(i)矩阵中的前K列子矩阵
Figure BDA00030968431700000614
和后M-K列子矩阵
Figure BDA00030968431700000615
分别对应信号子空间和噪声子空间。
进一步地，所述步骤S5过程如下：
利用公式(16)得到的噪声子空间
Figure BDA00030968431700000616
和包含真实信号角度信息的观测矩阵列向量，结合MUSIC算法的原理构造第i轮迭代的加权矩阵W(i)：
Figure BDA0003096843170000071
Figure BDA0003096843170000072
其中，
Figure BDA0003096843170000073
表示第i轮加权矩阵中的对角线上元素，
Figure BDA0003096843170000074
表示第i轮迭代中噪声子空间的共轭矩阵，a(θj),j＝1,2,…,N表示观测矩阵
Figure BDA0003096843170000075
中第j列的列向量，根据噪声子空间与真实信号矢量的正交性，在真实信号角度对应列向量的加权系数会增大其谱峰，而其它位置的列向量对应谱峰值则会减小，令i＝i+1，进入下一轮迭代，重复步骤S4和步骤S5，并将本轮迭代更新的加权矩阵W(i)代入下一轮迭代的凸优化方程中，直到满足迭代终止条件为止。
进一步地，所述步骤S5中迭代更新的终止条件包括以下两种：
(1)权重更新达到足够的迭代次数后终止；
(2)第i轮迭代中的凸优化方程解
Figure BDA0003096843170000076
与第i-1轮迭代中的凸优化方程解
Figure BDA0003096843170000077
之差的L2范数小于给定的误差参数时终止；
Figure BDA0003096843170000078
其中，ε为给定的误差参数，根据实际需要取值。
进一步地，所述步骤S6过程如下：
获取第i轮的凸优化方程解
Figure BDA0003096843170000079
根据
Figure BDA00030968431700000710
构建稀疏表示信号的空间谱密度函数Q：
Figure BDA00030968431700000711
其中，
Figure BDA00030968431700000712
表示
Figure BDA00030968431700000713
的第一列列向量，
Figure BDA00030968431700000714
表示该向量的L1范数，通过对Q进行谱峰搜索获取前K个峰值的位置，即对应虚拟的原始信号矢量
Figure BDA0003096843170000081
中真实信号入射角度的估计值
Figure BDA0003096843170000082
本发明相对于现有技术具有如下的优点及效果：
1、本发明公开的DOA估计方法基于压缩感知理论对接收信号进行压缩采样，相比传统的子空间类算法，该方法的采样频率不受奈奎斯特定理的限制，降低了采样后的数据存储和传输所消耗的资源，且通过稀疏重构即可恢复出原始信号。
2、本发明公开的DOA估计方法对接收信号矩阵进行了降维和降噪处理，在降低方法的运算复杂度的同时也提高了对噪声信号的抑制效果，在低信噪比的白噪声和色噪声环境下的DOA估计都具有更好的降噪性能。
3、本发明公开的DOA估计方法相比现有的加权凸优化算法具有精度更高的加权系数，通过迭代更新权重缩减了凸优化方程最优解的估计误差，进一步提高了算法的估计精度。
4、本发明公开的DOA估计方法对于多信源情况下的相邻入射角估计时的分辨力更好，相比传统的压缩感知类算法在抗模糊性上具有更高的性能。
附图说明
图1是本发明实施例提出的一种适用于低信噪比环境下的一维DOA估计方法中均匀线阵模型示意图；
图2是本发明实施例提出的一种适用于低信噪比环境下的一维DOA估计方法中接收信号模型示意图；
图3是本发明实施例提出的一种适用于低信噪比环境下的一维DOA估计方法的流程图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
实施例一
本实施例公开了一种适用于低信噪比环境下的一维DOA估计方法，所述估计方法包括以下步骤：
S1、建立一维的均匀线阵的阵列信号模型，经过L次采样得到接收信号矩阵Y，并根据阵列流型矩阵构建凸优化过程所需的观测矩阵
Figure BDA0003096843170000091
本实施例步骤S1的过程如下：
采用均匀线阵为接收阵列和基于远场窄带信号模型进行DOA估计，均匀线阵模型如附图1所示，假设均匀线阵上排布M个接收阵元，相邻阵元的间距为d，设有K个远场窄带声波信号入射到该均匀线阵上，阵元个数M和信源个数K满足关系K<M，每个信号对应的波长为λk，k＝1,2,…,K，真实信号的波达方向角分别θ＝{θ1,θ2,…,θK}，定义为信号入射方向与线阵所在直线的法线的夹角，接收阵元的采样信号y以列向量形式表示为：
Figure BDA0003096843170000092
其中，
Figure BDA0003096843170000093
表示阵列流型矩阵，
Figure BDA0003096843170000094
表示K个信源对应的目标信号矢量，
Figure BDA0003096843170000095
表示M个阵元的噪声信号矢量，
Figure BDA0003096843170000096
表示复数域，经过L次采样后得到的接收信号矩阵Y表达式如公式(2)所示：
Y＝A(θ)S+G 公式(2)
其中，S表示包含K个信源信息的K×L维目标信号矩阵，G表示M×L维的噪声矩阵，A(θ)表示M×K维的阵列流型矩阵，其表达式为：
Figure BDA0003096843170000101
压缩感知类算法通过对接收信号矩阵Y求解凸优化方程获得的N×K维最优解
Figure BDA0003096843170000102
中包含K个真实信源的角度信息，其中，凸优化方程的构建需要根据阵列流型矩阵A(θ)中的元素构建M×N维的观测矩阵
Figure BDA0003096843170000103
该观测矩阵
Figure BDA0003096843170000104
的表示形式为：
Figure BDA0003096843170000105
其中，N表示接收阵元上空域划分的角度个数，满足M<<N，且N越大，角度估计的位置越精确，λ表示信号的中心频率对应的波长值；
根据压缩感知理论，通过对接收信号矩阵Y进行凸优化处理得到稀疏表示信号的最优解，该最优解中包含目标信号的入射角度，凸优化方程如公式(5)所示：
Figure BDA0003096843170000106
在凸优化方程中的观测矩阵
Figure BDA0003096843170000107
基于阵列流型矩阵A(θ)中原子的参数构建。
接收信号模型如附图2所示，针对该接收信号模型中各阵元的接收信号，将目标信号矢量
Figure BDA0003096843170000108
扩展成包含K个真实信号的N×1维虚拟原始信号矢量
Figure BDA0003096843170000109
角度集合Θ所对应的虚拟原始信号矢量
Figure BDA00030968431700001010
是一个K-稀疏信号，表明
Figure BDA0003096843170000111
是可压缩的，在原始信号的稀疏解中包含来波方向角信息的标量位于虚拟原始信号矢量
Figure BDA0003096843170000112
的非零元素位置。
S2、对接收信号矩阵Y进行奇异值分解，得到从大到小排序的M个奇异值，对后M-K个奇异值取平均值作为噪声信号的平均功率；
本实施例的步骤S2过程如下：
在凸优化之前需要对接收信号矩阵Y做降维和降噪处理，目的是为了降低运算复杂度的同时抑制噪声的干扰。首先对接收信号矩阵Y进行奇异值分解构建降噪的奇异值矩阵Σ′：
Y＝UΣVH 公式(6)
Σ′＝Σ-ΣN 公式(7)
Figure BDA0003096843170000113
其中，
Figure BDA0003096843170000114
表示主对角线元素为奇异值构成的矩阵，
Figure BDA0003096843170000115
表示实数域，
Figure BDA0003096843170000116
表示Σ的主对角线元素中第K+1个到第M个奇异值之和的平均值作为噪声信号的平均功率。
S3、降低接收信号和稀疏信号的矩阵维数，同时对接收信号矩阵Y做降噪处理，并在迭代开始前初始化加权矩阵为W(0)＝IN，令i＝1，开始迭代更新权重系数，IN为N阶单位矩阵，N表示空域划分的角度数；
本实施例的步骤S3过程如下：
对接收信号矩阵Y降维处理的同时，减去噪声平均功率
Figure BDA0003096843170000117
得到处理后的接收信号矩阵YSVW：
YSVW＝UΣ′DK 公式(9)
Figure BDA0003096843170000121
其中，
Figure BDA0003096843170000122
通过右乘DK将列数为N的原接收信号矩阵压缩到只包含前K列信息，且保证包含信源信息的奇异值保留到降维后的矩阵中。相比现有的压缩感知估计方法，步骤S2在对接收信号矩阵Y作降维处理的同时，增加了降噪处理的步骤，目的是为了减少凸优化过程中由于噪声影响导致的估计偏差；另外，凸优化方程中待估计信号的L1范数替代了公式(5)中的L0范数作为衡量最优结果的参数，将约束条件转换为误差项并加入加权矩阵，凸优化方程变为：
Figure BDA0003096843170000123
其中，
Figure BDA0003096843170000124
表示使凸优化方程解为最小值的XSV取值，
Figure BDA0003096843170000125
和
Figure BDA0003096843170000126
分别表示稀疏表示信号的初始值和凸优化方程解，
Figure BDA0003096843170000127
表示矩阵的Frobenius范数的平方，||·||2,1表示矩阵中各列向量L2范数组成向量的L1范数值，
Figure BDA0003096843170000128
表示加权矩阵。公式(11)对应的二阶锥规划问题表示为公式(12)：
min p+ηq 公式(12)
min表示求解p+ηq的最小值，其约束条件为存在K个N×1维的矢量[Z1,Z2,…,Zk,…,ZK]和N×1维的矢量[γ1,γ2,…,γj,…,γN]T满足关系
Figure BDA0003096843170000129
且||γ1,γ2,…,γj,…,γN||1≤q。Zk和γn的表达式为：
Figure BDA00030968431700001210
||WXSV(j,:)||2≤γj,j＝1,2,…,N 公式(14)
其中，p、q为约束条件中两个范数和的临时变量，η表示正则化参数因子，用来平衡误差和稀疏性之间的关系，
Figure BDA0003096843170000131
表示经过降噪和降维处理后的接收信号矩阵，YSVW(:,k)为YSVW的第k列列向量，XSV(:,k)表示XSV的第k列列向量，
Figure BDA0003096843170000132
表示
Figure BDA0003096843170000133
的第j行行向量，
Figure BDA0003096843170000134
表示加权矩阵，
Figure BDA0003096843170000135
表示的Zk共轭向量，||·||1、||·||2和
Figure BDA0003096843170000136
分别表示向量的L1范数和L2范数以及Frobenius范数的平方。同时开始迭代前的加权矩阵W(0)初始化为N×N维的单位矩阵IN，令i＝1，开始进入第一轮迭代。
S4、对凸优化方程求解得到最优解
Figure BDA0003096843170000137
并根据
Figure BDA0003096843170000138
构造虚拟的接收信号矩阵
Figure BDA0003096843170000139
对
Figure BDA00030968431700001310
做奇异值分解，得到第i轮的噪声子空间
Figure BDA00030968431700001311
本实施例的步骤S4过程如下：
将加权矩阵W(i-1)代入公式(11)中进行凸优化方程的求解，得到第i轮凸优化方程解
Figure BDA00030968431700001312
根据其估计值构建出虚拟的接收信号矩阵
Figure BDA00030968431700001313
并对
Figure BDA00030968431700001314
进行奇异值分解：
Figure BDA00030968431700001315
Figure BDA00030968431700001316
其中，U(i)、Σ(i)、V(i)分别为第i轮迭代中
Figure BDA00030968431700001317
的左奇异矩阵、奇异值矩阵和右奇异矩阵，U(i)矩阵中的前K列子矩阵
Figure BDA00030968431700001318
和后M-K列子矩阵
Figure BDA00030968431700001319
分别对应信号子空间和噪声子空间；由于加权系数构成需要由接收信号矩阵，因此通过迭代的方式利用上一轮的凸优化方程解构建虚拟的接收信号矩阵可以作为本轮的加权系数构成，与加权的L1-SVD估计方法相比，能够实现迭代地更新加权系数，进一步提高估计的精度。
S5、利用噪声子空间与观测矩阵中目标信号对应矢量的正交关系更新第i轮加权矩阵W(i)，令i＝i+1，重复步骤S4到S5，直到收敛条件终止迭代；
利用式(16)得到的噪声子空间
Figure BDA0003096843170000141
和包含真实信号角度信息的观测矩阵列向量，结合MUSIC算法的原理构造第i轮迭代的加权矩阵W(i)：
Figure BDA0003096843170000142
Figure BDA0003096843170000143
其中，
Figure BDA0003096843170000144
表示第i轮加权矩阵中的对角线上元素，
Figure BDA0003096843170000145
表示第i轮迭代中噪声子空间的共轭矩阵，a(θj),j＝1,2,…,N表示观测矩阵
Figure BDA0003096843170000146
中第j列的向量，根据噪声子空间与真实信号矢量的正交性，在真实信号角度对应列向量的加权系数会增大其谱峰，而其它位置的列向量对应谱峰值则会减小，令i＝i+1，进入下一轮迭代，重复步骤S4和步骤S5，并将本轮迭代更新的加权矩阵W(i)代入下一轮迭代的凸优化方程中，直到满足迭代终止条件为止。
迭代更新的终止条件包括以下两种：
(1)权重更新达到足够的迭代次数后终止；
(2)第i轮迭代中的凸优化方程解
Figure BDA0003096843170000147
与第i-1轮迭代中的凸优化方程解
Figure BDA0003096843170000148
之差的L2范数小于给定的误差参数时终止；
Figure BDA0003096843170000149
其中，ε为给定的误差参数，根据实际需要取值。
S6、输出最后一次的凸优化方程解
Figure BDA00030968431700001410
通过谱峰搜索获取前K个峰值相应的位置即为角度估计值。
本实施例的步骤S6过程如下：
获取最后一次的凸优化方程解
Figure BDA0003096843170000151
根据
Figure BDA0003096843170000152
构建稀疏表示信号的空间谱密度函数Q：
Figure BDA0003096843170000153
其中，
Figure BDA0003096843170000154
表示
Figure BDA0003096843170000155
的第一列列向量，
Figure BDA0003096843170000156
表示该向量的L1范数。通过对Q进行谱峰搜索获取前K个峰值的位置，即对应原始信号矢量
Figure BDA0003096843170000157
中真实信号入射角度的估计值
Figure BDA0003096843170000158
与加权的L1-SVD算法相比，加权系数对角度位置加成更加准确，在步骤S1的条件下经过仿真对比，加权L1-SVD估计方法的估计偏差为0.18°，本算法的估计偏差为0.07°。
实施例二
本实施例提出一种基于迭代更新权重和一维均匀线阵的水下一维DOA估计方法，具体步骤如下：
T1、建立一维的均匀线阵的阵列信号模型。采用一维均匀线阵为接收阵列，阵元数M＝8，设有K个非相关的远场窄带声波信号入射到该均匀线阵上，信源个数K＝3，入射角度分别为[-30°，0°，30°]，信号的传播速度为c＝1500m/s，阵元间距d＝0.05m，为信源信号的中心频率15kHz对应波长的一半，阵列单次采样的接收信号形式为：
y＝[y1,y2,…,y8]T 公式(21)
经过快拍数L＝100次采样得到的接收信号矩阵为：
Figure BDA0003096843170000159
根据阵列流型矩阵A(θ)表达式构建凸优化方程所需的观测矩阵
Figure BDA00030968431700001510
其表达式如公式(6)所示。其中
Figure BDA0003096843170000161
N表示在接收阵元上空域划分的角度数量，且满足关系M<<N。这里根据附图2规定角度搜索的空域范围为[-90°～90°]，以均匀间隔为0.1°划分，共划分网格的角度数量N＝1801。
T2、对接收信号矩阵Y做奇异值分解得到奇异值矩阵Σ，得到共8个奇异值，从大到小对所有奇异值排列，取后5个奇异值的平均值作为噪声信号的平均功率，并根据公式(8)、公式(9)、公式(10)构建降噪后的奇异值矩阵Σ′。
T3、对原接收信号矩阵Y做降噪处理，减去噪声平均功率
Figure BDA0003096843170000162
并根据公式(11)对接收信号矩阵Y进行降维，得到处理后的接收信号矩阵YSVW。
T4、初始化加权矩阵W(0)＝IN。令i＝1，开始对加权矩阵进行更新。
T5、根据公式(13)求解凸优化方程得到第i轮凸优化方程解
Figure BDA0003096843170000163
并依此构建第i轮迭代中虚拟的接收信号矩阵
Figure BDA0003096843170000164
对
Figure BDA0003096843170000165
做奇异值分解得到
Figure BDA0003096843170000166
的噪声子空间
Figure BDA0003096843170000167
T6、根据公式(17)、(18)更新第i轮迭代的加权矩阵W(i)。令i＝i+1，进入下一轮迭代，重复S4到S6的步骤，并将第i-1轮迭代更新的加权矩阵W(i-1)代入下一轮迭代的凸优化方程中，直到迭代轮数大于等于6。
T7、经过6次权重更新后结束迭代过程，根据第6轮的凸优化方程解
Figure BDA0003096843170000168
构建稀疏表示信号的空间谱密度函数Q，并进行谱峰搜索，得到前3个谱峰值对应的角度为[-29.84°，0.09°，30.02°]，对目标估计达到了预期精度，说明估计结果正确，本发明方法可行。
上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (8)

1.一种适用于低信噪比环境下的一维DOA估计方法，其特征在于，所述估计方法包括以下步骤：
S1、建立一维的均匀线阵的阵列信号模型，经过L次采样得到接收信号矩阵Y，并根据阵列流型矩阵A(θ)构建凸优化过程所需的观测矩阵
Figure FDA0003096843160000017
S2、对接收信号矩阵Y进行奇异值分解，得到从大到小排序的M个奇异值，对后M-K个奇异值取平均值作为噪声信号的平均功率；
S3、降低接收信号和稀疏信号的矩阵维数，同时对接收信号矩阵Y做降噪处理，并初始化加权矩阵为W(0)＝IN，令i＝1，i表示迭代轮数，开始迭代更新权重系数，IN为N阶单位矩阵，N为空域划分的角度数；
S4、对凸优化方程求解得到最优解
Figure FDA0003096843160000011
并根据
Figure FDA0003096843160000012
构造虚拟的接收信号矩阵
Figure FDA0003096843160000013
对
Figure FDA0003096843160000014
做奇异值分解，得到第i轮的噪声子空间
Figure FDA0003096843160000015
S5、利用噪声子空间与观测矩阵中目标信号对应矢量的正交关系更新第i轮迭代的加权矩阵W(i)，令i＝i+1，重复步骤S4到S5，直到收敛条件终止迭代；
S6、输出最后一次的凸优化方程解
Figure FDA0003096843160000016
通过谱峰搜索获取前K个峰值相应的位置即为角度估计值。
2.根据权利要求1所述的一种基于迭代更新权重的一维DOA估计方法，其特征在于，所述步骤S1的过程如下：假设均匀线阵上排布M个接收阵元，相邻阵元的间距为d，设有K个远场窄带声波信号入射到该均匀线阵上，阵元个数M和信源个数K满足关系K<M，每个信号对应的波长为λk，k＝1,2,…,K，真实信号的波达方向角分别θ＝{θ1,θ2,…,θK}，定义为信号入射方向与线阵所在直线的法线的夹角，接收阵元的采样信号y以列向量形式表示为：
Figure FDA0003096843160000021
其中，
Figure FDA0003096843160000022
表示阵列流型矩阵，
Figure FDA0003096843160000023
表示K个信源对应的目标信号矢量，
Figure FDA0003096843160000024
表示M个阵元的噪声信号矢量，
Figure FDA0003096843160000025
表示复数域，经过L次采样后得到的接收信号矩阵Y表达式如公式(2)所示：
Y＝A(θ)S+G 公式(2)
其中，S表示包含K个信源信息的K×L维目标信号矩阵，G表示M×L维的噪声矩阵，A(θ)表示M×K维的阵列流型矩阵，其表达式为：
Figure FDA0003096843160000026
通过压缩感知类算法求解凸优化方程获得的N×K维稀疏解，该稀疏解中包含K个真实信源的角度信息，其中，凸优化方程的构建需要根据阵列流型矩阵中的元素构建M×N维的观测矩阵
Figure FDA0003096843160000027
该观测矩阵
Figure FDA0003096843160000028
的表示形式为：
Figure FDA0003096843160000029
其中，N表示接收阵元上空域划分的角度个数，满足M<<N，λ表示信号的中心频率对应的波长值；
根据压缩感知理论，通过对接收信号矩阵Y进行凸优化处理得到稀疏表示信号的最优解，该最优解中包含目标信号的入射角度，凸优化方程如公式(5)所示：
Figure FDA0003096843160000031
在凸优化方程中的观测矩阵
Figure FDA0003096843160000032
基于阵列流型矩阵A(θ)中原子的参数构建；
针对接收信号模型中第一个阵元的接收信号，将目标信号矢量
Figure FDA0003096843160000033
扩展成包含K个真实信号的N×1维虚拟原始信号矢量
Figure FDA0003096843160000034
N表示空域划分的角度数，角度集合Θ所对应的虚拟原始信号矢量
Figure FDA0003096843160000035
是一个K-稀疏信号，表明
Figure FDA0003096843160000036
是可压缩的，在原始信号的稀疏解中包含来波方向角信息的标量位于虚拟原始信号矢量
Figure FDA0003096843160000037
的非零元素位置。
3.根据权利要求2所述的一种基于迭代更新权重的一维DOA估计方法，其特征在于，所述步骤S2过程如下：
首先对接收信号矩阵Y进行奇异值分解构建降噪的奇异值矩阵Σ′：
Y＝UΣVH 公式(6)
Σ′＝Σ-ΣN 公式(7)
Figure FDA0003096843160000038
其中，
Figure FDA0003096843160000039
表示主对角线元素为奇异值构成的矩阵，
Figure FDA00030968431600000310
表示实数域，
Figure FDA00030968431600000311
表示Σ的主对角线元素中第K+1个到第M个奇异值之和的平均值作为噪声信号的平均功率。
4.根据权利要求3所述的一种基于迭代更新权重的一维DOA估计方法，其特征在于，所述步骤S3过程如下：
对接收信号矩阵Y降维处理的同时，减去噪声平均功率
Figure FDA0003096843160000041
得到处理后的接收信号矩阵YSVW：
YSVW＝UΣ′DK 公式(9)
Figure FDA0003096843160000042
其中，
Figure FDA0003096843160000043
凸优化方程中待估计信号的L1范数替代公式(5)中的L0范数作为衡量最优结果的参数，将约束条件转换为误差项并加入加权矩阵，凸优化方程变为：
Figure FDA0003096843160000044
其中，
Figure FDA0003096843160000045
表示使凸优化方程解为最小值的XSV取值，
Figure FDA0003096843160000046
和
Figure FDA0003096843160000047
分别表示稀疏表示信号的初始值和凸优化方程解，
Figure FDA0003096843160000048
表示矩阵的Frobenius范数的平方，||·||2,1表示矩阵中各列向量L2范数组成向量的L1范数值，
Figure FDA0003096843160000049
表示加权矩阵，公式(11)对应的二阶锥规划问题表示为公式(12)：
min p+ηq 公式(12)
min表示求解p+ηq的最小值，其约束条件为存在K个N×1维的矢量[Z1,Z2,…,Zk,…,ZK]和N×1维的矢量[γ1,γ2,…,γj,…,γN]T满足关系
Figure FDA0003096843160000051
且||γ1,γ2,…,γj,…,γN||1≤q，Zk和γj的表达式为：
Figure FDA0003096843160000052
||WXSV(j,:)||2≤γj,j＝1,2,…,N 公式(14)
其中，p、q为约束条件中两个范数和的临时变量，η表示正则化参数因子，用来平衡误差和稀疏性之间的关系，
Figure FDA0003096843160000053
表示经过降噪和降维处理后的接收信号矩阵，YSVW(:,k)为YSVW的第k列列向量，XSV(:,k)表示XSV的第k列列向量，
Figure FDA0003096843160000054
表示
Figure FDA0003096843160000055
的第j行行向量，
Figure FDA0003096843160000056
表示加权矩阵，
Figure FDA0003096843160000057
表示的Zk共轭向量，||·||1、||·||2和
Figure FDA0003096843160000058
分别表示向量的L1范数和L2范数以及Frobenius范数的平方，同时开始迭代前的加权矩阵W(0)初始化为N×N维的单位矩阵IN，重新令i＝1，开始进入第一轮迭代。
5.根据权利要求4所述的一种基于迭代更新权重的一维DOA估计方法，其特征在于，所述步骤S4过程如下：
将第i-1轮迭代的加权矩阵W(i-1)代入公式(11)中进行凸优化方程的求解，得到第i轮凸优化方程解
Figure FDA0003096843160000059
根据其估计值构建出虚拟的接收信号矩阵
Figure FDA00030968431600000510
并对
Figure FDA00030968431600000511
进行奇异值分解：
Figure FDA00030968431600000512
Figure FDA00030968431600000513
其中，U(i)、Σ(i)、V(i)分别为第i轮迭代中
Figure FDA00030968431600000514
的左奇异矩阵、奇异值矩阵和右奇异矩阵，U(i)矩阵中的前K列子矩阵
Figure FDA00030968431600000515
和后M-K列子矩阵
Figure FDA00030968431600000516
分别对应信号子空间和噪声子空间。
6.根据权利要求5所述的一种基于迭代更新权重的一维DOA估计方法，其特征在于，所述步骤S5过程如下：
利用式(16)得到的噪声子空间
Figure FDA00030968431600000517
和包含真实信号角度信息的观测矩阵列向量，结合MUSIC算法的原理构造第i轮迭代的加权矩阵W(i)：
Figure FDA0003096843160000061
Figure FDA0003096843160000062
其中，
Figure FDA0003096843160000063
表示第i轮加权矩阵中的对角线上元素，
Figure FDA0003096843160000064
表示第i轮迭代中噪声子空间的共轭矩阵，a(θj),j＝1,2,…,N表示观测矩阵
Figure FDA0003096843160000065
中第j列的列向量，根据噪声子空间与真实信号矢量的正交性，在真实信号角度对应列向量的加权系数会增大其谱峰，而其它位置的列向量对应谱峰值则会减小，令i＝i+1，进入下一轮迭代，重复步骤S4和步骤S5，并将本轮迭代更新的加权矩阵W(i)代入下一轮迭代的凸优化方程中，直到满足迭代终止条件为止。
7.根据权利要求6所述的一种基于迭代更新权重的一维DOA估计方法，其特征在于，所述步骤S5中迭代更新的终止条件包括以下两种：
(1)权重更新达到足够的迭代次数后终止；
(2)第i轮迭代中的凸优化方程解
Figure FDA0003096843160000066
与第i-1轮迭代中的凸优化方程解
Figure FDA0003096843160000067
之差的L2范数小于给定的误差参数时终止；
Figure FDA0003096843160000068
其中，ε为给定的误差参数，根据实际需要取值。
8.根据权利要求6所述的一种基于迭代更新权重的一维DOA估计方法，其特征在于，所述步骤S6过程如下：
获取第i轮的凸优化方程解
Figure FDA0003096843160000069
根据
Figure FDA00030968431600000610
构建稀疏表示信号的空间谱密度函数Q：
Figure FDA0003096843160000071
其中，
Figure FDA0003096843160000072
表示
Figure FDA0003096843160000073
的第一列列向量，
Figure FDA0003096843160000074
表示该向量的L1范数，通过对Q进行谱峰搜索获取前K个峰值的位置，即对应虚拟的原始信号矢量
Figure FDA0003096843160000075
中真实信号入射角度的估计值
Figure FDA0003096843160000076
CN202110614512.7A 2021-06-02 2021-06-02 一种适用于低信噪比环境下的一维doa估计方法 Active CN113504505B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110614512.7A CN113504505B (zh) 2021-06-02 2021-06-02 一种适用于低信噪比环境下的一维doa估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110614512.7A CN113504505B (zh) 2021-06-02 2021-06-02 一种适用于低信噪比环境下的一维doa估计方法

Publications (2)

Publication Number Publication Date
CN113504505A true CN113504505A (zh) 2021-10-15
CN113504505B CN113504505B (zh) 2023-11-03

Family

ID=78009296

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110614512.7A Active CN113504505B (zh) 2021-06-02 2021-06-02 一种适用于低信噪比环境下的一维doa估计方法

Country Status (1)

Country Link
CN (1) CN113504505B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116009058A (zh) * 2022-11-30 2023-04-25 杭州交大仪器设备有限公司 一种基于多探头传感器数据的地下管道定位方法
CN116879835A (zh) * 2023-07-25 2023-10-13 安徽大学 一种投影最小最大凹函数波达方向估计方法和装置

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106772225A (zh) * 2017-01-20 2017-05-31 大连大学 基于压缩感知的波束域doa估计
CN107544052A (zh) * 2017-08-07 2018-01-05 大连大学 一种基于矩阵补全的二阶统计量重构doa估计方法
WO2018094565A1 (zh) * 2016-11-22 2018-05-31 深圳大学 脉冲噪声下的波束成形方法及装置
CN110109053A (zh) * 2019-04-02 2019-08-09 华南理工大学 一种未知声速环境下快速doa估计方法
CN110308417A (zh) * 2019-05-30 2019-10-08 电子科技大学 基于矩阵填充的嵌套阵阵元失效下的波达方向估计方法及装置
CN111190136A (zh) * 2020-01-08 2020-05-22 华南理工大学 一种基于特定频率组合信号的一维doa估计方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2018094565A1 (zh) * 2016-11-22 2018-05-31 深圳大学 脉冲噪声下的波束成形方法及装置
CN106772225A (zh) * 2017-01-20 2017-05-31 大连大学 基于压缩感知的波束域doa估计
CN107544052A (zh) * 2017-08-07 2018-01-05 大连大学 一种基于矩阵补全的二阶统计量重构doa估计方法
CN110109053A (zh) * 2019-04-02 2019-08-09 华南理工大学 一种未知声速环境下快速doa估计方法
CN110308417A (zh) * 2019-05-30 2019-10-08 电子科技大学 基于矩阵填充的嵌套阵阵元失效下的波达方向估计方法及装置
CN111190136A (zh) * 2020-01-08 2020-05-22 华南理工大学 一种基于特定频率组合信号的一维doa估计方法

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116009058A (zh) * 2022-11-30 2023-04-25 杭州交大仪器设备有限公司 一种基于多探头传感器数据的地下管道定位方法
CN116009058B (zh) * 2022-11-30 2024-04-19 杭州交大仪器设备有限公司 一种基于多探头传感器数据的地下管道定位方法
CN116879835A (zh) * 2023-07-25 2023-10-13 安徽大学 一种投影最小最大凹函数波达方向估计方法和装置

Also Published As

Publication number Publication date
CN113504505B (zh) 2023-11-03

Similar Documents

Publication Publication Date Title
CN107290709B (zh) 基于范德蒙分解的互质阵列波达方向估计方法
CN110109051B (zh) 基于频控阵的互耦阵列doa估计方法
CN111337893B (zh) 一种基于实值稀疏贝叶斯学习的离格doa估计方法
CN107092004B (zh) 基于信号子空间旋转不变性的互质阵列波达方向估计方法
CN106093921B (zh) 基于稀疏分解理论的声矢量阵宽带测向方法
CN109298383B (zh) 一种基于变分贝叶斯推断的互质阵波达方向角估计方法
CN108872926B (zh) 一种基于凸优化的幅相误差校正及doa估计方法
CN113504505B (zh) 一种适用于低信噪比环境下的一维doa估计方法
CN109375154B (zh) 一种冲击噪声环境下基于均匀圆阵的相干信号参数估计方法
CN107544051A (zh) 嵌套阵列基于k‑r子空间的波达方向估计方法
CN104777491A (zh) 一种盲波束宽带干扰抑制方法及装置
CN112731273A (zh) 一种基于稀疏贝叶斯的低复杂度信号波达方向估计方法
CN111580042B (zh) 一种基于相位优化的深度学习测向方法
CN114726385B (zh) 基于功率估计的卫星导航接收机空域抗干扰方法
CN110749855B (zh) 一种基于协方差域零化的均匀线阵波达方向估计方法
CN114879133A (zh) 一种多径和高斯色噪声环境下的稀疏角度估计方法
CN109298382A (zh) 一种基于期望极大算法的非均匀直线阵波达方向角估计方法
CN115563444A (zh) 一种信号重构方法、装置、计算机设备及存储介质
CN113466782B (zh) 一种基于深度学习(dl)的互耦校正doa估计方法
KR20190001170A (ko) 신호의 도래각을 추정하는 방법 및 장치
CN114417917A (zh) 一种未知互耦条件下直接定位方法
WO2010066306A1 (en) Apparatus and method for constructing a sensor array used for direction of arrival (doa) estimation
CN112399366A (zh) 基于Hankel矩阵及WKNN方差提取的室内定位法
CN112087235A (zh) 基于伪逆感知字典的稀疏度自适应doa估计方法及系统
CN116226611A (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