Matched-field geoacoustic inversion using propagation invariant in a range-dependent waveguide

Good performance has been achieved with matched-field geoacoustic inversion (MFI) using controlled sound sources and a vertical line array (VLA) of receivers spanning a large portion of the water column.^1,2 Building on these efforts, other configurations with towed^3–5 and bottom-moored^6–8 horizontal line arrays (HLAs) have also been widely studied. Because MFI is a model-based inversion method, the key to successfully estimating the geoacoustic parameters of marine sediments is to develop accurate forward propagation models that can reasonably reproduce the measured acoustic fields.

In a real ocean environment, the bathymetry and sediment parameters of the seabed can be quite range dependent, and the sound speed profiles (SSPs) of the water column can vary in both space and time. These variabilities may strongly degrade the performance of MFI. First, it has been shown that in a range-dependent waveguide, the uncertainty of the inversion results increases with growing separation between the source and the VLA,^9,10 and it is recommended that this separation should be less than 2 km.¹⁰ Second, conventional range-independent MFI can yield only range-averaged properties of the sea bottom, which are not the results that are needed, especially when distinct bottom types are present along the sound path; instead, the exact “local” information of the sea bottom is needed for typical applications of matched-field processing (MFP). In this paper, the term “local” refers to results relating to a subarea, in contrast to results obtained by averaging over the whole propagation path. Third, unknown changes in the water column sound speed that may be caused by features such as internal waves give rise to uncertainty in the inverted results for the bottom sediments due to the coupling between the parameters of the water column and the sea bottom.¹¹ 

To address these problems, the forward model must be able to account for range-related variations, and this involves using methods such as adiabatic approximation,¹² mode coupling theory, or even more complicated theories, e.g., transport equations. These range dependent models are computationally time consuming, and their increased complexity poses a risk of unstable inversion results.¹³ By contrast, in this paper, instead of trying to describe this propagation process with complicated propagation models, we seek a way to eliminate the influence of the inhomogeneities in the intermediate zone, a concept that we refer to as propagation invariant (PI).¹⁴ 

In this paper, PI is used for geoacoustic inversion based on an experimental configuration consisting of a moving source and a fixed VLA that spans the whole water column. In Sec. 2, the basic theory of PI and the relevant MFI process are introduced. In Sec. 3, simulations are presented to illustrate the feasibility of this method. Summary and conclusions are given in Sec. 4.

In this subsection, a concrete realization of the concept of PI is presented in a two-dimensional (axially symmetric) acoustic waveguide. A typical experimental setup consisting of a fixed VLA of $NZ$ hydrophones and a moving source continually transmitting signals when moving away from the VLA, as shown in Fig. 1, is considered. The hydrophones of the VLA are located at depths of $z1,z2,…,zNz$⁠, respectively, while the source is assumed to remain at a fixed depth of $zs$⁠. We prefer a segment by segment inversion scheme, i.e., dividing the source path into many segments and performing the inversions segment by segment. For each segment, there are $Ns$ groups of signals to be processed, and the corresponding source-receiver distances are denoted by $r1,r2,r3,…,rNs$ as shown in Fig. 1. Our purpose is to invert the seabed acoustic parameters of each segment, without referring to the environments of previous segments.

For simplicity, in this section, it is assumed that the bathymetry, seabed type, and SSPs vary slowly, such that the adiabatic approximation can be applied in the waveguide.

For each segment, we can obtain an $Nz×Ns$ complex matrix of data at the frequency $f0$⁠,

where the frequency dependence is suppressed for notational simplicity. Based on the adiabatic approximation, the elements of the matrix P may be written as follows:

where $q=1,2,…,Nz$⁠, $l=1,2,…,Ns$⁠, $M$ denotes the effective number of trapped modes, $km(r)$ denotes the local complex horizontal wavenumber of the $mth$ mode at range r, and $ϕm$ denotes the corresponding orthogonal modal depth function. For the trapped modes, the following orthogonal relation is approximately satisfied:

Now, we rewrite Eq. (1) in the matrix product form

where

