⇦Chinese Physics Letters, 2018, Vol. 35, No. 4, Article code 044301⇨ Sequential Parameter Estimation Using Modal Dispersion Curves in Shallow Water * Xue-Dong Zhang(张雪冬)1,2, Li-Xin Wu(吴立新)1,2**, Hai-Qiang Niu(牛海强)1, Ren-He Zhang(张仁和)1,2 Affiliations 1State Key Laboratory of Acoustics, Institute of Acoustics, Chinese Academy of Sciences, Beijing 100190 2University of Chinese Academy of Sciences, Beijing 100049 Received 28 November 2017, online 13 March 2018 ^*Supported by the National Natural Science Foundation of China under Grant Nos 11434012, 11774374, 11404366 and 41561144006.
^**Corresponding author. Email: wlx@mail.ioa.ac.cn Citation Text: Zhang X D, Wu L X, Niu H Q and Zhang R H 2018 Chin. Phys. Lett. 35 044301 Abstract Existing sequential parameter estimation methods use the acoustic pressure of a line array as observations. The modal dispersion curves are employed to estimate the sound speed profile (SSP) and geoacoustic parameters based on the ensemble Kalman filter. The warping transform is implemented to the signals received by a single hydrophone to obtain the dispersion curves. The experimental data are collected at a range-independent shallow water site in the South China Sea. The results indicate that the SSPs are well estimated and the geoacoustic parameters are also well determined. Comparisons of the observed and estimated modal dispersion curves show good agreement. DOI:10.1088/0256-307X/35/4/044301 PACS:43.30.Pc, 43.60.Pt, 43.60.Wy © 2018 Chinese Physics Society Article Text Underwater acoustic waveguides, which consist of sound speed profiles (SSPs), surface boundary, geoacoustic properties, topography, etc., are important environmental components affecting underwater acoustic propagation. Acoustic waves can cover a large area and can carry large amounts of oceanic information, making environmental estimation feasible. Environmental parameter estimation was initially proposed by Munk and Wunsch.[1] It related the sound ray arrival time disturbances with sound speed fluctuations, which was referred to as tomography. After that, Tolstoy pioneered the matched-field processing into ocean acoustic tomography to estimate the environmental parameters.[2] However, there is also substantial disagreement when predicted and observed measurements are compared directly at specific times and locations. Sequential estimation can provide a state-space model for predicting and correcting the state vector of a system as measurement data become available. It is demonstrated to be a powerful estimation tool, employing predictions from previous estimates and corrects stemming from sound propagation models that relate acoustic measurements to the environment parameters. The sequential estimation method maintains a small calculation and a robust behavior.[3] The application of the sequential method to parameter estimation was initially presented in 1993,[4] and it was based on the state-space representation of the normal-mode propagation model, which was expanded in a Taylor series to employ the extended Kalman filter. Henceforth, in an environmental parameter estimation problem, several works used more advanced single sequential filters,[5,6] and some hybrid methods, such as the weighted ensemble Kalman filter (EnKF)[7] and the ensemble Kalman-particle filter.[8] All of the sequential estimation works mentioned above utilized the acoustic pressure collected from a line array as the observed measurements. There is little discussion concerning other acoustical measurements. This work proposes a sequential parameter estimation method using the modal dispersion curves of the received signals from one hydrophone. The single receiver method reduces the experimental cost and complexity of the recording systems, and it does not consider the fluctuation of the receiving array in inversion problems.[9-11] The EnKF,[12] which treats the probability distribution of the states and measurements as the distribution of ensemble members and performs well in highly nonlinear systems, is employed as an example of sequential methods. The modal dispersion curves are extracted by warping transform, which is demonstrated to be a reliable tool to separate the modes.[9-11,13,14] The sequential parameter estimation is a state estimation problem of a dynamic system. It is described by the state-space model and requires two equations: the state equation and the measurement equation. The state equation models the evolution of the state vector over time. The measurement equation performs the mapping from state vectors to observation vectors. This work is concerned with estimating the evolution of environmental parameters by time-varying acoustic observations. In this case, the state-space model is represented by $$\begin{align} {\boldsymbol x}_k =& f({\boldsymbol x}_{k-1})+{\boldsymbol v}_{k-1},~~ \tag {1} \end{align} $$ $$\begin{align} {\boldsymbol y}_k =& h({\boldsymbol x}_k)+{\boldsymbol w}_k,~~ \tag {2} \end{align} $$ where Eq. (1) is the state equation and Eq. (2) is the measurement equation. The state vector ${\boldsymbol x}_k$ contains the time-dependent SSPs in the water column at discrete time $k$ and the time-independent geoacoustic parameters, and $f(\cdot)$ is the state prediction operator, which models the evolution of the state vector and represents an interdisciplinary ocean physics model. In this work, $f(\cdot)$ is taken as the unit matrix ${\boldsymbol I}$ because the update rate of the state vector is higher than the change of the study area. The state noise ${\boldsymbol v}_k$ denotes the error between the utilized prediction operator $f(\cdot)$ and the true one, and it is usually treated as the Gaussian white noise. The measurement vector ${\boldsymbol y}_k$ is the arrival times of different modes at different frequencies. The function $h(\cdot)$ is related to the state vector ${\boldsymbol x}_k$ and the measurement vector ${\boldsymbol y}_k$. It is treated as the numerical acoustic propagation model, Kraken,[15] as the study area is considered to be range-independent. Here ${\boldsymbol w}_k$ is the measurement noise, which describes the error between the observed measurement vector ${\boldsymbol y}_k$ and the predicted measurement vector $h({\boldsymbol x}_k)$. The EnKF, in principal one sequential Monte Carlo approach, relies on the ensemble of samples to describe the probability density of the state and measurement vectors. It takes the mean of the state vector ensemble as the values of the estimated parameters. The initial ensemble $\{{\boldsymbol x}_1(i),i=1, 2,\ldots,N\}$ is assumed to be taken from a normal distribution, where $i$ refers to the $i$th ensemble member, and $N$ is the number of ensemble members. When assuming that the state noise ${\boldsymbol v}_k$ and the measurement noise ${\boldsymbol w}_k$ are uncorrelated Gaussian distributions and are parameterized by the zero-mean and covariance ${\boldsymbol Q}_k$ and ${\boldsymbol R}_k$, respectively, then the state-space model in Eqs. (1) and (2) can be rewritten as $$\begin{align} {\boldsymbol x}_k^{\rm c}=\,&{\boldsymbol x}_k^{\rm p}+{\boldsymbol K}[{\boldsymbol y}_k-h({\boldsymbol x}_k^{\rm p})],~~ \tag {3} \end{align} $$ $$\begin{align} {\boldsymbol K}=\,&{\boldsymbol P}_k^{\rm p} {\boldsymbol H}^{\rm T} ({\boldsymbol HP}_k^{\rm p} {\boldsymbol H}^{\rm T}+{\boldsymbol R}_k)^{-1},~~ \tag {4} \end{align} $$ where ${\boldsymbol x}_k^{\rm p}$ and ${\boldsymbol P}^{\rm p}_k$ denote the predicted state vector generated by the state equation and its covariance matrix, ${\boldsymbol x}_k^{\rm c}$ stands for the state vector corrected by the measurement equation, the superscript T is the transpose operator, ${\boldsymbol K}$ denotes the Kalman gain, and ${\boldsymbol H}$ refers to the linearized operator of $h(\cdot)$. The EnKF utilizes the statistics of the ensemble member to represent the probability density. Thus the empirical average $\bar{{\boldsymbol x}}_k^{\rm p}$ and sample covariance $\hat{{\boldsymbol P}}_k^{\rm p}$ of the predicted state vector ensemble ${\boldsymbol x}_k^{\rm p}$ at time $k$ are respectively given by $$\begin{align} \bar{{\boldsymbol x}}_k^{\rm p}=\,&\frac{1}{N}\sum_{i=1}^N{{\boldsymbol x}_k^{{\rm p},i}},~~ \tag {5} \end{align} $$ $$\begin{align} \hat{{\boldsymbol P}}_k^{\rm p}=\,&\frac{1}{N-1}\sum_{i=1}^N({\boldsymbol x}_k^{{\rm p},i}-\bar{{\boldsymbol x}}_k^{\rm p})({\boldsymbol x}_k^{{\rm p},i}-\bar{{\boldsymbol x}}_k^{\rm p})^{\rm T}.~~ \tag {6} \end{align} $$ Obviously, when $N$ converges to infinity, $\hat{{\boldsymbol P}}_k^{\rm p}$ converges to the true error covariance ${\boldsymbol P}_k^{\rm p}$, as $\bar{{\boldsymbol x}}_k^{\rm p}$ converges to the true state vector. Since the measurement operator $h(\cdot)$ in the acoustic parameter estimation is a wave equation involving partial differential terms, it is complex to linearize this operator. However, the analytical expression of ${\boldsymbol H}$ is not required because the terms in the state-space model, which contain ${\boldsymbol H}$, can be replaced by $$\begin{alignat}{1} \!\!\!\!\!\!\!{\boldsymbol P}_k^{\rm p} {\boldsymbol H}^{\rm T}=\,&\frac{1}{N-1}\sum_{i=1}^N{({\boldsymbol x}_k^{{\rm p},i}-\bar{{\boldsymbol x}}_k^{\rm p})[h({\boldsymbol x}_k^{{\rm p},i})-h(\bar{{\boldsymbol x}}_k^{\rm p})]^{\rm T}},~~ \tag {7} \end{alignat} $$ $$\begin{align} {\boldsymbol HP}_k^{\rm p} {\boldsymbol H}^{\rm T}=\,&\frac{1}{N-1}\sum_{i=1}^N[h({\boldsymbol x}_k^{{\rm p},i})-h(\bar{{\boldsymbol x}}_k^{\rm p})][h({\boldsymbol x}_k^{{\rm p},i})\\ &-h(\bar{{\boldsymbol x}}_k^{\rm p})]^{\rm T}.~~ \tag {8} \end{align} $$ From the above discussion, every state ensemble member ${\boldsymbol x}_k^{{\rm p},i}$ generated by the state equation can be corrected as $$ {\boldsymbol x}_k^{{\rm c},i}={\boldsymbol x}_k^{{\rm p},i}+{\boldsymbol K}[{\boldsymbol y}_k^i-h({\boldsymbol x}_k^{{\rm p},i})],~i=1,2,\ldots,N .~~ \tag {9} $$ The corrected state ensemble ${\boldsymbol x}_k^{\rm c}$ will be used as the initial value to recursively generate the predicted state ensemble ${\boldsymbol x}_{k+1}^{\rm p}$ for the next time. In the state-space model, the measurement vector ${\boldsymbol y}_k$, which is the travel time for different modes at different frequencies, can be extracted via warping transform. According to the normal mode theory, the time-domain received signal $s(t)$ can be expressed as the summation of $M$ normal modes $$ s(t)=\sum_{m=1}^M{A_m(t)e^{j{\it \Phi}_m(t)}},~ m=1, 2,\ldots,M,~~ \tag {10} $$ where ${\it \Phi}_m(t)$ and $A_m(t)$ denote the phase and instantaneous amplitude of the $m$th mode at time $t$, respectively, and $M$ is the number of modes propagating in the waveguide. The ideal waveguide model consists of an isovelocity water column between a pressure-release upper boundary and a rigid bottom. For such a model with depth $D$, range $r$ and velocity $c$, the phase ${\it \Phi}_m(t)$ can be expressed as $$ {\it \Phi}_m(t)=2\pi f_{\rm cm}\sqrt{t^2-(r/c)^2},~~ \tag {11} $$ where the cutoff frequency for mode $m$ is determined by $f_{\rm cm}=(2m-1)c/4D$. Thus the received signal in Eq. (10) can be expressed as $$ s(t)=\sum_{m=1}^M{A_m(t)e^{i2\pi f_{\rm cm}\sqrt{t^2-(r/c)^2}}}.~~ \tag {12} $$ In this form, the phase is nonlinear with time, but is linear with the time related operator $g(t)=\sqrt{t^2+(r/c)^2}$. The warping transform utilizes this operator to map the time domain and the warping domain, and transforms the received signal to the warped signal $$ W_{\rm g}s(t)=\sum_{m=1}^M{\sqrt{g'(t)}A_m[g(t)]e^{j2\pi f_{\rm cm}t}},~~ \tag {13} $$ where $\sqrt{g'(t)}$ ensures the energy balance after the warping transform. Thus the warped signal is represented by the summation of the $M$ single-frequency signals with frequency $f_{\rm cm}$ and amplitude $\sqrt{g'(t)}A_m[g(t)]$. Next, filter the warped signal $W_{\rm g}s(t)$ to obtain the single-frequency signal in the warping domain. Then, employ the inverse warping operator $g^{-1}=\sqrt{t^2-(r/c)^2}$ to obtain the single mode of the received signal in the time domain. Thus the modal dispersion curves are able to be extracted and adopted as the measurement vector ${\boldsymbol y}_k$ in the EnKF method. Although the introduced warping transform is modeled by an ideal waveguide, the warping transform generally features robust behavior to environment mismatch and performs well in most of the shallow water environment.

Fig. 1. (a) Time domain pulse-compression signal and (b) its spectrogram.

Fig. 2. The spectrograms of (a) the warped signal and (b) the extracted modal dispersion curves, where the modes are labeled with white numbers.

To validate the sequential parameter estimation method, an experiment was conducted in the South China Sea in June 2010. The waveguide was chosen to be range-independent, with the source/receiver range $r$ of 39.3 km and the water depth $D$ of 98.2 m. The stationary source, which was 7.4 m above the seafloor, emitted 160–250 Hz linear frequency modulation signals with a 1 s time length and 6 times a minute. The acoustic signal was collected by a horizontal line array laid on the bottom; the 61 signals received by one hydrophone from 8 June, 13:00 to 14:00 are used in this work. Pulse compression is implemented to enhance the signal-to-noise ratio and an example of the pulse-compression signal and its wavelet spectrogram are shown in Fig. 1. One can note that modes 1, 2 and 3 overlap and are indistinguishable. Figure 2(a) illustrates the short-time Fourier spectrogram of the warped signal. It is obvious that the warped modes are not perfectly horizontal because of the slight mismatch between the ideal waveguide model and the experimental environment, but modes 1–5 are clearly separated. Since the energy of mode 3 is not strong enough to be extracted, the warping transform is performed on modes 1, 2, 4 and 5. The modal dispersion curves are extracted as illustrated in Fig. 2(b), and the modes are labeled to ease the reading.

Fig. 3. The range-independent environmental model used in the estimations. The source is 7.4 m above the seafloor and the receiver is on the seabed.

Fig. 4. (a) The first three EOFs computed from the SSP database, weighted by their respective eigenvalue (black solid curve: EOF$_1$, blue dashed curve: EOF$_2$, and red dash-dot curve: EOF$_3$) and (b) the percentage of the first 15 EOF energy to the total terms.

The environmental model given in Fig. 3 is used for each sequential estimation time step. As mentioned above, since the waveguide is range-independent, the normal-mode numerical model, Kraken, is considered as the measurement operator to calculate the predicted arrival times with the parameters of the state vector. The state vector consists of 7 parameters including the first three empirical orthogonal function (EOF) coefficients[16] of the SSPs, sediment thickness $h_{\rm sed}$, sediment sound speed $c_{\rm sed}$, bottom sound speed $c_{\rm bot}$, and water depth $D$. The water depth is also taken into account because of the possible existence of tides in the study area. In this case, the sediment density $\rho_{\rm sed}$ and the bottom density $\rho_{\rm bot}$ vary with their sound speed according to the Hamilton empirical function $$ c=2330.4-1257.0\rho+387.7\rho^2,~~ \tag {14} $$ and they are not included in the state vector (this hybrid estimation method was presented in Ref. [17]). The attenuation coefficients of the sediment and bottom layer are not considered due to the fact that they are not sensitive to the modal dispersion curves. The data from this waveguide were also estimated in Ref. [10]. An EOF analysis of the SSPs is carried out using 2 moored T-string systems deployed at the source and receiver, respectively, from 7 June, 14:30 to 8 June, 20:30, giving a mean profile together with the first three EOFs (accounting for 92% of the variance in SSPs), in Fig. 4. The SSPs are estimated by tracking the EOF coefficients. The average values of the two T-string data at the time when the source emits signals are taken as the true SSPs.

**Table 1.** The initial settings and results of the EnKF parameter estimation.
Parameter	Representation	Unit	Initial value	${\boldsymbol P}^{1/2}_0$	${\boldsymbol Q}_k^{1/2}$	Search bounds	Estimated results
EOF coefficients	$\alpha_{1:3}$		[1,1,1]	1	1	[$-$20 20]
Sediment thickness	$h_{\rm sed}$	m	10	1	1	[0 30]	8.2
Sediment sound speed	$c_{\rm sed}$	m/s	1530	0.01	20	[1530 1650]	1574
Sediment density	$\rho_{\rm sed}$	g/cm$^3$	1.43				1.62
Bottom sound speed	$c_{\rm bot}$	m/s	1600	0.01	20	[1600 1700]	1654
Bottom density	$\rho_{\rm bot}$	g/cm$^3$	1.69				1.81
Water depth	$D$	m	90	0.01	1	[90 110]	97.0

The SSP and geoacoustic parameters are estimated using a 200-point EnKF and their initial values, initial state vector covariances, state equation noises, search bounds and results are summarized in Table 1. Generally, the sequential method does not need to set the search bounds of the state vector. However, because of the large values of the state equation noises ${\boldsymbol Q}^{1/2}_k$ for sediment and bottom sound speed, the search bounds are employed to ensure a reasonable environment model and they have little effect on the estimation results. Without loss of generality, the initial values of the first three EOF coefficients $\alpha_{1:3}$ are initialized to 1. The first three panels of Fig. 5 compare the estimated SSPs with the true SSPs, in 2-time-step increments, in which every profile is shifted to the right by 10 m/s from the previous one. By comparison, the estimated SSPs agree well with the measured SSPs. The performance of the EnKF in tracking the SSPs is tested by root mean square error (RMSE) between the estimated and true SSPs, which is defined as $$ {\rm RMSE}=\sqrt{\frac{1}{N_Z}\sum\nolimits_{z=1}^{N_Z}{|{\boldsymbol c}_{\rm true}(z)-{\boldsymbol c}(z)|^2}},~~ \tag {15} $$ where ${\boldsymbol c}_{{\rm true}}$ are the true SSPs, ${\boldsymbol c}$ are the estimated SSPs, and $N_Z$ is the length of one profile. A small RMSE indicates an excellent performance of the estimation method. The evolution of the RMSE is shown in the last panel of Fig. 5 with the time-averaged RMSE of 0.73 m/s. It is clear that the SSPs are well estimated. The evolving of the marginal posterior probability density functions (PDF) at every time step for geoacoustic parameters are given in Fig. 6, which shows clearly that the parameter expectations converge gradually, from the pre-convergent (1–40) to the convergent (41–61) stage. The value of $h_{\rm sed}$ decreases from 10 to 7.8 m. The value of $c_{\rm sed}$ increases quickly from 1530 to 1582 m/s, $b_{\rm bot}$ from 1600 to 1654 m/s, and $D$ from 90 to 97.5 m, respectively, along the track. In addition, it is obvious that $h_{\rm sed}$ and $D$ do not touch the search boundary in the convergent stage, while $c_{\rm sed}$ and $c_{\rm bot}$ are expanded within the search range. The means of the track results in the convergent stage are taken as the geoacoustic parameters estimation results (in Table 1). Even though detailed range-dependent ground true measurements do not exist, the results compare favorably with the previous study which estimated data from the same area with $h_{\rm sed}$ of 6.3 m, $c_{\rm sed}$ of 1579 m/s, $c_{\rm bot}$ of 1650 m/s, and the measured water depth of 98.2 m.[10] In that work, the parameters are estimated by a genetic algorithm with a population of 140 and generations of 500, and repeated computations on one signal are performed to ensure the cost function value converges to the minimum. Thus compared with the genetic algorithm, the sequential estimation method has less calculation and more robust behavior.

Fig. 5. The first three panels show the estimated (black solid curves) and true (red dashed curves) SSPs every 2 time steps. The last panel shows the evolution of the RMSE.

To validate this parameter estimation method, the estimated SSPs and geoacoustic parameters at the pre-convergent and convergent stages are used to calculate the travel times, and the estimated dispersion curves for modes 1, 2, 4 and 5 are compared with the experimental ones in Fig. 7. Figure 7(a) shows that the estimated curves for the signal at the pre-convergent stage do not agree well with those of the experimental signal, especially for modes 2 and 5. However, as shown in Fig. 7(b), the modal dispersion curves of the signal at the convergent stage coincide well with the experimental ones. This indicates that the presented estimation method effectively decreases the discrepancies between the estimated and experimental travel times.

Fig. 6. Estimated results as evolving marginal posterior PDFs for (a) sediment thickness, (b) sediment sound speed, (c) bottom sound speed and (d) water depth.

Fig. 7. The extracted (black solid curves) and calculated (red open circles) modal dispersion curves for the signals at (a) the pre-convergent and (b) convergent stage, and the modes are labeled with numbers.

In conclusion, the capability of the sequential method to estimate SSP and geoacoustic parameters in a spatially constant environment using modal dispersion curves has been presented. The evolutions of the estimated parameters are expressed by the expectations of marginal posterior PDFs. The SSPs are estimated at every time step, while the geoacoustic parameters converge to the true values gradually, and the final results of the geoacoustic parameters are represented by means of the state vectors at the convergent stage. The presented estimation method is employed to process the South China Sea experimental data. The results show that the estimated SSPs agree well with the true SSPs and the geoacoustic parameters converge gradually to the true values along the track. The inversion results are validated by comparing the estimated and observed travel times. References Ocean acoustic tomography: a scheme for large scale monitoringAcoustic tomography via matched field processingSequential geoacoustic inversion at the continental shelfbreakSound velocity profile estimation: a system theoretic approachInversion for Time-Evolving Sound-Speed Field in a Shallow Ocean by Ensemble Kalman FilteringGeoacoustic and source tracking using particle filtering: Experimental resultsA method for tracking time-evolving sound speed profiles using Kalman filtersTracking of time-evolving sound speed profiles in shallow water using an ensemble Kalman-particle filterSequential inversion of modal data for sound attenuation in sediment at the New Jersey ShelfBayesian geoacoustic inversion of single hydrophone light bulb data using warping dispersion analysisThe Ensemble Kalman Filter: theoretical formulation and practical implementationTheoretical analysis of warping operators for non-ideal shallow water waveguidesGeoacoustic Inversion Based on Modal Dispersion Curve for Range-Dependent EnvironmentAn underwater acoustic sound velocity data modelGeoacoustic inversion for sediments in the South China Sea based on a hybrid inversion scheme

[1]	Munk W and Wunsch C 1979 Deep. Sea. Res. 26 123
[2]	Tolstoy A et al 1991 J. Acoust. Soc. Am. 89 1119
[3]	Yardim C et al 2012 J. Acoust. Soc. Am. 131 1722
[4]	Candy J and Sullivan E 1993 IEEE J. Oceanic Eng. 18 240
[5]	Carrière O et al 2009 IEEE J. Oceanic Eng. 34 586
[6]	Yardim C et al 2010 J. Acoust. Soc. Am. 128 75
[7]	Huang J M et al 2014 J. Acoust. Soc. Am. 136 EL129
[8]	Li J L and Zhou H 2013 J. Acoust. Soc. Am. 133 1377
[9]	Duan R et al 2016 J. Acoust. Soc. Am. 139 70
[10]	Niu H Q et al 2014 Acta Acust. 39 1
[11]	Bonnel J et al 2013 J. Acoust. Soc. Am. 134 120
[12]	Evensen G 2003 Ocean Dyn. 53 343
[13]	Niu H Q et al 2014 J. Acoust. Soc. Am. 136 53
[14]	Guo X L et al 2015 Chin. Phys. Lett. 32 124302
[15]	Porter M B 1992 The Kraken Normal Mode Program Tech. Rep. DTIC Document
[16]	LeBlanc L R and Middleton F H 1980 J. Acoust. Soc. Am. 67 2055
[17]	Li Z L and Li F H 2010 Chin. J. Oceanol. Limnol. 28 990