Modeling of 2D Acoustic Radiation Patterns as a Control Problem

Abstract

A problem of modeling radiation patterns of wave sources in two-dimensional acoustic tomography is considered. Each source has its own radiation patterns, and their modeling will be used to improve the solvability of inverse problems of recovering the acoustic parameters of human soft tissues and come closer to building a digital twin of acoustic tomography. The problem is considered as a control problem of the right side for the velocities by spatial variables. Two statements are investigated—control by time or space functions. A numerical solution method is implemented. The results of numerical calculations are presented.

1. Introduction

This article considers a problem of modeling radiation patterns of wave sources in two-dimensional acoustic tomography. One of the important applications of acoustic tomography is ultrasound tomography in medicine [1,2,3,4,5,6,7,8,9,10,11,12]. Inverse problems of ultrasound tomography can be considered as coefficient-inverse problems [13,14,15,16,17,18,19,20,21]. A radiation pattern (or antenna pattern or far-field pattern) refers to the directional (angular) dependence of the strength of waves from an antenna [22,23].

A mathematical model of acoustic tomography based on a hyperbolic system of PDE and numerical methods of its solution are presented in [24,25]. The well-posedness of direct problems for hyperbolic systems of equations are presented in [26,27]. Inverse problems for hyperbolic systems are investigated theoretically in [28].

Each source has its own radiation patterns [29], and their modeling will increase the solvability of inverse problems of restoring the acoustic parameters of human soft tissues and come closer to building a digital twin of acoustic tomography.

A possibility of designing the radiation pattern of a structural–acoustic communication system by predicting the pressure distribution was demonstrated in [30].

A finite element method (FEM) was implemented for calculating the far field directivity pattern of individual array elements in phased-array antennas [31].

In [32], the influence of a cross-coupling lattice caused by the filling material on radiation patterns was investigated.

To represent the essence of cross-coupling, the radiation pattern of an individual array element was calculated in pulse mode directly using FEM [33]. Several array design options have been studied to reduce the level of crosstalk.

A method proposed in [34] made it possible to weaken the influence of reflected signals and obtain frequency characteristics that practically coincide with those measured in free-field conditions.

A differential evolution algorithm was applied to synthesize unequally spaced linear antenna arrays [35]. Position-only, position-phase synthesis, and effect of angle resolution have been studied.

A full wave-equation reverse-time, prestack migration algorithm was used to investigate the potential to correct the multi-offset, ground-penetrating radar amplitudes for antenna radiation pattern effects [36].

The concept of directional patterns appeared in marine acoustics when modeling tsunami waves. The initial disturbance may have oscillations in some direction, then the main part of the wave will propagate in this direction [37]. An asymptotically numerical description of the propagation of tsunami waves in an area with an inhomogeneous depth in the vicinity of wave fronts, which may have caustics, has been proposed [38]. A long-wave approximation was used. It was assumed that the size of the initial disturbance area is small compared to the characteristic length of the interval of change in the depth of the bottom and with the distance from the observation point. The description was based on the Maslov canonical operator approach [39]. A potential tsunami model with a simple piston-type source was used, and estimates of the geographical coordinates of two characteristic parameters of the tsunami source and the epicenter were given [40,41].

A parametric approach for classifying the radiation pattern of an acoustic source taking into account the signals captured by several microphones was developed in [42].

A cuckoo search optimisation algorithm was applied to model electromagnetic antennas [43].

To control the radiation pattern, an observation of dipole radiation of sound with a wide bandwidth through a decorated plate with periodic two-dimensional Helmholtz resonators was used [44].

A radiation pattern method of the characteristic basis function was presented as a general modeling method for obtaining high-precision approximations of antenna radiation patterns using a very small number of measurements in [45].

The results of a proposed source centering algorithm [46], which works on plane sections and is aimed at achieving a convex objective function on each plane section, turned out to be more accurate than the previous algorithms for sources with a higher-order radiation pattern, usually at higher frequencies.

A backtracking search optimization algorithm was applied to synthesize the directional pattern [47] of a linear antenna array with specified nulls in the interference directions by controlling only the amplitude, position and phase of the antenna array elements.

