SciELO - Scientific Electronic Library Online

 
vol.32 número4Dynamic right coprimt factorization and observer design for nonlinear systems índice de autoresíndice de assuntospesquisa de artigos
Home Pagelista alfabética de periódicos  

Serviços Personalizados

Journal

Artigo

Indicadores

  • Não possue artigos citadosCitado por SciELO

Links relacionados

  • Não possue artigos similaresSimilares em SciELO

Compartilhar


Latin American applied research

versão impressa ISSN 0327-0793versão On-line ISSN 1851-8796

Lat. Am. appl. res. v.32 n.4 Bahía Blanca dez. 2002

 

Numerical experiments on parameter estimation in acoustic media using the adjoint method

E. M. Fernández-Berdaguer1, L. V. Perez2 and J. E. Santos3

1 CONICET - Institute de Cáculo, Facultad de Ciencias Exactas y Naturales, Universidad de Buenos Aires
Departamento de Matemática, Facultad de Ingeniería, Universidad de Buenos Aires, Buenos Aires, 1431, Argentina
efernan@ulises.ic.fcen.uba.ar

2 CONICET - Institute de Cálculo, Facultad de Ciencias Exactas y Naturales, Universidad de Buenos Aires, Buenos Aires, 1431, Argentina
lperez@ing.unrc.edu.ar

3 CONICET - Departamento de Geofísica Aplicada, Universidad Nacional de La Plata, La Plata, 1900, Argentina
Department of Mathematics, Purdue University, West Lafayette, IN 47907, USA
santos@fcaglp.fcaglp.unlp.edu.ar

Abstract — We present an algorithm to solve the problem of estimation of the wave speed in a layered acoustic medium. It is based on the adjoint method. Even though it is a fast algorithm, it is able to detect approximately the interfaces and estimate the amplitudes of the speed parameter. Numerical examples illustrate the application of the algorithm in seismic inversion.

Keywords — Inverse Problems. Seismic Inversion. Acoustic Equation. Finite Element Methods.

I. INTRODUCTION

The purpose of this work is to present a fast algorithm to estimate the wave speed in a layered acoustic bounded domain.

In previous papers we treated this estimation problem using a Newton-type method on a quasilinearization algorithm. Newton-type methods are time consuming. In order to improve the performance we present a steepest descent procedure. We use the gradient of the linearized cost functional which is computed from its continuous version. The gradient is calculated using the adjoint of the Gateaux derivative of the solution of the equation with respect to the parameter.

Most of the necessary mathematics is presented in a previous paper (Fernández-Berdaguer, 1998). There the proofs are done for absorbing boundary conditions at the top and bottom of the domain. Here we use more realistic ones: a Dirichlet boundary condition at the top of the boundary, i.e., the air-solid interface and an absorbing boundary condition at the bottom boundary.

The present work describes the implementation of the method for the one-dimensional case in Q = (0,1). This simple case shows the most important aspects of the scheme, for example, the detection of the interfaces and the estimation of the parameter amplitudes, with minor computational burden.

This paper presents a fast estimation algorithm to solve a generalization of the velocity inversion procedure first described in Fernández-Berdaguer et al., 1996. In that paper we considered a layered medium consisting of layers of known thickness. Our method here allows us to attack problems where no a priori information is available about the parameterization of the unknown parameter function. The theory developed in Fernández-Berdaguer (1998), allows for the wave speed c(x) to be an arbitrary function in L(Ω); thus at the discrete level, c(x) may take different values at each point of the spatial discretization.

It is worth noting that even when our work is intended for applications in geophysical exploration, the analysis and algorithms presented here are also of interest in other fields such as non-destructive material testing (Burk and Weiss, 1979), polymer physics (Ferry, 1970) and ocean acoustics (Stoll, 1977).

The organization of the paper is as follows: in section 2 we present the direct model, the inverse problem and the sensitivity equations. In section 3 we state the algorithm and its implementation. Finally, in section 4 we present some numerical examples.

