SciELO - Scientific Electronic Library Online

 
vol.8 issue2OPTIMAL CODING OF BLENDED SEISMIC SOURCES FOR 2D FULL WAVEFORM INVERSION IN TIME author indexsubject indexarticles search
Home Pagealphabetic serial listing  

Services on Demand

Journal

Article

Indicators

Related links

  • On index processCited by Google
  • Have no similar articlesSimilars in SciELO
  • On index processSimilars in Google

Share


CT&F - Ciencia, Tecnología y Futuro

Print version ISSN 0122-5383

C.T.F Cienc. Tecnol. Futuro vol.8 no.2 Bucaramanga July/Dec. 2018

https://doi.org/10.29047/01225383.87 

Artículos originales

DIRECT AND INDIRECT INVERSION AND A NEW AND COMPREHENSIVE PERSPECTIVE ON THE ROLE OF PRIMARIES AND MULTIPLES IN SEISMIC DATA PROCESSING FOR STRUCTURE DETERMINATION AND AMPLITUDE ANALYSIS

INVERSIÓN DIRECTA E INDIRECTA Y UNA PERSPECTIVA NUEVA Y COMPLETA DEL ROL DE LAS REFLEXIONES PRIMARIAS Y LOS MÚLTIPLES EN EL PROCESAMIENTO DE DATOS SÍSMICOS EN LA DETERMINACIÓN ESTRUCTURAL Y ANÁLISIS DE AMPLITUD

Arthur-B Wegleina  * 

a M-OSRP/Physics Dept./ University of Houston, Houston, TX, 77204, USA.


ABSTRACT

The removal and use of multiples have a single shared goal and objective: the imaging and inversion of primaries. There are two kinds of primaries: recorded primaries and unrecorded primaries. For imaging recorded primaries using an industry standard practice smooth velocity model, recorded multiples must be removed, to avoid false and misleading images due to the multiples. Similarly, to find an approximate image of an unrecorded primary, that is a subevent of a recorded multiple, unrecorded multiples that are subevents of the recorded multiple must be removed, for exactly the same problem and reason that recorded multiples are needed to be eliminated. Direct inverse methods are employed to derive this new comprehensive perspective on primaries and multiples. Direct inverse methods not only assure that the problem of interest is solved, but equally important, that the problem of interest is the relevant problem that we (the petroleum industry) need to be interested in.

KEYWORDS: Direct and indirect inversion; Removing and using multiples; Migration; Inverse scattering series; ISS Q compensation without knowing Q

RESUMEN

La remoción y el uso de reflexiones múltiples tienen una sola meta y objetivo común: la construcción de imágenes (imaging) e inversión de reflexiones primarias. Existen dos clases de reflexiones primarias: las registradas y las no registradas. Dentro de la práctica estándar de la industria, se construyen imágenes de las reflexiones primarias registradas, empleando un modelo de velocidad suave, donde se deben remover los múltiples registrados para evitar los eventos falsos y engañosos provenientes de éstos. De igual forma, para encontrar una imagen apropiada de una reflexión primaria no registrada, o sea un sub-evento de un múltiple registrado, las múltiples no registradas que son sub-eventos de la múltiple registrada se deben retirar debido a exactamente los mismos problemas y razones por los que es necesario eliminar las múltiples registradas. Los métodos inversos directos no solo aseguran que se resuelve el problema en cuestión, sino lo que es igualmente importante, que dicho problema es el problema relevante en el que nosotros (la industria del petróleo) debemos interesarnos.

PALABRAS CLAVE: Inversión directa e indirecta; Remoción y uso de reflexiones múltiples; Migración; Series inversas de dispersión (SID); Compensación del coeficiente de calidad Q

1. INTRODUCTION

In this paper, we provide a new and comprehensive perspective on primaries and multiples, that encompasses both removing multiples and using multiples. We describe the original motivation and objectives behind these two initiatives, viewed almost always as "remove multiples versus use multiples". The premise behind that "versus" phrasing speaks to a competing and adversarial relationship. A contribution in this paper is placing these two activities and interests within a single comprehensive framework and platform. That in turn reveals and demonstrates their complementary rather than adversarial nature and relationship.

They are in fact after the same single exact goal, that is, to image primaries: both recorded primaries and unrecorded primaries. There are circumstances where a recorded multiple can be used to find an approximate image of an unrecorded subevent primary of the recorded multiple.

All direct methods for imaging and inversion require only primaries as input. To image recorded primaries requires that recorded multiples must first be removed. To try to use a recorded multiple to find an approximate image of an unrecorded primary subevent of the recorded multiple requires that unrecorded multiple subevents of the recorded multiple be removed. All multiples, recorded multiples and unrecorded multiples need to be removed. Not removing those recorded and unrecorded multiples will produce imaging artifacts and false and misleading images, when seeking to image recorded and unrecorded primaries, respectively.

For indirect methods that either: [1] solve a forward problem in an inverse sense, like AVO or [2] are model matching methods like, e.g., FWI, any data can be forward modeled and solved in an inverse sense, or model matched, respectively. In contrast, for direct inverse methods the data required and the algorithms called upon are explicitly and unambiguously defined.

In our view, direct and indirect methods each have a role to play, the former where the assumed physics (and an assumed partial differential equation governs the wave phenomena) captures some component of reality and the latter (indirect methods) is the only possible choice for the part of reality that is beyond our physical models, equations and assumptions. Furthermore, it would be ideal if the indirect method and the direct method were cooperative and consistent. That cooperation can be arranged by choosing the objective function or sought after quantity to be satisfied (in the indirect solution) as a property of the direct solution [1].

In this paper, we depend upon the clarity of direct methods to provide a new perspective that advances our understanding of the role of primaries and multiples in seismic exploration. That, in turn, allows us to recognize the priority of developing more effective multiple removal capability within a comprehensive strategy of providing increased seismic processing and interpretation effectiveness. The unmatched clarity and definitiveness of direct inverse methods provides a new, unambiguous and clearer perspective with a timely message on the role of primaries and multiples in seismic processing for structural determination and amplitude analysis.

STATE OF THE TECHNIQUE

DIRECT AND INDIRECT METHODS FOR STRUCTURAL DETERMINATION AND AMPLITUDE ANALYSIS

The starting point of our discussion of primaries and multiples begins with the key definitions and classification of direct and indirect inversion methods.

