CN111665469B - 一种基于空间时频分布的水下多径信号参数估计方法 - Google Patents

一种基于空间时频分布的水下多径信号参数估计方法 Download PDF

Info

Publication number
CN111665469B
CN111665469B CN202010530393.2A CN202010530393A CN111665469B CN 111665469 B CN111665469 B CN 111665469B CN 202010530393 A CN202010530393 A CN 202010530393A CN 111665469 B CN111665469 B CN 111665469B
Authority
CN
China
Prior art keywords
time
frequency
self
matrix
source
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
CN202010530393.2A
Other languages
English (en)
Other versions
CN111665469A (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.)
Zhejiang University ZJU
Original Assignee
Zhejiang University ZJU
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 Zhejiang University ZJU filed Critical Zhejiang University ZJU
Priority to CN202010530393.2A priority Critical patent/CN111665469B/zh
Publication of CN111665469A publication Critical patent/CN111665469A/zh
Application granted granted Critical
Publication of CN111665469B publication Critical patent/CN111665469B/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/802—Systems for determining direction or deviation from predetermined direction
    • G01S3/8027—By vectorial composition of signals received by plural, differently-oriented transducers
    • 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
    • G01S5/00—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations
    • G01S5/18—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations using ultrasonic, sonic, or infrasonic waves
    • G01S5/20—Position of source determined by a plurality of spaced direction-finders
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02D—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN INFORMATION AND COMMUNICATION TECHNOLOGIES [ICT], I.E. INFORMATION AND COMMUNICATION TECHNOLOGIES AIMING AT THE REDUCTION OF THEIR OWN ENERGY USE
    • Y02D30/00—Reducing energy consumption in communication networks
    • Y02D30/70—Reducing energy consumption in communication networks in wireless communication networks

Abstract

本发明公开了一种基于空间时频分布的水下多径信号参数估计方法，该方法将传统的阵列信号处理技术与时频分布相结合以提高系统的信号处理能力。首先，采用二次时频分布来减少交叉项和提高时频分辨率。其次，通过构造扩展的二维空间时频分布矩阵实现在被动场景下未知多径信号数目的精确估计。再次，根据信号在时频域的能量分布，引入向前向后平滑的思想以保证脊线的局部平滑性，利用属于同一多径信号的自源时频点的主特征向量相同的特性，获得每条多径信号对应脊线的时频点集合。最后，根据脊线检测结果可在欠定条件下实现多径信号瞬时频率和方位的准确估计。

Description

