Tutorial: unified 1D inversion of the acoustic reflection response.

Acoustic inversion in one-dimension gives impedance as a function of travel time. Inverting the reflection response is a linear problem. Recursive methods, from top to bottom or vice versa, are known and use a fundamental wave field that is computed from the reflection response. An integral over the solution to the Marchenko equation, on the other hand, retrieves the impedance at any vertical travel time instant. It is a non-recursive method, but requires the zero-frequency value of the reflection response. These methods use the same fundamental wave field in different ways. Combining the two methods leads to a non-recursive scheme that works with finite-frequency bandwidth. This can be used for target-oriented inversion. When a reflection response is available along a line over a horizontally layered medium, the thickness and wave velocity of any layer can be obtained together with the velocity of an adjacent layer and the density ratio of the two layers. Statistical analysis over 1000 noise realizations shows that the forward recursive method and the Marchenko-type method perform well on computed noisy data.

INTRODUCTION

Backus (1959) showed that the reflections in a marine seismic trace can be understood as the subsurface reflection response filtered by a water layer reverberation operator that is independent of the source and receiver depth. He also showed how a simple three‐term deconvolution filter can be constructed and used to remove these water‐layer reverberations. Kunetz (1964) showed that both free‐surface and internal multiples from the shallow part of the reflection response can be filtered from the data. Application of this filter removes overlap from multiples and the first deeper primary reflection which can then be analysed. He showed that the filter consists of two parts of a fundamental wave field and gave their mutual recurrence relation. Combining the filter and the recurrence relations leads to a forward recursive inversion scheme. It successively filters from top to bottom the local reflection coefficients from the data. The same idea, but using an energy equation, was used by Robinson (1967) and Robinson and Treitel (1977) to construct a similar forward recursive scheme. They also showed that the filters at each step can be used to compute the up‐ and down‐going wave fields at depth from the surface reflection response. Robinson and Treitel (1978) showed how the energy relation can be used to construct a scheme that starts in the bottom of a model. It leads to a backward recursive scheme to obtain the local reflection coefficients. All the work was done with the Goupillaud's (1961) model of a discrete‐layered medium. In these three schemes, one equation is used that involves the data. Kunetz (1964) used the convolutional model, whereas Robinson and Treitel (1977, 1978) used an energy relation. The other equations are the recursive relations of the two parts of the fundamental wave field to make a forward or backward recursion step. The reflection coefficient is retrieved as the amplitude of an event at the proper two‐way travel time, which is then also found. From the refection coefficients, the layer impedances are found and this completes the inversion.

Many people worked on the exact inverse solution of the Schrödinger equation to resolve the scattering potential from the scattered field (Agranovich and Marchenko 1963; Lamb 1980). This solution was taken by Ware and Aki (1969) who showed that if depth was converted to vertical travel time and the fields were flux‐normalized, the one‐dimensional acoustic wave equation is equal to the Schrödinger equation. The inverse solution consists of a time integral from which the fundamental reflection solution can be computed. The acoustic impedance at any fixed vertical travel time is obtained by integrating the fundamental reflection solution over all times. It requires the zero‐frequency component to be recorded in the seismic record (Berryman and Greene 1980). Rose (2002) showed how we can understand the solution to the Marchenko equation as creating a focused wave field at depth. This led Broggini, Snieder and Wapenaar (2012) to the understanding that this focus leads to a response at the acquisition surface as coming from a virtual source at the focal point. By reciprocity, this is the same as the understanding of finding the wave field at depth from a source at the acquisition surface by Robinson and Treitel (1977). Slob et al. (2014) showed that the fundamental reflection solution of the Marchenko equation is a particular combination of the up‐ and down‐going parts of the fundamental wave field. They showed how the Marchenko equation can be seen as a sum of two other equations. One that represents a seismic experiment and one that represents a time‐reversed seismic experiment. They derived these two equations by applying acoustic reciprocity of the time‐convolution and time‐correlation types. They showed that the local reflection coefficient at any vertical travel time instant is obtained directly from the up‐going part of the fundamental wave field. Wapenaar et al. (2013) derived the single‐sided Marchenko equation for a three‐dimensional heterogeneous medium, and connected the concepts of Marchenko inverse scattering and Green's function retrieval. This triggered new research on redatuming, imaging and inversion with minimal model information and multiple elimination without model information. Review and tutorial papers have been published that help getting into the subject of Marchenko redatuming and imaging (Wapenaar et al. 2017; Nowack and Kiraz 2018; Lomas and Curtis 2019).

Bardan and Robinson (2018) connect the result of the forward recursive scheme of Kunetz (1964) to the solution of the discrete version of the Marchenko equation from Berryman and Greene (1980). They conclude that the discrete solution of Kunetz (1964) is the same as that of the discrete Marchenko equation. Here, we show that the two bodies of thought described above are connected through the fundamental wave field. The fundamental wave field as derived in discrete forms of the recursive schemes has up‐ and down‐going parts and together make up the fundamental reflection solution of the Marchenko equation. We show how the equations that are used for the forward and backward recursive schemes can be combined at any depth level to yield the two parts of the Marchenko equation. We then show how the combination of the two equations removes the need for recursive solutions. Only in this way the inversion can be done at any chosen vertical travel time in a target‐oriented way for finite‐frequency‐bandwidth data. We show how the reflection response along a line over a horizontally layered medium leads to the possibility of obtaining layer thickness and velocity of a target layer, together with the velocity of the two adjacent layers and the density ratios of these layers.

First, we give the expressions for a two‐sided experiment to build the necessary expressions for the impulse reflection and transmission responses of a layered medium. We establish the convolutional model and the time‐reversed experiment that lead to the energy relations. The subsurface impulse reflection response is represented by the fundamental wave field and we give the recurrence relations between the up‐ and down‐going parts. We derive Green's function representations in terms of the reflection response and the fundamental wave field. From these representations, we derive the Marchenko equation. Second, we derive four inversion methods and show that they are based on the same set of equations. We show how the inversions are carried out and which information is obtained. Third, we show how only the inversion with the Marchenko‐type method can be extended to target‐oriented inversion of finite‐frequency‐bandwidth data. Finally, we give a numerical example to illustrate the performance of the methods on noisy data and evaluate the results.

THE CONVOLUTIONAL MODEL OF ACOUSTIC EXPERIMENTS

In this section, we describe the discrete‐layered model and the associated wave fields that can be measured on both sides of the model given a source in the top or in the bottom. Without loss of generality, we model the acquisition surface as a transparent boundary in the upper half space and in the lower half space. Inclusion of a free surface would not add an unknown reflector as explained by Kunetz (1964). We use the physics of acoustic wave propagation in a linear model, which says that the response of a medium is the convolution of the time signature of the source and the earth impulse response. Later, we use this understanding to unravel the subsurface reflection information from the data. We do not use z‐transforms and the associated Goupillaud's (1961) model, but use expressions in the time domain, with time denoted by t.

