Applicability of the single-scattering approximation for the ocean acoustic sounding

A mathematical model in the work that describes the process of pulsed high-frequency acoustic sounding in a fluctuating ocean is considered. The inverse problem for the non-stationary equation of acoustic radiative transfer is investigated. It consists in finding the coefficient of volume scattering from the angular distribution of the radiation flux density at a given point in space. In the single scattering approximation, a formula is obtained to find the scattering coefficient. A series of computational experiments was carried out to substantiate its use in acoustic sensing of a fluctuating water medium. A numerical analysis has shown that the application of the single scattering approximation is justified at least at a sensing range of the order of 100 meters. Moreover, double and triple scattered fields have the main effect on the error.


Introduction
In this study we consider the non-stationary radiative transfer equation as applied to acoustic sensing of fluctuating inhomogeneous medium [1][2][3][4][5][6][7]. The inverse problem for the radiative transfer equation is studied, which consists in finding the volume scattering coefficient of sound if the angular distribution of the radiation intensity at a certain point in unbounded space is given. This problem can be solved using the single-scattering approximation [7].
The aim of this work is the numerical verification of the allowability of the single-scattering approximation to solve the problem of determining the volume scattering coefficient in the marine medium. Similar questions were considered by many experts. For example, in [2], a simplified analysis of the applicability of the single-scattering acoustic wave approximation in problems of hydrobionts accumulations in the stationary case was carried out. Using sufficiently strict limitations, it was shown that at a sufficiently high density of hydrobionts in water and a high frequency of probe radiation (of the order of 10 − 100 kHz), the single-scattering approximation is insufficient to solve direct and inverse problems.
In this paper, we consider the admissibility of the single-approximation in the non-stationary case. First of all, we will be interested not in the proximity of the direct problem solution and its approximation, but in the proximity of the solutions to the corresponding inverse problems.
To this end, Monte Carlo method is constructed for numerically finding a solution to a direct problem. Its distinguishing feature is the fact that the radiation flux density corresponding to a scattering multiplicity of at least two is calculated by a statistical method and a single-scattered radiation is analytically found. Computational experiments were performed demonstrating the parasitic effect of multiple scattering during the reconstruction of an unknown medium.

Statement of the inverse problem
The non-stationary process of high-frequency wave fields propagation in a scattering medium is described by the integro-differential equation [6]: where r ∈ R 2 , t ∈ [0, T ] and the wave vector k belongs to the unit circle Ω = {k ∈ R 2 : |k| = 1}. The function I(r, k, t) is interpreted as the energy flux density of a wave at the time instant t at the point r, which propagates in the direction k with the velocity c. The quantities µ and σ are called the attenuation and scattering coefficients, and the function J describes the sources of the sound field. Let the initial condition for equation (1) have the form and we will assume that the function J describes a point radiation source concentrated at the origin, emitting a unit power pulse at the time t = 0: where δ is the Dirac delta function and I ± (r, k, t) = lim Assumption (3) about the structure of the radiation source is not only justified from a physical point of view, but also allows us to significantly simplify the study of a wide range of problems in mathematical physics [7][8][9][10].
The direct problem for the transport equation (1) is the problem of determining the function I from equation (1) and the initial condition (2) for all given coefficients (c, µ, σ, J). It is assumed that c and µ are scalar quantities and σ(r) is a piecewise continuous function in R 2 and σ(r) ≤ µ.
We will consider the inverse problem, which consists in determining the function σ from ratio (1), (2), (3) and the additional condition in which the quantities c, µ and the function P are known.