II. THE DIRECT MODEL, THE ESTIMATION PROBLEM AND THE SENSITIVITY EQUATIONS

A. The direct model

Let Ω = (0,1), Γ = ∂Ω = {0,1}. The problem is to estimate the parameter c(x) in the usual model for wave propagation in acoustic media:

(1)

with initial condition

p(c,x,t = 0)= pt(c,x,t = 0) = 0, x ∈ Ω, (2)

and boundary conditions

p(c, x = 0, t) = 0, (3)
(4)

In seismic exploration p(c, x, t) represents the pressure and c(x) the wave speed of the medium. The function S(x, t) on the right hand side of (1) is a given external source function. In the one dimensional case, x ∈ Ω represents depth. Equations (3) and (4) are respectively a Dirichlet boundary condition at the surface and an absorbing boundary condition which makes the bottom boundary transparent to downward going waves (Stephen et al., 1985).

The set of admissible parameters is denoted by . It consists of the functions c(x) bounded below and above by positive constants c* and c*, respectively; we also assume that c(x) is a linear combination of a fixed set of uniformly continuous functions defined in a neighborhood of the boundary.

We denote by (Ω) the set of functions in H1(Ω) such that (0) = 0. Let (·, ·) denote the usual inner product in L2(Ω) and set ⟨f,g⟩ = f(0)g(0) + f(1)g(1). The weak form of the direct problem (1)-(4) is obtained as usual by multiplying (1) by a test function H1(Ω) and integrating over Ω. After applying integration by parts in the second term in the left-hand side of the resulting equation and using the boundary condition (4) we get the following weak formulation of our differential model: find p(c,x,t) ∈ (Ω) such that

(5)

Let = W1,∞(0, T, L2(Ω)) ∩ L(0, T, H1(Ω)). The following result on the well posedness of the direct problem (1)-(4) is proved in Fernández-Berdaguer et al., (1996) with different boundary conditions; the proof holds for the present case.

Theorem 1Assume that cand that for an integer q ≥ 1, L2(0, T, L2(Ω)). Then the solution p(c,x,t) of (1)-(4) such thatand satisfies the estimates

(6)

where the positive constant C depends only on the total time T and the upper and lower bounds for c(x).

Next we describe the inverse problem.

B. The estimation problem

We assume that the pressure observations pobs(x,t) are recorded for t ∈ (0, T) at points inside Ω. Our objective is to infer from these measurements the actual function c(x), x ∈ Ω, characterizing the medium.

The estimation problem is solved using the output least squares criterion. That is: minimize the following functional

(7)

over the set .

When the set of parameters is considered as a subset of the space L2(Ω) it can be proven that J is continuous as a function of c (see pp. 220-222 in Fernández-Berdaguer et al., 1996). We conclude that the minimization problem has a solution over compact subsets of .

B.1. The sensitivity equations

