搜索

x

留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

联合简正波水波和底波频散特性的贝叶斯地声参数反演

郝望 段睿 杨坤德

引用本文:
Citation:

联合简正波水波和底波频散特性的贝叶斯地声参数反演

郝望, 段睿, 杨坤德

Bayesian geoacoustic parameter inversion based on dispersion characteristics of normal mode water wave and ground wave

Hao Wang, Duan Rui, Yang Kun-De
PDF
HTML
导出引用
  • 大多数基于浅海简正波模态频散数据的地声参数反演方法无法对深层底质声学参数进行可靠估计, 究其原因是仅利用了简正波水波频散特征, 忽略了与深层底质声学参数密切相关的底波频散特征, 因此, 本文在分析了包含水波和底波的浅海宽带数据的基础上, 结合底波频散特征对深层底质声学参数变化更加敏感的物理特性, 实现了基于完整简正波频散特性的贝叶斯地声参数反演, 并针对简正波宽带声场模型计算复杂度较高的现实问题, 利用变分贝叶斯蒙特卡罗方法的推断优势, 完成了未知参数的可靠估计和快速后验分析. 仿真和海上实验结果表明: 联合简正波水波和底波频散特征数据的贝叶斯地声参数反演, 不仅可以有效估计深层底质声学参数, 而且降低了其他相关环境参数的估计不确定性.
    Most of shallow water geoacoustic inversions based on modal dispersion cannot reliably estimate the deep geoacoustic parameters. Because these studies focused on the dispersions of water waves but ignored the dispersions of ground waves. Therefore, in this paper a Bayesian geoacoustic inversion is studied based on wideband modal dispersions of water waves and ground waves. Firstly, the modal dispersion curves with Airy phase components are discussed. Secondly, the Bayesian inversion theory and a novel sample-efficient inference algorithm, namely Variational Bayesian Monte Carlo, are introduced briefly. In the Bayesian inversion, the posterior probability densities of unknown parameters are inferred, which can provide the prediction closest to the observation data and the uncertainty of the prediction. Considering that the forward acoustic model is computationally intensive, the posterior analysis is carried out by using the Variational Bayesian Monte Carlo method. It is performed by finding the variational distribution closest to the target distribution and requires less computation time than the Markov chain Monto Carlo method. In the simulation study, a range-independent two-layer seabed, including the sediment layer and basement layer, is modeled, on the assumption that the water column is homogeneous. The function of shear wave in waveguide is ignored. The compressional sound speed of the sediment layer varies linearly from c1U to c1L between 0 and h1, while other geoacoustic parameters are constant. By comparing the inversion results with and without the information of ground waves for different signal-to-noise ratios, it can be concluded that the deep geoacoustic parameters are more sensitive to the dispersions of ground waves. And then, a shallow-water experimental study is carried out in the Bohai Sea of China. The average water depth is about 20 m. The wideband pulse signals are recorded by a hydrophone with a sensitivity of –170 dB re 1 V/μPa. The received signals include well-defined Airy phase components, and the modal dispersion curves of water waves and ground waves are extracted accurately. The experimental results indicate that the Bayesian inversion combining water and ground wave dispersions can not only estimate the deep geoacoustic parameters reliably, but also reduce the inversion uncertainties of other model parameters, such as the shallow geoacoustic parameters, water depth, etc. The estimated source-receiver range and water sound speed are close to their measured values. The modal dispersion curves predicted by the posterior mean samples are in good consistence with those extracted from the experimental data in different ranges. In addition, the good forecast of transmission loss also demonstrates the reliability of the joint inversion.
      通信作者: 段睿, duanrui@nwpu.edu.cn
    • 基金项目: 国家自然科学基金(批准号: 12074315, 42076198)和中船726研究所重点实验室基金(批准号: JCKY2021207CH02)资助的课题
      Corresponding author: Duan Rui, duanrui@nwpu.edu.cn
    • Funds: Project supported by the National Natural Science Foundation of China (Grant Nos. 12074315, 42076198) and the Open Fund for Key Laboratory of Underwater Acoustic Countermeasure Technology, China (Grant No. JCKY2021207CH02).
    [1]

    尚尔昌 2019 应用声学 38 468Google Scholar

    Shang E C 2019 Appl. Acoust. 38 468Google Scholar

    [2]

    Chapman N, Shang E C 2021 J. Theor. Comp. Acout. 29 2130004Google Scholar

    [3]

    Dosso S E, Dettmer J 2011 Inverse Probl. 27 055009Google Scholar

    [4]

    Shen Y, Pan X, Zheng Z, Gerstoft P 2020 J. Acoust. Soc. Am. 148 3279Google Scholar

    [5]

    杨坤德, 马远良 2009 物理学报 58 1798Google Scholar

    Yang K D, Ma Y L 2009 Acta Phys. Sin. 58 1798Google Scholar

    [6]

    Chapman N 2016 J. Mar. Sci. Eng. 4 61Google Scholar

    [7]

    Bonnel J, Chapman N 2011 J. Acoust. Soc. Am. 130 EL101Google Scholar

    [8]

    郭晓乐, 杨坤德, 马远良 2015 物理学报 64 174302Google Scholar

    Guo X L, Yang K D, Ma Y L 2015 Acta Phys. Sin. 64 174302Google Scholar

    [9]

    Duan R, Chapman N, Yang K D, Ma Y L 2016 J. Acoust. Soc. Am. 139 70Google Scholar

    [10]

    李佳蔚, 鹿力成, 郭圣明, 马力 2017 物理学报 66 204301Google Scholar

    Li J W, Lu L C, Guo S M, Ma L 2017 Acta Phys. Sin. 66 204301Google Scholar

    [11]

    Lin Y T, Lynch J F, Chotiros N, Chen C F, Newhall A, Turgut A, Schock S G, Chiu C S, Bartek L, Liu C S 2004 IEEE J. Oceanic Engineer. 29 1231Google Scholar

    [12]

    Wan L, Badiey M, Knobles D P, Wilson P S 2018 J. Acoust. Soc. Am. 143 EL199Google Scholar

    [13]

    Blei D M, Kucukelbir A, McAuliffe J D 2017 J. Am. Stat. Assoc. 112 859Google Scholar

    [14]

    Acerbi L 2018 arXiv: 1810.05558 v2 [stat. ML]

    [15]

    Che Y F, Wu X, Pastore G, Li W, Shirvan K 2021 Ann. Nucl. Energy 153 108046Google Scholar

    [16]

    Jensen F B, Kuperman W A, Porter M B, Schmidt H 2011 Computational Ocean Acoustics (Vol. 2) (New York: Springer) pp337–452

    [17]

    Dosso S E, Nielsen P L, Wilmut M J 2006 J. Acoust. Soc. Am. 119 208Google Scholar

    [18]

    Porter M B 1991 The KRAKEN Normal Mode Program (La Spezia: SACLANT Undersea Research Center) Technical Report SM-245

    [19]

    Bonnel J, Thode A, Wright D, Chapman N 2020 J. Acoust. Soc. Am. 147 1897Google Scholar

    [20]

    Daubechies I, Lu J F, Wu H T 2011 Appl. Comput. Harmon. A. 30 243Google Scholar

    [21]

    Dosso S E, Wilmut M J, Lapinski A L S 2001 IEEE J. Oceanic Engineer. 26 324Google Scholar

    [22]

    王鹏, 贾凯, 吴建政, 胡日军 2015 海洋地质与第四纪地质 35 23

    Wang P, Jia K, Wu J Z, Hu R J 2015 Marine Geology Quaternary Geology 35 23

    [23]

    张剑, 李日辉, 王中波, 张训华, 黄龙, 孙荣涛 2016 海洋地质与第四纪地质 36 1

    Zhang J, Li R H, Wang Z B, Zhang X H, Huang L, Sun R T 2016 Marine Geology & Quaternary Geology 36 1

    [24]

    Li Z L, Zhang R H 2004 Chin. Phys. Lett. 21 1100Google Scholar

  • 图 1  地声参数反演流程图

    Fig. 1.  Flow chart of geoacoustic inversion.

    图 2  3层水平分层介质波导模型

    Fig. 2.  Three-layer waveguide model.

    图 3  1阶模态群速度曲线 (a)不同的$ {h_1} $; (b)不同的${c_{{\text{1 L}}}}$; (c)不同的${c_2}$; (d)不同的${\rho _2}$, 圆圈给出了艾里相频率的位置

    Fig. 3.  Group velocity curves of Mode 1: (a) Different values of $ {h_1} $; (b) different values of ${c_{{\text{1 L}}}}$; (c) different values of ${c_2}$; (d) different values of ${\rho _2}$. The circles indicate the Airy phase frequencies.

    图 4  仿真信号 (a)归一化时域波形; (b)时频图, 其中1阶模态的底波、水波和Airy相结构清晰可辩

    Fig. 4.  Simulation signal: (a) Normalized time domain waveform; (b) time-frequency diagram. The ground wave, water wave, and Airy phase component of Mode 1 are well-defined.

    图 5  不同信噪比条件下底波频散曲线的单次提取结果(a)和提取误差标准差(b)

    Fig. 5.  Single extraction results (a) and standard deviations of extraction errors (b) of ground wave dispersion curve under different SNRs.

    图 6  两种数据条件下部分参数的一维边缘后验概率密度函数 (a)$ {h_1} $; (b)${c_{{\text{1 L}}}}$; (c)${c_2}$; (d)${\rho _2}$

    Fig. 6.  1 D marginal posterior probability densities of some parameters for Data 1 and Data 2: (a) $ {h_1} $; (b) ${c_{{\text{1 L}}}}$; (c) ${c_2}$; (d) ${\rho _2}$.

    图 7  实验描述 (a)实验地点; (b)水体声速剖面

    Fig. 7.  Description of the experiment: (a) Experimental site; (b) sound speed profile in water.

    图 8  实验信号 (a)归一化时域波形; (b)时频图, 其中1阶模态的底波、水波和Airy相结构清晰可辩

    Fig. 8.  Experimental signal: (a) Normalized time domain waveform; (b) time-frequency diagram. The ground wave, water wave, and Airy phase component of Mode 1 are well-defined.

    图 9  频散曲线的实测结果和后验预测结果 (a)数据条件1, 水波和底波频散数据, 插图为低频部分的局部放大; (b)数据条件2, 水波频散数据

    Fig. 9.  Measurements and posterior predictions of dispersion curves: (a) Data 1, water waves and ground waves, the enlarged view of the low frequency part is inset; (b) Data 2, water waves.

    图 10  待反演参数的一维(数据条件1和条件2)和二维(数据条件1)边缘后验概率密度函数

    Fig. 10.  1D (Data 1 and Data 2) and 2D (Data 1) marginal posterior probability densities of unknown parameters.

    图 11  针对数据条件1的归一化后验参数协方差矩阵

    Fig. 11.  Normalized posterior covariance matrix of unknown parameters for Data 1.

    图 12  针对数据条件1的VBMC收敛情况 (a)证据下界ELBO; (b) KL散度

    Fig. 12.  Convergences of the VBMC method for Data 1: (a) ELBO; (b) KL divergence.

    图 13  针对数据条件1的压缩波声速和密度剖面后验估计结果 (a) 2层海底模型; (b) 3层海底模型

    Fig. 13.  Posteriori estimates of compressed-wave sound speed profiles and density profiles for Data 1: (a) Two-layer seabed model; (b) three-layer seabed model.

    图 14  另一组宽带脉冲信号时频图和频散曲线预测结果(实线) (a)接收距离3.42 km; (b)接收距离5.35 km; (c)接收距离6.51 km

    Fig. 14.  Time-frequency diagrams and dispersion curve predictions (solid lines) of other broadband pulse signals: (a) The range is 3.42 km; (b) the range is 5.35 km; (c) the range is 6.51 km.

    图 15  传播损失的理论计算结果和实测结果比较, 声源深度6 m, 接收深度11.5 m (a)频率200 Hz; (b)频率315 Hz

    Fig. 15.  Comparisons of the theoretical and experimental transmission loss at source depth 6 m and receiver depth 11.5 m: (a) The frequency is 200 Hz; (b) the frequency is 315 Hz.

    表 1  VBMC算法步骤

    Table 1.  Steps of the VBMC algorithm.

    算法 : 变分贝叶斯蒙特卡罗(Acerbi, 2018)
    1: $t \leftarrow 0$, 初始化GMM模型参数${\phi _0}$和${K_0}$
    2: When 收敛条件未满足, do
    3: $t \leftarrow t + 1$ //迭代
    4:  for $ 1 \ldots {n_{\text{active}} }$ do
    5:  采样新的样本点$\boldsymbol{\theta} _{\rm new} \leftarrow {\rm argmax}_{\boldsymbol{\theta}} [a({\boldsymbol\theta})]$ //$a({\boldsymbol{\theta }})$
       为采样函数
    6:  计算样本点的$ \log p(\boldsymbol{d}, {\boldsymbol{\theta} _{\text{new}} }) $值, 并将其加入训练集
    7: $\log p({\boldsymbol{d}}, {\boldsymbol{\theta}} ) \leftarrow$训练GP模型, 得到最优的模型超参数
    8: ${K_t} \leftarrow $更新GMM模型参数
    9: ${\phi _t} \leftarrow $利用随机梯度下降法对ELBO寻优, 更新GMM
      模型参数
    10: end
    下载: 导出CSV

    表 2  反演参数空间和先验区间

    Table 2.  Parameter spaces and prior bounds for inversion.

    待反演参数单位仿真参数值先验区间
    沉积层厚度$ {h_1} $m47[10, 80]
    沉积层密度${\rho _1}$${\text{g/c}}{{\text{m}}^{\text{3}}}$1.5[1, 2]
    基底层密度${\rho _2}$${\text{g/c}}{{\text{m}}^{\text{3}}}$2.1[1.5, 2.5]
    沉积层上部声速$ {c_{1{\text{U}}}} $${\text{m/s}}$1675[1600, 1700]
    沉积层下部声速$ {c_{1{\text{L}}}} $${\text{m/s}}$1675[1650, 1750]
    基底层声速$ {c_2} $${\text{m/s}}$1765[1700, 1800]
    水深dwm21[15, 25]
    收发距离$r$km6.24[4, 8]
    时间因子$ {\delta _{\text{t} }} $s[–0.5, 0.5]
    水体声速cw${\text{m/s}}$1511[1500, 1520]
    下载: 导出CSV

    表 3  针对数据条件1的后验统计结果

    Table 3.  Posterior statistical results for Data 1.

    待反演参数均值方差95%可信区间实测值
    $ {h_1} $46.653.77[39.69, 54.21]
    ${\rho _1}$1.520.09[1.35, 1.69]
    ${\rho _2}$2.120.19[1.75, 2.44]
    $ {c_{1{\text{U}}}} $1675.899.40[1654.35, 1690.32]
    $ {c_{1{\text{L}}}} $1673.528.04[1660.41, 1691.35]
    $ {c_2} $1765.0315.78[1731.16, 1789.77]
    $ {d_{\text{w}}} $20.690.45[19.70, 21.41]19—22
    $r$6.240.30[5.63, 6.83]6.07
    $ {\delta _{\text{t}}} $0.2800.002[0.276, 0.284]
    $ {c_{\text{w}}} $1510.984.13[1503.08, 1517.91]1511—1513
    下载: 导出CSV
  • [1]

    尚尔昌 2019 应用声学 38 468Google Scholar

    Shang E C 2019 Appl. Acoust. 38 468Google Scholar

    [2]

    Chapman N, Shang E C 2021 J. Theor. Comp. Acout. 29 2130004Google Scholar

    [3]

    Dosso S E, Dettmer J 2011 Inverse Probl. 27 055009Google Scholar

    [4]

    Shen Y, Pan X, Zheng Z, Gerstoft P 2020 J. Acoust. Soc. Am. 148 3279Google Scholar

    [5]

    杨坤德, 马远良 2009 物理学报 58 1798Google Scholar

    Yang K D, Ma Y L 2009 Acta Phys. Sin. 58 1798Google Scholar

    [6]

    Chapman N 2016 J. Mar. Sci. Eng. 4 61Google Scholar

    [7]

    Bonnel J, Chapman N 2011 J. Acoust. Soc. Am. 130 EL101Google Scholar

    [8]

    郭晓乐, 杨坤德, 马远良 2015 物理学报 64 174302Google Scholar

    Guo X L, Yang K D, Ma Y L 2015 Acta Phys. Sin. 64 174302Google Scholar

    [9]

    Duan R, Chapman N, Yang K D, Ma Y L 2016 J. Acoust. Soc. Am. 139 70Google Scholar

    [10]

    李佳蔚, 鹿力成, 郭圣明, 马力 2017 物理学报 66 204301Google Scholar

    Li J W, Lu L C, Guo S M, Ma L 2017 Acta Phys. Sin. 66 204301Google Scholar

    [11]

    Lin Y T, Lynch J F, Chotiros N, Chen C F, Newhall A, Turgut A, Schock S G, Chiu C S, Bartek L, Liu C S 2004 IEEE J. Oceanic Engineer. 29 1231Google Scholar

    [12]

    Wan L, Badiey M, Knobles D P, Wilson P S 2018 J. Acoust. Soc. Am. 143 EL199Google Scholar

    [13]

    Blei D M, Kucukelbir A, McAuliffe J D 2017 J. Am. Stat. Assoc. 112 859Google Scholar

    [14]

    Acerbi L 2018 arXiv: 1810.05558 v2 [stat. ML]

    [15]

    Che Y F, Wu X, Pastore G, Li W, Shirvan K 2021 Ann. Nucl. Energy 153 108046Google Scholar

    [16]

    Jensen F B, Kuperman W A, Porter M B, Schmidt H 2011 Computational Ocean Acoustics (Vol. 2) (New York: Springer) pp337–452

    [17]

    Dosso S E, Nielsen P L, Wilmut M J 2006 J. Acoust. Soc. Am. 119 208Google Scholar

    [18]

    Porter M B 1991 The KRAKEN Normal Mode Program (La Spezia: SACLANT Undersea Research Center) Technical Report SM-245

    [19]

    Bonnel J, Thode A, Wright D, Chapman N 2020 J. Acoust. Soc. Am. 147 1897Google Scholar

    [20]

    Daubechies I, Lu J F, Wu H T 2011 Appl. Comput. Harmon. A. 30 243Google Scholar

    [21]

    Dosso S E, Wilmut M J, Lapinski A L S 2001 IEEE J. Oceanic Engineer. 26 324Google Scholar

    [22]

    王鹏, 贾凯, 吴建政, 胡日军 2015 海洋地质与第四纪地质 35 23

    Wang P, Jia K, Wu J Z, Hu R J 2015 Marine Geology Quaternary Geology 35 23

    [23]

    张剑, 李日辉, 王中波, 张训华, 黄龙, 孙荣涛 2016 海洋地质与第四纪地质 36 1

    Zhang J, Li R H, Wang Z B, Zhang X H, Huang L, Sun R T 2016 Marine Geology & Quaternary Geology 36 1

    [24]

    Li Z L, Zhang R H 2004 Chin. Phys. Lett. 21 1100Google Scholar

  • [1] 刘勇, 涂国华, 向星皓, 李晓虎, 郭启龙, 万兵兵. 横向矩形微槽抑制高超声速第二模态扰动波的参数化研究. 物理学报, 2022, 71(19): 194701. doi: 10.7498/aps.71.20220851
    [2] 高飞, 徐芳华, 李整林, 秦继兴. 大陆坡内波环境中声传播模态耦合及强度起伏特征. 物理学报, 2022, 71(20): 204301. doi: 10.7498/aps.71.20220634
    [3] 魏广宇, 陈凝飞, 仇志勇. 高能量粒子测地声模与Dimits区漂移波相互作用. 物理学报, 2022, 71(1): 015201. doi: 10.7498/aps.71.20211430
    [4] 孙冠文, 崔寒茵, 李超, 林伟军. 火星大气频散声速剖面建模方法及其对声传播路径的影响. 物理学报, 2022, 71(24): 244304. doi: 10.7498/aps.71.20221531
    [5] 魏广宇, 陈凝飞, 仇志勇. 高能量粒子测地声模与Dimits区漂移波相互作用. 物理学报, 2021, (): . doi: 10.7498/aps.70.20211430
    [6] 江鹏飞, 林建恒, 孙军平, 衣雪娟. 考虑噪声源深度分布的海洋环境噪声模型及地声参数反演. 物理学报, 2017, 66(1): 014306. doi: 10.7498/aps.66.014306
    [7] 郭晓乐, 杨坤德, 马远良, 杨秋龙. 一种基于简正波模态消频散变换的声源距离深度估计方法. 物理学报, 2016, 65(21): 214302. doi: 10.7498/aps.65.214302
    [8] 郭晓乐, 杨坤德, 马远良. 一种基于简正波模态频散的远距离宽带海底参数反演方法. 物理学报, 2015, 64(17): 174302. doi: 10.7498/aps.64.174302
    [9] 屠惠琳, 肖绍球, 杨智杰, 王秉中. 基于时间反演电磁波的微结构天线的单频点超分辨力聚焦研究. 物理学报, 2014, 63(8): 084102. doi: 10.7498/aps.63.084102
    [10] 周天, 李海森, 朱建军, 魏玉阔. 利用多角度海底反向散射信号进行地声参数估计. 物理学报, 2014, 63(8): 084302. doi: 10.7498/aps.63.084302
    [11] 张强, 李宏宇, 张立阳, 岳平, 史晋森. 陇中黄土高原自然植被下垫面陆面过程及其参数对降水波动的气候响应. 物理学报, 2013, 62(1): 019201. doi: 10.7498/aps.62.019201
    [12] 毛杰健, 杨建荣. 大尺度浅水波方程中相互调制的非线性波. 物理学报, 2013, 62(13): 130205. doi: 10.7498/aps.62.130205
    [13] 张亮, 张立凤, 吴海燕, 王骥鹏. 黏性水波振荡型行波解的存在性. 物理学报, 2009, 58(2): 703-711. doi: 10.7498/aps.58.703
    [14] 杨坤德, 马远良. 利用海底反射信号进行地声参数反演的方法. 物理学报, 2009, 58(3): 1798-1805. doi: 10.7498/aps.58.1798
    [15] 钟兰花, 吴福根. 水波在周期性钻孔底部结构中的传播及其能带. 物理学报, 2009, 58(9): 6363-6368. doi: 10.7498/aps.58.6363
    [16] 李富才, 孟 光. 窄频带Lamb波频散特性研究. 物理学报, 2008, 57(7): 4265-4272. doi: 10.7498/aps.57.4265
    [17] 莫嘉琪, 林万涛. 一类大气浅水波方程的近似解. 物理学报, 2007, 56(7): 3662-3666. doi: 10.7498/aps.56.3662
    [18] 肖 夏, 尤学一, 姚素英. 表征超大规模集成电路互连纳米薄膜硬度特性的声表面波的频散特性. 物理学报, 2007, 56(4): 2428-2433. doi: 10.7498/aps.56.2428
    [19] 孙 江, 姜 谦, 米 辛, 俞祖和, 傅盘铭. 利用场关联效应抑制瑞利型非简并四波混频的热背底. 物理学报, 2004, 53(2): 450-455. doi: 10.7498/aps.53.450
    [20] 张解放. 长水波近似方程的多孤子解. 物理学报, 1998, 47(9): 1416-1420. doi: 10.7498/aps.47.1416
