Review of Some Exact Methods for the Solution of the One-Dimensional Inverse Problems

Copyright © IFAC Identification and System Parameter Estimation 1982 . Washington D.C .. USA 1982

REVIEW OF SOME EXACT METHODS FOR THE SOLUTION OF THE ONE-DIMENSIONAL INVERSE PROBLEMS F. Santosa Department of Theoretical and Applied Mechanics, Cornell University, Ithaca, New York 14853, USA

Abstract. The inverse problem in geophysics is the determination of the inhomogeneous structure of an elastic medium from experimental observations. In this presentation, we first briefly review the Gelfand-Levitan and the Marchenko methods for the solution of the related quantum inverse problem. In these approaches, as well as in the methods for the geophysical problem we will discuss, the desired information is recovered from the solution of certain integral equations. Next, we consider the construction of a computational scheme for the inversion of geophysical data based on one of the methods. We will further show that the resulting al~orithm is stable under small perturbations to the data. We will finally demonstrate the inversion scheme by inverting noiseless and noise-corrupted synthetic seismograms. Inverse problems; exploration seismology; inverse scattering.

~Ceywords.

individually from get). The quantity that we can recover is the acoustical im-

INTRODUCTION

nedance, in this case, A = iJ:,~) 1/2 , and it will be given as a function of travel time. To be more precise, let us define the travel time for a sif,nal originality at x = to arrive at a point x as

The inverse problem in geophysics is the determination of the inhomogeneous structure of a medium by probing it with waves and observing the reflected signals. The earth is modelled as an elastic half-space, x~ 0, which is smoothly stratified in the direction of depth. Moreover, we consider oly plane wave propagations. In this way, the governing equations can be reduced to the one-dimensional wave equation on the half-line, and is given by p (x)u tt = ( ~ (x)ux)x, t

>

0, x

~

0.

°

n(x) =

<

0, x)

=

°

An impulsive traction is applied at

ux(t, x=o) =

o (t)/ ~ (O)

(5)

ds

A( n) =

(p( n )~( n ))1/2,

where the bars denote that p and ~ are put as functions of travel time n . Making the transformation from the space coordinate x to travel time n in (1), we arrive at

(1)

A( n)u tt = (A( n)un) n ' t

>

0, n

~

0.

(6)

The equations (2) and (4) remain unaltered by the transformation while (3) becomes u (t, n=O) = o (t)/A(0)

(2)

n

(7)

We will assume that A(O) is known, and now the problem is to find A( n) from get). In the subsequent discussion, we will present some methods for solving A( n) exactly from get).

x = 0, (3)

and the response of the medium is measured at the same point, u(t, x=O) = get)

1/2

Denote the acoustical impedance by