Methods for solving the above minimization problem require the computation of the derivative of the cost functional with respect to the parameter. Here the parameter is a function, thus we need a functional derivative. We will denote by D(c)δc the Gateaux derivative of p in the direction of the perturbation δc of c. In addition, we assume that the vector space of perturbations δc of c is such that (δc = 0 in a neighborhood of the top boundary.

The function D(c)δc(x,t) can be computed as the solution of the following differential equation:

(8)

with initial conditions

(D(c)δc)(x,t = 0) = (D(c)δc)t(x,t = 0) = 0, x ∈ Ω, (9)

and boundary conditions

D(c)δc(x = 0,t) = 0, (10)
(11)

The functional J(c) has a Gâteaux derivative with respect to the parameter c given by

(12)

Note that the differential system (8)-(11) is very similar to the forward problem (1)- (4), the differences are the source function and the boundary conditions.

The algorithm described in the next section is based on locating the zeros of J'(c).

III. THE ESTIMATION ALGORITHM AND ITS IMPLEMENTATION

Using the above directional derivative we propose a steepest descent method to estimate the parameter.

The iterative algorithm starts with an initial guess c0(x). From Eqn. (12) it follows that the direction of steepest descent of J at the point ck(x) in the k-iteration is given by

(13)

where D*(c) indicates the adjoint operator of D(c). The length μk of the descending step is computed linearizing p around ck (see the appendix). Finally, the parameter ck(x) is updated by ck+1(x) = ck(x) + μkdk. The expression has the representation:

(14)

where the function u is defined as

u(c,x,t) = v(c,x,T-t) (15)

and v is the solution of:

(16)

with initial condition

v(c,x,t = 0)= vt(c,x,t = 0) = 0, x ∈ Ω, (17)

and boundary conditions

v(c, x = 0, t) = 0, (18)
(19)

The function f in (16) is given by

(20)

Note again that this problem is similar to the problem represented by (1)-(4).

Summarizing, the minimization procedure is implemented as follows:

Estimation Algorithm

  • 1. Set k = 0, give an initial guess c0.
  • 2. Compute the direction of steepest descent dk defined by
  • 3. Compute the length μk of the descending step using the formula
  • 4. Update the wave speed as follows:
    ck+1 = ck + μk + dk.
  • 5. Check for convergence. If not, k = k + 1, go to 2.

To carry out step 2 it is necessary to perform the following calculations

  • 2.a) solve the direct problem (1)-(4) for the current value of ck,
  • 2.b) compute the second partial derivative ptt,
  • 2.c) solve (16)-(19) to obtain u and
  • 2.d) compute the integral in (14).

In order to perform the computation, we need to obtain approximations to the solution p(c, ·, ·) of (1)-(4), to the Gateaux derivative D(c)δc(·,·) and to the solution of the adjoint problem (16)-(19). This is done using a finite element procedure to be described below.

Let

and

where P1 denotes the set of polynomials of degree not greater than 1. Let M be a positive integer, Δt = T/M, and un = u(nΔt). Set


The discrete-time Galerkin procedure to compute an approximation to the solution of (1)-(4) is defined as follows:

Find such that

(21)

Similarly, approximations to the solution D(c)δc(c, ·, ·) of (8)-(11) can be computed as follows:

Find such that

(22)

The same procedure was used to approximate the solution v of the adjoint problem (16)-(19).

Find such that

(23)

It has been shown in Santos et al.(1988), that the procedures (21)-(23) are second order correct both in the time and space discretizations and they are stable provided that the Courant-Friedrichs-Levy stability constraint:

ΔtC1(c)h (24)

is satisfied.

IV. NUMERICAL EXAMPLES

The algorithm proposed was implemented to estimate the wave speed of the medium c(x) in the problem (1)-(4). We consider two models of stratified media to show the behavior of the algorithm on different configurations. The first model (model 1) consists of 5 layers with a total size of 1040 m. The second model (model 2) consists of 6 layers with total size of 880 m. The source chosen for all the experiments has the form S(x,t) = S(x - xs)g(t), where g(t) is a given waveform. We employed (Stephen et al., 1985),

(25)

with ξ related to a desired pulse frequency f0 by the formula and K chosen so that max|g(t)| = 1 (see Figure 1). In all the examples, we used f0 = 10 Hz.


Figure 1: Source function.

For both models we considered the case of a single receiver, located close to the surface, i.e., near x = 0. In the examples shown xsource = xreceiver = 17m. Figure 2 displays the simulated trace of observations for model 1 at the receiver position.


Figure 2: Simulated traces at the receiver position for model 1.

Since it can be assumed that the value of the wave speed c at the surface is known, the initial guess c0(x) was taken constant, equal to the value c(0). This is in agreement with the previous assumption that δc = 0 in a neighborhood of the top boundary.

Figure 3 shows the true wave speed profile for model 1, the initial guess and the updated profiles after 150, 300, 450 and 600 iterations. The final estimated profile is quite accurate, eventhough the wave speed corresponding to the bottom layer has been slightly underestimated.