Figure 1 shows the model configuration with possible up‐ and down‐going wave fields at the top and bottom of the layered medium. The up‐ and down‐going parts of the acoustic pressure, and , respectively, are defined at and and just below . The discrete‐layered medium has reflecting boundaries at depth levels , . Each layer is characterized by constant velocity, , and density, . The depth axis points downwards, hence . The thickness of each layer is , for . For an interface at , the reflection coefficient for a plane pressure wave coming from above is indicated by , with the impedance given by . The vertical travel time in each layer is denoted for . We define the travel time and distance to z 0 as and . Cumulative travel time from to is denoted . The transmission coefficient is denoted and the cumulative transmission coefficient is , where the plus‐sign applies to transmission from to and the minus‐sign from to .

Figure 1

The model for a physical experiment in a one‐dimensional piecewise‐homogeneous medium with reflecting boundaries with transparent boundary at , has down‐ and up‐going fields, and at and and at , respectively.

Reflection and transmission experiments

In the model shown in Fig. 1, the up‐ and down‐going fields are related to each other through the impulse reflection and transmission responses that is expressed as where * denotes temporal convolution, and denote the impulse reflection and transmission responses, at and , respectively, in case the down‐going wave field at and the up‐going field at , whereas and denote the impulse reflection and transmission responses, at and , respectively, when the up‐going wave field just below and the down‐going field at . They are impulse responses, or Green's functions. The column vector in the left‐hand side of equation (1) is the down‐going wave field that could be measured at and the up‐going wave field, , that could be measured at . The column vector in the right‐hand side of the equation contains the initiating down‐going wave field that could be emitted at and up‐going wave field that could be emitted at . These wave fields are shown in Fig. 1 and are connected to each other through the impulse reflection and transmission response matrix as expressed in equation (1). This is the convolutional model of a physical acoustic experiment. In this paper, all equations that represent an experiment have this structure. We illustrate equation (1) with a numerical example. We take a model with five reflectors (), which we use throughout the paper. The medium parameters and the local reflection coefficients are given in Table 1. We emit a Ricker wavelet at 75 m above the top reflecting boundary. It is shown as the down‐going field, , in Fig. 2(a) and . The reflection response, , is shown in Fig. 2(b). The total wave field everywhere in the model is shown in Fig. 2(c) as a function of vertical travel time, ζ, and recording time, t. The reflectors are indicated by the horizontal black lines and the vertical black line indicates the time it takes the initial down‐going wave to propagate to z 2. We use this level later to introduce a truncated medium which facilitates our analysis. The slanted black line marks . Zero vertical travel time indicates the acquisition surface in the top. The down‐going field is shown in Fig. 2(d). The sum of the two traces in Fig. 2(a,b) is the same as the field at zero vertical travel time in Fig. 2(c). The wave field in Fig. 2(d) is the same as the field shown at the latest vertical travel time in Fig. 2(c). Some of these events are indicated in Fig. 2(c) with arrows to show in which direction they propagate.

Table 1

Values for velocity, density, layer thickness and local reflection coefficient in the layered model

Velocity (m/s)	Density (kg/m3)	Thickness (m)	reflection Coefficient (−)
1500	1000	∞	0.6364
3000	2250	117	−0.6364
1500	1000	99	0.4545
2000	2000	85	−0.2075
1750	1500	111	0.3538
2750	2000	∞	–

Figure 2

The down‐going (a) and up‐going (b) wave field at , the wave field inside a medium with three reflectors as a function of vertical travel time and recording time (c) and the down‐going wave field below (d).

Time‐reversed reflection and transmission experiments

In a lossless medium, the wave equation is symmetric in time, which we exploit in a time‐reversed experiment. A down‐going field at the acquisition surface in a physical experiment is time reversed and becomes an up‐going field in a time‐reversed experiment and an up‐going field becomes a down‐going field upon time reversal. This configuration is shown in Fig. 3. By carrying out a time‐reversed experiment using the responses measured in the physical experiment, we recover our original source time function as the response (Fink 1992). Interchanging the up‐ and down‐going fields in the column vectors of equation (1) and reversing their time dependency gives the expression for a time‐reversed experiment, which is then given by Equation (2) is the mathematical expression of the experiment depicted in Fig. 3. For later convenience, we take the time reverse of equation (2), reorder the fields to make an equation that resembles equation (1) and find From here onward, equation (3) is how we define a time‐reversed experiment. It shows that when the field responses, and , at and , respectively, generated by and in the physical experiment, are used as the emitted wave fields in a time‐reversed experiment, the response is the original emitted field and . This is the convolutional model of the time‐reversed experiment corresponding to that of the physical experiment. We interpret Fig. 2 using equation (3). Equation (3) shows in the column vector in the right‐hand side that the wave fields shown in Fig. 2(b,d) are the input wave fields, from above and from below, respectively. According to equation (3), the medium responds in reverse time. This means that the wave field propagates inside the layered medium in the direction of decreasing recording time. This can be seen in Fig. 2(c) when all arrows are rotated 180°, such that what was upward pointing is now downward pointing, and vice versa, and they all point towards decreasing recording time. Consequently, the only result is the Ricker wavelet shown in Fig. 2(a).

Figure 3

The time‐reversed experiment that corresponds to the physical experiment of Fig. 1 has down‐ and up‐going fields, and at , and and at , respectively.

We use the fact that the input fields , where superscript t denotes matrix transposition, in equation (3), are the responses in the left‐hand side of equation (1). Combining these two equations results in a statement that the input is the same as the output. Since this must be true for any input, we decide to use once and and once and , which results in the matrix expression which expresses the unitary property of the reflection and transmission impulse response matrix. It is the well‐known expression of the conservation of acoustic energy in the system. The two equations involving the first column of the second matrix in the right‐hand side of equation (4) are written as Equations (5) and (6) describe the time‐reverse experiment of equation (3), with the incident fields given by and . Equation (5) was given in Robinson and Treitel (1977) and forms the basis of their methods, which we discuss later. We use these relations when we want to interrogate the medium in its interior using only the reflection response .

Relations between the reflection and transmission responses and the fundamental wave field