Representation of the direct problem solution in the form of a Neumann series and its summation using Monte Carlo method
The solution of the Cauchy problem (1), (2) is equivalent to the equation solution of the integral type, which for an unbounded domain has the following form [7]: The equation solution (5) can be found in the form of the Neumann series [11] I(r, Because of the radiation source J has the form of (3), the singular component I 0 of the Neumann series (6) contains the multiplying delta functions, and the sum I 1 +I 2 +... represents the regular part.
Since I + 0 (0, k, t) = 0, the solution of the Cauchy problem at the point r = 0 can be found in the form of the Neumann series In (12) for n = 2, the integration is carried out only with respect to the variables k 1 , τ 1 in exactly the same way as is done in relation (10). We use Monte Carlo method [11][12][13][14] to summarize the Neumann series (12). Denote the mathematical expectation of a random variable Θ by E[Θ], then for the functions I n at the point (r 0 , k 0 , t 0 ), the following expressions can be written down as In expressions (13), (14), the trajectory points (r 0 , k 0 , t 0 ), (r 1 , k 1 , t 1 ), , ..., (r n−1 , k n−1 , t n−1 ) are defined as follows: where k i is a random vector uniformly distributed on the unit circle Ω, τ i is a random variable uniformly distributed on a segment [0, τ (r i−1 , k i−1 , t i−1 )]. Taking into account relations (13), (14) the expression for the function I + (0, k, t) takes the form where the random variables θ n are recurrently determined θ n = θ n−1 σ(r n )τ (r n−1 , k n−1 , t n−1 ), θ 0 = 1, n = 1, 2, ....
When calculating the sum in (17), we restrict ourselves to N terms, which corresponds to accounting all scattering events in a medium of order not higher than N . Then the scheme for calculating the flux density I, which corresponds to the shorten sum of the Neumann series (17) takes the following form. Trajectory (15) is constructed, and the random variable is calculated as: Repeating the described procedure M times, we obtain a sample of the size M for a random variable Θ N . The average value of the sample gives us an estimate for the mathematical expectation of the random variable Θ N and, therefore, the approximate value of I at (r 0 , k 0 , t 0 ).

Numerical analysis of the multiple scattering effect on the accuracy of solving the inverse problem in the single-scattering approximation
The function τ (r, k, t) = ct/2 at the point r = 0, therefore, formula (13) for the density of a single-scattered radiation flux takes the form Given I + 1 (0, k, t) from (24), the function σ can be found. Indeed, for k = − r |r| and t = 2|r| c from ratio (19) we obtain  (20) gives an explicit solution to the inverse problem in the single-scattering approximation [7].
According to condition (4) in the formulation of the inverse problem, the total radiation flux P (k, t) = I + (0, k, t) at the point r = 0 is considered to be known, therefore, the calculation of the function σ(r) by the formula will lead to the distortion of tomographic images of the scattering coefficient. Below we will carry out a quantitative and qualitative analysis of distortions in the calculation of the function σ of the initial data by formula (21), which corresponds to acoustic sensing of the ocean by a pulsed sound source at frequencies of the order of 100kHz [6]. For convenience of describing the numerical results, the volume scattering coefficient is denoted as σ n which is calculated by the formula where P n is the solution of the radiative transfer equation in the n multiple scattering approximation calculated by formulas (13), (14) at the point (0, 0). Thus, according to (20) σ n = σ only for n = 1 and the remaining functions σ n , n ≥ 2 are some approximations of the coefficient σ, which are distorted by the influence of multiple scattering in the medium. The sensing region with the sound source concentrated at the origin was 160m×80m (r 1 ∈ [−80, 80], r 2 ∈ [0, −80]). The speed of sound and the attenuation coefficient had the following values: c = 1500m/s, µ = 0.018m −1 . The ratio σ/µ in the main water medium was 0.1 [3] and in inclusions it changed in the range from 0 to 1. A map of the scattering coefficient values in the sensing region is presented in Figure 1   The results of restoring the functions σ n , n = 2, 3, 4, 5 by formula (22) are presented in Figures 2,3,4,5. The reconstructed structure of the medium can be quite clearly distinguished in all tomographic images despite of the appearance of "exposure" from two large strongly scattering inclusions (σ/µ = 1).    The normalized mean square error was chosen to quantify the quality of the images obtained by modeling the radiation process taking into account scattering of the n-th multiplicity. ε n = i j (σ n (r 1,i , r 2,j ) − σ(r 1,i , r 2,j )) 2 i j σ 2 (r 1,i , r 2,j ) , adequately reflecting the quality of the images as a whole. Table 1 contains the calculated values of the normalized mean square errors ε n depending on n. The table shows that the error grows with increasing n, reaching 0.634, and for n > 3 the growth of errors ε n stabilizes. This indicates that the main contribution to the error in reconstructing the function σ is exerted by the destructive effect of the double and triple scattered fields. Thus, the single-scattering approximation allows one to recreate the qualitative structure of the medium during acoustic sensing in the ocean at a range of about 100 meters, however, obtaining the ensured quantitative estimates of the volume scattering coefficient still remains in doubt.