A developed new radiation pattern model, in which part of the energy is radiated inside and behind the lobes, the program and output files can be easily implemented to analyze the phenomenon of radio wave propagation in any algorithms and numerical calculations [48].

The results of experimental studies of linear arrays for hydroacoustic applications with a low level of acoustic coupling between elements through the supporting structure of the lattice and a strong coupling of elements with an acoustic load and modeling of the angular and frequency behavior of linear arrays were presented in [49].

A database of acoustic radiation patterns was recorded, modeled and analyzed for musical instruments in [50]. A centering algorithm was used to analyze the complexity of the acoustic radiation pattern in terms of the number of excitation points.

In [51], a comparative analysis of an equivalent source method and a boundary element method was carried out. The equivalent source method is also known as the method of fundamental solutions, the method of superposition of waves, the method of equivalent source, etc.

A signal modeling method was presented that allows measurements to be made at acoustic frequencies below the limits set by time, without echo [52]. The signals received during calibration were modeled, first using a simple model of a stationary response of a transducer and then using an extended model consisting of terms that were used to describe both the stationary and resonant behaviors of the device.

An optimization method based on a flower pollination algorithm was presented [53] for designing a linear antenna array to reduce the maximum level of side lobes by controlling the amplitudes or positions of the array elements.

A simple electrical method for analyzing the physical behavior of piezoelectric transducer arrays has been proposed [54].

A fairly general methodology was proposed for evaluating the radiation pattern of acoustic sources for a wide range of sources without the need for anechoic operating conditions [55].

It was shown [56] that using an equivalent source method makes it possible to reconstruct the sound field from measurements of the sound field excited by this emitter in a tank with reflecting boundaries.

Reconstruction of the radiation pattern of a sound source in free space from measurements of the field excited by this source in a pool was discussed [57]. The reconstruction procedure was based on the use of a reference acoustic monopole.

The modelling of the radiation pattern of acoustic sources was proposed in [58]. The problem was formulated as a control problem of finding arbitrary functions in the right-hand side of the equations for velocities, assuming that the pressure has a given shape.

In this article, we will consider a problem of modeling the radiation pattern of an acoustic wave source, which is formulated as a problem of controlling the right part of a system of equations for the velocities by spacial variables for a hyperbolic system of 2D first-order PDE, which is connected directly with conservation laws of continuum mechanics. Two statements are investigated: in time and in space. Hölder stability and global Carleman estimates for an inverse boundary value problem for dynamic hyperbolic equations were considered [59].

The structure of this paper is as follows. Section 2 contains statements of two inverse problems under consideration. In this section, we reduce the inverse problems to minimizing of a cost functional by a gradient method. Section 3 demonstrates numerical results of the solution of both inverse problems. In the Conclusion, some numerical results are discussed.

2. Problem Formulation

Let us consider the following 2D direct problem of acoustic wave propagation in the form of a first-order system of hyperbolic equations [60]:

$$\begin{array}{c}\frac{\partial u}{\partial t}+\frac{1}{\rho }\frac{\partial p}{\partial x}=q_1(x,y)r_1\left(t\right),\end{array}$$

(1)

$$\begin{array}{c}\frac{\partial v}{\partial t}+\frac{1}{\rho }\frac{\partial p}{\partial y}=q_2(x,y)r_2\left(t\right),\end{array}$$

(2)

$$\begin{array}{c}\frac{\partial p}{\partial t}+\rho c^2\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)={\theta }_s(x,y)I\left(t\right),\end{array}$$

(3)

$$\begin{array}{c}{u,v,p|}_{(x,y)\in \partial \mathsf{\Omega }}=0,\end{array}$$

(4)

$$\begin{array}{c}{u,v,p|}_{t=0}=0.\end{array}$$

(5)

Here, $(x,y)\in \mathsf{\Omega }=\left\{(x,y)\in [0,L]\times [0,L]\right\}$, $t\in (0,T)$ and $u=u(x,y,t)$, $v=v(x,y,t)$ are the components of the velocity vector with respect to x and y; $p=p(x,y,t)$ is the acoustic pressure; $\rho (x,y)$ is the density of the medium; $c(x,y)$ is the speed of the wave propagation; ${\theta }_s$ is the source location; and $I\left(t\right)$ is a Ricker wavelet (${\nu }_0$ is a frequency):

$$I(t,{\nu }_0)=\left(1-2{\left(\pi {\nu }_0\left(t-\frac{1}{{\nu }_0}\right)\right)}^2\right){\mathrm{e}}^{-\pi {\nu }_0\left(t-\frac{1}{{\nu }_0}\right)}.$$

(6)

The Equations (1)–(3) are the conservation laws of momentum and mass [60].

Let us denote as inverse problem I the problem of the reconstruction of $r_1\left(t\right)$ and $r_2\left(t\right)$ from the inverse problem data:

$$\begin{array}{c}p(x,y,t)=f(x,y,t),(x,y)\in \mathsf{\Omega },t\in (0,T).\end{array}$$

(7)

Stability and global uniqueness for source and coefficients problems for wave equations were investigated in [61,62,63,64].

We should also mention [65,66], where the authors considered the problem of detecting the source of acoustical noise in the Helmholtz equation from measurements of the acoustical pressure field inside the domain, proved the uniqueness of the reconstruction and obtained a sharp stability estimate.

In addition, in [67], the authors proposed a general framework to study uniqueness, stability and reconstruction for infinite-dimensional inverse problems when only a finite-dimensional approximation of the measurements is available. It is shown that the same estimate holds for a large class of inverse problems satisfying a Lipschitz stability, even with a finite number of measurements. The authors also derived a globally convergent reconstruction algorithm based on the Landweber iteration.

We reduce the inverse problem (1)–(5), (7) to the problem of minimization of the cost functional:

$$J(r_1,r_2)={\int }_0^T{\int }_{\mathsf{\Omega }}{\left[p(x,y,t;r_1,r_2)-f(x,y,t)\right]}^2\mathrm{d}x\mathrm{d}y\mathrm{d}t\to \underset{r_1,r_2}{min}$$

(8)

by a gradient method.

The gradient of the cost functional is calculated as follows [68]:

$$\begin{array}{c}{J}^{\prime }\left(r_1\right)\left(t\right)={\int }_{\mathsf{\Omega }}\rho (x,y)q_1(x,y){\mathsf{\Psi }}_1(x,y,t)\mathrm{d}\mathsf{\Omega },\end{array}$$

(9)

$$\begin{array}{c}{J}^{\prime }\left(r_2\right)\left(t\right)={\int }_{\mathsf{\Omega }}\rho (x,y)q_2(x,y){\mathsf{\Psi }}_2(x,y,t)\mathrm{d}\mathsf{\Omega }.\end{array}$$

(10)

Here, the functions ${\mathsf{\Psi }}_i$, $i=1,2,3$ solve the following adjoint problem:

$$\begin{array}{c}{{\mathsf{\Psi }}_1}_t+\frac{1}{\rho }{{\mathsf{\Psi }}_3}_x=0;\end{array}$$

(11)

$$\begin{array}{c}{{\mathsf{\Psi }}_2}_t+\frac{1}{\rho }{{\mathsf{\Psi }}_3}_y=0;\end{array}$$

(12)

$$\begin{array}{c}\frac{1}{\rho c^2}{{\mathsf{\Psi }}_3}_t+({{\mathsf{\Psi }}_1}_x+{{\mathsf{\Psi }}_2}_y)=2\left[p(x,y,t)-f(x,y,t)\right];\end{array}$$

(13)

$$\begin{array}{c}{\mathsf{\Psi }}_i(x,y,T)=0;\end{array}$$

(14)

$$\begin{array}{c}{\mathsf{\Psi }}_i{|}_{(x,y)\in \partial \mathsf{\Omega }}=0.\end{array}$$

(15)

Let inverse problem II be the problem of reconstruction of $q_1(x,y)$ and $q_2(x,y)$ in (1)–(5) from the known inverse problem data (7).

Note that if the inverse problem data (7) are used in Equations (1)–(3) and all functions are smooth enough, we obtain the following equation:

$$\frac{\partial q_1(x,y)}{\partial x}{r}_1\left(t\right)+\frac{\partial q_2(x,y)}{\partial y}{r}_2\left(t\right)={\theta }_s(x,y)I^{\prime }(t,{\nu }_0)-\frac{1}{\rho c^2}\frac{{\partial }^{2}f}{\partial t^2}+\frac{\partial }{\partial x}\left(\frac{1}{\rho }\frac{\partial f}{\partial x}\right)+\frac{\partial }{\partial y}\left(\frac{1}{\rho }\frac{\partial f}{\partial y}\right).$$

In inverse problem I, we find the functions $r_1\left(t\right)$ and $r_2\left(t\right)$; in inverse problem II, the functions $q_1(x,y)$ and $q_2(x,y)$ are unknowns. So far, no theoretical results on the uniqueness and stability of the solution have been obtained.

We reduce the inverse problem (1)–(5), (7) to the problem of minimization of the following cost functional:

$$J(q_1,q_2)={\int }_0^T{\int }_{\mathsf{\Omega }}{\left[p(x,y,t;q_1,q_2)-f(x,y,t)\right]}^2\mathrm{d}x\mathrm{d}y\mathrm{d}t\to \underset{q_1,q_2}{min}.$$

(16)

by a gradient method.

The gradient is calculated as follows:

$$\begin{array}{c}{J}^{\prime }\left(q_1\right)(x,y)={\int }_0^T{r}_1\left(t\right)\rho (x,y){\psi }_1(x,y,t)dt,\end{array}$$

(17)

$$\begin{array}{c}{J}^{\prime }\left(q_2\right)(x,y)={\int }_0^T{r}_2\left(t\right)\rho (x,y){\psi }_2(x,y,t)dt.\end{array}$$

(18)

Here, the functions ${\mathsf{\Psi }}_i$, $i=1,2$ are the solution of the adjoint problem (11)–(15).

At each iteration of the gradient method, we solve both the direct and adjoint problems. A comparative analysis of two gradient calculation strategies has been carried out [68]. It was shown that calculating the gradient from all sources and receivers gives the best approximate solution at the next iteration, in contrast to the case when the gradient is calculated from a single source. A new approach to calculating the gradient of the functional has been developed and implemented [69]. It saves almost twice the RAM and speeds up the inverse problem solving code by up to 25%, depending on the number of CPU cores used. The convergence of the gradient methods for hyperbolic equations is investigated in [70]. It has been shown [71] that adding a priori information about the inverse problem solution to the gradient method can significantly decrease the number of iterations. In [72], a stopping criterion for the gradient method, consistent with the accumulation of machine round-off errors, is it presented.

3. Numerical Results

We apply the MUSCL–Hancock method for the direct and adjoint problem solution. This method was introduced in [73] as a variant of the MUSCL scheme. This scheme is full time and space second-order accurate. It differs from the other MUSCL variants [74,75,76,77] in its second-order time accuracy. Indeed, the time accuracy rises when using a step based on a central-differencing scheme [78,79]. The $L_1$-stability of the MUSCL–Hancock scheme is considered in [80]. The influence of the method of solving the direct and adjoint problem on the solution of the inverse problem was investigated in [81]. The accuracy of the numerical simulation of nonlinear acoustic wave propagation in a liquid medium based on the Navier–Stokes equations is investigated in [82].

We use a homogeneous medium with acoustic parameters of the medium corresponding to water: $\rho (x,y)=1000$ kg/m${}^3$, $c(x,y)=1400$ m/s. The computational domain is $[0,0.3]\times [0,0.3]$ m. The grid parameters are chosen as follows: $N_x=N_y=300$, the CFL number is equal to 0.5. The initial guess is zero.

We conducted two tests for the two inverse problem statements. In Test 1, we know the exact solution for some model data. In Test 2, no exact solution is known. The data of the inverse problem are obtained from solving the direct problem after zeroing the pressure outside the desired sector.

For inverse problem 1 (Test 1 and Test 2), the functions are $q_1(x,y)=q_2(x,y)=1$. For Test 1, $r_1\left(t\right)$, $r_2\left(t\right)$ are as follows:

$$r_1\left(t\right)=I(t,20),r_2\left(t\right)=I(t,30).$$

Figure 1 presents the exact solution (green line) and inverse problem solution after 2000 iterations (using exact data). Figure 2 presents the inverse problem solution after 2000 iteration using noisy data $f_{\delta }(x,y,t)$ with $\delta =5\%$ and $\delta =10\%$.

Let us consider the results of numerical calculations for a sector with additional information on the pressure obtained as a solution to the direct problem (1)–(5) and then set to zero outside the specified sector.

Let us present the results of numerical calculations for inverse problem II. In order to make dimensions of the equations coincide, we consider the following:

$$r_1\left(t\right)=r_2\left(t\right)=\underset{0}{\overset{t}{\int }}I\left(\tau \right)\mathrm{d}\tau .$$

We also studied the dependence of the algorithm in the cases of exact and noisy data. We added some noise to the inverse problem data as follows:

$$f_{\delta }(x,y,t)=f(x,y,t)\left(1+\alpha (x,y,t)\frac{\delta }{100}\right).$$

Here, $f(x,y,t)$ are the exact data, $\delta$ is the percentage of noise level in the data and $\alpha (x,y,t)$ is a random variable uniformly distributed on the interval $(-1,1)$ for each variable x, y and t.

Figure 3 presents the exact solution and approximate solution of the inverse problem functions $q_1(x,y)$ and $q_2(x,y)$ after 1500 iterations. Figure 4 presents the exact solution and approximate solution of the inverse problem with noise in data $\delta =5\%$ and $\delta =10\%$, respectively. The inverse problem is quite stable, even for noisy data, due to the following consideration: We recover two functions of two spatial variables, $q_1(x,y)$ and $q_2(x,y)$, from the known inverse problem data of a function of three variables (two spatial variables and one time variable) $f(x,y,t)$.

Next, we consider the results of two numerical calculations for wide and narrow sectors. In this case, for some arbitrary data, no inverse problem solution may exist.

Figure 5 shows additional information (a sector of 90 degrees). The top of the sector is located on the right boundary.

The left column of Figure 6 presents the pressure with functions $q_1(x,y)=q_2(x,y)=0$. The right column of the Figure 6 presents the pressure with restored functions $q_1(x,y)$, $q_2(x,y)$ (see Figure 7) after 900 iterations.

Figure 8 considers the inverse problem data for a narrow sector. The sector moves from the middle of the side to the corners of the opposite one. The top of the sector is right on the border.

In the right column of Figure 9, the pressure is presented with the restored functions $q_1(x,y)$ and $q_2(x,y)$ (see Figure 10). Figure 10 presents functions $q_1,q_2$, obtained from the additional information (narrow sector) when solving the inverse problem. On Figure 11 $logJ(r_1,r_2)$ on nth iteration is presented (green line is for Test 3, red line is for Test 2).

4. Concluding Remarks

Two ways to simulate radiation patterns have been presented. The modelling of radiation patterns is formulated as control problems, which are solved by a gradient method. A numerical method for solving the inverse problems was implemented, and the results of numerical experiments were presented.

The numerical experiments have shown that if we simulate the radiation patterns by controlling $r_1\left(t\right)$ and $r_2\left(t\right)$ we may not obtain suitable results because the solution of the inverse problem of finding $r_1$ and $r_2$ may not exist. For the case of control by $q_1(x,y)$, $q_2(x,y)$ it is possible to reduce the pressure outside the sector by 15–20%, which is considered to be a good result.

5. Discussion

This work aimed to investigate some aspects of mathematical models of ultrasound tomography related to the radiation pattern of a sounding pulse and study a way to bring the model closer to the physical process. Naturally, the physical structure of the source of a sounding wave provides the direction in which the main part of sounding impulse propagates. The goal is to reflect such behaviour when considering an acoustic approximation based on a first order system of equations in time domain. This leads us to the control problem considered in this paper. One formulation was aimed to find a suitable wave time form that would provide the desired behaviour. The other, when we use functions of the spatial variables as means of control, changes the properties of the medium in the area, which corresponds to the source location. A natural way to continue this study is to try to use several sources and to try to control the radiation pattern by using the locations of these sources. In addition, the spatial functions of the second inverse problem can be considered as changes of the absorption properties in the area of the source. This connection will be studied as the next step of our research by introducing absorption parameters into the model.