Impulse reflection and transmission responses are not independent from each other. Compact frequency domain expressions were given in terms of the up‐ and down‐going parts of a fundamental wave field in optics by Abeles (1946) and in acoustics by Goupillaud (1961). Kunetz and d'Erceville (1962) introduced the term fundamental polynomials and Kunetz (1964) gave coupled recurrence relations for the fundamental polynomials. Here we give a treatment in our notation to link their equations to the recent work on Marchenko‐type equations to obtain reflectivity from the reflection response (Slob et al. 2014). If a layered system has reflecting boundaries, with the bottom reflector at and , we call it a truncated medium. Its impulse reflection response from a source in the top is denoted and its impulse transmission response just below is denoted . They are given in terms of the two parts of the fundamental wave field related to the truncated medium, and , as where denotes the cumulative transmission coefficient, as defined in the paragraph describing Fig. 1, and is explained below. Equations (7) and (8) can be interpreted as the convolutional model of a seismic reflection and transmission experiment. According to equation (1), when is the down‐going wave field incident on the truncated medium and , is the reflection response and the transmission response. The latter is the physical direct arrival in the impulse transmission response, which we denote as as shown in equation (8). For this reason, we interpret as the internal multiple eliminator, or the anti‐reverberation filter, for the transmission response of the truncated medium. The expression of equation (7) was given in Kunetz (1964).

The special character of the two wave fields is captured in their coupled recurrence relations given by (Kunetz 1964) where are the up‐ and down‐going parts of the fundamental wave field for the truncated medium that has its bottom reflector at . Both equations are used in the method of Kunetz (1964). The forward recursion is initialized by The importance of the initial unit amplitude impulse lies in the fact that is a causal minimum delay function for arbitrary n, provided all . For this reason, the inverse exists, denoted , which is a causal minimum delay function as well, with . Because is understood as the multiple eliminator, is understood as the multiple generator, such that .

It is of interest to make three observations from equations (9) and (10). First, have the same finite number of events and their number is equal to 2. This can be understood from the fact that for both filters have one event, see equations (11) and (12), and every time we add a reflector, we double the number of events. Second, all primary reflections that occur in the impulse reflection response have the local reflection coefficient as amplitude in and arrive at the expected two‐way travel time. This can be seen from equation (10) as follows. Let us take , is then constructed from , which contains the first reflection event at its physical arrival time, and which will never change when we increase n. The second term in the right‐hand side has a unit amplitude impulse in , which results in the contribution and contains the local reflection coefficient as an event at the physical arrival time. Then we take and has the reflection events from the first two reflectors as present in and is obtained from the unit amplitude impulse in . The second event in will create a non‐physical event in . This will happen for any value of n where all coefficients , , are contained in and comes from the leading impulse in . Third, in no events are present outside the time window between and . All events in are primary reflections and the ones that can be traced back to the initial unit amplitude impulse at are the physical primary reflections. Emitting into the truncated medium leads to the minimum number of events possible in a reflection and transmission experiment for that truncated medium. The reflection response is , which has the same number of events as the input signal. In this sense, the functions are called the up‐ and down‐going parts of the fundamental wave field, . The transmission response has a single event, coming from the initial down‐going impulse in .

To illustrate this, we use a truncated part of the model as in Fig. 2, with . We show the wave fields in the truncated medium in Fig. 4. The down‐going field, , is shown in Fig. 4(a). The corresponding reflection response, , is shown in Fig. 4(b). How the waves in Fig. 4(a) propagate into the medium and lead to waves in the response shown in Fig. 4(b) is shown in Fig. 4(c). Similar to Fig. 2(c), here the plot shows the propagation throughout the medium as a function of vertical travel time, ζ, and recording time, t. The reflectors are indicated by the horizontal black lines and the vertical black line indicates the time it takes the initial down‐going wave to propagate to z 2. The slanted black line marks . Zero vertical travel time indicates the acquisition surface. In the figure, the four waves that are emitted into the medium are indicated by the downward pointing arrows. These are the same as the ones shown in Fig. 4(a). The four waves that arrive at the acquisition surface are indicated by the upward pointing arrows. These are the same as the ones shown in Fig. 4(b). Note that the acquisition surface is an acoustic transparent surface where no reflections occur. Just below the black line that marks the bottom reflector we see the single wave propagating down into the lower half space as expressed in the left‐hand side of equation (8) with and indicated by the downward pointing arrow in the bottom of the figure. This is the physical first arrival in an impulse transmission experiment and is the only down‐going wave as shown in Fig. 4(d). This confirms that is the anti‐reverberation filter for the transmission response of the truncated medium. It has the same meaning for the reflection response, albeit that in that case non‐physical primary reflections end up in as can be seen in Fig. 4(b) where the third event is a non‐physical primary reflection. This event arises from the third event in which is emitted to prevent a multiple to be generated at the second reflector. Only then a single wave travels down below the third reflector.

Figure 4

The down‐going anti‐reverberation filter (a) and up‐going truncated medium response (b) at , the wave field inside the truncated medium with three reflectors as a function of vertical travel time and recording time (c) and the down‐going transmission response of the truncated medium in the bottom below (d).

To understand the corresponding time‐reversed experiment, we replace the reflection and transmission responses in equations (5) and (6) by those of the truncated medium. We then convolve both sides by and use equations (7) and (8) to find where and are defined as and but now for the truncated medium. Let us interpret Fig. 4 as a time‐reversed experiment. By rotating all arrows 180° in Fig. 4(c), all propagation takes place in the reversed‐time direction. The waves shown in Fig. 4(b) are now sent into the medium at zero vertical travel time and propagate in the direction of decreasing recording time. The reflection response at is given by the first term in the right‐hand side of equation (13). The corresponding transmission response at is given by the first term in the right‐hand side of equation (14). The wave in the bottom of the model is now an up‐going wave incident on the three reflectors and it propagates through the model in the decreasing time direction. Its initial strength is and it starts at . The transmission response at is given by the second term in the right‐hand side of equation (13). The corresponding reflection response from below at is given by the second term in the right‐hand side of equation (14). The sum of the two terms in the right‐hand side of equation (13) forms the total response at the acquisition surface, , given by the left‐hand side of equation (13). These are the waves shown in Fig. 4(a). The sum of the two terms in the right‐hand side of equation (14) forms the total response just below . This total response is zero as can be seen in Fig. 4(c), where no waves exist below the bottom reflector as a function of decreasing recording time and increasing vertical travel time. Hence, the waves in Fig. 4(c) can be understood as propagating in increasing or decreasing recording time direction without changing anything in their amplitude and time behaviour.

If we convolve both sides of equation (13) with , use equation (7) and reorder the terms, we obtain the useful relation (Kunetz 1964) To obtain the right‐hand side of equation (15), we have used an equation similar to equation (8) but for transmission from bottom to top, with replaced by and replaced by . We use equation (15) when we discuss inversion methods in the next section.

Single‐sided time‐reversed experiment

