Multifrequency near-field acoustic tomography and holography of 3D subbottom inhomogeneities

Abstract

A method of coherent multifrequency acoustic tomography and holography of spatially localized subbottom inhomogeneities in shallow seas is proposed. This method is based on solving of the near-field inverse scattering problem that makes it possible to realize a subwavelength resolution. It involves the analysis of measurement data obtained by the 2D transversal scanning with the source-receiver system along the sea bottom, over the area of sounded inhomogeneities. The solution begins with the Born approximation, where the original 3D integral equation for the scattered field is reduced to the 1D Fredholm equation of the first kind relative to the depth profile of the lateral spectrum of inhomogeneities. When solving this integral equation for each pair of spectral components, the generalized discrepancy method is in use. Then, corrections to the Born approximation can be obtained in the proposed iterative procedure. For distributed inhomogeneities, the inverse Fourier transform of the retrieved spectrum gives their 3D distribution that can be visualized as tomography images. For solid targets, this spectrum is used to obtain their shape (i.e. to solve the problem of computer holography). Corresponding results of the numerical simulation are presented.

Keywords:

• Acoustics
• near-field diagnostics
• tomography
• holography
• inverse scattering problem

1. Introduction

Acoustic methods of remote diagnostics of the sea environment are now widely applied in multiple applications [1], including the acoustic remote sounding of the seabed [2–17]. Subbottom profilers [18–25] give valuable information for marine technology and engineering, exploration of objects in basic sediments, ecological monitoring in the vicinity of ports, marine pipelines and offshore oil production platforms, and in the analysis of sediments in the marine archaeology. These profiling systems comprise a source of sound towed behind a vessel or firmly mounted to it that produces an acoustic signal (typically, pulse). Some part of the acoustic signal is reflected from the seafloor, whereas the remaining part penetrates the seafloor and is reflected from layers with different acoustic impedances. High-frequency systems (from 2 to 20 kHz) produce a high resolution (tens of cms) of sediment layers below the seafloor. Low-frequency systems (Sparker systems) working in 50 Hz to 4 kHz range have much greater penetration with typical resolution more than 2 m.

It is worth mentioning that the penetration of sound is decreased with frequency, necessitating a variety of subbottom profiling tools specific to different marine environments [25]. The corresponding dependence is not linear; in the range 100–1000 Hz, the penetration depth changes dramatically – from 100 m at 1000 Hz – up to 2000 m at 100 Hz [3]. As a result, various types of subbottom profilers have been developed: high-frequency (>2–3 kHz) boomer, pinger and chirper systems are developed to give very detailed information about near surface features, with bottom penetration of about 100 m. Medium frequency (~1 kHz) profilers such as the sparker system can penetrate to a depth of approximately 500 m, and maintain relatively good resolution. Very low frequency and high energy systems (<1 kHz), such as airgun systems may penetrate a number of kilometres into the seabed, with a corresponding reduction in resolution [25]. However, any one of these systems gives no more than a visualization of the scattered acoustic field. To obtain the distribution of physical parameters of subbottom inhomogeneities, it is necessary to solve corresponding inverse scattering problems.

Solving such a problem is an uneasy task as it demonstrates a history of similar investigations in electrodynamics. The method of scanning tomography proposed in this paper is based on hydrodynamics equations for the sound propagation in layered media [26,27] with low-contrast subbottom inhomogeneities, such as mud-like depositions. Because the statement of this problem for above-mentioned pulse profilers is too complicated, our method involves the analysis of multifrequency measurements of complex amplitudes of harmonic signals obtained by 2D transversal scanning along the bottom in the near-field zone above inhomogeneities. Unlike pulse sounding, where the depth sensitivity is related to the pulse delay, in this case, the frequency dependence of the near-field penetration is employed. In this approach, developed earlier in electromagnetic sounding [28–30], a subwavelength resolution can be realized by measurements of targets in the near-field zone with small-aperture probes. It can be especially important in cases of a strong sound absorption in the sounded medium, when it impossible to use higher frequencies. Based on such measurements computer tomography provides the retrieved 3D distribution of the sound speed perturbation in subbottom inhomogeneities that can be visualized as tomography images. For solid targets with a homogeneous internal structure, the inverse problem consists in calculation of their geometric shape. Its visualization gives us images of such a computer holography.

Within the scope of this approach, the inverse scattering problem is considered on the base of the 3D integral equation of the first kind obtained from the relevant Helmholtz equation. We apply here, the approach that has been developed by us in the near-field subsurface electromagnetic tomography [28–30]. In these papers, 2D distributions of the scattered signal over transversal coordinates are involved in analysis. To provide the depth sensitivity, these data should be obtained in dependence on a third parameter (such as signal frequency, sensor elevation, or its size). Here, we consider a multifrequency method as the most simple in realization.

Such a 3D problem leads to strong limitations of the grid size of discretization used at numerical calculations and, hence, to limitations of the achievable resolution. To overcome this difficulty, we propose to use the scheme of measurements that has been realized in electromagnetic diagnostics [29,30] – the scanning at the condition of the unchanging source-receiver relative position. In analysis, this scheme enables one to make 2D Fourier transform of the corresponding 3D integral equation over transversal coordinates in the Born approximation and, thereby, reduce this problem to the multiple solution of obtained 1D integral equation. It was shown in these papers that for subsurface targets in the near-field zone, the property of small-aperture antennas to generate and receive evanescent waves made it possible to realize a subwavelength resolution of such a tomography that has been demonstrated with the multifrequency microwave scanning system [30]. Note, that it is the only way when it is impossible to use higher frequencies (or short pulses) – technically, or because of large extinction in media.

In this paper, we use this approach to develop a hydroacoustic method of subbottom near-field multifrequency tomography (retrieval of 3D structure of sound speed perturbations in distributed inhomogeneities) and holography (retrieval of the shape of inhomogeneities with sharp boundaries). Possibilities of this new diagnostics are demonstrated in numerical simulation.

2. Theory

In Figure 1, one can see the scheme of proposed near-field tomography based on data of 2D scanning in x-y plane along the bottom in the near zone, over the region of subbottom inhomogeneities.

Figure 1. Scheme of the main elements of the acoustic system for the subbottom tomography in the sea waveguide.

This method can be used in diagnostics of objects in the bottom sediments. In the formulation of this problem, we can use the three-layered model (1 – liquid half-space, such as mud-like sediments; 2 – water (for example, shallow sea); 3 – air). If the bottom represents as a more consolidated environment, a more appropriate seabed model would be an elastic half-space [26]. In the liquid half-space at z > h (medium 1), we assume the density ρ = ρ₁, and the velocity of sound c = c₁. In the water layer (0 < z ≤ h, medium 2) ρ = ρ₂, the velocity of sound c = c₂. The third medium in the half-space z ≤ 0 is air. In general case, observable objects could be heterogeneous both of the speed of sound and density. In view of the complexity of this formulation, in this paper, we consider these problems separately. At first, we consider the case of refractive inhomogeneities, i.e. of a small perturbations c_s in a medium with the constant velocity c = c₀ of sound, i.e. $c=c_0+c_s(\mathbf{r})$, where $c_0=c_1,c_2$, $\mathbf{r}=(x,y,z)$.

From Euler, continuity and state equations for the perfect incompressible liquid, the wave equation for the sound wave pressure P is written in the linear approximation [26] as:

(1) $$\frac{1}{c^2}\frac{{\mathit{\partial }}^{2}P(\mathbf{r},t)}{\mathit{\partial }{t}^2}-{\mathrm{\nabla }}^{2}P(\mathbf{r},t)=f\left(\mathbf{r},t\right),$$(1)

where $P(\mathbf{r},t)$ – pressure at the point $\mathbf{r}$ and time t, c² = ∂P/∂ρ is the velocity of sound, f is the source. The analysis is based on the solution of the Helmholtz equation obtained from (1) for the complex amplitudes $P(\mathbf{r})$ of the pressure wave $P(\mathbf{r},t)=P(\mathbf{r})e^{-i\mathit{\omega }t}$ (ω is cyclic frequency) and of the source $f(\mathbf{r},t)=f(\mathbf{r})e^{-i\mathit{\omega }t}$:

(2) $${\mathrm{\nabla }}^{2}P+k^{2}P=-f(\mathit{\omega },\mathbf{r}),$$(2)

where k = ω/c. Equation (2) can describe the propagation of waves in absorbing media, if we assume $k=\frac{\mathit{\omega }}{c}{[1-\frac{i\mathit{\omega }}{\mathit{\rho }}(\mathit{\varsigma }+\frac{4}{3}\mathit{\eta })]}^{-1/2}$ [26]; here ζ, η – shear and bulk viscosity, respectively. Absorption becomes significant as the frequency increases. For example, absorption in sand and marine sediments at 100 kHz amounts 61 and 26 dB/m respectively [26]. Absorption in water is much smaller, and for distances and frequencies considered in this paper, it can be neglected.