计量
  • 文章访问数:  711
  • PDF下载量:  19
  • 被引次数: 0
出版历程
  • 收稿日期:  2022-08-31
  • 修回日期:  2023-01-02
  • 上网日期:  2023-01-07
  • 刊出日期:  2023-03-05

联合简正波水波和底波频散特性的贝叶斯地声参数反演

  • 1. 西北工业大学航海学院, 西安 710072
  • 2. 西北工业大学海洋研究院, 太仓 215400
  • 3. 西北工业大学, 海洋声学信息感知工业和信息化部重点实验室, 西安 710072
  • 通信作者: 段睿, duanrui@nwpu.edu.cn
    基金项目: 国家自然科学基金(批准号: 12074315, 42076198)和中船726研究所重点实验室基金(批准号: JCKY2021207CH02)资助的课题

摘要: 大多数基于浅海简正波模态频散数据的地声参数反演方法无法对深层底质声学参数进行可靠估计, 究其原因是仅利用了简正波水波频散特征, 忽略了与深层底质声学参数密切相关的底波频散特征, 因此, 本文在分析了包含水波和底波的浅海宽带数据的基础上, 结合底波频散特征对深层底质声学参数变化更加敏感的物理特性, 实现了基于完整简正波频散特性的贝叶斯地声参数反演, 并针对简正波宽带声场模型计算复杂度较高的现实问题, 利用变分贝叶斯蒙特卡罗方法的推断优势, 完成了未知参数的可靠估计和快速后验分析. 仿真和海上实验结果表明: 联合简正波水波和底波频散特征数据的贝叶斯地声参数反演, 不仅可以有效估计深层底质声学参数, 而且降低了其他相关环境参数的估计不确定性.