The above time‐reversed experiments were carried out as two‐sided experiments. We want to be able to find the medium parameters using only the single‐sided reflection response. For this reason, we assume that we know only the reflection response at and write equation (13) as where is the Green's function describing the pressure at generated by an up‐going impulse just below as indicated by the minus‐sign in superscript. Equation (16) states that when the up‐going fundamental wave field is sent into the medium for a time‐reversed experiment, the reflection response is the down‐going fundamental wave field minus the time‐reversed Green's function convolved with the direct arrival in the transmission response. Equation (16) is of course the same equation as equation (13), but with different meaning and interpretation. We illustrate equation (16) as a single‐sided time‐reversed experiment in Fig. 5. In this situation, the field that is incident on the medium only comes from above and is , shown in Fig. 5(a). The time‐reversed reflection response of the medium is shown in Fig. 5(b) where the function minus the scaled and delayed Green's function are depicted. The part of the filter that does not overlap with minus the scaled and delayed Green's function is shown in black, whereas the part of minus the scaled and delayed Green's function that does not overlap with the filter is shown in green. It can be seen that the two functions have one overlapping event, which is depicted as a dashed black–green line, because it is the sum of the two terms. The difference between Figs. 4(c), interpreted as a time‐reversed experiment, and 5(c) is the up‐going incident wave from the bottom that is present in Fig. 4(c) and absent in 5(c). This absence can be understood as emitting that impulse together with the same impulse with opposite sign. Consequently, no up‐going wave is visible in the bottom of Fig. 5(c) and as a result there are propagating waves that arise from the extra negative amplitude impulse that leads to minus the scaled and delayed time‐reversed Green's function at the acquisition surface. With this understanding, we are ready to use the fundamental wave field of the truncated medium in physical and time‐reversed experiments using the impulse reflection response of the actual medium.

Figure 5

The time‐reversed filter responses in the truncated medium; (a) The incident wave field, ; (b) the reflected wave field composed of (black line), (green line) and the sum where they overlap (dashed black‐green) at ; (c) the wave field propagating in reverse time inside a medium with three reflectors as a function of vertical travel time and recording time.

Filtering the impulse reflection response with fundamental wave fields of a truncated medium

In this section, we derive expressions for physical and time‐reversed experiments when the incident field is coming from above only. We begin with writing the impulse reflection response as the sum of the impulse reflection response of the truncated medium and a scaled and delayed Green's function. This was derived by Goupillaud (1961) and we write his expression as where is a Green's function of the actual layered medium. It describes the acoustic pressure at , generated by a down‐going unit‐amplitude impulse just below as indicated by the plus‐sign in superscript. The convolution of and describes the part of the impulse reflection response that is not described by . We convolve all terms in equation (17) with , use equations (7) and (8) and write it as a single‐sided experiment given by Equation (18) states that when the incident field is the down‐going part of the fundamental wave field of the truncated medium, the reflection response is the up‐going part of the fundamental wave field of the truncated medium plus a scaled and delayed Green's function of the actual medium. There is no overlap in time between these two parts of the reflection response. To show that, we illustrate equation (18) with a numerical example with and . Figure 6(a) shows as incident field which is the same as in Fig. 4(a). Figure 6(b) shows the corresponding reflection response, the black line is and is the same as in Fig. 4(b), whereas the green line is the scaled and delayed Green's function. Figure 6(c) shows the wave field in the actual medium. The dotted line marks the vertical travel time where the actual medium is truncated. The black solid lines indicate the reflecting boundaries and the dashed lines indicate the time window of the fundamental wave fields. This means that the part inside the time window marked by the dashed lines is the same as Fig. 4(c). Consequently, all waves outside the time window marked by the dashed lines are part of the delayed and scaled Green's function. As long as we truncate the actual medium between the third and fourth reflector, the wave field does not change anywhere.

Figure 6

The response of the actual medium to ; (a) incident field , (b) reflection response composed of (black line) and (green line) at and (c) the wave field propagating in a medium with five reflectors as a function of vertical travel time and recording time.

Equation (18) was used by Kunetz to invert a trace with a forward recursion method using equations (9) and (10). Equation (18) is also one of the two equations used for the Marchenko method that makes non‐recursive use possible. Before we can demonstrate that, we must obtain an equation for the single‐sided time‐reversed experiment. We need the expression for the impulse transmission response, which is written as where is the pressure Green's function of the actual layered medium. It describes the impulse response just below generated by the unit amplitude down‐going impulse starting just below . When the actual impulse is not just below but at , we need to convolve the Green's function with the transmission response of the truncated medium. This follows directly from the convolutional model of an acoustic experiment. We write equation (19) with the aid of equation (8) as an experiment, given by The left‐hand side of equation (20) is the wave field just below the black line marking the bottom reflector in Fig. 6(c).

In the two‐sided time‐reversed experiment, the responses as given in the left‐hand sides of equations (18) and (20) are emitted into the medium from top and bottom, respectively, and propagation takes place in reversed time. This leads to as the response at the acquisition surface, given by Equation (21) is found by convolving all terms in equation (5) with and then using equations (18) and (20) in the resulting right‐hand side. We understand equation (21) from Fig. 6(c) in the way we used Fig. 4(c) to understand equation (13). In a time‐reversed experiment, there are many waves sent in from below and from above by the scaled and delayed Green's functions, together with the waves in that are sent in from above. They result in at the acquisition surface. Sending in the scaled and delayed Green's functions from the top and bottom leads to a single up‐going impulse with amplitude that arrives just below z 2 at . It will then continue to propagate in reverse time. At the acquisition surface, this would lead to the time‐reversed scaled and delayed Green's function . Just below z 4, it would be the time‐reversed scaled and delayed Green's function . However, because we send back into the medium as well, these scaled and delayed time‐reverse Green's functions are cancelled. Hence, no waves exist for recording times smaller than the vertical travel time as can be seen in Fig. 6(c). The sum of the two terms involving the Green's functions in the right‐hand side of equation (21) is equal to the scaled and delayed time‐reversed impulse response, or the convolution of the Green's function with the direct arrival of the transmission response in the truncated medium, hence , at the acquisition surface. This is true for any n and m with and being an arbitrary example.

We carry out a single‐sided time‐reversed experiment and send in only . In that case, this time‐reversed Green's function will be part of the response with a minus sign. For this reason, we rewrite equation (21) as with Equation (22) is interpreted in the same way as equation (16). The only difference is that now the Green's function belongs to the actual medium and not to the truncated medium. We illustrate equation (22) with a numerical example where the waves propagate in the same model as used for Fig. 6, but now in reversed time. We send in only and the reflection response is as given in the left‐hand side of equation (22). Equation (22) is of course the same as equation (21) but interpreted differently. The wave field corresponding to the interpretation of equation (22) is shown in Fig. 7. Figure 7(a) shows the emitted wave field and Fig. 7(b) shows the corresponding reflection response, which again consists of two terms that overlap at and their sum is shown as a dashed green–black line. The rest of is shown in black, whereas the rest of minus the scaled and delayed time‐reversed Green's function is shown in green. Figure 7(c) shows the waves in the entire medium and the arrows indicate the direction of propagation at the acquisition surface. Similar to the missing up‐going wave in the bottom of Fig. 5(c), here the up‐going wave where the dashed and dotted lines coincide is absent. Therefore, similar to what we saw in Fig. 5(b), there is an overlap at zero time of the two terms in the left‐hand side of equation (22), but otherwise the two parts of the reflection response are disjoint.