In the presence of inhomogeneity, it is possible to represent parameter k² as $k^2(\mathbf{r})={\mathit{\omega }}^2/c^2(\mathbf{r})=[1-C_s(\mathbf{r})]{\mathit{\omega }}^2/c_0^2$, where $C_s(\mathbf{r})=1-c_0^2/{[c_0+c_s(\mathbf{r})]}^2$ is a dimensionless parameter of inhomogeneity (similarly, for the case of absorbing inhomogeneities $k^2=\frac{{\mathit{\omega }}^2}{c^2}/[1-\frac{i\mathit{\omega }}{\mathit{\rho }}(\mathit{\varsigma }+\frac{4}{3}\mathit{\eta })]=[1-C_s(\mathbf{r})]{\mathit{\omega }}^2/c_0^2$ and $C_s(\mathbf{r})=1-c_0^2/\{{[c_0+c_s(\mathbf{r})]}^2[1-i\mathit{\omega }(\mathit{\varsigma }+\frac{4}{3}\mathit{\eta })/\mathit{\rho }]\}$).

Then, using the Green function formalism, it is possible to obtain the integral expression for the pressure in the considered waveguide:

(3) $$P(\mathbf{r})=P_0(\mathbf{r})+P_1(\mathbf{r})=P_0(\mathbf{r})-\frac{{\mathit{\omega }}^2}{c_1^2}\underset{V}{\int }{C}_s({\mathbf{r}}^{\prime })G(\mathbf{r},{\mathbf{r}}^{\prime })P({\mathbf{r}}^{\prime })\text{d}{\mathbf{r}}^{\prime },$$(3)

where $P_0(\mathbf{r})=\int f(\mathit{\omega },{\mathbf{r}}^{\prime })G(\mathbf{r},{\mathbf{r}}^{\prime })\text{d}{\mathbf{r}}^{\prime }$ is the probing field, $G(\mathbf{r},{\mathbf{r}}^{\prime })$ is the dyadic Green function, $P_1(\mathbf{r})$ is the field scattered by inhomogeneities. In solving the direct problem (to determine the field distribution in a medium with inhomogeneities), Equation (3) is considered as a Fredholm integral equation of the 2nd kind, which has the solution:

(4) $$P(\mathbf{r})=P_0(\mathbf{r})+\underset{V}{\int }R(\mathbf{r},{\mathbf{r}}^{\prime })P_0({\mathbf{r}}^{\prime })\text{d}{\mathbf{r}}^{\prime },$$(4)

where the resolvent R is represented by Neumann series:

(5) $$R=K^{(1)}+K_{}^{(2)}+\cdots +K_{}^{(n)}+\cdots ,$$(5)

$$K^{(1)}(\mathbf{r},{\mathbf{r}}^{\prime })=-\frac{{\mathit{\omega }}^2}{c_1^2}{C}_s({\mathbf{r}}^{\prime })G(\mathbf{r},{\mathbf{r}}^{\prime }),K_{}^{(n)}=\underset{V}{\int }K(\mathbf{r},\mathbf{s})K_{}^{(n-1)}(\mathbf{s},{\mathbf{r}}^{\prime })\text{d}\mathbf{s}.$$

Terms of the series in (5) are interpreted as the contribution of single, double, etc. scattering. Green’s functions can be calculated as the expansion of the probing field in the 2D spectrum of plane waves over lateral (transverse) coordinates with the summation of geometric series of multiple-reflected waves [26,27]. This approach allowed us to obtain explicit expressions for the lateral spectrum of Green functions, which are used in the subsequent analysis (see details in Appendix 1). It is worth mentioning that such Green functions can be derived for an arbitrary multilayer medium.

To solve the considered tomography problem, i.e. to retrieve the 3D distribution of $C_s(\mathbf{r})$ by measurements of the scattered field, one should consider expression for $P_1(\mathbf{r})$ in (3) as a non-linear integral equation of the first kind that can be solved iteratively, beginning with the Born approximation, making the substitution $P({\mathbf{r}}^{\prime })\to P_0(\mathbf{r})$ in (3):(6) $$P_1(\mathbf{r})=-\frac{{\mathit{\omega }}^2}{c_1^2}\underset{V}{\int }{C}_s({\mathbf{r}}^{\prime })G(\mathbf{r},{\mathbf{r}}^{\prime })P_0({\mathbf{r}}^{\prime })\text{d}{\mathbf{r}}^{\prime }=-\frac{{\mathit{\omega }}^2}{c_1^2}\underset{V}{\int }{C}_s({\mathbf{r}}^{\prime })G(\mathbf{r},{\mathbf{r}}^{\prime })\int f(\mathit{\omega },{\mathbf{r}}^{″})G({\mathbf{r}}^{\prime },{\mathbf{r}}^{″})\text{d}{\mathbf{r}}^{″}\text{d}{\mathbf{r}}^{\prime }.$$(6)

It is a Fredholm integral equation of the first kind, so one deals with a well-known ill-posed problem. The solving of such 3D problems is especially difficult. For example, to obtain the solution with the resolution of 100 points for each coordinate, it requires 10¹² numbers in the computer memory to store the 6D kernel of the Equation (6). It leads to inevitable limitations of the discretization and, hence, of the achievable resolution. To overcome this difficulty, we use the approach based on the reduction of this equation to a convolution equation over lateral coordinates because such equations can be reduced to 1D equations and solved relative to each pair of their Fourier components separately. Unfortunately, (6) isn’t a convolution equation, and such a reduction, in general, is impossible. Nevertheless, there are statements of this problem, where this transformation is still possible.

The first one is to regard the product $\tilde{C}({\mathbf{r}}^{\prime })=C_s({\mathbf{r}}^{\prime })P_0({\mathbf{r}}^{\prime })$ as the desired function. However, this desired function is frequency-dependent, and the most promising multifrequency method can’t be based on this approach. The second approach is to use a single plane wave as the probing field $P_0({\mathbf{r}}^{\prime })$ that can be realized in sounding with a remote source operating in a single-mode regime. In this case, the kernel in (6) is formed by the product of the Green function with an exponent that makes it possible to carry out the convolution over the lateral coordinates. Disadvantages of this method are related to frequency limitations of the single-mode operation when the receiver should be placed far away from the source – in the range where the field has a single-mode structure.

The third approach, which has been used in similar electrodynamics problems [29,30], is to use data of the scanning with a source-receiver system at the condition of the unchanging source-receiver relative position. In this case, making the Fourier transform over transversal coordinates in (6), we obtain a 1D integral equation (see details in Appendix 1):(7) $$\begin{array}{cc}{P}_1(k_x,k_y,z,\mathit{\delta }\mathbf{r})=\hfill & -16{\mathit{\pi }}^4\frac{{\mathit{\omega }}^2}{c_1^2}\underset{z^{\prime }}{\int }{C}_s(k_x,k_y,z^{\prime })\int \int e^{-i({\mathit{\kappa }}_x-k_x)\mathit{\delta }x-i({\mathit{\kappa }}_y-k_y)\mathit{\delta }y}\text{d}{\mathit{\kappa }}_x\text{d}{\mathit{\kappa }}_y\hfill \\ \hfill & \times \underset{z^{″}}{\int }f({\mathit{\kappa }}_x-k_x,{\mathit{\kappa }}_y-k_y,z^{″}-z-\mathit{\delta }z)G_{21}({\mathit{\kappa }}_x-k_x,{\mathit{\kappa }}_y-k_y,z^{\prime },z^{″})G_{12}({\mathit{\kappa }}_x,{\mathit{\kappa }}_y,z,z^{\prime })\text{d}{z}^{″}\text{d}{z}^{\prime },\hfill \end{array}$$(7)

where spectra of Green functions $G_{12},G_{21}$ are determined in Appendix 1, $\mathit{\delta }\mathbf{r}=(\mathit{\delta }x,\mathit{\delta }y,\mathit{\delta }z)$ is the offset of the receiver position relative to the position of the source. For brevity, here and below, we use for Fourier transforms the same notation, so values and their transforms are identified by their arguments. The derived 1D Fredholm integral equation of the first kind should be sequentially solved for each pair of spectral components $k_x,k_y$ relative to the depth profile of the lateral spectrum of inhomogeneities C_s(k_x, k_y, z′). From the practical point of view, a great advantage of this measurement scheme is the fact that all variations of the signal (received scattered field P₁) can be related only to subbottom inhomogeneities.

Within frameworks of this approach, it is possible to propose some corrections beyond the Born approximation, though a direct reduction of the problem to a 1D equation is no longer possible. Indeed, if the solution in the Born approximation is obtained, it can be used for the correction of the probing field in subsequent iterations:

(8) $$P_1(\mathbf{r},\mathit{\delta }\mathbf{r})=-\frac{{\mathit{\omega }}^2}{c_0^2}\underset{V}{\int }{C}_1^n({\mathbf{r}}^{\prime })G_{12}(\mathbf{r},{\mathbf{r}}^{\prime })[P_0({\mathbf{r}}^{\prime })+P_1(C_1^{n-1},{\mathbf{r}}^{\prime })]\text{d}{\mathbf{r}}^{\prime }.$$(8)