and $Λ=diag{e−i∫0r1km(x)dx/∫0r1km(x)dx}$ is the acoustic propagation matrix from r₁ to the VLA. The matrices in Eq. (4) have the following physical interpretations: the matrix $Φ$ is the receiver weight matrix, $Λ$ is the propagation matrix from the source position corresponding to the first signal to the VLA receivers, and V is the source weight matrix.

To define PI, we choose the first $(Ns−1)$ columns of P to form a new matrix $Pr$ and the last column to form $P1$⁠. It can be inferred that the difference between $Pr$ and $P1$ arises only from the source weight matrix or vector, denoted as $Vr$ and $V1$⁠, respectively, while $Φ$ and $Λ$ are common. PI can be achieved by means of the transformation $Pr+P1$⁠, and from the above equations, one can obtain

where the superscript “+” denotes the generalized inverse matrix,

the superscript “H” represents the complex conjugate transpose operation. The optimal regularization parameter $μ$ can be selected by means of the ordinary cross-validation (OCV) method or the L-curve criterion.¹⁵ 

In this paper, we refer to $Pr+P1$ as the PI array vector and $Vr+V1$ as the PI source vector; the PI array vector is calculated from the array data while the PI source vector is related to the relative source locations in the selected segment and the local acoustic mode wavenumbers at these source locations. From Eq. (6), it can be seen that the PI array vector is equal to the PI source vector, and the propagation matrix plays no role in this relationship. In this sense, the PI array vector is independent of the intermediate propagation process, thus giving rise to the term “propagation invariant.”¹⁴ PI is a kind of data transformation technique that can get rid of intermediate propagation and produce a measurement of the local “sound pressure field.”

It should be remarked that the adiabatic approximation is only adopted for simplicity of the derivation of PI; as a matter of fact, the concept of PI holds for the case when acoustic modal coupling occurs. In such cases, the pressure field can still be expressed in the matrix product form in Eq. (4). The difference lies in the matrix $Λ$⁠, whose specific expression is related to the environment from r₁ to the VLA. However, this does not destroy the concept of PI because $Λ$ is eliminated in the PI array vector.

In this subsection, we will show how PI can benefit MFP-based geoacoustic inversion in a range dependent waveguide. We will use the PI array vector for matched field inversion, and this can offer two advantages. First, the PI array vector has nothing to do with the intermediate propagation process in the waveguide, so that the forward model is only required to describe the propagation from the first source range r₁ to the last r_Ns in the selected segment, as in Eq. (6). This can highly reduce the complexity of the forward model. Second, the inversion result based on the PI array vector is the local one, that is, only related to the geoacoustic parameters of the selected segment, and the range variance of sediment property can be acquired segment by segment independently. This is why we need PI for MFP-based geoacoustic inversion.

The process of PI-based MFI can be summarized as follows:

1. Divide the path into different segments.

2. Calculate the PI array vector $X=Pr+P1$ at selected segment.

3. Use the forward model to calculate the replica, i.e., PI source vector: Y = V_r⁺V₁. The forward model is established using a range independent model. The input parameters of the forward model include: SSP of the water column, water depth, sediment properties of the local area, and the source ranges of the selected segment. The forward model first uses the local environment parameters to calculate the modal wavenumbers, then calculates the PI source vector using a range independent version of Eq. (6).

4. The cost function $Q=YHX/XHXYHY$ is used to generate the optimized inversion result.

Among the input parameters of the forward model, the source ranges must be known and can be provided by a Global Positioning System (GPS) instrument, while the water depth and SSP of the water column can be either provided by echo sounder and CTD instruments, or inverted together with sediment properties. For the multi-parameter inversion problems, it is expected that the global optimization methods, such as simulated annealing and genetic algorithms, can also benefit from the concept of PI.

In the simulations, as a proof of the concept, we intend to invert the range variance of the sediment property. We choose to invert only the sediment sound speed $cb$ with all other sediment parameters held fixed at their known values. The reasoning is that compared with sediment density and attenuation, the sediment sound speed is the most sensitive parameter.³ In the forward model, the SSP in the water column and water depth are also assumed to be known, taking the value at the mid-point of the segment.