Figure 7

The time‐reversed response of the actual medium to ;(a) the incident wave field, ; (b) the reflected wave field composed of (black line) and (green line) at of a time‐reversed experiment; (c) the wave field propagating in reverse time inside a medium with five reflectors as a function of vertical travel time and recording time.

To facilitate the derivation of the Marchenko equation and later analysis of non‐recursive inversion, we change the notation of the fundamental wave fields to allow positioning the down‐going impulse at an arbitrary vertical travel time ζ. We write them as with ζ as a free parameter that defines the time window in which the function exists. Hence, for and and for . With this new notation, we write equations (18) and (22) in the time domain as where for .

Derivation of the Marchenko equation

For the later purpose of Marchenko inversion, we derive the Marchenko equation. We define a shifted version of the fundamental wave field as for , where it is noted that the impulse at zero time, which is the first term in , is not included in . We use this wave field in equations (24) and (25), take the time‐reverse of the latter equation and restrict the time window such that the contributions from the Green's functions are excluded. This means we evaluate the equations in the interval and find where shows up as a separate term in the left‐hand side of equation (26), because the initial delta‐function of is not part of the function but it does contribute to the convolution in equation (24). We are free to add or subtract the two equations and decide to subtract. We introduce the fundamental wave field and end up with and . Equation (28) is the Marchenko equation. The functions are the same as the functions as defined in equations (21) and (22) in Slob et al. (2014), but they used z 0 instead of and used a depth level instead of vertical travel time ζ. With this result, the relation between the fundamental wave fields and the kernel of the Marchenko equation is established.

FOUR INVERSION ALGORITHMS

In this section, we describe the inverse filtering method of Kunetz (1964) who used a forward recursive scheme and of Robinson and Treitel (1978) who used a backward recursive scheme. We then use the Marchenko equation to directly compute the impedance for any vertical travel time ζ as given in Berryman and Greene (1980). We end the section with the non‐recursive scheme to obtain the local reflection coefficient at any vertical travel time. All four methods presented here are exact for infinite bandwidth data. As we show below, the direct Marchenko inversion computes the impedance from the zero‐frequency value in the data, which normally is not available. The other three methods offer a data filtering technique that can be used when the source wavelet is known from pre‐processing, because they do not rely on the zero‐frequency information to be present in the data. Presence of the source time signature in the data brings band‐limitation and the associated limited resolution in the proper retrieval of reflection coefficients. In one‐dimension (1D), the information available and retrievable from the reflection response is at best the acoustic impedance as a function of vertical travel time. All four schemes perform this task through direct data filtering methods. In this sense, full waveform inversion using only the subsurface reflection response of a discretely layered lossless 1D medium is a linear problem. The presence of the pressure‐free or rigid surface as acquisition surface does not increase the number of unknowns and slightly modified versions of the schemes presented here will remain valid and exact. Kunetz (1964) and Robinson and Treitel (1978) have included those surfaces in their analysis. Singh et al. (2015) has included it for the Marchenko scheme. It follows from the analysis in Zhang and Slob (2019) for the Marchenko‐type scheme.

Kunetz' inversion method by forward recursion

Kunetz and d'Erceville (1962) derived the recursive expressions for the fundamental wave field and included the effects of a pressure‐free or rigid surface. Kunetz (1964) recognized that the fundamental wave field can be used in a recursive manner to remove overlap from shallow multiples from the first primary that occurs below that shallow part. His algorithm includes the free surface, but here we assume it is not present. The reasoning is as follows. Because it is a forward recursive scheme, we compute filters from previous inversion results and find the reflection coefficient from the next primary in the data. Hence, all reflection coefficients are obtained from the data, not from the filter. The first reflection in the data is a primary with the local reflection coefficient as amplitude, hence r 0 is found directly from the first event in and we find it as where ε is a small number, but large enough to ensure that the reflection occurs in the non‐zero part of the windowed reflection response as expressed in the right‐hand side of the equation. In practice, half the time length of the source wavelet is used. We then know r 0 and t 0 with which the associated fundamental wave field is computed according to equations (11) and (12). We also know and perform the following recursive steps for . We compute the next reflection coefficient from equation (18) as where . Note that in the right‐hand side of equation (30) we convolve the data with , while the truncation is around . The convolution with removes all multiples generated at any of the reflectors from the surface to from the data. In the convolution result, the first event after is the primary reflection from the reflector at and this is the first event in the scaled and delayed Green's function in equation (18). This has the cumulative two‐way transmission coefficients in its amplitude and that is why this factor is present in the denominator of equation (30). Then is computed from equations (9) and (10), we add 1 to n and evaluate equation (30) again until the entire trace is predicted. Note that the denominator in the right‐hand side of equation (30) can be found directly using equation (15). The impedance of each layer is found from and this completes the inversion.

Based on this same idea, a slightly different version of the forward recursive scheme was implemented in Robinson (1967) and was derived in detail in Robinson and Treitel (1977) who also gave a numerical example. Kunetz (1964) showed that the acoustic response at depth, generated by an impulsive source at the free surface, is a function whose autocorrelation has a causal part that is the reflection response that would be measured at the surface. By reciprocity, as was already remarked by Robinson and Treitel (1977), this is the same as what Claerbout (1968) wrote.

Robinson and Treitel's inversion method by backward recursion

