CN111563233B - 基于非线性贝叶斯理论的水声信号频散特征被动反演方法 - Google Patents

Abstract

本发明涉及一种基于非线性贝叶斯理论的简正波频散曲线被动反演方法，从统计角度定性和定量的分析反演结果的不确定性，为反演结果的可信度提供参考。首先利用噪声干涉仪技术与简正波分离技术，被动提取各阶水声信号的频散特征，然后利用自适应单纯形法模拟退火算法开展海底参数的反演，给出各个反演参数的边缘概率分布，为反演结果的准确性提供评判参考。

Description

基于非线性贝叶斯理论的水声信号频散特征被动反演方法

技术领域

本发明属于水声信号频散特征被动反演方法，涉及一种基于非线性贝叶斯理论的水声信号频散特征被动反演方法，将贝叶斯理论应用于水声信号频散特征提取并被动反演中，从而为被动反演结果的准确性提供评判依据。

背景技术

文献《Characterizing the seabed by using noise interferometry and timewarping》给出了一种基于海洋环境噪声信号的水声信号频散特征被动提取方法，并将其应用于地声参数反演中。首先利用海洋环境噪声信号提取两点之间的时空互相关函数，通过对互相关函数做warping变换分离不同阶模态信号，进而得到水声信号各阶简正波的频散特征，最后利用频散特征对地声参数进行了反演。论文虽然给出了海底参数反演结果，但是反演的基底参数不理想，与主动反演结果差异较大，此外，该论文没有对参数的不确定性进行分析，无法评价反演结果的可信度。

发明内容

要解决的技术问题

为了避免现有技术的不足之处，本发明提出一种基于非线性贝叶斯理论的水声信号频散特征被动反演方法。

技术方案

一种基于非线性贝叶斯理论的水声信号频散特征被动反演方法，其特征在于步骤如下：

步骤1、利用噪声干涉仪技术与简正波分离技术，被动提取各阶水声信号的频散特征：

根据噪声干涉仪技术，空间两点噪声场互相关函数的导数用两点之间的因果格林函数与反因果格林函数之和表征：

其中C(τ,r₁,r₂)，G(τ,r₁,r₂)和G(-τ,r₁,r₂)分别为r₁,r₂两点间噪声场时空互相关函数、因果格林函数与反因果格林函数，τ为时延；

通过对互谱密度函数C(ω,r₁,r₂)做反傅里叶变换得到噪声场时空互相关函数C(τ,r₁,r₂)：

其中

为反傅里叶变换，C(ω,r₁,r₂)为r₁,r₂两点间噪声场互谱密度函数；

通过对两个水听器接收的声压信号求数学期望得到互谱密度函数：

C(ω,r₁,r₂)＝E[P(ω,r₁)P^*(ω,r₂)]

其中P(ω,r₁)和P^*(ω,r₂)分别为第一个水听器在角频率ω处的复声压以及第二个水听器在角频率ω处的复声压共轭，采用Warping变换来进行简正波分离，Warping变换之后的噪声场时空互相关函数W(τ,r₁,r₂)为：

其中h(t)为warping算子：

t_R为信号起始时间，t_R＝R/c，R和c分别为两水听器间距以及水中声速；

步骤2、利用自适应单纯形法模拟退火算法开展海底参数的反演，给出各个反演参数的边缘概率分布：

随机变量d和m分别表示观测数据向量和模型参数向量，N和M则分别表示向量d和m的变量个数；根据贝叶斯定理，d和m的关系表示为：

其中，P(m|d)为后验概率密度；P(d|m)为条件概率。P(m)和P(d)为模型参数向量和观测数据向量的先验概率；

构建代价函数：

其中d_k为第k阶模态的实测频散曲线，d_k(m)为参数向量下第k阶模态的理论频散曲线，N_k为第k阶模态频散曲线上取样的数量；

结合先验概率的代价函数φ(m)为：

φ(m)＝E(m)-log_eP(m)

最大后验概率分布

