In Beam-Displacement Ray-Mode (BDRM) theory, when a sound source with a small grazing angle propagates in a waveguide where the speed of sound of the seabed is larger than that in water, the dispersion curve of the SRBR normal mode can also be calculated with the dispersion formula [25]: f n t ≈ n π c 0 t 2 π H t 2 − r / c 0 2 − P c 0 4 π H − Pr 4 π H t 2 c 0 (5) where the P is the bottom reflection phase shift parameter, which is a quantity associated with the seabed parameter, expressed as: P ≈ 2 ρ 2 / ρ 1 1 − c 0 / c b 2 2 + c 0 / c b 2 α c b 2 π f + 1 − c 0 / c b 2 1 − c 0 / c b 2 2 + c 0 / c b 2 α c b 2 π f (6) where ρ ₁, ρ ₂, c _b, and α are the water density, seabed density, and seabed absorption, respectively. According to Eq. 6, the P is a frequency-dependent quantity. However, when the attenuation of the seabed is minor, and the frequency exceeds 10 Hz, the effect of frequency on P can be neglected, and Eq. 6 can be simplified to: P ≈ 2 ρ 2 ρ 1 1 − c 0 / c b 2 (7)

The inversion efficiency can be improved by reducing the number of inversion parameters by describing the bottom properties with parameter P. Moreover, the dispersion curve can be calculated more rapidly by using the dispersion formula compared with the sound field model (e.g., KRAKENC) [26].

According to the above, the dispersion curves of each order of simple positive waves can be calculated with Eq.5 and 6, when the parameters of the marine environment are known. KRAKENC is a more mature model for acoustic field calculation; therefore, the dispersion curves obtained by KRAKENC according to Eq. 5 are used as the reference curves.

However, the dispersion curves obtained with the dispersion formula Eq. 5 do not contain frequencies lower than Ariy frequency (Ariy frequency: The frequency corresponding to the minimum value of the group velocity of the normal mode). It is also an important part of our future work to obtain dispersion formulas that can calculate complete dispersion curves.

In practical applications, detailed environmental parameters are difficult to obtain, and the dispersion curves can be extracted through analysis of the received signals. The combination of time-frequency analysis techniques, warping transform, modal filtering, and ridge extraction techniques can achieve accurate extraction of dispersion curves in measured data.

The warping transform can transform the time domain signal pr e t r , t (the time domain result of pr e f f , r ) by a time transformation factor h t = t 2 + r / c 0 2 , thereby transforming the sound pressure signal into a superposition of multiple single frequency signals to achieve effective separation of different normal modes. The modal filtering technique can then be used to extract normal modes.

The extracted modes are inverted with the inverse warping operator h t = t 2 − r / c 0 2 , and finally the dispersion curve of each mode is extracted with techniques such as ridge extraction. In the above process, two parameters are used: the propagation distance r and the speed of sound in seawater c0. Because the warping transform is invertible, the two parameters are process quantities, and the final results are not dependent on the accuracy of the two parameters. However, the arrival time of each normal mode must be accurate, and great care must be taken in processing the received signal [27].

The above process can be described by the simulation results in the Pekeris waveguide, with the same simulation environment as in the previous section, taking the example of extracting the dispersion curve of normal mode x3.

In the spectrogram of the received signal shown in Figure 3A, the energy of mode ×3 is strong, thus facilitating its extraction. The received signal contains at least 5th-order normal modes, but the energy of different-order modes varies, and the energy of mode four is weak and appears blurry in Figure 3B. Mode x3 is extracted after modal filtering and unwarping transform, and the spectrogram and extracted dispersion curve of modal x3 are shown in Figures 3C, D, respectively.

The dispersion formula Eq. 5, data-extraction and theoretical calculation based on an acoustic model are three methods for calculating the dispersion curve. The accuracy of the calculated results is verified by comparing with the results from KRAKEN, used as reference value. The simulation analysis is performed in the Pekeris waveguide and the waveguide with a thermocline, respectively, whose simulation parameters are shown in Figure 4. The propagation range of source is r = 5 km, and the broadband of source is 100–200 Hz, the sampling rate is f _s = 1 kHz, the spectrum of the received signal is obtained by the short-time Fourier transform (STFT). The spectrum information is illustrated in Figure 3, where the information of RGB is given. The locations of the source and receiver, and the seabed environment are the same as in the Pekeris waveguide, but the speed of sound in water has a thermocline at 20–40 m. The difference between the speed of sound at the surface and the seabed is 9 m/s, and the source is located in the thermocline. Figure 5A and Figure 5B show the TFA results of the received signal in the two waveguides and the dispersion curves calculated with the three methods, respectively.

Figure 5 shows that the dispersion curves obtained through the three methods match well in the Pekeris waveguide, and are consistent with the mode energy change trend. In a thermocline waveguide, the dispersion curves extracted from data and calculated with Eq. 5 are consistent, but the low order normal modes do not fully conform to the characteristics of SRBR normal modes, and the consistency is poor. The results of the extracted dispersion curves, particularly the high energy modes, represent the energy variation trend of each mode and can be used in underwater applications.