Robinson and Treitel (1978) recognized that all local reflection coefficients are present at their correct two‐way travel time in the up‐going part of the fundamental wave field computed for the bottom reflector, . They write that if higher order products of reflection coefficients can be neglected, we need to solve only for . They called those corrupted primaries and we use their terminology. From what we have seen above, when a model has just 11 reflectors, the number of events in is 1024 and only 11 of them are the desired primaries. All the other 1013 will be present in the same time window, overlap the physical primaries and their large number will outweigh their individual small strength. Hence, it is useful to not neglect the non‐physical primaries. Robinson and Treitel (1978) gave a backward recursion scheme that uses the reverse of equations (9) and (10), given by With these equations, we need to start in the bottom and work our way up. Backward recursion implies that the data are used to compute the filter that corresponds to the bottom reflector and find the reflection coefficient directly from the latest event in the up‐going filter . We then use that inversion result to compute the filters that corresponds to next higher reflector and find the reflection coefficient in the same way as the previous one. Hence, the data are used only to compute the filters and each reflection coefficient is found directly from the up‐going part of the filter after each recursion step. To find an equation from which can be found from the reflection response, we convolve equation (5) with the fundamental wave field This leads to where denotes the autocorrelation of the reflection response. Note that equation (34) is equal to equation (22) for , but written such that is avoided. This is possible only for , because then . The second term in the right‐hand side of equation (34) is zero for positive times. We write the down‐going fundamental wave field as , where contains all unknown waves that are emitted after the initial impulse. Because this is a causal function, we evaluate equation (34) for positive times only. This leads to This equation can be solved for for all positive times available. Once this wave field is known, equation (7) is used to determine . Within that wave field, the last event contains the local reflection coefficient of the bottom reflector in the data. It is found as from which and the two‐way travel time, , are known. With this information, the scheme then determines from equations (32) and (33). The last event in contains the reflection coefficient at its two‐way travel time . In this way, we recursively move upward in the data until the local reflection coefficients and two‐way travel times of all reflectors are computed. Impedances are obtained using equation (31) and the inversion is complete.

Marchenko impedance inversion by non‐recursive filtering

Equation (35) resembles equation (28) remarkably well, even though the equations have quite different interpretations and meaning. Equation (28) uses only the reflection response and can be solved for the wave field within its time window of validity for any constant value of ζ. The impedance at that particular vertical travel time is directly obtained by evaluating (Berryman and Greene 1980) The value 1 in the right‐hand side comes from the fact that the initial impulse is not part of the wave field but contributes to the impedance. With this evaluation the inversion is complete. The integral can be seen as a Fourier transformation at zero frequency, which is the only frequency used to compute the impedance.

Marchenko‐type inversion by non‐recursive filtering

Equations (24) and (25) can be used in two ways to perform the inversion. For both, the first step is to use them together to simultaneously solve for for any fixed value of ζ, let us say , by limiting the evaluation of the equations to . These are given by Equations (38) and (39) are the same as equations (24) and (25) in Slob et al. (2014) with a time shift for the fundamental wave fields. The first method retrieves the reflection coefficient as the last event in at its physical two‐way travel time, , or, We observe that this is similar to the backward recursive scheme, but here it is not recursive and we can evaluate at any vertical travel time. Of course, from one reflection coefficient, only the impedance ratio can be obtained, cf. equation (31). The second method uses to evaluate equation (24) for larger values of t. The first event in will be the next primary reflection event in the data, it will have its physical amplitude, , and will be present at its two‐way travel time, . Equation (15) is evaluated to determine the factor after which the reflection coefficient is known. It is given by We observe that this step is similar to the forward recursive scheme, but here it is not recursive and we can perform this step at any vertical travel time.

In this method, we combine the ideas of the forward and backward recursive schemes and make it non‐recursive by using the idea behind the Marchenko equation. As a consequence, only this method can be performed with finite frequency bandwidth data in a target‐oriented manner. Once we have made the first non‐recursive step at a chosen vertical travel time, we can explore the target zone by using the recursive relations of equations (32) and (33) to move in upward direction, or by using the recursive relations of equations (9) and (10) to move in downward direction. The benefit of this possibility is illustrated in the next section for 3D wave fields in a horizontally layered medium.

TARGET‐ORIENTED INVERSION WITH 3D FINITE FREQUENCY BANDWIDTH WAVE FIELDS

When one shot gather of the reflection response of a horizontally layered medium is available, we transform it to the horizontal‐slowness intercept‐time domain, with radial slowness, s, and intercept time, τ. We assume that the up‐going pressure is known at the acquisition surface together with the source time signature, , which we assume to be zero phase for simplicity. The source time signature has a finite‐frequency bandwidth and a zero mean. The Ricker wavelet used in the above examples qualifies. Real source time signatures often have smaller bandwidths than the Ricker wavelet. The up‐going pressure is then given by The functions convolved with the wavelet are denoted . Equations (38) and (39) are now written as and the equations can be solved for each value of s separately. Note that the filters remain impulse response functions, within the bandwidth of the data. These two equations can be solved for at any chosen vertical travel time, ζ, at any value of s. Suppose we solve the equations for several values of s and for and . Then we know and the latest event in is the event with amplitude and intercept time , in which is the normalized vertical slowness in layer . Hence, the first inversion method obtains We can now move in upward direction using the backward recursion of equations (32) and (33), which remain valid by replacing by for each value of s for propagating waves.

For the second method, we use the band‐limited version of equation (24) for later times and obtain . The earliest event in is the reflection from and we know its intercept time , but not yet the local reflection coefficient, because all transmission effects are still present. We use the band‐limited version of equation (25) to compute and obtain the local reflection coefficient of the deeper reflector as (Slob et al. 2014; Wapenaar et al. 2014) We can move in downward direction using the forward recursion of equations (9) and (10), which remain valid by replacing by for each value of s for propagating waves .

We now outline a different target‐oriented inversion procedure than given in Slob, Wapenaar and Treitel (2018). After one Marchenko solution and one Green's function computation, we know the intercept times inside layer n as , and the reflection coefficients and . We first use the intercept time to recover the thickness of the layer. Suppose the ray parameter is sampled with samples as with and . The layer thickness is found from the intercept times and the slowness values. The general expression for the intercept time is given by with known and we rewrite equation (47) as which leads to with . We find the velocity in the layer from With this velocity we invert the reflection coefficient for and is obtained as final parameter. The reflection coefficient is written as where is ratio of the impedances at the two sides of the reflector and is obtained from as and is the ratio of the squared q‐factors given by Because is known, we find from Now both velocities are known and the density ratio can be computed from equation (52). Because is also known, we can repeat this analysis and find and . With this result the inversion is complete. This can be repeated by performing a forward or a backward recursion step, in which case we would continue with the method of Kunetz or Robinson and Treitel, respectively, without having to perform more Marchenko steps.

NUMERICAL EXAMPLE

The methods described above are exact methods for piecewise constant one‐dimensional layered models. With impulse reflection response data they return all layer impedance values without error. When we use finite frequency bandwidth reflection response data, by convolving the impulse reflection response with a Ricker wavelet, Marchenko impedance inversion does not work anymore. When noise is added to the data the backward recursion method soon does not work anymore. The primary reason is that the autocorrelation of the reflection response is used to compute the fundamental wave field. The overall noise amplitude in the down‐going part of the fundamental wave field that we obtain is too high to estimate the up‐going part of the fundamental wave field to start the inversion in the bottom of the model. We do not show results for these two methods. We illustrate how the method of Kunetz (KI), of equation (30), and the Marchenko‐type method (MI), of equation (40), perform on computed data with multiplicative noise in a model with 12 reflectors. Once the reflection coefficients and their two‐way travel times are obtained, we use equation (31) recursively to compute the impedance as a function of two‐way travel time. Using the same model, we show how well the target‐oriented Marchenko‐type inversion (TOMI) scheme works.