Here, u(t,x) is the particle displacement, p, the densit y , and ~ , the shear modulus. For the propagation of longitudinal waves, ~ is correspondingly replaced by ( A +2 ~ ). The system is assumed to be quiescent initially and hence u(t

x

Jo ( p (s)/~(s))

THE GELFAND-LEVITAN AND THE METHODS

(4)

MA~C~E~O

The Gelfand-Levitan and the Marchenko methods were developed for the solution of the inverse problem in quantum scattering theory (see e.g. Chadan and Sabatier, 1977). Their application to the geophysical inverse

The inverse problem is to obtain as much information as possible on the nature of p (x) and ~ (x) given the model in (1) (3), and the data in (4). As we shall see, we cannot reconstruct p (x) and ~(x) 951

F. Santosa

952

problem comes from the fact that (6) can be reduced to the Schroedinger equation. Assume that A is twice- continuously differentiable, and let W Al/2u. Substituting this expression in (6) we get w

- W +q( ll ) W =O 1111 q( ll ) = (A l / 2)"/A l / 2 .

(8)

tt

where Define the time-Fourier transform of w(k, n) as ~ (k, n ). ~ (k' ll ) satisfies 2' , (9) W" + k W = q( ll)W

The idea is to relate

S12(k)

to

q( ll).

For the case where q( ll ) = 0 for 11 < 0, that is, for the situation in which the inhomogeneous half-space 11 ~ 0 is joined with a homogeneous half-space n < 0, the unknown potential q( n) is recovered from a kernel M( n ,s) which satisfies M( 11 , z) + R( ll+Z) +

R(t)

=

f~ ooM( ll ,S)R( s+ z)ds

0,

~ f oo S (k) cp-ikt dk 21' - 00 12

(17)

A

Here, where primes denote differentiations in 11 . Equation (9) is precisely the Schroedin~er equation for the S-waves. The Gelfand-Levitan method considers the following Sturm-Liouville problem, 2 cp " + k CP = q( ll)

The problem is to relate the spectral function ~(k) to q( ll ). In physics, ~ (k) can be calculated using data from scattering experiments. The key in the method is substituting the representation for CP (k, ll ) in (9), which is obtained from the Paley-Wiener theorem, (12)

CP (k, ll ) = COSkll + f llK( ll ,S)cosksds o

21

=

OO

where

T( ll ,Z) =

A

f oCOSkll coskz(v(k)

The unknown potential recovered in

q( ll )

-

0, (13)

2

~)dk.

is finally

1 d 2" q( ll ) = dll K( ll , Z= ll)

Therefore, if one can find v(k) experiments, one can reconstruct using equation (13) - (14).

(14) from the q (11 )

=

OO

-2k/nf g(t)sinktdt o

M( ll , Z=ll )

THE METHODS OF SYMES, GOPINATH AND S0NDHI, AND CARROLL AND SANTOSA Symes (1979) constructed a nonlinear version of the Gelfand-Levitan method specifically for the geophysical problem. The construction, like that of Burrid~e, is in the time domain, uses Riemann function representations. He showed that if one solves L(y,z) from (see also Burridge, 1980)

t

[g( ll+Z)+g( !ll -Z !)] o

(15)

The Harchenko method consider (10) for _ 00 < 11 < 00 with the following asymptotic boundary conditions,

The first term corresponds to the incoming probe, and the second, the reflected waves.

(18)

We remark that the discussions here are much simplified and the reader is refered to Chadan and Sabatier (1977) for a detailed exposition of the full quantum scattering problem. The application of the Marchenko method for the geophysical problem can be found in Ware and Aki (1968) where they discussed the calculation of R(t) from get) in (4). The method is also rederived in a direct way in a recent review by Newton (1981). Both the Gelfand-Levitan and the Marchenko methods are presented in the time domain by Burridge (1980) and precise connections between the two methods are made. Time domain formulations have the advantaRe that they are more direct than the Fourier transform domain formulations.

L( ll ,Z)+fllL( ll ,S)L(z,s)ds

This method was first applied to the geophysical inverse problem by Alekseev (1962), who showed that the spectral function, v(k), can be calculated directly from the data get) in (4). The relation is given by V(k)

d dn

The main tool used is a Titchmarsh theorem for Fourier integrals. Thus, knowi~g S12(k), which can be calculated from exper~metaI data, q( ll ) can be found.

Here, K( ll ,S) is an unknown kernel which can be found by solving the Gelfand-Levitan equation T( ll ,Z) + K( ll ,Z) + f:K(ll,S)T(S,z)ds

q( ll ) =

where

get)

(19)

is directly from the data in

(4), one can recover the unknown potential

in

1 d 2"q( ll ) = dll [L( ll , z=ri)]

(21))

A generalization of this method can be found in (Symes, 1981) where the question of stability of the inversion to perturbations on the data is d~scussed. The method of Gopinath and Sondhi (1971) (see also Burridge, 1980) was derived for the transmission line inverse problem. The method cleverly convolves a known solution of (6) with the solution U(t, ll ) of the problem to obtain the equation

953

Methods for the Solut i on of One- d i mension a l Inve r se Problems 1 =

In - ng r' ( l t-T I )G( T , n )dT +

where

g r (t) = 1 + get)

for

t > O.

unknown i mpedance with A(O) = 1, ly rela t ed t o G(t, n ) in

The

is direct -

(22)

A(y) = [G(t = n , n )]2