In addition, Figure 5A shows that the mode order corresponding to the higher energy modes (modes x1, x2, and x3) in the received signal are 2, 3, and 5. The properties of the dispersion curves correspond to the mode order, and the accurate identification of the extracted modes is key to the passive ranging of sources by using the dispersion curves in the case of unknown marine environmental parameters. However, the identification of the modes in the case of unknown a priori knowledge presents the following difficulties.

1) Because of the weak intensity of the first-order mode, the first mode that can be detected in the TFA results of the received signal is not necessarily the true first-order normal mode.

The posterior probability density (PPD) of the parameters to be inverted is used as the estimation result, in an important feature of Bayesian theory. The PPD is associated with the conditional probability density of measured data P r o d m and the prior distribution of m and d, P r o m and P r o d , and the PPD can be written as [28]: P r o m d = P r o d m P r o m / P r o d (8) where m is a vector of inversion parameters with length M, and d is the parameter vector of measurement data. In general, P r o m d a is positively correlated with the likelihood function L m , so P r o m d is: P r o m d ∝ L m P r o m (9)

The likelihood function L m can be represented by the misfit function E m , as shown in the following equation: L m = P r o d m ∝ exp − E m (10)

Furthermore, where E(m) is the data matching function, a generalized misfit combining data and prior can be defined as φ m = E m − ln P r o m (11)

The φ(m) is also known as the cost function. According to Eq.10 and Eq.11, the P r o m d can also be expressed as P r o m d = exp − φ m ∫ exp − φ m ′ , d d m ′ (12)

To decrease the difficulty in calculating the integral in the above equation, the principle for parameter estimation is the maximum a posteriori (MAP). The inversion parameter values can initially be determined through the optimization search algorithm. The optimal value of a set of inversion parameters is obtained when the maximum likelihood function is maximal. m ^ = A r g min L m = A r g min φ m (13)

The uncertainty of the inversion parameters is characterized by the one-dimensional marginal probability distribution of the parameters P r o m i d = ∫ δ m i ′ − m i P r o m ′ d d m ′ (14) Where m _i refers to the ith inversion parameter。

As described above, the cost function can be derived from a likelihood function, which is key to achieving parameter estimation. The input function and the replica are two important components of the likelihood function. The dispersion curve F n m extracted from the data is taken as the input function, and the replica F n c is calculated with Eq. 5. The results in Figure 5 confirm the feasibility of calculating the replica with Eq. 5. If the data error satisfies a Gaussian distribution, the likelihood function is: L m , r = ∏ n = 1 N exp − 1 2 F n m − F n r c T C n − 1 F n m − F n r c 2 π M n / 2 C n 1 / 2 (15) where ∑ n = 1 n = N is the number of extracted dispersion curves, does not represent the exact order n of normal mode, and M n is the length of the data. C n = υ n I is the covariance matrix, and υ n and I are the unknown variance and the identity matrix, respectively. Therefore, the likelihood function is: L m , r = ∏ n = 1 N 1 2 π υ n 1 / 2 exp − 1 2 υ n F n m − F n r c 2 (16) let ∂ L ∂ υ n = 0 , then υ n is υ ^ n = F n m − F n r c 2 / M n (17) the cost function is: φ m , n r = 1 2 ∑ n = 1 N M n ln F n − F n r c 2 (18)

In the absence of prior knowledge of the environment, the true order of mode n must be determined simultaneously, on the basis of the assumption that the true order is nr Therefore, each set of data also must satisfy: φ n n r = 1 2 ∑ n x = 1 N x M n ln F n m − F n x c 2 (19) where n x ∈ 1 , N x is the range of normal mode numbers, and N _x is sufficiently large to exceed the maximum order of modes obtained from the received data, thus ensuring the accuracy of the numerical determination.

Eq. 18 and Eq. 19 affect each other, and the order determination takes precedence over the parameter inversion in the optimization search process. After the order of normal mode numbers of the nth data set is determined by Eq. 19, the nth data set is then used to perform the parameter search and inversion according to Eq. 18. The optimal solution of the parameters is substituted into Eq. 12 to obtain the PPD of each inversion parameter, and then the effective inversion of the parameters is achieved. The Genetic Algorithm is used as a method of parameter search for optimization.

According to Eq. 5, the calculation of the replica in an unknown ocean environment requires the depth H, average speed of sound c₀, bottom reflection phase shift parameter P, and range r. Herein, the above four parameters are estimated jointly with the mode order.

The sensitivity of the four parameters to be inverted is analyzed in the waveguide with a thermocline. The simulation environment is shown in Figure 4. To fully verify the validity of the method for the mode order identification, the depth of the receiver is at the node of mode 5, which is also near the node of mode 2; i.e., theoretically the received signal has low energy of modes 5 and 2.