The model information is given in Table 2. The source emits a 30 Hz Ricker wavelet with which we compute the reflection response. For illustrating TOMI, we compute the response for 10 different slowness values and use multiplicative noise. We compute the noise in the frequency domain with a base amplitude of 0.1 and a random phase after which it is multiplied by and added to the reflection response. This is equivalent to creating a random white noise trace in time domain, convolving it with the reflection response and then adding it to the reflection response. The time domain noise has extreme values of ±0.009. We compute 1000 different noise realizations and add each to the reflection response. We have tested the inversion schemes with only additive noise, and with both additive and multiplicative noise, and the results statistics are the same as the results shown here for multiplicative noise.

Table 2

Values for the medium parameters in the model for the numerical example

Layer number	Velocity (m/s)	Density (kg/m3)	Thickness (m)	Z/Z_0 (−)
0	1700	1430	∞	1.000
1	2300	2750	64.0	2.602
2	1900	2000	49.2	1.563
3	1700	1430	153.6	1.000
4	2100	1750	216.4	1.512
5	3200	2930	335.9	3.857
6	2000	1750	55.7	1.440
7	2100	2110	154.8	1.823
8	3300	1970	85.4	2.674
9	2500	2110	199.0	2.170
10	3000	2110	162.6	2.604
11	2500	2250	147.7	2.314
12	2900	2300	∞	2.744

For each of these 1000 datasets, we compute inversion results with KI and MI. In the inversion, we first compute the reflection coefficients for each method. For both methods, we perform blind inversion as a fully automated process. The only model assumption is that all layers have constant parameters. We do not assume a predetermined number of layers, but we assume that the layers are thick enough to prevent resolution problems to occur. The 30 Hz Ricker wavelet is 62 ms long to an amplitude level of 0.0023. This means we assume that the two‐way travel time inside each layer is 62 ms or more. In the model we use in this example, four layers have a two‐way travel time between 51 ms and 56 ms (layer numbers 1, 2, 6 and 8). In the 51 ms window, the smallest wavelet amplitude is 0.04, which is more than four times the noise level.

The KI method is implemented as follows. We search the extreme value in the reflection response within the first 62 ms time window. When this value is at or below the noise level, we assume it is noise and search again by shifting the search window by one window length. When this value is above the noise level, we assume it is part of a reflection and search for another extreme value within the same window length but with half a window length shift. When that extremum is larger than the one found earlier, its location and value provide the time location and reflection coefficient value. Otherwise the location and value of the previously found extremum provide the time location and reflection coefficient value. We then apply the forward recursion step and search the resulting trace by shifting the time search window by one window length from the detected reflector location. This is repeated until the end time is reached. With this method, the number of detected reflectors has a mean value of 12.1, a maximum of 19 and a standard deviation of 0.4.

The MI method is implemented in two steps. We first solve the coupled Marchenko‐type equations to obtain the primary reflectivity as a function of travel time. The resulting trace contains only primary reflections convolved with the Ricker wavelet. On that trace we proceed as we did in KI but now without the need to use the recursive steps, because all multiples have been eliminated already. We use a threshold of 0.04 to decide whether an extremum is part of a reflection or noise. This higher level can be used because in the trace to be analysed the transmission effects have been eliminated and true reflectivity amplitudes are present. The number of detected reflectors has a mean of 12.5, a maximum of 21 and a standard deviation of 1.1. For both results, we compute the constant impedance values for each layer for all travel times within that layer to plot the results.

Figure 8 shows the mean normalized impedance values and the value at one standard deviation away from the mean as a function of travel time as obtained from KI in Fig. 8(a) and from MI in Fig. 8(b). The impedance of the upper half space is used as normalization factor. The true values are shown in black solid lines. The KI retrieved values are shown in red with the mean by a dashed line and the values at one standard deviation by dotted lines. The MI results are shown in a similar way in blue colours. It can be seen that the mean values are quite accurately obtained for both methods. The mean value errors are below 0.5% for KI and MI. The mean timing errors are below ±4 samples for both methods. In Fig. 8(a), we can see that the standard deviation of the impedance values suddenly jumps up in layer with number 6, which is where the first error in the arrival time estimations occur. The standard deviation grows in the next two layers but then stabilizes until the end of the trace. In Fig. 8(b), we can see that the standard deviation of the impedance values slowly grows and continues to grow after the last reflector is found. From reflector 6 onward in this model, KI shows larger standard deviations than MI but at the end of the trace the standard deviations are the same.

Figure 8

The retrieved mean impedance value and at one standard deviation for KI (a) and MI (b). The mean values are normalized to the impedance of the upper half space. The true normalized impedance is shown by solid black lines, the mean impedances by dashed lines and at one standard deviation by dotted lines, in red for KI and in blue for MI.

Both methods retrieve the impedance of each layer with acceptable errors in the mean values of the retrieved results. This is illustrated in Fig. 9 where it can be seen that the erratic behaviour of the error in the impedance between 1.5 seconds and 2 seconds is the effect of the noise that causes the methods to detect a reflector. It looks so erratic because at every noise realization those non‐physical reflectors are detected at slightly different times. They do not create large changes in the impedance. Note that the error axis is bounded by ±0.5% error. The horizontal lines are spiky errors that occur after 0.8 seconds. They all coincide with small timing errors of the reflector locations. These errors are similar in both methods, below 10% in amplitude and less than 4 ms in time location.

Figure 9

The error in the mean values of the normalized impedance as a function of travel time for KI (red solid) and MI (blue dashed). The horizontal lines in both results indicate large errors (<10%) arising from small errors in the retrieved travel times of the reflectors.

The coupled Marchenko equations can be used at any two‐way travel‐time value for TOMI. We assume that the reflection response is obtained such that we have data at zero incidence, , and at nine angles of incidence, , in the upper half space from 15° to 28°, such that the radial slowness is given by . We choose the same noise levels and compute 1000 realization of the reflection responses at these angles of incidence. For each dataset, we perform the inversion as outlined in equations (49)–(54) by computing the mean value using all values for the radial slowness. We choose the truncation time values inside layer number 6 (). The latest event in has amplitude and the time of arrival is . As explained in the text above equation (46), we then use the band‐limited version of equation (24) for and find the time of arrival of the first event in , which is equal to . We now know the interval intercept time . The arrival times are not significantly influenced by the noise after one Marchenko step. This means that the layer thickness, d 6 and the corresponding velocity c 6 are found with mean value errors well below 1% and very small standard deviations. The true and mean retrieved reflection amplitude r 5 from is shown as a function of incidence angle in Fig. 10. The vertical lines indicate the retrieved results within one standard deviation. The figure shows that the errors in the mean retrieved value and the standard deviation increase with increasing angle of incidence. The maximum error in the mean value is below 1% and occurs at 28°. The values at one standard deviation have an error below 3%. The numerical values are summarized in Table 3.