The method works for A' ( n ) continuous as opposed t o the methods in the previous sec is assumed to be cont i nt ion where A" ( n ) uous . The au t ho r, i n collaboration with Carroll (Ca rr oll and Santosa , 1981) , generalized the Gelfand-Levi t an me t hod for the case where A' (n ) is cont i nuous . He r e , we pose a simila r Sturm- Li ouville problem ,

(A(n)~ ' ) ' /A(n ) ~(k,O)

= 1

,

+

k2~

~ ' (k , O)

0 (23)

0

= 6 (x- y)/A( n )

(24)

The idea i s to relate v (k) to A( n ). By applying the Paley - Wiener theorem , we were able to show that the unknown impedance A(n) can be recovered from a kernal H( n ,s) wh i ch sa ti sf i es T( n , z) + H( n , z) - f nH( n ,s)T (s , z)ds o s where

T( n , z) =

f:sink ncoskz/k( ~ (k)

v(k) = - 2k/ lT f g(t)sinktdt o

(29)

Next , define the TH = -

T( n , z) = Ts(s , z) = -

21

T( n , z)

and

-

j

Ti(j) =

L ~=l

[ 6 . ~ -T . ~ ]h.(j) ~

~

(31)

~

Thus , if we t hink of the data p,(t) as given by a time series at time increments of tl , we can construct an approximate solution of H( n ,x) d i scretely as hi(j) using (31) . The unknown impedance can be calculated from ,

(1 _ h . (j)) - 2 J

J

J

(32)

Therefore , if we are given a time series {get = j tl ) } for j = 1, 2 , 3 , ... , 25 , we can solve the j x j matrix equa t ion (31) for j = 1 , 2,3 , ... , J for hi(j) , from which we can find

a . " A( n =H), J

for

1 ~ j ~ J . Moreover , we can show that this approximat i on is conve r ~ent , t hat is as tl ~ 0 and j tl = y , a . ~ A( n) . The J

Ts(S ' z)

g~( l s - z l )]

(30)

Then , equation (25) is replaced by the approximate version ,

convergence is pointwise and the er r or committed by the discretization is of the order tl . Next , we study the behavior of this approximation to perturbation on the data get) . Let the noise corrupted data g*(t) = get) + h(t) , where het) and h'(t) are small perturbations . The re covered impedance from g* will be denoted

[gr( n+ z ) - gr( n- z )]

t [g~(s+z)

~ [g r' ( ~tl+i tl ) -g r' ( I~tl - i tll )]

a.

In a subsequent paper (Santosa (1982)) , we showed that this formulation can be given time - doma i n meaning . If we set, as before , g r (t) = 1 + get), and A(O) = 1 , we can find the kernels from gr(t) ,

matrix

j x j

2/ lT )dk .

(27)

PROBLEMS

hi (j) " H(y = H , x = ill), i < j

A( n = j tl )

(26)

I~VERSE

[gr(jMitl ) - gr(j tl - i tl )] , i < j

a.

Mo r eove r, using the same arguements as Alekseev , we related the spectral function v(k) t o the data get) , OO

t

Ti (j)

(25)

o

The unknown i mpedance is reconstructed using

-2 A( n ) /A(O) = (1 - H( n , z = n )) .

O~

The following discussion is a summary of the results present in Santosa (1932) . We are basically interested in solving equation (25) numerically. We propose to do so by r eplacing the integral i n the equation with a finite same approximation . Define the j - dimensional vectors

We also consider the completeness relat i on

Joo~ (k , x) ~ (k , n ) ~ (k)dk o

ASPECTS

NUMERICA~

(21)

G(t , n )

(28)

Thus , the inversion is made more direct . Stability results for this method can be found in Carroll and Santosa (1982). Although all t hese methods are exact, the actual solu t ion of the inverse problem itself in volves solv i ng integrals equations whose general solutions are not known . Therefore , an important question is can one construct an accurate numerical scheme based on these me t hods for the inversions of the data . Below we d i scuss a computational method based on the t echnique of Carroll and Santosa .

by and

j

a~ .

Let

J h

j

o

L

L

= max 1 ( 6 i ~ -Ti £ ) 0 £ i =l Ihi(j) 1 (cf . (31)) , then s