和一维边缘概率密度分布(MPD)P(m_i|d)用下式表征：

P(m_i|d)＝∫δ(m_i-m′_i)P(m′|d)dm′

其中

就是使()内达到最小值时的变量的取值。

所述互谱密度函数：对两个水听器接收的噪声信号做傅里叶变换得到各个频率点的声压值，第一个水听器处的声压值与第二个水听器处的声压值的复共轭二者相乘即得到两点间互谱密度：

其中p_j(ω,r₁)和p_j(ω,r₂)分别为第一个和第二个水听器接收的第j个快拍信号的频域声压解，J为快拍的总数，p_j(ω,r₁)可以通过对第j个快拍信号x_j(t)做傅里叶变换得到。

有益效果

本发明提出的一种基于非线性贝叶斯理论的简正波频散曲线被动反演方法，从统计角度定性和定量的分析反演结果的不确定性，为反演结果的可信度提供参考。首先利用噪声干涉仪技术与简正波分离技术，被动提取各阶水声信号的频散特征，然后利用自适应单纯形法模拟退火算法开展海底参数的反演，给出各个反演参数的边缘概率分布，为反演结果的准确性提供评判参考。

附图说明

图1：基于非线性贝叶斯理论的水声信号频散特征被动反演实施方案图

图2：实验站点与实验期间海水声速

图3：(a)噪声场时空互相关函数及其(b)时频图

图4：(a)warping变换之后的噪声场时空互相关函数及其(b)时频图

(c)单阶噪声场互相关信号波形(d)实测数据提取的各阶模态频散曲线

图5：各个参数边缘概率分布

具体实施方式

现结合实施例、附图对本发明作进一步描述：

图1给出的是基于非线性贝叶斯理论的水声信号频散特征被动反演流程，具体实施如下：先对空间两点噪声信号做互相关得到两点间噪声场互谱密度函数，进而得到两点间时空互相关函数；再对时空互相关函数做warping变换分离不同阶模态信号，提取各阶模态的频散曲线；最后利用水声信号频散特征构建贝叶斯反演框架，对地声参数进行反演，最终给出反演的各参数的一维边缘概率分布，评价反演结果是否可信。

发明通过以下几个步骤完成：

(a)利用噪声干涉仪与简正波分离技术提取水声信号频散特征

根据噪声干涉仪技术可知，空间两点噪声场互相关函数的导数可以用两点之间的因果格林函数与反因果格林函数之和表征：

其中C(τ,r₁,r₂)，G(τ,r₁,r₂)和G(-τ,r₁,r₂)分别为r₁,r₂两点间噪声场时空互相关函数、因果格林函数与反因果格林函数。事实上，通常空间两点的互相关函数比其导数能够更好的表征声场格林函数，所以通常更多采用两点噪声场互相关函数，而非其导数来近似声场格林函数。

噪声场时空互相关函数C(τ,r₁,r₂)可以通过对互谱密度函数C(ω,r₁,r₂)做反傅里叶变换得到：

其中

为反傅里叶变换，C(ω,r₁,r₂)为r₁,r₂两点间噪声场互谱密度函数，互谱密度函数可以通过对两个水听器接收的声压信号求数学期望得到：

C(ω,r₁,r₂)＝E[P(ω,r₁)P^*(ω,r₂)] (3)

P(ω,r₁)和P^*(ω,r₂)分别为第一个水听器在角频率ω处的复声压以及第二个水听器在角频率ω处的复声压共轭。实际数据处理中，互谱密度函数可以通过如下步骤获得：对两个水听器接收的噪声信号做傅里叶变换得到各个频率点的声压值，第一个水听器处的声压值与第二个水听器处的声压值的复共轭二者相乘即可得到两点间互谱密度：

其中p_j(ω,r₁)和p_j(ω,r₂)分别为两个水听器接收的第j个快拍信号的在角频率ω处的复声压，p_j(ω,r₁)可以通过对第j个快拍信号做短时傅里叶变换得到，J为快拍的总数。