Figure 10

The true and mean retrieved reflection amplitude r 5 as a function of incidence angle; the vertical lines indicate the standard deviation of the retrieved results.

Table 3

Values for the medium parameters in the model and for the TOMI result

	d 6 (m)	c 6 (m/s)	c 5 (m/s)	ρ_6/ρ_5 (‐)
True	55.7	2000	3200	0.594
Mean	55.9	2007	3209	0.595
Std	0.2	8	86	0.023

DISCUSSION

The method of Kunetz (KI) is a forward recursive method and involves only convolutions of modelled fundamental wave fields with the data. The fundamental wave fields are modelled using equations (9) and (10). The only influence of the noisy data on the fundamental wave fields are the amplitude and arrival time of the estimated reflection coefficient. For this reason, this method is quite well behaved with regard to noise. Estimated reflection coefficients at shallower depth are used to compensate for the transmission effects for each new estimated reflection coefficient. Small amplitude errors lead to increasingly larger errors because of increasingly incorrect correction for transmission effects. Those errors and errors in the estimated arrival times of the reflections lead to incorrect times of the events in the fundamental wave fields. These lead in turn to incomplete cancelation of the multiple reflections. When the amplitudes of these remnant multiple refections are above the noise level, they lead to the estimation of non‐physical reflectors. The method suffers from increasing errors in estimating reflection coefficients and arrival times, because these are used to model fundamental wave fields, which in turn are used to detect and estimate the arrival time and reflection coefficient of a deeper reflector. These effects are minimal in the example shown here. The advantage of KI is that no equation has to be solved. It is consists of one convolution and modelling the fundamental wave fields at each step. This makes it a very fast method.

The Marchenko‐type inversion is a non‐recursive method that computes the fundamental wave fields by filtering the reflection response. Filtering involves a convolution and a correlation of the fundamental wave fields with the data. As a consequence, the noise in the data enters into the fundamental wave fields. The correlation of two noisy traces increases the noise in the result. This can be seen by the increase of the relative standard deviation in the impedance values obtained as a function of travel time with this method. An advantage of Marchenko‐type method (MI) is that it can be adapted to accommodate an unknown source wavelet, but then the down‐going part of the field must be known (Ravasi 2017). Another advantage of MI is that higher thresholds can be used to detect a reflector, because the local reflection coefficient values are present in the up‐going part of the fundamental wave field. This is necessary because the correlations make the noise levels grow with increasing travel time for the inversion. The larger standard deviation in the number of detected reflectors in MI than in KI is because the noise levels increase with increasing inversion times. Most of the reflectors detected in addition to the reflectors in the model occur after the last reflector of the model has been detected. Another implication of higher noise levels is that weaker reflectors whose reflections have large two‐way travel times will be missed by MI. This does not happen in the example shown here.

We evaluate the target oriented inversion method on noisy reflection responses for normal incidence and nine non‐zero values for the radial slowness. We use the interval intercept times at the available angles of incidence to find the value for the layer thickness for each noise realization. Then the velocity inside the layer is obtained from the layer thickness and the normal incidence intercept time. Estimating the intercept times of the reflections from the top and bottom boundaries of the layer can be done with high accuracy. This means that the estimate of the intercept times of the layer are very accurately obtained. For this reason, both the layer thickness and the velocity in the layer are obtained with high accuracy and small standard deviations. Figure 10 shows that the mean values of the retrieved reflection strength is very accurate for all values of the incidence angle and the largest standard deviation is below 3% of the reflection coefficient value. For this reason and the fact that c 6 is accurately obtained from the intercept times, the mean value of the velocity in the layer number 5 and the density ratio of the layers 5 and 6 are very accurately obtained as well. Because of the increased standard deviation in the reflection coefficient, the standard deviation in the retrieved values of c 5 is higher, but still acceptable at less than 3%. Given the lower sensitivity to noise of the Kunetz method, it seems a good idea to continue with this method recursively after one Marchenko step. The Marchenko method is already adapted to retrieve local reflectivity in three‐dimensional (3D) that could be used for inversion (Wapenaar et al. 2014). Another interesting option seems to investigate the method of Kunetz further to see whether it can be extended for use in a general 3D setting.

CONCLUSIONS

We have discussed two recursive and two non‐recursive methods for computing the reflection coefficients of a discrete layered model. The recursive methods use one equation that involves the data and two model equations that are recursively used to compute the fundamental wave fields. The forward recursive scheme computes the reflection coefficient from the data. The backward recursive scheme computes the reflection coefficients from the up‐going part of the fundamental wave field. The non‐recursive methods use only equations that involve the data to compute the fundamental wave fields at any chosen travel time. The first non‐recursive method is the Marchenko impedance inversion method. The method first computes a particular sum of the up‐ and down‐going parts of the fundamental wave fields. It then computes the impedance directly from the zero‐frequency value of that sum of the fundamental wave fields. The second is a Marchenko‐type method that first computes the up‐ and down‐going parts of the fundamental wave field from the data and then computes the reflection coefficient. The first three methods were developed independently almost 60 years ago. We have shown that these were revived briefly several times during the six following decades, but remained seen as independent. We have shown from the underlying physics that all four methods are based on the concept of the fundamental wave field.

In the Marchenko‐type method, the two equations are solved as a coupled set of equations for a chosen vertical travel time. We can then proceed in two ways. The first is to find the local reflection coefficient directly in the up‐going part of the fundamental wave field. It is the coefficient of the reflector that has a vertical travel time less than but closest to the chosen time instant. The second is to compute the local reflection coefficient from the first event in the Green's function. The fundamental wave fields are used to correct its amplitude for two‐way transmission effects. Either way, this scheme leads to the possibility of performing target‐oriented inversion. We have shown for a one‐dimensional model that when data are available for several ray parameters, the non‐recursive target‐oriented inversion can provide layer velocity, layer thickness and the ratio of densities of two adjacent layers. The non‐recursive nature of this method makes it a good candidate for extension to laterally varying media. Its sensitivity to noise can create problems at large times. The method of Kunetz is not very sensitive to noise (additive and/or multiplicative) and is very fast compared with the Marchenko‐type method, although it can be implemented only in a recursive way. The theory of the Marchenko‐type method is available for three‐dimensional (3D) inversion. The concept of forward recursive filtering seems worth investigating to see whether it can be adapted for use in a general 3D setting.

DATA SHARING

The data that are used to make the plots of Figs. 8, 9, 10 can be obtained from the first author on request.