I .

i

2j

+ Ll h(p tl ) l}

(33)

p =O Thus , if we allow the data to contain a small amount of perturbations het) and h ' (t), the reconstructed profile a ~ J

ISPE - 2- E*

-1

F. Santosa

954 differs from

a

j

by a small amount also.

REFERENCES

This is the stability. In the figures below, we show some results of a numerical experiment on the developed inversion method. In all the inversions, we take 6 = 0.01 and we took J = 50. In Fig. 1, we show the original profile. In Fig. 2, the reflection data g (t) = 1 + g(t), calculated using

Alekseev, A. (1962). Some inverse problems of the theory of wave propagation, Pt. I and 11, Izv. Akad. Nauk. SSSR, Ser. Geofiz., 11, 1515-1531. Burridge, R. (1980). ThelGelfand-Levitan, The Marchenko, and the Gopinath-Sondhi integral equations of inverse scattering theory, regarded in the context of inverse impulse-response problems, Wave Motion, 2, 305-323. Carroll, R. and Santosa, F. (1981). Inverse scattering techniques for a onedimensional inverse problem in geophysics, Math. Meth. Appl. Sci., 1, l45-l7l. Carroll, R. and Santosa, F. (1982). Stability for the one-dimensional inverse problem via the Gelfand-Levitan equation, to appear in Applicable Analysis. Chadan, K. and Sabatier, P. (1977). Inverse Problems in Quantum Scattering Theory, Springer-Verlag, New York. Gopinath, B. and Sondhi, M. (1971). Inversion of the telegraph equation and the synthesis of nonuniform lines, Proc. IEEE, 59, 3R3-392. Newton, R. (1981):- Inversion of reflection data for layered media: review of exact methods, Geophys. J.R. astr. Soc., ~, 191-215. Santosa, F. (1982). Numerical scheme for the inversion of acoustical impedance profile based on the Ge1fand-Levitan method, to appear in Geophys. J.R. astr. Soc. Symes, W. (1979). Inverse boundary value problems and a theorem of Gelfand and Levitan, J. Math. Anal. App1., 11, 379-402. Symes, W. (1981). Stable solutions of the inverse reflection problem for a smoothly stratified elastic medium, SIAM J. Math. Anal., 12, 421-453. Ware, J. and Aki, K. (1969~ Continuous and discrete inverse scattering problems in a stratified elastic medium, J. Acous. Soc. Am., ~, 911-921.

r

the method of characteristics is shown. Fig. 3 shows the reconstruction from the synthetic data. In Fig. 4, we added to the signal gr(t), noise of the amount 100 n

=

100

L Ih'(p6) 11 L Ig'(p 6)1=0.30.

In Fig. o p=l p=l r 4, we display the result of inverting the noise-corrupted synthetic data. In conclusion, we remark that it is impossible to include in this presentation the other methods for the one-dimensional inverse problem. The methods discussed here were chosen because they display the common feature, which is that the desired information is recovered in the solution of certain integral equations. It should be noted that many other results, theoretical and applicable, exist in the inverse literature. ACKNOWLEDGEMENT This work was supported by the Material Science Cent er at Cornell University.

3.

2.

A (y) I.~---

O. L-_ _. L -_ _ 0.1 0.2

~

_ _..I....-_ _-'--_----J

0.3

0.4

0.5

Y Fig. 1.

The original acoustical impedance function.

955

Methods fo r the So lution of One-dimensional Inverse Problems

0.2

-0.2 Fig . 2 .

The theo r e tical ref l ection data for the impedance gi ven in

3. 2.

A (y) I . ~----

O . ~---L---~--~---~--~

0.1

0.2

0.3

0.4

0.5

Y Fi g . 3 . The r econs truc ted function from the theoretical data.

0 .2

-0 .2 Fi g . 4 . The no i s e-corrupted reflection data.

~ ig .

1.

956

F. Santosa

3.

2.

A (y) 1. 1-----OL-----~-------L------~------L---

01

02

0.3

__~

0.4

05

Y Fig . 5. The reconstructed function from noisy data .