Performing 2D Fourier transform over the lateral coordinates, obtain the k-space representation of (8):

(9) $$\begin{array}{cc}{P}_1(k_x,k_y,z,\mathit{\delta }\mathbf{r})=\hfill & -\frac{{\mathit{\omega }}^2}{c_0^2}16{\mathit{\pi }}^4\underset{z^{\prime }}{\int }\{C_1^{(n)}(k_x,k_y,z^{\prime })\underset{-\infty }{\overset{\infty }{\int }}\underset{-\infty }{\overset{\infty }{\int }}{e}^{-i({\mathit{\kappa }}_x-k_x)\mathit{\delta }x-i({\mathit{\kappa }}_y-k_y)\mathit{\delta }y}\text{d}{\mathit{\kappa }}_x\text{d}{\mathit{\kappa }}_y\hfill \\ \hfill & \times \underset{z^{″}}{\int }f({\mathit{\kappa }}_x-k_x,{\mathit{\kappa }}_y-k_y,z^{″}-z-\mathit{\delta }z)G_{21}({\mathit{\kappa }}_x-k_x,{\mathit{\kappa }}_y-k_y,z^{\prime },z^{″})G_{12}({\mathit{\kappa }}_x,{\mathit{\kappa }}_y,z,z^{\prime })\text{d}{z}^{″}\hfill \\ \hfill & -\frac{{\mathit{\omega }}^2}{c_0^2}16{\mathit{\pi }}^4\underset{-\infty }{\overset{\infty }{\int }}\underset{-\infty }{\overset{\infty }{\int }}{C}_1^{(n)}({\mathit{\kappa }}_x,{\mathit{\kappa }}_y,z^{\prime })\underset{z^{″}}{\int }{C}_1^{(n-1)}(k_x-{\mathit{\kappa }}_x,k_y-{\mathit{\kappa }}_y,z^{″})\underset{z^{″\prime }}{\int }\underset{-\infty }{\overset{\infty }{\int }}\underset{-\infty }{\overset{\infty }{\int }}{e}^{-i({\mathit{\kappa }}_x+{\mathit{\kappa }}_x^{\prime }-k_x)\mathit{\delta }x-i({\mathit{\kappa }}_y+{\mathit{\kappa }}_y^{\prime }-k_y)\mathit{\delta }y}\hfill \\ \hfill & \times f({\mathit{\kappa }}_x+{\mathit{\kappa }}_x^{\prime }-k_x,{\mathit{\kappa }}_y+{\mathit{\kappa }}_y^{\prime }-k_y,z^{″\prime }-z-\mathit{\delta }z)G_{21}({\mathit{\kappa }}_x+{\mathit{\kappa }}_x^{\prime }-k_x,{\mathit{\kappa }}_y+{\mathit{\kappa }}_y^{\prime }-k_y,z^{″},z^{″\prime })\hfill \\ \hfill & \times G_{11}({\mathit{\kappa }}_x^{\prime },{\mathit{\kappa }}_y^{\prime },z^{\prime },z^{″})G_{12}({\mathit{\kappa }}_x+{\mathit{\kappa }}_x^{\prime },{\mathit{\kappa }}_y+{\mathit{\kappa }}_y^{\prime },z,z^{\prime })\hfill \\ \hfill & \text{d}{\mathit{\kappa }}_x^{\prime }\text{d}{\mathit{\kappa }}_y^{\prime }\text{d}{\mathit{\kappa }}_x\text{d}{\mathit{\kappa }}_y\text{d}{z}^{″\prime }\text{d}{z}^{″}\}\text{d}{z}^{\prime }.\hfill \end{array}$$(9)

As you can see, this equation is 3D – in contrast to (7) obtained in the Born approximation. Nevertheless, it can be used for corrections as a 1D equation, if one assumes $C_1^{(n)}=C_1^{(n-1)}$ in the second term on (9). In such a statement, corrections can be obtained with available computers.

It is suitable to note that for inhomogeneities in the far-field zone the solution of the inverse scattering problem leads to fundamental Rayleigh limitations of the achievable resolution. But in the proposed method, where we consider the scheme with near-field measurements, a subwavelength resolution can be realized. As it was mentioned above, to retrieve the desired 3D distribution, data of 2D scanning over transversal coordinates (x, y) should be measured in dependence on a third parameter that provides the depth sensitivity. For that, it is possible to use the multilevel scanning at a single frequency (i.e. the scanning in x-y plane at several altitudes above the bottom) or multifrequency measurements at a single level. The latter method is simpler in use, and, in this statement, the integral Equation (7) can be rewritten in a compact form as:

(10) $$P_1(k_x,k_y,\mathit{\omega })=\underset{z^{\prime }}{\int }{C}_s(k_x,k_y,z^{\prime })K(k_x,k_y,z^{\prime },\mathit{\omega })\text{d}{z}^{\prime }.$$(10)

In solving this equation, we use the algorithm [29] based on the Tikhonov’s generalized discrepancy principal [32] for complex-valued functions in Hilbert space ${\mathbf{W}}_{\mathbf{2}}^{\mathbf{1}}$ (Sobolev’s space). Detailed comments to Tikhonov’s method in application to various problems of physical diagnostics can be found in the book [31]. Let us represent (10) in the operator form:

(11) $$\mathbf{K}{C}_s=P_1^{\mathit{\delta }},$$(11)

where $P_1^{\mathit{\delta }}$ is a data vector that includes measurement errors. The measure of data errors is determined as

(12) $${\mathit{\delta }}_P^2=sup{∥\mathbf{K}{C}_s-P_1^{\mathit{\delta }}∥}_{L_2}^2=sup\left(\frac{1}{\mathrm{\Delta }\mathit{\omega }}{\int }_{\mathrm{\Delta }\mathit{\omega }}{\left|P_1(\mathit{\omega })-P_1^{\mathit{\delta }}(\mathit{\omega })\right|}^2\text{d}\mathit{\omega }\right)$$(12)

where $P_1(\mathit{\omega })$ in (10) corresponds to the exact solution C_s(z), Δω is the frequency band of analysis, $∥x∥$ is the norm of a function as an element of Hilbert space L₂ (space of complex-valued functions, summable with the square). In frameworks of this method, it is necessary to take into account errors in the kernel of this equation that include discretization errors at the numerical solution as well as errors of the Born approximation:

(13) $${\mathit{\delta }}_h^2=sup{∥{\mathbf{K}}_h{C}_s-\mathbf{K}{C}_s∥}_{L_2}^2=sup{∥{\mathbf{K}}_h{C}_s-P_1∥}_{L_2}^2=h^2||C_s{||}_{W_2^1}^2,$$(13)

where h is a constant, ${\mathbf{K}}_h$ is the numerical approximation of $\mathbf{K}$ in (10),

(14) $${∥C_s∥}_{W_2^1}^2=\frac{1}{\mathrm{\Delta }z}{\int }_{z_{min}}^{z_{max}}\left\{{\left|C_s\right|}^2(z)+{\left|\mathrm{\Delta }z\frac{\mathrm{d}{C}_s}{\mathrm{d}z}\right|}^2\right\}\mathrm{d}z.$$(14)

where Δz = z_max − z_min is the depth region of analysis. Note, that unlike the error δ_P, this measure of kernel errors is proportional to the norm of the desired solution. Both kinds of errors can lead to incompatibility of data with the equation to be solved. The corresponding measure of errors related to this incompatibility δ_μ can’t exceed the sum of two previous error components [32], so it is reasonable to determine it as

(15) $${\mathit{\delta }}_{\mathit{\mu }}^2=inf{∥{\mathbf{K}}_{\mathbf{h}}{C}_s-P_1^{\mathit{\delta }}∥}_{L_2}^2={({\mathit{\delta }}_P+{\mathit{\delta }}_h^{})}^2.$$(15)

In this method, the approximate solution minimizes the smoothing functional

(16) $${\mathbf{M}}^{\mathit{\alpha }}(C_s)={∥{\mathbf{K}}_{\mathbf{h}}{C}_s-P_1^{\mathit{\delta }}∥}_{L_2}^2+\mathit{\alpha }{∥C_s∥}_{W_2^1}^2.$$(16)

In the numerical algorithm, the conjugate gradient method was applied to minimize this functional. The regularization parameter α that determines the extent of the solution smoothing is determined from the 1D equation of generalized discrepancy:

(17) $$\mathit{\rho }(\mathit{\alpha })={∥{\mathbf{K}}_{\mathbf{h}}{C}_s^{\mathit{\alpha }}-P_1^{\mathit{\delta }}∥}_{L_2}^2-{\mathit{\delta }}^2=0.$$(17)

where $C_s^{\mathit{\alpha }}$ is the function that minimizes the functional (16), δ² is the sum of all three sources discrepancy that includes all three error components [32]:

(18) $${\mathit{\delta }}^2={\left({\mathit{\delta }}_P+{\mathit{\delta }}_h\right)}^2+{\mathit{\delta }}_{\mathit{\mu }}^2.$$(18)