由于海洋波导中接收的信号是典型的多分量信号，如果传播距离不够远的话不同阶模态在时频域会混叠，无法直接对其分离。所以需要借助信号处理手段将不同阶模态分离，这里采用Warping变换来进行简正波分离。

Warping变换之后的噪声场时空互相关函数W(τ,r₁,r₂)为：

其中warping算子h(t)为：

t_R为信号起始时间，t_R＝R/c，R和c分别为两水听器间距以及水中声速。Warping变换后的各阶简正波变成了以该阶简正波的截止频率为信号频率的单频信号，信号转换成了若干先后到达的单频信号的组合。对W(τ,r₁,r₂)做时频分析，利用窄带滤波即可得到噪声场时空互相关信号各单模分量W_m(τ,r₁,r₂)。Warping是一种可逆变换，对W_m(τ,r₁,r₂)做逆变换即可恢复噪声时空互相关函数单阶信号C_m(τ,r₁,r₂)。

其中warping逆变换的算子为：

对C_m(τ,r₁,r₂)做时频分析处理，即可得到第m阶简正波的频散特征。

(b)基于贝叶斯理论的水声信号频散特征被动反演

反演结果的不确定性是评价反演结果可信度的重要指标。假设随机变量d和m分别表示观测数据向量和模型参数向量，N和M则分别表示向量d和m的变量个数。根据贝叶斯定理，d和m的关系可表示为：

其中，P(m\d)为后验概率密度；P(d\m)为条件概率。P(m)和P(d)为模型参数向量和观测数据向量的先验概率。最大似然函数L(d|m)与代价函数E(m)之间的关系为L(d|m)∝exp(-E(m))，本发明采用如下形式的代价函数：

其中d_k为第k阶模态的实测频散特征，d_k(m)为参数向量下第k阶模态的理论频散特征。N_k为第k阶模态频散特征上取样的数量。更一般的结合先验概率的代价函数φ(m)定义为：

φ(m)＝E(m)-log_eP(m) (11)

最大后验概率分布

和一维边缘概率密度分布(MPD)P(m_i|d)可以用下式表征：

P(m_i|d)＝∫δ(m_i-m′_i)P(m′|d)dm′ (13)

本发明利用实测的海洋环境噪声数据验证算法的性能。实验站点以及实验期间的海水声速剖面见图2，水听器布放站点处海深为100m，海底地形变换缓慢，可近似认为海底是水平不变的。两个接收水听器都布放在95m深度处，两个水听器水平间距为5km。需要说明的是：考虑到水听器在布放时下沉的过程中位置会有所改变，所以利用GPS数据估计的水听器间距可能存在一定的误差，所以这里将水听器间距也作为未知量来反演。待反演的参数集为Ω＝[R,H,C_s,ρ_s,C_b,ρ_b,Δt]，其中各参数的含义分别为两水听器间水平距离、沉积层厚度、沉积层纵波声速、沉积层密度、基底纵波声速和基底密度，以及相对时间漂移量。

图3(a)给出的是实测数据提取的噪声场互相关函数正半轴部分，图3(b)给出的是噪声场互相关函数的时频图，从图中可以看到有四阶模态存在，但是彼此在时频域上有重叠，很难直接提取到各阶模态的频散特征。图4(a)和(b)分别给出的是warped变换之后的噪声场互相关函数及其时频图。从图中图4(b)可以看到warped变换之后的各阶噪声互相关信号在时频域上彼此不重叠，通过窄带滤波就能很容易地得到单阶模态信号。图4(c)给出的是1-4阶模态时域波形。对单阶模态信号做时频分析就能很容易获得各阶模态的频散曲线。为了获得最高分辨率的时频图，本发明采用Wigner-Ville Distribution来对各阶模态时域信号做时频分析，图4(d)中四条白色实线即为提取的4阶简正波的频散曲线。

表1各参数搜索区间以及反演得到的最优值