The geometry used in the simulations is shown in Fig. 1. The bathymetry and seabed sediment properties both vary with range. There is a gentle slope in the waveguide, which starts at a range of 4 km, with a sea depth of 140 m, and ends at a range of 12.5 km, with a sea depth of 105 m; the other areas of the waveguide have a flat bottom. The bottom sediment is set to be a half-space liquid medium. For the first flat area (from 0 to 4 km), the properties of the sediment are constant, with values of $cb=1650$ m/s, $ρb=1.4$ g/cm³, and $αb=0.25$ dB/λ, and for the second flat area (from 12.5 to 15 km), the properties of the sediment are also constant, with values of $cb=1561$ m/s, $ρb=1.25$ g/cm³, and $αb=0.1$ dB/λ. In the slope area, the properties of the sediment vary linearly from those of the first flat area to those of the second flat area. For simplicity, it is assumed that the sound speed in the water column is 1500 m/s.

In the simulation, the 70 hydrophones of the VLA are equally spaced to cover the whole water column, as shown in Fig. 1. As the source moves away from the VLA at a speed of 10 kn and a fixed depth of 5 m, it emits a linear frequency modulated (LFM) signal (from 150 to 250 Hz) every 10 s during its movement from 1 to 15 km.

The sound pressure fields are calculated based on the kraken program model using the adiabatic approximation, and the signal waveforms are calculated via inverse fast Fourier transform (IFFT) transformations. The spectrum levels of the source (denoted by SL) are 195, 185, 175, and 165 dB, respectively, while the spectrum level of the ambient ocean noise (assumed to be Gaussian white noise) is 95 dB.

Figure 2 shows the sound transmission loss (in dB) from the acoustic source to the hydrophones in the VLA for the setup depicted in Fig. 1. Two examples of the waveform and time-frequency spectrogram of the received signal at a depth of 100 m when the source-receiver distance is 14.5 km are also given in Fig. 2. The transmission loss is approximately 70 dB when the source is at 14.5 km; thus, a high signal-to-noise ratio (SNR) is observed in the received signal when SL ≥ 185 dB, while the SNR is low when SL = 165 dB.

In the following, we show the results of MFI based on the concept of PI under different SNR conditions.

We use the segment by segment inversion scheme to examine the evolution of the sediment sound speed with range. The effect of segment aperture length will be discussed below. In this simulation the aperture length is set as 1000 m with an element space of $Δr=50 m$⁠. There is a total of 14 segments along the whole path with no override, labeled with Arabic numerals, see Table 1. For each data segment, there are 20 groups of signals; the first 19 VLA signals are used to form the matrix $Pr$⁠, while the last VLA signal is treated as $P1$⁠. The PI array vector is calculated by $X=Pr+P1$⁠, while the replica PI vector $Y=Vr+V1$ is calculated using the forward model, as described in Sec. 2.2.

Before performing MFI, we first examine the validity of PI for various SNR conditions. For the purpose of validity, the PI source vector Y is calculated by the forward model with the sediment sound speed known and the validity of PI can be measured by the cost function Q. Figure 3(a) shows the cost function values of Q for all 14 segments with different SL values. High values of Q are observed for high SL values, while when SL = 165 dB Q is small for segments of large distance. This implies the validity of the PI on the condition of high SNR. It should be added that when calculating the generalized inverse matrix using Eq. (7), the value of μ is highly related to SNR conditions. When SNR is high, μ takes the value of 0.0033λ₀, where λ₀ is the maximum of the singular values of $Pr+P1$⁠; while when SNR is very low, μ takes the value of 0.1622λ₀.

We also discussed the effect of segment aperture length, shown in Fig. 3(b). One can observe that Q > 0.98 for all aperture lengths when SNR is high. Larger aperture length gives smaller Q because we use a range independent version of Eq. (6). From Fig. 3(b) one can conclude that the choice of aperture length is very flexible as long as SNR is high. However, when the aperture length is small the sidelobes are high (not shown here) and the relatively large aperture length can give more robust results. Thus, the aperture length we use in the simulation is 1000 m.

Now, PI-based MFI can be performed. In the inversion process, the bottom sound speed is inverted through an exhaustive search in the interval from 1510 to 1800 m/s. We use signals at 11 distinct frequencies (in 10 Hz intervals) in the band from 150 to 250 Hz to separately invert the sediment sound speed, and the average and uncertainties of the inversion results are then obtained. The inversion results for the bottom sound speed $cb$ obtained by following the steps in Sec. 2.2 when SL = 165, 175, 185, and 195 dB are shown in Fig. 4(a).