The value of the regularization parameter (and the extent of smoothing) is determined by the level of this error parameter. As δ² tends to zero, the approximate solution converges to the exact solution in the norm ${\mathbf{W}}_{\mathbf{2}}^{\mathbf{1}}$ and, according to Sobolev’s including theorem, it converges uniformly in the space C with the maximum modulus as a norm – that is the main advantage of this Tikhonov’s approach. Parameters δ_h and δ_μ can be determined in the numerical simulation. Finally, the desired distribution is determined from the solution of (10) in k-space as:

(19) $$C_s(x,y,z)=\int \int C_s(k_x,k_y,z)exp(i{k}_{x}x+i{k}_{y}y)\text{d}{k}_x\text{d}{k}_y.$$(19)

Then, the distribution of the sound speed inhomogeneity is easily calculated as $c_s(\mathbf{r})=c_0^2[(1-\sqrt{1-C_s(\mathbf{r})})/\sqrt{1-C_s(\mathbf{r})}]$. Based on the numerical simulation, it is possible to obtain the estimation of the solution accuracy in the C-norm that is a stronger estimation when compared to L₂ or mean square estimations. So, to obtain true information about the solution accuracy, it is enough to study the method for typical or expected targets. The numerical simulation includes:

(a)	Generation of real-valued non-correlated random gauss-distributed errors δp(x, y, ω)with rms = σ in each point (x,y) of analysis; transformation of errors in the k-space to obtain complex-valued data errors δp(kx, ky, ω)
(b)	Production of input data in (11): $P_1^{\mathit{\delta }}(k_x,k_y,\mathit{\omega })=P_1(k_x,k_y,\mathit{\omega })+\mathit{\delta }(k_x,k_y,\mathit{\omega })$;
(c)	Solving the inverse problem to retrieve the 3D distribution of inhomogeneity;
(d)	Comparison of initial and retrieved distributions.

In the solving algorithm, the first term of $\mathit{\rho }(\mathit{\alpha })$ in (17) should be minimized up to the value of δ² that includes three components (18). The first is ${\mathit{\delta }}_P^2(k_x,k_y)=sup(\frac{1}{\mathrm{\Delta }\mathit{\omega }}{\int }_{\mathrm{\Delta }\mathit{\omega }}{\left|\mathit{\delta }(k_x,k_y,\mathit{\omega })\right|}^2\text{d}\mathit{\omega })$. It is possible to estimate this parameter as the mean value over multiple realizations of random errors, or to use its calculation for a single realization. However, to build a universal algorithm, it is reasonable to use a universal parameter for each pair (k_x, k_y). We choose this parameter as the mean value over the k-space region of analysis:

(20) $${\mathit{\delta }}_P^2=\frac{1}{\mathrm{\Delta }{k}_x\mathrm{\Delta }{k}_y\mathrm{\Delta }\mathit{\omega }}\underset{\mathrm{\Delta }\mathit{\omega }}{\int }\int \int {\left|\mathit{\delta }(k_x,k_y,\mathit{\omega })\right|}^2\text{d}{k}_x\text{d}{k}_y\text{d}\mathit{\omega }$$(20)

that is practically independent on the realization of random errors. In our simulation, we used the pseudorandom number generator with the fixed rms = σ in every random sequence, so ${\mathit{\sigma }}^2=\frac{1}{\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }\mathit{\omega }}\underset{\mathrm{\Delta }\mathit{\omega }}{\int }\int \int {\left|\mathit{\delta }(x,y,\mathit{\omega })\right|}^2\text{d}x\text{d}x\text{d}\mathit{\omega }$. Using the Plancherel theorem, one obtains ${\mathit{\delta }}_P^2=4{\mathit{\pi }}^2{\mathit{\sigma }}^2$. In real measurements, it would be reasonable to estimate these noise parameters by data of measurements in a region without subsurface inhomogeneities – outside the studied place.

In our case, when we know the exact formula for the kernel K of the integral equation, the second parameter ${\mathit{\delta }}_h^2$ is related to the level of discretization, and we minimize this parameter gradually increasing the discretization level to satisfy the condition ${\mathit{\delta }}_h^2<<{\mathit{\delta }}_p^2$. The necessary level of discretization depends on target parameters (its permittivity, sizes, shape, depth). At a scarce discretization, retrieved images can be smoothed, so we used the numerical simulation to satisfy the mentioned condition for our test targets. As a result, it was obtained that the discretization 80 × 80 × 80 in the region of analysis 40 × 40 × 20 m³ gives about the same results as that at 40 × 40 × 40 m³. According to [32], the third part of δ² (measure of incompatibility) satisfies ${\mathit{\delta }}_{\mathit{\mu }}^2\le {({\mathit{\delta }}_p+{\mathit{\delta }}_h)}^2$. So, neglecting the contribution of δ_h, we obtain the estimation ${\mathit{\delta }}^2=2{\mathit{\delta }}_p^2$ that has been used in analysis.

The frequency range of analysis has been determined so that the near-field distribution changed significantly over a sounded object. This condition provides the depth sensitivity of measurements. The number of frequencies has been chosen so that the difference of signal variations at neighbour frequencies exceeds the error level. It is clear that this number increases with the decrease of errors.

Here it should be noted that the inverse problem is simplified for quite a typical case of solid targets composed of a homogeneous material – in this case it is enough to determine the shape of their boundaries. However, it is impossible to obtain the exact positions of boundaries for such targets from (19) – even neglecting retrieval errors – because it is well-known that step-functions can’t be expressed exactly nearby points of gradient discontinuity by arbitrary long Fourier series that represent (19). In such points, the Dini test for pointwise convergence of Fourier series is unfilled, and results are distorted by the Gibbs effect. So, if it is known, a priori that we deal with such a target, the problem can be solved rigorously using our approach.

When it is known a priori that $C_s=C_s^0=const$, the holography problem (i.e. the problem of the target shape retrieval) can be solved using k-space solution C_s(κ_x, κ_y, z) of (7). For that, let us represent the shape of a simple connected target by two single-valued functions x₁(y, z), x₂(y, z) (see in Figure 2).

Figure 2. Boundaries of a target in the section z = const represented by functions x₁(y, z), x₂(y, z).

In this case, retrieved k-space spectrum of a target can be written as a Fourier transform with finite limits:

(21) $$C_s(k_x,k_y,z)=\frac{1}{4{\mathit{\pi }}^2}\underset{y_1}{\overset{y_2}{\int }}\underset{x_1(y)}{\overset{x_2(y)}{\int }}{C}_s^0{e}^{-i{k}_{x}x-i{k}_{y}y}\text{d}x\text{d}y=\frac{C_s^0}{4{\mathit{\pi }}^2}\underset{y_1}{\overset{y_2}{\int }}exp(-i{k}_{y}y)\frac{1}{i{k}_x}(e^{-i{k}_x{x}_1(y)}-e^{-i{k}_x{x}_2(y)})\text{d}y.$$(21)

Then, making the inverse Fourier transform of (21) over k_y

(22) $$C_s(k_x,y^{\prime },z)=\underset{-\infty }{\overset{\infty }{\int }}{C}_s(k_x,k_y,z)exp(i{k}_y{y}^{\prime })\text{d}{k}_y,$$(22)

obtain (after change y′ → y) the complex-valued transcendent equation:

(23) $$C_s(k_x,y,z)=\frac{C_s^0}{2\mathit{\pi }i{k}_x}(e^{-i{k}_x{x}_1(y,z)}-e^{-i{k}_x{x}_2(y,z)})$$(23)

that is equivalent to the system of two real-valued equations. Then, the desired shape of the target expressed by two functions x₁(y, z), x₂(y, z) is obtained from the solution of (23). It should be mentioned that this equation is overdetermined: it can be solved at each value of k_x that makes it possible to optimize the solution. This optimization could be based on the optimal choice of k_x or on the solving (23) as a system of equations with different values of k_x. Depending on available conditions, the value of $C_s^0$ may be known a priori, or obtained from tomography results (19), or from the solution of (23) as an expanded system with some extra equations with different values of k_x. It is clear that similar equation can be written for the shape expressed by functions y₁(x,z), y₂(x,z).

3. Numerical simulation

In numerical simulation, we consider the case of the point source $f(\mathit{\omega },{\mathbf{r}}^{\prime })=\mathit{\delta }({\mathbf{r}}_{\mathbf{0}}-{\mathbf{r}}^{\prime })$ of the probing field because small aperture sources (that should be used in near-field sounding) have about the same k-space spectrum of the emitted field as the point source. The probing field of the point source is expressed by the Green function $P_0(\mathbf{r})=\int f(\mathit{\omega },{\mathbf{r}}^{\prime })G(\mathbf{r},{\mathbf{r}}^{\prime })\text{d}{\mathbf{r}}^{\prime }=G(\mathbf{r},{\mathbf{r}}_0)$. In this case, the transversal spectrum (7) of the scattered field measured in the point $\mathbf{r}={\mathbf{r}}_{\mathbf{0}}-\mathit{\delta }\mathbf{r}$ is expressed as:

(24) $$\begin{array}{cc}{P}_1(k_x,k_y,z,\mathit{\delta }\mathbf{r})=\hfill & -4{\mathit{\pi }}^2{f}_0\frac{{\mathit{\omega }}^2}{c_0^2}\underset{z^{\prime }}{\int }{C}_s(k_x,k_y,z^{\prime })\int \int e^{-i({\mathit{\kappa }}_x-k_x)\mathit{\delta }x-i({\mathit{\kappa }}_y-k_y)\mathit{\delta }y}\text{d}{\mathit{\kappa }}_x\text{d}{\mathit{\kappa }}_y\hfill \\ \hfill & \times G_{21}({\mathit{\kappa }}_x-k_x,{\mathit{\kappa }}_y-k_y,z^{\prime },z+\mathit{\delta }z)G_{12}({\mathit{\kappa }}_x,{\mathit{\kappa }}_y,z,z^{\prime })\text{d}{z}^{\prime },\hfill \end{array}$$(24)

In our simulation, we consider a waveguide water layer with depth h = 20 m and parameters ${\mathit{\rho }}_2^{}=1{\text{g/cm}}^3,c_2=1470\text{m/s}$ above the liquid half-space with ${\mathit{\rho }}_1^{}=1.92{\text{g/cm}}^3,c_1\text{= 1675 ; m/s}$ that simulates mud bottom deposits (see in Figure 1). In Figure 3, distributions of the probing field amplitude (relative to its maximum value P_max) are shown in the vertical section at two frequencies in the band of analysis.

Figure 3. Distributions of the probing field $\text{Re}{P}_0(x,z)/P_{\text{max}}$ for the source at $x_0=12.5\text{m},y_0=0,z_0=19\text{m}$ in the vertical section $y=y_0+0.1\text{m}$ at frequencies 100 (left) and 200 Hz (right).

From these images of two field distributions, it is seen that the near zone penetrates much deeper into the bottom for the source at the frequency 100 Hz than for the source at 200 Hz. This frequency dependence of the near field penetration explains the depth sensitivity of multifrequency sounding of targets in the near zone. It is also seen that in both cases, the range of penetration of the strong field is much shorter than the wavelength. This property provides a subwavelength sensitivity of near-field sounding.

The tomography method for distributed inhomogeneities is simulated with two 3D gauss functions $C_s(x,y,z)=C_{m}exp[-(\frac{{(x-x_m)}^2}{{\mathit{\sigma }}_x^2}+\frac{{(y-y_m)}^2}{{\mathit{\sigma }}_y^2}+\frac{{(z-z_m)}^2}{{\mathit{\sigma }}_z^2})]$ with parameters: C_m1 = C_m2 = 0.25, z_m = 30 m, σ_x1 = 7.7 m, σ_x2 = 13.7 m, σ_y1 = 13.0 m, σ_y2 = 6.0 m, σ_z1 = 6.0 m, σ_z2 = 3.0 m at 12 m distance between centres. The source was located at z₀ = 19.9 m; the receiver – at z_r = 19 m (at 0.1 m and 1 m above the media interface, respectively). Data of ‘measurements’ include calculated complex amplitudes of the scattered field (with added random errors) that simulate results of the scanning over lateral coordinates in the range $40\times 40{\text{m}}^2$ at five frequencies in the region 33–165 Hz. In Figure 4, corresponding images of calculated 2D distributions of $\left|P_1(x,y,z=z_r)\right|/P_{\text{max}}$ (the relation of the scattered field amplitude to its maximum value) are demonstrated.

Figure 4. Distributions of the scattered field $\left|P_1(x,y,z=z_r)\right|/P_{\text{max}}$ at five frequencies used in analysis: f₁ = 33 Hz, f₂ = 66 Hz, f₃ = 99 Hz, f₄ = 132 Hz, and f₅ = 165 Hz.

One can see that the image of the scattered field is rather smoothed at the lowest frequency f₁. Then, at the frequency f₂, it resembles simulated targets; at higher frequencies the structure of images becomes more complicated. To simulate the described tomography method, ‘data’ at five frequencies are used to solve the inverse problem in the k-space (24); then its solution C_s(k_x, k_y, z) has been transferred in (19) to C_s(x, y, z) (in Cartesian coordinates). In Figure 5, results of the tomography simulation are demonstrated for simulated gauss inhomogeneities in the vertical section C_s(x, y = 0, z) at the level of random errors $\mathit{\delta }/∥P_1∥$ of 7%.

Figure 5. Measurements scheme and simulation results of multifrequency tomography (vertical section). Initial distribution of inhomogeneities is shown in the left part of bottom layer, tomogram – in the right part.

In Figure 6, tomography results obtained at two values of simulated errors are shown in the horizontal section at the depth of maxima (at z_m = 30 m).

Figure 6. Simulation of tomography (horizontal section in the range $40\times 40{\text{m}}^2$).

Note: (a) initial distribution C_s(x, y, z = 30 m); (b) tomogram retrieved at error level 7%; (c) tomogram retrieved at error level 35%.

It is seen in Figures 5 and 6 that the applied algorithm works quite good. Results demonstrate that the left inhomogeneity (broader in z-direction) is retrieved better than the right one because corresponding small-scale near-field components of the probing field (that provided the depth resolution) are significantly faded at depth on maxima. The comparison of tomogram details to the shortest wavelength ${\mathit{\lambda }}_{\text{min}}$ in analysis (shown by the arrow in Figure 5) demonstrates the subwavelength resolution of proposed method. In Figures 5 and 6, it is possible to see a good retrieval quality at the data error level of 7%, however, the general shape of inhomogeneities is still seen even at errors of 35% (in Figure 6). At the error of 7%, the value of C_s in the left maximum amounts about 0.8 C_m. It is seen that retrieval errors decrease with simulated data errors; however, we have obtained also that the further convergence continues only at the simultaneous increase of the number of frequencies in analysis that should be chosen in a broader frequency band.

The case of solid objects with sharp boundaries is especially interesting because it seems as more typical. In this case, the solution of the inverse problem in k-space obtained from integral Equation (24) can be used to determine also the shape of such targets from the obtained non-linear Equation (23). In Figure 7, one can see results of the numerical simulation, where the k-space distribution C_s(k_x, k_y, z′) has been retrieved from the solution of (24) for such an object (parallelepiped with C_s = const = 0.25 with sizes ∆x = 5 m, ∆y = 12 m, ∆z = 6 m at z = 25 m). Results are given for the relation of real part of the complex-valued spectrum to its maximal value C_Smax at z = 25 m (the depth position of the rectangular target).

Figure 7. Numerical simulation of k-space tomography method for a solid object: retrieval of the transversal spectrum of the parallelepiped target from the solution of (24) at the error level of 2%. Left, transversal spectrum ($\text{Re}{C}_{\text{s}}(k_{\text{x}},k_{\text{y}},z=25\text{m})/C_{S\text{max}}$) of simulated target shown in the spectral range $0\le k_x\le 2{\text{m}}^{-1},0\le k_y\le 2{\text{m}}^{-1}$; right, its retrieved spectrum.

Comparing initial and retrieved images in Figure 7, one can see the influence of random errors on the retrieval quality of $C_s(k_x,k_y,z=25\text{m})$ that is used in the further analysis to obtain holography images demonstrated below in Figures 8 and 9. The error level in the integral metric L₂ amounts about 10%, so we have some amplification of the error’s level when compared to that of input data. It is a common property of ill-posed problems [31]. Using k-space distributions, shown in Figure 7, holography images of the parallelepiped target have been obtained from the Equation (23) for this target at z = 25 m and, similarly, for the target at z = 30 m. In Figure 8, the shape of these targets is visualized as two functions y₁(x, z), y₂(x, z).

Figure 8. Numerical simulation of holography method for a parallelepiped object.

Notes: Top row: target at z = 25 m. Left, images of the target and its retrieval (part described by function y₁(x, z); right, target and its retrieval described by function y₂(x, z). Bottom row: target at z = 30 m.

Figure 9. Numerical simulation of holography method for a parallelepiped object (another viewpoint). Top row: target at z = 25 m.

Notes: Left, images of target and its retrieval described by function x₁(y, z); right, target and its retrieval described by function x₂(y, z). Bottom row: target at z = 30 m.

In Figure 9, holography images of these targets are shown from a different viewpoint.

Results demonstrate a good quality of the retrieval – one can easily identify the target’s shape and determine its sizes and depth position with reasonable accuracy. Based on such numerical simulations, we found a way for the optimization of the solution of (23) using the freedom in the choice of the k_x value. It was obtained that the optimal value is k_x ≅ 2π/L, where L is a typical scale of a target. At the optimal value, one obtains the most sharp holography images, so the process of optimization is similar with the optical focusing.

In Figure 10, possibilities of the method are demonstrated for spherical targets.

Figure 10. Numerical simulation of holography method for spherical objects.

Notes: (a) spherical target with radius r = 3 m at z = 30 m; (b) spherical target with r = 6 m at z = 30 m; (c) spherical target with r = 6 m at z = 27 m with a spherical cavity; (d) spherical target with r = 6 m at z = 30 m in the absorbing medium. Here, to place images described by functions x₁(y, z), x₂(y, z) in the single picture, these images have been shifted to the left and to the right from the central position of targets.

Results demonstrate somewhat better quality for spherical targets in Figure 10(a) and (b) when compared to the case of parallelepiped targets in Figures 8 and 9. Also, the smoothed spherical surface seems less distorted than sharp-edged shapes. In Figure 10(c), the retrieved image of the spherical surface with a cavity is shown. The retrieved image of this cavity is somewhat distorted that demonstrates a typical property of ill-posed problems that it is more difficult to retrieve more complicated distributions [31].

In Figure 10(d), one can see retrieved images in an absorbing medium with $k_1=\frac{\mathit{\omega }}{c}{[1-0.01i\mathit{\omega }]}^{-1/2}$ for the same target as in Figure 10(b). At that, the skin-depth $d=c/[\mathit{\omega }\text{Im}(k_1)]$ amounts 7.5 m at the highest frequency 165 Hz; at the lowest frequency 33 Hz this medium is almost transparent. Comparing these images with results for the medium without absorption in Figure 10(b), one can see a smoothing effect related to the signal fading.

It is important to note that similar results have been obtained for tomography and holography in waveguides with depth from 20 cm up to 200 m if all sizes (depths of targets, their sizes, wavelengths of analysis) change proportionally. So, for waveguide depth h = 20 cm, similar results can be obtained using measurements in 3.3–16.5 kHz frequency band.

Finally, we note a possibility to modify this method that is important to practical applications. As it was found in similar electromagnetic diagnostics [30], it appeared difficult to discern sounded subsurface objects in images of the scattered signal at separate frequencies (such, as shown in Figure 4) against the noise produced by the surface scattering. Nevertheless, it was found that it is possible to obtain much better images of subsurface targets, using the transformation of multifrequency data to the synthesized pseudopulse.

(25) $$P_1(x_r,y_r,t)={\int }_{\mathrm{\Delta }\mathit{\omega }}^{}{P}_1(x_r,y_r,\mathit{\omega })exp(i\mathit{\omega }t)\text{d}\mathit{\omega }.$$(25)

It can be represented in dependence on the effective depth parameter z_s as P₁(x, y, z_s) = P₁(x, y, t = −2z_s/c₁) (taking into account the sound velocity in a medium and the signal path to and from a scattering element). At that, subsurface inhomogeneities could be clearly seen in visualised 2D images of $\left|P_1(x,y,z_s)\right|$. The strong maximum of $\left|P_1(z_s)\right|$ observed in every point of the scanning region will mark the value of z_s that corresponds to the surface scattering. The same transformation $P_1(k_x,k_y,t)={\int }_{\mathrm{\Delta }\mathit{\omega }}^{}{P}_1(k_x,k_y,\mathit{\omega })exp(i\mathit{\omega }t)\text{d}\mathit{\omega }$ can be applied to the Equation (10) (compact form of (7)):

(26) $$P_1(k_x,k_y,z_s)=\underset{z^{\prime }}{\int }{C}_s(k_x,k_y,z^{\prime })K(k_x,k_y,z^{\prime },z_s)\text{d}{z}^{\prime },$$(26)

$$K(k_x,k_y,z^{\prime },z_s)={\int }_0^{\infty }K(k_x,k_y,z^{\prime },\mathit{\omega })exp(-i\mathit{\omega }2z_s/c_1)\text{d}\mathit{\omega },$$

where P₁(k_x, k_y, z_s) = P₁(k_x, k_y, t = −2z_s/c₁). The same method, as above, can be applied to solve this equation.

In the case of diagnostics that uses real pulse measurements, one has the problem related to calculations of necessary reflection coefficients and Green functions that are, in general, unknown for time-dependent signals. It is possible overcome this problem by the Fourier transformation of the signal to frequency domain. Then, one can apply the theory for multifrequency diagnostics, developed here.

4. Tomography and holography problems for inhomogeneities of density

In the case of inhomogeneities of density, the wave equation for complex amplitudes of a harmonic signal is written [26] as:

(27) $$\mathrm{\Delta }P+k^{2}P-grad(ln\mathit{\rho })gradP=0.$$(27)

It is possible to apply the same method, as above, considering the inhomogeneity as an effective source of sound. Then, using the Green function for homogeneous media, instead of (6), obtain the generalized equation for the scattered field in the Born approximation:

(28) $$P_1(\mathbf{r})=-\frac{{\mathit{\omega }}^2}{c_0^2}\underset{V}{\int }\left[C_s({\mathbf{r}}^{\prime })+\frac{1}{{\mathit{\rho }}_0}grad{\mathit{\rho }}_s({\mathbf{r}}^{\prime })grad{P}_0({\mathbf{r}}^{\prime })\right]G_{12}(\mathbf{r},{\mathbf{r}}^{\prime })P_0({\mathbf{r}}^{\prime })\text{d}{\mathbf{r}}^{\prime }.$$(28)

where ${\mathit{\rho }}_s({\mathbf{r}}^{\prime })$ is the perturbation of density $\mathit{\rho }({\mathbf{r}}^{\prime })={\mathit{\rho }}_0+{\mathit{\rho }}_s({\mathbf{r}}^{\prime })$. It is possible to take by parts the term related to the perturbation:

(29) $$\begin{array}{cc}{P}_1(\mathbf{r})=\hfill & -\frac{{\mathit{\omega }}^2}{c_0^2}\underset{V}{\int }\{C_s({\mathbf{r}}^{\prime })G_{12}(\mathbf{r},{\mathbf{r}}^{\prime })P_0({\mathbf{r}}^{\prime })\hfill \\ \hfill & -{\mathit{\rho }}_s({\mathbf{r}}^{\prime })\frac{1}{{\mathit{\rho }}_0}([\frac{\mathit{\partial }{P}_0({\mathbf{r}}^{\prime })}{\mathit{\partial }{x}^{\prime }}{G}_{12}(\mathbf{r},{\mathbf{r}}^{\prime })P_0({\mathbf{r}}^{\prime }){]}_x^{\prime }+[\frac{\mathit{\partial }{P}_0({\mathbf{r}}^{\prime })}{\mathit{\partial }{y}^{\prime }}{G}_{12}(\mathbf{r},{\mathbf{r}}^{\prime })P_0({\mathbf{r}}^{\prime }){]}_y^{\prime }+[\frac{\mathit{\partial }{P}_0({\mathbf{r}}^{\prime })}{\mathit{\partial }{z}^{\prime }}{G}_{12}(\mathbf{r},{\mathbf{r}}^{\prime })P_0({\mathbf{r}}^{\prime }){]}_z^{\prime })\}\text{d}{\mathbf{r}}^{\prime }.\hfill \end{array}$$(29)

Of course, it is difficult to solve problem for inhomogeneities of sound and density simultaneously. However, it is still possible when these inhomogeneities are not independent, for example, when we know the target material or components of a mixture. This is just the proper condition for the above-considered holography diagnostics of solid targets. For the scheme of sounding with the fixed source-receiver system, (29) is a convolution equation over transversal coordinates that it is possible to solve in the same way as (6) – transforming it to the 1D integral equation in k-space.

Let us assume for a solid target $C_s({\mathbf{r}}^{\prime })=s{\mathit{\rho }}_s(\mathbf{r})$, where the material parameter s is a constant. Then it is possible to rewrite (29) as

(30) $$\begin{array}{cc}{P}_1(\mathbf{r},\mathit{\omega })\hfill & =-\frac{{\mathit{\omega }}^2}{c_0^2}\underset{V}{\int }{\mathit{\rho }}_s({\mathbf{r}}^{\prime })\{s{G}_{12}(\mathbf{r},{\mathbf{r}}^{\prime })P_0({\mathbf{r}}^{\prime })-\frac{1}{{\mathit{\rho }}_0}([\frac{\mathit{\partial }{P}_0({\mathbf{r}}^{\prime })}{\mathit{\partial }{x}^{\prime }}{G}_{12}(\mathbf{r},{\mathbf{r}}^{\prime })P_0({\mathbf{r}}^{\prime }){]}_x^{\prime }+[\frac{\mathit{\partial }{P}_0({\mathbf{r}}^{\prime })}{\mathit{\partial }{y}^{\prime }}{G}_{12}(\mathbf{r},{\mathbf{r}}^{\prime })P_0({\mathbf{r}}^{\prime }){]}_y^{\prime }\hfill \\ \hfill & +[\frac{\mathit{\partial }{P}_0({\mathbf{r}}^{\prime })}{\mathit{\partial }{z}^{\prime }}{G}_{12}(\mathbf{r},{\mathbf{r}}^{\prime })P_0({\mathbf{r}}^{\prime }){]}_z^{\prime })\}{\mathbf{dr}}^{\prime }\hfill \\ \hfill & =\underset{V}{\int }{\mathit{\rho }}_s(x^{\prime },y^{\prime },z^{\prime })K(x-x^{\prime },y-y^{\prime },z,z^{\prime },\mathit{\omega })\text{d}{x}^{\prime }\text{d}{y}^{\prime }\text{d}{z}^{\prime }.\hfill \end{array}$$(30)

The Fourier transform over transversal coordinates x and y reduces (30) to a 1D integral equation of the corresponding inverse problem in the k-space:

(31) $$P_1(k_x,k_y,z,\mathit{\omega })=\underset{z}{\int }{\mathit{\rho }}_s(k_x,k_y,z^{\prime })K(k_x,k_y,z,z^{\prime },\mathit{\omega })\text{d}{z}^{\prime }.$$(31)

Then, in the same way as in (21)–(23), obtain the equation of holography that determines the shape of a target by two functions x₁(y, z) and x₂(y, z):

(32) $${\mathit{\rho }}_s(k_x,y,z)=\frac{{\mathit{\rho }}_s^0}{2\mathit{\pi }i{k}_x}(e^{-i{k}_x{x}_1(y,z)}-e^{-i{k}_x{x}_2(y,z)}),$$(32)

where ${\mathit{\rho }}_s^0$ is the a priori known density of the target, ${\mathit{\rho }}_s(k_x,y^{\prime },z)={\int }_{-\infty }^{\infty }{\mathit{\rho }}_s(k_x,k_y,z)exp(i{k}_y{y}^{\prime })\text{d}{k}_y$. Therefore, holography results demonstrated in Figures 7–10 are valid for the proposed near-field diagnostics even in this, more generous case.

5. Conclusion

A new method of near-field multifrequency tomography and holography of subbottom inhomogeneities is developed and studied in the numerical simulation. The developed approaches can be applied to diagnostics of low-contrast under-bottom targets, such as, for example, archaeology artefacts. Similar near-field approaches can be developed in sound defectoscopy and diagnostics of living tissues.

Funding

This work was supported by the Russian Foundation for Basic Research [grant number 17-07-00488]; partly by grant of Ministry of education and science of Russian Federation under Agreement on August 27, 2013 No. 02.B.49.21.0003 between Ministry of education and science and Lobachevsky State University of Nizhny Novgorod and within frameworks of the state plan in the region of scientific activity (the project part code No. 1727).

Disclosure statement

No potential conflict of interest was reported by the authors.

Appendix 1

Direct problem

To solve the direct problem, it is necessary to find the Green function $G(\mathbf{r},{\mathbf{r}}^{\prime })$ in the Equation (2) with taking into account, the boundary conditions in the waveguide, i.e. to solve the problem for the point source located in the point $\mathbf{r}={\mathbf{r}}_{\mathbf{0}}$:

(1A) $$\mathrm{\Delta }P+k^{2}P=-\mathit{\delta }(\mathbf{r}-{\mathbf{r}}_{\mathbf{0}}),$$(1A)

Then, the solution of (2) with the source ($-f(\mathit{\omega },\mathbf{r})$) is given by

(2A) $$P_0(\mathbf{r})=\int f(\mathit{\omega },{\mathbf{r}}^{\prime })G(\mathbf{r},{\mathbf{r}}^{\prime })\text{d}{\mathbf{r}}^{\prime }.$$(2A)

Green functions in any multilayer medium are obtained from the Green function in a homogeneous medium that is expressed by the spherical wave:

(3A) $$G_0(\mathbf{r},{\mathbf{r}}_{\mathbf{0}})=\frac{1}{4\mathit{\pi }}\frac{exp(ik\left|\mathbf{r}-{\mathbf{r}}_{\mathbf{0}}\right|)}{\left|\mathbf{r}-{\mathbf{r}}_{\mathbf{0}}\right|}.$$(3A)

Following [26,27], we use the plane wave expansion of (3A):

(4A) $$G_0(\mathbf{r},{\mathbf{r}}_{\mathbf{0}})=\frac{i}{8{\mathit{\pi }}^2}\int \int e^{i[{\mathit{\kappa }}_x(x-x_0)+{\mathit{\kappa }}_y(y-y_0)+\sqrt{k^2-{\mathit{\kappa }}_x^2-{\mathit{\kappa }}_y^2}\left|z-z_0\right|]}\frac{\text{d}{\mathit{\kappa }}_x\text{d}{\mathit{\kappa }}_y}{\sqrt{k^2-{\mathit{\kappa }}_x^2-{\mathit{\kappa }}_y^2}},$$(4A)

This Green function can be considered as the 2D inverse Fourier transform over differences $(x-x_0),(y-y_0)$ of its lateral spectrum:

(5A) $$G_0({\mathit{\kappa }}_x,{\mathit{\kappa }}_y)=\frac{i}{8{\mathit{\pi }}^2}\frac{e^{i\sqrt{k^2-{\mathit{\kappa }}_x^2-{\mathit{\kappa }}_y^2}\left|z-z_0\right|}}{\sqrt{k^2-{\mathit{\kappa }}_x^2-{\mathit{\kappa }}_y^2}}=\frac{i}{8{\mathit{\pi }}^2}\frac{e^{i{k}_z\left|z-z_0\right|}}{k_z}.$$(5A)

To obtain Green functions in a multilayer medium, it is necessary to match fields on boundaries of layers taking into account boundary conditions. For that, it is possible to express the field in each layer as a sum of multiply reflected plane waves [26,27]. For the considered problem of tomography in the three-layer medium (Pekeris waveguide), it is necessary to obtain fields in the medium 1 from the source located in the medium 2 and fields in the medium 2 from the source located in the medium 1, i.e. Green functions G₂₁, G₁₂:$$\begin{array}{cc}{G}_{21}(\mathbf{r},{\mathbf{r}}_{\mathbf{0}})\hfill & =\frac{i}{8{\mathit{\pi }}^2}\int \int e^{i[{\mathit{\kappa }}_x(x-x_0)+{\mathit{\kappa }}_y(y-y_0)]}{T}_{21}{e}^{i{k}_{1z}(z-h)}(e^{i{k}_{2z}(h-z_0)}+R_{23}{e}^{i{k}_{2z}(d+z_0)})\times \sum _{l=0}^{\infty }{(R_{23}{R}_{21})}^l{e}^{2i{k}_{2z}hl}\frac{\text{d}{\mathit{\kappa }}_x\text{d}{\mathit{\kappa }}_y}{k_{2z}}=\hfill \\ \hfill & =\frac{i}{8{\mathit{\pi }}^2}\int \int T_{21}{e}^{i[{\mathit{\kappa }}_x(x-x_0)+{\mathit{\kappa }}_y(y-y_0)+k_{1z}(z-h)]}\frac{e^{i{k}_{2z}(h-z_0)}+R_{23}{e}^{i{k}_{2z}(h+z_0)}}{1-R_{23}{R}_{21}{e}^{2i{k}_{2z}h}}\frac{\text{d}{\mathit{\kappa }}_x\text{d}{\mathit{\kappa }}_y}{k_{2z}},\hfill \end{array}$$

(6A) $$G_{21}(\mathbf{r},{\mathbf{r}}_{\mathbf{0}})=\frac{i}{8{\mathit{\pi }}^2}\int \int T_{21}\frac{{\text{e}}^{i{\mathit{\kappa }}_{2z}(h-z_0)}+R_{23}{\text{e}}^{i{\mathit{\kappa }}_{2z}(h+z_0)}}{(1-R_{23}{R}_{21}{\text{e}}^{2i{\mathit{\kappa }}_{2z}h})}{e}^{i[{\mathit{\kappa }}_x(x-x_0)+{\mathit{\kappa }}_y(y-y_0)+{\mathit{\kappa }}_{1z}(z-h)]}\frac{\text{d}{\mathit{\kappa }}_x\text{d}{\mathit{\kappa }}_y}{\sqrt{k_2^2-{\mathit{\kappa }}_{\perp }^2}},$$(6A) $$\begin{array}{cc}{G}_{12}(\mathbf{r},{\mathbf{r}}_{\mathbf{0}})\hfill & =\frac{i}{8{\mathit{\pi }}^2}\int \int e^{i[{\mathit{\kappa }}_x(x-x_0)+{\mathit{\kappa }}_y(y-y_0)]}{T}_{12}{e}^{i{k}_{1z}(z_0-h)}(e^{i{k}_{2z}(h-z)}+R_{23}{e}^{i{k}_{2z}(h+z)})\times \sum _{l=0}^{\infty }{(R_{23}{R}_{21})}^l{e}^{2i{k}_{2z}hl}\frac{\text{d}{\mathit{\kappa }}_x\text{d}{\mathit{\kappa }}_y}{k_{1z}}\hfill \\ \hfill & =\frac{i}{8{\mathit{\pi }}^2}\int \int T_{21}{e}^{i[{\mathit{\kappa }}_x(x-x_0)+{\mathit{\kappa }}_y(y-y_0)+k_{1z}(z_0-h)]}\frac{e^{i{k}_{2z}(h-z)}+R_{23}{e}^{i{k}_{2z}(h+z)}}{1-R_{23}{R}_{21}{e}^{2i{k}_{2z}h}}\frac{\text{d}{\mathit{\kappa }}_x\text{d}{\mathit{\kappa }}_y}{k_{1z}},\hfill \end{array}$$

(7A) $$G_{12}(\mathbf{r},{\mathbf{r}}_{\mathbf{0}})=\frac{i}{8{\mathit{\pi }}^2}\int \int T_{12}\frac{e^{-i{k}_{2z}(z-h)}+R_{23}{e}^{i{k}_{2z}(z+h)}}{1-R_{23}{R}_{21}{e}^{2i{k}_{2z}h}}{e}^{i[{\mathit{\kappa }}_x(x-x_0)+{\mathit{\kappa }}_y(y-y_0)+\sqrt{k_1^2-{\mathit{\kappa }}_{\perp }^2}(z_0-h)]}\frac{\text{d}{\mathit{\kappa }}_x\text{d}{\mathit{\kappa }}_y}{\sqrt{k_1^2-{\mathit{\kappa }}_{\perp }^2}},$$(7A)

where ${\mathit{\kappa }}_{1z}=\sqrt{k_1^2-{\mathit{\kappa }}_{\perp }^2},{\mathit{\kappa }}_{2z}=\sqrt{k_2^2-{\mathit{\kappa }}_{\perp }^2}$, and reflection and transmission coefficients $R_{ij},T_{ij}$ from ith layer to jth layer are determined by the method of input impedances in [26,27] as(8A) $$R_{ij}=\frac{m_{ij}\sqrt{k_i^2-{\mathit{\kappa }}_{\perp }^2}-\sqrt{n_{ij}^2{k}_i^2-{\mathit{\kappa }}_{\perp }^2}}{m_{ij}\sqrt{k_i^2-{\mathit{\kappa }}_{\perp }^2}+\sqrt{n_{ij}^2{k}_i^2-{\mathit{\kappa }}_{\perp }^2}},T_{ij}=1+R_{ij},$$(8A)

where $m_{ij}={\mathit{\rho }}_j/{\mathit{\rho }}_i,n_{ij}=c_i/c_j$. In the case considered in this paper, the reflection coefficient on the water–air interface is R₂₃ = −1. Expressions (6A), (7A) are also expressed as inverse Fourier transforms of their transversal spectra, used in formulas (7), (24) as:

(9A) $$G_{21}({\mathit{\kappa }}_x,{\mathit{\kappa }}_y,z,z^{\prime })=\frac{i}{8{\mathit{\pi }}^2}{T}_{21}\frac{e^{i{k}_{2z}(h-z^{\prime })}-e^{i{k}_{2z}(h+z^{\prime })}}{(1+R_{21}{e}^{2i{k}_{2z}h})k_{2z}}{e}^{i{k}_{1z}(z-h)},$$(9A)

(10A) $$G_{12}({\mathit{\kappa }}_x,{\mathit{\kappa }}_y,z,z^{\prime })=\frac{i}{8{\mathit{\pi }}^2}{T}_{12}\frac{e^{-i{k}_{2z}(z-h)}-e^{i{k}_{2z}(z+h)}}{(1+R_{21}{e}^{2i{k}_{2z}h})k_{1z}}{e}^{i{k}_{1z}(z^{\prime }-h)}.$$(10A)

At k_1,2 > κ_⊥ these spectra describe the wave (propagating) part of a field; at k_1,2 < κ_⊥ they correspond to near-field (evanescent) components that are sharply fading from a source providing a subwavelength resolution of the proposed diagnostics in the near zone. In this diagnostics, the frequency dependence of the extinction depth of evanescent waves is in use to obtain a subwavelength resolution – unlike the far-field technique [18–25] in hydroacoustic sounding and in the seismic sounding (for example, in the wave migration technique [33]) that uses the time delay as the depth sensitivity parameter.

Inverse problem

Let us consider the case when the scattering inhomogeneities are placed under the bottom in the medium 1, whereas the source-receiver system is located in the medium 2 (water). Then, the scattered field measured in the medium 2 is determined as:

(11A) $$P_1(\mathbf{r})=-\frac{{\mathit{\omega }}^2}{c^2}\underset{V}{\int }{C}_s({\mathbf{r}}^{\prime })G_{12}(\mathbf{r},{\mathbf{r}}^{\prime })P_0({\mathbf{r}}^{\prime })\text{d}{\mathbf{r}}^{\prime }.$$(11A)

It is possible to simplify the problem using the 2D Fourier transform to obtain the 1D integral equation for the depth profile of the transversal spectrum

(12A) $$C_s(k_x,k_y,z)=\frac{1}{4{\mathit{\pi }}^2}\int \int C_s(x,y,z)exp(-i{k}_{x}x-i{k}_{y}y)\text{d}x\text{d}y$$(12A)

However, in general, (11A) isn’t a convolution equation. For that, the probing field should depend on differences of transversal coordinates, i.e. it can be expressed as P₀ = P(x − x′, y − y′, z, z′). This condition can be fulfilled if the relative source-receiver position is unchanged in the scanning, i.e. the offset of the receiver position relative to the position of the source $\mathit{\delta }\mathbf{r}=(\mathit{\delta }x,\mathit{\delta }y,\mathit{\delta }z)=const$. Then, the probing field can be expressed as

(13A) $$P_0({\mathbf{r}}^{\prime })=\underset{V^{″}}{\int }f({\mathbf{r}}^{″}-\mathbf{r}-\mathit{\delta }\mathbf{r})G_{21}({\mathbf{r}}^{\prime },{\mathbf{r}}^{″})\text{d}{\mathbf{r}}^{″}$$(13A)

or, using the k-space representation of the Green function (9A), as

(14A) $$P_0({\mathbf{r}}^{\prime })=\underset{V^{″}}{\int }f({\mathbf{r}}^{″}-\mathbf{r}-\mathit{\delta }\mathbf{r})\int \int G_{21}({\mathit{\kappa }}_x,{\mathit{\kappa }}_y,z^{\prime },z^{″})e^{i[{\mathit{\kappa }}_x(x^{\prime }-x^{″})+{\mathit{\kappa }}_y(y^{\prime }-y^{″})]}\text{d}{\mathit{\kappa }}_x\text{d}{\mathit{\kappa }}_y\text{d}{\mathbf{r}}^{″}.$$(14A)

After the change of variables $x^{″}-x-\mathit{\delta }x\to x^{″},y^{″}-y-\mathit{\delta }y\to y^{″}$, obtain the expression(15A) $$P_0({\mathbf{r}}^{\prime })=\underset{V^{″}}{\int }f(x^{″},y^{″},z^{″}-z-\mathit{\delta }z)\int \int G_{21}({\mathit{\kappa }}_x,{\mathit{\kappa }}_y,z^{\prime },z^{″})e^{-i[{\mathit{\kappa }}_x(x^{″}+\mathit{\delta }x)+{\mathit{\kappa }}_y(y^{″}+\mathit{\delta }y)]}{e}^{i[{\mathit{\kappa }}_x(x^{\prime }-x)+{\mathit{\kappa }}_y(y^{\prime }-y)]}\text{d}{\mathit{\kappa }}_x\text{d}{\mathit{\kappa }}_y\text{d}{\mathbf{r}}^{″},$$(15A)

where the dependence on differences of transversal coordinates is given explicitly.

Finally, making the 2D inverse Fourier transform over transversal coordinates in (11A)

(16A) $$P_1(k_x,k_y,z,\mathit{\delta }\mathbf{r})=\frac{1}{4{\mathit{\pi }}^2}\int \int e^{-i{k}_{x}x-i{k}_{y}y}[-\frac{{\mathit{\omega }}^2}{c^2}\underset{V}{\int }{C}_s({\mathbf{r}}^{\prime })G_{12}(\mathbf{r},{\mathbf{r}}^{\prime })P_0({\mathbf{r}}^{\prime })\text{d}{\mathbf{r}}^{\prime }]\text{d}x\text{d}y,$$(16A)

obtain the desired 1D integral equation:

(17A) $$\begin{array}{cc}{P}_1(k_x,k_y,z,\mathit{\delta }\mathbf{r})=\hfill & -16{\mathit{\pi }}^4\frac{{\mathit{\omega }}^2}{c^2}\underset{z^{\prime }}{\int }{C}_s(k_x,k_y,z^{\prime })\int \int e^{-i({\mathit{\kappa }}_x-k_x)\mathit{\delta }x-i({\mathit{\kappa }}_y-k_y)\mathit{\delta }y}\text{d}{\mathit{\kappa }}_x\text{d}{\mathit{\kappa }}_y\hfill \\ \hfill & \times \underset{z^{″}}{\int }f({\mathit{\kappa }}_x-k_x,{\mathit{\kappa }}_y-k_y,z^{″}-z-\mathit{\delta }z)G_{21}({\mathit{\kappa }}_x-k_x,{\mathit{\kappa }}_y-k_y,z^{\prime },z^{″})G_{12}({\mathit{\kappa }}_x,{\mathit{\kappa }}_y,z,z^{\prime })\text{d}{z}^{″}\text{d}{z}^{\prime }.\hfill \end{array}$$(17A)