表1给出的是各个反演参数的搜索区间以及反演得到的最优值，并与ZangXiaoQin的博士论文《Acoustic Green's Function Extraction in the Ocean》给出的反演结果进行了对比，可以看到二者吻合的非常好，证明了反演结果的准确性。图5给出的是被反演的各个参数的一维边缘概率分布。其中竖线对应的是反演得到的最优值，曲线对应的被反演参数在反演区间范围内的概率分布。从图中可以看到：时间漂移量、水听器间距、沉积层厚度以及沉积层声速对简正波频散特征很敏感，反演得到的最优值基本都出现在边缘概率分布的最大值附近，证明上述几个参数反演的结果可信度很高。相对于上述几个参数，沉积层密度和基底密度对简正波频散特性敏感度相对较弱，但是这两个参数反演过程中出现在最优值的概率仍然相对较高，证明这两个参数的反演结果具有一定的可信度。基底声速对反演敏感度较差，相对于其他参数，反演得到的该参数可信度最低。

Claims (2)

1.一种基于非线性贝叶斯理论的水声信号频散特征被动反演方法，其特征在于步骤如下：

步骤1、利用噪声干涉仪技术与简正波分离技术，被动提取各阶水声信号的频散特征：

根据噪声干涉仪技术，空间两点噪声场互相关函数的导数用两点之间的因果格林函数与反因果格林函数之和表征：

其中C(τ,r₁,r₂)，G(τ,r₁,r₂)和G(-τ,r₁,r₂)分别为r₁,r₂两点间噪声场时空互相关函数、因果格林函数与反因果格林函数，τ为时延；

通过对互谱密度函数C(ω,r₁,r₂)做反傅里叶变换得到噪声场时空互相关函数C(τ,r₁,r₂)：

其中

为反傅里叶变换，C(ω,r₁,r₂)为r₁,r₂两点间噪声场互谱密度函数；

通过对两个水听器接收的声压信号求数学期望得到互谱密度函数：

C(ω,r₁,r₂)＝E[P(ω,r₁)P^*(ω,r₂)]

其中P(ω,r₁)和P^*(ω,r₂)分别为第一个水听器在角频率ω处的复声压以及第二个水听器在角频率ω处的复声压共轭，采用Warping变换来进行简正波分离，Warping变换之后的噪声场时空互相关函数W(τ,r₁,r₂)为：

其中h(t)为warping算子：

t_R为信号起始时间，t_R＝R/c，R和c分别为两水听器间距以及水中声速；

步骤2、利用自适应单纯形法模拟退火算法开展海底参数的反演，给出各个反演参数的边缘概率分布：

随机变量d和m分别表示观测数据向量和模型参数向量，N和M则分别表示向量d和m的变量个数；根据贝叶斯定理，d和m的关系表示为：

其中，P(m|d)为后验概率密度；P(d|m)为条件概率。P(m)和P(d)为模型参数向量和观测数据向量的先验概率；

构建代价函数：

其中d_k为第k阶模态的实测频散曲线，d_k(m)为参数向量下第k阶模态的理论频散曲线，N_k为第k阶模态频散曲线上取样的数量；

结合先验概率的代价函数φ(m)为：

φ(m)＝E(m)-log_e P(m)

最大后验概率分布

和一维边缘概率密度分布(MPD)P(m_i|d)用下式表征：

P(m_i|d)＝∫δ(m_i-m′_i)P(m′|d)dm′

其中

就是使()内达到最小值时的变量的取值。

2.根据权利要求1所述基于非线性贝叶斯理论的水声信号频散特征被动反演方法，其特征在于：所述互谱密度函数：对两个水听器接收的噪声信号做傅里叶变换得到各个频率点的声压值，第一个水听器处的声压值与第二个水听器处的声压值的复共轭二者相乘即得到两点间互谱密度：

其中p_j(ω,r₁)和p_j(ω,r₂)分别为第一个和第二个水听器接收的第j个快拍信号的频域声压解，J为快拍的总数，p_j(ω,r₁)可以通过对第j个快拍信号x_j(t)做傅里叶变换得到。

Images