| Citation: | Xiaoman Li, Biao Wang, Xuejie Bi, Hong Wu. A fast inversion method for ocean parameters based on dispersion curves with a single hydrophone[J]. Acta Oceanologica Sinica, 2022, 41(9): 71-85. doi: 10.1007/s13131-022-1999-z
|
Propagation characteristics of underwater sound sources are influenced greatly by the marine environmental characteristics. The inversion of ocean properties in situ is the basis of exploring marine structure and analyzing marine characteristics, and it is also a valuable field in underwater acoustic application research (Bonnel et al., 2013). Various inversion methods of ocean parameters have been developed (Li et al., 2019b; Bao et al., 2017), and the use of an array of hydrophones is the most common method. In methods with a hydrophone array, the water column can be covered with hydrophones (vertical array or horizontal array), the measured data contains the information of the whole sound field space, so marine environment parameters can be accurately estimated. However, in practical applications, the cost of hydrophone array is high, what’s more, the array form of hydrophone array can be affected greatly by ocean internal motion (such as ocean current), which will greatly reduce the reliability of measured data (Walker et al., 2005; Gingras and Gerstoft, 1995; Fallat and Dosso, 1999). Therefore, passive positioning or ranging methods of underwater targets with a hydrophone have become a major research trend in the field of underwater acoustics (Le Gac et al., 2003; Heaney, 2004; Bonnel et al., 2012). The above methods have been greatly improved in terms of cost saving and application systems.
Compared with the technology in a hydrophone array, the information received by single hydrophone is lower in quantity. Therefore, more prior knowledge of sound fields and more accurate information about sound field extraction are necessary. In addition, a mismatch of environmental parameters will greatly affect the accuracy of underwater acoustic parameter inversion. The Bayesian methodology can realize the geoacoustic inversion, and can get more accurate inversion results (Le Touze et al., 2009). The Bayesian inversion methodology can take unknown ocean parameters as inversion parameters, which requires less prior knowledge of environment and has good environmental tolerance (Dosso and Wilmut, 2008).
The Bayesian inversion theory is applied in underwater acoustic research. In an uncertain ocean environment, Bayesian focalization and marginalization method is applied to localization of multiple sound sources (Dosso and Wilmut, 2011). The multi-source localization and continuous localization of moving sound sources in an uncertain marine environment can also be realized based on Bayesian theory (Li et al., 2016, 2018). The inversion of ocean parameters is realized by matching the extracted dispersion curves with the result calculated by the acoustic model (Bonnel et al., 2013). The marine parameters can also be inverted by matching mode signals extracted from a received signal with calculated replicas from an acoustic model (Li et al., 2019a). However, there are two shortcomings in the application of Bayesian estimation in underwater acoustic parameter inversion. On the one hand, when there are many inversion parameters, the correlation between different parameters will lead to the multiple values in inversion results; on the other hand, when the acoustic model is used to calculate the replica, the final estimation results will be affected by the accuracy and completeness of the model, and the calculation process is complicated and uncontrollable.
In this paper, a fast inversion method comprising marine environment parameters based on Bayesian methodology with single hydrophone in shallow water is presented. The input function of Bayesian methodology is the dispersion curves, which are extracted from the received signal by a signal processing tool named warping operator in time-frequency domain (Bonnel et al., 2020). Regarding the surface-reflected bottom reflected (SRBR) modes, when the sound source is incidental at a small glancing angle, the bottom characteristics can be described by a parameter-the bottom reflection phase shift parameter P based on the Beam-Displacement Ray-Mode (BDRM) theory (Shang et al., 2012; Niu et al., 2014). Therefore, a new dispersion formula can be used to calculate the replica, which is different from the use of an acoustic model. The dispersion formula is composed of four marine parameters: depth, propagation range, average velocity of sound in water and P. The parameters are optimized by establishing an effective cost function, and the maximum value of the one-dimensional marginal probability density (MPD) of the parameters is taken as the inversion result, and a fast inversion method of the marine environment parameters based on Bayesian methodology with single hydrophone is thus realized. In this way, only three parameters are involved in the inversion process, and the inversion speed and efficiency are both improved.
This paper has five parts: Section 2 describes the extraction and estimation methodology for the dispersion curves. Section 3 presents the inversion scheme, while Section 4 describes the inversion results of measured data. At last, Section 5 elaborates the conclusions.
Normal mode theory is the basic theory of acoustic signal propagation, when a broadband impulse sound source S(f) is propagated in shallow water waveguide which is independent of range, the measured signal P(f,r) is the sum of normal modes, which can be written as (Jensen et al., 2011)
| $$ P(f,r) = S(f)\sum\limits_{n = 1}^N {\frac{{{{\rm{e}}^{{\rm{j}}\pi /4}}}}{{\sqrt {8\pi } \rho ({z_{\rm{s}}})}}} {U_n}({z_{\rm{s}}},f){U_n}({z_{\rm{r}}},f)\frac{{{{\rm{e}}^{{\rm{j}}{\mu _n}{(}f{)}r{ - }{\beta _n}{(}f{)}r}}}}{{\sqrt {{\mu _n}r} }} , $$ | (1) |
where zs, zr, ρ and r are the depth of sound source, received hydrophone, density in water and propagation range; the maximum number of normal modes is N; Eigen-function Un, a function related to depth; μn(f) is the horizontal wavenumber, while βn(f) is the attenuation coefficient.
In received signal P(f,r), each normal mode signal carries a large amount of marine environmental information in the process of propagation. The propagation speed (group speed) of each normal mode is related to the frequency, which makes the arrival time of each normal mode varies with frequency. The energy of each normal mode is concentrated around the dispersion curve, which describes the relationship between arrival time of each mode and frequency. The dispersion curve of the nth normal mode is given by:
| $$ {t_n}(f) = \frac{r}{{{v_{{\rm{g}}n}}(f)}}{\text{ = }}\frac{{r\partial {\mu _n}}}{{\partial (2\pi f)}} , $$ | (2) |
where vgn(f ) is the group speed of the nth normal mode. The energy of the received signal are concentrated on different curves that change with frequency in the time-frequency (TF) domain after signal processing analysis (Cai et al., 2021). Figure 1a shows the spectrogram of the measured signal after time-frequency analysis in TF domain and the dispersion curves (black solid line), the time-frequency analysis technology named short time Fourier transform (STFT) is used here. It can be seen from Fig. 1a that the energy of each mode changes along the curve, and the change trend of different modes is different, which fully reflects the dispersion characteristics.
For inversion system, dispersion curves are good choices to be used as an input function (Bonnel et al., 2013). To extract the dispersion curves, time-frequency analysis (TFA) technology is used as a tool for signal preprocessing. In time-frequency domain, the normal mode signals can be separated and extracted by warping operator that acts on the nonlinear phase of the normal mode. The dispersion curves can then also be extracted by a data extraction technique in the TF domain after the warping operator (Bonnel et al., 2020). The received sound pressure field signal p(t) can be transformed into a new signal
| $$ {\mathbf{W}}p(t) = \sqrt {\left| {\frac{{\partial h(t)}}{{\partial t}}} \right|} p\left[ {h(t)} \right] , $$ | (3) |
where
| $$ h(t) = \sqrt {{t^2} + t_r^2} , $$ | (4) |
| $$ {h^{ - 1}}(t) = \sqrt {{t^2} - t_r^2} , $$ | (5) |
where
The whole procedure for extracting mode curves is described in Fig. 1. The spectrum of the received signal for a broadband impulse source with a frequency band of 100 Hz to 300 Hz propagating in shallow water is shown in Fig. 1a. when the parameters r = 20 km and c0 = 1 500 m/s are used, the warped result in TF domain is presented in Fig. 1b. The 3rd mode extracted from the received signal is taken as an example, and the spectrum and the dispersion curves are both shown in Figs 1c and d. The results have the same trajectory. According to the above method, all normal mode curves of the received signal can be extracted.
As stated above, dispersion curves contain a substantial amount of environmental information, thus providing a basis for the inversion of environmental parameters. This method performs well in calculating dispersion curves with an acoustic model and detailed knowledge of the environment. However, for onsite applications, acoustic parameters are difficult to obtain, and unknown parameters are used as inversion parameters. The use of too many inversion parameters not only requires extensive calculation but also leads to interactions between parameters, thus decreasing the accuracy of the replica. The accuracy of the acoustic model also affects the final inversion results. To address the above problems, the mode dispersion curves are calculated with a dispersion formula instead of the acoustic model herein.
In terms of SRBR modes, the sound speed profile in water can be expressed as follows (Wang et al., 2016):
| $$ c(z) = \bar c[1 - a(z)] , $$ | (6) |
where
| $$ k(z) = \omega /c(z) \approx \bar k[1 + a(z)] , $$ | (7) |
where
In BDRM theory (Niu et al., 2014), the mode is propagated by an eigen-ray, which is shown in Fig. 2.
The received signal is the sum of WKBZ modes (Niu et al., 2014), and the Eigen-equation (EE) can be written as:
| $$ 2\int_{{z_{n1}}}^{{z_{n2}}} {\sqrt {{k^2}(z) - \mu _n^2} {\rm{d}}z + } {\varphi _1}(\mu _n^{}) + {\varphi _2}(\mu _n^{}) = 2(n - 1)\pi ,{\text{ }}n = 1,2,3,... $$ | (8) |
where φ1 and φ2 are the reflection phases on surface and seabed boundaries, respectively. In Eqs (9) and (10),
| $$ {\text{δ}_{n1}}({\mu _n}) = - \left. \frac{\partial {\varphi _1}(\mu _n)}{\partial \mu _n} \right|_{\omega} ,\quad {\text{δ} _{n2}}({\mu _n}) =-\left. \frac{\partial {\varphi_2}(\mu_n)}{\partial \mu _n} \right|_{\omega}, $$ | (9) |
| $$ {\tau _{n1}}({\mu _n}) = \left. \frac{{\partial {\varphi _1}({\mu _n})}}{{\partial \omega }}\right|_{{\mu _n}},\quad {\tau _{n2}}({\mu _n}) = \left. \frac{\partial {\varphi _2}({\mu _n})}{\partial \omega } \right|_{{\mu _n}} . $$ | (10) |
The cycle distance of eigen-ray is:
| $$ {S_n} \approx \frac{{2H{\mu _n}}}{{\sqrt {{{\bar k}^2} - \mu _n^2} }} . $$ | (11) |
And the cycle time is:
| $$ {T_n} \approx \frac{{2H\bar k}}{{\bar c\sqrt {{{\bar k}^2} - \mu _n^2} }} . $$ | (12) |
The approximate expression of the horizontal wave number
| $$ {\mu _n} \approx \frac{{\bar k{r_0}}}{{\bar ct'}} , $$ | (13) |
where
Considering the small incident angle (
| $$ 2\int_0^H {\sqrt {{k^2}(z) - u_n^2} {\rm{d}}z + {\phi _b} - \pi = 2n\pi } , $$ | (14) |
where θ is the incident angle, and
| $$ {\phi _b} \cong - \pi + P\theta . $$ | (15) |
In Eq. (15), P is independent of frequency. For a waveguide with a liquid, uniform, half-infinite bottom, P is (Wang and Shang, 2013)
| $$ P \approx \dfrac{{{{\sqrt 2 {\rho _2}} \mathord{\left/ {\vphantom {{\sqrt 2 {\rho _2}} {{\rho _1}}}} \right. } {{\rho _1}}}\sqrt {\sqrt {{{\left( {1 - {{\left( {{{{c_{\text{0}}}} \mathord{\left/ {\vphantom {{{c_{\text{0}}}} {{c_b}}}} \right. } {{c_b}}}} \right)}^2}} \right)}^2} + {{\left( {{{{c_0}} \mathord{\left/ {\vphantom {{{c_0}} {{c_b}}}} \right. } {{c_b}}}} \right)}^2}\dfrac{{\alpha {c_b}}}{{2\pi f}}} + 1 - {{\left( {{{{c_0}} \mathord{\left/ {\vphantom {{{c_0}} {{c_b}}}} \right. } {{c_b}}}} \right)}^2}} }}{{\sqrt {{{\left( {1 - {{\left( {{{{c_0}} \mathord{\left/ {\vphantom {{{c_0}} {{c_b}}}} \right. } {{c_b}}}} \right)}^2}} \right)}^2} + {{\left( {{{{c_0}} \mathord{\left/ {\vphantom {{{c_0}} {{c_b}}}} \right. } {{c_b}}}} \right)}^2}\dfrac{{\alpha {c_b}}}{{2\pi f}}} }} , $$ | (16) |
where ρ1 and ρ2 are the water density and the bottom density; c0 and cb are the sound speed; α is the seabed absorption parameter. When the source is incidental at a small grazing angle θ, and the seabed absorption coefficient is small, the influence of frequency on P can be ignored in a certain high-frequency range (greater than tens of Hertz). P can be simplified as
| $$ P \approx \frac{{2{\rho _2}}}{{{\rho _1}\sqrt {1 - {{\left( {{{{c_{\text{0}}}} \mathord{\left/ {\vphantom {{{c_{\text{0}}}} {{c_b}}}} \right. } {{c_b}}}} \right)}^2}} }} . $$ | (17) |
When Eq. (15) is substituted into Eq. (14), a Eigen-equation can be obtained,
| $$ 2{\displaystyle {\int }_{ 0}^{H}\sqrt{{k}^{2}(z)-{\mu }_{n}^{2}}{\rm{d}}z+P{\theta }_{n}=2n\pi ,\;\;}n=1,2,3,\mathrm{...} $$ | (18) |
where
| $$ 2\int_0^H {k(z)\sin {\theta _n}{\rm{d}}z + P\sin {\theta _n} = 2n\pi } ,{\text{ }}n = 1,2,3,... $$ | (19) |
Eq. (19) can also be expressed as
| $$ 2\int_0^{H + \Delta H} {\sqrt {{k^2}(z) - u_n^2} {\rm{d}}z = 2n\pi } , $$ | (20) |
where
| $$\left(H + \frac{P}{{2\bar k}}\right)\bar k\sqrt {1 - \left( {\frac{{{r_0}}}{{\bar ct(1 + \varepsilon )}}} \right)_{}^2} = n\pi . $$ | (21) |
The instantaneous frequency expression of the nth normal mode is obtained,
| $$ f(t') = \frac{{n\bar {ct}}}{{2H\sqrt {t'^2 - {{\left( {\dfrac{r}{{\bar c}}} \right)}^2}} }} - \frac{{P\bar c}}{{4\pi H}}, $$ | (22) |
where
| $$ t' = t(1 + \varepsilon ). $$ | (23) |
Therefore, the incident angle can also be written approximately as
| $$ \theta \approx \sin \theta \approx \frac{{\sqrt {{{\bar k}^2} - \mu _n^2} }}{{\bar k}} \approx \frac{{\sqrt {{t^2} - \left( {{r \mathord{\left/ {\vphantom {r {\bar c}}} \right. } {\bar c}}} \right)_{}^2} }}{t} . $$ | (24) |
The bottom reflection phase shift φb is:
| $$ {\varphi _b} \approx - \pi + P \cdot \frac{{\sqrt {{t^2} - \left( {{r \mathord{\left/ {\vphantom {r {\bar c}}} \right. } {\bar c}}} \right)_{}^2} }}{t} . $$ | (25) |
Submitting the Eq. (25) into Eqs (9) and (10), the beam displacement and the time delay of the seabed is
| $$ {\text{δ} _2}({\mu _n}) = \frac{{{{\varphi '}_b}t}}{{\mu _n^{}}}{\text{,}}\quad {\tau _2}({\mu _n}) = \frac{{{{\varphi '}_b}r}}{{\mu _n^{}{{\bar c}^2}}} . $$ | (26) |
After combining Eqs (26), (11), and (12), the expression of ε becomes
| $$ \varepsilon = \frac{{P({t^2} - \left( {{r \mathord{\left/ {\vphantom {r {\bar c}}} \right. } {\bar c}}} \right)_{}^2)}}{{2\bar kH{t^2}}} . $$ | (27) |
According to Eq. (23), Eq. (22) can be expanded with the Taylor series at t, and when the higher-order minor term is ignored, Eqs (27) and (17) are submitted to Eq. (22), and the dispersion formula of the nth normal mode is
| $$ f(t) \approx \frac{{n\pi \bar ct}}{{2\pi H\sqrt {{t^2} - \left( {{r \mathord{\left/ {\vphantom {r {\bar c}}} \right. } {\bar c}}} \right)_{}^2} }} - \frac{{P\bar c}}{{4\pi H}} - \frac{{Pr}}{{4\pi H{t^2}\bar c}}, $$ | (28) |
| $$\begin{split} f(t) \approx &\frac{{n\pi \bar ct}}{{2\pi H\sqrt {{t^2} - \left( {{r \mathord{\left/ {\vphantom {r {\bar c}}} \right. } {\bar c}}} \right)_{}^2} }} - \frac{{{{2{\rho _2}} \mathord{\left/ {\vphantom {{2{\rho _2}} {{\rho _1}}}} \right. } {{\rho _1}}}\sqrt {1 - {{\left( {{{\bar c} \mathord{\left/ {\vphantom {{\bar c} {{c_b}}}} \right. } {{c_b}}}} \right)}^2}} \bar c}}{{4\pi H}} -\\ &\frac{{{{2{\rho _2}} \mathord{\left/ {\vphantom {{2{\rho _2}} {{\rho _1}}}} \right. } {{\rho _1}}}\sqrt {1 - {{\left( {{{\bar c} \mathord{\left/ {\vphantom {{\bar c} {{c_b}}}} \right. } {{c_b}}}} \right)}^2}} r}}{{4\pi H{t^2}\bar c}}.\end{split} $$ | (29) |
When the dispersion Eq. (28) is used to calculate dispersion curves, some marine parameters are used. The bottom reflection phase shift parameter P can describe the seabed characteristics, thus making subsequent inversion work more convenient. However, when more detailed marine environmental parameters are required, Eq. (29) can also be used to invert the seabed parameters.
The average speed of sound in seawater is an important parameter in the dispersion formula. In practice, the sound speed profile of waveguide changes with depth. To verify the validity and accuracy of Eqs (28) and (29) in calculating the dispersion curves, the simulation analysis process was performed in a shallow water waveguide with thermocline and isovelocity sound speed profiles respectively. The simulation parameters are shown in Fig. 3, and the depth is H=100 m. A thermocline is present in the water, and the average sound speed is c0=1 531.5 m/s, which is the sound speed of the isovelocity waveguide. A chirp sound source with a bandwidth of 50–200 Hz is emitted at a depth of zs=90 m, and the receiver is located on the seabed.
Three methods are used to calculate the dispersion cures concurrently. KRAKEN is a widely used acoustic model; therefore, the results calculated with KRAKEN are used as a reference value (Porter, 1992). Dispersion curves are calculated with the KRAKEN program (method 1), Eq. (29) (method 2), and the method described in Section 2.2 (method 3). The calculation results of methods 2 and 3 are compared with the results of the first method. The curves obtained with the three methods are represented by black solid lines, blue solid lines with circles, and red dots, respectively. The calculation results of the dispersion curves for an isovelocity waveguide are shown in Fig. 4, and the results for a thermocline waveguide are shown in Fig. 5.
As shown in Fig. 4, the results obtained with the three methods are in good agreement, thus indicating that the dispersion formula is effective for the isovelocity waveguide. In the results of the waveguide with thermocline shown in Fig. 5, the dispersion curves of normal mode with the number greater than two obtained with the three methods tend to be consistent. However, for the first mode, the data extraction results (method 3) are essentially consistent with the KRAKEN calculation results (method 1), whereas the results obtained with the dispersion formula (method 3) appear to deviate (as shown in the enlarged part of the figure). This result is because the derivation of the dispersion formula is based on SRBR normal modes. In shallow water waveguides with a thermocline, the low-order mode (particularly the first mode) does not meet the conditions of an SRBR normal mode in many cases; therefore the dispersion formula herein would not be applicable. However, for shallow water waveguides, the received signal would contain many order normal modes. Therefore, if the first mode is discarded, Eqs (28) and (29) can be used instead of the acoustic model to calculate the dispersion curves of the SRBR modes in shallow water, which in turn can be used to calculate a replica of the dispersion curves.
As stated above, dispersion curves contain a significant amount of environmental information; therefore, they can be used as an input signal to an inversion methodology of ocean parameters. Bayesian inverse theory is used as the inverse algorithm in this paper.
The parameter vector m is inverted when d represents measured data and E(m) is the data misfit function. The generalized misfit function is usually defined by combining data and a prior information,
| $$ \varphi ({\boldsymbol{m}}){{ = }}E({\boldsymbol{m}}){{ - }}\ln Pro({\boldsymbol{m}}) , $$ | (30) |
where Pro(m) is the prior probability density function of the model parameters, and is generally assumed to satisfy a uniform distribution. The posterior probability density (PPD) is usually used to describe parameter estimates and uncertainty. The conditional probability density Pro(m|d) can be expressed as (Dosso, 2002)
| $$ Pro({\boldsymbol{m}}\left| {\boldsymbol{d}} \right.) = \frac{{\exp [ - \varphi ({\boldsymbol{m}})]}}{{\int\limits_{} {\exp [ - \varphi ({\boldsymbol{m'}},{\boldsymbol{d}})]{\rm{d}}{\boldsymbol{m'}}} }} , $$ | (31) |
where φ(m) is the cost function. The integration interval is an M-D parameter space. To obtain the probability distribution of the parameters in model i, the energy values of each parameter can be weighted according to the Boltzmann function,
| $$ Pro({{\boldsymbol{m}}_i}|{\boldsymbol{d}}) = \frac{{\exp ( - [\varphi ({{\boldsymbol{m}}_i}) - \varphi ({{\boldsymbol{m}}_{i0}})]/[t - \varphi ({{\boldsymbol{m}}_{i0}})])}}{{\sum\limits_{k = i}^M {\exp ( - [\varphi ({{\boldsymbol{m}}_i}) - \varphi ({{\boldsymbol{m}}_{i0}})]/[t - \varphi ({{\boldsymbol{m}}_{i0}})])} }} , $$ | (32) |
where M is the number of variables, mi0 is the optimum value of mi, t is the minimum value of φ(m), and the value of the posterior probability tends to be optimal. t generally takes the average of the minimum values of the first 50 cost functions in the optimization process. The one-dimensional marginal probability densities of parameters can be written as
| $$ Pro({{\boldsymbol{m}}_i}|{\boldsymbol{d}}) = \int { \text{δ} ({{\boldsymbol{m}}_i} - {{\boldsymbol{m}}_i}^\prime )Pro({\boldsymbol{m'}}|{\boldsymbol{d}})} {\rm{d}}{\boldsymbol{m'}} . $$ | (33) |
The δ function is a unit impulse function, which is equal to 0 at all points except 0, and its integral over the entire domain is equal to 1. In this inversion scheme, the maximum approximation probability (MAP) is used to replace the posterior mean as the estimation result of the model parameters in this paper. When the MAP estimation reaches the maximum value, there is
| $$ {{\hat {\boldsymbol{m}}}} = {\rm{Ar}}{{\rm{g}}_{\max}}\{ Pro({\boldsymbol{m}}\left| {\boldsymbol{d}} \right.)\} . $$ | (34) |
On the base of Eq. (30), Eq. (34) can be written as
| $$ {{\hat {\boldsymbol{m}}}} = {\rm{Ar}}{{\rm{g}}_{\min }}\{ \varphi ({\boldsymbol{m}})\} . $$ | (35) |
In Eq. (35),
In the Bayesian inversion theory, calculating the replica is important. As stated above, the replica can be calculated with using Eq. (29). In this way, the inversion parameters are: seabed density ρb, seabed sound velocity cb, propagation distance and sea depth, water depth H and propagation range r. Only the average sound velocity in water
In this method, the replica does not involve information such as the frequency spectrum of the sound source, therefore, sound source information need not to be estimated invisibly in the cost functions, thus, greatly improving the accuracy. What’s more, only five parameters are required to calculate the replica, and the calculation speed and inversion efficiency are substantially improved. The cost function is established as follows.
The curve with frequency f>fairy (fairy is the airy frequency, and the group speed is the minimum value at fairy) has a higher amount of energy and is easy to extract, in agreement with the dispersion results in the TF domain, moreover, the extracted value is effective. To obtain better inversion results, the extracted dispersion curves from data on the frequency band f>fairy, Fn, is used as the input function in Bayesian inversion theory. The dispersion curves
Under the assumption that the data errors satisfy a Gaussian distribution, the likelihood function is
| $$ L({\boldsymbol{m}},r) = \prod\limits_{n = 1}^N {\frac{1}{{{{(2\pi )}^{{M_n}/2}}{{\left| {{{\boldsymbol{C}}_n}} \right|}^{1/2}}}}\exp \left[ - \frac{1}{2}{{({{\boldsymbol{F}}_n} - {\boldsymbol{F}}_n^c)}^{\rm{T}}}{\boldsymbol{C}}_n^{ - 1}({{\boldsymbol{F}}_n} - {\boldsymbol{F}}_n^c)\right] } . $$ | (36) |
There is a covariance matrix
| $$ L({\boldsymbol{m}},r) = \prod\limits_{n = 1}^N {\frac{1}{{{{(2\pi {\upsilon _k})}^{{M_n}/2}}}}\exp \left( - \frac{1}{{2{\upsilon _k}}}{{\left| {{{\boldsymbol{F}}_n} - {\boldsymbol{F}}_n^c} \right|}^2}\right)} . $$ | (37) |
The
| $$ \frac{{\partial L}}{{\partial \upsilon _n^{}}} = 0 . $$ | (38) |
Therefore
| $$ \hat \upsilon _n^{} = {\left| {{{\boldsymbol{F}}_n} - {\boldsymbol{F}}_n^c} \right|^2}/{M_n} . $$ | (39) |
According to the relationship between the likelihood function and the cost function,
| $$ L({\boldsymbol{m}},r) = \exp [- \varphi ({\boldsymbol{m}})] . $$ | (40) |
Therefore, the cost function is
| $$ \varphi ({\boldsymbol{m}}) = \frac{1}{2}\sum\limits_{n = 1}^N {{M_n}\ln } ({\left| {{{\boldsymbol{F}}_n} - {\boldsymbol{F}}_n^c} \right|^2}) . $$ | (41) |
The parameters are then optimized, and the obtained optimal values are substituted into Eq. (31) to solve the posterior probability distribution of each inversion parameter. The Genetic Algorithm (GA) (Whitley, 1994) is used to calculate the optimal solution. When the search result finally tends to a stable and convergent value, this value is the optimal value. The initial population size, generation number, crossover probability, selection probability, and mutation probability are 1 000, 3 000, 0.8, 0.5, and 0.05, respectively, and a global optimum will be ensured by 30 sets of programs.
The sensitivity of all inversion parameters to the cost function in Eq. (41) is analyzed. The sound speed is a parameter that can be easily obtained through experience. The analysis process is performed in a shallow water waveguide. The simulation parameters are shown in Fig. 3, and the depth is H=100 m. The average sound speed is c0=1 531.5 m/s. A chirp sound source with a bandwidth of 50–200 Hz is emitted at a depth of zs=90 m, the hydrophone is installed at depth zr=100 m, the propagation range r=10 km, the signal-to-noise ratio is 15 dB, the seabed sound speed cb=1 700 m/s, the seabed density ρb=1.6 g/cm3 and the bottom reflection phase shift parameter P=7.355. The replica is calculated by the Eq. (29).
When other parameters are the true values and the analyzed parameter changes within a certain range, the analysis results are as seen in Fig. 6, it shows that, near the reference value, the cost function curve becomes very sharp and gradually reaches the minimum value. Therefore, the cost function is valid for all inversion parameters. What is more, compared with the cb and ρb, the cost function change most dramatically near to the true value of r and H, especially the range r, which shows that the cost function is most sensitive to r and H, indicating that the cost function is more sensitive to them. Therefore, when the cost function Eq. (41) is used to invert the parameters, the probability of obtaining an accurate r value is higher.
In order to obtain accurate inversion results, it is necessary to analyze the correlation among the inversion parameters. Correlation results are used to analyze whether there is coupling between parameters, and then the accuracy of inversion results is estimated. The two-dimensional sensitivity analysis can be used to study the parameter inter-relationship. In the process of two-dimensional sensitivity analysis, the two parameters to be analyzed will change within the search interval, while the other parameters are true values. The location where the cost function theoretically minimizes (both parameters to be analyzed are true value) is indicated by a white asterisk. After the cost function value is normalized, the analysis results are shown in Fig. 7.
The coupling between the inversion parameters will lead to errors in the inversion results, and reducing the inversion parameters is an effective means to solve this problem. As mentioned above, seabed characteristics can be described by bottom reflection phase shift parameter P, so cb and ρb can be combined into one parameter P. when the replica is calculated by the Eq. (28), the number of inversion parameter will be resuced to 3, the two-dimension sensitivity analysis result in this case is shown in Fig. 8.
It can be seen that the propagation range r is in dependent of other parameters in Figs 8a and b. Figure 8c shows that the minimum value of the cost function changes with the P and H approximately as a line, but the line can be approximately equal to the true value of H, indicating that the H is more sensitive to the cost function than the inversion of four parameters. In order to reduce the the influence of parameter correlation and improve the accuracy of inversion, both Eqs (28) and (29) are used to the inversion process. Equation (28) is first used to calculate the replica, the seabed characteristics is decribed by the P. In the case of three inversion parameters, the propagation range r and the depth H are estimated. Then , the r and H will be submitted to Eq. (29) to calculate the replica again as a priori parameters, the seabed parameters cb and ρb can be inversed. The flowchart of parameter inversion is shown in Fig. 9.
The inversion parameters of the waveguide in Section 3.3 were searched over a relatively wide interval by GA, and the optimization interval of the parameters is shown in Table 1. When the optimal values of parameters are searched, the PPD and MPD of the parameters can be estimated by Eqs (31) and (33).
| Parameter | Reference value | Search bound | Inversion value | Error/% |
| r/km | 10 | [5, 15] | 9.99 | 0.1 |
| −cb/(m·s−1) | −1 700 | −[1 600, 2 000] | −1 690.83 | −0.54 |
| $\rho $b/(g·cm−3) | 1.6 | [1.2, 2] | 1.56 | 2.5 |
| H/m | 100 | [50, 180] | 100.78 | 0.78 |
DownLoad:
CSV
Depending on the flowchart of parameter inversion shown in Fig. 9, the range r and the depth H are inversed firstly, the seabed parameters cb and ρb are estimated next. The variation of inversion parameters with the iterations number of 30 sets is shown in Figs 10 and 11. It can be seen that the convergence rate of the three inversion parameters is fast. The propagation range is the first to converge to the true value. The convergence rate of the parameter is often in direct proportion to its sensitivity to the cost function. What’s more, fast convergence rate also shows the high efficiency of the inversion method. When GA is used for optimization, the optimization time of 5 000 generations is less than 3 min.
The one-dimensional marginal probability of inversion parameters is shown in Fig. 11. The maximum value of the one-dimensional marginal probability densities corresponds to the final inversion result of the parameter.
It can be seen from Fig. 11 and Table 1 that the convergence probability of propagation r near the true value is high, the global optimal value is 9.99 km, the error percentage of which is 0.1%, the global optimal value of cb is 1 690.83 m/s, the error percentage of which being 0.54%, and the global optimal value of H is 100.31 m, the error percentage of which is 0.78%. The global optimal value of ρb is 2.5%. Figure 7d shows that there is a sub-maximum near the global optimal value, which is related to the correlation between ρb and cb. In order to avoid this problem, Hamilton empirical formulas
The inversion errors of all parameters are within 3%, which proves the accuracy and effectiveness of the method in this paper. The inverted parameters were substituted into Eq. (29) to calculate the dispersion curves, and the results were compared with the extracted results from the data. Figure 12 shows that the results were in good agreement, which proves that the estimated values of the three inversion parameters are accurate and effective.
Through the above simulation results, these conclusions can be obtained: (1) the inversion method is suitable for the SRBR normal modes in shallow water waveguide, and (2) the replica calculated by the dispersion formula is simple, effective, and fast.
The dispersion formulas shown in Eqs (28) and (29) contain several parameters, but the average sound speed in water is not used as an inversion parameter for the following two reasons: (1) the sound speed is a parameter that can be easily obtained through experience, (2) fewer inversion parameters can avoid errors caused by correlation and improve the inversion efficiency. However, inaccurate sound speed estimation will affect the entire inversion process. Now the influence of sound speed mismatch on inversion results is analyzed.
The simulation environment is shown in Fig. 3, the average sound speed is c0=1 531.5 m/s. when the sound speed in the range of 1 500 m/s to 1 560 m/s, the errors of the three inversion parameters are analyzed under different sound velocity values. According to the inversion process described above, with 5 m/s as the interval, the final inversion results are shown in Fig. 13.
It can be seen from Fig. 13 that the relative errors of r, H are both less than 3%, which indicates that the inversion parameters are robust when the average sound speed changes within the range of ±30 m/s. The relative errors of cb and ρb are within 10%. In addition to the coupling of these two seabed parameters, the reason for the larger error is that these two parameters and the speed in water have a greater correlation in describing the dispersion characteristics of the waveguide. However, when the average sound speed in water changes within the range of ±10 m/s, the error of each parameter can be controlled within 5%. Therefore, it is feasible to use the average speed of sound in seawater to calculate the replica.
The data is obtained by an experiment performed in an area of the Yellow Sea in China. Experimental sea area is a shallow water waveguide with a liquid bottom, the schematic diagram of the experimental device is shown in Fig. 14. The explosion sound source is emitted, the water depth is H=30 m, the sound velocity is approximately equal to isovelocity with c0=1 461 m/s (the sound speed profile is shown in Fig. 14). The ocean bottom seismometer (OBS) was placed on the bottom. The OBS is a vector hydrophone with four signal receiving channels, X, Y, Z and P. The signal processing frequency band is 100–300 Hz, the propagation range measured by a global positioning system (GPS) is r=29.52 km.
The received signal in time domain of four channels is shown in Fig. 15. The STFT is used to preprocess the signals, the spectrum of the received signal after STFT in TF domain are shown in Fig. 16. It can be seen from Fig. 16 that the received signals in X, Y and P contain three normal modes with higher energy. The energy concentration of the signal in the Z channel is poor, but three modes can also be obtained after warping. Therefore, all four sets of data are used for inversion analysis. Then the dispersion curve extraction method mentioned above is used to process the data. The extracted curves (red dotted line) are also shown in Fig. 16. The extracted dispersion curve is consistent with the energy change curve of each normal mode in the time-frequency domain. The extracted curves are used as an input of inversion methodology. Different mode has different cutoff frequencies (the lowest frequency that can excite mode), so the dispersion curves of modes have different frequency range. Different modes will use different frequency band for inversion. The inversion results of the four sets of data from the same waveguide should be consistent. The four sets of data are used for parameters inversion, and then the inversion results are compared and analyzed.
According to the inversion process in Fig. 9, the propagation range r, the depth H, the seabed cb and the seabed density ρb will be inversed. The inversion values are shown in Table 2. The inverted results were compared with the reference value. According to the inversion results, the results inversed from the four sets of data are similar. The dispersion curves are calculated by substituting four sets of inversion parameters into Eq. (29), and compared with the spectrum of received signal in TF domain and the extracted curves, the results are shown in Fig. 17. In each channel, the dispersion curves obtained from the inversion parameters are in good agreement with the data extraction results. Dispersion is a characteristic of the waveguide and has nothing to do with the receiving channel, so the final inversion value is the average of the four sets of results.
| Parameter | Reference value | Search bound | Inversion value | ||||
| X | Y | Z | P | Mean value | |||
| r/km | 29.52 | [10, 20] | 29.21 | 29.21 | 29.22 | 29.21 | 29.21 |
| cb/(m·s−1) | − | [1 500, 2 000] | 1 561.43 | 1 607.13 | 1 627.22 | 1 564.45 | 1 590.06 |
| $\rho $b/(g·cm−3) | − | [1.2, 2] | 1.55 | 1.56 | 1.62 | 1.52 | 1.55 |
| H/m | 30 | [20, 100] | 29.02 | 30.14 | 31.65 | 29.71 | 30.11 |
| Note: − represents no data. | |||||||
DownLoad:
CSV
Compared with the GPS measurement result obtained during the experiment (29.52 km), the inversion propagation range r=29.21 km is a reliable value. Furthermore, the inversion depth is H=30.11 m. The errors of r and H were both less than 1%. The inversion results of cb and ρb are 1 590.6 m/s and 1.55 g/cm3, respectively, and the accuracy of values can be verified by the dispersion curves calculated from the inversion parameters. When the calculated result of the dispersion curve is consistent with the extracted result, the accuracy of the inversion result can be proved. The comparison result of the dispersion curves calculated by inversion parameters of four sets and the dispersion curves calculated by the final inversion value is shown in Fig. 18. After comparison and analysis, it is found that all the dispersion curves are in good agreement, which proves the validity and accuracy of the inversion results.
Bayesian inversion theory is a good way to inverse marine environmental parameters when there is little prior knowledge. A fast inversion method of marine environment parameters based on dispersion characteristics with single hydrophone is presented in this paper. This method is suitable for broadband pulsed sound source signals that propagate over long distances in shallow water waveguides. The normal modes in received signal could be described by dispersion curves in the TF domain, which can be extracted by a signal processing method known as the warping operator. In this paper, the curves are used as the input of inversion methodology. In this method, there are many advantages of acquiring input functions, only a single hydrophone is used, multiple groups are acquired at the same time, the extraction process is simple and easy to operate, and the cost is low.
For Bayesian methodology, the biggest advantage is to realize the inversion for ocean parameters based on less prior knowledge. Sensitive cost function and efficient replica are the keys to accurate inversion. A new dispersion formula is used to calculate the replica in this paper. The dispersion formula is composed by four parameters, which are all sensitive to cost function.
In shallow water, as for SRBR modes, the bottom characteristics can be described by the bottom reflection phase shift parameter P when the sound source is incidental at a small glancing angle based on the BDRM theory. Therefore, the number of inversion parameters was reduced because of the application of P, so the propagation range and the water depth can be inversed quickly and accurately. Then the inversion results can be used ton inverse other seabed parameters as a prior knowledge, which improves the inversion accuracy and efficiency again.
What is more, compared with the sound field model, the advantages of using a dispersion formula to calculate the replica include requiring less prior knowledge and a faster calculation speed.
The error of the inversion results of the simulation and experimental data is within 5%, which is in good agreement with the real values. The validity and accuracy of the inversion method have been proved. However, the applicability of this method in complex ocean waveguides needs further research and improvement.
| [1] |
Bao Qingliu, Zhang Youguang, Lin Mingsen, et al. 2017. An ocean current inversion accuracy analysis based on a Doppler spectrum model. Acta Oceanologica Sinica, 36(9): 101–107. doi: 10.1007/s13131-017-1115-y
|
| [2] |
Bonnel J, Dosso S E, Chapma N R. 2013. Bayesian geoacoustic inversion of single hydrophone light bulb data using warping dispersion analysis. The Journal of the Acoustical Society of America, 134(1): 120–130. doi: 10.1121/1.4809678
|
| [3] |
Bonnel J, Gervaise C, Nicolas B, et al. 2012. Single-receiver geoacoustic inversion using modal reversal. The Journal of the Acoustical Society of America, 131(1): 119–128. doi: 10.1121/1.3664083
|
| [4] |
Bonnel J, Nicolas B, Mars J I, et al. 2010. Estimation of modal group velocities with a single receiver for geoacoustic inversion in shallow water. The Journal of the Acoustical Society of America, 128(2): 719–727. doi: 10.1121/1.3459855
|
| [5] |
Bonnel J, Thode A, Wright D, et al. 2020. Nonlinear time-warping made simple: a step-by-step tutorial on underwater acoustic modal separation with a single hydrophone. The Journal of the Acoustical Society of America, 147(3): 1897–1926. doi: 10.1121/10.0000937
|
| [6] |
Cai Haiyan, Jiang Qingtang, Li Lin, et al. 2021. Analysis of adaptive short-time Fourier transform-based synchrosqueezing transform. Analysis and Applications, 19(1): 71–105. doi: 10.1142/S0219530520400047
|
| [7] |
Dosso S E. 2002. Quantifying uncertainty in geoacoustic inversion. I. A fast Gibbs sampler approach. The Journal of the Acoustical Society of America, 111(1): 129–142. doi: 10.1121/1.1419086
|
| [8] |
Dosso S E, Wilmut M J. 2008. Uncertainty estimation in simultaneous Bayesian tracking and environmental inversion. The Journal of the Acoustical Society of America, 124(1): 82–97. doi: 10.1121/1.2918244
|
| [9] |
Dosso S E, Wilmut M J. 2011. Bayesian multiple-source localization in an uncertain ocean environment. The Journal of the Acoustical Society of America, 129(6): 3577–3589. doi: 10.1121/1.3575594
|
| [10] |
Fallat M R, Dosso S E. 1999. Geoacoustic inversion via local, global, and hybrid algorithms. The Journal of the Acoustical Society of America, 105(6): 3219–3230. doi: 10.1121/1.424651
|
| [11] |
Gingras D F, Gerstoft P. 1995. Inversion for geometric and geoacoustic parameters in shallow water: experimental results. The Journal of the Acoustical Society of America, 97(6): 3589–3598. doi: 10.1121/1.412442
|
| [12] |
Heaney K D. 2004. Rapid geoacoustic characterization using a surface ship of opportunity. IEEE Journal of Oceanic Engineering, 29(1): 88–99. doi: 10.1109/JOE.2003.823286
|
| [13] |
Jensen F B, Kuperman W A, Porter M B, et al. 2011. Computational Ocean Acoustics. 2nd ed. New York: Springer, 354–356
|
| [14] |
Le Gac J C, Asch M, Stephan Y, et al. 2003. Geoacoustic inversion of broad-band acoustic data in shallow water on a single hydrophone. IEEE Journal of Oceanic Engineering, 28(3): 479–493. doi: 10.1109/JOE.2003.816689
|
| [15] |
Le Touze G, Nicolas B, Mars J I, et al. 2009. Matched representations and filters for guided waves. IEEE Transactions on Signal Processing, 57(5): 1783–1795. doi: 10.1109/TSP.2009.2013907
|
| [16] |
Li Xiaoman, Piao Shengchun, Zhang Minghui, et al. 2019a. A passive source location method in a shallow water waveguide with a single sensor based on Bayesian theory. Sensors (Basel), 19(6): 1452. doi: 10.3390/s19061452
|
| [17] |
Li Qianqian, Shi Juan, Li Zhenglin, et al. 2019b. Acoustic sound speed profile inversion based on orthogonal matching pursuit. Acta Oceanologica Sinica, 38(11): 149–157. doi: 10.1007/s13131-019-1505-4
|
| [18] |
Li Qianqian, Yang Fanlin, Zhang Kai, et al. 2016. Moving source parameter estimation in an uncertain environment. Acta Physica Sinica, 65(16): 164304. doi: 10.7498/aps.65.164304
|
| [19] |
Li Qianqian, Yang Fanlin, Zhang Kai. 2018. Multiple source localization using Bayesian theory in an uncertain environment. Haiyang Xuebao (in Chinese), 40(1): 39–46
|
| [20] |
Niu Haiqiang, Zhang Renhe, Li Zhenglin. 2014. Theoretical analysis of warping operators for non-ideal shallow water waveguides. The Journal of the Acoustical Society of America, 136(1): 53–65. doi: 10.1121/1.4883370
|
| [21] |
Porter M B. 1992. The KRAKEN normal mode program. Washington, DC: Naval Research Lab
|
| [22] |
Shang E C, Wu J R, Zhao Z D. 2012. Relating waveguide invariant and bottom reflection phase-shift parameter P in a Pekeris waveguide. The Journal of the Acoustical Society of America, 131(5): 3691–3697. doi: 10.1121/1.3699242
|
| [23] |
Walker S C, Roux P, Kuperman W A. 2005. Data-based mode extraction with a partial water column spanning array. The Journal of the Acoustical Society of America, 118(3): 1518–1525. doi: 10.1121/1.1993149
|
| [24] |
Wang Dong, Guo Lianghao, Liu Jianjun, et al. 2016. Passive impulsive source range estimation based on warping operator in shallow water. Acta Physica Sinica, 65(10): 104302. doi: 10.7498/aps.65.104302
|
| [25] |
Wang Dezhao, Shang Erchang. 2013. Underwater Acoustics. 2rd ed. Beijing: Science Press, 158–164
|
| [26] |
Whitley D. 1994. A genetic algorithm tutorial. Statistics and Computing, 4(2): 65–85
|
| 1. | Xiaoman Li, Hongyun Chen, Hongyu Lu, et al. A method of underwater sound source range estimation without prior knowledge based on single sensor in shallow water. Frontiers in Physics, 2023, 11 doi:10.3389/fphy.2023.1109220 | |
| 2. | Mian Wei, Shuqiao Wang, Yucheng Lin. Systematic notes on three new Luthela (Mesothelae, Heptathelidae) spiders from China, with their descriptions. ZooKeys, 2023, 1159: 151. doi:10.3897/zookeys.1159.90120 | |
| 3. | Varat Sivayyapram, Chawakorn Kunsete, Xin Xu, et al. Molecular phylogeny, biogeography, and species delimitation of segmented spider genus Liphistius (Araneae: Liphistiidae) in Thailand. Zoological Journal of the Linnean Society, 2023. doi:10.1093/zoolinnean/zlad149 | |
| 4. | Yan Zhang, Zhaoyang Chen, Daiqin Li, et al. Three new species of the primitively segmented spider genus Songthela (Mesothelae, Liphistiidae, Heptathelinae) from Hunan Province, China. ZooKeys, 2023, 1154: 17. doi:10.3897/zookeys.1154.98273 |
Figures(18) / Tables(2)
Supported by:
Beijing Renhe Information Technology Co. Ltd
Xiaoman Li, Biao Wang, Xuejie Bi, Hong Wu. A fast inversion method for ocean parameters based on dispersion curves with a single hydrophone[J]. Acta Oceanologica Sinica, 2022, 41(9): 71-85. doi: 10.1007/s13131-022-1999-z
| Parameter | Reference value | Search bound | Inversion value | Error/% |
| r/km | 10 | [5, 15] | 9.99 | 0.1 |
| −cb/(m·s−1) | −1 700 | −[1 600, 2 000] | −1 690.83 | −0.54 |
| $\rho $b/(g·cm−3) | 1.6 | [1.2, 2] | 1.56 | 2.5 |
| H/m | 100 | [50, 180] | 100.78 | 0.78 |
DownLoad:
CSV
| Parameter | Reference value | Search bound | Inversion value | ||||
| X | Y | Z | P | Mean value | |||
| r/km | 29.52 | [10, 20] | 29.21 | 29.21 | 29.22 | 29.21 | 29.21 |
| cb/(m·s−1) | − | [1 500, 2 000] | 1 561.43 | 1 607.13 | 1 627.22 | 1 564.45 | 1 590.06 |
| $\rho $b/(g·cm−3) | − | [1.2, 2] | 1.55 | 1.56 | 1.62 | 1.52 | 1.55 |
| H/m | 30 | [20, 100] | 29.02 | 30.14 | 31.65 | 29.71 | 30.11 |
| Note: − represents no data. | |||||||
DownLoad:
CSV
| Parameter | Reference value | Search bound | Inversion value | Error/% |
| r/km | 10 | [5, 15] | 9.99 | 0.1 |
| −cb/(m·s−1) | −1 700 | −[1 600, 2 000] | −1 690.83 | −0.54 |
| $\rho $b/(g·cm−3) | 1.6 | [1.2, 2] | 1.56 | 2.5 |
| H/m | 100 | [50, 180] | 100.78 | 0.78 |
| Parameter | Reference value | Search bound | Inversion value | ||||
| X | Y | Z | P | Mean value | |||
| r/km | 29.52 | [10, 20] | 29.21 | 29.21 | 29.22 | 29.21 | 29.21 |
| cb/(m·s−1) | − | [1 500, 2 000] | 1 561.43 | 1 607.13 | 1 627.22 | 1 564.45 | 1 590.06 |
| $\rho $b/(g·cm−3) | − | [1.2, 2] | 1.55 | 1.56 | 1.62 | 1.52 | 1.55 |
| H/m | 30 | [20, 100] | 29.02 | 30.14 | 31.65 | 29.71 | 30.11 |
| Note: − represents no data. | |||||||
DownLoad:
DownLoad:
