A fast algorithm for solving a three-dimensional inverse multiple frequency problems of scalar acoustics in a cylindrical region
aa r X i v : . [ m a t h . NA ] F e b A fast algorithm for solving a three-dimensional inverse multiple frequencyproblems of scalar acoustics in a cylindrical regionAnatoly B. Bakushinsky , Alexander S. Leonov Abstract
A new algorithm for the stable solution of a three-dimensional scalar inverseproblem of acoustic sounding of an inhomogeneous medium in a cylindrical regionis proposed. The data of the problem is the complex amplitude of the wave field,measured outside the region of acoustic inhomogeneities in a cylindrical layer. Us-ing the Fourier transform and Fourier series, the inverse problem is reduced tosolving a set of one-dimensional Fredholm integral equations of the first kind, tothe subsequent calculation of the complex amplitude of the wave field in the regionof inhomogeneity, and then to finding the required sound velocity field in this re-gion. The algorithm allows solving the inverse problem on a personal computer ofaverage performance for sufficiently fine three-dimensional grids in tens of seconds.A numerical study of the accuracy of the proposed algorithm for solving model in-verse problems at various frequencies is carried out, and the issues of stability ofthe algorithm with respect to data perturbations are investigated.
Keywords:
Three-dimensional wave equation, inverse coefficient problem, regulariz-ing algorithm, fast Fourier transform.
AMS Mathematics Subject Classification:
Let the scalar function p ( x , t ) define an acoustic wave field depending on the coordinates x = ( x, y, z ) and time t ≥ Q ⊂ R . The domain is an infinite cylinder ofthe form Q = { ( x, y, z ) : x + y ≤ b , z ∈ R } . The field is created by sources localizedin a known region S . The medium is characterized by the local phase velocity of sound c ( x ) and has a constant density. Moreover, it is known that c ( x ) = c = const outsidethe given region X, X ⊂ Q , X ∩ S = ∅ . In the region X , the function c ( x ) can vary,and this is interpreted as the location of acoustic inhomogeneities there. In this case, thefield p ( x , t ) for a harmonic source of the form f ( x , ω ) e iωt with a known frequency ω canbe found within the linear acoustics approximation as p ( x , t ) = u ( x , ω ) e iωt , where thecomplex amplitude u ( x , ω ) satisfies the equation∆ u ( x , ω ) + k u ( x , ω ) = f ( x , ω ) + ω ξ ( x ) u ( x , ω ) , x ∈ Q. (1)Here k = ωc , and ξ ( x ) = c − − c − ( x ). We also assume that the boundary condition (cid:18) ∂u∂ n + σ ( x ) u (cid:19) ∂Q = 0 Federal Research Center Computer Science and Control of Russian Academy of Sciences, Institutefor Systems Analysis, 117312 Moscow; Mari State University, Lenin Sqr. 1, 424000 Yoshkar-Ola Department of Mathematics, National Nuclear Research University ’MEPHI’, Kashirskoe Shosse 31,115409 Moscow
1s satisfied with known function σ ( x ), and the radiation condition is fulfilled in the coor-dinate z . Here n is the outer normal to the surface ∂Q . Without going into details of theconditions for the coefficients σ ( x ) and f ( x , ω ), we make the following assumption. Assumption 1.
The function ξ ( x ) is continuous with compact support in X and thecorresponding problem (1) with the indicated additional conditions has a unique solution u ( x , ω ) ∈ H ( Q ) for each considered ω .Finding such a function u ( x , ω ) is direct problem . We are interested in the following inverse problem for the equation with the indicated additional conditions: knowing thecomplex amplitude of the field u ( x , ω ) for some set of frequencies ω in the domain Y ( Y ⊂ Q, Y ∩ X = ∅ , Y ∩ S = ∅ ) find the coefficient ξ ( x ), i.e. the function c ( x ) definingacoustic inhomogeneities in the region X . Introducing the Green’s function G ( x , x ′ , ω ) forthe Helmholtz equation in the domain Q , we can reduce the inverse problem under certainassumptions about the smoothness of the functions u ( x , ω ) , f ( x , ω ) , c ( x ) to a nonlinearsystem of integral equations for the unknowns u ( x ′ , ω ) and ξ ( x ′ ) , x ′ ∈ X : u ( x , ω ) = u ( x , ω ) + ω Z X G ( x , x ′ , ω ) ξ ( x ′ ) u ( x ′ , ω ) d x ′ , x ∈ X, (2) ω Z X G ( x , x ′ , ω ) ξ ( x ′ ) u ( x ′ , ω ) d x ′ = w ( x , ω ) , x ∈ Y (see, for example, [1, 2, 3], etc.). The functions u ( x , ω ) = Z X G ( x , x ′ , ω ) f ( x ′ , ω ) d x ′ , w ( x , ω ) = u ( x , ω ) − u ( x , ω ) , x ∈ Y, included in (2) are known (computable) functions, and the values u ( x ′ , ω ) and ξ ( x ′ ) , x ′ ∈ X must be defined.The problems similar to (2) has been sufficiently well investigated theoretically. Inparticular, questions of the existence and uniqueness of their solution for various domains Q, X, Y have been studied (see, for example, [2]–[4] and etc.). We emphasize that thepurpose of this article is not to study the properties of the inverse problem (2). We onlypropose an effective numerical method for solving it in the case of cylindrical regions
X, Y . That is why we do not carry out a detailed reduction of the inverse problem to thesystem (2) indicating all the requirements for the coefficients, but make only one moreassumption.
Assumption 2.
Reduction of the inverse problem to the (2) system is possible forgiven u ( x , ω ) ∈ L ( Q ). The system (2) is solvable and defines, for all considered ω , somesolution of the inverse problem, that is a continuous function ξ ( x ) with compact supportin X and the function u ( x , ω ) ∈ L ( X ).Note that under such assumption the system (2) may have more than one solution.In many papers, various numerical methods for solving problems similar to (2) intwo-dimensional and three-dimensional formulations are considered. For example, in thework [3] the system (2) is reduced to a nonlinear operator equation that is then solvedby a special iterative method. This approach does not take into account explicitly thatthis equation is an ill-posed problem, and the iterative method used is not regularized.2evertheless, the method works for model scatterers of medium strength [3, p.101-105].In the work [4], the regularized Gauss – Newton method was used to solve the system ofequations (2) in the three-dimensional axially symmetric case. A special gradient methodand the Fletcher-Reeves method were used for a similar problem in the paper [5]. Othergradient techniques have been used in [6],[7].Note that there are also alternative approaches to solving the inverse problem ofacoustic sounding that are not related to a system of the type (2). In particular, theoriginal method of boundary control was proposed and developed in the works [8],[9] thatallows solving three-dimensional problems (2). In [10], the well-known method of R.G.Novikov [11] was applied to solve the two-dimensional inverse acoustic scattering problem.In the subsequent work [12], a comparative analysis of a variant of this method andsome other functional-analytical methods for solving two-dimensional inverse problemsof acoustic scattering was carried out. Methods by M.V. Klibanov summarized in themonograph [13], as well as methods from the monograph [14], turned out to be verypromising in processing real experimental data. Also of great interest are recent works[15],[16].All the above methods for solving the inverse problem of acoustic sounding requiresignificant computational resources in a three-dimensional formulation. In this regard,we note the articles [17],[18], in which the original inverse coefficient problem for thewave equation, from which the equation (1) is actually obtained, is reduced to a three-dimensional Fredholm integral equation of the first kind with the right-hand side con-taining special integrals of the recorded field. The proposed method for solving thisequation turned out to be very effective numerically, and allows one to quickly solvethree-dimensional inverse problems for sufficiently fine grids on a personal computer (PC)even without parallelization. This method is also very useful for the (2) problem, and wewill demonstrate it below.In this article, we adhere to the following scheme for solving a nonlinear system(2) foreach considered ω .1) Introducing the notation v ( x ′ , ω ) = ξ ( x ′ ) u ( x ′ , ω ) , x ′ ∈ X , solve the second equation(linear integral Fredholm equation of the first kind), written in the form ω Z X G ( x , x ′ , ω ) v ( x ′ , ω ) d x ′ = w ( x , ω ) , x ∈ Y, (3)for the functions v ( x ′ , ω ).2) With found function v ( x ′ , ω ), calculate the function u ( x , ω ) , x ∈ X, from the firstequality of the system (2), written in the form u ( x , ω ) = u ( x , ω ) + ω Z X G ( x , x ′ , ω ) v ( x ′ , ω ) d x ′ , x ∈ X. (4)3) Find the function ξ ( x ) from the equation v ( x , ω ) = ξ ( x ) u ( x , ω ) , x ∈ X usingcomputed functions v ( x , ω ) , u ( x , ω ).A similar scheme has been used before. For example, in [19], it was used to solve anill-posed problem, a three-dimensional equation similar to (3). However, regularizationmethods were not used there, and the solution of the inverse problem was sought for arather narrow class of functions with a piecewise constant current.3elow we will show that under some not very burdensome special assumptions aboutthe set of data registration Y and about the set X , where the solution of the inverseproblem is sought, scheme 1) - 3) is effectively implemented numerically and allows solvingthe corresponding three-dimensional inverse problem for fairly fine grids even on a personalcomputer of average performance in a few tens of seconds without parallelization. Theproposed algorithm for solving the inverse problem and its numerical study are the mainresults of this work.The article is organized as follows. In Sect.2, we consider the geometric scheme usedand reduction of the inverse problem to a final system of equations. In Sect.3, a methodfor obtaining model data for the inverse problem (Algorithm 1) is given, and Algorithm 2for solving our inverse problem in a cylindrical domain is presented and discussed. Sect.4is devoted to a finite-dimensional approximation and solution of model inverse problems.Here we consider some details of algorithms and present results of numerical experiments.The properties of the algorithm, such as the speed of finding a solution and sensitivityof solutions to input data errors, are discussed in Sect.5. Finally, in Sect.6 we formulatemain conclusions. Everywhere below, the domains
X, Y have the form of cylinders: the solution domain is X = { ≤ r ≤ a } × R z , r = p x + y , and the observation region is Y = { r ≤ r ≤ b } × R z . Figure 1 shows the disposition of these regions. The possible position of thesources is also conventionally shown there. We consider the cylinders to be infinite in thevariable z : z ∈ ( −∞ , + ∞ ), bearing in mind that the sought-for function ξ is compactlysupported and in particular with respect to this variable. In what follows we will denoteas X xy , Y xy sections of cylindrical areas X, Y by plane, perpendicular to the Oz axis.Green’s function for the problem (1) with the indicated additional conditions in thecylindrical domain Q is well known, and here we do not write it out explicitly. Note animportant feature of this function: in cylindrical coordinates x = ( r, ϕ, z ) , x ′ = ( r ′ , ϕ ′ , z ′ )it has the form: G ( x , x ′ , ω ) = G (cid:18)q r + r ′ − rr ′ cos( ϕ − ϕ ′ ) + ( z − z ′ ) ; ω (cid:19) == G ( r, r ′ , cos( ϕ − ϕ ′ ) , z − z ′ ; ω ) , and the expression for the function G can be found, for example, in [20]. In the samecoordinates, we represent the functions u ( x , ω ) = u ( r, ϕ, z ) , v ( x , ω ) = v ( r, ϕ, z ). Thenthe equations (3), (4) can be written as follows: ω a Z π Z ∞ Z −∞ G ( r, r ′ , cos( ϕ − ϕ ′ ) , z − z ′ ; ω ) v ( r ′ , ϕ ′ , z ′ ; ω ) r ′ dr ′ dϕ ′ dz ′ == w ( r, ϕ, z ; ω ) , r ∈ [ r , b ] , ϕ ∈ [0 , π ] , z ∈ R , (5)4 ba r * *** * * * *Y Y X Figure 1:
Geometric scheme of data registration for the inverse problem: X is the domain ofwave field scatterers, Y is the region of data registration, asterisks are the conditional positionsof the sources. and u ( r, ϕ, z ; ω ) − u ( r, ϕ, z ; ω ) == ω a Z π Z ∞ Z −∞ G ( r, r ′ , cos( ϕ − ϕ ′ ) , z − z ′ ; ω ) v ( r ′ , ϕ ′ , z ′ ; ω ) r ′ dr ′ dϕ ′ dz ′ ,r ∈ [0 , a ] , ϕ ∈ [0 , π ] , z ∈ R . (6)Here v ( r ′ , ϕ ′ , z ′ ; ω ) = ξ ( r ′ , ϕ ′ , z ′ ) u ( r ′ , ϕ ′ , z ′ ; ω ). Next, we introduce the Fourier transforms F z [ · ](Ω) = R + ∞−∞ [ · ] e i Ω z dz with respect to the variable z (or z ′ ) for functions G , v, u, u , w as elements of space L (see Assumption 2):˜ G ( r, r ′ , cos ϕ, Ω; ω ) = F z [ G ( r, r ′ , cos ϕ, z ; ω )] (Ω) , ˜ w ( r, ϕ, Ω; ω ) = F z [ w ( r, ϕ, z ; ω )] (Ω) , ˜ u ( r, ϕ, Ω; ω ) = F z [ u ( r, ϕ, z ; ω )] (Ω) , ˜ u ( r, ϕ, Ω; ω ) = F z [ u ( r, ϕ, z ; ω )] (Ω) . Then, by the convolution theorem, the equality (5) and (6) can be written as ω a Z π Z ˜ G ( r, r ′ , cos( ϕ − ϕ ′ ) , Ω; ω )˜ v ( r ′ , ϕ ′ , Ω; ω ) r ′ dr ′ dϕ ′ == ˜ w ( r, ϕ, Ω; ω ) , r ∈ [ r , b ] , (7)5 u ( r, ϕ, Ω; ω ) − ˜ u ( r, ϕ, Ω; ω ) == ω a Z π Z ˜ G ( r, r ′ , cos( ϕ − ϕ ′ ) , Ω; ω )˜ v ( r ′ , ϕ ′ , Ω; ω ) r ′ dr ′ dϕ ′ , r ∈ [0 , a ] , (8)with ˜ v ( r ′ , ϕ ′ , Ω; ω ) = F z [ v ( r ′ , ϕ ′ , z ′ ; ω )] (Ω) = F z [ ξ ( r ′ , ϕ ′ , z ′ ) u ( r ′ , ϕ ′ , z ′ ; ω )]. We also intro-duce expansions of the functions ˜ G, ˜ v, ˜ u, ˜ u , ˜ w in the basis system of functions { e inϕ } , n ∈ Z , in the space L (0 , π ):˜ G ( r, r ′ , cos ϕ, Ω; ω ) = P n G n ( r, r ′ , Ω; ω ) e inϕ , ˜ v ( r ′ , ϕ ′ , Ω; ω ) = P n v n ( r ′ , Ω; ω ) e inϕ ′ , ˜ u ( r, ϕ, Ω; ω ) = P n u n ( r ′ , Ω; ω ) e inϕ , ˜ u ( r, ϕ, Ω; ω ) = P n u n ( r ′ , Ω; ω ) e inϕ ;˜ w ( r, ϕ, Ω; ω ) = P n w n ( r ′ , Ω; ω ) e inϕ ; ϕ, ϕ ′ ∈ [0 , π ]with coefficients G n ( r, r ′ , Ω; ω ) = π R π ˜ G ( r, r ′ , cos ϕ, Ω; ω ) e − inϕ dϕ,v n ( r ′ , Ω; ω ) = π R π ˜ v ( r ′ , ϕ ′ , Ω; ω ) e − inϕ ′ dϕ ′ ,w n ( r, Ω; ω ) = π R π ˜ w ( r, ϕ, Ω; ω ) e − inϕ dϕ,u n ( r, Ω; ω ) = π R π ˜ u ( r, ϕ, Ω; ω ) e − inϕ dϕ,u n ( r, Ω; ω ) = π R π ˜ u ( r, ϕ, Ω; ω ) e − inϕ dϕ. Then the relations (7), (8) can be reduced to the following system of equalities that istrue for any n ∈ Z , Ω ∈ R and all considered ω : ω Z a G n ( r, r ′ , Ω; ω ) v n ( r ′ , Ω; ω ) r ′ dr ′ = 12 π w n ( r, Ω; ω ) , r ∈ [ r , b ] , (9)12 π ( u n ( r, Ω; ω ) − u n ( r, Ω; ω )) = ω Z a G n ( r, r ′ , Ω; ω ) v n ( r ′ , Ω; ω ) r ′ dr ′ , r ∈ [0 , a ] . (10)Moreover, the equality v n ( r ′ , Ω; ω ) = 12 π Z π F z [ v ( r ′ , ϕ, z ; ω )] e − inϕ dϕ == 12 π Z π F z [ ξ ( r ′ , ϕ, z ) u ( r ′ , ϕ, z ; ω )] e − inϕ dϕ == 12 π Z π F z " ξ ( r ′ , ϕ, z ) X m u m ( r ′ , Ω; ω ) e imϕ e − inϕ dϕ (11)is satisfied. The relations (9), (10) are equations for the unknown functions v n ( r ′ , Ω; ω )and u n ( r, Ω; ω ), depending on one variable r ′ or r . Other arguments of these functions,i.e. Ω , ω , are parameters. 6 Algorithms for solving direct and inverse problemsin a cylindrical domain
The direct problem considered below is to find the function w ( r, ϕ, z, ω ) , r ∈ [ r , b ] , ϕ ∈ [0 , π ] , z ∈ R , from the equalities (5), (6) using the known finite function ξ ( r ′ , ϕ ′ , z ′ ) andthe given function of sources u ( r, ϕ, z ; ω ) for a set of frequencies ω under consideration.For this, the equalities (5), (6) are reduced to the system (7),(8), and then to the systemof relations (9) – (11). As a result, the calculation of the function w can be representedas the following algorithm. Algorithm 1
Step 1). For the set of frequencies ω under consideration, calculate the Fourier trans-forms in z :˜ G ( r, r ′ , cos ϕ, Ω; ω ) = F z [ G ( r, r ′ , cos ϕ, z ; ω )] (Ω) , ˜ u ( r, ϕ, Ω; ω ) = F z [ u ( r, ϕ, z ; ω )] (Ω)and expand the resulting functions into Fourier series in the variable ϕ :˜ G ( r, r ′ , cos ϕ, Ω; ω ) = X n G n ( r, r ′ , Ω; ω ) e inϕ , ˜ u ( r, ϕ, Ω; ω ) = X n u n ( r, r ′ , Ω; ω ) e inϕ . Both of these procedures can be implemented using the Fast Discrete Fourier Transform(FFT).Step 2). For each parameters ω, Ω, we implement the following iterative process ofsolving the equations (10), (11) with respect to the set of function { u n ( r, Ω; ω ) } , n ∈ Z : v ( k ) n ( r ′ , Ω; ω ) = 12 π Z π F z " ξ ( r ′ , ϕ, z ) X m u ( k ) m ( r ′ , Ω; ω ) e imϕ e − inϕ dϕ, r ′ ∈ [0 , a ] , (12) u ( k +1) n ( r, Ω; ω ) = u n ( r, Ω; ω ) + 2 πω Z a G n ( r, r ′ , Ω; ω ) v ( k ) n ( r ′ , Ω; ω ) r ′ dr ′ ,r ∈ [0 , a ] , k = 0 , , , ..., (13)with an initial guess { u (0) n ( r, Ω; ω ) } = { u n ( r, Ω; ω ) } .Step 3). Stop the process by some rule at iteration number ν and obtain an approxi-mate solution { u ( ν ) n ( r, Ω; ω ) } of the system of equations (10), (11), and related functions { v ( ν ) n ( r, Ω; ω ) } . 7tep 4). Calculate an approximate set of values { w n ( r, Ω; ω ) } from (9) : w ( ν ) n ( r, Ω; ω ) = 2 πω Z a G n ( r, r ′ , Ω; ω ) v ( ν ) n ( r ′ , Ω; ω ) r ′ dr ′ , r ∈ [ r , b ] (14)and accept the function ˜ w ( ν ) ( r, ϕ, Ω; ω ) = X n w ( ν ) n ( r ′ , Ω; ω ) e inϕ (15)or its inverse Fourier transform as an approximate solution of the direct problem.We will not carry out a theoretical analysis of the convergence of Algorithm 1 here,since we do not formally use it in solving the inverse problem. It is only needed togenerate the model data w . We only note that, as follows from the general theory ofsolving integral equations of the second kind (see, for example, [21]), the algorithm willrapidly converge, at least for small ω . Below we will demonstrate numerical examplesconfirming this statement. The task is to find for each frequency ω the solution ξ ( r ′ , ϕ ′ , z ′ ) of the system (5), (6)using given function w ( r, ϕ, z, ω ) , r ∈ [ r , b ] , ϕ ∈ [0 , π ] , z ∈ R , and given source function u ( r, ϕ, z ; ω ). To do this, we reduce the problem (5), (6) to the system (9) – (11), assumingthe sets of data functions { w n ( r, Ω; ω ) } and { u n ( r, Ω; ω ) } are calculated. Then from thesystem and these sets we find the function ξ . The solution procedure is presented in theform of the following algorithm. Algorithm 2
Step 1). For all considered parameters ω, Ω and all used n we solve one-dimensionalintegral equations of the first kind corresponding to the equalities (9): ω Z a G n ( r, r ′ , Ω; ω ) v n ( r ′ , Ω; ω ) r ′ dr ′ = 12 π w n ( r, Ω; ω ) , r ∈ [ r , b ] . (16)Here we use a suitable regularization method (regularizing algorithm, RA) for these ill-posed problems. The result is a set of approximate solutions { v n ( r ′ , Ω; ω ) } .Step 2). Using the found functions { v n ( r ′ , Ω; ω ) } , we calculate from the equalities (10)the set of functions { u n ( r, Ω; ω ) } : u n ( r, Ω; ω ) = u n ( r, Ω; ω ) + 2 πω Z a G n ( r, r ′ , Ω; ω ) v n ( r ′ , Ω; ω ) r ′ dr ′ , r ∈ [0 , a ] . (17)Step 3). We restore the functions v ( r, ϕ, z, ω ) and u ( r, ϕ, z, ω ) for r ∈ [0 , a ] , ϕ ∈ [0 , π ] , z ∈ R using the sets { v n ( r, Ω; ω ) } and { u n ( r, Ω; ω ) } , summing the correspondingFourier series in the variable ϕ and then calculating inverse Fourier transform in thevariable z , F − [ · ]: v ( r, ϕ, z ; ω ) = F − "X n v n ( r, Ω; ω ) e inϕ ( z ) ,u ( r, ϕ, z ; ω ) = F − "X n u n ( r, Ω; ω ) e inϕ ( z ) , ( r, ϕ, z ) ∈ X, ξ ( r, ϕ, z ) from the equation u ( r, ϕ, z, ω ) ξ ( r, ϕ, z ) = v ( r, ϕ, z, ω )for each point ( r, ϕ, z ) ∈ X . This can be done for each ω under consideration, and theresult will generally depend on ω . Further, for example, one can average the results byone or another method over the value of ω .Step 5). Finally, we calculate the function c ( r, ϕ, z ) from the equality ξ ( r, ϕ, z ) = c − − c − ( r, ϕ, z ).Let’s make some comments on Algorithm 2.a) To implement Step 1 it is necessary to clarify the properties of the functions v, w .Using Assumption 2, we presume that the inclusions v ( x , ω ) ∈ L ( X ), w ( x , ω ) ∈ L ( Y )are valid for each considered frequency ω . The first inclusion follows from the compact-ness of the support of the function ξ , and the second is postulated. In this case, theknown methods for solving linear ill-posed problems in Hilbert spaces are applicable tothe equations (16) (see, for example, [4],[5],[22],[23],[24], etc.).b) In solving our inverse problem, the equation (3), i.e. (5) and the equations (16)generated by it, may have more than one solution for the finite set of frequencies ω used.Therefore, it is important to establish a connection between the solutions of all theseequations. The connection is substantiated by the following statements. Theorem 1
1) Let w ( r, ϕ, z ; ω ) ∈ L ( Y ) for every ω . Then any solution v ( r, ϕ, z ; ω ) ∈ L ( X ) to the equation (5) can be represented in the form v ( r, ϕ, z ; ω ) = F − "X n v n ( r, Ω; ω ) e inϕ ( z ) , ( r, ϕ, z ) ∈ X, (18) where the functions v n ( r, Ω; ω ) ∈ L { [0 , a ] × R Ω } satisfy the integral equations (9) foreach ω . Conversely, if v n ( r, Ω; ω ) ∈ L { [0 , a ] × R Ω } are solutions of equations (9) suchthat P n k v n ( r, Ω; ω ) k L { [0 ,a ] × R Ω } < ∞ for each ω , then a function of the form (18) is thesolution to the equation (5) .2) For each ω the equality k v ( r, ϕ, z ; ω ) k L ( X ) = X n k v n ( r, Ω; ω ) k L { [0 ,a ] × R Ω } (19) holds. Corollary 1
Let the functions ¯ v n ( r, Ω; ω ) ∈ L { [0 , a ] × R Ω } be normal solutions to equa-tions (9) (solutions with minimal norm). Then the function ¯ v ( r, ϕ, z ; ω ) = F − "X n ¯ v n ( r, Ω; ω ) e inϕ ( z ) , ( r, ϕ, z ) ∈ X, is the unique normal solution to the equation (5) . Everywhere below, it is assumed that the equations (3), (4) and their consequences arewritten in dimensionless form with c = 1, so that k = ω . Model domains are as follows: Q = (cid:8) ( x, y, z ) : x + y ≤ , | z | ≤ (cid:9) , X = (cid:8) ( x, y, z ) : x + y ≤ , | z | (cid:9) ,Y = (cid:8) ( x, y, z ) : 3 ≤ x + y , | z | ≤ (cid:9) . It is also supposed that σ ( x ) = 0 and model sources are given in the form f ( x ) = M s X m =1 A m δ ( x − x m ) , where x m are coordinates of δ -shaped point sources. Then u ( x , ω ) = M s X m =1 A m G ( x , x ′ m , ω ) , and the Fourier transform of this function, ˜ u ( r, ϕ, Ω; ω ), can be calculated in advance.We did not set ourselves the goal of optimizing the number, positions and amplitudes ofsources. In all calculations, it was assumed that M s = 8 , A m = 1 , x m = ( r m , ϕ m , z m ) with r m = 4 .
01 and ϕ m = h , π , − π , π, , π , − π , π i , z m = [ − , − , − , − , , , , . For the numerical study of the proposed algorithms, two direct and inverse model problemsare considered.
The first problem has a solution ξ ( x, y, z ) of the form ξ ( x, y, z ) = ξ ( x, y, z ) + ξ ( x, y, z ) , ( x, y, z ) ∈ Q ; ξ ( x, y, z ) = { A exp {− R ( x, y, z ) } , ( x, y, z ) ∈ Q ; 0 , ( x, y, z ) ∈ Q \ Q } ,ξ ( x, y, z ) = { A exp {− R ( x, y, z ) } , ( x, y, z ) ∈ Q ; 0 , ( x, y, z ) ∈ Q \ Q } , R ( x, y, z ) = 5( x − . + 5 y + 0 . z + 0 . ,R ( x, y, z ) = 5( x + 0 . + 5( y − . + 0 . z − . and Q = (cid:8) ( x, y, z ) : ( x − . + y + 0 . z + 0 . ≤ . (cid:9) Q = (cid:8) ( x, y, z ) : ( x + 0 . + ( y − . + 0 . z − . ≤ . (cid:9) This function simulates small local inhomogeneities of the medium, the position of whichand the corresponding velocity distributions must be found. Algorithm 2 is tuned specif-ically to search for such inhomogeneities. The value A determines the contrast∆ cc = max x ( p − c ξ ( x ) ) − A = 0 . l ∼ . c = 1 , ω = 3: ∆ cc ≫ c lω . Such scatterers are quite common inpractice. The second model task will be described below.The equations (5), (6) were approximated in the domains X, Y by the finite-differencemethod on uniform grids of variables r, r ′ , ϕ, z . The sizes of the grids are determinedby the numbers N r , N r ′ , N ϕ , N z . In the region X , the grid has size N r ′ × N ϕ × N z , andin the region Y , the size is N r × N ϕ × N z . Specific dimensions will be given belowfor each example. Discrete analogues of the functions ˜ G ( r, r ′ , cos ϕ, Ω; ω ), ˜ u ( r, ϕ, Ω; ω ),used in Algorithms 1 and 2, were calculated for the considered frequencies ω from theknown values G ( r, r ′ , cos ϕ, z ; ω ), u ( r, ϕ, z ; ω ) using the fast Fourier transform with grid { Ω ( m ) } N z m =1 in the variable Ω. Fourier series expansions were also implemented using theFFT with n ∈ [0 , N ϕ − Now we present typical results of a numerical study of the iterative process (12),(13) toobtain data for the first inverse problem on grids of size N r = 32 , N r ′ = 33 , N ϕ = 90 , N z =64. Figure 2 shows a comparison of the convergence rate for the process with variousquantities ω = k .The iterations were stopped by the value∆ k ( ω ) = (cid:26)P n (cid:13)(cid:13)(cid:13) u ( k ) n ( r, Ω; ω ) − u ( k − n ( r, Ω; ω ) (cid:13)(cid:13)(cid:13) L { Π } (cid:27) / (cid:26)P n (cid:13)(cid:13)(cid:13) u (0) n ( r, Ω; ω ) (cid:13)(cid:13)(cid:13) L { Π } (cid:27) / ,
10 20 30 40 50 60 70 80 90 100 k -15 -10 -5 ∆ k ( ω ) ω =3 ω =2 ω =110 -13 Figure 2:
Convergence rate of iterations (12), (13) for different ω = k . r φ R e ( w ( ν ) )
601 42 32 0 4 r φ I m ( w ( ν ) ) Figure 3:
Typical data w ( ν ) ( r, ϕ, ω = 3) for solving the inverse problem. when the condition ∆ ν ( ω ) − was satisfied for the iteration number ν . Here Π = { ( r, Ω) ∈ [0 , a ] × R Ω } . After that, using the found set of functions n u ( ν ) n ( r, Ω; ω ) o thefunctions n v ( ν ) n ( r ′ , Ω; ω ) o were calculated by the formula (12) with k = ν . Then, using theformula (14), the functions n w ( ν ) n ( r, Ω; ω ) o were found. Further, they were transformedaccording to (15), in the function ˜ w ( ν ) ( r, ϕ, Ω; ω ). The inverse Fourier transform of thelast function with respect to the variable z , w ( ν ) ( r ′ , ϕ, z ; ω ), represents the data for solvingthe inverse problem. The form of this function, found for ω = k = 3, is shown in Fig.3for z = 0.Further, the data for solving the inverse problem were specified with some perturba-tions that can be interpreted as a measurement error. In our calculations, this was mod-eled by imposing an additive normally distributed pseudo-random noise with zero meanon the function w ( ν ) ( r, ϕ, z ; ω ) so that the resulting approximate function w ( ν ) δ ( r, ϕ, z ; ω )12ould satisfy the condition (cid:13)(cid:13)(cid:13) w ( ν ) δ ( r, ϕ, z ; ω ) − w ( ν ) ( r, ϕ, z ; ω ) (cid:13)(cid:13)(cid:13) L ( Y ) ≤ δ (cid:13)(cid:13) w ( ν ) ( r, ϕ, z ; ω ) (cid:13)(cid:13) L ( Y ) . This corresponds to approximate data with relative error δ . The first step of Algorithm 2, i.e. solving equations of the first kind (16) by regularizationmethods was discussed in the papers [17],[18] in connection with the solution of anotherinverse problem, formally similar to the considered one and differing from it only in theform of the kernel and the right-hand side. In these papers, it was noted that for eachconsidered frequency ω , the discretization used reduces the equations (16) to a system oflinear algebraic equations (SLAE) of the form A ( m ) n V ( m ) n = W ( m ) n for each n ∈ [0 , N ϕ − m ∈ [1 , N z ]. Here A ( m ) n = (cid:2) µ ij G n ( r i , r ′ j , Ω ( m ) ; ω ) (cid:3) N r ,N r ′ i =1 ,j =1 is the matrix of the system ob-tained by discretizing the kernel of the equation (16) on the considered grid of size N r × N r ′ ,quantities Ω ( m ) are grid points along Ω, and µ ij are quadrature coefficients for calculatingintegrals in (16). The right sides of the system, W ( m ) n = πω (cid:2) w n (cid:0) r i , Ω ( m ) ; ω (cid:1)(cid:3) N r i =1 , arecolumn vectors of height N r , and the column vector V ( m ) n of height N ′ r contains unknowns v n ( r ′ j , Ω ( m ) ; ω ). Thus, when performing Step 1 of Algorithm 2 for each ω it is necessaryto solve N ϕ × N z systems of linear equations with matrices of size N r × N r ′ .We solved the indicated SLAEs using various versions of the Tikhonov regularizationmethod [22],[23] and using the TSVD method [24]. The justification for these methodsis given in [17],[18] for similar problem. The best calculation results were obtained usingthe TSVD method. We present them further.Step 2 of Algorithm 2 does not cause any difficulties for the discretized problem,since is reduced to matrix multiplication of discrete quantities G n , v n and addition of theresult with a discrete analogue of the function u n . Step 3 was performed using the inverseFFT. Finally, Step 4 was implemented for each considered frequency ω using the followingprocedure for finding the normal pseudosolution of the equation uξ = v by the TSVDmethod at each point ( r ′ , ϕ ′ , z ) ∈ X : ξ = (cid:8) vu , | u | > tol; 0 , | u | ≤ tol (cid:9) with tol = 10 − .Further, it is easy to recalculate the function ξ ( x ) into c ( x ). For the sake of brevity, wedo not do this in the examples below, presenting the value ξ ( x ) directly in the figures. We start by applying Algorithm 2 in numerical experiments of solving the first modelinverse problem on grids of size N r = 32 , N r ′ = 33 , N ϕ = 90 , N z = 64 for ω = k = 3.The problem was solved with exact data, more precisely, calculated by Algorithm 1, andapproximate data with different error levels δ . For a qualitative comparison, Fig.4 showsthe exact and approximate solutions of the inverse problem, ξ exact ( x, y, z ) and ξ appr ( x, y, z ),in different sections with z = const. The exact solution is shown in the first line of thefigure. The second line contains approximate solution for exact data. The third line showsthe solution for the perturbed data with δ = 10 − . The figure demonstrates a fairly high13igure 4: Qualitative comparison of the exact solution ξ exact ( x, y, z ) and the calculated approx-imate solutions ξ appr ( x, y, z ) of the inverse problem in different sections z = const. A) exactsolution; B) approximate solution obtained for exact data by Algorithm 2; C) approximatesolution for disturbed data with δ = 10 − . -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 z ∆ L ( z ) δ =0 δ =10 -10 δ =10 -9 δ =10 -8 Figure 5:
The first model problem for ω = 3. Relative error ∆ L ( z ) of approximate solutionsfor different z at different levels of data perturbation δ . z ∆ L ( z ) δ = 0 δ = 10 -11 δ = 10 -10 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 z ∆ L ( z ) δ =0 δ =10 -10 δ =10 -9 δ =10 -8 Figure 6:
The first model problem. Relative error ∆ L ( z ) of approximate solutions for different z at different levels of data perturbation δ . Left: for ω = 1; right: for ω = 2. Figure 7:
Second model problem. A qualitative comparison of the positions and geometry ofthe reconstructed inhomogeneity ξ ( x ) for different δ . N r T I P ( N r , N φ ) , s e c . N φ =45N φ =90N φ =180 Figure 8:
Time T IP ( N r , N ϕ ) of solving the inverse problem for different N r , N ϕ . sensitivity of solutions to data disturbances. More detailed information on the accuracyof solving the inverse problem, i.e. on the relative error∆ L ( z ) = k ξ appr ( x, y, z ) − ξ exact ( x, y, z ) k L ( X xy ) max z k ξ exact ( x, y, z ) k L ( X xy ) of the approximate solution ξ appr ( x, y, z ) in the sections z = const, is presented in Fig.5for different levels δ of data error.For comparison, Fig.6 shows the errors of approximate solutions obtained for different δ in solving model problems with ω = 1 and ω = 2. The improvement in accuracy isobvious with increasing number ω .When solving the inverse problem under consideration, it is very important to knowhow accurately the algorithm allows one to determine the positions of the investigatedlocal scatterers. For illustration, the second model problem was solved, which differs fromthe first one only in expressions for ξ ( x, y, z ) and ξ ( x, y, z ): ξ ( x, y, z ) = { A ; ( x, y, z ) ∈ Q ; 0 , ( x, y, z ) ∈ Q \ Q } ξ ( x, y, z ) = { A ; ( x, y, z ) ∈ Q ; 0 , ( x, y, z ) ∈ Q \ Q } , A = 0 . . Such a solution corresponds to two ellipsoidal scatterers lying in the X region and filledwith substances having different refractive indices. Fig.7 shows qualitatively the influenceof the problem data perturbation on the determination of the position and geometry ofinhomogeneities. The contrast of the exact solution is 0.291. It is seen that the positionscan be determined quite accurately in the case of δ = 10 − . This is also possible for δ = 10 − when using appropriate noise filtering in the found solution.16 Some properties of Algorithm 2
All calculations were carried out in MATLAB on a PC with an Intel (R) Core (TM) i7-7700CPU 3.60 GHz, 16GB RAM without parallelization. Algorithm 2 for solving the inverseproblem turned out to be fast enough. We present the results of corresponding numericalexperiments for solving the first inverse problem with ω = 2. In the experiments, itwas assumed that the z grid is fixed ( N z = 64), and only the sizes N r , N ϕ of the r and ϕ grids change. In addition, it was assumed that N r ′ = N r + 1. Then the time forsolving the inverse problem is a function of the form T IP ( N r , N ϕ ). This dependence isshown in Fig.8. When the grid size in the variable z is changed, the time T IP ( N r , N ϕ )changes proportionally to the number N z , since the time is determined by solving N z × N ϕ equations of the form (16).Note once again that the inverse problem being solved is very sensitive to input dataerrors. When solving it with double precision, introducing random errors with an ampli-tude of the order of 10 − into the right-hand sides of Eqs.(16) leads to serious distortions ofthe solution. This happens when using the TSVD method and the regularization methodtoo. This sensitivity is associated with a very fast decrease in the singular numbers of thematrices A ( m ) n for SLAEs solved in Step 1 of Algorithm 2, and this is a specific feature ofthe inverse problem being solved. A similar property of the inverse coefficient problem forthe wave equation was noted before in the works [17, 18]. The corresponding theoreticalestimates of the error under various a priori assumptions on the exact solution can befound in [4, 5]. From the numerical experiments carried out in this work, the following conclusions canbe drawn.1. The considered three-dimensional inverse problem of scalar acoustics in a cylindricalregion can be solved using Algorithm 2 for sufficiently fine grids in a few tens of secondson a PC of average performance even without parallelization. For this, one should usethe scheme for recording data of the inverse problem in a cylindrical layer indicated inthe article. The proposed algorithm can be easily parallelized.2. The inverse problem under consideration is in itself very sensitive to data pertur-bations; to obtain a detailed approximate solution, data measured with high accuracy arerequired. This feature of the problem does not depend on the used algorithm.3. Algorithm 2 makes it possible to reliably determine the position and shape of smalllocal inhomogeneities of the acoustic medium with data having small errors.
Acknowledgments
This work was supported by the Russian Science Foundation (project 20-11-20085)for the first author in part of substantiating numerical algorithms and the Programm ofCompetitiveness Increase of the National Research Nuclear University MEPhI (MoscowEngineering Physics Institute); contract no. 02.a03.21.0005, 27.08.2013 for the secondauthor. 17 eferences [1] Ramm A.G.
Multidimensional Inverse Scattering Problems,
Pitman Monogr. Surv.Pure Appl. Math. 51. Harlow: Longman Scientific & Technical, 1992.[2] Colton D., Kress R.
Inverse Acoustic and Electromagnetic Scattering Theory,
Inverse scattering problems in acoustics,
M., Publish-ing House of Moscow State University, 1989 (in Russian).[4] Bakushinsky A., Goncharsky A.
Ill-Posed Problems: Theory and Applications.
Dor-drecht: Kluwer Academic Publishers, 1994.[5] Bakushinsky A.B., Kokurin M.Yu.
Iterative methods for approximate solution of in-verse problems,
Mathematics and Its Applications. Dordrecht: Kluwer Academic Pub-lishers, 2004.[6] A. V. Goncharsky and S. Y. Romanov, On two approaches to the solution of coefficientinverse problems for wave equations, Zh. Vychisl. Mat. Mat. Fiz. 52 (2012), no.2,pp.263-269.[7] A. V. Goncharsky and S. Y. Romanov, Supercomputer technologies in inverse prob-lems of ultrasound tomography, Inverse Problems 29 (2013), no.7, Article ID 075004.[8] Belishev M.I. Recent progress in the boundary control method, Inverse Problems.2007. V.23. N5. P.1-67.[9] Pestov L.N., Bolgova V.M., Danilin A.N. Numerical reconstruction of the threedi-mensional speed of sound by the method of boundary control, Bulletin of Ugra StateUniversity, 2011, Issue 3, pp.92-98 (in Russian).[10] Burov V.A., Alekseenko N.V., Rumyantseva O.D. Multifrequency Generalization ofthe Novikov Algorithm for the Two-Dimensional Inverse Scattering Problem, AcousticJournal, V.55, No.6, 2009, pp.784-798.[11] Novikov P. G. Reconstruction of the two-dimensional Schrodinger operator from thescattering amplitude at a fixed energy, Funktsional. analysis and its adj., T.20, No.3,1986, pp.90-91.[12] Burov V.A., Vecherin S.N., Morozov S.A., Rumyantseva O.D. Modeling of the ExactSolution of the Inverse Scattering Problem by Functional Methods. Acoustic Journal,Vol.56, No.4, 2010, pp.516-536.[13] Beilina L., Klibanov M.V.,
Approximate Global Convergence and Adaptivity for Co-efficient Inverse Problems,
New York: Springer, 2012.[14] Kabanikhin S.I., Satybaev A.D., Shishlenin M.A.
Direct Methods of Solving Multidi-mensional Inverse Hyperbolic Problems,
Utrecht: VSP, 2004.1815] Klibanov M.V., Kolesov A.E. Convexification of a 3-D coefficient inverse scatteringproblem // Computers and Mathematics with Applications. 2019. V.77. P.1681-1702.[16] Klibanov M.V., Kolesov A.E., Nguyen Dinh-Liem. Convexification method for aninverse scattering problem and its performance for experimental backscatter data forburied targets, SIAM J. Imaging Sciences. 2019. V.12, N.1. P.576-603.[17] Bakushinsky A.B., Leonov A.S. Fast numerical method of solving 3D coefficientinverse problem for wave equation with integral data, Journal of Inverse and Ill-PosedProblems. 2018. V.26. Issue 4. P.477-492.[18] A.B. Bakushinskii, and A.S. Leonov, Low-Cost Numerical Method for Solving a Co-efficient Inverse Problem for the Wave Equation in Three-Dimensional Space. Comp.Math. and Math. Phys., 2018, Vol. 58, No. 4, pp.548-561.[19] Evstigneev R.O., Medvedik M.Yu., Smirnov Yu.G., Tsupak A.A. The inverse problemof body’s heterogeneity recovery for early diagnostics of diseases using microwave to-mography, University proceedings, Volga region, Physical and Mathematical Sciences,2017, No.4 (44), pp. 3-17 (in Russian).[20] B.M. Budak, A.A. Samarskii and A.N. Tikhonov,
A Collection of Problems on Math-ematical Physics,
Pergamon Press, Oxford, 1964[21] F. Riesz and B. Sz.-Nagy,
Functional Analysis,
Frederick Ungar Publishing Co., N.Y.,1955.[22] A. N. Tikhonov, A. V. Goncharsky, V. V. Stepanov and A. G. Yagola,
NumericalMethods for the Solution of Ill-Posed Problems , Math. Appl. 328, Kluwer AcademicPublishers, Dordrecht, 1995.[23] A. S. Leonov,
Solution of Ill-Posed Inverse Problems. Theory Review, Practical Algo-rithms and MATLAB Demonstrations , Librokom, Moscow, 2010, 2013 (in Russian).[24] Engl H.W., Hanke M, Neubauer A.