- Email: [email protected]
A BEM based methodology to solve inverse problems considering fictitious background media
Engineering Analysis with Boundary Elements 66 (2016) 109–118
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
A BEM based methodology to solve inverse problems considering fictitious background media Markcilei Lima Dan a,n, Webe João Mansur b,c, Franciane Conceição Peters b,c a
Federal Institute of Espírito Santo, Mechanical Engineering Department, 29300-970 Cachoeiro de Itapemirim, Espírito Santo, Brazil Federal University of Rio de Janeiro, COPPE/Department of Civil Engineering, 21945-970 Rio de Janeiro, Brazil c Lamemo Laboratory/COPPE/Federal University of Rio de Janeiro, Av. Pedro Calmon s/n, Anexo ao Centro de Tecnologia, Ilha do Fundão, 21941-596 Rio de Janeiro, RJ, Brazil b
art ic l e i nf o
a b s t r a c t
Article history: Received 26 March 2015 Received in revised form 29 November 2015 Accepted 24 January 2016 Available online 2 March 2016
This paper proposes a new BEM based inversion methodology named Fictitious Background Media Inverse Formulation (FBMIF) to solve inverse problems described in terms of the Helmholtz equation for inhomogeneous media. The BEM based approach considers a Green’s function related to a homogeneous background medium, a fictitious medium, which works as a transfer function. Thereafter, considering a nonrestrictive linearization scheme, the FBMIF is completely defined. Three examples are used to test the efficacy of the proposed model of inversion, which can be extended to other mathematical models different from those governed by the Helmholtz equation. & 2016 Elsevier Ltd. All rights reserved.
Keywords: Boundary Element Method Helmholtz equation Inverse problems Inhomogeneous media Variable acoustic velocity profile
1. Introduction The problem of imaging internal body structures from external measurements often arises in some branches of science such as medicine, geology and engineering. Medical radiography has been one of the pioneer techniques used. This technique consists of estimating the internal structure of a body from the projection of its shade on a special photosensitive film. To illuminate the internal structure of opaque objects some kind of penetrating radiation, such as X-Ray, has to be used. Although somewhat archaic, this technique still has extensive use in both medical and industrial applications. The natural evolution of radiography, a pioneer inversion technique, led to the development of the Computerized Tomography (CT). CT is an effective imaging technique; however, it makes use of considerable charges of radiation. Thus, developing imaging techniques based on the use of some type of less aggressive energy becomes necessary. In this sense, some researchers have dedicated their efforts to establish imaging techniques based on the measurement of the scattered acoustical field generated when beams of ultrasonic waves reach regions of non-homogeneities inside the studied bodies. As a result, the so-called Diffraction Tomography (DT) was developed [1–3]. However, most of the studies developed n
Corresponding author. E-mail addresses: [email protected] (M.L. Dan), [email protected] (W.J. Mansur), [email protected] (F.C. Peters). http://dx.doi.org/10.1016/j.enganabound.2016.01.011 0955-7997/& 2016 Elsevier Ltd. All rights reserved.
in this sense were devoted to the medical area. Aiming to expand the field of application of the DT to geophysical problems, in 1984, Devaney [4] presented an elegant and low computational cost formulation, named by the author Geophysical Diffraction Tomography (GDT). Although Devaney's proposition has been successfully applied to imaging geophysical profiles, it is mainly applicable to geophysical problems characterized by the existence of an inner, non-homogeneous region embedded in a homogeneous infinite background medium. This restriction is intrinsic to the GDT, as its proposition was built based on the Lippmann–Schwinger Integral Equation (LSIE). The LSIE models the scattered acoustical field due to the incidence of wave fronts on non-homogeneous regions embedded in an infinite homogeneous background medium [5]. In addition, LSIE becomes nonlinear when it is applied to the solution of inverse problems. Thus, it is also necessary to apply some linearization scheme to LSIE before using the GDT algorithm. The most common schemes used to linearize LSIE are the Born and the Rytov approximations. The Born approximation represents well problems whose velocity fields are characterized by weak scatterers; whereas the Rytov approximation, although allowing scatterers stronger than the Born's approximation, is indicated for problems whose velocity fields are characterized by low gradients [1,6]. Besides the need of having to consider an infinite homogeneous background medium, the application of linearization schemes brings extra restrictions to GDT’s fields of application. To deal with the problem of the existence of a homogeneous background medium in geophysical applications, some researchers have proposed variations of
110
M.L. Dan et al. / Engineering Analysis with Boundary Elements 66 (2016) 109–118
the GDT capable of consider a stratified background medium [7,8]. However, these variations still employ Born or Rytov linearization schemes. A robust inversion approach based on LSIE should make use of nonrestrictive linearization methods. In addition to this, it has to overcome the requirement of the existence of a homogeneous real background medium. Considering geophysical applications, a robust inversion scheme based on LSIE can be developed in terms of iterative methods. However, the use of iterative methods (i.e., optimization methods) to solve inverse problems also requires the construction of an effective solution scheme to solve direct problems. In this sense and observing that LSIE is closely related to the integral equation of the Boundary Element Method (BEM) for the Helmholtz equation for inhomogeneous media, Fu [9] proposed an integral sentence named by the author Generalized Lippmann– Schwinger Integral Equation (GLSIE). The GLSIE proposed by Fu was obtained through physical considerations concerning scattered waves occurring in an inhomogeneous region that is embedded in a semi-infinite homogeneous background medium and it also incorporated boundary effects occurring on the border of the nonhomogeneous region. However, GLSIE was applied only to solve direct problems. Thus, there is a vast field of research in the sense of establishing inversion schemes based on GLSIE. However, such propositions have also to overcome the problems of nonlinearity and the one related to the existence of a homogeneous real background medium. Nonlinearity of inverse problems approaches can be well dealt with by iterative solution schemes [10–19]. However, the difficulty related to the consideration of homogeneous background medium remains. With respect to the proposition of direct inversion schemes based on GLSIE, the application of some linearization scheme, as in the case of the Diffraction Tomography, is necessary. However, this linearization scheme has to be nonrestrictive. In this sense, this paper proposes a direct BEM inversion approach called Fictitious Background Media Inverse Formulation (FBMIF). This proposition overcomes the problem of considering a homogeneous background medium and proposes a nonrestrictive linearization scheme able to deal with the nonlinearity of the problem. Applied to the Helmholtz equation for inhomogeneous media, the FBMIF leads to a BEM sentence equivalent to that of GLSIE. However, FBMIF states that it is possible to consider a fictitious background medium in order to obtain a numerical solution for a given problem and that this numerical solution has to be independent of the fictitious background medium chosen. To deal with the nonlinearities, FBMIF proposes a simple and nonrestrictive change of variables. The collocation method is used to obtain a linear system that is assumed to have full row rank. Then, this linear system is solved considering the Singular Value Decomposition of the coefficient matrix. The mathematical background that enables the use of FBMIF to solve inverse problems is also explored. Finally, the performance of the proposed methodology is tested using three examples.
2. Coefficient inverse problem based on the Helmholtz equation The inverse problems treated in this paper are related to the Helmholtz equation for nonhomogeneous media and can be classified in the theory of the inverse problems as coefficient inverse problems or as inverse medium problems [20]. Before defining such an inverse problem, it is didactic to first define its analog direct problem. To establish a well-posed direct problem (a stationary one), one must define the following:
the domain in which the process is studied; the equation that describes the process; the conditions on the boundary of the domain.
c(x,y)
y
x
Fig. 1. Schematic representation of an acoustic domain Ω and its boundary Γ.
For a nonhomogeneous medium exhibiting a non-homogeneous propagation velocity distribution and a constant density, the Helmholtz equation reads: 2
∇2 p þ k ðx; yÞp ¼ 0
ð1Þ
where kðx; yÞ represents a wavenumber defined as: kðx; yÞ ¼
ω
ð2Þ
cðx; yÞ
ω being an angular frequency. An integration domain Ω, limited by a boundary Γ and characterized by an acoustic velocity field cðx; yÞ is shown in Fig. 1: Hence, a well-posed direct problem is defined if the boundary conditions (natural, essential or mixed ones) are defined on the boundary Γ . However, when the acoustic wave velocity field represented by the function cðx; yÞ is not known, a coefficient inverse problem is defined. In this case, the acoustic velocity field cðx; yÞ forms an extra set of unknowns that has also to be found at the end of the process of solving the inverse problem. In order to determine these unknowns, the formulation of the inverse problem has to be supplied with some additional information about the solution of the problem that, in the present approach comprises both the natural and the essential boundary conditions. Hence, in this paper all of the boundary conditions listed below are considered as known: dp ð3Þ ¼ q ðω; Γ Þ ðNatural Boundary ConditionsÞ dn Γ
pΓ ¼ p ðω; Γ Þ
ðEssential Boundary ConditionsÞ
ð4Þ
In the form described above, the coefficient inverse problem is well defined. However, even if more restrictions were imposed to this problem, it will not imply that the problem is well posed in the Hadamard sense, since almost all inverse problems do not comply with the Hadamard stability condition [20]. In addition, the coefficient inverse problems are usually nonlinear. Hence, as in the case of the inversion procedures based on LSIE, some linearization scheme has to be applied in order to make it possible to solve the problem in a direct form. On the other hand, linearization schemes are generally restrictive because of the approximations introduced by them. Alternatively, linearization schemes can be avoided solving the coefficient inverse problems through an iterative process [21]. In turn, an iterative solution scheme of an inverse problem starts with an initial guess of the velocity field cðx; yÞ in Ω, establishing an iterative process that should converge to the solution of the inverse problem but for a predetermined tolerance. This process requires solving linear systems of equations
M.L. Dan et al. / Engineering Analysis with Boundary Elements 66 (2016) 109–118
at each iteration. Therefore, iterative solution schemes for inverse problems become very expensive in terms of computational cost. Thus, there is a gap to overcome in the sense of developing direct solution schemes for inverse problems that can be used as an alternative to the iterative schemes. Hence, this paper proposes a new inversion technique that is able to solve the coefficient inverse problems in a direct form. The proposed methodology leads to a system of linear equations able to recover the pressure field as well as to estimate the acoustic velocity profile inside the acoustical domain Ω. The proposed method is based on a BEM integral sentence called, in the context of this article, Fictitious Background Media Inverse Formulation (FBMIF).
3. The FBMIF for Helmholtz coefficient inverse problem
2
2
2
ð5Þ
where k0 corresponds to the wavenumber of some homogeneous medium characterized by a constant acoustic velocity profile c0: k0 ¼
ω
ð6Þ
c0
It is worth noting that the homogeneous medium corresponding to the wavenumber k0 does not exist in fact. Thus, this homogeneous medium is fictitious. Eq. (5) can be easily rewritten as shown below: h i 2 2 2 ∇2 p þ k0 p ¼ k ðx; yÞ k0 p ð7Þ Hence, making: h i 2 2 kv ¼ k ðx; yÞ k0
2
∇ p þ k0 p ¼ kv p
where k0 corresponds to the wavenumber of the fictitious medium. As can be extensively found in the literature, Eq. (10) is satisfied by the following function: - i Ψ ¼ H20 k0 j r r ξ j ð11Þ 4 where H 20 is the Hankel function of the second kind and zero order. Eq. (11) corresponds to the Green’s function of the fictitious background medium considered in this paper.
Now, considering Eq. (11) and following common procedures used to establish integral sentences for BEM, the following integral equation can be obtained from Eq. (9): Z Z Z c ξ p ξ qΨ dΓ þ pq dΓ ¼ kv pΨ dΩ ð12Þ Γ
ð9Þ
Eq. (9) corresponds exactly to Eq. (1). However, in the form presented above, Eq. (9) leads to some important conclusions about the physical problem analyzed. First, the pressure field modeled by Eq. (1) can be mathematically represented as the response of a homogeneous fictitious medium subjected to a domain source given by kvp, where kv is given by Eq. (8). Second, the pressure field that obeys Eq. (9) is independent of the fictitious background medium considered. The third conclusion is that the fictitious background medium is actually responsible for transmitting the effects of the domain sources ( kvp) to all points of the domain Ω. These three observations, together with the consideration of a Green’s function related to the fictitious background medium, lead to the proposition of the FBMIF. 3.1. Green’s function of the fictitious background medium Any integral sentence of the BEM is directly related to a Green’s function, which is directly related to the physical problem under consideration. However, taking into account the second observation made in the last paragraph of Section 2, one can state that it is not necessary to relate the Green's function to the properties of a real physical medium. Thus, the Green's function used here corresponds to a Green's function of a fictitious medium. Thereby,
Γ
Ω
∂p Ψ and c ξ are constants that depend on the , q ¼ ∂∂n where q ¼ ∂n location of the source point ξ [22]. For kðx; yÞ known (or, in an equivalent way, cðx; yÞ), Eq. (12) can be directly applied to solve direct boundary value problems. However, in the case of inverse problems, kðx; yÞ is not known “a priori”. Thus, kv becomes an additional unknown of the inverse problem and Eq. (12) becomes non-linear. To eliminate the nonlinearity in Eq. (12), a simple and non-restrictive change of variable is proposed as follows: p^ ¼
kv p
ð13Þ
α
Substituting (13) in Eq. (12), the linearity of the integral Eq. (12) ^ is obtained in terms of the new unknown p: Z Z Z αp^ Ψ dΩ ð14Þ c ξ p ξ qΨ dΓ þ pq dΓ ¼ Γ
ð8Þ
one can write: 2
considering an infinite medium characterized by the acoustic wave velocity c0 and subjected to a domain point source δ ξ; X , the correlated Helmholtz problem is defined by the equation below: 2 ∇2 Ψ þ k0 Ψ ¼ δ ξ; r ð10Þ
3.2. Integral formulation
The methodology presented in this paper can be extended to other mathematical models. Thus, it is convenient to explain in detail the procedure considered to obtain the integral sentence used to establish the proposed FBMIF. The main idea considered here consists of modifying the original mathematical model of the problem under study by adding to and subtracting from Eq. (1) the 2 term k0 p . So, one obtains: ∇2 p þ k ðx; yÞp k0 p þ k0 p ¼ 0
111
Γ
Ω
Eq. (14) corresponds to the integral sentence for the proposed FBMIF. The parameter α was introduced in Eq. (14) because, α o1 has a regularization effect on the inverse solution obtained by FBMIF. 3.3. Matrix formulation Since there is a domain integral in Eq. (14), the FBMIF has to consider some BEM technique able to deal with domain integrals. As p^ carries valuable information on the medium properties, it is advisable to use some discretization technique that permits to recover the unknown p^ at discrete points inside the domain. Thus, the method of cell domain discretization [23–25] is considered here, the domain Ω (shown in Fig. 1) being discretized by linear cells. The boundary discretization is carried out using quadratic boundary elements, the discretized mesh being such that each functional cell node that falls on the boundary Γ coincides with one functional node of a boundary element. Hence, if a total of N functional nodes are used in the boundary discretization, a total of N þ M functional nodes results from the cell domain discretization, with a total of M nodes inside domain Ω. Hence, considering the collocation points coincident with the N þ M functional nodes resulting from the problem discretization, Eq. (14) gives rise to the following linear system: HpΓ GqΓ þ IpΩ FFp^ ¼ 0
ð15Þ
where the matrices H and G are the typical ones of the BEM presenting orders ðN þ MÞxN. Matrix FF is obtained from the
112
M.L. Dan et al. / Engineering Analysis with Boundary Elements 66 (2016) 109–118
domain integral that appears in Eq. (14) and presents order ðN þ MÞxðN þ MÞ. As one has M collocation points inside the domain, the identity matrix I that corresponds to the free term c ξ , first term on the l.h.s. of Eq. (14), presents order MxM. Now, applying the boundary conditions to Eq. (15), and remembering that for the coefficient inverse problem defined in Section 2 both the natural and the essential boundary conditions are considered as known, one obtains: IpΩ FFp^ ¼ b1
ð16Þ
where b1 ¼ HpΓ þ GqΓ
ð17Þ
Eq. (16) corresponds to a system of N þ M equations and N þ 2M unknowns. There is a total of N unknowns p^ corresponding to the cell nodes that fall on the boundary and a total of M unknowns p^ corresponding to the cell nodes inside domain Ω. The remainder M unknowns correspond to the nodal values of pΩ for the collocation points inside domain Ω. Thus, if the values of kv on boundary Γ are known, a total of N unknowns corresponding to the values of p^ on boundary Γ can be also eliminated (p^ Γ ¼ kv p Γ ), resulting in the following linear system: Ax ¼ b2
ð18Þ
where the matrix A is a rectangular matrix of order ðN þ MÞ 2M and the vector b2 now incorporates the values of p^ Γ . 3.4. Inverse solution by FBMIF The inverse solution procedure considered in this paper was developed for a problem discretization that obeys the following condition: 2M 4 ðN þ MÞ (it generally occurs). Such discretization condition obviously implies that the linear system (18) admits infinite solutions. In this situation, it is also reasonable to admit that the rows of the matrix A are linearly independent. It means that the matrix A is a full row rank matrix and admits a right inverse that leads to a vector solution xR with special features. Details of the properties of the vector solution xR will be discussed in a more convenient moment. At this moment, we only illustrate how such a solution can be obtained. Thus, considering the following change of variable: xR ¼ AT y
ð19Þ
Eq. (18) can be written in terms of the new unknown vector y as: AAT y ¼ b2
ð20Þ
Considering that matrix A is a full row rank matrix, matrix AAT represents an invertible matrix of order ðN þMÞxðN þMÞ. Hence, Eq. (20) can be promptly solved for the unknowns y: h i1 b2 ð21Þ y ¼ AAT Finally, considering Eq. (19), the solution xR is obtained multiplying Eq. (21) by AT : h i1 b2 ð22Þ xR ¼ AT AAT Vector xR corresponds to the inverse solution obtained by FBMIF. However, the solution scheme shown above has to be avoided because matrix A is generally ill conditioned, although this fact does not invalidate the hypothesis of full row rank. Notice that the same vector solution xR can be obtained by a more efficient solution scheme that considers the generalized inverse A þ obtained from the Singular Value Decomposition (SVD) of the
matrix A: xR ¼ A þ b2
ð23Þ
given that: h i1 A þ ¼ AT AAT
ð24Þ
for a full row rank matrix A [26]. In addition, the generalized inverse A þ can properly deal with some unexpected linear dependence that can occur because of numerical issues, and it still allows the application of some regularization techniques. Finally, one should notice that the vector solution xR contains nodal values of the unknowns pΩ and p^ inside domain. Hence, it is possible to recover discrete values of cðx; yÞ by the following equation:
ω
ci ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 αp^ p i þ k0
ð25Þ
where Eq. (25) was obtained from Eq. (13) by considering Eq. (8) and the wavenumber definition. The index i is used to refer to ^ and cðx; yÞ inside domain. Hence, the the nodal values of pΩ , p, inverse solution of FBMIF is completely defined by Eqs. (23) (or (22)) and (25). 3.5. Mathematical interpretation of the FBMIF solution As mentioned before, it is possible to find an infinite set of vectors x that fits exactly to the linear system represented by Eq. (18). This infinite set of solutions is represented by the following general form: xG ¼ xRow þ xNull
ð26Þ
where the general solution xG represents all vectors x that exactly fit to the linear system (18); the vector xRow represents the row space solution admitted by the linear system (18) and xNull represents all vectors belonging to the null space of the matrix A. Notice that any vector xNull belonging to the null space of the matrix A exactly fits to equation: AxNull ¼ 0
ð27Þ
Then, multiplying Eq. (26) by A, one obtains: AxRow ¼ b2
ð28Þ
given that AxG ¼ b2 . The row space and the null space being mutually orthogonal, xRow actually represents the component of the vector xG that exactly fits to the vector b2 . By definition, it is known that the row space of the matrix A is obtained as a linear combination of the rows of the matrix A. Thus, it is clear that Eq. (19) defines a vector xR belonging to the row space of the matrix A. In addition, the vector xR is uniquely defined by Eq. (19), since AT represents a full column rank matrix of order 2MxðN þ MÞ, with 2M 4 ðN þ MÞ. Hence, it can be concluded that FBMIF establishes a methodology to find the unique solution xR that belongs to the row space solution of the matrix A and that exactly fits the problem conditions represented by the vector b2 . However, the exact solution of the inverse problem can have some components in the null space of A. This fact can reduce drastically the efficiency of the inverse methodology proposed. But this inconvenience is solved when, artificially, one reduces the importance of the null space of the matrix A. In this sense, we note that the dimension of the null space of A can be obtained by the following equation: Null Space Dimension of A ¼ No of Column of A Rank of A ð29Þ
M.L. Dan et al. / Engineering Analysis with Boundary Elements 66 (2016) 109–118
Hence, for a full row rank matrix A one obtains: Null Space Dimension of A ¼ 2M ðN þ MÞ ¼ ðM NÞ
ð30Þ
Eq. (30) suggests that the null space of A becomes more significant as the number M of domain nodes increases over the number N of functional nodes used on the boundary discretization. However, issues related to the boundary or domain discretization are of minor importance. Most important factors are the ones related to the domain geometry. Hence, Eq. (30) actually suggests that the null space of A becomes less significant as slimmer domains are considered. Agreeing with the comment above, it was observed that for large domains the efficiency of the proposed methodology deteriorates, although it has efficiently solved some inverse problems defined in large domains for a very low frequency range. However, it has not been possible yet to define geometrical parameters able to characterize domains that are well dealt with by the proposed inversion methodology, as extensive tests have to be conducted in order to achieve such purpose. This subject, as well as the application of the proposed methodology to analyze practical problems, will be treated in future works. Now, the main objective is to depict the performance of the proposed methodology considering only inverse problems defined in slim domains.
4. Numerical examples
^ On the (and undesirable) signal change to the auxiliary variable p. other hand, the proposition of Eq. (31) can be justified by the third observation made in the last paragraph presented just before the Subsection 3.1, as the consideration of bigger fictitious acoustical velocities, i.e., c0 Z cMax , means that a fictitious background medium presenting suitable acoustical properties will be efficient in transmitting information from point to point inside the domain. 4.1. Example 1 The first example examined corresponds to a Helmholtz problem defined in an inhomogeneous medium characterized by the following linear varying acoustic velocity profile: cðxÞ ¼ x þ 1:5
ð32Þ
The analytical solution considered for this direct Helmholtz problem is: ffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffi 1 2 pðxÞ ¼ x þ1:5e i ω 4 lnðx þ 1:5Þ ð33Þ Eq. (33) represents a propagating harmonic wave for ω 4 1=2. To define an analog inverse problem, we only need to consider the acoustic velocity profile given by Eq. (32) as unknown , the boundary conditions necessary to correctly model the inverse problem being promptly obtained from Eq. (33):
pffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffi 1 2 pΓ ¼ x þ 1:5e i ω 4 lnðx þ 1:5Þ ðEssential BCÞ ð34aÞ Γ
This section presents three examples that illustrate FBMIF performance when solving inverse problems. Each one of them has an analog direct problem with either analytical or numerical solutions known. These direct solutions are used to generate the total set of boundary conditions needed to correctly model the inverse problems under consideration. Then, the simulations of the inverse problems are performed considering that the acoustic velocity field cðx; yÞ inside domain is not known. In the first example, a continuous velocity profile with linear variation is assumed. This example has an analytical solution known and the problem is characterized by a smooth gradient of velocity. Aiming to illustrate that the proposed inversion methodology is not restrictive concerning the acoustic velocity profile considered, the second example uses a discontinuous acoustic velocity field that characterizes a two-layered stratified medium. These examples are defined in a regular integration domain, as shown in Fig. 2, and the acoustical velocity profile varies only with the x coordinate (c ¼ cðxÞ). In the third example, an irregular Lshaped integration domain and a two-dimensional variable velocity profile (c ¼ cðx; yÞ) are considered. Finally, it has to be pointed out that the following criterion was used to set values for the fictitious wavenumber k0: k0 rkMin
113
ð31Þ
ω is the smallest physical wavenumber present in where kMin ¼ cMax the hmedium considered. In this sense and recalling that i 2 2 kv ¼ k ðx; yÞ k0 , it is observed that Eq. (31) prevents artificial
Fig. 2. Regular acoustic domain considered in the first and the second examples.
qffiffiffiffiffiffiffiffiffiffiffiffiffi 1 3 1 1 pffiffiffiffiffiffiffiffiffiffi 2 1 ∂p 2 2 þi ω 4 i ω 4 lnðx þ 1:5Þ A dx 5 ¼ 4@ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e ∂nΓ dn x þ 1:5 20
ðNatural BCÞ Γ
ð34bÞ The inverse problem considered in this subsection is completely defined. For numerical simulations, the integration domain presented in Fig. 2 was discretized using 44 quadratic boundary elements and 160 linear internal cells. Fig. 3 shows the mesh used in the numerical simulation. The numerical results shown below were obtained considering two different situations. The first one is represented by Eq. (18) and corresponds to the situation where nodal values of cðxÞ (or, equivalently, kv) are previously known on the boundary Γ . The second one seeks for an inverse problem solution, considering only the natural and the essential boundary conditions, as stated in Eq. (16). For both cases, the parameter α (introduced by Eq. (13)) is equal to 0:01 and the following values of angular frequency were considered: ω ¼ 6rad=s (Fig. 4) and ω ¼ 2rad=s (Fig. 5). The results presented in Figs. 4 and 5 correspond to the nodal values of acoustic velocity recovered inside the domain obtained from Eqs. (18) or (25), considering the application of the zero order Tikhonov regularization technique for a regularization parameter αT ¼ 10 5 . Figs. 6 and 7 show the nodal values of acoustic velocity recovered inside the domain (by Eq. (25)) considering numerical simulations performed from Eq. (16). In this case, it was possible to estimate the existent acoustic velocity field in Ω without considering the application of any regularization technique. The results shown above demonstrate the efficiency of FBMIF in solving this inverse problem, as the acoustic velocity profile could be efficiently estimated even without considering extra “boundary conditions” represented by prescribed values of p^ on the boundary Γ (p^ Γ ). Figs. 4 and 5 show an excellent agreement between the estimated and the exact acoustic velocity profile, although they represent regularized solutions obtained considering prescribed values of p^ Γ as it is stated in Eq. (18). On the other hand, maintaining the coefficients relative to the unknowns p^Γ in the coefficient matrix, as considered in Eq. (16), the conditioning of the coefficient matrix is favored. Thus, it was possible to solve the
114
M.L. Dan et al. / Engineering Analysis with Boundary Elements 66 (2016) 109–118
Fig. 3. Boundary and domain discretization considered in the first and the second examples.
Acoustic Velocity (m/s)
14.0 Calc.
Exact
10.5
7.0
3.5
0.0
88
127
166
205
Internal Node
Fig. 4. Nodal values of acoustic velocity recovered inside domain Ω (see Fig. 3) obtained from Eq. (18) (with p^ prescribed on the boundary Γ) considering a fictitious wavenumber k0 ¼ 0:40m 1 and an angular frequency ω ¼ 6rad=s.
14.0
Acoustic Velocity (m/s)
Calc.
Exact
10.5
7.0
3.5
0.0 88
127
166
205
Internal Node
Fig. 5. Nodal values of acoustic velocity recovered inside domain Ω (see Fig. 3) obtained from Eq. (18) (with p^ prescribed on the boundary Γ) considering a fictitious wavenumber k0 ¼ 0:15m 1 and an angular frequency ω ¼ 2rad=s.
14.0
Acoustic Velocity (m/s)
Calc.
Exact
10.5
7.0
3.5
0.0 88
127
166
205
Internal Node
Fig. 6. Nodal values of acoustic velocity recovered inside domain Ω (see Fig. 3) obtained from Eq. (16) (without prescribing p^ on the boundary Γ) considering a fictitious wavenumber k0 ¼ 0:40m 1 and an angular frequency ω ¼ 6rad=s.
M.L. Dan et al. / Engineering Analysis with Boundary Elements 66 (2016) 109–118
115
14.0 Acoustic Velocity (m/s)
Calc.
Exact
10.5
7.0
3.5
0.0 88
127
166
205
Internal Node
Fig. 7. Nodal values of acoustic velocity recovered inside domain Ω (see Fig. 3) obtained from Eq. (16) (without prescribing p^ on boundary Γ) considering a fictitious wavenumber k0 ¼ 0.15 m 1 and an angular frequency ω ¼ 2rad=s.
inverse problem without considering the application of any regularization technique as shown in Figs. 6 and 7. It also means that Eq. (16) meets the hypothesis of full row rank in a stronger form than Eq. (18), although it was not observed that this condition was violated by Eq. (18). Despite being possible to observe that the results obtained without considering prescribed values of p^ on the boundary Γ are slightly worse than the ones obtained considering p^Γ prescribed, it is noted that the FBMIF efficiently solved this inverse problem, even when the natural boundary conditions and the essential boundary conditions are all the information available. In the same way, significant influence of the angular frequency ω in the inverse solution obtained was not observed in the frequency range tested. 4.2. Example 2 In the second example, we analyze a Helmholtz problem defined in an inhomogeneous medium characterized by a discontinuous acoustic velocity profile as shown in Fig. 8. For numerical simulations the domain discretization was performed considering the mesh shown in Fig. 3, since the integration domain is the same shown in Fig. 2. This Helmholtz problem has no known analytical solution. Thus, the boundary conditions necessary to correctly define the inverse problem have to be numerically generated. In this sense, a direct problem is defined by arbitrarily imposing the following boundary conditions: dp ¼ 20 þ 5 i ð35aÞ dnΓ 1 dp ¼ 100 10 i dnΓ 2
ð35bÞ
dp ¼ 30 þ 40 i dnΓ 3
ð35cÞ
dp ¼ 80 10 i dnΓ 4
ð35dÞ
Hence, assuming an angular frequency ω ¼ 50rad=s, the essential boundary conditions can be generated. Again, the numerical simulations were performed considering Γ . As preor not considering prescribed values of p^ on boundary viously done in the case of prescribed values of p^ Γ , the numerical solution was regularized considering the application of the zero order Tikhonov regularization technique for a regularization parameter αT ¼ 10 5 . In the other case, no regularization technique was used. In both cases, the parameter α was equal to 0:01.
Fig. 8. Integration domain and discontinuous acoustic velocity profile considered in the second example.
Fig. 9 presents the nodal values of the acoustic velocity field recovered inside domain Ω considering prescribed values of p^ on boundary Γ . In Fig. 10 we observe again that, even without considering prescribed values of p^ on boundary Γ , it was possible to estimate the nodal values cðxÞ inside domain Ω with reasonable accuracy. However, the results presented in Fig. 9 are higher. This is because the prescription of p^ on boundary Γ favors the representation of the inverse problem solution, mainly concerning the recovered values of the variable p^ inside the domain. On the other hand, this gain of representativeness causes the impoverishment of the coefficient matrix, thus the application of regularization techniques is needed. When the coefficients relative to the unknown p^ are kept in the coefficient matrix, the system conditioning is favored, allowing to obtain a non-regularized solution for the inverse problem. However, the information omitted ^ resulting in a decreases the representativeness of the unknown p, slightly worse, although still consistent, estimative for the acoustical velocity profile. In this sense, we observe that FBMIF arises as an efficient alternative tool to simulate inverse problems for slim domains. Furthermore, it has to be noted that this example considers a two layered physical medium characterized by a discontinuous acoustic velocity profile. It is important to notice that a nonrestrictive linearization scheme has to be able to recover efficiently acoustic profiles characterized by severe gradients, or discontinuous ones in extreme cases, as well as for acoustic fields characterized by smooth behavior like the one shown in the previous example. Hence, the results shown above lead to the conclusion that the linearization scheme represented by Eq. (13) is nonrestrictive, as it was possible to recover the nodal values of velocity with excellent agreement, even for a discontinuous velocity profile. 4.3. Example 3 This example considers a Helmholtz problem defined in an Lshaped domain with an acoustical velocity profile cðx; yÞ varying in both x and y directions in a quadratic form as shown in Fig. 11.
116
M.L. Dan et al. / Engineering Analysis with Boundary Elements 66 (2016) 109–118 350.0
Acoustic Velocity (m/s)
Exact
Calc.
300.0
250.0
200.0
150.0 88
127
166
205
Internal Node
Fig. 9. Nodal values of acoustic velocity recovered inside domain Ω (see Fig. 3) for prescribed values of p^ on boundary Γ, ω ¼ 50rad=s and k0 ¼ 0:09m 1 .
350.0
Acoustic Velocity (m/s)
Exact
Calc.
300.0
250.0
200.0
150.0
88
127
166
205
Internal Node
Fig. 10. Nodal values of acoustic velocity recovered inside domain Ω (see Fig. 3) without considering prescribed values of p^ on boundary Γ for ω ¼ 50rad=s and k0 ¼ 0:09m 1 .
Again, it has to be noted that this Helmholtz problem has no known analytical solution. Thus, the data is the numerical solution of Helmholtz problem for an angular frequency ω ¼ 50 rad=s , considering the following natural boundary conditions: dp ¼ 100 þ 200 i ð36aÞ dnΓ1 dp ¼ 150 þ100 i dnΓ 2
ð36bÞ
dp ¼ 300 90 i dnΓ 3
ð36cÞ
dp dn
¼ 100 þ 150 i
ð36dÞ
dp ¼ 200 100 i dnΓ 5
ð36eÞ
dp dn
Γ4
Γ6
¼ 150 þ200 i
ð36f Þ
Following the same idea outlined in the previous examples, the results presented below consider situations where nodal values of p^ were prescribed and not prescribed on boundary Γ . In the case of prescribed values of p^ Γ , a zero order Tikhonov regularization for a regularization parameter αT ¼ 10 4 was applied. In the other case, no regularization technique was used. Here, the parameter α is equal to 0:01. The boundary and the domain discretization were performed respectively considering 80 quadratic boundary elements and 304 linear internal cells. The node numbering sequence follows the same pattern presented in Fig. 3.
Fig. 11. L-shaped integration domain and acoustical velocity profile cðx; yÞ considered in the third example.
The results presented below show the recovered values of the acoustic velocity field cðx; yÞ at each internal node used in the discretization of the integration domain shown in Fig. 11. The results depicted in Figs. 12 and 13 above correspond respectively to the values of acoustic velocity recovered inside domain Ω with or without considering prescribed values of p^ on boundary Γ . Once more, it is observed that, even in the case of not prescribed values of p^ on boundary Γ, it was possible to estimate the nodal values cðx; yÞ inside domain Ω with reasonable accuracy. Different from the previous two examples, this problem considers both: an acoustical wave propagation velocity field varying as a
M.L. Dan et al. / Engineering Analysis with Boundary Elements 66 (2016) 109–118
117
500.0
Acoustic Velocity (m/s)
Calc.
Exact
400.0
300.0
199
349
275
200.0 319
100.0 244
167
0.0 160
235
310
385
Internal Node
Fig. 12. Nodal values of acoustic velocity recovered inside domain Ω considering prescribed values of p^ on boundary Γ for ω ¼ 50rad=s and k0 ¼ 0:1m 1 . 500.0
Acoustic Velocity (m/s)
Calc.
Exact
400.0
300.0 199
349
275
200.0 319
100.0 244
167
0.0 160
235
310
385
Internal Node
Fig. 13. Nodal values of acoustic velocity recovered inside domain Ω without considering prescribed values of p^ on boundary Γ for ω ¼ 50rad=s and k0 ¼ 0:1m 1 .
quadratic function of both, x and y coordinates and an L-shaped irregular integration domain. Nevertheless, it was observed that these factors have not greatly influenced the performance of FBMIF. For a quantitative analysis, it is pointed out that the results presented in Fig. 12 produce a mean error percentage of 0.32% with peaks of nodal errors of 5.32%, 5.90% and 9.01%, occurring at the nodes 167, 244 and 319 respectively. As the knee of the L occurs at the nodes 199, 275 and 349, these values of error peaks are not due to the domain geometry, but due to a change of sign (introduced by the unknown p) occurring in both real and ima^ With respect to Fig. 13, the mean ginary parts of the variable p. error percentage produced was of 4.90%. The same effects of the change of sign occurring at the nodes 167, 244 and 319 are observed, as well as a slight sensibility to the knee of the L occurring in the neighborhood of the nodes 199, 275 and 349. The effect observed in the neighborhood of the nodes 199, 275 and 349 can also be due to a change of sign occurring in the real part of p^ at the nodes 194, 271 and 346. Anyway, it is clear that this matter is not strongly influencing the value of the mean error observed, as well as the inverse solution obtained by FBMIF was not significantly influenced by the domain geometry considered in this example, even when not prescribed values of p^ on boundary Γ were considered.
5. Conclusion This paper proposes a new BEM based inversion methodology named FBMIF. This proposition was only possible because it was noted that a numerical solution of an inverse problem could be obtained by considering a Green’s function related to fictitious background media. Thereafter, the proposal of a nonrestrictive linearization scheme led to the establishment of the FBMIF. To test
the efficacy of the model of inversion proposed, three examples of inverse problems were considered. The first one was defined from a direct problem with a known analytical solution. This example is mainly characterized by a low gradient of velocity. In the second example, a discontinuous case of the acoustic velocity profile was considered. This example can be thought as an extreme case for high gradient of velocity. Then, it was possible to confirm that the proposed linearization scheme is really a nonrestrictive one. In this sense, there is no evidence that some velocity profile represented by a continuous function with smooth behavior could introduce restriction to the application of the inversion scheme proposed. This observation can be extended even to discontinuous functions. However, the same observation cannot be applied to the integration domain because the matrix obtained by the FBMIF admits a nontrivial null space. As it was noted, the influence of this null space on the FBMIF solution is closely related to the domain geometry. Although the FBMIF has been proposed with a clear mathematical background, it has not yet been possible to define geometrical parameters able to characterize a domain geometry that can be well treated by the proposed inversion methodology. Therefore, in the absence of more conclusive criteria, it can be primarily affirmed that the proposed methodology is indicated to inverse problems defined on slim domains. To enforce this conclusion, a slim L-shaped integration domain was considered. It was not observed that the geometry of this L-shaped slim domain has influenced the FBMIF inverse solution. Finally, it has to be pointed out that FBMIF also admits a nonhomogeneous fictitious background medium, as long as the Green’s function related to the medium considered can be determined. However, this matter as well as issues related to the geometrical restriction to the application of the FBMIF will be treated in future works.
118
M.L. Dan et al. / Engineering Analysis with Boundary Elements 66 (2016) 109–118
Acknowledgments To the National Council for Scientific and Technological Development (CNPq) (Grant no. 306933/2014-4), the Research Foundation of the State of Rio de Janeiro (FAPERJ) (Grant no. E/26/ 201.453/2014), and PETROBRAS for their support (Grant no. 0050.0070743.11.9).
References [1] Kaveh M, Mueller RK, Iverson RD. Ultrasonic tomography based on perturbation solutions of the wave equation. Comput. Gr. Image Process. 1979;9 (2):105–16. [2] Mueller RK. Diffraction tomography I: the wave-equation. Ultrason. Imaging 1980;2(3):213–22. [3] Devaney AJ. A filtered backpropagation algorithm for diffraction tomography. Ultrason. Imaging 1982;4:336–50. [4] Devaney AJ. Geophysical diffraction tomography. IEEE Trans. on Geosci. Remote Sens. 1984;GE-22(n. 1):3–13. [5] Tien-When L, Inderwiesen PL. Fundamentals of Seismic Tomography. Tulsa: Society of Exploration Geophysicists; 1994. [6] Devaney AJ. Generalized projection-slice theorem for fan beam diffraction tomography. Ultrason. Imaging 1985;7(3):264–75. [7] Devaney AJ, Zhang D. Geophysical diffraction tomography in a layered background. Wave Motion 1991;14:243–65. [8] Harris JM, Wang GY. Diffraction tomography for inhomogeneities in layered background medium. Geophysics 1996;61(2):570–83. [9] Fu LY. Numerical study of generalized lipmann-schwinger integral equation including surface topography. Geophysics 2003;68(2):665–71. [10] Tarantola A. Inversion of seismic reflection data in the acoustic approximation. Geophysics 1984;49(8):1259–66.
[11] Kolb P, Collino F, Lailly P. Pre-stack inversion of a 1-D medium. Proc. IEEE 1986;74(3):498–508. [12] Virieux J. P-SV wave propagation in heterogeneous media: velocity-stress finite-difference method. Geophysics 1986;51(4):889–901. [13] Mora P. Nonlinear two-dimensional elastic inversion of multioffset seismic data. Geophysics 1987;52(9):1211–28. [14] Crase E, Pica A, Noble M, et al. Robust Elastic nonlinear waveform inversion: application to real data. Geophysics 1990;55(5):527–38. [15] Pica A, Diet JP, Tarantola A. Nonlinear inversion of seismic reflection data in a laterally invariant medium. Geophysics 1990;55(3):284–92. [16] Xu T, McMechan GA, Sun R. 3-D prestack full-wavefield inversion. Geophysics 1995;60(6):1805–18. [17] Djikpéssé HA, Tarantola A. Multiparameter ℓ1 norm waveform fitting: interpretation of gulf of mexico reflection seismograms. Geophysics 1999;64 (4):1023–35. [18] Kleefeld A, Reißel M. The Levenberg–Marquardt method applied to a parameter estimation problem arising from electrical resistivity tomography. Appl. Math. Comput. 2011;217(9):4490–501. [19] Tao Y, Sen MK. Frequency-domain full waveform inversion with a scatteringintegral approach and its sensitivity analysis. J. Geophys. Eng. 2013;10(6). [20] Kabanikhin SI. Definitions and examples of inverse and Ill-posed problems. J. Inverse Ill-Posed Probl. 2008;16(4):317–57. [21] Pratt RG, Shin C, Hicks GJ. Gauss–Newton and Full Newton Methods in frequency-space seismic waveform inversion. Geophys. J. Int. 1998;133 (2):341–62. [22] Wu TW. Boundary Element Acoustics: Fundamentals and Computer Codes. Boston: WIT Press; 2000. [23] Brebbia CA, Skerget P. Diffusion-convection problems using boundary elements. Adv. Water Resour. 1984;7(2):50–7. [24] Onishi K, Kuroki T, Tanaka M. An application of a boundary element method to natural convection. Appl. Math. Model. 1984;8(6):383–90. [25] Takhteyev V, Brebbia CA. Analytical integrations in boundary elements. Eng. Anal. Bound Elem. 1990;7(2):95–100. [26] Aster RC, Borchers B, Thurber CH. Parameter Estimation and Inverse Problems. London, UK: Elsevier Academic Press; 2005.
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
A BEM based methodology to solve inverse problems considering fictitious background media Markcilei Lima Dan a,n, Webe João Mansur b,c, Franciane Conceição Peters b,c a
Federal Institute of Espírito Santo, Mechanical Engineering Department, 29300-970 Cachoeiro de Itapemirim, Espírito Santo, Brazil Federal University of Rio de Janeiro, COPPE/Department of Civil Engineering, 21945-970 Rio de Janeiro, Brazil c Lamemo Laboratory/COPPE/Federal University of Rio de Janeiro, Av. Pedro Calmon s/n, Anexo ao Centro de Tecnologia, Ilha do Fundão, 21941-596 Rio de Janeiro, RJ, Brazil b
art ic l e i nf o
a b s t r a c t
Article history: Received 26 March 2015 Received in revised form 29 November 2015 Accepted 24 January 2016 Available online 2 March 2016
This paper proposes a new BEM based inversion methodology named Fictitious Background Media Inverse Formulation (FBMIF) to solve inverse problems described in terms of the Helmholtz equation for inhomogeneous media. The BEM based approach considers a Green’s function related to a homogeneous background medium, a fictitious medium, which works as a transfer function. Thereafter, considering a nonrestrictive linearization scheme, the FBMIF is completely defined. Three examples are used to test the efficacy of the proposed model of inversion, which can be extended to other mathematical models different from those governed by the Helmholtz equation. & 2016 Elsevier Ltd. All rights reserved.
Keywords: Boundary Element Method Helmholtz equation Inverse problems Inhomogeneous media Variable acoustic velocity profile
1. Introduction The problem of imaging internal body structures from external measurements often arises in some branches of science such as medicine, geology and engineering. Medical radiography has been one of the pioneer techniques used. This technique consists of estimating the internal structure of a body from the projection of its shade on a special photosensitive film. To illuminate the internal structure of opaque objects some kind of penetrating radiation, such as X-Ray, has to be used. Although somewhat archaic, this technique still has extensive use in both medical and industrial applications. The natural evolution of radiography, a pioneer inversion technique, led to the development of the Computerized Tomography (CT). CT is an effective imaging technique; however, it makes use of considerable charges of radiation. Thus, developing imaging techniques based on the use of some type of less aggressive energy becomes necessary. In this sense, some researchers have dedicated their efforts to establish imaging techniques based on the measurement of the scattered acoustical field generated when beams of ultrasonic waves reach regions of non-homogeneities inside the studied bodies. As a result, the so-called Diffraction Tomography (DT) was developed [1–3]. However, most of the studies developed n
Corresponding author. E-mail addresses: [email protected] (M.L. Dan), [email protected] (W.J. Mansur), [email protected] (F.C. Peters). http://dx.doi.org/10.1016/j.enganabound.2016.01.011 0955-7997/& 2016 Elsevier Ltd. All rights reserved.
in this sense were devoted to the medical area. Aiming to expand the field of application of the DT to geophysical problems, in 1984, Devaney [4] presented an elegant and low computational cost formulation, named by the author Geophysical Diffraction Tomography (GDT). Although Devaney's proposition has been successfully applied to imaging geophysical profiles, it is mainly applicable to geophysical problems characterized by the existence of an inner, non-homogeneous region embedded in a homogeneous infinite background medium. This restriction is intrinsic to the GDT, as its proposition was built based on the Lippmann–Schwinger Integral Equation (LSIE). The LSIE models the scattered acoustical field due to the incidence of wave fronts on non-homogeneous regions embedded in an infinite homogeneous background medium [5]. In addition, LSIE becomes nonlinear when it is applied to the solution of inverse problems. Thus, it is also necessary to apply some linearization scheme to LSIE before using the GDT algorithm. The most common schemes used to linearize LSIE are the Born and the Rytov approximations. The Born approximation represents well problems whose velocity fields are characterized by weak scatterers; whereas the Rytov approximation, although allowing scatterers stronger than the Born's approximation, is indicated for problems whose velocity fields are characterized by low gradients [1,6]. Besides the need of having to consider an infinite homogeneous background medium, the application of linearization schemes brings extra restrictions to GDT’s fields of application. To deal with the problem of the existence of a homogeneous background medium in geophysical applications, some researchers have proposed variations of
110
M.L. Dan et al. / Engineering Analysis with Boundary Elements 66 (2016) 109–118
the GDT capable of consider a stratified background medium [7,8]. However, these variations still employ Born or Rytov linearization schemes. A robust inversion approach based on LSIE should make use of nonrestrictive linearization methods. In addition to this, it has to overcome the requirement of the existence of a homogeneous real background medium. Considering geophysical applications, a robust inversion scheme based on LSIE can be developed in terms of iterative methods. However, the use of iterative methods (i.e., optimization methods) to solve inverse problems also requires the construction of an effective solution scheme to solve direct problems. In this sense and observing that LSIE is closely related to the integral equation of the Boundary Element Method (BEM) for the Helmholtz equation for inhomogeneous media, Fu [9] proposed an integral sentence named by the author Generalized Lippmann– Schwinger Integral Equation (GLSIE). The GLSIE proposed by Fu was obtained through physical considerations concerning scattered waves occurring in an inhomogeneous region that is embedded in a semi-infinite homogeneous background medium and it also incorporated boundary effects occurring on the border of the nonhomogeneous region. However, GLSIE was applied only to solve direct problems. Thus, there is a vast field of research in the sense of establishing inversion schemes based on GLSIE. However, such propositions have also to overcome the problems of nonlinearity and the one related to the existence of a homogeneous real background medium. Nonlinearity of inverse problems approaches can be well dealt with by iterative solution schemes [10–19]. However, the difficulty related to the consideration of homogeneous background medium remains. With respect to the proposition of direct inversion schemes based on GLSIE, the application of some linearization scheme, as in the case of the Diffraction Tomography, is necessary. However, this linearization scheme has to be nonrestrictive. In this sense, this paper proposes a direct BEM inversion approach called Fictitious Background Media Inverse Formulation (FBMIF). This proposition overcomes the problem of considering a homogeneous background medium and proposes a nonrestrictive linearization scheme able to deal with the nonlinearity of the problem. Applied to the Helmholtz equation for inhomogeneous media, the FBMIF leads to a BEM sentence equivalent to that of GLSIE. However, FBMIF states that it is possible to consider a fictitious background medium in order to obtain a numerical solution for a given problem and that this numerical solution has to be independent of the fictitious background medium chosen. To deal with the nonlinearities, FBMIF proposes a simple and nonrestrictive change of variables. The collocation method is used to obtain a linear system that is assumed to have full row rank. Then, this linear system is solved considering the Singular Value Decomposition of the coefficient matrix. The mathematical background that enables the use of FBMIF to solve inverse problems is also explored. Finally, the performance of the proposed methodology is tested using three examples.
2. Coefficient inverse problem based on the Helmholtz equation The inverse problems treated in this paper are related to the Helmholtz equation for nonhomogeneous media and can be classified in the theory of the inverse problems as coefficient inverse problems or as inverse medium problems [20]. Before defining such an inverse problem, it is didactic to first define its analog direct problem. To establish a well-posed direct problem (a stationary one), one must define the following:
the domain in which the process is studied; the equation that describes the process; the conditions on the boundary of the domain.
c(x,y)
y
x
Fig. 1. Schematic representation of an acoustic domain Ω and its boundary Γ.
For a nonhomogeneous medium exhibiting a non-homogeneous propagation velocity distribution and a constant density, the Helmholtz equation reads: 2
∇2 p þ k ðx; yÞp ¼ 0
ð1Þ
where kðx; yÞ represents a wavenumber defined as: kðx; yÞ ¼
ω
ð2Þ
cðx; yÞ
ω being an angular frequency. An integration domain Ω, limited by a boundary Γ and characterized by an acoustic velocity field cðx; yÞ is shown in Fig. 1: Hence, a well-posed direct problem is defined if the boundary conditions (natural, essential or mixed ones) are defined on the boundary Γ . However, when the acoustic wave velocity field represented by the function cðx; yÞ is not known, a coefficient inverse problem is defined. In this case, the acoustic velocity field cðx; yÞ forms an extra set of unknowns that has also to be found at the end of the process of solving the inverse problem. In order to determine these unknowns, the formulation of the inverse problem has to be supplied with some additional information about the solution of the problem that, in the present approach comprises both the natural and the essential boundary conditions. Hence, in this paper all of the boundary conditions listed below are considered as known: dp ð3Þ ¼ q ðω; Γ Þ ðNatural Boundary ConditionsÞ dn Γ
pΓ ¼ p ðω; Γ Þ
ðEssential Boundary ConditionsÞ
ð4Þ
In the form described above, the coefficient inverse problem is well defined. However, even if more restrictions were imposed to this problem, it will not imply that the problem is well posed in the Hadamard sense, since almost all inverse problems do not comply with the Hadamard stability condition [20]. In addition, the coefficient inverse problems are usually nonlinear. Hence, as in the case of the inversion procedures based on LSIE, some linearization scheme has to be applied in order to make it possible to solve the problem in a direct form. On the other hand, linearization schemes are generally restrictive because of the approximations introduced by them. Alternatively, linearization schemes can be avoided solving the coefficient inverse problems through an iterative process [21]. In turn, an iterative solution scheme of an inverse problem starts with an initial guess of the velocity field cðx; yÞ in Ω, establishing an iterative process that should converge to the solution of the inverse problem but for a predetermined tolerance. This process requires solving linear systems of equations
M.L. Dan et al. / Engineering Analysis with Boundary Elements 66 (2016) 109–118
at each iteration. Therefore, iterative solution schemes for inverse problems become very expensive in terms of computational cost. Thus, there is a gap to overcome in the sense of developing direct solution schemes for inverse problems that can be used as an alternative to the iterative schemes. Hence, this paper proposes a new inversion technique that is able to solve the coefficient inverse problems in a direct form. The proposed methodology leads to a system of linear equations able to recover the pressure field as well as to estimate the acoustic velocity profile inside the acoustical domain Ω. The proposed method is based on a BEM integral sentence called, in the context of this article, Fictitious Background Media Inverse Formulation (FBMIF).
3. The FBMIF for Helmholtz coefficient inverse problem
2
2
2
ð5Þ
where k0 corresponds to the wavenumber of some homogeneous medium characterized by a constant acoustic velocity profile c0: k0 ¼
ω
ð6Þ
c0
It is worth noting that the homogeneous medium corresponding to the wavenumber k0 does not exist in fact. Thus, this homogeneous medium is fictitious. Eq. (5) can be easily rewritten as shown below: h i 2 2 2 ∇2 p þ k0 p ¼ k ðx; yÞ k0 p ð7Þ Hence, making: h i 2 2 kv ¼ k ðx; yÞ k0
2
∇ p þ k0 p ¼ kv p
where k0 corresponds to the wavenumber of the fictitious medium. As can be extensively found in the literature, Eq. (10) is satisfied by the following function: - i Ψ ¼ H20 k0 j r r ξ j ð11Þ 4 where H 20 is the Hankel function of the second kind and zero order. Eq. (11) corresponds to the Green’s function of the fictitious background medium considered in this paper.
Now, considering Eq. (11) and following common procedures used to establish integral sentences for BEM, the following integral equation can be obtained from Eq. (9): Z Z Z c ξ p ξ qΨ dΓ þ pq dΓ ¼ kv pΨ dΩ ð12Þ Γ
ð9Þ
Eq. (9) corresponds exactly to Eq. (1). However, in the form presented above, Eq. (9) leads to some important conclusions about the physical problem analyzed. First, the pressure field modeled by Eq. (1) can be mathematically represented as the response of a homogeneous fictitious medium subjected to a domain source given by kvp, where kv is given by Eq. (8). Second, the pressure field that obeys Eq. (9) is independent of the fictitious background medium considered. The third conclusion is that the fictitious background medium is actually responsible for transmitting the effects of the domain sources ( kvp) to all points of the domain Ω. These three observations, together with the consideration of a Green’s function related to the fictitious background medium, lead to the proposition of the FBMIF. 3.1. Green’s function of the fictitious background medium Any integral sentence of the BEM is directly related to a Green’s function, which is directly related to the physical problem under consideration. However, taking into account the second observation made in the last paragraph of Section 2, one can state that it is not necessary to relate the Green's function to the properties of a real physical medium. Thus, the Green's function used here corresponds to a Green's function of a fictitious medium. Thereby,
Γ
Ω
∂p Ψ and c ξ are constants that depend on the , q ¼ ∂∂n where q ¼ ∂n location of the source point ξ [22]. For kðx; yÞ known (or, in an equivalent way, cðx; yÞ), Eq. (12) can be directly applied to solve direct boundary value problems. However, in the case of inverse problems, kðx; yÞ is not known “a priori”. Thus, kv becomes an additional unknown of the inverse problem and Eq. (12) becomes non-linear. To eliminate the nonlinearity in Eq. (12), a simple and non-restrictive change of variable is proposed as follows: p^ ¼
kv p
ð13Þ
α
Substituting (13) in Eq. (12), the linearity of the integral Eq. (12) ^ is obtained in terms of the new unknown p: Z Z Z αp^ Ψ dΩ ð14Þ c ξ p ξ qΨ dΓ þ pq dΓ ¼ Γ
ð8Þ
one can write: 2
considering an infinite medium characterized by the acoustic wave velocity c0 and subjected to a domain point source δ ξ; X , the correlated Helmholtz problem is defined by the equation below: 2 ∇2 Ψ þ k0 Ψ ¼ δ ξ; r ð10Þ
3.2. Integral formulation
The methodology presented in this paper can be extended to other mathematical models. Thus, it is convenient to explain in detail the procedure considered to obtain the integral sentence used to establish the proposed FBMIF. The main idea considered here consists of modifying the original mathematical model of the problem under study by adding to and subtracting from Eq. (1) the 2 term k0 p . So, one obtains: ∇2 p þ k ðx; yÞp k0 p þ k0 p ¼ 0
111
Γ
Ω
Eq. (14) corresponds to the integral sentence for the proposed FBMIF. The parameter α was introduced in Eq. (14) because, α o1 has a regularization effect on the inverse solution obtained by FBMIF. 3.3. Matrix formulation Since there is a domain integral in Eq. (14), the FBMIF has to consider some BEM technique able to deal with domain integrals. As p^ carries valuable information on the medium properties, it is advisable to use some discretization technique that permits to recover the unknown p^ at discrete points inside the domain. Thus, the method of cell domain discretization [23–25] is considered here, the domain Ω (shown in Fig. 1) being discretized by linear cells. The boundary discretization is carried out using quadratic boundary elements, the discretized mesh being such that each functional cell node that falls on the boundary Γ coincides with one functional node of a boundary element. Hence, if a total of N functional nodes are used in the boundary discretization, a total of N þ M functional nodes results from the cell domain discretization, with a total of M nodes inside domain Ω. Hence, considering the collocation points coincident with the N þ M functional nodes resulting from the problem discretization, Eq. (14) gives rise to the following linear system: HpΓ GqΓ þ IpΩ FFp^ ¼ 0
ð15Þ
where the matrices H and G are the typical ones of the BEM presenting orders ðN þ MÞxN. Matrix FF is obtained from the
112
M.L. Dan et al. / Engineering Analysis with Boundary Elements 66 (2016) 109–118
domain integral that appears in Eq. (14) and presents order ðN þ MÞxðN þ MÞ. As one has M collocation points inside the domain, the identity matrix I that corresponds to the free term c ξ , first term on the l.h.s. of Eq. (14), presents order MxM. Now, applying the boundary conditions to Eq. (15), and remembering that for the coefficient inverse problem defined in Section 2 both the natural and the essential boundary conditions are considered as known, one obtains: IpΩ FFp^ ¼ b1
ð16Þ
where b1 ¼ HpΓ þ GqΓ
ð17Þ
Eq. (16) corresponds to a system of N þ M equations and N þ 2M unknowns. There is a total of N unknowns p^ corresponding to the cell nodes that fall on the boundary and a total of M unknowns p^ corresponding to the cell nodes inside domain Ω. The remainder M unknowns correspond to the nodal values of pΩ for the collocation points inside domain Ω. Thus, if the values of kv on boundary Γ are known, a total of N unknowns corresponding to the values of p^ on boundary Γ can be also eliminated (p^ Γ ¼ kv p Γ ), resulting in the following linear system: Ax ¼ b2
ð18Þ
where the matrix A is a rectangular matrix of order ðN þ MÞ 2M and the vector b2 now incorporates the values of p^ Γ . 3.4. Inverse solution by FBMIF The inverse solution procedure considered in this paper was developed for a problem discretization that obeys the following condition: 2M 4 ðN þ MÞ (it generally occurs). Such discretization condition obviously implies that the linear system (18) admits infinite solutions. In this situation, it is also reasonable to admit that the rows of the matrix A are linearly independent. It means that the matrix A is a full row rank matrix and admits a right inverse that leads to a vector solution xR with special features. Details of the properties of the vector solution xR will be discussed in a more convenient moment. At this moment, we only illustrate how such a solution can be obtained. Thus, considering the following change of variable: xR ¼ AT y
ð19Þ
Eq. (18) can be written in terms of the new unknown vector y as: AAT y ¼ b2
ð20Þ
Considering that matrix A is a full row rank matrix, matrix AAT represents an invertible matrix of order ðN þMÞxðN þMÞ. Hence, Eq. (20) can be promptly solved for the unknowns y: h i1 b2 ð21Þ y ¼ AAT Finally, considering Eq. (19), the solution xR is obtained multiplying Eq. (21) by AT : h i1 b2 ð22Þ xR ¼ AT AAT Vector xR corresponds to the inverse solution obtained by FBMIF. However, the solution scheme shown above has to be avoided because matrix A is generally ill conditioned, although this fact does not invalidate the hypothesis of full row rank. Notice that the same vector solution xR can be obtained by a more efficient solution scheme that considers the generalized inverse A þ obtained from the Singular Value Decomposition (SVD) of the
matrix A: xR ¼ A þ b2
ð23Þ
given that: h i1 A þ ¼ AT AAT
ð24Þ
for a full row rank matrix A [26]. In addition, the generalized inverse A þ can properly deal with some unexpected linear dependence that can occur because of numerical issues, and it still allows the application of some regularization techniques. Finally, one should notice that the vector solution xR contains nodal values of the unknowns pΩ and p^ inside domain. Hence, it is possible to recover discrete values of cðx; yÞ by the following equation:
ω
ci ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 αp^ p i þ k0
ð25Þ
where Eq. (25) was obtained from Eq. (13) by considering Eq. (8) and the wavenumber definition. The index i is used to refer to ^ and cðx; yÞ inside domain. Hence, the the nodal values of pΩ , p, inverse solution of FBMIF is completely defined by Eqs. (23) (or (22)) and (25). 3.5. Mathematical interpretation of the FBMIF solution As mentioned before, it is possible to find an infinite set of vectors x that fits exactly to the linear system represented by Eq. (18). This infinite set of solutions is represented by the following general form: xG ¼ xRow þ xNull
ð26Þ
where the general solution xG represents all vectors x that exactly fit to the linear system (18); the vector xRow represents the row space solution admitted by the linear system (18) and xNull represents all vectors belonging to the null space of the matrix A. Notice that any vector xNull belonging to the null space of the matrix A exactly fits to equation: AxNull ¼ 0
ð27Þ
Then, multiplying Eq. (26) by A, one obtains: AxRow ¼ b2
ð28Þ
given that AxG ¼ b2 . The row space and the null space being mutually orthogonal, xRow actually represents the component of the vector xG that exactly fits to the vector b2 . By definition, it is known that the row space of the matrix A is obtained as a linear combination of the rows of the matrix A. Thus, it is clear that Eq. (19) defines a vector xR belonging to the row space of the matrix A. In addition, the vector xR is uniquely defined by Eq. (19), since AT represents a full column rank matrix of order 2MxðN þ MÞ, with 2M 4 ðN þ MÞ. Hence, it can be concluded that FBMIF establishes a methodology to find the unique solution xR that belongs to the row space solution of the matrix A and that exactly fits the problem conditions represented by the vector b2 . However, the exact solution of the inverse problem can have some components in the null space of A. This fact can reduce drastically the efficiency of the inverse methodology proposed. But this inconvenience is solved when, artificially, one reduces the importance of the null space of the matrix A. In this sense, we note that the dimension of the null space of A can be obtained by the following equation: Null Space Dimension of A ¼ No of Column of A Rank of A ð29Þ
M.L. Dan et al. / Engineering Analysis with Boundary Elements 66 (2016) 109–118
Hence, for a full row rank matrix A one obtains: Null Space Dimension of A ¼ 2M ðN þ MÞ ¼ ðM NÞ
ð30Þ
Eq. (30) suggests that the null space of A becomes more significant as the number M of domain nodes increases over the number N of functional nodes used on the boundary discretization. However, issues related to the boundary or domain discretization are of minor importance. Most important factors are the ones related to the domain geometry. Hence, Eq. (30) actually suggests that the null space of A becomes less significant as slimmer domains are considered. Agreeing with the comment above, it was observed that for large domains the efficiency of the proposed methodology deteriorates, although it has efficiently solved some inverse problems defined in large domains for a very low frequency range. However, it has not been possible yet to define geometrical parameters able to characterize domains that are well dealt with by the proposed inversion methodology, as extensive tests have to be conducted in order to achieve such purpose. This subject, as well as the application of the proposed methodology to analyze practical problems, will be treated in future works. Now, the main objective is to depict the performance of the proposed methodology considering only inverse problems defined in slim domains.
4. Numerical examples
^ On the (and undesirable) signal change to the auxiliary variable p. other hand, the proposition of Eq. (31) can be justified by the third observation made in the last paragraph presented just before the Subsection 3.1, as the consideration of bigger fictitious acoustical velocities, i.e., c0 Z cMax , means that a fictitious background medium presenting suitable acoustical properties will be efficient in transmitting information from point to point inside the domain. 4.1. Example 1 The first example examined corresponds to a Helmholtz problem defined in an inhomogeneous medium characterized by the following linear varying acoustic velocity profile: cðxÞ ¼ x þ 1:5
ð32Þ
The analytical solution considered for this direct Helmholtz problem is: ffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffi 1 2 pðxÞ ¼ x þ1:5e i ω 4 lnðx þ 1:5Þ ð33Þ Eq. (33) represents a propagating harmonic wave for ω 4 1=2. To define an analog inverse problem, we only need to consider the acoustic velocity profile given by Eq. (32) as unknown , the boundary conditions necessary to correctly model the inverse problem being promptly obtained from Eq. (33):
pffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffi 1 2 pΓ ¼ x þ 1:5e i ω 4 lnðx þ 1:5Þ ðEssential BCÞ ð34aÞ Γ
This section presents three examples that illustrate FBMIF performance when solving inverse problems. Each one of them has an analog direct problem with either analytical or numerical solutions known. These direct solutions are used to generate the total set of boundary conditions needed to correctly model the inverse problems under consideration. Then, the simulations of the inverse problems are performed considering that the acoustic velocity field cðx; yÞ inside domain is not known. In the first example, a continuous velocity profile with linear variation is assumed. This example has an analytical solution known and the problem is characterized by a smooth gradient of velocity. Aiming to illustrate that the proposed inversion methodology is not restrictive concerning the acoustic velocity profile considered, the second example uses a discontinuous acoustic velocity field that characterizes a two-layered stratified medium. These examples are defined in a regular integration domain, as shown in Fig. 2, and the acoustical velocity profile varies only with the x coordinate (c ¼ cðxÞ). In the third example, an irregular Lshaped integration domain and a two-dimensional variable velocity profile (c ¼ cðx; yÞ) are considered. Finally, it has to be pointed out that the following criterion was used to set values for the fictitious wavenumber k0: k0 rkMin
113
ð31Þ
ω is the smallest physical wavenumber present in where kMin ¼ cMax the hmedium considered. In this sense and recalling that i 2 2 kv ¼ k ðx; yÞ k0 , it is observed that Eq. (31) prevents artificial
Fig. 2. Regular acoustic domain considered in the first and the second examples.
qffiffiffiffiffiffiffiffiffiffiffiffiffi 1 3 1 1 pffiffiffiffiffiffiffiffiffiffi 2 1 ∂p 2 2 þi ω 4 i ω 4 lnðx þ 1:5Þ A dx 5 ¼ 4@ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e ∂nΓ dn x þ 1:5 20
ðNatural BCÞ Γ
ð34bÞ The inverse problem considered in this subsection is completely defined. For numerical simulations, the integration domain presented in Fig. 2 was discretized using 44 quadratic boundary elements and 160 linear internal cells. Fig. 3 shows the mesh used in the numerical simulation. The numerical results shown below were obtained considering two different situations. The first one is represented by Eq. (18) and corresponds to the situation where nodal values of cðxÞ (or, equivalently, kv) are previously known on the boundary Γ . The second one seeks for an inverse problem solution, considering only the natural and the essential boundary conditions, as stated in Eq. (16). For both cases, the parameter α (introduced by Eq. (13)) is equal to 0:01 and the following values of angular frequency were considered: ω ¼ 6rad=s (Fig. 4) and ω ¼ 2rad=s (Fig. 5). The results presented in Figs. 4 and 5 correspond to the nodal values of acoustic velocity recovered inside the domain obtained from Eqs. (18) or (25), considering the application of the zero order Tikhonov regularization technique for a regularization parameter αT ¼ 10 5 . Figs. 6 and 7 show the nodal values of acoustic velocity recovered inside the domain (by Eq. (25)) considering numerical simulations performed from Eq. (16). In this case, it was possible to estimate the existent acoustic velocity field in Ω without considering the application of any regularization technique. The results shown above demonstrate the efficiency of FBMIF in solving this inverse problem, as the acoustic velocity profile could be efficiently estimated even without considering extra “boundary conditions” represented by prescribed values of p^ on the boundary Γ (p^ Γ ). Figs. 4 and 5 show an excellent agreement between the estimated and the exact acoustic velocity profile, although they represent regularized solutions obtained considering prescribed values of p^ Γ as it is stated in Eq. (18). On the other hand, maintaining the coefficients relative to the unknowns p^Γ in the coefficient matrix, as considered in Eq. (16), the conditioning of the coefficient matrix is favored. Thus, it was possible to solve the
114
M.L. Dan et al. / Engineering Analysis with Boundary Elements 66 (2016) 109–118
Fig. 3. Boundary and domain discretization considered in the first and the second examples.
Acoustic Velocity (m/s)
14.0 Calc.
Exact
10.5
7.0
3.5
0.0
88
127
166
205
Internal Node
Fig. 4. Nodal values of acoustic velocity recovered inside domain Ω (see Fig. 3) obtained from Eq. (18) (with p^ prescribed on the boundary Γ) considering a fictitious wavenumber k0 ¼ 0:40m 1 and an angular frequency ω ¼ 6rad=s.
14.0
Acoustic Velocity (m/s)
Calc.
Exact
10.5
7.0
3.5
0.0 88
127
166
205
Internal Node
Fig. 5. Nodal values of acoustic velocity recovered inside domain Ω (see Fig. 3) obtained from Eq. (18) (with p^ prescribed on the boundary Γ) considering a fictitious wavenumber k0 ¼ 0:15m 1 and an angular frequency ω ¼ 2rad=s.
14.0
Acoustic Velocity (m/s)
Calc.
Exact
10.5
7.0
3.5
0.0 88
127
166
205
Internal Node
Fig. 6. Nodal values of acoustic velocity recovered inside domain Ω (see Fig. 3) obtained from Eq. (16) (without prescribing p^ on the boundary Γ) considering a fictitious wavenumber k0 ¼ 0:40m 1 and an angular frequency ω ¼ 6rad=s.
M.L. Dan et al. / Engineering Analysis with Boundary Elements 66 (2016) 109–118
115
14.0 Acoustic Velocity (m/s)
Calc.
Exact
10.5
7.0
3.5
0.0 88
127
166
205
Internal Node
Fig. 7. Nodal values of acoustic velocity recovered inside domain Ω (see Fig. 3) obtained from Eq. (16) (without prescribing p^ on boundary Γ) considering a fictitious wavenumber k0 ¼ 0.15 m 1 and an angular frequency ω ¼ 2rad=s.
inverse problem without considering the application of any regularization technique as shown in Figs. 6 and 7. It also means that Eq. (16) meets the hypothesis of full row rank in a stronger form than Eq. (18), although it was not observed that this condition was violated by Eq. (18). Despite being possible to observe that the results obtained without considering prescribed values of p^ on the boundary Γ are slightly worse than the ones obtained considering p^Γ prescribed, it is noted that the FBMIF efficiently solved this inverse problem, even when the natural boundary conditions and the essential boundary conditions are all the information available. In the same way, significant influence of the angular frequency ω in the inverse solution obtained was not observed in the frequency range tested. 4.2. Example 2 In the second example, we analyze a Helmholtz problem defined in an inhomogeneous medium characterized by a discontinuous acoustic velocity profile as shown in Fig. 8. For numerical simulations the domain discretization was performed considering the mesh shown in Fig. 3, since the integration domain is the same shown in Fig. 2. This Helmholtz problem has no known analytical solution. Thus, the boundary conditions necessary to correctly define the inverse problem have to be numerically generated. In this sense, a direct problem is defined by arbitrarily imposing the following boundary conditions: dp ¼ 20 þ 5 i ð35aÞ dnΓ 1 dp ¼ 100 10 i dnΓ 2
ð35bÞ
dp ¼ 30 þ 40 i dnΓ 3
ð35cÞ
dp ¼ 80 10 i dnΓ 4
ð35dÞ
Hence, assuming an angular frequency ω ¼ 50rad=s, the essential boundary conditions can be generated. Again, the numerical simulations were performed considering Γ . As preor not considering prescribed values of p^ on boundary viously done in the case of prescribed values of p^ Γ , the numerical solution was regularized considering the application of the zero order Tikhonov regularization technique for a regularization parameter αT ¼ 10 5 . In the other case, no regularization technique was used. In both cases, the parameter α was equal to 0:01.
Fig. 8. Integration domain and discontinuous acoustic velocity profile considered in the second example.
Fig. 9 presents the nodal values of the acoustic velocity field recovered inside domain Ω considering prescribed values of p^ on boundary Γ . In Fig. 10 we observe again that, even without considering prescribed values of p^ on boundary Γ , it was possible to estimate the nodal values cðxÞ inside domain Ω with reasonable accuracy. However, the results presented in Fig. 9 are higher. This is because the prescription of p^ on boundary Γ favors the representation of the inverse problem solution, mainly concerning the recovered values of the variable p^ inside the domain. On the other hand, this gain of representativeness causes the impoverishment of the coefficient matrix, thus the application of regularization techniques is needed. When the coefficients relative to the unknown p^ are kept in the coefficient matrix, the system conditioning is favored, allowing to obtain a non-regularized solution for the inverse problem. However, the information omitted ^ resulting in a decreases the representativeness of the unknown p, slightly worse, although still consistent, estimative for the acoustical velocity profile. In this sense, we observe that FBMIF arises as an efficient alternative tool to simulate inverse problems for slim domains. Furthermore, it has to be noted that this example considers a two layered physical medium characterized by a discontinuous acoustic velocity profile. It is important to notice that a nonrestrictive linearization scheme has to be able to recover efficiently acoustic profiles characterized by severe gradients, or discontinuous ones in extreme cases, as well as for acoustic fields characterized by smooth behavior like the one shown in the previous example. Hence, the results shown above lead to the conclusion that the linearization scheme represented by Eq. (13) is nonrestrictive, as it was possible to recover the nodal values of velocity with excellent agreement, even for a discontinuous velocity profile. 4.3. Example 3 This example considers a Helmholtz problem defined in an Lshaped domain with an acoustical velocity profile cðx; yÞ varying in both x and y directions in a quadratic form as shown in Fig. 11.
116
M.L. Dan et al. / Engineering Analysis with Boundary Elements 66 (2016) 109–118 350.0
Acoustic Velocity (m/s)
Exact
Calc.
300.0
250.0
200.0
150.0 88
127
166
205
Internal Node
Fig. 9. Nodal values of acoustic velocity recovered inside domain Ω (see Fig. 3) for prescribed values of p^ on boundary Γ, ω ¼ 50rad=s and k0 ¼ 0:09m 1 .
350.0
Acoustic Velocity (m/s)
Exact
Calc.
300.0
250.0
200.0
150.0
88
127
166
205
Internal Node
Fig. 10. Nodal values of acoustic velocity recovered inside domain Ω (see Fig. 3) without considering prescribed values of p^ on boundary Γ for ω ¼ 50rad=s and k0 ¼ 0:09m 1 .
Again, it has to be noted that this Helmholtz problem has no known analytical solution. Thus, the data is the numerical solution of Helmholtz problem for an angular frequency ω ¼ 50 rad=s , considering the following natural boundary conditions: dp ¼ 100 þ 200 i ð36aÞ dnΓ1 dp ¼ 150 þ100 i dnΓ 2
ð36bÞ
dp ¼ 300 90 i dnΓ 3
ð36cÞ
dp dn
¼ 100 þ 150 i
ð36dÞ
dp ¼ 200 100 i dnΓ 5
ð36eÞ
dp dn
Γ4
Γ6
¼ 150 þ200 i
ð36f Þ
Following the same idea outlined in the previous examples, the results presented below consider situations where nodal values of p^ were prescribed and not prescribed on boundary Γ . In the case of prescribed values of p^ Γ , a zero order Tikhonov regularization for a regularization parameter αT ¼ 10 4 was applied. In the other case, no regularization technique was used. Here, the parameter α is equal to 0:01. The boundary and the domain discretization were performed respectively considering 80 quadratic boundary elements and 304 linear internal cells. The node numbering sequence follows the same pattern presented in Fig. 3.
Fig. 11. L-shaped integration domain and acoustical velocity profile cðx; yÞ considered in the third example.
The results presented below show the recovered values of the acoustic velocity field cðx; yÞ at each internal node used in the discretization of the integration domain shown in Fig. 11. The results depicted in Figs. 12 and 13 above correspond respectively to the values of acoustic velocity recovered inside domain Ω with or without considering prescribed values of p^ on boundary Γ . Once more, it is observed that, even in the case of not prescribed values of p^ on boundary Γ, it was possible to estimate the nodal values cðx; yÞ inside domain Ω with reasonable accuracy. Different from the previous two examples, this problem considers both: an acoustical wave propagation velocity field varying as a
M.L. Dan et al. / Engineering Analysis with Boundary Elements 66 (2016) 109–118
117
500.0
Acoustic Velocity (m/s)
Calc.
Exact
400.0
300.0
199
349
275
200.0 319
100.0 244
167
0.0 160
235
310
385
Internal Node
Fig. 12. Nodal values of acoustic velocity recovered inside domain Ω considering prescribed values of p^ on boundary Γ for ω ¼ 50rad=s and k0 ¼ 0:1m 1 . 500.0
Acoustic Velocity (m/s)
Calc.
Exact
400.0
300.0 199
349
275
200.0 319
100.0 244
167
0.0 160
235
310
385
Internal Node
Fig. 13. Nodal values of acoustic velocity recovered inside domain Ω without considering prescribed values of p^ on boundary Γ for ω ¼ 50rad=s and k0 ¼ 0:1m 1 .
quadratic function of both, x and y coordinates and an L-shaped irregular integration domain. Nevertheless, it was observed that these factors have not greatly influenced the performance of FBMIF. For a quantitative analysis, it is pointed out that the results presented in Fig. 12 produce a mean error percentage of 0.32% with peaks of nodal errors of 5.32%, 5.90% and 9.01%, occurring at the nodes 167, 244 and 319 respectively. As the knee of the L occurs at the nodes 199, 275 and 349, these values of error peaks are not due to the domain geometry, but due to a change of sign (introduced by the unknown p) occurring in both real and ima^ With respect to Fig. 13, the mean ginary parts of the variable p. error percentage produced was of 4.90%. The same effects of the change of sign occurring at the nodes 167, 244 and 319 are observed, as well as a slight sensibility to the knee of the L occurring in the neighborhood of the nodes 199, 275 and 349. The effect observed in the neighborhood of the nodes 199, 275 and 349 can also be due to a change of sign occurring in the real part of p^ at the nodes 194, 271 and 346. Anyway, it is clear that this matter is not strongly influencing the value of the mean error observed, as well as the inverse solution obtained by FBMIF was not significantly influenced by the domain geometry considered in this example, even when not prescribed values of p^ on boundary Γ were considered.
5. Conclusion This paper proposes a new BEM based inversion methodology named FBMIF. This proposition was only possible because it was noted that a numerical solution of an inverse problem could be obtained by considering a Green’s function related to fictitious background media. Thereafter, the proposal of a nonrestrictive linearization scheme led to the establishment of the FBMIF. To test
the efficacy of the model of inversion proposed, three examples of inverse problems were considered. The first one was defined from a direct problem with a known analytical solution. This example is mainly characterized by a low gradient of velocity. In the second example, a discontinuous case of the acoustic velocity profile was considered. This example can be thought as an extreme case for high gradient of velocity. Then, it was possible to confirm that the proposed linearization scheme is really a nonrestrictive one. In this sense, there is no evidence that some velocity profile represented by a continuous function with smooth behavior could introduce restriction to the application of the inversion scheme proposed. This observation can be extended even to discontinuous functions. However, the same observation cannot be applied to the integration domain because the matrix obtained by the FBMIF admits a nontrivial null space. As it was noted, the influence of this null space on the FBMIF solution is closely related to the domain geometry. Although the FBMIF has been proposed with a clear mathematical background, it has not yet been possible to define geometrical parameters able to characterize a domain geometry that can be well treated by the proposed inversion methodology. Therefore, in the absence of more conclusive criteria, it can be primarily affirmed that the proposed methodology is indicated to inverse problems defined on slim domains. To enforce this conclusion, a slim L-shaped integration domain was considered. It was not observed that the geometry of this L-shaped slim domain has influenced the FBMIF inverse solution. Finally, it has to be pointed out that FBMIF also admits a nonhomogeneous fictitious background medium, as long as the Green’s function related to the medium considered can be determined. However, this matter as well as issues related to the geometrical restriction to the application of the FBMIF will be treated in future works.
118
M.L. Dan et al. / Engineering Analysis with Boundary Elements 66 (2016) 109–118
Acknowledgments To the National Council for Scientific and Technological Development (CNPq) (Grant no. 306933/2014-4), the Research Foundation of the State of Rio de Janeiro (FAPERJ) (Grant no. E/26/ 201.453/2014), and PETROBRAS for their support (Grant no. 0050.0070743.11.9).
References [1] Kaveh M, Mueller RK, Iverson RD. Ultrasonic tomography based on perturbation solutions of the wave equation. Comput. Gr. Image Process. 1979;9 (2):105–16. [2] Mueller RK. Diffraction tomography I: the wave-equation. Ultrason. Imaging 1980;2(3):213–22. [3] Devaney AJ. A filtered backpropagation algorithm for diffraction tomography. Ultrason. Imaging 1982;4:336–50. [4] Devaney AJ. Geophysical diffraction tomography. IEEE Trans. on Geosci. Remote Sens. 1984;GE-22(n. 1):3–13. [5] Tien-When L, Inderwiesen PL. Fundamentals of Seismic Tomography. Tulsa: Society of Exploration Geophysicists; 1994. [6] Devaney AJ. Generalized projection-slice theorem for fan beam diffraction tomography. Ultrason. Imaging 1985;7(3):264–75. [7] Devaney AJ, Zhang D. Geophysical diffraction tomography in a layered background. Wave Motion 1991;14:243–65. [8] Harris JM, Wang GY. Diffraction tomography for inhomogeneities in layered background medium. Geophysics 1996;61(2):570–83. [9] Fu LY. Numerical study of generalized lipmann-schwinger integral equation including surface topography. Geophysics 2003;68(2):665–71. [10] Tarantola A. Inversion of seismic reflection data in the acoustic approximation. Geophysics 1984;49(8):1259–66.
[11] Kolb P, Collino F, Lailly P. Pre-stack inversion of a 1-D medium. Proc. IEEE 1986;74(3):498–508. [12] Virieux J. P-SV wave propagation in heterogeneous media: velocity-stress finite-difference method. Geophysics 1986;51(4):889–901. [13] Mora P. Nonlinear two-dimensional elastic inversion of multioffset seismic data. Geophysics 1987;52(9):1211–28. [14] Crase E, Pica A, Noble M, et al. Robust Elastic nonlinear waveform inversion: application to real data. Geophysics 1990;55(5):527–38. [15] Pica A, Diet JP, Tarantola A. Nonlinear inversion of seismic reflection data in a laterally invariant medium. Geophysics 1990;55(3):284–92. [16] Xu T, McMechan GA, Sun R. 3-D prestack full-wavefield inversion. Geophysics 1995;60(6):1805–18. [17] Djikpéssé HA, Tarantola A. Multiparameter ℓ1 norm waveform fitting: interpretation of gulf of mexico reflection seismograms. Geophysics 1999;64 (4):1023–35. [18] Kleefeld A, Reißel M. The Levenberg–Marquardt method applied to a parameter estimation problem arising from electrical resistivity tomography. Appl. Math. Comput. 2011;217(9):4490–501. [19] Tao Y, Sen MK. Frequency-domain full waveform inversion with a scatteringintegral approach and its sensitivity analysis. J. Geophys. Eng. 2013;10(6). [20] Kabanikhin SI. Definitions and examples of inverse and Ill-posed problems. J. Inverse Ill-Posed Probl. 2008;16(4):317–57. [21] Pratt RG, Shin C, Hicks GJ. Gauss–Newton and Full Newton Methods in frequency-space seismic waveform inversion. Geophys. J. Int. 1998;133 (2):341–62. [22] Wu TW. Boundary Element Acoustics: Fundamentals and Computer Codes. Boston: WIT Press; 2000. [23] Brebbia CA, Skerget P. Diffusion-convection problems using boundary elements. Adv. Water Resour. 1984;7(2):50–7. [24] Onishi K, Kuroki T, Tanaka M. An application of a boundary element method to natural convection. Appl. Math. Model. 1984;8(6):383–90. [25] Takhteyev V, Brebbia CA. Analytical integrations in boundary elements. Eng. Anal. Bound Elem. 1990;7(2):95–100. [26] Aster RC, Borchers B, Thurber CH. Parameter Estimation and Inverse Problems. London, UK: Elsevier Academic Press; 2005.