Inverse methods are either direct or indirect (see, e.g., the definition and examples of direct and indirect inversion in e.g., [2], [3]. Direct methods provide assurance and confidence, that we are solving the problem of interest. For example the direct solution of the quadratic equation ax 2+bx+c=0 has roots x=(-b±√(b2-4ac)/2a). Nobody would consider an indirect solution of the quadratic equation. The clear logic and reasoning behind choosing a direct solution of the quadratic equation [where indirect solutions of guessing roots, and matching and minimizing cost functions, e.g., to seek values of x to minimize (ax 2+bx+c)n n or integrals of such an expression] carries over to all math and math-physics problems (including inverse seismic problems) wherever a direct solution exists.

HOW TO DETERMINE WHETHER A PROBLEM OF INTEREST IS THE PROBLEM WE (THE PETROLEUM INDUSTRY) NEED TO BE INTERESTED IN?

In addition to knowing that we are solving the problem of interest, and equally important, direct solutions communicate whether the problem of interest is the problem that we (the petroleum industry) need to be interested in. When a direct solution does not result in an improved drill success rate, we know that the problem we have chosen to solve is not the right problem since the solution is direct and cannot be the issue. On the other hand with an indirect method, if the result is not an improved drill success rate, then the issue can be either the chosen problem, or the particular choice within the plethora of indirect solution methods, or both.

Among key aspects in effectively designing and managing an industrial or academic research program, are: (1) to be able to identify and select the problems and challenges that need to be addressed and (2) what benefit would derive from a new and effective method that addresses a specific challenge. From our perspective, benefit is measured by an increase in successful exploration drilling and optimizing appraisal and development drill placement. Challenges arise when the assumptions behind current seismic methods are not satisfied. As noted in the previous section, direct methods play a unique role in problem identification, a critical aspect of defining research objectives and programs.

THE DISCONNECT BETWEEN THE SUCCESS RATE OFTEN CLAIMED BY RESEARCHERS AND THE REALITY OF THE DRILL SUCCESS RATE IN DEEP WATER FRONTIER EXPLORATION

In our experience, the most important ingredient in defining challenges, and prioritized open issues that need to be addressed is asking the end-user [e.g., the seismic interpreter, the drill decision makers in operating business units] what methods are [and are not] working, and under what circumstances. These individuals by and large only have the single interest and objective in their focus on effectiveness, avoiding dry-holes and making successful drill decisions. On the other hand, there is too often a serious disconnect between, e.g., the frontier drill success rate of 1 in 10 in the deep water Gulf of Mexico and the 100% success rate typically reported by researchers at international professional conferences and in societal publications. Too often, researchers will simply refuse to recognize major problems and challenges that exist, unless and until they might feel comfortable they can address them. That avoidance of recognizing major challenges and obstacles is what we call "the disconnect". Within that "disconnect" resides a tremendous positive set of opportunities to define major E and P challenges and to address them. Seeking funding can be tricky, since we are directed to research departments for that support. There will often be one individual (or at most a very small number) within a research organization who is (are) able to recognize prioritized and significant pressing challenges and problems and is fascinated rather than frightened by new visions of what might be possible. We have been enormously fortunate for the funding and support we have received and do receive, and we are enormously grateful and deeply appreciative.

RELEVANT RESEARCH PROGRAMS BEGIN BY DEFINING CURRENT SHORTCOMINGS AND ADDRESSING ACTUAL SEISMIC CHALLENGES

In our view a research program needs to begin by defining the actual real world seismic challenges and pressing issues, and the current method shortcomings being addressed | and then developing and delivering methods that address those challenges. A relevant research program must begin with the problem, and then seeks a solution; it does not involve a method looking for a problem.

We encourage and welcome and need new seismic methods and capability. However, in our view, all new ideas for imaging and inversion (e.g., interferometry and Marchenko and virtual sources), need to begin by clearly defining the shortcoming and limitation of current capability that they are addressing. And they need to specifically define the challenge and potential added value relative to the current high water mark of imaging and inversion methods that is, Stolt Claerbout III migration-inversion, for automatically and simultaneously imaging and inverting specular and non-specular reflectors (curved reflectors, diffractors and pinch-outs) see, e.g., [4], [5] and the inverse scattering series (ISS) task specific subseries for depth imaging and inversion, [6], where the former (Stolt CIII) require and the latter (ISS methods) do not require subsurface information, respectively.

There is too often an insular inward looking aspect to research projects that do not deem it necessary to show relevant differential added value relative to current capability. In our view, that is essential | and it is incumbent upon us as researchers to explain where a new advance sits within the seismic toolbox, and the circumstances when it will be (and will not be) the appropriate and indicated choice among method options. The research objective is to increase tool-box options.

All of our research reporting needs to have the method assumptions clearly spelled out in the conclusions, and to delineate the remaining open issues that need to be addressed and will require new concepts and future contributions.

MULTI-D DIRECT INVERSION

The inverse scattering series (ISS) is the only direct inversion method for a multidimensional subsurface. Solving a forward problem in an inverse sense is not equivalent to a direct inverse solution [2] (see Appendix A). Many methods for parameter estimation, e.g., AVO, are solving a forward problem in an inverse sense and are indirect inversion methods. The direct ISS method for determining earth material properties, defines both the precise data required and the algorithms that directly output earth mechanical properties. For an elastic earth model of the subsurface the required data for parameter estimation and amplitude analysis is a matrix of multi-component data, and a complete set of shot records, with only primaries see e.g. [2], [3], [7].

With indirect methods any data can be matched: one trace, one or several shot records, one component, multicomponent data, with primaries only or primaries and multiples, with pressure measurements and displacement and spatial derivatives of these quantities, and stress, or only just multiples. Added to that are the innumerable choices of cost functions, generalized inverses, the often ill-posed nature of indirect methods, and local and global search engines.

Direct and indirect parameter inversion have been compared for a normal incident plane wave on a 1D acoustic model, and full bandwidth analytic data [2], [8], [9]. In the latter example and comparison [even when the iterative linear inverse has a linear approximate provided by the analytic linear ISS parameter estimate], the direct ISS method has more rapid convergence and a broader region of convergence. The difference in clarity and effectiveness between indirect methods and direct methods (where the direct methods specify the data requirements, provide well defined algorithms that produce the linear and explicit and unique higher order contribution to the sought after earth mechanical properties) increases as subsurface circumstances become more realistic and complex, and, in particular with an elastic or anelastic subsurface and with band-limited noisy data [2], [3].

There are two categories of direct methods for imaging and inversion: (1) those that require subsurface information, and (2) those that do not require subsurface information. For Stolt CIII migration, see [5], [10] the most general and effective imaging principle and migration method, a smooth velocity model will suffice for structural determination and reflector location. For more ambitious objectives (using Stolt CIII migration) beyond structural determination, such as amplitude analysis for target identification, all elastic and inelastic subsurface properties need to be provided above the target. For all migration methods, e.g., Stolt CIII and CII RTM or Kirchhoff, in practice a smooth velocity is employed and all recorded multiples must first be removed, to avoid false and misleading images from recorded multiples, before imaging and inverting recorded primaries. Stolt CIII imaging, the current high water mark of migration and migration-inversion capability requires recorded primaries as input.

The ISS is the only direct inversion method for a multidimensional earth (see, e.g., [6]). It can be applied with or without subsurface information. The direct ISS depth imaging subseries reduces to the single term Stolt CIII migration algorithm for the case of adequate subsurface information above the structure to be imaged ([6], [12] and http://mosrp.uh.edu/news/invited-presentation-petrobras-workshop-aug-2016). In the case where there is adequate velocity information above a reflector, all the ISS imaging subseries terms beyond the linear first term (the linear first term in the ISS depth imaging subseries corresponds to Stolt CIII migration) will vanish for imaging that reflector. See the 2017 executive summary video http://mosrp.uh.edu/news/executive-summary-progress-2017. However, the fact that all the ISS derived methods can be formulated and applied directly and without any subsurface information has been its unique strength and advantage, distinguishing it from all other seismic processing approaches and methods.

Once we recognize that "absolutely no subsurface information is required" property of the entire ISS, and every individual term in the series, then that leads to the idea of locating isolated task subseries of the ISS that can perform and achieve every seismic objective directly in terms of data and without subsurface information. There are isolated task subseries that perform free surface multiple removal, then internal multiple removal, followed by distinct series that migrate and invert primaries, and perform Q compensation directly and without subsurface elastic or inelastic information. [See e.g. [13] for Q compensation without knowing, estimating or determining Q.]

The ISS is the only direct inversion methodology for a multidimensional subsurface, it does not require subsurface information and multiples are removed prior to performing the tasks of structural determination and amplitude analysis, the latter inputting only primaries. The only direct inversion method for a multi-dimensional subsurface without subsurface information treat multiples as coherent noise that needs to be removed. If ISS depth imaging and inversion subseries needed multiples it would not have distinct ISS subseries that remove free surface and internal multiples. When the velocity information above a reflector is known, the ISS depth imaging reduces to Stolt CIII migration, and for a smooth velocity model, multiples will cause artifacts and must be removed. All direct imaging and inversion methods [with or without subsurface information] call for an adequate set of primaries, and require as a prerequisite that all multiples be removed.

All direct imaging and inversion methods [with or without subsurface information] call for an adequate set of primaries, and require as a prerequisite that all multiples be removed.

USING MULTIPLES

There has been considerable literature recently on using multiples. We will show (below) the consistency, interrelated nature and precisely aligned objectives of the remove multiples and the use multiples activity. Within that new framework we explain based on the need for imaging recorded primaries the need to remove recorded multiples. That reality drives and defines the need and priority of effective multiple removal. And Stolt CIII stands alone in capability and beyond all migration methods including all RTM methods for delivering the maximally resolved and delineated structure and the most effective amplitude analysis at both simple planar and complex corrugated and diffractive structure (see e.g. [5], [11]). The use of multiples to provide an approximate image of an unrecorded primary, cannot produce a Stolt CIII image of the unrecorded primary, instead it provides a weaker and approximate RTM imaging result. Approximate images of unrecorded primaries extracted as subevents of a recorded multiple do not deliver Stolt CIII imaging and inversion capability that Stolt CIII migration imaging and inversion delivery can only be achieved with recorded primaries. In the sections that follow, we review and exemplify the recent advances in the arenas of removing and using multiples, and describe open issues and challenges that need to be addressed.

A NEW AND COMPREHENSIVE PERSPECTIVE ON THE ROLE OF PRIMARIES AND MULTIPLES IS SEISMIC PROCESSING FOR STRUCTURAL DETERMINATION AND AMPLITUDE ANALYSIS

A major activity within M-OSRP has been and remains the development and delivery of fundamentally new and more effective methods for removing free surface and internal multiples, for offshore and on-shore plays, without damaging proximal or interfering events. That is, the current focus is on removing multiples that interfere with target or reservoir identifying primaries, without damaging the primaries. More effective multiple removal remains an active and priority seismic research topic. That is an essential requirement to be able to derive full benefit from the new Stolt CIII migration-inversion methods that are the currently most effective method at imaging and inverting primaries.

We recognize that there is considerable attention and communication these days on "using multiples". In the note below and in the executive summary video http://mosrp.uh.edu/news/executive-summary-progress-2017 we present a new perspective on the removal and using of multiples.

As we noted, all direct methods for imaging and inversion require a complete set of primaries. However due to limits in acquisition some primaries are recorded and others are not recorded. Primaries are therefore classified as recorded and unrecorded.

To image recorded primaries, with a smooth velocity model, recorded multiples need first to be removed. If not removed, each multiple will produce a false and misleading structural image. For unrecorded primaries, the idea begins with the assumption that there are recorded multiples that consist of two subevents, one recorded and one unrecorded and the latter being an unrecorded primary [14]- [17]. Then the recorded multiple and the recorded subevent of the multiple are used to find an approximate image of the unrecorded primary that is a subevent of the recorded multiple. However the unrecorded subevent of the recorded multiple that we assume is an unrecorded primary, might in fact be an unrecorded multiple, and not as assumed an unrecorded primary. Any unrecorded multiple that is a subevent of the recorded multiple must be removed to avoid it producing a false and misleading structural image. See Figure 1.

Figure 1 Using a recorded multiple to find an approximate image of an unrecorded primary of the multiple: illustrate the need to remove unrecorded multiples. A solid line () is a recorded event, and a dashed line (- - -) connotes an unrecorded event. 

Hence, to image recorded primaries recorded multiples must first be removed, and to find an approximate image of an unrecorded primary requires unrecorded multiples to be removed.

The very use of multiples speaks to the primacy of primaries. Multiples are only useful if it contains as a subevent an unrecorded primary. A multiple that has all its subevents recorded has absolutely no use or value. All primaries are useful | and there is no substitute for a complete set of recorded primaries. Multiples can at times be useful, but are not in any sense the "new primary".

The recorded multiple event that can be used (at times) to find an approximate image of an unrecorded primary, must as an event be removed in order to image recorded primaries. The removing and using of multiples are always about our interest in primaries, both recorded and unrecorded primaries that we seek and require, and hence removing and using multiples are not adversarial, they serve the same single purpose and objective.

In the Executive Summary Presentation in the link there is a detailed discussion on this new perspective regarding removing and using multiples.

Here is the link with the executive summary video: http://mosrp.uh.edu/news/executive-summary-progress-2017.

Basically: (1) to image recorded primaries, with a smooth velocity model, recorded multiples must be removed and (2) for unrecorded primaries, to use a recorded multiple and a recorded subevent of the multiple to find an approximate image of an unrecorded primary subevent of the recorded multiple, any unrecorded multiple that is a subevent of the recorded multiple must be removed.

Hence, to image recorded primaries recorded multiples must be removed, and to find an approximate image of an unrecorded primary requires unrecorded multiples to be removed. The recorded multiple event that can be used (at times) to find an approximate image of an unrecorded primary must as an event be removed in order to image recorded primaries.

The key point is that it is primaries, both recorded and unrecorded primaries that we seek and require, and removing and using multiples are not adversarial, they serve the same single purpose and objective: the imaging of primaries. Multiples (recorded and unrecorded) need to be removed in order to image primaries (recorded and unrecorded, respectively).

If the multiple does not contain an unrecorded primary subevent, then is has no use. Whether or not an unrecorded primary is within the recorded multiple determines whether the recorded multiple is or is not useful. The use or lack of use of the multiple depends on whether or not a specific and particular primary has not been or has been recorded. What use is a multiple where all primary subevents of the multiple have been recorded. The answer: absolutely no use or value, none whatsoever the only interest for us in such a multiple is (as always) to remove that recorded multiple in order to not produce false, misleading and injurious images when migrating recorded primaries.

Hence multiples are NOT now rehabilitated events on equal footing with recorded primaries. They are NOT the new primaries and multiples are never migrated (That idea and thought of "migrating multiples" has no meaning. See [18], but as events themselves must always be removed. For those pursuing the use of multiples, it is of interest to know how unrecorded multiples will be removed.

The use of multiples is worthwhile to pursue, and to develop and deliver. Their value directly depends on the lack of adequate recorded primaries. There is no substitute for recorded primaries for the extraction of complex structural information and subsequent amplitude analysis. The high water mark of migration capability, Stolt CIII migration for heterogeneous media see e.g., [5], [10], requires recorded primaries. Methods that use a recorded multiple to obtain an approximate image of an unrecorded primary subevent, cannot achieve a Stolt CIII migration delivery and resolution effectiveness under any circumstances. The greatestdi_erential added-value [compared to all other migration methods] derives from Stolt CIII migration for complex structure determination and subsequent amplitude analysis, see [5], [11]. The priority of recorded primaries drives the priority of removing recorded multiples. In the next several sections we review the status, recent advances and open issues in removing multiples.

MULTIPLE REMOVAL: A BRIEF HISTORIC OVERVIEW AND UPDATE ON RECENT PROGRESS AND OPEN ISSUES

Multiple removal has a long history in seismic exploration. Among early and effective methods for removing multiples are CMP stacking, deconvolution, FK and Radon Filtering. These methods made assumptions about either: (1) the statistical, random and periodic nature of seismic events, (2) the ability to determine an accurate velocity model, (3) the assumed move-out differences between primaries and multiples, and (4) subsurface information including knowledge about the reflectors that generate the multiples.

However, as the industry trend moved to deep water and ever more complex offshore and on-shore plays, the assumptions behind those methods often could not be satisfied and therefore these methods were frequently unable to be effective and failed.

Methods that sought to avoid those limiting assumptions include SRME for free surface multiples and the distinct inverse scattering subseries (ISS) for removing free surface and internal multiples. SRME did not require subsurface information but only predicted the approximate time and approximate amplitude of first order free surface multiples at all offsets. In contrast, the ISS free surface multiple removal algorithm does not require subsurface information and predicts the exact time and exact amplitude of all orders of free surface multiples at all offsets. A quantitative comparison of SRME and the ISS free surface multiple elimination [ISS FSME] algorithm can be found in [19], [20]. That analysis helps to define when SRME and ISS free surface multiple elimination are the appropriate and indicated choice within the multiple removal seismic toolbox. SRME relies on an energy minimization adaptive subtraction to fill the gap between its amplitude and time prediction and the amplitude and time of the free surface multiple. That adaptive energy criteria assumes that there is less energy, in an interval of time, when a multiple is removed compared to when it is present. That assumption can and will fail for interfering or proximal events. The ISS free surface multiple elimination method does not require an adaptive energy subtraction, and, hence, is effective whether or not the multiple is isolated or if it is proximal or interfering with other events.

A key and central objective in multiple removal is not to damage target and reservoir primaries. For isolated free surface multiples SRME can at times be a reasonable tool box option. However, for free surface multiples that are proximal to (or interfering with other events), e.g., primaries, the ISS free surface multiple elimination algorithm is an important option and could be the appropriate and indicated choice. The ISS free surface multiple elimination requires the direct wave and source and receiver ghosts to be removed. Later in this paper, we will provide references that utilize variants of Green's theorem to remove the direct wave and ghosts, without damaging the reflection data.

The inverse scattering series internal multiple attenuation algorithm [6], [21], [22] is, at this time, the only internal multiple algorithm that does not require any subsurface information, no knowledge of the multiple generators and no seismic interpreter intervention. It is a multidimensional method that predicts the exact time and approximate amplitude of all internal multiples at all offsets. It is the current high water mark of internal multiple capability in the petroleum industry. In the sections below we review and exemplify the free surface and internal multiple removal status and describe recent advances in internal multiple elimination, and open issues.

THE ISS FSME AND SRME EQUATIONS

[6], [22] and[23] developed the multi-D ISS FSME algorithm from the Inverse Scattering Subseries for removing free-surface multiples (See Equations 1 and 2). For a 2D subsurface and towed streamer data, the ISS FSME algorithm for data without free surface multiples is

In Equations 1 and 2 D'1 is the input deghosted reflection data containing primaries, and free surface and internal multiples and D' is the output with primaries and only internal multiples. The quantities A(ω),ϵg and ϵs, in Equation 2 are the source signature, receiver depth and source depth, respectively; kg, ks are W Q wavenumbers of receivers and sources, ω is the temporal frequency and the obliquity 200 factor q is q=√(ω2/co 2) -k2).

The first term in this algorithm is the input data, D'1 in a 2D case, which is the Fourier transform of the deghosted prestack reflection data, with the direct wave and its ghost removed. The subsequent prediction terms, represented by D'2 , D'3 ,..., provide predictions of free-surface multiples of different orders. Specifically, each term in D'n (where n=2,3,4,...) performs two functions: (1) it predicts the nth order free-surface multiple and (2) it alters all higher order free-surface multiples to be prepared to be removed by higher-order D'j terms, where n=n+1,n+2,... The order of a free surface multiple is defined by the number of times the multiple has a downward reflection at the free surface.

The sum of these predictions (D 2'+D 3'+...+D n+1') provide free-surface-multiple predictions with accurate time and accurate amplitude (in opposite polarity) for free-surface multiples up to n-th order [6], [24].

The data, D' with free-surface multiples eliminated is obtained by Equation 1.

For the SRME approach, the method begins with the removal of: (1) the direct wave and (2) source and receiver ghosts, and then the SRME approximate free surface multiple prediction [25]- [27], M, can be expressed as follows

The input, P, is the prestack data for one temporal frequency, ω and where x g and x s are the location of the receivers and sources, respectively. Notice that, the input P for SRME and the input D 1 ' for ISS FSME are the same and both assume the removal of the direct wave and the source and receiver ghosts.

The output M in Equation 3 is the time and amplitude approximate free surface multiple prediction, provided within the approximations and assumptions in the SRME derivation and algorithm. The difference between the SRME approximate free surface multiple prediction, Equation 3 and the ISS exact free surface multiple predictor Equations 1 and 2 resides in the obliquity factor q, a function of frequency and wavenumber, and hence it causes an error in the SRME amplitude and phase prediction of the free surface multiple at all offsets. This SRME approximate free surface multiple is then energy minimization adaptive subtracted from the data in an attempt to match the amplitude and phase of the free surface multiple and thereby obtain the data without free-surface multiples. That lack of an accurate time and amplitude prediction in SRME is explicitly recognized by the energy minimization adaptive subtraction as a necessary and intrinsic part of the algorithm.

A quantitative comparison of SRME and the ISS Free Surface Multiple Elimination (FSME) algorithm can be found in [19], [20], see Figure 2 and 3. That analysis helps to define when SRME and ISS free surface elimination are the appropriate and indicated choice within the free-surface multiple removal seismic toolbox.

Figure 2 Model used to generated synthetic data. Two primaries (Blue) and one free-surface multiple (Red) are generated. 

Figure 3 (a) Input data generated using model shown in Figure 2. Two primaries are pointed by the blue arrows, one free-surface multiple is pointed by the red arrow, (b) ISS free-surface multiple prediction (c) SRME free-surface multiple prediction (d) Actual primaries in the data (e) Result after ISS FSME (f) Result after SRME + Adaptive subtraction. The free-surface multiple is interfering with the recorded primary. The SRME + Adaptive damages the primary that interferes with the free surface multiple. The ISS free-surface algorithm effectively removes the free surface multiple without damaging the primary. 

The result shows SRME + adaptive subtraction can be an effective and appropriate choice to remove isolated free-surface multiples, but can be injurious when applied to remove a FS multiple that is proximal or interfering with other events. The ISS FSME is effective and the appropriate choice whether or not the FS multiple is isolated or interfering with other events. The ISS FSME can surgically remove free-surface multiple that interfere with primaries or other events, and without damaging primaries.

There are many off-shore and on-shore plays where it is not clear, a priori, whether there are (or are not) free surface multiples that interfere with other events. The ISS free surface multiple eliminator is always a prudent choice.

THE CURRENT HIGH WATER MARK OF FREE SURFACE AND INTERNAL MULTIPLE REMOVAL

The ISS free surface multiple elimination algorithm (see e.g., [6], [22], [28]) predicts both the exact time and amplitude of all orders of free surface multiples at all offsets. It is effective with either isolated or interfering free surface multiples.

The ISS internal multiple attenuation algorithm attenuates internal multiples predicting the exact time and approximate amplitude of internal multiples is the only internal multiple algorithm that requires absolutely no subsurface information and often will be applied along with an energy minimization adaptive subtraction to remove an internal multiple that is not proximal to other events. To remove an internal multiple that is proximal to or interferes with other events (and therefore cannot rely on energy minimization, since the energy minimization criteria itself can fail under those circumstances), we need a more capable prediction, to surgically remove the multiple without damaging a nearby or interfering event. ISS internal multiple elimination had its origins in [29], discussion in [30], and an initial algorithm development in [31] and a fuller development and multidimensional algorithm in [32]- [34]

The ISS internal multiple attenuation algorithm is model type independent. That is, one absolutely unchanged algorithm (and with no change whatsoever in the computer code) predicts the precise time and approximate amplitude of all internal multiples independent of whether the subsurface is acoustic, elastic, anisotropic or anelastic. Filling the gap between the SS internal multiple attenuation and the elimination of the internal multiples is currently assuming an acoustic medium However a major contributor to the ISS internal multiple eliminator is the ISS internal multiple attenuator and the latter is model type independent. As we noted the gap filling part of the latter ISS internal multiple elimination algorithm [35] is based on an acoustic medium, and the effectiveness under different circumstances for acoustic, elastic and an-elastic media is evaluated in [34], [36] and [37]

THE ISS INTERNAL-MULTIPLE ATTENUATION ALGORITHM AND THE ID ISS INTERNAL MULTIPLE ELIMINATION ALGORITHM

The ISS internal-multiple attenuation algorithm was pioneered and developed in [21] and [22], The ID normal incidence version of the algorithm is presented in Equation 4 below (The 2D version is given in [6], [21] and [22] and the 3D version is a straightforward extension).

In Equation 4b 1 (z) is the constant velocity Stolt migration of the reflection data resulting from a ID normal incidence spike plane wave. ϵ1 and ϵ2 are two small positive numbers introduced to strictly maintain a lower-higher-lower relationship between the three water speed images and to avoid two water speed images at the same depth. b3 (k) is the predicted first order internal multiples in the vertical wavenumber domain. This algorithm can predict the correct time and approximate amplitude of all first-order internal multiples at once without any subsurface information.

Innanen and colleagues (e.g., [38]) have investigated the sensitivity of the choice of epsilon in Equation 4 in terms of the required lower higher lower pseudo depth relation that the subevents need to satisfy in order to combine to predict an internal multiple. They have suggested and have exemplified a non-stationary epsilon strategy, that navigate the issues between a too small (predictor becomes a "primary-like" artifact) and too large (missing predicting some internal multiples) epsilon value, and they propose that a priori geologic information can assist. Our view is that the very meaning of a primary and an internal multiple is a bandwidth dependent concept, and hence, e.g., there are events that we consider to be primaries that in fact under broader bandwidth would be a superposition of subresolution internal multiples. The ISS internal multiple attenuation and elimination algorithms assume definitions of primaries and internal multiples that are defined and have meaning within the bandwidth of the recorded data set.

The ISS internal-multiple attenuation algorithm automatically uses three primaries in the data to predict a first-order internal multiple. (Note that this algorithm is model type independent and it operates by taking into account all possible combinations of primaries that can be combined in a lower-higher-lower sense to predict internal multiples.).

The following equations are the 1D pre-stack ISS internal multiple elimination algorithm, see [32] and [33] for details and 1D examples].

where

The data without internal multiples, D is provided by Equations 5-7 and b 1 +b E =-2iqD" (see the discussion later in this paper and Table 1 for details and a broader perspective on the processing chain and steps in multiple removal).

Table 1 Steps that towed streamer data goes through (in 2D, it is similar in 3D) in the removal of free surface and internal multiples. D, D' and D" are the recorded data, data with free surface multiples removed and data without free surface and internal multiples, respectively. 

THE FIRST INVERSE-SCATTERING-SERIES INTERNAL MULTIPLE ELIMINATION METHOD FOR A MULTIDIMENSIONAL SUBSURFACE

The multi-dimensional ISS internal multiple elimination algorithm, see e.g., [35], [39] is provided in Equations 8, 9 and 10. This elimination formula is for all first order internal multiples from all reflectors at once, and without subsurface information. A first order internal multiple has one downward reflection in its history.

Similar to the 1D ISS internal multiple elimination algorithm [40] it is useful to introduce two intermediate functions F(k1,k2,z) and g(k 1 ,k 2 ,z) as follows:

Once again b 1 +b E =-2iqD", where D" is the data without first order internal multiples. The generalization for eliminating higher order internal multiples follows from the corresponding higher order ISS internal multiple attenuation algorithm in [6] and [21].

We thought that it might be useful at this point to remind ourselves of the data processing steps that towed streamer data goes through from the time it is recorded to the removal of all multiples and consists only of primaries. See Table 1.

SYNTHETIC DATA EXAMPLE OF MULTI D INTERNAL MULTIPLE ELIMINATION

The left part of the Figure 4 shows the 2D model. The data is generated by a finite difference method. The acoustic model is designed so that the base salt primary is negatively interfering with a first order internal multiple whose downward reflection is at the water bottom.

Figure 4 The model and zero offset traces of data. The base salt is almost invisible because the primary generated by the base salt is negatively interfering with an internal multiple. 

The base salt is almost invisible because the primary from base salt is negatively interfering with the water bottom downward reflected internal multiple. The left hand side of Figure 5 shows that the ISS internal multiple attenuation plus energy minimization adaptive subtraction does not recover the base salt image. The right hand side of Figure 5 shows the result after ISS internal-multiple elimination, The base salt is recovered. It demonstrates that the elimination algorithm can predict both the correct time and amplitude and can eliminate internal multiples without damaging an interfering or proximal primary.

Figure 5 The model and zero offset traces of data. The base salt is almost invisible because the primary generated by the base salt is negatively interfering with an internal multiple. 

CONCLUSION ON INTERNAL MULTIPLE REMOVAL

The ISS multi-dimensional internal-multiple-elimination algorithm that removes internal multiples is one part of a three-pronged strategy that is a response to current seismic processing and interpretation challenge that occurs when primaries and internal multiples are proximal to and/or interfere with each other. That can frequently occur in on-shore and off-shore plays.

The other two parts of the three part strategy involve: (1) preprocessing for on-shore plays and (2) developing a new adaptive criteria for the internal multiple elimination algorithm. Recent progress in preprocessing non-horizontal undulating off-shore cables and on-shore acquisition can be found in the following references: [36], [41]- [50]. An example of a new adaptive criteria for the case of the ISS free surface elimination is provided in [1] We are pursuing a similar criteria (that derives as a property of the SS predictor) for internal multiple elimination.

The ISS internal multiple elimination is a direct solution for the removal of multiples within the assumed physics and acquisition requirements. The adaptive step is indirect and is designed for addressing the parts of reality and e.g. linear wave propagation assumption and acquisition limitations that are outside and beyond our assumed physics.

This new internal multiple elimination algorithm addresses the prediction shortcoming of the current most capable internal-multiple-removal method (ISS internal-multiple-attenuation algorithm plus adaptive subtraction). Meanwhile, this elimination algorithm retains the strengths of the ISS internal-multiple-attenuation algorithm that can predict all internal multiples at once and requiring no subsurface information. This ISS internal-multiple-elimination algorithm is more effective and more compute-intensive than the current industry-standard most capable internal-multiple-removal method, i.e., the ISS internal multiple attenuator. Within the three-pronged strategy, our plans include developing an alternative adaptive-subtraction criteria for internal-multiple elimination derived from, and always aligned with the ISS elimination algorithm. That would be analogous to the new adaptive criteria for free-surface-multiple removal proposed by [1], as a replacement for internal multiple elimination for the energy-minimization criteria for adaptive subtraction. We provide this new multi-dimensional internal-multiple-elimination method as a new internal-multiple-removal capability in the multiple-removal toolbox that can remove internal multiples that interfere with primaries without damaging the primary and without subsurface information.

Various strategies to provide an effective eliminator in anelastic media include: (1) developing a modeltype independent ISS internal multiple eliminator and (2) employ an ISS subseries that inputs data that has experienced absorption and dispersion and outputs the data as though it had only experienced an elastic subsurface, without knowing or needing to know or to estimate or determine the absorption and dispersion mechanism [13]

The capability and potential of the recently developed inverse scattering subseries that performs Q compensation without knowing, estimating or determining Q [13] and without low frequency or zero frequency data is illustrated in Figures 6, 7 and 8. The same ISS subseries that performs Q compensation without needing to know or determine Q, can be easily adjusted to provide a subsurface map of Q. That advance has implications in many seismic and non-seismic applications. Among applications are the ability to avoid the need for low and zero frequency data in all amplitude analysis methods, including all indirect model-matching and updating methods. There are very significant applications of this new Q compensation method to electromagnetic prospecting and data analysis. These two references [51] and [52] show how the role of Q in seismic wave propagation corresponds to conductivity in electromagnetic propagation. The latter represents the practical potential of producing a subsurface conductivity map, and therefore a way to distinguish water from oil in the surface.

Figure 6 Two-reflector model for Q compensation without Q. [13]

Figure 7 Left: Data generated by the model with Q. Middle: The data (with Q) after ISS Q compensation without Q Right: Data generated by the same model but without Q. [13]

Figure 8 One trace comparison magnifying the event in the previous slide between 3.2 s -3.5 s. Red line: Data with Q. Green line: Data with Q after Q compensation. Blue line: Data without Q [34]

COMMENTS ON DIRECT INVERSION AND INDIRECT INVERSION (MODEL MATCHING):

The only direct inverse methods for parameter estimation the parameter estimation subseries of the inverse scattering series, pioneered by [7], [53], [54], (see [5]) specify the data and algorithms. The required data is a complete set of shot records with multi-component primaries. In these references it is shown that the elastic inhomogeneous isotropic elastic wave equation becomes a matrix operator identity in terms of a data matrix (in 2D),

The inverse solution for V is a generalized Geometric series in the data matrix D, where the V PP , V PS , V SP and V SS contain the sought after mechanical properties of the subsurface.

The forward series for D, in terms of V can be solved for one component of D, say, for example D PP. If one were to consider solving the latter forward problem for D PP in "an inverse sense", one would incorrectly deduce that D PP is an adequate data for an inverse solution. That thinking would violate the basic operator identity relationship (the Lippmann Schwinger equation) that solves for V (or any one component of V) in terms of the data matrix D. Please also see Appendix A for the detail relationship between D and V and why the entire data matrix D is required for a direct elastic parameter inversion, Equations A-28-A-32 and in addition a simple analytic example demonstrating that solving a forward problem in an inverse sense is not the same as solving the inverse problem directly.

In contrast, to the specificity in terms of data and inverse algorithms provided by direct solutions, with model matching methods, e.g., n the recent model matching approach, FWI there is no guide, no underlying theory or conceptual platform for what data is adequate, in principle, one trace, many traces, multi-component traces, and horizontal and vertical derivatives of displacement and pressure, and stress measurements and gravity data in fact, absolutely any data can be chosen to be model matched, including only one trace, or traces with only multiples. It seems reasonable that adding more data and data types would provide more constraints to search algorithms that might benefit and assist the parameter identification objective and reduce ill posedness however while including free surface multiples with primaries is often viewed as helpful, with added data constraints for the modeling to match, the addition of internal multiples seems in practice to be "too fuir model matching with too many complicated constraints to satisfy Under most circumstances internal multiples are attenuated before a FWI model matching begins. It seems that model matching with only primaries is viewed as not "full" enough, with primaries and free surface multiples that feels just right and perfectly full, and with the addition of internal multiples, apparently a little "too full", We are back to the lack of an underlying theory and framework, Why would a so called "full" wavefield inversion need to exclude internal multiples?

A great pedagogical advantage of indirect model methods is they are conceptually simple and readily understandable. Take a modeled trace and an actual trace and try to adjust the model parameters so the two traces match. Not hard to follow and understand. Indirect model matching methods also require a great deal of computer power and investment for search algorithms, and that expenditure "must" be based on a firm scientific foundation. In contrast, direct methods require an investment in understanding the physics and math-physics behind forward and inverse scattering. The first mention and derivation of the Lippmann Schwinger equation often has many in the audience dreaming of and pining for the simplicity of model matching concepts. Direct methods are often applied without any understanding of derivations behind the equations being coded and the services provided. For example, every major service company today and many oil companies provide a service based on the ISS internal multiple attenuator. It is extremely rare to find an individual who understands the underlying math-physics message and promise of the ISS series and how the ISS derives the isolated task subseries that attenuates internal multiples.

Indirect methods have a useful role and place within the seismic toolbox, and as with all seismic methods (including all migration, Green's theorem and ISS methods), we welcome and encourage a balanced view of the benefits, shortcomings and open issues.

CONCLUSION

Multiple removal and using multiples have one single exact goal: imaging primaries, recorded and unrecorded primaries. To be effective at reaching that objective recorded and unrecorded multiples must be removed. Since recorded primaries have the greatest potential (via Stolt CIII migration and migration-inversion and ISS depth imaging and inversion) for delivering structure and amplitude analysis, the removal of recorded multiples as a concomitant high priority and interest.

The confusion over "using" multiples is not a harmless misunderstanding without consequences because if multiples were in fact the new signal and the equivalent of primaries then we should no longer remove multiples, no more than we remove primaries that is the danger that derives from a misinformed premise and conclusion in thinking that removing and using multiples are adversarial.

In the history of useful methods and contributions that seek to accommodate limited data acquisition, like DMO, and 2D and 2.5D processing with asymptotic techniques in the cross line direction, eventually data acquisition advances to provide the data necessary to reach processing and interpretation goals and methods that seek to accommodate limited data become less interesting and less relevant.

Multiple removal is a permanent issue, whereas multiple usage is transient, and the latter will eventually be replaced by a more complete recording of primaries. In the interim, advances in both removing and using multiples are welcome and needed.

ACKNOWLEDGEMENTS

The M-OSRP sponsors are gratefully appreciated for their support. The author is deeply grateful to Ecopetrol for the invitation to write this overview paper. Jim Mayhan is thanked for assisting in the preparation of this paper.

REFERENCES

[1] Weglein, A. B., Short note: An alternative adaptive subtraction criteria (to energy minimization) for free surface multiple removal, 2012, M-OSRP 2011-2012 Annual Report, 375. [ Links ]

[2] Weglein, A. B., A direct inverse method for subsurface properties: The conceptual and practical benefit and added value in comparison with all current indirect methods, for example, amplitude-variation-with-offset and full-waveform inversion, 2017, Interpretation, 5, SL89-SL105. [ Links ]

[3] Weglein, A. B., A timely and necessary antidote to indirect methods and so-called P-wave FWI, 2013, The Leading Edge, 32, 1192-1204. [ Links ]

[4] Stolt, R. H., and A. B. Weglein, (2012). Seismic imaging and inversion: Application of linear inverse theory, Cambridge University Press. [ Links ]

[5] Weglein, A. , J. Mayhan, Y. Zou, Q. Fu, F. Liu, J. Wu, C. Ma, X. Lin, and R. Stolt, The first migration method that is equally effective for all acquired frequencies for imaging and inverting at the target and reservoir, 2016, 86 th Annual International Meeting, SEG, Expanded Abstracts, 4266-4272. [ Links ]

[6] Weglein, A. B., F. V. Araújo, P. M. Carvalho, R. H. Stolt, K. H. Matson, R. T. Coates, D. Corrigan, D. J. Foster, S. A. Shaw, and H. Zhang, Inverse scattering series and seismic exploration, 2003, Inverse Problems, 19, R27-R83. [ Links ]

[7] Zhang, H., (2006). Direct non-linear acoustic and elastic inversion: Towards fundamentally new comprehensive and realistic target identification, PhD thesis, University of Houston, Houston, USA. [ Links ]

[8] Yang, J., (2014). Extending the Inverse Scattering Series free-surface multiple elimination and internal multiple attenuation algorithms by incorporating the source wavelet and radiation pattern: Examining and evaluating the benefit and added-value, PhD thesis, University of Houston, Houston, USA. [ Links ]

[9] Yang, J., and A. B. Weglein, Incorporating the source wavelet and radiation pattern into the ISS internal multiple attenuation algorithm: theory and examples, 2014, M-OSRP 2013-2014 Annual Report, 63-78. [ Links ]

[10] Zou, Y., Q. Fu, and A. B. Weglein, A wedge resolution comparison between RTM and the first migration method that is equally effective at all frequencies at the target: tests and analysis with both conventional and broadband data, 2017a, 87 th Annual International Meeting SEG, Expanded Abstracts, 4468-4472. [ Links ]

[11] Zou, Y., C. Ma, and A. B. Weglein. (2017b). A multidimensional method that eliminates internal multiples: a new tool box option for removing multiples that interfere with primaries, without damaging the primary, and without any knowledge of subsurface properties, presented at 87 th Annual International Meeting, SEG, Expanded Abstracts. (Submitted). [ Links ]

[12] Weglein, A. B., F. Liu, X. Li, P. Terenghi, E. Kragh, J.D. Mayhan, Z. Wang, J. Mispel, L. Amundsen, H. Liang, L. Tang, and S.-Y. Hsu, Inverse scattering series direct depth imaging without the velocity model: First field data examples, 2012, Journal of Seismic Exploration, 21, 1-28. [ Links ]

[13] Zou, Y., and A. Weglein, An inverse scattering subseries that performs Q compensation without knowing, estimating or determining Q, and without low and zero frequency data, 2018, Journal of Seismic Exploration . (Submitted May 2018). [ Links ]

[14] Whitmore, N. D., A. Valenciano, S. Lu, and Chemingui, N. (2011a). Imaging of Primaries and Multiples with Image Space Surface Related Multiple Elimination, presented at the 73rd EAGE Conference & Exhibition, Extended Abstracts. [ Links ]

[15] Whitmore, N. D., A. A. Valenciano, S. Lu, and Chemingui, N. (2011b). Deep water prestack imaging with primaries and multiples, presented at the Twelfth International Congress of the Brazilian Geophysical Society, Sociedade Brasileira de Geofísica. [ Links ]

[16] Liu, Y., X. Chang, D. Jin, R. He, H. Sun, and Y. Zheng, Reverse time migration of multiples for subsalt imaging, 2011, Geophysics, 76, WB209-WB216. [ Links ]

[17] Lu, S., N. D. Whitmore, A. A. Valenciano, and Chemingui, N.(2011). Imaging of primaries and multiples with 3D SEAM synthetic: 81 st Annual International Meeting, SEG, Expanded Abstracts, 3217-3221. [ Links ]

[18] Weglein, A. B., Multiples: Signal or noise?, 2016, Geophysics , 81, V283-V302. [ Links ]

[19] Ma, C., Q. Fu, and A. B. Weglein, Analysis, (2018a). testing and comparison of the inverse scattering series (ISS) free-surface multiple-elimination (FSME) algorithm, and the industry standard SRME plus energy minimization adaptive subtraction, presented at the M-OSRP 2017-2018 Annual Report. [ Links ]

[20] Ma, C., Q. Fu, and A. B. Weglein. (2018b). Analysis, testing and comparison of the inverse scattering series (ISS) free-surface multiple-elimination (FSME) algorithm, and the industry-standard SRME plus energy minimization adaptive subtraction, presented at the 88th Annual International Meeting, SEG, Expanded Abstracts. (Submitted). [ Links ]

[21] Araújo, F. V., A. B. Weglein, P. M. Carvalho, and R. H. Stolt, 1994, Inverse scattering series for multiple attenuation: An example with surface and internal multiples: 64th Annual International Meeting, SEG, Expanded Abstracts, 1039-1041. [ Links ]

[22] Weglein, A. B., F. A. Gasparotto, P. M. Carvalho, and R. H. Stolt, An inverse-scattering series method for attenuating multiples in seismic reflection data, 1997, Geophysics , 62, 1975-1989. [ Links ]

[23] Carvalho, P. M., and A. B. Weglein, Examples of a nonlinear inversion method based on the T matrix of scattering theory: Application to multiple suppression, 1991, 61 st Annual International Meeting, SEG, Expanded Abstracts, 1319-1322. [ Links ]

[24] Ma, C., and A. Weglein, Examining the interdependence and cooperation of the terms in the distinct inverse-scattering subseries for free-surface multiple and internal multiple removal, 2016, 86 th Annual International Meeting, SEG, Expanded Abstracts, 4561-4565. [ Links ]

[25] Berkhout, A. J., 1985, Seismic migration: Theoretical aspects, third edition, Elsevier Publishing Co. [ Links ]

[26] Verschuur, D. J., (1991). Surface-related multiple elimination, an inversion approach, PhD thesis, Delft Univ. of Technology, Netherlands. [ Links ]

[27] Verschuur, D., A. J. Berkhout, and C. P. A. Wapenaar, Adaptive surface-related multiple elimination, 1992, Geophysics , 57, 1166-1177. [ Links ]

[28] Carvalho, P. M., A. B. Weglein, and R. H. Stolt, Nonlinear inverse scattering for multiple suppression: Application to real data. Part I, 1992, 62 nd Annual International Meeting, SEG, Expanded Abstracts, 1093-1095. [ Links ]

[29] Weglein, A. B., and K. Matson, Inverse-scattering interval multiple attenuation: An analytic example and subevent interpretation, in Mathematical methods in geophysical imaging, 1998, SPIE, 1008-1017. [ Links ]

[30] Ramírez, A. C., and A. B. Weglein, An inverse scattering internal multiple elimination method: Beyond attenuation, a new algorithm and initial tests, 2005, 75th Annual International Meeting, SEG, Expanded Abstracts, 2115-2118. [ Links ]

[31] Herrera, W., and A. B. Weglein, Eliminating first-order internal multiples with downward reflection at the shallowest interface: Theory and initial examples, 2013, 83 rd Annual International Meeting, SEG, Expanded Abstracts, 4131-4135. [ Links ]

[32] Zou, Y., (2017). The first multi-dimensional inverse-scattering-series internal-multiple-elimination method: A new toolbox option for removing internal multiples that interfere with a primary, without damaging the primary, and without any knowledge of subsurface properties II. Tests and analysis for resolution comparisons between reverse time migration (RTM) and the first migration method that is equally effective at all frequencies at the target, PhD thesis, University of Houston, Houston, USA. [ Links ]

[33] Zou, Y., C. Ma, and A. Weglein, The first inverse-scattering-series internal multiple elimination method for a multidimensional subsurface, 2016, 86th Annual International Meeting, SEG, Expanded Abstracts, 4550-4554. [ Links ]

[34] Zou, Y., and A. Weglein, An inverse scattering subseries that performs Q compensation without knowing, estimating or determining Q, and without low and zero frequency data, 2018, Journal of Seismic Exploration . (Submitted May 2018). [ Links ]

[35] Zou, Y., C. Ma, and A. B. Weglein, (2018b). A multidimensional method that eliminates internal multiples: a new tool box option for removing multiples that interfere with primaries, without damaging the primary, and without any knowledge of subsurface properties, presented at the 88 th Annual International Meeting, SEG, Expanded Abstracts. (Submitted). [ Links ]

[36] Wu, J., and A. B. Weglein, A new method for deghosting data collected on a depth-variable acquisition surface by combining Green's theorem wave separation followed by a Stolt extended Claerbout III wave prediction for oneway propagating waves, 2017, 87th Annual International Meeting, SEG, Expanded Abstracts, 4859-4864. [ Links ]

[37] Fu, Q., Y. Zou, J. Wu, and A. B. Weglein, Analysis of the inverse scattering series (ISS) internal multiple attenuation and elimination algorithms as effective tool box choices for absorptive and dispersive media with interfering events, presented at the 88 th Annual International Meeting, SEG, Expanded Abstracts, 2018. (Submitted). [ Links ]

[38] Innanen, K. A., 2017, Time- and offset-domain internal multiple prediction with nonstationary parameters: Geophysics , 82, V105-V116. [ Links ]

[39] Zou, Y., C. Ma, and A. Weglein, 2018a, The first multidimensional inverse-scattering-series internal-multiple-elimination method: a new toolbox option for removing internal multiples that interfere with a primary, without damaging the primary, and without any knowledge of subsurface properties: Geophysics . (In preparation). [ Links ]

[40] Zou, Y., and A. B. Weglein, 2014, An internal-multiple elimination algorithm for all reflectors for 1D earth Part I: Strengths and limitations: Journal of Seismic Exploration , 23, 393-404. [ Links ]

[41] Weglein, A. B., J. D. Mayhan, L. Amundsen, and H. Liang, Green's theorem de-ghosting algorithms in the k, io (e.g., P-Vz de-ghosting) as a special case of x, io algorithms (based on Green's theorem) with: (1) significant practical advantages and disadvantages of algorithms in each domain, and (2) a new message, implication and opportunity for marine towed streamer, ocean bottom and on-shore acquisition and applications, 2013a, M-OSRP 2012-2013 Annual Report, 7-19. [ Links ]

[42] Weglein, A. B., J. D. Mayhan, L. Amundsen, H. Liang, J. Wu, L. Tang, Y. Luo, and Q. Fu, Green's theorem de-ghosting algorithms in k, ω (e.g., P-Vz deghosting) as a special case of x, ω algorithms (based on Green's theorem) with: (1) significant practical advantages and disadvantages of algorithms in each domain, and (2) a new message, implication and opportunity for marine towed streamer, ocean bottom and on-shore acquisition and applications, 2013b, Journal of Seismic Exploration , 22, 389-412. [ Links ]

[43] Weglein, A. B., and B. G. Secrest, Wavelet estimation for a multidimensional acoustic or elastic earth, 1990, Geophysics , 55, 902-913. [ Links ]

[44] Osen, A., L. Amundsen, and A. Reitan, Removal of water-layer multiples from multicomponent sea-bottom data: Geophysics , 64, 838-851. [ Links ]

[45] Osen, A., B. G. Secrest, L. Amundsen, and A. Reitan, 1998, Wavelet estimation from marine pressure measurements, 1999, Geophysics , 63, 2108-2119. [ Links ]

[46] Tan, T. H., Wavelet spectrum estimation, 1999, Geophysics , 64, 1836-1846. [ Links ]

[47] Zhang, J., (2007). Wave theory based data preparation for inverse scattering multiple removal, depth imaging and parameter estimation: analysis and numerical tests of Green's theorem deghosting theory, PhD thesis, University of Houston, Houston, USA. [ Links ]

[48] Mayhan, J. D., and A. B. Weglein, First application of Green's theorem-derived source and receiver deghosting on deep-water Gulf of Mexico synthetic (SEAM) and field data, 2013, Geophysics , 78, WA77-WA89. [ Links ]

[49] Shen, Y., Z. Zhang, and A. B. Weglein, Initial tests of the new two-step strategy method for deghosting on a depth variable acquisition surface by combining Greens theorem and a Stolt extended Claerbout III wave prediction for one way propagating waves, 2016, M-OSRP 2015-2016 Annual Report, 121-142. [ Links ]

[50] Amundsen, L., A. Reitan, A. B. Weglein, and BØrn Ursin, On seismic deghosting using Green's theorem, 2016, Geophysics , 81, V317-V325. [ Links ]

[51] Lax, M., and D. F. Nelson, Maxwell equations in material form, 1976, Physical Review B, 13, 1777-1784. [ Links ]

[52] Wolski, A., PHYS370 Advanced Electromagnetism, 2014, [Online]. Available at Available at http://pcwww.liv.ac.uk/~awolski/main_teaching_Liverpool_PHYS370.htm . (Accessed 28 May 2018; see Part 3: Electromagnetic waves in conducting media). [ Links ]

[53] Li, X., I. (2011). Multi-parameter depth imaging using the inverse scattering series; II.- Multi-component direct non-linear inversion for elastic earth properties using the inverse scattering series, PhD thesis, University of Houston, Houston, USA. [ Links ]

[54] Liang, H., (2013). Addressing several key outstanding issues and extending the capability of the inverse scattering subseries for internal multiple attenuation, depth imaging, and parameter estimation, PhD thesis, University of Houston, Houston, USA. [ Links ]

[55] Taylor, J. R., (1972). Scattering theory: The quantum theory of nonrelativistic collisions, John Wiley & Sons, Inc. [ Links ]

[56] Claerbout, J. F., Toward a unified theory of reflector mapping,1971, Geophysics , 36, 467-481. [ Links ]

[57] Stolt, R. H., Migration by Fourier transform, 1978, Geophysics , 43, 23-48. [ Links ]

[58] Schneider, W. A., Integral formulation for migration in two and three dimensions, 1978, Geophysics , 43, 49-76. [ Links ]

[59] Ware, J. A., and K. Aki, Continuous and discrete inverse-scattering problems in a stratified elastic medium. I. Plane waves at normal incidence, 1969, The Journal of the Acoustical Society of America, 45, 911-921. [ Links ]

[60] Beylkin, G., and R. Burridge, Linearized inverse scattering problem of acoustics and elasticity, 1990, Wave Motion, 12, 15-52. [ Links ]

[61] Boyse, W. E., and J. B. Keller, Inverse elastic scattering in three dimensions, 1986, The Journal of the Acoustical Society of America, 79, 215-218. [ Links ]

[62] Burridge, R., M. de Hoop, D. Miller, and C. Spencer, Multiparameter inversion in anisotropic elastic media, 1998, Geophysical Journal International, 134, 757-777. [ Links ]

[63] Castagna, J. P., and S. W. Smith, Comparison of AVO indicators: A modeling study, 1994, Geophysics , 59, 1849-1855. [ Links ]

[64] Clayton, R. W., and R. H. Stolt, A Born-WKBJ inversion method for acoustic reflection data, 1981, Geophysics , 46, 1559-1567. [ Links ]

[65] Foster, D. J., R. G. Keys, and F. D. Lane, Interpretation of AVO anomalies, 2010, Geophysics , 75, 75A3-75A13. [ Links ]

[66] Goodway, B., The magic of Lamé, 2010, The Leading Edge, 29, 14-32. [ Links ]

[67] Goodway, B., T. Chen, and J. Downton, Improved AVO fluid detection and lithology discrimination using Lamé petrophysical parameters; " λp", " μp", and "λ/μ fluid stack", from P and S inversions, 1997, 67 th Annual International Meeting, SEG, Expanded Abstracts, 183-186. [ Links ]

[68] Shuey, R. T., A simplication of the Zoeppritz equations, 1985, Geophysics , 50, 609-614. [ Links ]

[69] Smith, G. C., and M. Gidlow, A comparison of the fluid factor with i and u in AVO analysis, 2000, 70 th Annual International Meeting, SEG, Expanded Abstracts, 122-125. [ Links ]

[70] Stolt, R. H., Seismic inversion revisited: Geophysical Inversion, Society for Industrial and Applied Mathematics, 3-19. (Proceedings of the Geophysical Inversion Workshop, September 27-29, 1989, Houston, Texas). [ Links ]

[71] Stolt, R. H., and A. B. Weglein, Migration and inversion of seismic data, 1985, Geophysics , 50, 2458-2472. [ Links ]

[72] Weglein, A. B., H. Zhang, A. C. Ramírez, F. Liu, and J. E. M. Lira, Clarifying the underlying and fundamental meaning of the approximate linear inversion of seismic data, 2009, Geophysics , 74, WCD1-WCD13. [ Links ]

[73] Stolt, R. H., and A. B. Weglein, (2012). Seismic imaging and inversion: Application of linear inverse theory , Cambridge University Press. [ Links ]

APPENDIX

WHY SOLVING A FORWARD PROBLEM IN AN INVERSE SENSE IS NOT THE SAME AS SOLVING AN INVERSE PROBLEM DIRECTLY

Scattering theory is a form of perturbation theory, starting with a medium described by L0 and a perturbation in that medium described by L.

and seek the relationship between L0-L = V and G-G0=I|JS.

The operator identity

[for a fixed source function] is the exact relationship between changes in a medium and changes in the wavefield. The operator identity Equation A-1 can be solved for G as

and expanded as

Equation A-3 is called the Born or Neumann series in scattering theory literature (see, e.g., [55]). Equation A-3 has the form of a generalized geometric series

where we identify a=G0 and r=VG0 in Equation A-3, and

where the portion of S that is linear, quadratic, ... in r is:

and the sum is

Solving Equation A-6 for r, in terms of S/a produces the inverse geometric series,

where rn is the portion of r that is nth order in S/a. When S is a geometric power series in r, then r is a geometric power series in S. The former is the forward series and the latter is the inverse series. That is exactly what the inverse scattering series represents, the inverse geometric series of the forward series Equation A-3. This is the simplest prototype of an inverse series for r, i.e., the inverse of the forward geometric series for S.

SOLVING A FORWARD PROBLEM IN AN INVERSE SENSE IS NOT THE SAME AS SOLVING AN INVERSE PROBLEM DIRECTLY

We will show that in general solving a forward problem in an inverse sense is not the same as solving an inverse problem directly. The exception is when the exact direct inverse is linear, as e.g. in the theory of wave equation migration (see, e.g. [4], [5], [56]- [58];). For wave equation migration, given a velocity model, the direct and exact structure map output relationship is a linear function of the input recorded reflection data. Hence, (the direct and exact linear inverse represented by) wave equation migration can be viewed as wave equation modeling (a forward problem) run backwards in either depth or time.

To help explain the latter statement, if we assume S=ar (that is, that there is an exact linear forward relationship between S and r) then r=S/a is solving the inverse problem directly. In that case, solving the forward problem in an inverse sense is the same as solving the inverse problem directly, that is, it provides a direct inverse solution.

However, if the forward exact relationship is non-linear, for example, a geometric polynomial or series (for |r|<1)

and solving the forward problem (Equation A-8) in an inverse sense for r will have n roots, r1,r2....,rn. As n →∞, the number of roots →∞. However, from the direct nonlinear forward problem S=ar/(1-r), we found the direct inverse solution r=S/(a+S) has one real root.

This discussion above provides an extremely simple, transparent and compelling illustration of how solving a forward problem in an inverse sense is not the same as solving the inverse problem directly when there is a non-linear forward and non-linear inverse problem. The difference between solving a forward problem in an inverse sense (for example using Equation A-3 to solve for V) and solving an inverse problem directly (for example, Equations A-10-A-12 is much more serious, substantive and practically significant the further we move away from a scalar single component acoustic framework. For example, it is hard to overstate the differences when examining the direct and indirect inversion of the elastic heterogeneous wave equation for earth mechanical properties, and the consequences for structural and amplitude analysis and interpretation. This is a central flaw in many inverse approaches, including AVO and FWI [3].

The inverse scattering series[6] corresponding to the forward series Equation A-3 and generalizing the scalar form Equation A-7

where Vn is the portion of V that is nth order in measured data, D. The expansion of V in Equation A-9, in terms of G0 and D=(G-G0)ms, the inverse scattering series [6] can be obtained as

To illustrate how to solve Equations A-10, A-11, A-12, ... for V1, V2, V3, ... consider the marine case with L0 corresponding to a homogeneous reference medium of water. G0 is the Green's function for propagation in water. D is the data measured for example, with towed streamer acquisition, G is the total field the hydrophone records on the measurement surface, and G0 is the field the reference wave (due to L0) would record at the receiver. V then represents the difference between earth properties L and water properties L0. The solution for V is found from Equation A-9. Substituting Equation A-9 into the forward series Equation A-3, then evaluating Equation A-3 on the measurement surface and setting terms that are equal order in the data to be equal we find equations on each side of the Equation A-10, A-11, A-12... Solving Equation A-10 for V1 involves the data D and G0 (water speed propagator) and solving for V1 is analytic, and corresponds to a prestack water-speed Stolt FK migration of the data D.

Hence, solving for V1 involves an analytic water speed FK migration of the data D. Solving for V2 from Equation A-11 involves the same water-speed analytic Stolt FK migration of -G0V1G0V1G0, a quantity that depends on V1 and G0, where V1 depends on data and water speed, and G0 is the water speed Green's function. Each term in the series produces Vn as an analytic Stolt FK migration of a new "effective data", where the effective data, the right-hand side of Equations A-10-A-12, are multiplicative combinations of factors that only depend on the data, D, and G0. Hence, every term in the ISS is directly computed in terms of data and water speed. That's the direct non-linear inverse solution.

There are closed form inverse solutions for a 1D earth and a normal incident plane wave (see, e.g., [59]) but the inverse scattering series is the only direct inverse method for a multi-dimensional subsurface.

The inverse scattering series provides a direct method for obtaining the subsurface properties contained within the differential operator L, by inverting the series order-by-order to solve for the perturbation operator V, using only the measured data D and a reference Green's function G0, for any assumed earth model type. Equations A-10-A-12 provide V in terms of V1, V2..., and each of the Vi is computable directly in terms of D and G0. There is one equation, Equation A-10, that exactly produces V1, and V1 is the exact portion of V that is linear in the measured data, D. The inverse operation to determine V1, V2, V3,..., is analytic, and never is updated with a bandlimited data, D. The band-limited nature of D never enters an updating process as occurs in iterative linear inversion, non-linear AVO and FWI.

Examples of the fallacy in thinking that solving a forward problem in an inverse sense (order by order or otherwise) is equivalent to a direct inverse solution are in several of our papers, for example, Equations 14-27 (from [7]) in the first reference are the direct elastic heterogeneous isotropic inverse, and LS requires a matrix solution, in terms of a set of multicomponent data, since the heterogeneous elastic LS equation is a matrix equation.

ONE CAN SOLVE THE FORWARD PROBLEM AS A FORWARD SERIES FOR PP DATA ALONE

(or for SS, SP... alone) each separately in terms of VPP, VPS, VSS... one could take the latter PP data and term by term "invert it" solving a forward problem in an inverse sense and erroneously conclude that the direct non-linear inverse can be solved in terms of only PP data, term by term order by order, and that is not a direct inverse solution; that PP "solution" violates the LS operator identity for the elastic wave equation and its matrix inverse solution, where all data components are necessary.

DIRECT INVERSE AND INDIRECT INVERSE

Since iterative linear inversion is the concept and thinking behind many inverse approaches we thought to make explicit the difference between that approach and a direct inverse method. The direct 2D elastic isotropic inverse solution described in Appendix A is not iterative linear inversion. Iterative linear inversion starts with Equation A-10. In that approach, we solve for V1 and then change the reference medium iteratively. The new differential operator L0' and the new reference medium G0' satisfy

In the indirect iterative linear approach, all steps basically relate to the linear relationship Equation A-10 with a new reference background medium, with differential operator L0' and a new reference Green's function G0' where in terms of the new updated reference, L0', Equation A-10 becomes

where V1' is the portion of V linear in data (G-G0')ms. We can continue to update L0' and G0', and hope that indirect procedure is solving for the perturbation operator V. In contrast, the direct inverse solution Equations A-9 and A-12 calls for a single unchanged reference medium, for computing V1,V2,... For a homogeneous reference medium, V1,V2,. are each obtained by a single unchanged analytic inverse. We remind ourselves that the inverse to find V1 from data, is the same exact unchanged analytic inverse operation to find V2,V3..., from Equations A-10, A-11,...,which is completely distinct and different from Equations A-13 and A-14 and higher iterates.

For ISS direct inversion, there are no numerical inverses, no generalized inverses, no inverses of matrices that are computed from and contain noisy band-limited data. The latter issue is terribly troublesome and difficult and a serious practical problem which simply does not exist or occur with direct ISS methods. The inverse of operators that contain and depend on band-limited noisy data is a central and intrinsic characteristic and practical pitfall of indirect methods, model matching, updating, iterative linear inverse approaches (e.g. AVO and FWI).

FURTHER SUBSTANTIVE DIFFERENCES BETWEEN ITERATIVE LINEAR MODEL MATCHING INVERSION AND DIRECT INVERSION FROM THE LIPPMANN-SCHWINGER EQUATION AND THE INVERSE SCATTERING SERIES

The difference between iterative linear and the direct inverse of Equation A-10 is much more substantive and serious than merely a different way to solve G0V1G0 = D [Equation A-10], for V1. If Equation A-10 is someone's entire basic theory, you can mistakenly think that

is sufficient to update

(generalizing Equations A-13 and A-14). Note that ^ indicates variables are transformed to PS space. This step loses contact with and violates the basic operator identity G = G0 + G0VG for the elastic wave equation. The fundamental identity G=G0+G0VG for the elastic wave equation is a non-linear multiplicative matrix relationship. For the forward and inverse series the input and output variables are matrices. The inverse solution for a change in an earth mechanical property has a nonlinear coupled dependence on \underline(all] the data components

in 2D and the P, SH, SV 3 x 3 generalization in 3D.

A unique expansion of VG0 in orders of measurement values of (G-G0) is

The scattering-theory equation allows that forward series form the opportunity to find a direct inverse solution. Substituting Equation A-17 into Equation A-3 and setting the terms of equal order in the data to be equal, we have D=G0V1G0, where the higher order terms are V2, V3...,as given in [6] page R33 Equation 7-14.

For the elastic equation, V is a matrix and the relationship between the data and V1 is

where V1, V2 are linear, quadratic contributions to V in terms of the data,

The changes in elastic properties and density are contained in

and that leads to direct and explicit solutions for the changes in mechanical properties in orders of the data,

where γ, μ and ρ are the bulk modulus, shear modulus and density, respectively.

The ability of the forward series to have a direct inverse series derives from (1) the identity among G, G0, V provided by the scattering equation and then (2) the recognition that the forward solution can be viewed as a geometric series for the data, D, in terms of VG0. The latter derives the direct inverse series for VG0 in terms of the data.

Viewing the forward problem and series as the Taylor series

in which the derivatives are Frechet derivatives, in terms of Δm does not offer a direct inverse series, and hence there is no choice but to solve the forward series in an inverse sense. It is that fact that results in all current AVO and FWI methods being modeling methods that are solved in an inverse sense. Among references that solve a forward problem in an inverse sense in P-wave AVO are [60]- [71]. The intervention of the explicit relationship among G, G0, and V (the scattering equation) in a Taylor series-like form produces a geometric series and a direct inverse solution.

The linear equations are:

where ay [1], aμ [1], and ap [1] are the linear estimates of the changes in bulk modulus, shear modulus, and density, respectively. kg is the Fourier conjugate to the receiver position xg and vg and ηg are the vertical wavenumbers for the P and S reference waves, respectively, where

and α0 and β0 are the P and S velocities in the reference medium, respectively. The direct quadratic non-linear equations are

Because V̂ 1 PP relates to D̂PP , V̂1 PS relates to D̂PS, and so on, the four components of the data will be coupled in the nonlinear elastic inversion. We cannot perform the direct nonlinear inversion without knowing all components of the data. Thus, the direct nonlinear solution determines the data needed for a direct inverse. That, in turn, defines what a linear estimate means. That is, a linear estimate of a parameter is an estimate of a parameter that is linear in data that can directly invert for that parameter. Since DPP, DPS, DSP, and DSS are needed to determine ay, aμ, and ap directly, a linear estimate for any one of these quantities requires simultaneously solving Equations A-24-A-27. See, e.g., [72] for further detail.

Those direct nonlinear formulas are like the direct solution for the quadratic equation mentioned above and solve directly and nonlinearly for changes in the velocities, α, β and the density p in a 1D elastic Earth. [73], present the linear equations for a 3D Earth that generalize Equations A-24-A-27. Those formulas prescribe precisely what data you need as input, and they dictate how to compute those sought-after mechanical properties, given the necessary data. There is no search or cost function, and the unambiguous and unequivocal data needed are full multicomponent data PP, PS, SP, and SS for all traces in each of the P and S shot records. The direct algorithm determines first the data needed and then the appropriate algorithms for using those data to directly compute the sought-after changes in the Earth's mechanical properties. Hence, any method that calls itself inversion (let alone full-wave inversion) for determining changes in elastic properties, and in particular the P-wave velocity a, and that inputs only P-data, is more off base, misguided, and lost than the methods that sought two or more functions of depth from a single trace. You can model-match P-data until the cows come home, and that takes a lot of computational effort and people with advanced degrees in math and physics computing Frechet derivatives, and requires sophisticated LP norm cost functions and local or global search engines, so it must be reasonable, scientific, and worthwhile. Why can not we use just PP-data to invert for changes in VP, VS, and density, because Zoeppritz says that we can model PP from those quantities, and because we have, using PP-data with angle variation, enough dimension? As stated above, data dimension is good, but not good enough for a direct inversion of those elastic properties.

Adopting Equations A-15 and A-16 as in AVO and FWI, there is a violation of the fundamental relationship between changes in a medium and changes in a wavefield, G=G0+G0VG, which is as serious as considering problems involving a right triangle and violating the Pythagorean theorem. That is, iteratively updating PP data with an elastic model violates the basic relationship between changes in a medium, V, and changes in the wavefield, G-G0, for the simplest elastic earth model.

This direct inverse method for parameter estimation provides a platform for amplitude analysis, and a solid framework and direct methodology for the goals and objectives of indirect methods like AVO and FWI. A direct method for the purposes of amplitude analysis provides a method that derives from, respects and honors the fundamental identity and relationship G = G0 + G0VG. Iteratively inverting multi-component data has the correct data but does not correspond to a direct inverse algorithm. To honor G=G0+G0VG, you need both the data and the algorithm that direct inverse prescribes. Not recognizing the message that an operator identity and the elastic wave equation unequivocally communicate is a fundamental and significant contribution to the gap in effectiveness in current AVO and FWI methods and application. This analysis generalizes to 3D with P, SH, and SV data ([4]).

Received: May 29, 2018; Revised: June 13, 2018; Accepted: August 23, 2018

Creative Commons License This is an open-access article distributed under the terms of the Creative Commons Attribution License