English Abstract

    • 浅海波导环境复杂多样, 其声学特性参数主要包括水体和海底参数, 例如声速和密度等, 这些参数对于声场预报和目标定位等研究有着重要意义. 水体声学参数可以通过直接测量的方式快速准确地获得, 但是对于海底声学参数, 也称地声参数, 直接测量起来较为困难, 需要耗费较大的人力和物力, 且通常只能获得实验站位处的浅层测量结果, 因此, 利用水声观测数据进行地声参数反演的研究就显得尤为重要[1], 地声参数反演具有成本低、效率高和范围广等众多优势, 现有的地声参数反演方法有传统的匹配反演和线性扰动反演, 以及近些年研究的热点, 如贝叶斯反演和神经网络等[2].

      在浅海地声参数反演研究中, 经常使用的水声观测数据有声压场数据[3,4]、海底反射数据[5,6]和模态数据[7-10]等, 其中, 利用模态频散数据进行地声参数反演的研究已经发展较为成熟, 该方法的主要特点是仅利用单个接收水听器就可以实现相关环境参数的可靠估计, 无需已知声源级, 例如: 文献[7]提出利用Warping变换对一定距离上的宽带脉冲接收信号进行模态分离, 基于分离的水波模态频散特征数据进行地声参数反演; 另外, 针对远距离海底参数反演性能下降的问题, 文献[8]提出采用双声源模态相对到达时间差数据进行地声参数反演, 消除了距离累积带来的误差, 实现了远距离地声参数的可靠估计. 虽然, 上述研究已经充分证明利用浅海模态频散数据可以有效地估计相关地声参数, 但是不难发现, 几乎所有的研究都是基于模态的水波频散特征, 而对于底波频散特征的研究极少, 所能利用的实验数据也较少, 导致反演结果仅能给出海底浅层地声参数的有效估计, 对于深层地声参数基本无效, 因此, 本文联合简正波的水波和底波频散特性进行地声参数反演, 一方面实现对深层地声参数的可靠估计, 另一方面降低其他模型输入参数的反演不确定性.

      值得注意的是, 2004年和2018年, Lin等[11]和Wan等[12]分别提出利用包含Airy相结构的模态频散特征数据对海底的压缩波声速进行反演, 但是对于其他类型的环境参数文中没有进行过多讨论, 除此之外, 前者没有进行参数反演的不确定性分析, 后者采用效率较低的基于采样策略的推断方法进行了简单的不确定性分析. 随着非线性贝叶斯近似推断理论的发展, 地声参数反演不仅可以给出模型输入参数的后验估计结果, 而且可以定量地描述参数估计的不确定性以及参数之间的相关性, 对于大多数地声参数反演问题, 前向模型的巨大计算代价导致基于采样策略的近似推断方法无法在短时间内完成参数的后验分析, 比如遗传算法(genetic algorithm, GA)和马尔科夫链蒙特卡罗(Markov chain Monto Carlo, MCMC)采样方法等, 相比之下, 变分推断(variational inference, VI)[13]作为另一类近似推断方法, 因其计算效率更高, 且计算结果较为可靠, 将有望代替MCMC采样, 成为一种更加实用的后验分析工具. 本文使用一种数据高效的近似推断方法, 也称变分贝叶斯蒙特卡罗(variational Bayesian Monte Carlo, VBMC)[14,15], 完成地声参数反演的后验分析工作.

      本文第2节为基本理论, 分别介绍了包含Airy相结构的频散曲线、贝叶斯地声参数反演理论、VBMC方法以及地声参数反演流程; 第3节为仿真数据验证; 第4节为实验数据验证; 第5节为结论.

    • 根据简正波声场理论, 当声源发射深度为$ {z_{\text{s}}} $, 水听器接收深度为$ {z_{\text{r}}} $, 收发距离为$r$时, 接收复声压谱$ Y(f) $可以表示为[16]

      $ Y(f)\approx AS(f){\displaystyle \sum _{m=1}^{M}{\psi }_{m}(f,{z}_{\text{s}}){\psi }_{m}(f,{z}_{\text{r}})}\times \frac{{\text{e}}^{\text{j}{k}_{rm}\left(f\right)r}}{\sqrt{{k}_{rm}(f)r}}\text{, } $

      其中, $ f $为信号频率, $ A $为常数因子, $ S(f) $为声源的发射谱, $M$为总的传播模态数, $ {\psi _m} $$ {k_{rm}} $分别为第$ m $阶简正波模态的深度函数和水平波数.

      由(1)式可得, 复声压谱$Y$是由各阶简正波叠加组成, 对于第$ m $阶简正波, 不同频率的信号能量具有不同的传播速度, 也称群速度$ v_{\text{g}}^{(m)} $, 其表达式为

      $ v_{\text{g}}^{\left( m \right)}(f) = \frac{{2{\text{π d}}f}}{{{\text{d}}{k_{rm}}(f)}}. $

      理论上, 随频率变化的群速度曲线被Airy相频率$f_{\text{a}}^{(m)}$分成两个部分, 分别称为简正波传播中的底波和水波部分, 特别地, 群速度取值在越靠近$f_{\text{a}}^{(m)}$时越小, 并在该点处取得极小值, 此时$ v_{\text{g}}^{(m)} $关于$ f $的一阶导函数等于0, 可以写作[11]

      $ {\left. {\frac{{{\text{d}}v_{\text{g}}^{\left( m \right)}(f)}}{{{\text{d}}f}}} \right|_{f = f_a^{(m)}}} = 0. $

      各阶简正波模态信号具有随频率变化的群速度, 因此, 当收发距离$r$足够远时, 宽带脉冲信号会表现出明显的模内频散和模间频散现象, 假设在$ {t_0} $时刻发射一个宽带脉冲信号, 第$m$阶模态信号到达接收点的绝对时间可以表示为

      $ {t^{(m)}}(f,r) = {t_0} + \frac{r}{{v_{\text{g}}^{\left( m \right)}(f)}}, $

      然而, 在实际情况中, 因为收发系统非严格同步, 所以通常获得第$m$阶模态信号的相对到达时间, 本文称为模态频散曲线, 可以表示为

      $ \begin{split} t_{{\text{re}}}^{(m)}(f,r) = \;&{t_0} + \frac{r}{{v_{\text{g}}^{(m)}(f)}} - t' \\ =\;& \frac{r}{{v_{\text{g}}^{(m)}(f)}} - \left( {t' - {t_0}} \right) \\ =\;& \frac{r}{{v_{\text{g}} ^{(m)}(f)}} - \left( {\frac{r}{{{c_{\text{w}}}}} - {\delta _{\text{t}}}} \right), \end{split} $

      式中, $ t' $为实验数据处理的参考时间, $ {c_{\text{w}}} $为海水的平均声速, $ {\delta _{\text{t}}} $为未知的常数时间因子, 该时间因子在反演时作为辅助变量可以快速收敛. 根据(5)式可知, 模态相对到达时间$ t_{{\text{re}}}^{(m)}(f, r) $和群速度$ v_{\text{g}}^{(m)}(f) $密切相关, 而群速度由环境参数(水体声学参数和地声参数)决定, 因此, 利用$ t_{{\text{re}}}^{(m)}(f, r) $的完整观测数据, 可以实现对相关环境参数的反演.

    • 在贝叶斯反演框架内, 未知的模型输入参数被描述为多维随机变量, 并从后验概率的角度对其进行估计和不确定性分析. 结合本文的具体问题, 采用随机变量${\boldsymbol{\theta}}$表示未知的环境参数, 包含部分水体参数和地声参数, ${\boldsymbol{d}} = [{{\boldsymbol{d}}_1}, {{\boldsymbol{d}}_2}, \cdots , {{\boldsymbol{d}}_M}]$表示频散曲线的观测向量集合, 其中, ${{\boldsymbol{d}}_m}(m = 1, 2, \cdots , M)$为第$m$阶模态的观测结果, 频点数为$ {N_m} $, 通常情况下, ${{\boldsymbol{d}}_m}$可以表示为正演模型预测项${g_m}({\boldsymbol{\theta}} )$和误差项${{\boldsymbol{\varepsilon}} _m}$之和, 即

      $ {{\boldsymbol{d}}_m} = {g_m}({\boldsymbol{\theta }}) + {{\boldsymbol{\varepsilon}} _m}. $

      $ {g_m}({\boldsymbol{\theta}} ) $由(5)式给出, ${{\boldsymbol{\varepsilon}} _m}$包含模型计算误差和实际测量误差.

      贝叶斯地声参数反演的核心是贝叶斯公式, 随机变量${\boldsymbol{\theta}}$的后验概率密度函数$p({\boldsymbol{\theta}} |{\boldsymbol{d}})$可以写作[3]:

      $ p({\boldsymbol{\theta}} |{\boldsymbol{d}}) = {{p({\boldsymbol{\theta}} )p({\boldsymbol{d}}|{\boldsymbol{\theta}} )} /{p({\boldsymbol{d}})}}, $

      $ p({\boldsymbol{d}}) = \int {p({\boldsymbol{\theta}} )p({\boldsymbol{d}}|{\boldsymbol{\theta}} )} {\text{d}}{\boldsymbol{\theta}} , $

      式中, $p({\boldsymbol{\theta}} )$表示先验概率密度函数, $p({\boldsymbol{d}})$为归一化因子, 与参数${\boldsymbol{\theta}}$无关, 似然函数$p({\boldsymbol{d}}|{\boldsymbol{\theta}} )$$L({\boldsymbol{\theta}} )$表示误差项的统计分布, 本文假设不同阶模态的频散曲线数据误差$ {{\boldsymbol{\varepsilon}} _m} $之间相互独立, 服从高斯分布$N({\boldsymbol{0}}, {{\boldsymbol{C}}_m})$, 因此, 似然函数$L({\boldsymbol{\theta}} )$可以表示为

      $ \begin{split} \;& L({\boldsymbol{\theta}} ) = \prod\limits_{m = 1,2,\cdots,M} {{L_m}({\boldsymbol{\theta}} )} \\ =\;& \prod\limits_{m = 1,2,\cdots,M} \frac{1}{{{{(2{\text{π }})}^{{N_m}/2}}|{{\boldsymbol{C}}_m}{|^{1/2}}}} \times \\ & \exp \left[ - \frac{1}{2}{{({{\boldsymbol{d}}_m} - {g_m}({\boldsymbol{\theta}} ))}^ \top }{\boldsymbol{C}}_m^{ - 1}({{\boldsymbol{d}}_m} - {g_m}({\boldsymbol{\theta}} )) \right] , \end{split} $

      其中, $^ \top $表示向量转置, ${{\boldsymbol{C}}_m}$$ {N_m} \times {N_m} $维的数据误差协方差矩阵, 通常假设${{\boldsymbol{C}}_m} = \sigma _m^2 {\boldsymbol{I}}$, 在地声参数反演中, 可以利用多次迭代的数据残差分析计算得到[17].

      获得环境参数的后验概率密度函数$p({\boldsymbol{\theta}} |{\boldsymbol{d}})$后, 其他统计结果可以利用(10)—(12)式计算得到, 包括一维边缘后验概率密度函数$p({\theta _i}|{\boldsymbol{d}})$、后验均值估计$ {{\boldsymbol{\theta}} _0} $和参数后验协方差矩阵${\boldsymbol{D}}$:

      $ p({\theta _i}|{\boldsymbol{d}}) = \int {\delta (\theta _i' - {\theta _i})p({{\boldsymbol{\theta}} '}|{\boldsymbol{d}}){\text{d}}{\boldsymbol{\theta}} '} , $

      $ {{\boldsymbol{\theta}} _0} = \int {{\boldsymbol{\theta}} 'p({{\boldsymbol{\theta}} '}|{\boldsymbol{d}}){\text{d}}{{\boldsymbol{\theta}} '}} , $

      $ {\boldsymbol{D}} = \int {({{\boldsymbol{\theta}} '} - {{\boldsymbol{\theta}} _0}){{({{\boldsymbol{\theta}} '} - {{\boldsymbol{\theta}} _0})}^ \top }p({{\boldsymbol{\theta}} '}|{\boldsymbol{d}}){\text{d}}{{\boldsymbol{\theta}} '},} $

      其中, $\delta $为狄拉克函数, 除上式外, 还可以采用最高概率密度可信区间来定量描述参数的不确定性, 例如, 95%可信区间表示包含边缘概率密度函数95%覆盖面积的最小宽度区间.

    • 对于非线性贝叶斯推断问题, 通常采用MCMC或VI方法近似求解, 当前向模型的计算复杂度较高时, 相比于MCMC方法(通常需要数以万计次的前向模型计算), 基于寻优策略的VI方法可以更加高效地完成后验分析工作, 本文不再赘述其基本理论, 感兴趣的读者可以参考相关文献[13], 这里仅给出散度(Kullback-Leibler divergence, KL)和证据下界(evidence lower bound, ELBO)的计算表达式, 分别为

      $ \text{KL}[{q}_{\phi }({\boldsymbol{\theta}} )||p({\boldsymbol{\theta}} |{\boldsymbol{d}})]={E}_{\phi }\left[{\rm{log}}\frac{{q}_{\phi }({\boldsymbol{\theta}} )}{p({\boldsymbol{\theta}} |{\boldsymbol{d}})}\right]\geqslant 0\text{, } $

      $\begin{split} {\text{ELBO}} =\;& {E_\phi }\left[ {\log \frac{{p({\boldsymbol{d}}|{\boldsymbol{\theta}} )p({\boldsymbol{\theta}} )}}{{{q_\phi }({\boldsymbol{\theta}} )}}} \right] \\ =\;& {E_\phi }[\log p({\boldsymbol{d}},{\boldsymbol{\theta}} )] + {E_\phi }[\log q_\phi ^{ - 1}({\boldsymbol{\theta}} )], \end{split}$

      其中, ${q_\phi }({\boldsymbol{\theta}} )$为参数化的分布族, 也称变分族, $\phi $表示变分参数, ${E_\phi }$表示关于$ {q_\phi } $求期望, $p({\boldsymbol{d}}, {\boldsymbol{\theta}} )$为联合概率密度函数.

      VI方法的核心思想是利用优化算法求出使ELBO取得最大值时的变分参数ϕ*, 此时后验概率密度函数$p({\boldsymbol{\theta}} |{\boldsymbol{d}})$可以被近似为$q_\phi ^*({\boldsymbol{\theta}} )$, 为了计算(14)式中非线性函数的复杂期望, Acerbi[14]在2018年提出了一种数据高效的近似推断方法, 也称VBMC方法, 该方法将参数化的分布族${q_\phi }({\boldsymbol{\theta}} )$建立为高斯混合模型(Gaussian mixed model, GMM), 并使用高斯过程模型(Gaussian process, GP)拟合$\log p({\boldsymbol{d}},{\boldsymbol{ \theta}} )$, 同时, 将基于主动采样的贝叶斯正交(Bayesian quadrature and active sampling)以及蒙特卡罗采样(Monte Carlo sampling)和VI方法相结合, 利用随机梯度下降法对证据下界函数进行寻优, 使其非常适合具有复杂似然函数的非线性贝叶斯推断问题, 大量仿真和实验结果表明, 该方法使用较小的计算代价可以获得和MCMC类方法几乎相同的后验分析准确度[15], 主要步骤如表1所示.

      算法 : 变分贝叶斯蒙特卡罗(Acerbi, 2018)
      1: $t \leftarrow 0$, 初始化GMM模型参数${\phi _0}$和${K_0}$
      2: When 收敛条件未满足, do
      3: $t \leftarrow t + 1$ //迭代
      4:  for $ 1 \ldots {n_{\text{active}} }$ do
      5:  采样新的样本点$\boldsymbol{\theta} _{\rm new} \leftarrow {\rm argmax}_{\boldsymbol{\theta}} [a({\boldsymbol\theta})]$ //$a({\boldsymbol{\theta }})$
         为采样函数
      6:  计算样本点的$ \log p(\boldsymbol{d}, {\boldsymbol{\theta} _{\text{new}} }) $值, 并将其加入训练集
      7: $\log p({\boldsymbol{d}}, {\boldsymbol{\theta}} ) \leftarrow$训练GP模型, 得到最优的模型超参数
      8: ${K_t} \leftarrow $更新GMM模型参数
      9: ${\phi _t} \leftarrow $利用随机梯度下降法对ELBO寻优, 更新GMM
        模型参数
      10: end

      表 1  VBMC算法步骤

      Table 1.  Steps of the VBMC algorithm.

    • 地声参数反演流程如图1所示, 假设模型的输入参数已知, 利用Kraken声场模型[18]计算出模态群速度$ v_{\text{g}}^{(m)} $后, 代入(5)式求出模态频散曲线, 即为地声参数反演中的正演模型$ {g_m}({\boldsymbol{\theta }}) $. 对于实际接收信号, 利用Warping变换[19]和同步压缩小波变换[20]分别对水波和底波信号的频散曲线进行分离和提取, 即得到地声参数反演中的观测数据$ {{\boldsymbol{d}}_m} $.

      图  1  地声参数反演流程图

      Figure 1.  Flow chart of geoacoustic inversion.

      在贝叶斯反演框架内(图1中虚线框所示), 首先, 设定待反演参数${\boldsymbol{\theta}}$的先验区间, 假设在区间上服从均匀分布; 然后, 基于正演模型${g_m}({\boldsymbol{\theta}} )$和观测数据${{\boldsymbol{d}}_m}$建立贝叶斯反演模型, 利用多次迭代的数据残差分析方法估计误差协方差矩阵${{\boldsymbol{C}}_m}$, 图中${{\boldsymbol{\theta}} _{{\text{ML}}}}$为参数的最大似然估计值, 在本文中使用自适应单纯形模拟退火算法求得[21]; 最后, 基于稳定的$ {{\boldsymbol{C}}_m} $估计结果, 采用VBMC方法推断后验概率密度函数$p({\boldsymbol{\theta}} |{\boldsymbol{d}})$, 完成参数后验分析. 在实际应用中, 对于最终获得的反演结果, 通常还需要采用多种直接或间接的方式对其进行有效性验证, 例如原位测量对比和传播损失对比等. 值得注意的是, 变分推断类方法对初始状态的选取较为敏感, 因此, 本文将${{\boldsymbol{\theta}} _{{\text{ML}}}}$作为VBMC算法的初始值.

    • 图2给出了一个具有3层水平分层介质的波导模型, 包括水体层、沉积层和半空间基底层: 在水体层中, 海水声速$ {c_{\text{w}}} $和海深$ {d_{\text{w}}} $均为常数; 在沉积层中, 压缩波声速随着沉积厚度的增大线性变化, 在$[0, {h_1}]$的范围内从$ {c_{1{\text{U}}}} $变化到$ {c_{1{\text{L}}}} $, 其密度${\rho _1}$为一个常数; 在基底半空间中, 压缩波声速$ {c_2} $和密度${\rho _2}$均为常数, 值得注意的是, 对于大多数浅海粉砂质砂或砂质粉砂海底, 忽略剪切波的影响是合理的, 所有的模型输入参数总结如表2所示.

      图  2  3层水平分层介质波导模型

      Figure 2.  Three-layer waveguide model.

      待反演参数单位仿真参数值先验区间
      沉积层厚度$ {h_1} $m47[10, 80]
      沉积层密度${\rho _1}$${\text{g/c}}{{\text{m}}^{\text{3}}}$1.5[1, 2]
      基底层密度${\rho _2}$${\text{g/c}}{{\text{m}}^{\text{3}}}$2.1[1.5, 2.5]
      沉积层上部声速$ {c_{1{\text{U}}}} $${\text{m/s}}$1675[1600, 1700]
      沉积层下部声速$ {c_{1{\text{L}}}} $${\text{m/s}}$1675[1650, 1750]
      基底层声速$ {c_2} $${\text{m/s}}$1765[1700, 1800]
      水深dwm21[15, 25]
      收发距离$r$km6.24[4, 8]
      时间因子$ {\delta _{\text{t} }} $s[–0.5, 0.5]
      水体声速cw${\text{m/s}}$1511[1500, 1520]

      表 2  反演参数空间和先验区间

      Table 2.  Parameter spaces and prior bounds for inversion.

      表2中的仿真参数值代入Kraken声场模型, 计算1阶模态群速度曲线, 如图3所示, 图中的圆圈给出Airy相频率的位置, 从整体上看, 群速度曲线被Airy相频率分成了左右两部分, 分别称为简正波传播中的底波和水波部分. 当$ {h_1} $, ${c_{{\text{1 L}}}}$, ${c_2}$${\rho _2}$取不同值时, 群速度曲线在小于Airy相频率的部分变化明显, 而在大于Airy相频率的部分基本不变, 说明底波的群速度特征对深层的地声参数变化较为敏感, 充分利用这一特性, 可以实现相关地声参数的可靠反演, 另外, 从图3的结果可以看出, $ {h_1} $${c_{{\text{1 L}}}}$的变化会引起Airy相频率的变化, 而${c_2}$${\rho _2}$的变化对Airy相频率无影响.

      图  3  1阶模态群速度曲线 (a)不同的$ {h_1} $; (b)不同的${c_{{\text{1 L}}}}$; (c)不同的${c_2}$; (d)不同的${\rho _2}$, 圆圈给出了艾里相频率的位置

      Figure 3.  Group velocity curves of Mode 1: (a) Different values of $ {h_1} $; (b) different values of ${c_{{\text{1 L}}}}$; (c) different values of ${c_2}$; (d) different values of ${\rho _2}$. The circles indicate the Airy phase frequencies.

    • 在仿真接收信号的频散特性时, 假设声源深度zs = 6 m, 接收深度zr = 11.5 m, 收发距离r = 6.24 km, 海底的衰减系数为0 $ {\text{dB/m}} $, 在不考虑环境噪声的情况下, 宽带脉冲信号的接收时域波形和时频分析结果如图4所示, 在0—500 Hz频带范围内, 接收信号包含前3阶简正波模态, 其中, 1阶和3阶模态信号能量较强, 2阶模态信号能量较弱. 利用图1所示的方法对1阶和3阶模态频散曲线进行分离和提取, 结果如图4(b)中的“×”符号所示, 可以看出, 在理想情况下, 实验提取的频散曲线与时频谱图显示的频散特征基本一致, 并且与理论计算出的频散曲线(图中的实线)基本吻合.

      图  4  仿真信号 (a)归一化时域波形; (b)时频图, 其中1阶模态的底波、水波和Airy相结构清晰可辩

      Figure 4.  Simulation signal: (a) Normalized time domain waveform; (b) time-frequency diagram. The ground wave, water wave, and Airy phase component of Mode 1 are well-defined.

      在实际中, 底波信号的幅度通常远小于水波信号的幅度, 因此, 实验中声源级大小、沉积物吸收以及环境噪声等因素都会影响底波的观测, 导致其频散曲线的提取存在误差, 综合考虑上述因素, 在仿真的底波信号中加入不同量级的高斯白噪声后, 进行频散曲线的提取, 分析不同信噪比(signal to noise ratio, SNR)条件下的提取误差, 结果如图5所示, 图5(a)对比了不同信噪比条件下的单次提取结果和理论模型计算结果; 假设提取误差在不同频点上服从独立的零均值高斯分布, 图5(b)给出了不同信噪比条件下误差标准差的估计结果, 可以看出, 随着信噪比逐渐降低, 误差标准差逐渐增大.

      图  5  不同信噪比条件下底波频散曲线的单次提取结果(a)和提取误差标准差(b)

      Figure 5.  Single extraction results (a) and standard deviations of extraction errors (b) of ground wave dispersion curve under different SNRs.

    • 为了讨论引入简正波的底波频散特性对地声参数反演的影响, 本文针对两种数据条件, 分别进行仿真研究: 数据条件1(Data 1), 包括1阶模态的底波和1阶, 3阶模态的水波频散特征数据; 数据条件2(Data 2), 包括1阶, 3阶模态的水波频散特征数据, 特别地, 对于数据条件1, 考虑不同的底波信噪比情况, 即SNR = 5, 0, –5和–10 dB.

      根据图1所示的流程进行地声参数反演, 部分深层参数的一维边缘后验概率密度函数如图6所示, 从图6(a)(d)分别为$ {h_1}, {c_{{\text{1 L}}}}, {c_2}, {\rho _2}$, 图中的点划线表示参数真值, 可以发现: 1)当底波信号信噪比较高时(SNR = 5 dB, 0 dB), 数据条件1的边缘后验概率密度函数较为尖锐, 即反演的不确定性较小, 且最大后验概率估计值(概率密度函数的最大值点)与参数真值较为接近, 相比之下, 数据条件2的边缘后验概率密度函数非常平坦, 几乎等于先验概率密度函数, 因此, 无法对该参数进行可靠估计; 2)随着底波信噪比的降低, 其频散曲线的提取误差标准差逐渐增大(如图5(b)所示), 包含的有用信息越来越少, 数据条件1的地声参数反演效果将逐渐退化为数据条件2的地声参数反演效果. 总之, 当底波信号的信噪比足够高时, 基于简正波水波和底波完整频散特性的贝叶斯地声参数反演可以准确估计出深层地声参数.

      图  6  两种数据条件下部分参数的一维边缘后验概率密度函数 (a)$ {h_1} $; (b)${c_{{\text{1 L}}}}$; (c)${c_2}$; (d)${\rho _2}$

      Figure 6.  1 D marginal posterior probability densities of some parameters for Data 1 and Data 2: (a) $ {h_1} $; (b) ${c_{{\text{1 L}}}}$; (c) ${c_2}$; (d) ${\rho _2}$.

    • 图7给出某次渤海实验的实验地点和实测的水体声速剖面(sound speed profile, SSP), 在图7(a)中, 符号▲表示11号自容式水听器的具体位置, 深度传感器显示其布放深度约为11.5 m, 其灵敏度约为$- 170{\text{ dB re }}1\; {{\rm{V}}/\text{μ} {\rm{Pa}}}$, 采集系统的采样率为30 kHz, 实验过程中, 宽带脉冲声源在图中★符号所示位置处发射信号, 发射深度约为6 m, 发射距离约为6.07 km, 利用温盐深仪对水体声速剖面进行测量, 结果如图7(b)所示, 测点处的海深约为20 m, 海面附近的声速约为1512.5 m/s, 海底附近的声速约为1511 m/s, 海水混合程度较高, 介质较均匀.

      图  7  实验描述 (a)实验地点; (b)水体声速剖面

      Figure 7.  Description of the experiment: (a) Experimental site; (b) sound speed profile in water.

    • 自容式水听器接收信号的时域波形和时频分析结果如图8所示, 信号包含两阶简正波模态, 通过仿真对比可以确定分别为1阶和3阶模态, 其中, 1阶模态包含水波和底波两组成分, 而3阶模态仅包含水波成分, 1阶模态的底波成分到达接收点的时间最早, 准正弦的底波信号能量主要集中在35—74 Hz之间, 大幅值的水波信号紧随其后, 与底波信号相比, 水波信号的波形更加尖锐, 频散特征更弱, 1阶模态信号中到达时间最晚的部分(图8(b)中方框区域)即为Airy相结构, 对应的Airy相频率约为74 Hz(点划线). 实验中, 接收信号的信噪比较高, 有利于进行模态频散特性分析.

      图  8  实验信号 (a)归一化时域波形; (b)时频图, 其中1阶模态的底波、水波和Airy相结构清晰可辩

      Figure 8.  Experimental signal: (a) Normalized time domain waveform; (b) time-frequency diagram. The ground wave, water wave, and Airy phase component of Mode 1 are well-defined.

      利用Warping变换和同步压缩小波变换对水波信号和底波信号的频散曲线分别进行提取, 结果如图8(b)中的“×”符号所示, 频散曲线和时频图中的频散特征一致性较好, 值得注意的是, 宽带声源信号在频谱上存在弱干涉, 导致接收信号的模态能量在感兴趣的频段范围内出现强弱交替的现象, 但是, 其对频散特征的影响较小, 实际处理时可以忽略.

    • 与仿真研究相同, 针对两种实验数据条件分别进行地声参数反演, 如图9所示, 图中使用误差条给出了实验数据误差标准差${\sigma _m}$的估计结果.

      图  9  频散曲线的实测结果和后验预测结果 (a)数据条件1, 水波和底波频散数据, 插图为低频部分的局部放大; (b)数据条件2, 水波频散数据

      Figure 9.  Measurements and posterior predictions of dispersion curves: (a) Data 1, water waves and ground waves, the enlarged view of the low frequency part is inset; (b) Data 2, water waves.

      待反演参数的边缘后验概率密度函数如图10所示, 对角线图形给出了一维估计结果, 实线和虚线分别表示数据条件1和条件2, 在先验区间内, 不同数据条件下的10个参数结果均表现为单峰结构, 但是, 因为引入了底波频散特征, 数据条件1的参数估计不确定性(概率密度函数的分布范围)明显小于数据条件2, 特别是对于深层地声参数, 这一结果尤为明显, 除此之外, 联合水波和底波完整频散特性进行地声参数反演, 可以准确地分辨沉积层和基底半空间层, 即$ {h_1} $的估计不确定性非常小, 综上所述, 底波频散特征数据可以提供更多关于海底声学特性的信息, 使相关参数的反演更加准确.

      图  10  待反演参数的一维(数据条件1和条件2)和二维(数据条件1)边缘后验概率密度函数

      Figure 10.  1D (Data 1 and Data 2) and 2D (Data 1) marginal posterior probability densities of unknown parameters.

      利用待反演参数的后验估计结果对模态频散曲线进行预测, 如图9所示, 虚线给出后验均值样本的预测值, 细实线给出100个后验随机样本的预测值, 在不确定性范围内, 预测结果和实验数据的一致性较好, 充分说明图2所示的浅海波导模型可以有效模拟此次海上实验的真实环境.

      针对数据条件1的地声参数反演给出了待反演参数的可靠估计, 后面重点对其进行讨论. 参数的后验统计结果总结在表3中, 包括均值、方差和95%可信区间, 其中, 沉积层参数结果表明, 实验海域的海底覆盖一层中高声速沉积物, 平均厚度约为(46.65±3.77) m, 上部的压缩波声速约为(1675.89±9.40) m/s, 下部的压缩波声速约为(1673.52±8.04) m/s, 平均密度约为(1.52±0.09)${\text{g/c}}{{\text{m}}^{\text{3}}}$, 沉积层中介质的低频声学特性随着沉积厚度的增大变化不大. 在基底层中, 压缩波声速约为(1765.03±15.78)${\text{m/s}}$, 其95%的可信区间为[1731.16, 1789.77] m/s, 密度约为(2.12±0.19) ${\text{g/c}}{{\text{m}}^{\text{3}}}$, 其95%可信区间为[1.75, 2.44] ${\text{g/c}}{{\text{m}}^{\text{3}}}$. 本文假设环境参数与距离无关, 因此海水深度$ {d_{\text{w}}} $和声速$ {c_{\text{w}}} $代表了实验海域的平均估计结果, 分别为(20.69$ \pm $0.45) m (实测结果为19—22 m)和(1510.98$ \pm $4.13) m/s (实测结果为1511—1513 m/s), 收发距离$r$的估计结果为(6.24±0.3) km(实测结果为6.07 km), 时间因子$ {\delta _{\text{t}}} $作为一个辅助变量, 在反演过程中快速收敛, 并且估计不确定性较小, 其均值为0.280 s, 方差为0.002 s.

      待反演参数均值方差95%可信区间实测值
      $ {h_1} $46.653.77[39.69, 54.21]
      ${\rho _1}$1.520.09[1.35, 1.69]
      ${\rho _2}$2.120.19[1.75, 2.44]
      $ {c_{1{\text{U}}}} $1675.899.40[1654.35, 1690.32]
      $ {c_{1{\text{L}}}} $1673.528.04[1660.41, 1691.35]
      $ {c_2} $1765.0315.78[1731.16, 1789.77]
      $ {d_{\text{w}}} $20.690.45[19.70, 21.41]19—22
      $r$6.240.30[5.63, 6.83]6.07
      $ {\delta _{\text{t}}} $0.2800.002[0.276, 0.284]
      $ {c_{\text{w}}} $1510.984.13[1503.08, 1517.91]1511—1513

      表 3  针对数据条件1的后验统计结果

      Table 3.  Posterior statistical results for Data 1.

      待反演参数之间的相关性可以通过二维边缘后验概率密度函数(图10中的非对角线图形)和归一化后验参数协方差矩阵(图11中的热图)进行描述, 结果显示, 不同参数之间存在一定的相关性, 例如: 水平距离$r$和平均水深$ {d_{\text{w}}} $之间存在明显的正相关; 水平距离$r$和时间因子$ {\delta _{\text{t}}} $之间存在明显的负相关, 因此, 当其中一个参数的估计出现偏差时, 与之相关的其他参数的估计也会出现偏差.

      图  11  针对数据条件1的归一化后验参数协方差矩阵

      Figure 11.  Normalized posterior covariance matrix of unknown parameters for Data 1.

      图12(a), (b)分别给出VBMC推断过程中ELBO和KL散度的收敛情况, 可以看出, 在经过约80次迭代之后, 该方法已经基本收敛, 此时前向模型的计算次数仅为400余次, 与MCMC方法相比, 计算次数显著减少, 虽然VBMC推断还包括有梯度下降算法等额外的计算损耗, 但是其总的后验分析时间是远小于MCMC方法的. 利用VBMC方法求出近似的变分后验概率密度函数后, 通过简单的采样算法获得其有效样本, 并进行统计分析.

      图  12  针对数据条件1的VBMC收敛情况 (a)证据下界ELBO; (b) KL散度

      Figure 12.  Convergences of the VBMC method for Data 1: (a) ELBO; (b) KL divergence.

      进一步讨论海底的分层情况, 针对数据条件1, 假设在沉积层和基底层之间还存在一个中间层, 该层的厚度、压缩波声速和密度均为未知常数, 同样, 利用2.4节所述方法, 可以获得如图13(b)所示的波导环境声学参数剖面, 包括压缩波声速(左)和密度(右), 图中粗实线表示后验均值结果, 细实线表示若干个后验随机结果, 图13(a)给出没有中间层的波导环境声学参数剖面, 对比发现, 两种海底分层结构下的波导环境几乎相同, 海水的实测声速剖面(虚线)和平均估计结果(实线)均拟合较好, 海底中间层的出现等效于在50 m深度附近将沉积层分成上下两部分, 除了下部的密度外, 两部分的底质声学特性基本一致.

      图  13  针对数据条件1的压缩波声速和密度剖面后验估计结果 (a) 2层海底模型; (b) 3层海底模型

      Figure 13.  Posteriori estimates of compressed-wave sound speed profiles and density profiles for Data 1: (a) Two-layer seabed model; (b) three-layer seabed model.

    • 表3给出的参数均值代入Kraken声场模型, 计算模态群速度, 结合另一组已知收发距离的宽带脉冲信号进行频散曲线预测, 图14给出了不同距离上接收信号的时频分析结果, 从图14(a)(c)分别为3.42, 5.35和6.51 km, 实线为频散曲线的预测结果, 两者的一致性较好, 充分证明上述地声参数反演结果的有效性和可靠性.

      图  14  另一组宽带脉冲信号时频图和频散曲线预测结果(实线) (a)接收距离3.42 km; (b)接收距离5.35 km; (c)接收距离6.51 km

      Figure 14.  Time-frequency diagrams and dispersion curve predictions (solid lines) of other broadband pulse signals: (a) The range is 3.42 km; (b) the range is 5.35 km; (c) the range is 6.51 km.

      历史调查和研究显示, 实验海域的海深在15—25 m之间, 海底沉积物多为粉砂质砂或砂质粉砂[22,23], 平均粒径的变化范围为[2$\varPhi $, 5$\varPhi $], 文献[10, 24]给出了相似浅海环境下粉砂质沉积物低频衰减系数$\alpha $的取值, 当频率为200 Hz时, $\alpha \approx 0.02{\kern 1 pt} {\kern 1 pt} {\kern 1 pt} {\text{dB/m}}$, 当频率为315 Hz时, $\alpha \approx 0.06{\kern 1 pt} {\kern 1 pt} {\kern 1 pt} {\kern 1 pt} {\text{dB/m}}$. 声源深度为6 m, 接收深度为11.5 m, 结合其他地声参数反演结果进行传播损失预报, 如图15所示, 传播损失的距离变化范围为[0, 10] km, 图中圆圈为实测传播损失, 距离分别为3.42, 5.35和6.51 km, 可以看出, 利用反演得到的地声参数(图中实线)可以较好地预测实验海域的传播损失.

      图  15  传播损失的理论计算结果和实测结果比较, 声源深度6 m, 接收深度11.5 m (a)频率200 Hz; (b)频率315 Hz

      Figure 15.  Comparisons of the theoretical and experimental transmission loss at source depth 6 m and receiver depth 11.5 m: (a) The frequency is 200 Hz; (b) the frequency is 315 Hz.

    • 针对基于简正波水波频散特性的地声参数反演方法无法有效估计深层地声参数的现实问题, 本文联合底波和水波频散特性进行贝叶斯地声参数反演, 创新地采用变分贝叶斯蒙特卡罗算法对反演结果进行快速后验分析, 不仅实现了深层地声参数的可靠估计, 而且大大降低了浅层地声参数和水体参数的估计不确定性, 仿真和实验结果显示: 待反演参数的边缘后验概率密度函数表现为先验区间内的尖峰结构, 频散曲线的预测结果与测量结果拟合较好, 说明通过此方法反演出的波导环境模型比较可靠. 在实验数据验证中, 收发距离和海水平均声速两个参数的估计结果与实际测量结果基本一致, 反演得到的地声参数不但可以用于解释本组实验数据, 而且可以较好地用于不同距离上另一组实验数据的模态频散特征预测, 除此之外, 传播损失的准确预报也充分证明了联合反演方法的有效性.

      在本文的地声参数反演研究中, 假设海洋环境参数不随距离发生变化, 然而, 实际的浅海环境复杂多样, 在水平和垂直方向上均表现出非均匀性, 因此, 如何实现距离相关浅海波导环境下地声参数的可靠反演是下一步需要研究的工作.

参考文献 (24)

目录

    /

    返回文章
    返回