一种基于空间时频分布的水下多径信号参数估计方法
技术领域
本发明涉及水下多径信号参数估计技术领域，特别涉及一种基于空间时频分布的水下多径信 号参数估计方法。
背景技术
水下方位估计技术的研究在海洋信息的获取、环境监测、资源开发、水下灾害预警、海洋权 益维护与安全防卫等方面都发挥着重要作用。
由于海洋水声信道是一个时变-空变-频变的随机多径传输信道，其存在背景噪声强、干扰源多、 可用带宽窄、传输时延大、传播损失大、多途效应以及多普勒效应等严重问题，是目前所知开发难 度最大的无线通信信道。水声信号经过海洋水声信道的传输之后会受到大量的反射和折射，接收端 除了能接收到信源发射的直达路径的信号外，还会接收到其它经过反射和折射的多径信号，这给水 下目标方位估计带来了巨大的挑战。因此，研究如何仅根据阵列接收信号在被动场景下实现多径信 号的方位估计一直是该领域的研究热点和难点，对水下通信、遥测遥控、声呐等海洋信息应用领域 都有极为重要的意义。
方位估计是阵列信号处理领域长期以来的研究热点之一，目前国内外已经涌现出大量优秀的研 究成果，例如比较典型的MUSIC算法、ESPRIT算法、最大似然估计、子空间算法及其各种突破 性改进算法，并取得了非常好的效果，已经逐渐建立了一套完整全面的理论体系。虽然当前对平稳 窄带信号的方位估计研究已经取得了丰硕的成果，但是对非平稳信号的方位估计却鲜有突破，成果 相对较少。传统的方位估计方法在处理被动场景下的多径信号的方位估计存在以下不足：
(1)传统的方位估计方法都假设信号是统计特性不随时间变化的平稳高斯过程，即线性、平 稳、高斯性。然而，水下多径信号本质上是由多个不同频率、不同幅值、不同时延的单模态信号叠 加而成的复杂多模态、非平稳信号。所有的多径信号到达接收端的时延都是不同的，多径信号的频 率成分是随时间而变化的。因此，传统的方位估计方法并不适用于多径信号的方位估计。
(2)传统的方位估计方法利用的是阵列信号的空时统计信息，没有充分挖掘利用信号本身的 时频信息，这必将会损失很多有用的信息，导致其在处理非平稳信号的方位估计问题上存在难以回 避的不足。
(3)传统的方位估计方法通常要求阵元的数目大于源信号数目，即信号环境应满足超定条件。 但是，由于实际的水下环境复杂多变，在未知信源辐射信号和水声传播信道的情况下，信源辐射信 号经多径传播到达接收端后，阵列接收到的多径信号的数目往往是未知且时变的，而接收阵元的数 目有限，因而常常会导致阵列接收到的多径信号的数目大于阵元数目的情况，即欠定条件下的方位 估计，导致传统的方位估计方法并不适用。
此外，实现多径信号参数的准确估计对反演水声传播信道以及协助目标完成定位、跟踪都具有 很强的理论及现实意义。因此，亟需解决在被动场景中实现欠定条件下的多径信号参数估计问题。
发明内容
针对现有方位估计技术在处理被动场景中欠定条件下的多径信号方位估计中存在的不足，本 发明提供一种基于空间时频分布的水下多径信号参数估计方法。本发明采用了时频分析技术，不同 于传统的傅里叶变换，时频分析方法不仅能够有效提取信号的频域信息，同时还能获得频率随时间 的变化特性，进而对信号的时域和频域进行联合分析。本发明将传统的基于空时处理的阵列信号处 理技术与时频分析相结合以获得信号的空间时频分布矩阵。该方法将时域研究对象拓展至时频域， 即将信号在时频域通过时频滤波选取有效的时频点集合，经过时频滤波后的时频点集合可简化为平 稳的窄带信号，利用这些时频点构成的空间时频分布矩阵代替传统的接收信号协方差矩阵，再利用 空间谱估计方法完成方位估计。该方法能够分离在时域和频域难以分离的信号，具有非常好的空间 分辨力和系统的信号处理能力，同时对噪声和干扰也具有非常好的鲁棒性，不仅可以在被动场景下 多径信号数目未知且时变的情况下实现欠定条件下方位的准确估计，而且可以准确估计多径信号的 数目和瞬时频率，易于推广至更简单的主动场景下的多径信号参数估计。
为了解决上述技术问题，本发明采用以下技术方案予以实现：
一种基于空间时频分布的水下多径信号参数估计方法，包括以下步骤：
步骤(1)：在t时刻，使用由M个全向阵元组成的均匀线阵接收水下多径信号，获得t时刻 的阵列输出向量Y(t)＝[y1(t) y2(t) … yM(t)]T，其中，ym(t)表示第m个阵元的输出；
步骤(2)：采用平滑伪魏格纳分布求解阵列输出向量Y(t)的空间时频分布矩阵DYY(t,f)；
步骤(3)：针对DYY(t,f)中的每个采样时间切片设置一个噪声阈值，滤除时频平面上的所有 噪声时频点，将经噪声阈值选择后的时频点保存至集合Ω；
步骤(4)：针对DYY(t,f)进行预白化处理，得到白化空间时频分布矩阵
Figure BDA0002534991740000031
根据
Figure BDA0002534991740000032
设计筛选准则，从集合Ω中筛选出自源时频点，保存至集合ΩA；
步骤(5)：将集合ΩA内的所有自源时频点对应的空间时频分布矩阵组成三阶张量
Figure BDA0002534991740000033
L表示集合ΩA内自源时频点的数量；然后将三阶张量P转化成一个扩展的二维空间时频分布矩阵
Figure BDA0002534991740000034
并通过对矩阵R进行奇异值分解，求解矩阵R的秩，得出信号源的数目N；
步骤(6)：脊线检测：采用前向搜索和后向搜索提取每一条多径信号对应的时频脊线，获得 每一条脊线包含的自源时频点、自源时频点数目、自源时频点对应的瞬时频率；
步骤(7)：根据脊线检测的结果实现多径信号瞬时频率
Figure BDA0002534991740000035
和方位
Figure BDA0002534991740000036
的估计。
与现有技术相比，本发明有益效果是：
本发明将传统的阵列信号处理技术与时频分布相结合以进一步提高系统的信号处理能力。
1)本发明将传统的基于空时处理的阵列信号处理技术与时频分析相结合以获得信号的空间时 频分布矩阵。时频分析技术不同于传统的傅里叶变换，时频分析方法不仅能够有效提取信号的频域 信息，同时还能获得频率随时间的变化特性，进而对信号的时域和频域进行联合分析。该方法将时 域研究对象拓展至时频域，能够分离在时域和频域难以分离的信号，具有非常好的空间分辨力和系 统的信号处理能力，同时对噪声和干扰也具有非常好的鲁棒性。
2)本发明采用二次时频分布中的平滑伪魏格纳分布来减少交叉项，并获得较高的时频分辨率， 同时通过设置噪声阈值和选择自源时频点来降低噪声和交叉项对估计性能的影响。通过构造扩展的 二维空间时频分布矩阵和奇异值分解在被动场景下仅根据阵列接收信号估计出未知多径信号的数 目。
3)本发明充分利用了信号时频表示提供的信号在时频域的能量分布，引入向前向后平滑的思 想以保证脊线的局部平滑性，并利用属于同一多径信号的所有自源时频点具有相同主特征向量的特 性对不同的多径信号分量进行分离，获得每条多径信号对应脊线的时频点集合。
4)本发明根据脊线检测结果实现多径信号瞬时频率和方位的准确估计。该方法以时频分布矩 阵代替传统的协方差矩阵，对具有不同时频特性的多径信号分量进行逐个处理，可在被动场景下多 径信号数目未知且时变的情况下实现欠定条件下的方位估计，有效提高方位估计的精度，并准确估 计多径信号的数目和瞬时频率，而且很容易推广至更简单的主动场景下的多径信号参数估计。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述 中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例， 对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附 图。
图1是本发明所提出的基于空间时频分布的水下多径信号参数估计方法的流程图；
图2是本发明实施例中参考阵元接收到的水下多径信号示意图；
图3是本发明实施例中水下多径信号的时频分布示意图。
图4是本发明实施例中检测到的水下多径信号的脊线示意图。
图5是本发明实施例中基于空间时频分布的水下多径信号瞬时频率估计性能对比。
图6是本发明实施例中基于空间时频分布的水下多径信号方位估计空间谱。
图7是本发明实施例中基于空间时频分布的水下多径信号方位估计性能对比。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述， 显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例， 本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的 范围。
一种基于空间时频分布的水下多径信号参数估计方法，其基本流程如图1所示，主要包括以 下步骤：
S1：在t时刻，使用由M个全向阵元组成的均匀线阵接收水下多径信号，获得t时刻的阵列输 出向量Y(t)＝[y1(t) y2(t) … yM(t)]T，其中，ym(t)表示第m个阵元的输出；
S2：采用平滑伪魏格纳分布求解阵列输出向量Y(t)的空间时频分布矩阵DYY(t,f)；
S3：针对DYY(t,f)中的每个采样时间切片设置一个噪声阈值，滤除时频平面上的所有噪声时 频点，将经噪声阈值选择后的时频点保存至集合Ω；
S4：针对DYY(t,f)进行预白化处理，得到白化空间时频分布矩阵
Figure BDA0002534991740000041
根据
Figure BDA0002534991740000042
设 计筛选准则，从集合Ω中筛选出自源时频点，保存至集合ΩA；
S5：将集合ΩA内的所有自源时频点对应的空间时频分布矩阵组成三阶张量
Figure BDA0002534991740000043
L表 示集合ΩA内自源时频点的数量；然后将三阶张量P转化成一个扩展的二维空间时频分布矩阵
Figure BDA0002534991740000044
并通过对矩阵R进行奇异值分解，求解矩阵R的秩，得出信号源的数目N；
S6：脊线检测：采用前向搜索和后向搜索提取每一条多径信号对应的时频脊线，获得每一条 脊线包含的自源时频点、自源时频点数目、自源时频点对应的瞬时频率；
S7：根据脊线检测的结果实现多径信号瞬时频率
Figure BDA0002534991740000051
和方位
Figure BDA0002534991740000052
的估计。
在本发明的一个具体实施中，对多径信号的接收过程进行了介绍。
在t时刻，使用由M个全向阵元组成的均匀线阵接收多径信号获得阵列输出向量 Y(t)＝[y1(t) y2(t) … yM(t)]T，其中，第m个阵元的输出为ym(t)。具体步骤为：
S11、假设空间中存在K个不相关的远场窄带信源sk(t),(k＝1,2,…,K)，其中，第k个信源对 应的中心频率为fk c。由于水下环境中多径传播条件的影响，假设第k个信源发射的信号经过Nk条 路径入射到由M个全向阵元组成的均匀线阵上，第k个信源对应的Nk条多径信号的入射角分别为
Figure BDA0002534991740000053
其中θk,n∈[-π,π]，(1≤n≤Nk)。假设不同信源的所有多径信号之间都 是不相交的，并且每条多径信号都可以看作是一个虚拟远场窄带信源产生的信号。则K个信源对 应的总的多径信号的数目N表示为：
Figure BDA0002534991740000054
S12、设阵列中阵元的间距为d，以阵元1为参考阵元，考虑到加性高斯白噪声的影响，在t时 刻第m个阵元的输出ym(t)表示为：
Figure BDA0002534991740000055
其中，ak,n表示第k个信源的第n条多径信号对应的衰减因子；Δfk,n表示多普勒频移；Δtk,n表示 传播延迟；
Figure BDA0002534991740000056
表示第m个阵元相对于参考阵元的接收延迟；c是声波在水中的 传播速度；nm(t)是第m个阵元接收到的平稳零均值加性高斯白噪声，即满足E[nm(t)]＝0，信号与 噪声之间不相关。
S13、在t时刻，M×1维的阵列输出向量Y(t)可以表示为：
Figure BDA0002534991740000057
其中，A(Θ)＝[A1(Θ) A2(Θ) … AM(Θ)]T表示阵列的方向矩阵，
其中，Am(Θ)＝[α(θ1) α(θ2) … α(θK))，
Figure BDA0002534991740000061
是第k个 信源对应的阵列方向矢量，
Figure BDA0002534991740000062
表示第k个信源的第n条多径信号对应的中心频率；
Γ＝diag[Γ1 Γ2 … ΓK]，其中，
Figure BDA0002534991740000063
B＝diag[B1 B2 … BK]，其中，
Figure BDA0002534991740000064
S(t)＝[S1(t) S2(t) … SK(t)]T为信号源矢量，其中第k个信源对应的多径信号分别表示为
Figure BDA0002534991740000065
其中，
Figure BDA0002534991740000066
T表示矩阵的转置；
N(t)＝[n1(t) n2(t) … nM(t)]T为阵列的噪声矢量。
在本发明的一个具体实施中，对空间时频分布矩阵DYY(t,f)的求解过程进行了介绍。
采用平滑伪魏格纳分布求解阵列输出向量Y(t)的空间时频分布矩阵DYY(t,f)，具体步骤为：
将阵列输出向量Y(t)的空间时频分布矩阵表示为：
Figure BDA0002534991740000067
其中，t和f分别表示时间和频率索引，(t,f)表示时频点，每一个时频点对应一个的空间时频分布 矩阵DYY(t,f)，l和τ分别表示时移和频移，φ(l,τ)为时频分布的核函数，采用不同的φ(l,τ)得到不 同的时频分布。(·)H表示向量的复共轭转置。
将式(3)带入式(4)可得阵列输出向量Y(t)的空间时频分布矩阵表示为：
DYY(t,f)＝ΨDSS(t,f)ΨH+ΨDSN(t,f)+DNS(t,f)ΨH+DNN(t,f) (5)
其中，对于每一个时频点(t,f)都有DYY(t,f)∈CM×M，其中
Figure BDA0002534991740000068
Ψ＝A(Θ)ΓB包含了信号的空间信息以及传播时延、幅度、衰减等信息，并将信号的时频分布矩阵 中的自源时频分布矩阵和互源时频分布矩阵映射到阵列输出向量的空间时频分布矩阵，空间时频分 布矩阵同时反映了信号的空间和时频信息；DSS(t,f)是多径信号的时频分布矩阵，本发明假设共有 N条多径信号入射到均匀线阵，因此，对于每一个时频点(t,f)，DSS(t,f)为N×N的矩阵，表示为：
Figure BDA0002534991740000069
其中，
Figure BDA0002534991740000071
其对角线元素
Figure BDA0002534991740000072
表示多径信号的自源时频分 布矩阵，非对角线元素
Figure BDA0002534991740000073
表示为互源时频分布矩阵；DSN(t,f)和DNS(t,f)分别表示 信源和噪声之间的互源时频分布矩阵，DNN(t,f)表示噪声的时频分布矩阵。假设信号和噪声分量 之间不相关，且噪声均值为零，因此可以得出信号和噪声向量之间的互源时频分布矩阵为零，因此， 阵列输出的空间时频分布矩阵表示为：
Figure BDA0002534991740000074
其中，δ2为白噪声方差，IM表示单位矩阵。
在本发明的一个具体实施中，对时频点的去噪过程和自源时频点的筛选过程进行了介绍。
S3、设置噪声阈值。
为降低噪声的影响，本实施例通过设置噪声阈值，滤除DYY(t,f)中噪声对应的时频点，并将 DYY(t,f)中剩余信号能量足够大的时频点保存在集合Ω。具体步骤为：
信号在时频平面中的能量相对集中，而噪声的能量较低，且通常均匀分布在整个时频平面上。 因此，可以通过为DYY(t,f)中的每个采样时间切片ti设置一个噪声阈值来判定每个时频点是否为噪 声点，进而滤除时频平面上的所有噪声时频点(t,f)，仅保留信号能量足够大的时频点(t,f)。
Figure BDA0002534991740000075
其中，||·||F为Frobenius范数，i和j分别是时间和频率索引，ε1为门限值，与噪声水平有关(在信 噪比SNR＝10dB时，一般取0.05)，集合Ω表示经噪声阈值选择后的时频点(t,f)的集合，通过去 除噪声点可以有效增强方位估计算法对噪声的鲁棒性。
S4、自源时频点的选择。
非平稳信号方位估计的性能在很大程度上取决于如何准确地选取信号时频点，阵列接收信号 的空间时频分布中包括由信号的自源时频分布产生的自源时频点和信号间的互源时频分布产生的 互源时频点，互源时频点会严重影响信号方位估计的性能。因此，本发明需要先对信号的时频分布 进行优化，以去除信号间的互源时频点，选取信号对应的自源时频点。
首先，本实施例通过定义一个白化矩阵W对阵列输出向量Y(t)的空间时频分布矩阵DYY(t,f) 进行预白化处理，然后通过设定特定的准则从集合Ω中选择自源时频点(t,f)a，去除互源时频点 (t,f)c，并将选择后的自源时频点保存至集合ΩA。具体步骤为：
S41、阵列接收信号Y(t)的协方差矩阵表示为：
Figure BDA0002534991740000081
其中，RS表示信号源矢量的协方差矩阵。定义一个N×M维的白化矩阵W＝∑-1/2VH。其中，∑和V 分别表示对RY进行特征分解后与信号相对应的特征值和特征向量。对于任何一个白化矩阵W都存 在一个N×N维的酉矩阵U使得U＝WΨ，即(WΨ)(WΨ)H＝UUH＝I，其中，I表示单位矩阵。
S42、如式(10)所示，通过在DYY(t,f)左右两侧分别乘以白化矩阵W得到阵列输出向量Y(t) 的白化空间时频分布矩阵
Figure BDA0002534991740000082
Figure BDA0002534991740000083
S43、对于任意的互源时频点(t,f)c，其相应的白化空间时频分布矩阵满足如下条件。
Figure BDA0002534991740000084
其中，下标c表示该时频点为互源时频点，tr{·}表示矩阵的迹。
S44、按照以下准则从集合Ω中选择自源时频点(t,f)a，下标a表示该时频点为自源时频点。
Figure BDA0002534991740000085
其中，阈值ε2为小于1经验阈值(通常当SNR＝10dB时，ε2＝0.85)，选择后的自源时频点保存至 集合ΩA。
在本发明的一个具体实施中，对多径信号数目的估计过程进行了介绍。
S5、多径信号数目的估计。
本实施例将集合A内的所有自源时频点对应的空间时频分布矩阵组成三阶张量
Figure BDA0002534991740000086
然后将三阶张量P转化成一个扩展的二维空间时频分布矩阵
Figure BDA0002534991740000087
并通过对矩阵R进行奇异 值分解估计出信号源的数目N。具体步骤为：
S51、假设自源时频点集合ΩA内总共包含L个时频点，即：
ΩA＝{(t,f)l|1≤l≤L} (14)
其中，下标l表示时频点索引，L个自源时频点对应的L个空间时频分布矩阵表示分别为：
DYY(t,f)l＝ΨDSS(t,f)lΨH,(1≤l≤L) (15)
S52、利用L个自源时频点对应的L个空间时频分布矩阵构造一个三阶张量
Figure BDA0002534991740000088
其中， P的第(i,j,l)个元素Pijl表示为：
Pijl＝[DYY(k,f)l]ij,(1≤i,j≤M,1≤l≤L) (16)
S53、根据三阶张量P构造一个扩展的二维空间时频分布矩阵
Figure BDA0002534991740000091
且满足R(i-1)M+j,l＝Pijl。
S54、定义矩阵
Figure BDA0002534991740000092
Cl,n＝[DSS(k,f)l]nn，则矩阵R可转化为：
R＝(Ψ⊙Ψ*)CT (17)
式中，
Figure BDA0002534991740000093
其中，(·)*、⊙和
Figure BDA0002534991740000094
分别表示复数的共 轭、Khatri-Rao乘积以及Kronecker乘积。
S55、求解矩阵R的秩估计接收信号Y(t)中总的多径信号的数目N。
在本发明的一个具体实施中，对脊线检测过程进行了介绍。
S6、脊线检测。
实际上，在时频聚集性较高的时频平面上，信号的能量主要集中在被称为脊线的区域。因此， 时频分布的脊线不仅能够体现信号瞬时频率的变化情况，而且包含了表征信号特征的关键信息。而 且，本发明假设所有的多径信号在时频域上都是不相交的，即每一条多径信号分量都对应着一条脊 线。因此，只要检测出属于每一条多径信号的脊线，就可以获得属于每个多径信号的时频点集合， 进而实现多径信号瞬时频率和方位的准确估计。
首先，该发明充分考虑了信号时频表示所提供的信号在时频域的能量分布，并引入向前向后 平滑的思想，在保证探测到的脊线在相邻时刻的局部平滑性的基础上寻找所有的局部极值点。其次， 利用属于同一多径信号的所有自源时频点应具有相同主特征向量的性质，获得每一个多径信号对应 脊线的时频点集合。具体步骤为：
S61、参数设置。本发明假设系统的采样频率为fs，在时间轴，采样间隔为Δt＝1/fs，信号总 的采样数目为T，在频率轴，谱线间隔为Δf＝fs/N，最大搜索频率Fmax＝fs/2，假设在t时刻包含 的多径信号的数目为Nt，其中，0≤t≤T。
构建一个零矩阵N_tf＝zeros(1,N)用来保存每条脊线含有的自源时频点的个数；
构建N个空集合Ωn＝[],(1≤n≤N)用来保存每条脊线包含的所有自源时频点；
构建N个空集合Fn＝[],(1≤n≤N)用来保存每条脊线包含的所有自源时频点对应的瞬时频率；
构建N个空集合En＝[],(1≤n≤N)用来保存集合Fn中包含的所有自源时频点(t,f)对应的单位 化特征向量
Figure BDA0002534991740000095
S62、初始化。
S621、本发明初始化时间切片t＝T/2，并假设多径信号在中间时刻t能够取到最大的源数目N， 即Nt＝N；
S622、求解DYY(t,f)在t时刻的Nt个极值点对应的瞬时频率ft,n，其中(0≤f≤Fmax)：
Figure BDA0002534991740000101
保存相应的时频点和瞬时频率，即(t,ft,n)∈Ωn和ft,n∈Fn，并使N_tf(n)＝N_tf(n)+1；
S623、分别对Nt个极值点对应的自源时频点(t,f)的空间时频分布矩阵DYY(t,ft,n)，(1≤n≤Nt) 进行特征分解：
DYY(t,ft,n)＝V(t,ft,n)Λ(t,ft,n)V(t,ft,n)H (19)
其中，V(t,ft,n)是一个M×M的酉矩阵，它的列向量表示特征向量，Λ(t,ft,n)为M×M的对角矩阵， 其对角元素为与所有特征向量对应的特征值。设最大的特征值对应的特征向量分别表示为v(t,ft,n)， 对每一个特征向量进行单位化处理。
Figure BDA0002534991740000102
并保存单位化之后的特征向量，即
Figure BDA0002534991740000103
S63、前向搜索。
S631、对每一个时刻h，其中h＝t-1,t-2,…,1，首先判断h+1时刻所有极值点对应的自源时 频点是否为边缘时频点，进而估计h时刻对应的多径信号的数目Nh。
若任意的f∈(fh+1,n-ΔF,fh+1,n+ΔF)，其中，1≤n≤Nh+1，ΔF为允许搜索的频率范围，对应时 频点(h+1,f)都满足：
DYY(h+1,f)＝0 (21)
则认为该自源时频点(h+1,fh+1,n)为左边缘时频点，h时刻的目标数Nh＝Nh+1-1。
S632、根据式(18)，求解在h时刻DYY(h,f)，(1≤f≤Fmax)的Nh个极值点对应的频率fh,j，(1≤j≤Nh)。
S633、根据式(19)和(20)求在h时刻第j(1≤j≤Nh)个极值点对应的自源时频点的空间时 频分布矩阵DYY(h,fh,j)对应的单位化特征向量为
Figure BDA0002534991740000111
与h+1时刻的每一个特征向量
Figure BDA0002534991740000112
进行比较，其中(1≤n≤Nh+1)，若满足：
Figure BDA0002534991740000113
则该极值点对应的自源时频点属于第n条脊线，即(h,fh,j)∈Ωn和fh,j∈Fn，保存单位化后的特征 向量
Figure BDA0002534991740000114
且第n条脊线中包含的自源时频点数目N_tf(n)＝N_tf(n)+1。
S64、后向搜索。
S641、对每一个时刻h，其中h＝t+1,t+2,…,T，根据步骤S631判断h-1时刻所有极值点对 应的时频点是否为边缘时频点，并估计h时刻对应的多径信号的数目Nh。
若任意的f∈(fh-1,n-ΔF,fh-1,n+ΔF)，其中，1≤n≤Nh-1，对应时频点(h-1,f)都满足：
DYY(h-1,f)＝0 (23)
则认为该自源时频点(h-1,fh-1,n)为右边缘时频点，h时刻的目标数Nh＝Nh-1-1。
S642、根据式(18)，求解在h时刻DYY(h,f)，(1≤f≤Fmax)的Nh个极值点对应的频率fh,j，(1≤j≤Nh)。
S643、根据式(19)和(20)求在h时刻第j(1≤j≤Nh)个极值点对应的自源时频点的空间时 频分布矩阵DYY(h,fh,j)对应的单位化特征向量为
Figure BDA0002534991740000115
与h-1时刻的每一个特征向量
Figure BDA0002534991740000116
进行比较，其中(1≤n≤Nh-1)，若满足：
Figure BDA0002534991740000117
则该极值点对应的自源时频点属于第n条脊线，即(h,fh,j)∈Ωn和fh,j∈Fn，保存单位化后的特征 向量
Figure BDA0002534991740000118
且第n条脊线中包含的自源时频点数目N_tf(n)＝N_tf(n)+1。
S65、参数输出。输出N_tf、Ωn,(1≤n≤N)以及Fn,(1≤n≤N)。
在本发明的一个具体实施中，对多径信号参数估计过程进行了介绍。
S7、多径信号参数估计。
本发明根据脊线检测的结果进行多径信号的参数估计，主要包括信号的瞬时频率估计
Figure BDA0002534991740000119
以及 方位估计
Figure BDA00025349917400001110
具体步骤为：
S71、瞬时频率估计。根据每条脊线包含的所有自源时频点对应的瞬时频率集合Fn,(1≤n≤N)， 可估计第n条多径信号对应的瞬时频率为：
Figure BDA0002534991740000121
S72、方位估计。计算每条脊线对应的时频点集合Ωn(1≤n≤N)中所有时频点所对应的空间时 频分布矩阵的平均值，然后使用TF-music算法实现多径信号的方位估计。具体步骤为：
S721、构建与每一条脊线相对应的空间时频分布矩阵。计算第n条脊线对应的时频点集合 Ωn(1≤n≤N)中所有时频点对应的空间时频分布矩阵的平均值。
Figure BDA0002534991740000122
其中，N_tf(n)表示第n条脊线对应的时频点集合Ωn中的时频点的个数。
S722、对第n条脊线对应的平均空间时频分布矩阵
Figure BDA0002534991740000123
进行特征分解。
Figure BDA0002534991740000124
其中，∑S为大特征值组成的对角阵，∑N为小特征值组成的对角阵，US为信号子空间，UN为噪声 子空间。
S723、定义空间谱函数：
Figure BDA0002534991740000125
其中，a(θ)是阵列的方向矢量，与信号的频率和波达方向有关，信号频率确定后，就只与波达方 向有关；
第n条多径信号的入射角的估计值为：
Figure BDA0002534991740000126
其中，argmax(·)是argmax函数。
实施例
为验证所提基于空间时频分布的水下多径信号参数估计方法的有效性，本发明进行了仿真验 证。使用均方根误差来评估参数的估计性能。估计误差表示为：
Figure BDA0002534991740000131
其中，S表示蒙特卡罗仿真的次数，P表示参数的个数，ξs,p表示第p个参数在第s次仿真中的实 际值，
Figure BDA0002534991740000132
表示ξs,p的估计值，RMSE仿真结果均通过500次蒙特卡罗实验得到。
仿真条件：假设均匀线阵的阵元数M＝8，阵元间距d＝λ/2，λ为波长，系统的采样频率 fs＝1024Hz，采样数目T＝1024。假设水中共存在两个不相关的窄带信源，信源的中心频率分别 为200Hz和400Hz。在所有仿真中仅考虑水底和水面的一次反射和二次反射，即每个信源经过五 条路径传播入射到阵列上，所有多径的参数设置如表1所示。
表1多径信号参数设置
参数 信源1 信源1
入射角度 -60,-50,-40,-30,-20 30,40,50,60,70
多普勒频移 0,-15,15,-30,30 0,-15,15,-30,30
相对传播时延 0,0.0135,0.0227,0.0750,0.1116 0,0.0265,0.0100,0.1141,0.0495
信道衰减系数 1,0.9980,0.9883,0.9808,0.9716 1,0.9937,0.9976,0.9731,0.9882
如图2所示为在参考阵元处接收到的多径信号，图中的横坐标纵坐标分别表示时间和信号幅 值。图3为经过去燥和单源点选择后水下多径信号对应的时频分布图，图中的横坐标纵坐标分别表 示时间和频率，该图可以清楚地反映多径信号的频率成分随时间的变化情况。
如图4所示为采用脊线检测算法检测到的水下多径信号的脊线示意图，图中的横坐标纵坐标 分别表示时间和频率，从图中可以看出本发明提出的脊线检测方法可以准确地提取信号的时频脊 线，为水下多径信号的瞬时频率估计和方位估计提供保障。
如图5所示为本发明所提方法的瞬时频率估计与基于短时傅里叶变换的瞬时频率估计的 RMSE随信噪比的变化情况，图中的横坐标纵坐标分别表示时间和频率估计误差，从图中可以看出， 本发明所提的方法能够非常准确地估计信号的瞬时频率，且随着信噪比的增加，估计精度越来越高。
如图6所示比较了在角度间隔为0.1°时本发明所提方法分别采用估计瞬时频率和真实瞬时频 率进行方位估计时的空间谱。图中的横坐标纵坐标分别表示时间和空间谱，从图中可以看出，本发 明所提方法能够准确地识别出每条多径信号的角度，并且在使用估计瞬时频率进行方位估计时也能 实现非常好的估计性能。
图7比较了本发明所提方位估计方法和基于短时傅里叶变换的方位估计方法的RMSE随信噪 比的变化情况，图中的横坐标纵坐标分别表示时间和方位估计误差，与其它方位估计算法相比，本 发明所提方法的估计精度有明显提升，且所需的阵元数目大幅度减少，即可实现在欠定条件下的方 位估计，这是因为本发明所提的方法利用脊线检测算法分别检测出属于每一条多径信号的脊线，获 得属于每个多径信号的时频点集合，然后对每一个多径信号分别进行方位估计。综上，本发明提出 的基于空间时频分布的水下多径信号参数估计方案可以实现在被动场景中的欠定条件下的多径信 号参数的准确估计，并且该方法对噪声和干扰具有较强的鲁棒性。
步骤S11中，第k个远场窄带信源sk(t),(k＝1,2,…,K)通常表示为一个单频脉冲信号：
Figure BDA0002534991740000141
其中，Ak表示信号的幅度(不失一般性，一般设其幅度为1)，
Figure BDA0002534991740000142
表示信号的中心频率。
步骤S2中，信号s(t)对应的自平滑伪魏格纳分布的离散形式定义为：
Figure BDA0002534991740000143
其中，(·)*表示复数的共轭，l和τ分别表示时间和频率索引，φ(l,τ)为时频分布的核函数，给出不 同的φ(l,τ)得到不同的时频分布。
信号s1(t)和s2(t)之间的互平滑伪魏格纳分布的离散形式定义为：
Figure BDA0002534991740000144
平滑伪魏格纳分布是一种二次时频分布，可以通过对时间和频率进行同时加窗进行平滑处理，可以 在抑制交叉项的同时获得较好的时频聚集性。
步骤S4中，若一个时频点(t,f)a仅聚集了信号自源时频分布矩阵
Figure BDA0002534991740000145
的能量，则时频点 (t,f)a称为信号si(t)的自源时频点，下标a表示该时频点为自源时频点。自源时频点所对应的信号 空间时频分布矩阵DSS(t,f)为一个对角矩阵。
若一个时频点(t,f)c仅聚集了信号互源时频分布矩阵
Figure BDA0002534991740000146
(i≠j)的能量，则时频点(t,f)c称为信号si(t)和sj(t)的互源时频点，下标c表示该时频点为互源时频点。互源时频点所对应的信号 空间时频分布矩阵DSS(t,f)为一个对角线元素为零、非对角线元素不为零的矩阵。
步骤S5中，当多径信号的数目满足N≤min(M2,L)时，Ψ⊙Ψ*和C是列满秩的，即满足 rank(Ψ⊙Ψ*)＝N和rank(CT)＝N。因此，矩阵R也是列满秩的，矩阵R的秩等于多径信号的数目 N，即rank(R)＝N。
步骤S6中，本发明假设所有信号源对应的所有多径信号在时频域都是不重合的，即一个自源 时频点仅与一个多径信号相对应。假设所有多径信号对应的自源时频点集合表示为ΩA，第n个多 径信号对应的自源时频点集合表示为Ωn(1≤n≤N)，其中，
Figure BDA0002534991740000151
且满足
Figure BDA0002534991740000152
(m≠n)，其中m和n分别表示不同的多径信号的索引。因此，对于时频平面中的任何一个自源时 频点对应的时频分布矩阵都满足：
Figure BDA0002534991740000153
其中，
Figure BDA0002534991740000154
表示第n个多径信号在自源时频点(t,f)i处的空间时频分布矩阵，它可以看作是矩 阵DYY(t,f)i进行特征分解后的主特征值，αn可以看作是主特征值对应的特征向量，而其它特征值 为零，因此，DYY(t,f)i是一个秩为1的矩阵，其特征分解可以表示为：
DYY(t,f)i＝V(t,f)iΛ((t,f)iV(t,f)i (36)
其中，V(t,f)i是一个M×M的酉矩阵，它的列向量v(t,f)i表示特征向量，Λ((t,f)i为M×M的对角 矩阵，其对角元素为与所有特征向量对应的特征值。
由以上分析可知，对于集合Ωn中的任何两个自源时频点(t,f)i和(t,f)j，与矩阵DYY(t,f)i和 DYY(t,f)j相对应的主特征向量是相同的，这一性质可以为寻找属于每个多径信号的所有自源时频 点提供有力的工具。
当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的普通 技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。

Claims (6)

1.一种基于空间时频分布的水下多径信号参数估计方法，其特征在于，包括以下步骤：
步骤(1)：在t时刻，使用由M个全向阵元组成的均匀线阵接收水下多径信号，获得t时刻的阵列输出向量Y(t)＝[y1(t) y2(t)…yM(t)]T，其中，ym(t)表示第m个阵元的输出；
步骤(2)：采用平滑伪魏格纳分布求解阵列输出向量Y(t)的空间时频分布矩阵DYY(t,f)；
步骤(3)：针对DYY(t,f)中的每个采样时间切片设置一个噪声阈值，滤除时频平面上的所有噪声时频点，将经噪声阈值选择后的时频点保存至集合Ω；
步骤(4)：针对DYY(t,f)进行预白化处理，得到白化空间时频分布矩阵
Figure FDA0003716689830000011
根据
Figure FDA0003716689830000012
设计筛选准则，从集合Ω中筛选出自源时频点，保存至集合ΩA；
步骤(5)：将集合ΩA内的所有自源时频点对应的空间时频分布矩阵组成三阶张量
Figure FDA0003716689830000013
L表示集合ΩA内自源时频点的数量；然后将三阶张量P转化成一个扩展的二维空间时频分布矩阵
Figure FDA0003716689830000014
并通过对矩阵R进行奇异值分解，求解矩阵R的秩，得出信号源的数目N；
步骤(6)：脊线检测：采用前向搜索和后向搜索提取每一条多径信号对应的时频脊线，获得每一条脊线包含的自源时频点、自源时频点数目、自源时频点对应的瞬时频率；
步骤(7)：根据脊线检测的结果实现多径信号瞬时频率
Figure FDA0003716689830000015
和方位
Figure FDA0003716689830000016
的估计。
2.根据权利要求1所述的基于空间时频分布的水下多径信号参数估计方法，其特征在于所述的步骤(1)具体为：
在t时刻，第m个阵元的输出ym(t)表示为：
Figure FDA0003716689830000017
其中，ak,n表示第k个信源的第n条多径信号对应的衰减因子；Δfk,n表示多普勒频移；Δtk,n表示传播延迟；τk,n表示第m个阵元相对于参考阵元的接收延迟；c是声波在水中的传播速度；nm(t)是第m个阵元接收到的平稳零均值加性高斯白噪声；sm(t)是第m个阵元接收到的多径信号；K是远场窄带信源的数量；Nk表示第k个信源对应的多径信号的数目；
Figure FDA0003716689830000018
是第k个信源对应的中心频率；
阵列输出向量Y(t)表示为：
Y(t)＝[y1(t) y2(t)…yM(t)]T
＝A(Θ)ΓBS(t)+N(t)
其中，A(Θ)表示阵列的方向矩阵，Γ和B表示对角矩阵，S(t)为t时刻的信号源矢量，T表示矩阵的转置，N(t)为阵列的噪声矢量。
3.根据权利要求1所述的基于空间时频分布的水下多径信号参数估计方法，其特征在于所述的步骤(2)具体为：
将阵列输出向量Y(t)进行平滑伪魏格纳变换，得到空间时频分布矩阵DYY(t,f)，表示为：
DYY(t,f)＝ΨDSS(t,f)ΨH+δ2IM
其中，t和f分别表示时间和频率的索引，(t,f)表示时频点，每一个时频点对应一个M*M维的空间时频分布矩阵DYY(t,f)；Ψ＝A(Θ)ΓB，A(Θ)表示阵列的方向矩阵，Γ和B表示对角矩阵；(·)H表示向量的复共轭转置，δ2为白噪声方差，IM表示单位矩阵；DSS(t,f)是多径信号的时频分布矩阵，表示为：
Figure FDA0003716689830000021
对角线元素
Figure FDA0003716689830000022
表示多径信号的自源时频分布矩阵，非对角线元素
Figure FDA0003716689830000023
表示为互源时频分布矩阵。
4.根据权利要求1所述的基于空间时频分布的水下多径信号参数估计方法，其特征在于所述的步骤(4)具体为：
对空间时频分布矩阵DYY(t,f)进行预白化处理，得到白化空间时频分布矩阵
Figure FDA0003716689830000024
公式为：
Figure FDA0003716689830000025
其中，W为白化矩阵，U＝WΨ表示酉矩阵，Ψ＝A(Θ)ΓB，A(Θ)表示阵列的方向矩阵，Γ和B表示对角矩阵，(WΨ)(WΨ)H＝UUH＝I，I表示单位矩阵；
根据白化空间时频分布矩阵
Figure FDA0003716689830000026
按照以下准则从集合Ω中选择自源时频点：
Figure FDA0003716689830000027
其中，ε2为阈值，tr{·}表示矩阵的迹，选择后的自源时频点(t,f)保存至集合ΩA。
5.根据权利要求1所述的基于空间时频分布的水下多径信号参数估计方法，其特征在于所述的步骤(6)具体为：
步骤(6.1)：参数设置：
构建一个零矩阵N_tf＝zeros(1,N)用来保存每条脊线含有的自源时频点的个数；
构建N个空集合Ωn,(1≤n≤N)用来保存每条脊线包含的所有自源时频点；
构建N个空集合Fn,(1≤n≤N)用来保存每条脊线包含自源时频点对应的瞬时频率；
初始化时间切片t＝T/2，T表示信号总的采样数目；
步骤(6.2)：求解空间时频分布矩阵DYY(t,f)在t时刻的Nt个极值点对应的瞬时频率ft,n，其中(0≤f≤Fmax)：
Figure FDA0003716689830000031
保存相应的时频点和瞬时频率，即(t,ft,n)∈Ωn和ft,n∈Fn，并使N_tf(n)＝N_tf(n)+1；Nt表示在t时刻包含的多径信号的数目，初始化Nt＝N；
分别对Nt个极值点对应的自源时频点的空间时频分布矩阵DYY(t,ft,n)，(1≤n≤Nt)进行特征分解：
DYY(t,ft,n)＝V(t,ft,n)Λ(t,ft,n)V(t,ft,n)H
其中，V(t,ft,n)是一个M×M的酉矩阵，其列向量表示特征向量，Λ(t,ft,n)为M×M的对角矩阵，对角元素为与所有特征向量对应的特征值，将最大的特征值对应的特征向量表示为v(t,ft,n)，将Nt个最大的特征值对应的特征向量进行单位化处理，得到单位化之后的特征向量，即
Figure FDA0003716689830000032
步骤(6.3)：前向搜索：
对每一个时刻h，其中h＝t-1,t-2,…,1，首先判断h+1时刻所有极值点对应的自源时频点是否为边缘时频点，若任意的f∈(fh+1,n-ΔF,fh+1,n+ΔF)，其中，1≤n≤Nh+1，ΔF为允许搜索的频率范围，对应时频点(h+1,f)都满足：
DYY(h+1,f)＝0
则该自源时频点(h+1,fh+1,n)为左边缘时频点，h时刻的多径信号目标数Nh＝Nh+1-1；
求解在h时刻DYY(h,f)，(1≤f≤Fmax)的Nh个极值点对应的频率fh,j，(1≤j≤Nh)，根据h时刻第j(1≤j≤Nh)个极值点对应的自源时频点的空间时频分布矩阵DYY(h,fh,j)，求解对应的单位化特征向量为
Figure FDA0003716689830000041
与h+1时刻的每一个特征向量
Figure FDA0003716689830000042
进行比较，其中(1≤n≤Nh+1)，若满足：
Figure FDA0003716689830000043
则该极值点对应的自源时频点属于第n条脊线，即(h,fh,j)∈Ωn和fh,j∈Fn，且第n条脊线中包含的自源时频点数目N_tf(n)＝N_tf(n)+1；
步骤(6.4)：后向搜索：
对每一个时刻h，其中h＝t+1,t+2,…,T，判断h-1时刻所有极值点对应的时频点是否为边缘时频点，若任意的f∈(fh-1,n-ΔF,fh-1,n+ΔF)，其中，1≤n≤Nh-1，对应时频点(h-1,f)都满足：
DYY(h-1,f)＝0
则该自源时频点(h-1,fh-1,n)为右边缘时频点，h时刻的目标数Nh＝Nh-1-1；
求解在h时刻DYY(h,f)，(1≤f≤Fmax)的Nh个极值点对应的频率fh,j，(1≤j≤Nh)，根据h时刻第j(1≤j≤Nh)个极值点对应的自源时频点的空间时频分布矩阵DYY(h,fh,j)，求解对应的单位化特征向量为
Figure FDA0003716689830000044
与h-1时刻的每一个特征向量
Figure FDA0003716689830000045
进行比较，其中(1≤n≤Nh-1)，若满足：
Figure FDA0003716689830000046
则该极值点对应的自源时频点属于第n条脊线，即(h,fh,j)∈Ωn和fh,j∈Fn，且第n条脊线中包含的自源时频点数目N_tf(n)＝N_tf(n)+1；
步骤(6.5)：输出Ωn、N_tf(n)、以及Fn,(1≤n≤N)，得到每一条脊线包含的自源时频点、自源时频点数目、自源时频点对应的瞬时频率。
6.根据权利要求1所述的基于空间时频分布的水下多径信号参数估计方法，其特征在于所述的步骤(7)具体为：
步骤(7.1)：瞬时频率
Figure FDA0003716689830000047
估计：
根据每条脊线包含的所有自源时频点对应的瞬时频率集合Fn,(1≤n≤N)，得到第n条多径信号对应的瞬时频率的估计值
Figure FDA0003716689830000048
为：
Figure FDA0003716689830000049
其中，N_tf(n)为第n条脊线包含的自源时频点数目，fi是第n条脊线包含的第i个自源时频点对应的瞬时频率，Fn是第n条脊线包含的所有自源时频点对应的瞬时频率集合；
步骤(7.2)：方位
Figure FDA0003716689830000051
估计：
计算每条脊线包含的所有自源时频点对应的空间时频分布矩阵的平均值，公式为：
Figure FDA0003716689830000052
其中，
Figure FDA0003716689830000053
是第n条脊线对应的平均空间时频分布矩阵，Ωn是第n条脊线包含的所有自源时频点集合；
对
Figure FDA0003716689830000054
进行特征分解，公式为：
Figure FDA0003716689830000055
其中，ΣS为大特征值组成的对角阵，ΣN为小特征值组成的对角阵，US为信号子空间，UN为噪声子空间；(·)H表示向量的复共轭转置；
采用TF-music算法进行方位估计，定义空间谱函数：
Figure FDA0003716689830000056
其中，a(θ)是阵列的方向矢量；
第n条多径信号的入射角的估计值
Figure FDA0003716689830000057
为：
Figure FDA0003716689830000058
其中，argmax(·)是argmax函数。
CN202010530393.2A 2020-06-11 2020-06-11 一种基于空间时频分布的水下多径信号参数估计方法 Active CN111665469B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010530393.2A CN111665469B (zh) 2020-06-11 2020-06-11 一种基于空间时频分布的水下多径信号参数估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010530393.2A CN111665469B (zh) 2020-06-11 2020-06-11 一种基于空间时频分布的水下多径信号参数估计方法

Publications (2)

Publication Number Publication Date
CN111665469A CN111665469A (zh) 2020-09-15
CN111665469B true CN111665469B (zh) 2022-08-23

Family

ID=72386683

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010530393.2A Active CN111665469B (zh) 2020-06-11 2020-06-11 一种基于空间时频分布的水下多径信号参数估计方法

Country Status (1)

Country Link
CN (1) CN111665469B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114866159B (zh) * 2022-04-01 2023-04-21 华南理工大学 一种多分量线性调频信号时频分析方法

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2008094701A1 (en) * 2007-01-31 2008-08-07 Signal Labs, Inc. System and methods for multistep target detection and parameter estimation
CN105842656A (zh) * 2016-05-31 2016-08-10 黑龙江工程学院 基于雅克比旋转联合对角化的空时频方位估计方法
WO2018049595A1 (zh) * 2016-09-14 2018-03-22 深圳大学 一种基于交替方向乘子法的稳健稀疏恢复stap方法及其系统
CN108680910A (zh) * 2018-05-15 2018-10-19 河海大学 基于波形认知的调频广播外辐射源雷达目标检测方法
CN109633558A (zh) * 2018-10-25 2019-04-16 上海无线电设备研究所 一种基于极化时频分布的波达方向估计算法
CN110243605A (zh) * 2019-06-10 2019-09-17 苏州大学 多源时频脊线提取方法
CN110764053A (zh) * 2019-10-22 2020-02-07 浙江大学 一种基于水下传感器网络的多目标被动定位方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR101413229B1 (ko) * 2013-05-13 2014-08-06 한국과학기술원 방향 추정 장치 및 방법

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2008094701A1 (en) * 2007-01-31 2008-08-07 Signal Labs, Inc. System and methods for multistep target detection and parameter estimation
CN105842656A (zh) * 2016-05-31 2016-08-10 黑龙江工程学院 基于雅克比旋转联合对角化的空时频方位估计方法
WO2018049595A1 (zh) * 2016-09-14 2018-03-22 深圳大学 一种基于交替方向乘子法的稳健稀疏恢复stap方法及其系统
CN108680910A (zh) * 2018-05-15 2018-10-19 河海大学 基于波形认知的调频广播外辐射源雷达目标检测方法
CN109633558A (zh) * 2018-10-25 2019-04-16 上海无线电设备研究所 一种基于极化时频分布的波达方向估计算法
CN110243605A (zh) * 2019-06-10 2019-09-17 苏州大学 多源时频脊线提取方法
CN110764053A (zh) * 2019-10-22 2020-02-07 浙江大学 一种基于水下传感器网络的多目标被动定位方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
A Two-step 2-D DOA Estimation Algorithm Via Underwater Acoustic Vector-Sensor Array;Xueyan Han 等;《2019 IEEE International Conference on Signal Processing, Communications and Computing (ICSPCC)》;20191231;全文 *
Time-frequency Diffusion Analysis for DOA Estimation of Tropospheric Scatter Signal;Jing Wang 等;《2016 Progress In Electromagnetic Research Symposium (PIERS)》;20160811;全文 *
基于时频子空间的DOA估计;应武等;《航天电子对抗》;20040630(第03期);全文 *

Also Published As

Publication number Publication date
CN111665469A (zh) 2020-09-15

Similar Documents

Publication Publication Date Title
Krishnaveni et al. Beamforming for direction-of-arrival (DOA) estimation-a survey
CN108375763B (zh) 一种应用于多声源环境的分频定位方法
CN110045323B (zh) 一种基于矩阵填充的互质阵稳健自适应波束形成算法
CN109459744B (zh) 一种实现多干扰抑制的稳健自适应波束形成方法
CN110320490B (zh) 一种无直达信号条件下的无线电波达方向估计方法
CN109541548B (zh) 一种基于匹配场的空气声呐定位方法
CN109765521B (zh) 一种基于子阵划分的波束域成像方法
CN107576931B (zh) 一种基于协方差低维度迭代稀疏重构的相关/相干信号波达方向估计方法
CN109407055A (zh) 基于多径利用的波束形成方法
CN113064147A (zh) 一种低信噪比下新型匹配场被动定位方法
CN114726385B (zh) 基于功率估计的卫星导航接收机空域抗干扰方法
Durrani et al. Eigenfilter approaches to adaptive array processing
CN111665469B (zh) 一种基于空间时频分布的水下多径信号参数估计方法
CN114884841A (zh) 基于高阶统计和非均匀阵列的欠定参数联合估计方法
Liu et al. Target localization in high-coherence multipath environment based on low-rank decomposition and sparse representation
Okane et al. Resolution improvement of wideband direction-of-arrival estimation" Squared-TOPS"
Yu et al. A robust minimum variance beamformer with new constraint on uncertainty of steering vector
Lefort et al. Sub-antenna processing for coherence loss in underwater direction of arrival estimation
CN113805139A (zh) 基于聚焦变换的宽带信号稀疏表示波达方向估计方法
CN110320492B (zh) 一种基于导向矢量对称特性的低复杂度doa估计方法
Harrison et al. Robust source localization in an acoustic waveguide
CN108181611B (zh) 基于子空间的压缩感知高分辨阵列处理方法
Sun et al. Robust adaptive beamforming for coprime array based on interference power estimation via subspace
Jiang et al. The DOA estimation for a mixture of uncorrelated and coherent sources via decomposing the coprime array
Mao et al. An Effective Algorithm for Direction-of-Arrival Estimation of Coherent Signals with ULA

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