Figure 3: Initial, estimated and true velocity profiles for model 1.

Figure 4 displays the true profile for model 2, the initial guess and the updated profiles after 50, 200 and 350 iterations. In this example, although the final estimated profile is not so accurate, the algorithm was able to find a good approximation of the stratified medium.


Figure 4: Initial, estimated and true velocity profiles for model 2.

Figure 5 shows the behavior of the estimation algorithm for model 2 but in this case the simulated observations were corrupted by zero mean gaussian white noise. To choose the noise variance, three different time windows were considered, according to the amplitude of the trace values. In each window the variance value was chosen in such a way that the ratio of the maximun noise to the maximun signal did not exceed 15 percent. It can be seen that in spite of the presence of noise in the observed data the estimation procedure yields a very good estimate of the true wave speed profile, with an accuracy similar to that in the noise-free case (see Figure 4).


Figure 5: Initial, estimated and true velocity profiles for model 2 in the case of noise corrupted observations.

V. CONCLUSIONS

We presented an iterative method to estimate the wave speed in one dimensional acoustic layered media.

The inverse problem is posed as a functional optimization problem and the gradient of the cost functional was computed by solving the associate adjoint problem. This allows for the formulation of the steepest descent algorithm to solve the estimation problem independently of the discretization scheme used to solve the associated partial differential equations. Here we employed a standard Galerkin procedure using C0-piecewise linear functions as approximating functions. Higher order polynomials or other type of numerical methods could be used to discretize the continuos estimation algorithm. Using higher order polynomials may yield faster and better estimates for the wave speed profiles.

The proposed algorithm requires the solution of three forward problems per iteration, as opposed to full Newton schemes as presented in Fernández-Berdaguer et al., (1996), where we had to solve as many forward problems as parameters. On the other hand, full Newton schemes have much better resolution, so the steepest descent procedure presented here may be used as a first step in a method combining both approaches. This idea is currently being tested.

VI. APPENDIX

The computation of μk proceeds as follows: considering that dk is the direction of steepest descent of J at the point ck(x) find nk that minimizes J(ck + μkdk). We look for a local minimun of J;

Linearizing p around the current value ck and neglecting second order terms, it follows that:

Solving the above equation for μk, the expression used in the algorithm is obtained.

REFERENCES
1. Burk, J.J. and W. Weiss, "Non destructive evaluation of materials", Plenum Press, 265-268 (1979).         [ Links ]
2. Fernández-Berdaguer, E.M., "Parameter Estimation in Acoustic Media using the Adjoint Method", SIAM J. Control Optim. 36, 1315-1330 (1998).         [ Links ]
3. Fernández-Berdaguer, E.M., J.E. Santos and D. Sheen, "An Iterative Procedure for Estimation of Variable Coefficients in a Hyperbolic System", Applied Mathematics and Computation, 76, 213-250 (1996).         [ Links ]
4. Ferry, J.D., "Viscoelastic properties of polymers", John Wiley & Sons, (1970).         [ Links ]
5. Santos, J.E., J. Douglas Jr., M.E. Morley and O.M. Lovera, "Finite element method for a model for full waveforn acoustic logging", IMA J. Numer. Anal. 8, 415-433 (1988).         [ Links ]
6. Stephen, R.A., F. Cardo Casas and C. H. Cheng, "Finite-difference acoustic logs", Geophysics, 50, 1588-1609 (1985).         [ Links ]
7. Stoll, R.D., "Acoustic waves in ocean sediments", Geophysics 8, 715-725 (1977).
        [ Links ]

Received: May 9, 2001.
Accepted for publication: March 18, 2002.
Recommended by Subject Editor R. Sánchez Peña.

Creative Commons License Todo o conteúdo deste periódico, exceto onde está identificado, está licenciado sob uma Licença Creative Commons