The input function is the dispersion curves extracted from the data, and the rest of the parameters are taken as true values when sensitivity analysis is performed for one parameter. The analysis results are shown in Figures 6A–E.

Figure 6 shows that the four inversion parameters are sensitive to the cost function. Compared with other parameters, the value cost function curve of r varies more sharply around the true value, thus fully demonstrating the sensitivity of the range to the cost function and the accuracy of inversion. The dispersion curve data of mode 5 from the received signal are selected to analyze the sensitivity of the normal mode order n. Figure 6E indicates that the cost function has the smallest value at the true order of the mode, thus indicating the validity of the cost function.

In addition, for a more adequate analysis of the validity of the cost function, two-dimensional correlations of inversion parameters are analyzed. The two parameters to be analyzed vary within the search interval, and the remaining parameters are taken as true values. The true values of the two parameters are marked with a white star. Figure 7 shows that the correlation between the range r and other parameters is low, with a nearly straight line at r ₀, thus indicating that r plays a major role in the influence of the cost function with high sensitivity. However, the correlation between P and H is high, the correlation between c ₀ and H gradually decreases with depth, and the speed of sound gradually dominates. Parameter coupling caused by correlation between parameters affects the accuracy of the inversion results. However, as shown in Figure 7, the parameter coupling does not affect the sensitivity of r. Comprehensively, the four inversion parameters and mode order have high sensitivity to the joint cost function.

The depth of the experimental sea H = 24.5 m, and the SSP in the water has a thermocline, as shown in Figure 11, the average speed of sound is c0 = 1497 m/s, the broadband source with the broadband 200–500 Hz is located at z_s = 10 m, and p = 6.33. The signal is received by a single hydrophone at z_r = 9 m, and the receiving depth is near the node of mode 3. r ₀ = 4.96 km, as measured by GPS. The source is a long pulse, and the received signal is processed by pulse compression for the next step of analysis.

The received signal is shown in Figure 12A, and the signal to noise ratio is high. The dispersion curves are extracted after warping transform. Figure 12B shows the time-frequency analysis result of the received signal. The normal modes with stronger energy are selected for dispersion curve extraction, and the extraction results are shown in Figure 12B. The received signal contains four orders modes with strong energy, and the position of mode 1 is clear, but several modes are missing after mode 2, and the number of the missing modes cannot be known, according to the figure. The part of the dispersion curve with high energy and strong aggregation is used as the input function. To ensure accuracy of the inversion, the frequency range matching the different orders of the simple positive waves is selected, and the final dispersion curve as the input function is shown in Figure12B. The search intervals of the parameters to be inverted and the reference values are shown in Table 2.

Parameter inversion and modal identification are performed according to the procedure described above, and the inversion results are shown in Table 2. As shown in Table 2 and Figure 13, the results of normal mode number determination are accurate, and the results indicate that the missing mode is mode 3. Because of the location of the receiver, the energy of mode 3 is low and it is unobserved in the time-frequency domain of the received signal.

The errors of c0 and H are less than 5%, on the basis of comparison of the inversion results with the reference values, which are 2.04% and 0.002%, respectively. Although the error of P is larger, the estimated value of the propagation range is r = 5.00 km, which is a very accurate value with respect to the GPS result of r ₀ = 4.96 km, and the error is only 0.008%. The inversion results fully illustrate that the method has strong stability with respect to marine environmental parameters. The inversed dispersion curves match well with the data extraction results and the spectrogram of the received data, thus demonstrating the superiority of the proposed method in passive ranging and normal mode order determination.

The depth of the experiment sea is H = 30 m, the SSP is isovelocity with the speed c0 = 1460 m/s, and p = 5.43. An explosion sound signal was received by the ocean bottom seismometer located at the seabed after propagating 35.12 km, and four channels (P, X, Y, and Z) of data were obtained. The P-channel data were selected for inversion, and the received signal is shown in Figure 14A, which indicated that the signal to noise ratio of the measured data is lower than that in data 1.

Figure 14B shows the TFA result of the received signal, indicating four normal modes with strong energy. However, the first mode in the Figure is not the true mode 1 in the received signal; that is, the first few orders of the normal mode in the signal are missing. We therefore must determine the number of normal mode to achieve the inversion of parameters. The dispersion curves in Figure 14B are the parts after band selection and serve as the input signal. The search intervals, reference values, and inversion results of the parameters to be inverted are shown in Table 3.

The comparison result of dispersion curves obtained with the different methods is shown in Figure 15, which shows that the inversed dispersion curves are in good agreement with those from the other methods. From Figure 15 and Table 3, the orders of normal modes used for inversion are 3, 4, 5, and 6, and the mode identification results are accurate. Consequently, the first two modes did not appear in the time-frequency domain of the received signal and warped signal.

The error of the distance estimation results with respect to those obtained by GPS is only 0.26%, thus further validating the accuracy of the source ranging. In addition, the estimation errors of other parameters are within 5%.