From Fig. 4(a), one can see that the inversion results are highly accurate when SL = 185 and 195 dB, while the inversion results show considerable deviations when SL = 165 dB. Additionally, the uncertainty of the inversion results increases as the SNR decreases. This is because the generalized inverse operation used in the PI method depends on the SNR, and a high SNR is usually required. Besides, it is seen that the inverted values are systematically a little larger than the true values. This is logical because we use a range independent model to describe the sound propagation through the slope area.

For comparison with PI-based MFI, two other inversion schemes are also considered. The first is conventional MFI. In conventional MFI, the simulated sound pressure fields are directly matched with the replica fields. In the conventional MFI scheme, the sediment density, attenuation and bathymetry, which are taken to be known, are range dependent, while the sediment sound speed is range independent. To calculate the replica pressure field, we apply the kraken model using the adiabatic approximation. Note that the inverted sediment sound speed is the range-averaged along the sound path in the conventional MFI scheme.

The second scheme is MFI based on the concept of phase conjugation (denoted as PC-MFI). Phase conjugation or backpropagation has been applied in source location, geoacoustic inversion, and underwater communication.^16–18 In the PC-MFI scheme, we use $PrHP1$⁠, similar to the PI array vector, to provide the local measurement which is then used for local sediment property inversion; the replica is calculated by $VrHV1$⁠, similar to Eq. (6). Combined with the segment by segment approach, PC-MFI can also predict the range variance of sediment property. What we do in the PC-MFI is technically different from that in Ref. 15; we use phase conjugation to provide the local measurement for inversion, while Thode used this to produce a virtual receiver array with the help of a guide source; however, the idea of using phase conjugation to compensate for the intermediate propagation process is common.

The PI-based, conventional and PC-based MFI results are plotted in Fig. 4(b) and numerically listed in Table 1. It is seen that PI-based MFI enables good estimation of the local sediment sound speed along the whole path, while PC-based MFI achieves similar performance but with generally larger deviations. By contrast, the “conventional” MFI method can produce good estimates at close distances but shows an increasing deviation at larger distances. This is because conventional MFI can estimate only the average sediment properties along the sound path, while the PI-MFI and PC-MFI method can invert the local properties.

The reason why PI-based MFI outperforms PC-based MFI lies in the generalized inverse operation. The generalized inverse and complex conjugate transpose operations can both provide the local “sound pressure field,” but the complex conjugate transpose operation applied in PC-based MFI can only compensate for the phase terms of the acoustic modes, while the generalized inverse operation applied in PI-based MFI can also compensate for geometric spreading and bottom attenuation in addition. Consequently, in the PI-based MFI scheme, there is more energy in higher acoustic modes, which usually contain more information about the seabed sediment, and consequently, PI-based MFI offers higher accuracy than PC-based MFI does.

In this paper, a scheme for matched-field geoacoustic inversion (MFI) based on the concept of propagation invariant (PI) in a range-dependent waveguide is presented. An experimental geometry consisting of a fixed vertical line array (VLA) of hydrophones and a moving source ship is considered. By the segment by segment inversion scheme, PI-MFI can offer the range variance of bottom sediment property. The proposed inversion scheme is verified by simulations of a case in which there is a gentle slope in the waveguide.

The PI approach can be used to eliminate the effects of intermediate propagation by taking the (generalized) inverse of the VLA data segment. This can be done because the information of the propagation matrix is contained in the received data matrix and can be canceled out by multiplying the inverse of one data matrix by another. The advantage of using the PI concept in MFI is that it can reduce the complexity of the forward propagation model and provide the local sediment property. The method of phase conjugation can also achieve similar results, although it shows inferior performance.

The proposed PI concept serves as the basis for an efficient method of geoacoustic inversion, and we believe that it can also be beneficial for source location and underwater communication in complex environments.

Thanks are due to Professor Ning Wang in the Department of Marine Technology, Ocean University of China, for valuable discussions and suggestions. This research was supported by the National Natural Science Foundation of China (Grant No. 11674294) and by the Fundamental Research Funds for the Central Universities (Grant No. 201713056).