# Closed-loop multiple-scattering imaging with sparse seismic measurements

Closed-loop multiple-scattering imaging with sparse seismic measurements Summary In the theoretical situation of noise-free, complete data volumes (‘perfect data’), seismic data matrices are fully filled and multiple-scattering operators have the minimum-phase property. Perfect data allow direct inversion methods to be successful in removing surface and internal multiple scattering. Moreover, under these perfect data conditions direct source wavefields realize complete illumination (no irrecoverable shadow zones) and, therefore, primary reflections (first-order response) can provide us with the complete seismic image. However, in practice seismic measurements always contain noise and we never have complete data volumes at our disposal. We actually deal with sparse data matrices that cannot be directly inverted. The message of this paper is that in practice multiple scattering (including source ghosting) must not be removed but must be utilized. It is explained that in the real world we badly need multiple scattering to fill the illumination gaps in the subsurface. It is also explained that the proposed multiple-scattering imaging algorithm gives us the opportunity to decompose both the image and the wavefields into order-based constituents, making the multiple scattering extension easy to apply. Last but not least, the algorithm allows us to use the minimum-phase property to validate and improve images in an objective way. Seismic interferometry, Wave propagation, Wave scattering and diffraction INTRODUCTION Multiple scattering events (‘multiples’) have always played a vital role in seismic imaging. Traditionally, they are considered as unwanted events (being referred to as ‘shot-generated noise’) that interfere with—or may even overshadow—the primary reflections (the latter being considered as the ‘shot-generated signal’). A representative example of scientific arguments to consider multiple scattering as shot-generated noise can be found in Weglein (2014, 2015). In the past 50 yr enormous investments have been made to remove multiples, making the seismic reflection response linear in reflectivity. However, recently the seismic community started to realize that lack of understanding was the reason of this linearization strategy. There is increasing theoretical and practical evidence that multiples contain subsurface information that cannot be extracted from primary reflections, see for instance Whitmore et al. (2010) and Davidenko et al. (2015). Hence, multiples should not be removed as noise but they should be utilized as signal. The only noise left in seismic data is background noise (Berkhout 2014). The author considers utilization of multiples the next big step forward in the practice of seismic imaging (often referred to as ‘seismic migration’). We will see that this also applies to the utilization of source ghosts (the interaction of near-surface sources with the surface). Feedback model at the surface $${(z_{0}^{+})}$$ In marine data—and in hard-surface land data—multiples that are generated by the surface reflector are strong and must be addressed. Fig. 1(a) shows the feedback model of the surface-related multiple generation process (Berkhout 1982). Figure 1. View largeDownload slide (a) Feedback model at the surface, showing the generation of the surface-related multiples $$({{\bf M}}_0^ - = {{\bf P}}_0^ - {{{\bf A}}^ \cap }{{\bf P}}_{}^ - )$$. Note the important role of multiple generation operator [I + A∩P − ] in the downgoing (Q + ) and upgoing (P − ) wavefields at $${\rm z}_0^ +$$, where A∩ = (S + ) − 1R∩. (b) Feedback model at the surface, showing the generation of both the source ghost response and the surface-related multiples in the data (‘surface multiple scattering’). Note the extended multiple generation operator $$[{{\bf I}} + {{{\bf A}}^ \cap }{{{\underline{\bf P}}}^ - }]$$ in the downgoing (Q + ) and upgoing (P − ) wavefields. Extended multiple removal has a double function, i.e., simultaneously removal of the source ghost response and the surface-related multiples $$({{\bf M}}_0^ - )$$. Figure 1. View largeDownload slide (a) Feedback model at the surface, showing the generation of the surface-related multiples $$({{\bf M}}_0^ - = {{\bf P}}_0^ - {{{\bf A}}^ \cap }{{\bf P}}_{}^ - )$$. Note the important role of multiple generation operator [I + A∩P − ] in the downgoing (Q + ) and upgoing (P − ) wavefields at $${\rm z}_0^ +$$, where A∩ = (S + ) − 1R∩. (b) Feedback model at the surface, showing the generation of both the source ghost response and the surface-related multiples in the data (‘surface multiple scattering’). Note the extended multiple generation operator $$[{{\bf I}} + {{{\bf A}}^ \cap }{{{\underline{\bf P}}}^ - }]$$ in the downgoing (Q + ) and upgoing (P − ) wavefields. Extended multiple removal has a double function, i.e., simultaneously removal of the source ghost response and the surface-related multiples $$({{\bf M}}_0^ - )$$. Using the matrix operator notation (see also Appendix  D), the expression of the upgoing wavefields $$({{\bf P}}_{}^ - )$$ arriving at the surface (z0) can be written as:   $${{\bf P}}_{}^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = {{\bf X}}_{}^ \cup ({\rm z}_0^ + ,{\rm z}_0^ + ){{{\bf S}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}})$$ (1a)or   $${{{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = {{\bf X}}_0^ \cup ({\rm z}_0^ + ,{\rm z}_0^ + ){{{\bf Q}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}})$$ (1b)with $${\rm z}_0^ +$$ being the lower part of the surface and   $${{{\bf Q}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = {{{\bf S}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) + {{{\bf R}}^ \cap }({\rm z}_0^ + ,{\rm z}_0^ + ){{{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}}),$$ (1c)where $${{{\bf S}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = {{{\bf W}}^*}({\rm z}_0^ + ,{\rm z}_{\rm s}^{}){{{\bf S}}^ + }({{\rm z}_{\rm s}})$$ with W being the scattering-free propagation matrix and the complex conjugate version W* being its inverse. Note that source depth level zs may be generalized to zs(x, y), meaning that each column of matrix S + (zs) may have its own depth level. In the above expressions each column of wavefield matrix Q + represents the total downgoing wavefield that leaves the surface, each column of the source matrix S + represents the downgoing source (array) wavefield that excludes the source ghost (will be included later), one column of the reflectivity matrix $${{\bf R}}_{}^ \cap$$ represents the angle-dependent reflection operator for upgoing wavefields at each gridpoint of the surface and one column of the wavefield matrix $${{\bf P}}_{}^ -$$ represents the upgoing wavefields in one shot record. Finally, the surface-included earth's impulse response matrix $${{\bf X}}_{}^ \cup$$ transforms the downgoing source wavefields $${{\bf S}}_{}^ +$$ into the upgoing reflected wavefields $${{\bf P}}_{}^ -$$ and the surface-excluded earth's impulse response matrix $${{\bf X}}_0^ \cup$$ transforms the downgoing total wavefields $${{\bf Q}}_{}^ +$$ into the upgoing reflected wavefields $${{\bf P}}_{}^ -$$. In the ideal situation S + (zs) equals the unity matrix (S + = I) with a unit source element at each gridpoint of depth level zs. Note also that from the physics point of view, $${{\bf X}}_0^ \cup$$ represents one roundtrip (from surface into the subsurface and back to the surface) and $${{\bf X}}_{}^ \cup$$ represents many roundtrips. This can be easily seen from the feedback model in Fig. 1(a). For the response of a single, possibly blended source wavefield—blending meaning the creation of a decomposable areal source wavefield—matrix S + is reduced to column vector $$\vec{{\rm S}}_{\rm j}^ + = {{\bf S}}_{}^ + {\vec{\rm I}_{\rm j}}$$ and expression (1c) becomes a vector-matrix equation:   $$\vec{\rm Q}_{\rm j}^ + ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = \vec{{\rm S}}_{\rm j}^ + ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) + {{{\bf R}}^ \cap }({\rm z}_0^ + ,{\rm z}_0^ + )\vec{\rm P}_{\rm j}^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}),$$ (1d)showing that each gridpoint at the surface functions as one secondary source element of the naturally blended source array $${{{\bf R}}^ \cap }\vec{\rm P}_{\rm j}^ -$$. It is the response of this naturally blended source array that results in the surface-related multiples: $$({{\bf M}}_0^ - {\vec{\rm I}_{\rm j}} = {{\bf X}}_0^ \cup [{{{\bf R}}^ \cap }{{\bf P}}_{}^ - ]{\vec{\rm I}_{\rm j}} = {{\bf X}}_0^ \cup \delta {{\bf P}}_{}^ + {\vec{\rm I}_{\rm j}} = {{\bf X}}_0^ \cup \delta \vec{\rm P}_{\rm j}^ + )$$. If $${{\bf D}}_{}^ - ({{\rm z}_{\rm d}};{\rm z}_0^ + )$$ represents the detector matrix that selects the upgoing wavefields at detector level zd, where D − includes the detector (array) directivity as well, then the wavefield matrix $${{\bf P}}_{}^ -$$ is given by:   \begin{eqnarray} {{{\bf P}}^ - }({{\rm z}_{\rm d}};{{\rm z}_{\rm s}}) = {{{\bf D}}^ - }({{\rm z}_{\rm d}};{\rm z}_0^ + ){{\bf P}}_{}^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = {{{\bf D}}^ - }({{\rm z}_{\rm d}};{\rm z}_0^ + ){{\bf X}}_{}^ \cup ({\rm z}_0^ + ,{\rm z}_0^ + ){{{\bf S}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}})\end{eqnarray} (2a)or   $${{{\bf P}}^ - }({{\rm z}_{\rm d}};{{\rm z}_{\rm s}}) = {{{\bf D}}^ - }({{\rm z}_{\rm d}};{\rm z}_0^ + ){{\bf X}}_0^ \cup ({\rm z}_0^ + ,{\rm z}_0^ + ){{{\bf Q}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}})$$ (2b)with   $${{\bf D}}_{}^ - ({{\rm z}_{\rm d}};{\rm z}_0^ + ) = {{\bf D}}_{}^ - ({{\rm z}_{\rm d}}){{\bf W}}_{}^*({{\rm z}_{\rm d}},{\rm z}_0^ + ).$$ (2c)Similar to source level zs, detector level zd may be generalized to zd(x, y), meaning that each column of matrix D − (zd) may have its own depth level. Expressions (2a)–(2c) quantify the influence of acquisition matrices, $${{\bf S}}_{}^ + ({{\rm z}_{\rm s}})$$ and D − (zd) at the start and the end of the wavefield life cycle. In the 3-D practice both matrices are far from complete and, therefore, their influence is significant on the data space. This influence occurs in the common source gathers (columns of $${{\bf P}}_{}^ -$$) and common detector gathers (rows of $${{\bf P}}_{}^ -$$) as spatial aliasing and spatial band limitation. In industrial acquisition design, the elements of $${{\bf S}}_{}^ +$$ and D − are chosen such that aliasing and band limitation are minimized under pre-specified economical and logistical constraints. These constraints unavoidably lead to sparse 3-D acquisition matrices that cannot be directly inverted. We will see that this practical reality has an important influence on the design of effective multiple removal algorithms. The same is true for full wavefield redatuming algorithms and full wavefield imaging algorithms. Surface scattering removal for perfect data volumes Using expressions (1b) and (1c) and (2b) and (2c), we may write (see Fig. 1a):   \begin{eqnarray} {{{\bf P}}^ - }({{\rm z}_{\rm d}};{{\rm z}_{\rm s}}) &=& {{{\bf D}}^ - }({{\rm z}_{\rm d}};{\rm z}_0^ + ){{{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}})\nonumber\\ &=& {{{\bf D}}^ - }({{\rm z}_{\rm d}};{\rm z}_0^ + ){{\bf X}}_0^ \cup ({\rm z}_0^ + ,{\rm z}_0^ + ){{{\bf Q}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) \nonumber\\ &=& {{{\bf D}}^ - }({{\rm z}_{\rm d}};{\rm z}_0^ + ){{\bf X}}_0^ \cup ({\rm z}_0^ + ,{\rm z}_0^ + )[{{{\bf S}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) + {{{\bf R}}^ \cap }({\rm z}_0^ + ,{\rm z}_0^ + ){{{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}})] \nonumber\\ &=& {{\bf P}}_0^ - ({{\rm z}_{\rm d}};{\rm z}_{\rm s}^{})[{{\bf I}} + {{{\bf A}}^ \cap }({\rm z}_{\rm s}^{},{\rm z}_{\rm d}^{}){{{\bf P}}^ - }({\rm z}_{\rm d}^{};{{\rm z}_{\rm s}})], \end{eqnarray} (3a)where   $${{\bf P}}_0^ - ({{\rm z}_{\rm d}};{{\rm z}_{\rm s}}) = {{{\bf D}}^ - }({{\rm z}_{\rm d}};{\rm z}_0^ + ){{\bf X}}_0^ \cup ({\rm z}_0^ + ,{\rm z}_0^ + ){{{\bf S}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}}),$$ (3b)  $${{\bf A}}_{}^ \cap ({{\rm z}_{\rm s}},{{\rm z}_{\rm d}}) = {[{{{\bf S}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}})]^{ - 1}}{{{\bf R}}^ \cap }({\rm z}_0^ + ,{\rm z}_0^ + ){[{{\bf D}}_{}^ - ({{\rm z}_{\rm d}};{\rm z}_0^ + )]^{ - 1}}.$$ (3c)Note also that:   $${{{\bf Q}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = {{\bf S}}_{}^ + ({\rm z}_0^ + ;{{\rm z}_{\rm s}})[{{\bf I}} + {{{\bf A}}^ \cap }({{\rm z}_{\rm s}},{{\rm z}_{\rm d}}){{{\bf P}}^ - }({\rm z}_{\rm d}^{};{{\rm z}_{\rm s}})].$$ (3d)If we leave the source and detector level out of the notation, we may simplify expressions (3d) and (3a):   $${{{\bf Q}}^ + } = {{\bf S}}_{}^ + [{{\bf I}} + {{{\bf A}}^ \cap }{{{\bf P}}^ - }] \quad{\rm at}\quad {\rm z}_0^ +$$ (4a)  $${{{\bf P}}^ - } = {{\bf P}}_0^ - [{{\bf I}} + {{{\bf A}}^ \cap }{{{\bf P}}^ - }]\quad{\rm at}\quad {\rm z}_{\rm d}^{},$$ (4b)showing that [I + A∩P − ] is the surface multiple generation operator, transforming $${{\bf S}}_{}^ +$$ into $${{\bf Q}}_{}^ +$$ and $${{\bf P}}_0^ -$$ into P − . From expressions (4a) and (4b) it directly follows that   $${{\bf S}}_{}^ + = {{\bf Q}}_{}^ + {[{{\bf I}} + {{{\bf A}}^ \cap }{{{\bf P}}^ - }]^{ - 1}} \quad{\rm at}\quad {\rm z}_0^ + ,$$ (5a)  $${{\bf P}}_0^ - = {{{\bf P}}^ - }{[{{\bf I}} + {{{\bf A}}^ \cap }{{{\bf P}}^ - }]^{ - 1}} \quad{\rm at}\quad {\rm z}_{\rm d}^{},$$ (5b)showing that [I + A∩P − ] − 1 equals the surface multiple removal operator, transforming $${{\bf Q}}_{}^ +$$ into $${{\bf S}}_{}^ +$$ and P − into $${{\bf P}}_0^ -$$. Making use of the series expansion of the removal operator, direct inversion expression (5b) may also be written as:   \begin{eqnarray} {{\bf P}}_0^ - &=& {{{\bf P}}^ - } - {{{\bf P}}^ - }[{{{\bf A}}^ \cap }{{{\bf P}}^ - }] + {{{\bf P}}^ - }{[{{{\bf A}}^ \cap }{{{\bf P}}^ - }]^2} - {{{\bf P}}^ - }{[{{{\bf A}}^ \cap }{{{\bf P}}^ - }]^3} + \cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\nonumber\\ &=& {{{\bf P}}^ - } - [{{{\bf P}}^ - }{{{\bf A}}^ \cap }]{{{\bf P}}^ - } + {[{{{\bf P}}^ - }{{{\bf A}}^ \cap }]^2}{{{\bf P}}^ - } - {[{{{\bf P}}^ - }{{{\bf A}}^ \cap }]^3}{{{\bf P}}^ - } + \cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\quad {\rm at}\quad {\rm z}_{\rm d}^{}. \end{eqnarray} (5c)Expression (5c) represents the well-known inverse scattering series at the surface (Berkhout 1982; Araujo et al.1994; Verschuur & Berkhout 1997). If all elements of surface operator A∩ exist and are known (‘perfect data’), we have two ways of surface-multiple removal by applying direct inversion: Eq. (5b) is used, meaning that multiple-generating matrix operator [I + A∩P − ] need be inverted for; Inverse scattering series (5c) is used, meaning that each previous term need be multiplied with [A∩P − ] at the source side or [P − A∩] at the detector side. Now, let us return to expressions (4a) and (4b). If we pose that the determination of impulse response (IR) matrix $${{\bf X}}_0^ \cup$$ is the ultimate aim of surface-related pre-processing (leading to a subsurface response that has been made independent of the conditions at the surface), then we may write:   $${{\bf X}}_0^ \cup = {({{{\bf D}}^ - })^{ - 1}}{{{\bf P}}^ - }{({{{\bf Q}}^ + })^{ - 1}} = {({{{\bf D}}^ - })^{ - 1}}\frac{{{{\bf P}}_0^ - [{{\bf I}} + {{{\bf A}}^ \cap }{{{\bf P}}^ - }]}}{{{{\bf S}}_{}^ + [{{\bf I}} + {{{\bf A}}^ \cap }{{{\bf P}}^ - }]}} = {({{{\bf D}}^ - })^{ - 1}}{{\bf P}}_0^ - {({{\bf S}}_{}^ + )^{ - 1}},$$ (6)assuming again that the elements of A∩ exist at every gridpoint of the source and detector depth levels. Expression (6) reveals that in the situation of a perfect data volume (no noise, no missing data), surface multiples do not provide extra information for the determination of subsurface IR matrix $${{\bf X}}_0^ \cup$$! Expression (6) also explains that for the determination of $${{\bf X}}_0^ \cup$$ the complex division with double illumination matrix $$({{\bf Q}}_{}^ + )$$ can be avoided by (i) surface multiple removal using the inverse scattering series, followed by (ii) source-detector deconvolution:   $$({\rm i})\,\,{{\bf P}}_0^ - = {{{\bf P}}^ - }{({{\bf I}} + {{{\bf A}}^ \cap }{{{\bf P}}^ - })^{ - 1}} = {{{\bf P}}^ - } - {{{\bf P}}^ - }({{{\bf A}}^ \cap }{{{\bf P}}^ - }) + {{{\bf P}}^ - }{({{{\bf A}}^ \cap }{{{\bf P}}^ - })^2} - {{{\bf P}}^ - }{({{{\bf A}}^ \cap }{{{\bf P}}^ - })^3} + \cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot$$ (7a)or   $${{\bf P}}_0^ - = {({{\bf I}} + {{{\bf P}}^ - }{{{\bf A}}^ \cap })^{ - 1}}{{{\bf P}}^ - } = {{{\bf P}}^ - } - ({{{\bf P}}^ - }{{{\bf A}}^ \cap }){{{\bf P}}^ - } + {({{{\bf P}}^ - }{{{\bf A}}^ \cap })^2}{{{\bf P}}^ - } - {({{{\bf P}}^ - }{{{\bf A}}^ \cap })^3}{{{\bf P}}^ - } + \cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot$$ (7b)and   $$({\rm ii})\,\,{{\bf X}}_0^ \cup = {({{{\bf D}}^ - })^{ - 1}}{{\bf P}}_0^ - {({{{\bf S}}^ + })^{ - 1}}$$ (7c)at each gridpoint on detector level zd and source level zs respectively. Surface scattering removal for sparse data volumes Due to economical and logical constraints in practice, D − and S + are sparse matrices, meaning that [D − ] − 1 and [S + ] − 1 do not exist and therefore operator A∩ does not exist (from theoretical to real data). Moreover, $${{\bf P}}_{}^ -$$ contains always noise and full knowledge of true source and detector behaviour is not available (D − and S + are only approximately known). This means that direct inversion as formulated by expression (6) causes significant problems in practice. In the following, an alternative way of thinking is introduced with respect handling multiples in sparse data volumes. First, let us reveal a fundamental property of the multiple-generating system. From expression (4b) it follows that:   $$\begin{array}{@{}l@{}} {{\bf P}}_0^ - = {{{\bf P}}^ - } - {{\bf P}}_0^ - {{{\bf A}}^ \cap }{{{\bf P}}^ - } = [{{\bf I}} - {{\bf P}}_0^ - {{{\bf A}}^ \cap }]{{{\bf P}}^ - }\quad {\rm at} \quad {\rm z}_{\rm d}^{}, \end{array}$$ (8a)showing that $${{\bf P}}_0^ - {{{\bf A}}^ \cap }$$ functions as a prediction operator. Using expression (7b), we may also write:   $${{\bf P}}_0^ - = {[{{\bf I}} + {{{\bf P}}^ - }{{{\bf A}}^ \cap }]^{ - 1}}{{{\bf P}}^ - }\quad {\rm at} \quad {\rm z}_{\rm d}^{}.$$ (8b)If we compare expression (8a) with expression (8b) it follows that:   $${[{{\bf I}} + {{{\bf P}}^ - }{{{\bf A}}^ \cap }]^{ - 1}} = [{{\bf I}} - {{\bf P}}_0^ - {{{\bf A}}^ \cap }]\quad {\rm at} \quad {\rm z}_{\rm d}^{}$$ (9a)or   $${[{{\bf I}} - {{\bf P}}_0^ - {{{\bf A}}^ \cap }]^{ - 1}} = [{{\bf I}} + {{{\bf P}}^ - }{{{\bf A}}^ \cap }]\quad {\rm at} \quad {\rm z}_{\rm d}^{}.$$ (9b)Expressions (9a) and (9b) tell us that if surface operator A∩ exists, meaning that at every gridpoint wavefields are available, we may conclude that: The inverse of causal operator [I + P − A∩] exists The inversion result $$[{{\bf I}} - {{\bf P}}_0^ - {{{\bf A}}^ \cap }]$$ is again causal. These two conditions are exactly the two necessary and sufficient conditions for an operator to have the minimum-phase property (Berkhout 1971). Hence, expression (9b) extends the classical minimum-phase property of a time-series to the 3-D seismic surface multiple-generating system. Note that the minimum-phase condition for a wavefield is significantly more than the causality condition. In all areas of science (think particularly at today's climate sciences) the practical relevance of minimum phase is that the future can be fully predicted by knowledge of the past! Later in this paper, it will be shown that this fundamental property also applies to the Earth's internal multiple-generating system. Minimum-phase test for multiple removal solutions From expression (9b) it follows that solution $${{\bf P}}_0^ -$$ must obey the condition:   $$[{{\bf I}} - {{\bf P}}_0^ - {{{\bf A}}^ \cap }][{{\bf I}} + {{{\bf P}}^ - }{{{\bf A}}^ \cap }] = {{\bf I}} \qquad \hbox{(minimum-phase condition)}$$ (10a)or, simplifying the notation by writing $${{{\bf Y}}_0} = {{\bf P}}_0^ - {{{\bf A}}^ \cap }$$ and Y = P − A∩,   $$$$[{{\bf I}} - {{\bf Y}}_0^{}][{{\bf I}} + {{\bf Y}}] = {{\bf I}}.$$$$ (10b)If $${{{\bf \hat{Y}}}_0}$$ represents some initial estimate of Y0, then error matrix E0 is given by:   $${{{\bf E}}_0} = [{{\bf I}} - {{\bf \hat{Y}}}_0^{}][{{\bf I}} + {{\bf Y}}] - {{\bf I}}.$$ (11a)Note that E0 is a measure for the violation of the minimum-phase property due to a wrong prediction result. It is also a measure for the difference $$[{{{\bf Y}}_0} - {{{\bf \hat{Y}}}_0}]$$. If we write $${{{\bf Y}}_0} = {{{\bf \hat{Y}}}_0} + \Delta {{{\bf Y}}_0}$$, then it follows from expression (11a) that   $$[{{\bf I}} - {{\bf \hat{Y}}}_0^{}][{{\bf I}} + {{\bf Y}}] = {{\bf I}} + \Delta {{\bf Y}}_0^{}[{{\bf I}} + {{\bf Y}}] = {{\bf I}} + {{{\bf E}}_0}.$$ (11b)Using expressions (11a) and (11b) we may conclude that at $${\rm z}_0^ +$$:   $${\rm (1)}\,\,{{{\bf E}}_0} = [{{\bf I}} - {{\bf \hat{Y}}}_0^{}][{{\bf I}} + {{\bf Y}}] - {{\bf I}} = {{\bf Y}} - {{\bf \hat{Y}}}_0^{} - {{\bf \hat{Y}}}_0^{}{{\bf Y}}$$ (12a)  $${\rm (2)}\,\,\Delta {{\bf Y}}_0^{} = {{{\bf E}}_0}{[{{\bf I}} + {{\bf Y}}]^{ - 1}} \approx {{{\bf E}}_0}[{{\bf I}} - {{{\bf \hat{Y}}}_0}].$$ (12b)In Appendix  A the iterative solution scheme is shown, allowing also an error in input matrix Y. In the following, another application of the minimum-phase property is proposed. Today, many surface multiple removal algorithms are in use. Most of them belong to the data-driven SRME family, based on the feedback model (Fig. 1a). Although different in-house versions show different performances, at the end of the day all multiple scattering solutions must obey the minimum-phase condition (10a). This does not only apply to the output data $${{\bf P}}_0^ -$$ but also to the interpolated data at the missing data points of input data P − (Appendix  A). Let us assume that the performance of a specific multiple removal algorithm need be objectively tested. Then input-output product $$[{{\bf I}} - {{\bf \hat{Y}}}_0^{}][{{\bf I}} + {{\bf \hat{Y}}}]$$ can be easily computed (product of two data matrices). If we exclude the local data around the origin (‘offset = 0, time = 0’), each column of this product represents an error shot record $$({{{\bf E}}_0}{\vec{\rm I}_{\rm j}})$$. The smaller the samples in the error records, the better the performance of the algorithm. Expression (12b) shows that the error records can be used to improve the solution (see Appendix  A):   $${({{\bf P}}_0^ - {{{\bf A}}^ \cap })^{\rm (i)}} = {({{\bf P}}_0^ - {{{\bf A}}^ \cap })^{({\rm i} - 1)}} + {{\bf E}}_0^{({\rm i} - 1)}[{{\bf I}} - {({{\bf P}}_0^ - {{{\bf A}}^ \cap })^{({\rm i} - 1)}}]$$ (13a)  $${({{{\bf P}}^ - }{{{\bf A}}^ \cap })^{(\rm i)}} = {({{{\bf P}}^ - }{{{\bf A}}^ \cap })^{({\rm i} - 1)}} + [{{\bf I}} + {({{\bf P}}_{}^ - {{{\bf A}}^ \cap })^{({\rm i} - 1)}}{{{\bf E}}^{({\rm i} - 1)}}],$$ (13b)where i = 1 represents the first improvement step. For algorithms with a very good performance, one step may be already sufficient. Note that $${{{\bf P}}_0}{{{\bf A}}^ \cap } = {{\bf X}}_0^ \cup {{{\bf R}}^ \cap }$$. Hence, by using knowledge of the surface reflectivity the IR matrix $${{\bf X}}_0^ \cup$$ is known (of course within the constraint of the limited seismic bandwidth). If we apply the Fourier transform with respect to the source and detector coordinates (meaning multiplication by Fourier matrix F and its Hermitian version FH), minimum-phase condition (10a) becomes a condition for a plane wave response.   $$[{{\bf I}} - ({{\bf F}}{{\bf P}}_0^ - {{{\bf F}}^{\rm H}})({{\bf F}}{{{\bf A}}^ \cap }{{{\bf F}}^{\rm H}})][{{\bf I}} + ({{\bf F}}{{{\bf P}}^ - }{{{\bf F}}^{\rm H}})({{\bf F}}{{{\bf A}}^ \cap }{{{\bf F}}^{\rm H}})] = {{\bf I}}$$ (14a)or   $$[{{\bf I}} - {{\bf \tilde{P}}}_0^ - {{{\bf \tilde{A}}}^ \cap }][{{\bf I}} + {{{\bf \tilde{P}}}^ - }{{{\bf \tilde{A}}}^ \cap }] = {{\bf I}}.$$ (14b)For a 1-D subsurface the data matrices $${{\bf \tilde{P}}}_{}^ - {\rm and}\, {{\bf \tilde{P}}}_0^ -$$ become diagonal matrices, each diagonal element representing one frequency component of a plane wave response, and minimum phase condition (14b) can be simplified to the minimum-phase condition of a time-series for each plane wave (Berkhout 1971, 1974). Including the source ghost in surface multiple scattering In marine data we deal with prominent ghosts at the source and the detector side:   $${{\bf S}}({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = {{{\bf W}}^*}({\rm z}_0^ + ,{\rm z}_{\rm s}^{}){{{\bf S}}^ + }({{\rm z}_{\rm s}}) + {{{\bf R}}^ \cap }({\rm z}_0^ + ,{\rm z}_0^ + ){{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm s}^{}){{{\bf S}}^ - }({{\rm z}_{\rm s}})$$ (15a)  $${{\bf D}}({{\rm z}_{\rm d}};{\rm z}_0^ + ) = {{\bf D}}_{}^ - ({{\rm z}_{\rm d}}){{\bf W}}_{}^*({{\rm z}_{\rm d}},{\rm z}_0^ + ) + {{\bf D}}_{}^ + ({{\rm z}_{\rm d}}){{\bf W}}({{\rm z}_{\rm d}},{\rm z}_0^ + ){{{\bf R}}^ \cap }({\rm z}_0^ + ,{\rm z}_0^ + ),$$ (15b)leading to generalization of expression (2a): $${{\bf P}}({{\rm z}_{\rm d}};{{\rm z}_{\rm s}}) = {{{\bf D}}^ - }({{\rm z}_{\rm d}};{\rm z}_0^ + ){{\bf X}}_{}^ \cup ({\rm z}_0^ + ,{\rm z}_0^ + ){{\bf S}}({\rm z}_0^ + ;{{\rm z}_{\rm s}})$$. It is common practice to apply detector deghosting to the data. Looking at expressions (15b), this implies ‘division’ by the operator $${{\bf D}}({{\rm z}_{\rm d}};{\rm z}_0^ + )$$. Although it is for the theoretical framework not required, let us assume that this detector deghosting process has been successfully applied. For source deghosting, the situation may be much more complex due to the nonlinear interference of the source with the surface during firing time. It leads to different values of operator R∩ in (15a) and (15b), meaning that reciprocity does not hold anymore. An alternative solution is proposed by utilizing the source ghost (R∩W − S − ) for extra illumination (extended double illumination), by including it in the Q-matrix (see Fig. 1b):   \begin{eqnarray} {{{\bf Q}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) &=& {{{\bf W}}^*}({\rm z}_0^ + ,{\rm z}_{\rm s}^{}){{{\bf S}}^ + }({{\rm z}_{\rm s}}) + {{{\bf R}}^ \cap }({\rm z}_0^ + ,{\rm z}_0^ + )[{{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm s}^{}){{{\bf S}}^ - }({{\rm z}_{\rm s}}) + {{{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}})] \nonumber\\ &=& {{{\bf W}}^*}({\rm z}_0^ + ,{\rm z}_{\rm s}^{}){{{\bf S}}^ + }({{\rm z}_{\rm s}}) + {{{\bf R}}^ \cap }({\rm z}_0^ + ,{\rm z}_0^ + ){{{\underline{\bf P}}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) \nonumber \\ & =& {{{\bf S}}^ + }({\rm z}_0^ + ;{\rm z}_{\rm s}^{}) + \delta {{{\underline{\bf P}}}^ + }({\rm z}_0^ + ;{\rm z}_{\rm s}^{}) \end{eqnarray} (16)with $${{\underline{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = {{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm s}^{}){{{\bf S}}^ - }({{\rm z}_{\rm s}}) + {{{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}})$$ and $$\delta {{\underline{\bf P}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = {{{\bf R}}^ \cap }({\rm z}_0^ + ,{\rm z}_0^ + ){{\underline{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}})$$. The rationale for including this option is that the ghost source and the multiple source are both coded by the surface reflectivity operator R∩. It is also interesting to realize that the expression of the data with surface-related multiples and source ghost response (‘surface multiple scattering’) is given by (see Fig. 1b):   $${{\bf P}}_{}^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = {{\bf P}}_0^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) + {{\bf P}}_0^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}){{{\bf A}}^ \cap }({{\rm z}_{\rm s}};{\rm z}_0^ + ){\underline{\bf P }}_{}^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}})$$ (17a)with   $${{\underline{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = {{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm s}^{}){{{\bf S}}^ - }({{\rm z}_{\rm s}}) + {{{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}}).$$ (17b)Note that the extended multiple removal process removes simultaneously the multiples of both the direct source and the ghost source: output $${{\bf P}}_0^ -$$ is multiple-free and ghost-free. Imaging with surface multiple scattering The expression of the upgoing wavefields that arrive at the surface $$({\rm z}_0^ + )$$ is given by:   \begin{eqnarray}{{\underline{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = {{\bf X}}_0^ \cup ({\rm z}_0^ + ,{\rm z}_0^ + ){{{\bf Q}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) + {{{\bf W}}^ - }({\rm z}_0^ + ,{{\rm z}_{\rm s}}){{{\bf S}}^ - }({{\rm z}_{\rm s}}),\end{eqnarray} (18a)where   \begin{eqnarray} {{{\bf Q}}^ + }({\rm z}_0^ + ;{\rm z}_{\rm s}^{}) &=& {{{\bf S}}^ + }({\rm z}_0^ + ;{\rm z}_{\rm s}^{}) + {{{\bf R}}^ \cap }({\rm z}_0^ + ,{\rm z}_0^ + ){{{\underline{\bf P}}}^ - }({\rm z}_0^ + ;{\rm z}_{\rm s}^{})\nonumber\\ &=& {{{\bf S}}^ + }({\rm z}_0^ + ;{\rm z}_{\rm s}^{})[{{\bf I}} + {{{\bf A}}^ \cap }({\rm z}_{\rm s}^{},{\rm z}_0^ + ){{{\underline{\bf P}}}^ - }({\rm z}_0^ + ;{\rm z}_{\rm s}^{})] \end{eqnarray} (18b)represents extended double illumination. If we neglect the internal multiples for the moment, we may write for the Earth's first-order transfer matrix (WRW-model):   $${{\bf X}}_0^ \cup ({\rm z}_0^ + ,{\rm z}_0^ + ) = \sum\limits_{\rm n = 1}^{{\rm m} - 1} {[{{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm n}^ - )} {{{\bf R}}^ \cup }({\rm z}_{\rm n}^ - ,{\rm z}_{\rm n}^ - ){{\bf W}}({\rm z}_{\rm n}^ - ,{\rm z}_0^ + )] + {{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm m}^ - ){\underline{\bf {X}}}_{}^ \cup ({\rm z}_{\rm m}^ + ,{\rm z}_{\rm m}^ + ){{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_0^ + ),$$ (18c)where matrix W equals the wavefield propagation matrix and where one column of reflectivity matrix R∪ equals the angle-dependent reflection operator for the downgoing wavefields at one gridpoint of depth level zm. Let us look at perfect data volumes first. Imaging with surface multiple scattering for perfect data In the seismic migration algorithm the downward travelling wavefields (Q + ) that leave $${\rm z}_0^ +$$ are numerically transported into the subsurface by forward extrapolation. It yields the full illuminating wavefields (P + ) that arrive at the upper part of depth level $${\rm z}_{\rm m}^{}$$, being indicated by $${\rm z}_{\rm m}^ -$$ (m = 1, 2, ……). Again, if we neglect the internal multiples for the moment (see Fig. 2):   \begin{eqnarray} {{{\bf P}}^ +}({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}) &=& {{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_0^ + ){{{\bf Q}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) \nonumber \\ &=& {{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_0^ + ){{{\bf S}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}})[{\bf I} + {{{\bf A}}^ \cap }({\rm z}_{\rm s}^{},{\rm z}_0^ + ){{{\underline{\bf P}}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}})] \nonumber \\ &=& {{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_{\rm s}^{}){{{\bf S}}^ + }({{\rm z}_{\rm s}})[{\bf I} + {{{\bf A}}^ \cap }({\rm z}_{\rm s}^{},{\rm z}_0^ + ){{{\underline{\bf P}}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}})]\nonumber \\ &=& {{\bf P}}_0^ + ({\rm z}_{\rm m}^ - ;{\rm z}_{\rm s}^ + )[{\bf I} + {{{\bf A}}^ \cap }({\rm z}_{\rm s}^{},{\rm z}_0^ + ){{{\underline{\bf P}}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}})], \end{eqnarray} (19a)where $${{\bf P}}_0^ + = {{\bf W}}_{}^ + {{\bf S}}_{}^ +$$ represents the illuminating wavefield without surface ghosts and without multiples. Similarly, the subsurface response—being represented by the upward travelling wavefields P − arriving at $${\rm z}_0^ +$$—are numerically transported into the subsurface by inverse extrapolation, yielding the response wavefields Q − that leave $${\rm z}_{\rm m}^ -$$ (m = 1, 2, ……):   \begin{eqnarray} {{{\bf Q}}^ - }({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}) &=& {{{\bf W}}^*}({\rm z}_{\rm m}^ - ,{\rm z}_0^ + )\left[ {{{{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) - \sum\limits_{\rm n = 1}^{{\rm m} - 1} {{{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm n}^ - )} {{\bf R}}_{}^ \cup ({\rm z}_{\rm n}^ - ,{\rm z}_{\rm n}^ - ){{{\bf P}}^ + }({\rm z}_{\rm n}^ - ;{{\rm z}_{\rm s}})} \right] \nonumber\\ &=& {{{\bf W}}^*}({\rm z}_{\rm m}^ - ,{\rm z}_0^ + )\left[ {{{\bf P}}_0^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) - \sum\limits_{\rm n = 1}^{{\rm m} - 1} {{{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm n}^ - )} {{\bf R}}_{}^ \cup ({\rm z}_{\rm n}^ - ,{\rm z}_{\rm n}^ - ){{\bf P}}_0^ + ({\rm z}_{\rm n}^ - ;{{\rm z}_{\rm s}})} \right]\left[ {{{\bf I}} + {{\bf A}}_{}^ \cap ({\rm z}_{\rm s}^{},{\rm z}_0^ + ){{{\underline{\bf P}}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}})} \right] \nonumber\\ &=& {{\bf Q}}_0^ - ({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})\left[ {{{\bf I}} + {{\bf A}}_{}^ \cap ({\rm z}_{\rm s}^{},{\rm z}_0^ + ){{{\underline{\bf P}}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}})} \right], \end{eqnarray} (19b)where $${{\bf Q}}_0^ -$$ represents the response wavefields at $${\rm z}_{\rm m}^ -$$without surface ghosts and without multiples, the term ∑W − R∪P + represents the response wavefields from shallower depth levels $${\rm z}_{\rm n}^ -$$ (zn < zm) that need be removed and $${{\bf P}}_0^ + = {{{\bf W}}^ + }{{{\bf S}}^ + }$$ being the direct source wavefields that arrive at $${\rm z}_{\rm n}^ -$$. Note that subtraction term ∑W − R∪P + is needed to keep Q − causal. At level $${\rm z}_{\rm m}^ -$$ full response wavefields Q − and full incident wavefields P + are related according to (see Fig. 2):   $${{{\bf Q}}^ - }({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}) = {{{\underline{\bf X}} }}_{}^ \cup ({\rm z}_{\rm m}^ - ,{\rm z}_{\rm m}^ - ){{{\bf P}}^ + }({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}).$$ (20a) Figure 2. View largeDownload slide Feedback model for the wavefields at depth level $${\rm z}_{\rm m}^ -$$ for the simplified situation that the internal multiples in the overburden can be neglected. Note that in practice, surface operator R∩ may be a time-variant operator due to nonlinearity at the source and a changing sea state. Figure 2. View largeDownload slide Feedback model for the wavefields at depth level $${\rm z}_{\rm m}^ -$$ for the simplified situation that the internal multiples in the overburden can be neglected. Note that in practice, surface operator R∩ may be a time-variant operator due to nonlinearity at the source and a changing sea state. Applying direct inversion and using expressions (19a) and (19b), transfer operator $${\underline{\bf {X}}}_{}^ \cup$$ is given by:   \begin{eqnarray} {\underline{\bf {X}}}_{}^ \cup ({\rm z}_{\rm m}^ - ,{\rm z}_{\rm m}^ - ) &=& {{\bf Q}}_{}^ - ({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}){[{{\bf P}}_{}^ + ({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}})]^{ - 1}}\nonumber\\ &=& {{\bf Q}}_0^ - ({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}){[{{\bf P}}_0^ + ({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}})]^{ - 1}}. \end{eqnarray} (20b)Expression (20b) reveals that if we deal with perfect data (no noise, no missing data), direct inversion can be applied and the result is insensitive to the surface multiple scattering (source ghost response + the surface-related multiples). In other words, with perfect data the surface multiple scattering does not give an extra contribution to the estimate of transfer operator $${\underline{\bf {X}}}_{}^ \cup$$. Note that $${\underline{\bf {X}}}_{}^ \cup$$ contains the reflectivity at $${\rm z}_{\rm m}^ -$$, i.e., $${\underline{\bf {X}}}_{}^ \cup = {{\bf R}}_{}^ \cup + {{\bf X}}_{}^ \cup$$. This means that for perfect data R∪ can be completely determined by surface-free response $${{\bf P}}_0^ -$$. If we apply Fourier transformation with respect to the detector coordinates on depth level $${\rm z}_{\rm m}^ -$$ (multiplication by Fourier matrix F to expression (20b) on the detector side), reflection coefficient matrix $${{{\bf \tilde{R}}}^ \cup }$$ can be found at zero intercept time:   $${{\bf \tilde{R}}}_{}^ \cup ({\rm z}_{\rm m}^ - ,{\rm z}_{\rm m}^ - ) = {{\bf F}}({\rm z}_{\rm m}^ - ){{\bf Q}}_{}^ - ({\rm z}_{\rm m}^ - ;{{\rm z}_0}){[{{\bf P}}_{}^ + ({\rm z}_{\rm m}^ - ;{{\rm z}_0})]^{ - 1}} = {{\bf F}}({\rm z}_{\rm m}^ - ){{\bf Q}}_0^ - ({\rm z}_{\rm m}^ - ;{{\rm z}_0}){[{{\bf P}}_0^ + ({\rm z}_{\rm m}^ - ;{{\rm z}_0})]^{ - 1}} \quad {\rm at} \quad \tau = 0,$$ (20c)where $${{\bf \tilde{R}}}_{}^ \cup = {{\bf F}}{{\bf R}}_{}^ \cup$$. Bear in mind that the columns of $${{{\bf \tilde{R}}}^ \cup }$$ represent the angle-dependent plane-wave reflection coefficients at each gridpoint of depth level $${\rm z}_{\rm m}^ -$$. Note also that expression (20c) is consistent with expression (20b) and confirms that in the situation of perfect data surface multiple scattering does not contribute to the redatuming result and to the image. Next, we will focus on the realistic situation of sparse data volumes. Imaging with surface multiple scattering for sparse data In the situation of sparse data, direct inversion combined with applying the imaging principle is replaced by ‘reflectivity estimation by minimization’ (m = 1, 2, …….):   $${\rm E}({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}) = \sum\limits_j {\sum\limits_\omega {{{\left\| {[{{\bf Q}}_{}^ - ({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}) - {{\bf R}}_{}^ \cup ({\rm z}_{\rm m}^ - ,{\rm z}_{\rm m}^ - ){{\bf P}}_{}^ + ({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}})]{{\vec{{\rm I}}}_{\rm j}}} \right\|}^2}} } = {\rm minimum},$$ (21a)where j indicates the shot number, $${\vec{\rm I}_{\rm j}}$$ selects the right column,   $${{{\bf Q}}^ - }({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}) = {{{\bf W}}^*}({\rm z}_{\rm m}^ - ,{\rm z}_0^ + )\left[ {{{{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) - \sum\limits_{\rm n = 1}^{{\rm m} - 1} {{{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm n}^ - )} {{\bf R}}_{}^ \cup ({\rm z}_{\rm n}^ - ,{\rm z}_{\rm n}^ - ){{{\bf P}}^ + }({\rm z}_{\rm n}^ - ;{{\rm z}_{\rm s}})} \right]$$ (21b)represents the causal upgoing wavefields that leave depth level $${\rm z}_{\rm m}^ -$$ and   $${{\bf P}}_{}^ + ({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}) = {{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_{\rm s}^{}){{\bf S}}_{}^ + ({{\rm z}_{\rm s}}) + {{\bf W}}({\rm z}_{\rm m}^ - ;{\rm z}_0^ + )[{{\bf R}}_{}^ \cap ({\rm z}_0^ + ,{\rm z}_0^ + ){\underline{\bf P}}_{}^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}})]$$ (21c)represents the extended double illuminating wavefields that arrive at depth level $${\rm z}_{\rm m}^ -$$. Note that (21a) physically means ‘optimized double focusing’ (Berkhout 1982). Using expression (21c), minimization by (21a) can be refined to (notation is simplified by omitting the depth levels):   $$E = \sum\limits_{\rm j} {\sum\limits_\omega {{{\left\| {[{{\bf Q}}_{}^ - - {{\bf R}}_{\rm a}^ \cup {{{\bf W}}^ + }{{\bf S}}_{}^ + - {{\bf R}}_{\rm b}^ \cup {{{\bf W}}^ + }\delta {\underline{\bf P}}_{}^ + ]{{\vec{{\rm I}}}_{\rm j}}} \right\|}^2}} } = {\rm minimum\,at}\ {\rm z}_{\rm m}^ - ,$$ (22)where $$\delta {{\underline{\bf {P} }}^ + } = {{{\bf R}}^ \cap }{{\underline{\bf {P} }}^ - }$$ at $${\rm z}_0^ +$$ and $${{\bf R}}_{\rm a}^ \cup {{{\bf W}}^ + }{{\bf S}}_{}^ + + {{\bf R}}_{\rm b}^ \cup {{{\bf W}}^ + }\delta {\underline{\bf {P} }}_{}^ + = {{\bf R}}_{}^ \cup {{\bf P}}_{}^ +$$ represents extended wide-angle reflection process at depth level $${\rm z}_{\rm m}^ -$$. From minimization expression (22) we see the benefit of reflectivity estimation by utilizing the surface multiple scattering: Realistic data volumes, represented by parse data matrices can be handled (from direct inversion to iterative minimization); Images can be computed per shot record, or per receiver gather, meaning that only one domain need be properly sampled (impossible with data-driven multiple removal); By utilizing surface multiple scattering two images are obtained (‘image decomposition’), one from the primary source response $$({\rm yielding}\,\,{{\bf R}}_{\rm a}^ \cup )$$ and one from the secondary source response $$({\rm yielding}\,\,{{\bf R}}_{\rm b}^ \cup )$$. Note the significant advantage: instead of applying an elaborate removal process, the ghost source response and the surface multiples are utilized. From the above we may conclude that surface multiple scattering (source ghost response + surface multiples) should not be removed but it should be utilized. In addition, we may conclude that image decomposition$$({{\bf R}}_{}^ \cup\, {\rm into}\,{{\bf R}}_{\rm a}^ \cup ,{{\bf R}}_{\rm b}^ \cup )$$ also allows us to separate primaries and surface multiples by using the WRW-model (‘wavefield decomposition’):   $$$${{\bf P}}_{\rm pr}^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = \sum\limits_{\rm m = 1}^{\rm M} {{{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm m}^ - ){{\bf R}}_{\rm a}^ \cup ({\rm z}_{\rm m}^ - ,{\rm z}_{\rm m}^ - ){{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_{\rm s}^{}){{\bf S}}_{}^ + ({{\rm z}_{\rm s}})}\quad \left( {\rm primaries} \right)$$$$ (23a)  $${{\bf M}}_{\rm sfc}^ - ({\rm z}_0^ + ;{{\rm z}_0}) = \sum\limits_{\rm m = 1}^{\rm M} {{{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm m}^ - ){{\bf R}}_{\rm b}^ \cup ({\rm z}_{\rm m}^ - ,{\rm z}_{\rm m}^ - ){{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_0^ + )\delta {\underline{\bf {P} }}_{}^ + ({\rm z}_0^ + ;{{\rm z}_0})}\quad \left( {\rm surface\, multiples} \right)$$ (23b)with $$\delta {{\underline{\bf {P} }}^ + } = {{{\bf R}}^ \cap }({{{\bf W}}^ - }{{{\bf S}}^ - } + {{{\bf P}}^ - })$$. Note also that all wavefield components include the missing data as well (wavefield reconstruction) and superposition leads to the full response: $${{\bf P}}_{}^ - = {{\bf P}}_{\rm pr}^ - + {{\bf M}}_{\rm sfc}^ -$$ at the surface, where $${{\bf M}}_{\rm sfc}^ -$$ includes the primaries from the ghost source. Of course, both images $$({{\bf R}}_{\rm a}^ \cup ,{{\bf R}}_{\rm b}^ \cup )$$ can be combined into one image $$({{\bf R}}_{}^ \cup )$$ at each subsurface gridpoint k by a data-driven illumination-weighted addition at each subsurface gridpoint k:   $$\vec{\rm R}_{}^ \cup ({\rm k}) = [{\rm A}_{\rm k}^2\vec{\rm R}_{\rm a}^ \cup ({\rm k}) + {\rm B}_{\rm k}^2\vec{R}_{\rm b}^ \cup ({\rm k})]/[{\rm A}_{\rm k}^2 + {\rm B}_{\rm k}^2]$$ (24)with $${\rm A}_{\rm k}^2 = \sum\nolimits_{\rm j} {{\rm a}_{\rm kj}^2}\, {\rm and}\, {\rm B}_{\rm k}^2 = \sum\nolimits_{\rm j} {\rm b_{\rm kj}^2}$$, where akj equals the measured energy of the primary incident wavefield at gridpoint k and $${\rm b_{\rm kj}}$$ equals the measured energy of the secondary incident wavefield at gridpoint k, both being initially generated by primary source j. Using the prior information that reflectivity is an angle-dependent scalar at each depth level (R does not contain traveltimes), wavefield separation and wavefield reconstruction should be an integral part of the closed-loop migration algorithm. In the following we will see that the same conclusion holds for the internal multiple generating system. Feedback model below the surface $$({\rm z}_1^ + ,{\rm z}_2^ + ,{\rm z}_3^ + ,........)$$ Fig. 3(a) visualizes the feedback model at $${\rm z}_{\rm m}^ +$$. Note that IR matrix $${{\underline{\bf {X}} }}_{}^ \cap$$ contains all multiples in the overburden $$(z < {\rm z}_{\rm m}^ + )$$. If we would neglect the internal multiples, we may write:   $$\delta {\underline{\bf {X}}}_{}^ \cap ({\rm z}_{\rm m}^ + ,{\rm z}_{\rm m}^ + ) = {{\bf W}}({\rm z}_{\rm m}^ + ,{\rm z}_0^ + ){{\bf R}}_{}^ \cap ({\rm z}_0^ + ,{\rm z}_0^ + ){{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm m}^ + ),$$ (25a)leading to Fig. 2. An interesting refinement of expression (25a) is given by:   $${\underline{\bf {X}}}_{\rm pr}^ \cap ({\rm z}_{\rm m}^ + ,{\rm z}_{\rm m}^ + ) = \sum\limits_{\rm n = 0}^{\rm m} {[{{\bf W}}({\rm z}_{\rm m}^ + ,{\rm z}_{\rm n}^ + ){{\bf R}}_{}^ \cap ({\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ + ){{\bf W}}({\rm z}_{\rm n}^ + ,{\rm z}_{\rm m}^ + )]} ,$$ (25b)being the WRW-model of $${{\underline{\bf {X}}}^ \cap }$$ (linearity in R∩). It generates the third-order multiples in the data. In the following we will use all orders of multiples. At $${\rm z}_{\rm m}^ +$$ the feedback equations can be written as:   \begin{eqnarray} {{{\bf P}}^ - }({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}}) &=& {{\bf X}}_0^ \cup ({\rm z}_{\rm m}^ + ,{\rm z}_{\rm m}^ + ){{{\bf Q}}^ + }({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}}) \nonumber\\ &=& {{\bf X}}_0^ \cup ({\rm z}_{\rm m}^ + ,{\rm z}_{\rm m}^ + )[{{\bf W}}({\rm z}_{\rm m}^ + ,{{\rm z}_{\rm s}}){{{\bf S}}^ + }({{\rm z}_{\rm s}}) + {{{\underline{\bf {X} }}}^ \cap }({\rm z}_{\rm m}^ + ,{\rm z}_{\rm m}^ + ){{{\bf P}}^ - }({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})] \nonumber\\ &=& {{\bf X}}_0^ \cup ({\rm z}_{\rm m}^ + ,{\rm z}_{\rm m}^ + ){{\bf Q}}_0^ + ({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})[{{\bf I}} + {{{\bf A}}^ \cap }({{\rm z}_{\rm s}},{\rm z}_{\rm m}^ + ){{{\bf P}}^ - }({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})] \nonumber\\ &=& {{\bf P}}_0^ - ({{\rm z}_0};{{\rm z}_{\rm s}})[{{\bf I}} + {{{\bf A}}^ \cap }({{\rm z}_{\rm s}},{\rm z}_{\rm m}^ + ){{{\bf P}}^ - }({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})], \end{eqnarray} (26a)where (m = 1, 2, ………):   $${{\bf Q}}_0^ + ({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}}) = {{\bf W}}({\rm z}_{\rm m}^ + ,{{\rm z}_{\rm s}}){{\bf S}}_{}^ + ({{\rm z}_{\rm s}}),$$ (26b)  $${{\bf A}}_{}^ \cap ({{\rm z}_{\rm s}},{\rm z}_{\rm m}^ + ) = {[{{\bf W}}({\rm z}_{\rm m}^ + ,{{\rm z}_{\rm s}}){{{\bf S}}^ + }({{\rm z}_{\rm s}})]^{ - 1}}{{\underline{\bf {X} }}^ \cap }({\rm z}_{\rm m}^ + ,{\rm z}_{\rm m}^ + )$$ (26c)and where $${{\underline{\bf {X} }}^ \cap }({\rm z}_{\rm m}^ + ,{\rm z}_{\rm m}^ + )$$ equals the full up-down IR matrix of the overburden $$(z \le {\rm z}_{\rm m}^ + )$$ that includes all internal multiples. Similar to the expression of the down-up IR matrix $${{\underline{\bf {X} }}^ \cup } = {{{\bf R}}^ \cup } + {{{\bf X}}^ \cup }$$ at $${\rm z}_{\rm m}^ -$$, we may write for the expression of the up-down IR matrix at $${\rm z}_{\rm m}^ +$$: $${{\underline{\bf {X} }}^ \cap } = {{{\bf R}}^ \cap } + {{{\bf X}}^ \cap }$$. We will use this property to image R∩. From expression (26a) it follows that:   $${{{\bf Q}}^ + }({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}}) = {{\bf Q}}_0^ + ({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})[{{\bf I}} + {{{\bf A}}^ \cap }({{\rm z}_{\rm s}},{\rm z}_{\rm m}^ + ){{{\bf P}}^ - }({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})].$$ (26d) Figure 3. View largeDownload slide (a) Feedback model in the subsurface (m > 0), showing the generation of the internal multiples in the data. Note again the important role of multiple generation operator [I + A∩P − ] at $${\rm z}_{\rm m}^ +$$. Compare with Fig. 1(a). (b) Extension of (a), showing also the propagation operators in layer $$({\rm z}_{\rm m}^ + ,{\rm z}_{{\rm m} + 1}^ - )$$. Compare with Fig. 2. Figure 3. View largeDownload slide (a) Feedback model in the subsurface (m > 0), showing the generation of the internal multiples in the data. Note again the important role of multiple generation operator [I + A∩P − ] at $${\rm z}_{\rm m}^ +$$. Compare with Fig. 1(a). (b) Extension of (a), showing also the propagation operators in layer $$({\rm z}_{\rm m}^ + ,{\rm z}_{{\rm m} + 1}^ - )$$. Compare with Fig. 2. Note the resemblance between eqs (3a)–(3d) and (26a)–(26d). Similar to what we saw at the surface $$({\rm z}_0^ + )$$, expression (26a) shows that $$[{{\bf I}} + {{{\bf A}}^ \cap }{{\bf P}}_{}^ - ]$$ is the overburden multiple generation operator at $${\rm z}_{\rm m}^ +$$, transforming $${{\bf P}}_0^ -$$ into $${{\bf P}}_{}^ -$$ and $${{\bf Q}}_0^ +$$ into $${{\bf Q}}_{}^ +$$. From expression (26a) it also directly follows that (omitting the depth level indication in the notation):   $${{\bf P}}_0^ - = {{{\bf P}}^ - }{[{{\bf I}} + {{{\bf A}}^ \cap }{{{\bf P}}^ - }]^{ - 1}}\quad {\rm at} \quad{\rm z}_{\rm m}^ +$$ (27a)  $${{\bf Q}}_0^ + = {{{\bf Q}}^ + }{[{{\bf I}} + {{{\bf A}}^ \cap }{{{\bf P}}^ - }]^{ - 1}} \quad {\rm at} \quad {\rm z}_{\rm m}^ + ,$$ (27b)showing that [I + A∩P − ] − 1 equals the overburden multiple removal operator, transforming $${{\bf P}}_{}^ -$$ into $${{\bf P}}_0^ -$$ and $${{\bf Q}}_{}^ +$$ into $${{\bf Q}}_0^ +$$ at $${\rm z}_{\rm m}^ +$$. Considering the large similarity between the feedback model at the surface $$({\rm z}_0^ + )$$ and in the subsurface $$({\rm z}_{\rm m}^ + )$$, we can derive in exactly the same way ($${\underline{\bf {X} }}_{}^ \cap$$ plays the same role as R∩ and W + S + replaces S + ):   $${[{{\bf I}} + {{{\bf P}}^ - }{{{\bf A}}^ \cap }]^{ - 1}} = [{{\bf I}} - {{\bf P}}_0^ - {{{\bf A}}^ \cap }]\quad {\rm at} \quad {\rm z}_{\rm m}^ +$$ (28a)or   $${[{{\bf I}} - {{\bf P}}_0^ - {{{\bf A}}^ \cap }]^{ - 1}}[{{\bf I}} + {{{\bf P}}^ - }{{{\bf A}}^ \cap }] = {{\bf I}}\quad {\rm at} \quad {\rm z}_{\rm m}^ +$$ (28b)where   {{\bf A}}_{}^ \cap = $${{\bf A}}_{}^ \cap ({{\rm z}_{\rm s}},{\rm z}_{\rm m}^ + ) = {[{{\bf W}}({\rm z}_{\rm m}^ + ,{{\rm z}_{\rm s}}){{{\bf S}}^ + }({{\rm z}_{\rm s}})]^{ - 1}}{{\underline{\bf {X} }}^ \cap }({\rm z}_{\rm m}^ + ,{\rm z}_{\rm m}^ + )$$. (28c) Eq. (28a) shows that not only at the surface (m = 0) but also in the subsurface (m = 1, 2, …. .) the multiple-generating operator [I + P − A∩] has the minimum-phase property. This means that at each depth level the wavefields can be validated by testing whether they fulfil the minimum-phase condition (28b) and, optionally, the error shot records can be visualized at each depth level as a QC for the full wavefield imaging system. We consider the minimum-phase based validation and updating capability a new opportunity to improve our migration algorithms on field data in an objective way. Imaging with all multiple scattering events for perfect data If we make use of expressions (27a) and (27b) and assume perfect data (no noise, no missing data), then direct inversion leads to (m = 1, 2, …….):   $${{\bf X}}_0^ \cup ({\rm z}_{\rm m}^ + ,{\rm z}_{\rm m}^ + ) = {{{\bf P}}^ - }({\rm z}_{\rm m}^ + ;{\rm z}_{\rm s}^{}){[{{{\bf Q}}^ + }({\rm z}_{\rm m}^ + ;{\rm z}_{\rm s}^{})]^{ - 1}} = {{\bf P}}_0^ - ({\rm z}_{\rm m}^ + ;{\rm z}_{\rm s}^{}){[{{\bf Q}}_0^ + ({\rm z}_{\rm m}^ + ;{\rm z}_{\rm s}^{})]^{ - 1}}.$$ (29) Hence, we see that overburden operator [I + A∩P − ] drops out and, therefore, overburden multiple scattering (surface-related + internal) does not contribute to redatuming result $${{\bf X}}_0^ \cup ({\rm z}_{\rm m}^ + ,{\rm z}_{\rm m}^ + )$$. This means that under the ideal condition of perfect data the combination of removing all multiples from the overburden (surface and internal) followed by traditional redatuming (‘first-order redatuming’) will reveal the desired transfer operator $${{\bf X}}_0^ \cup$$. It also explains why the concept of Marchenko redatuming (Broggini & Snieder 2012) is valid in the situation of perfect data (Wapenaar et al.2014). Using expression (29) and moving inside layer $$({\rm z}_{\rm m}^ + ,{\rm z}_{{\rm m} + 1}^ - )$$ from $${\rm z}_{\rm m}^ +$$ to $${\rm z}_{{\rm m} + 1}^ -$$ by wavefield extrapolation, then we obtain (see Fig. 3b):   \begin{eqnarray} {{\bf Q}}_{}^ - ({\rm z}_{{\rm m} + 1}^ - ;{{\rm z}_{\rm s}}) &=& {{{\bf W}}^*}({\rm z}_{{\rm m} + 1}^ - ,{\rm z}_{\rm m}^ + ){{\bf P}}_{}^ - ({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}}) \nonumber\\ &=& {{{\bf W}}^*}({\rm z}_{{\rm m} + 1}^ - ,{\rm z}_{\rm m}^ + ){{\bf P}}_0^ - ({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})[{{\bf I}} + {{{\bf A}}^ \cap }({\rm z}_{\rm s}^{},{\rm z}_{\rm m}^ + ){{\bf P}}_{}^ - ({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})] \nonumber\\ &=& {{\bf Q}}_0^ - ({\rm z}_{{\rm m} + 1}^ + ;{{\rm z}_{\rm s}})[{{\bf I}} + {{{\bf A}}^ \cap }({\rm z}_{\rm s}^{},{\rm z}_{\rm m}^ + ){{\bf P}}_{}^ - ({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})], \end{eqnarray} (30a)  \begin{eqnarray} {{\bf P}}_{}^ + ({\rm z}_{{\rm m} + 1}^ - ;{{\rm z}_{\rm s}}) &=& {{\bf W}}({\rm z}_{{\rm m} + 1}^ - ,{\rm z}_{\rm m}^ + ){{\bf Q}}_{}^ + ({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}}) \nonumber\\ &=& {{\bf W}}({\rm z}_{{\rm m} + 1}^ - ,{\rm z}_{\rm m}^ + ){{\bf Q}}_0^ + ({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})[{{\bf I}} + {{{\bf A}}^ \cap }({\rm z}_{\rm s}^{},{\rm z}_{\rm m}^ + ){{\bf P}}_{}^ - ({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})] \nonumber\\ &=& {{\bf P}}_0^ + ({\rm z}_{{\rm m} + 1}^ + ;{{\rm z}_{\rm s}})[{{\bf I}} + {{{\bf A}}^ \cap }({\rm z}_{\rm s}^{},{\rm z}_{\rm m}^ + ){{\bf P}}_{}^ - ({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})]. \end{eqnarray} (30b) Using (30a) and (30b), expression (29) can be extended to:   $${{\bf {X} }}_{}^ \cup ({\rm z}_{{\rm m} + 1}^ - ,{\rm z}_{{\rm m} + 1}^ - ) = {{\bf Q}}_{}^ - ({\rm z}_{{\rm m} + 1}^ - ;{\rm z}_{\rm s}^{}){[{{\bf P}}_{}^ + ({\rm z}_{{\rm m} + 1}^ - ;{\rm z}_{\rm s}^{})]^{ - 1}} = {{\bf Q}}_0^ - ({\rm z}^-_{\rm m+1}{\rm z}_{\rm s}^{}){[{{\bf P}}_0^ + ({\rm z}_{\rm m+1}^-{\rm z}_{\rm s}^{})]^{ - 1}}.$$ (31a) Taking into account that $${{\underline{\bf {X} }}^ \cup } = {{{\bf X}}^ \cup } + {{{\bf R}}^ \cup }$$, we may write (compare with expression 20c):   $${{\bf \tilde{R}}}_{}^ \cup ({\rm z}_{{\rm m} + 1}^ - ,{\rm z}_{{\rm m} + 1}^ - ) = {{\bf F}}({\rm z}_{{\rm m} + 1}^ - ){{\bf Q}}_0^ - ({\rm z}_{{\rm m} + 1}^ - ;{\rm z}_{\rm s}^{}){[{{\bf P}}_0^ + ({\rm z}_{{\rm m} + 1}^ - ;{\rm z}_{\rm s}^{})]^{ - 1}}\quad {\rm at} \quad \tau = 0,$$ (31b)yielding the angle-dependent reflection coefficients at each gridpoint of depth level $${\rm z}_{{\rm m} + 1}^ -$$. Expression (31b) is consistent with redatuming output (29): ‘in the ideal situation of perfect data multiples drop out and they do not contribute to the image’. This means that in the situation of perfect data volumes, it is justified to follow the traditional approach and apply multiple removal prior to imaging. It also explains why the concept of direct inversion works well in the theoretical situation of perfect data (Weglein et al. 2010). Imaging with all multiple scattering events for sparse data As was already mentioned at the start, in practical situations the theoretical situation of perfect data does not exist. Similar to what we saw when dealing with surface multiple scattering, in practical situations (i) multiple scattering should not be removed but utilized and (ii) direct inversion must be replaced by minimization. Now, if we look again at expressions (21a)–(21d) and we include internal multiples, then these expressions can be extended to (see Fig. 4):   $$E({\rm z}_{{\rm m} + 1}^ - ;{{\rm z}_{\rm s}}) = \sum\limits_{\rm j} {\sum\limits_\omega {{{\left\| {[{{\bf Q}}_{}^ - ({\rm z}_{{\rm m} + 1}^ - ;{{\rm z}_{\rm s}}) - {{\bf R}}_{}^ \cup ({\rm z}_{{\rm m} + 1}^ - ,{\rm z}_{{\rm m} + 1}^ - ){{\bf P}}_{}^ + ({\rm z}_{{\rm m} + 1}^ - ;{{\rm z}_{\rm s}})]{{\vec{{\rm I}}}_{\rm j}}} \right\|}^{2}}} } = {\rm minimum},$$ (32a)where   $${{{\bf Q}}^ - }({\rm z}_{{\rm m} + 1}^ - ;{{\rm z}_{\rm s}}) = {{{\bf W}}^*}({\rm z}_{{\rm m} + 1}^ - ,{\rm z}_0^ + )\left[ {{{{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) - \sum\limits_{\ell = 1}^{\rm m} {{{\bf W}}({\rm z}_0^ + ,{\rm z}_\ell ^ - )} {{\bf R}}_{}^ \cup ({\rm z}_\ell ^ - ,{\rm z}_\ell ^ - ){{{\bf P}}^ + }({\rm z}_\ell ^ - ;{{\rm z}_{\rm s}})} \right]$$ (32b)represents the upgoing wavefields that leave depth level $${\rm z}_{{\rm m} + 1}^ -$$,   $${{\bf P}}_{}^ + ({\rm z}_{{\rm m} + 1}^ - ;{{\rm z}_{\rm s}}) = {{\bf W}}({\rm z}_{{\rm m} + 1}^ - ,{\rm z}_{\rm s}^{}){{\bf S}}_{}^ + ({{\rm z}_{\rm s}}) + {{\bf W}}({\rm z}_{{\rm m} + 1}^ - ,{\rm z}_0^ + )\delta {\underline{\bf {P} }}_{}^ + ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) + \sum\limits_{\rm n = 1}^{\rm m} {{{\bf W}}({\rm z}_{{\rm m} + 1}^ - ,{\rm z}_{\rm n}^ + )\delta {{\bf P}}_{}^ + ({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}})}$$ (32c)represents the triple downgoing wavefields that arrive at depth level $${\rm z}_{{\rm m} + 1}^ -$$ (triple illumination) and where $$\delta {\underline{\bf P}}_{}^ + ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = {{\bf R}}_{}^ \cap ({\rm z}_0^ + ,{\rm z}_0^ + ){\underline{\bf {P} }}_{}^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}})\, {\rm and}\, \delta {{\bf P}}_{\!\!\!\!\!\!-}^ + ({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}}) = {{\bf R}}_{}^ \cap ({\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ + ){{\bf P}}_{}^ - ({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}})$$ represent the secondary illuminating wavefields that leave the surface $${\rm z}_0^ +$$ and depth level $${\rm z}_{\rm n}^ +$$ (n = 1, 2, …. . m) respectively. Simplifying the notation by omitting the depth levels, minimization expression (32a) may also be written as (m = 1, 2, ……):   $${{\rm E}_{{\rm m + 1}}} = \sum\limits_{\rm j} {\sum\limits_\omega {{{\left\| {\left[ {{{\bf Q}}_{}^ - - {{\bf R}}_{\rm a}^ \cup {{{\bf W}}^ + }{{\bf S}}_{}^ + - {{\bf R}}_{\rm b}^ \cup {{{\bf W}}^ + }\delta {\underline{\bf {P} }}_{}^ + - {{\bf R}}_{\rm c}^ \cup \sum\limits_{\rm n = 1}^{\rm m} {{{{\bf W}}^ + }\delta {{\bf P}}_{}^ + } } \right]{{\vec{{\rm I}}}_{\rm j}}} \right\|}^{2}}} } = {\rm minimum\, at}\ {\rm z}_{\rm m}^ - .$$ (33) Figure 4. View largeDownload slide Wavefield propagation in the subsurface, showing the generation of multiple scattering (surface + internal) in the downgoing and upgoing wavefields during propagation. It explains how full wavefield extrapolation need be applied. In our full wavefield migration algorithm (FWM) this is implemented in a closed-loop. Bear in mind that if the source ghost is included in the surface multiple scattering, we need to take $${{{\bf R}}^ \cap }{{\underline{\bf {P} }}^ - }$$ for n = 0, where $${{\underline{\bf {P} }}^ - } = {{{\bf W}}^ - }{{{\bf S}}^ - } + {{{\bf P}}^ - }$$ at the surface. Figure 4. View largeDownload slide Wavefield propagation in the subsurface, showing the generation of multiple scattering (surface + internal) in the downgoing and upgoing wavefields during propagation. It explains how full wavefield extrapolation need be applied. In our full wavefield migration algorithm (FWM) this is implemented in a closed-loop. Bear in mind that if the source ghost is included in the surface multiple scattering, we need to take $${{{\bf R}}^ \cap }{{\underline{\bf {P} }}^ - }$$ for n = 0, where $${{\underline{\bf {P} }}^ - } = {{{\bf W}}^ - }{{{\bf S}}^ - } + {{{\bf P}}^ - }$$ at the surface. From expression (33) we see the benefit of triple reflectivity estimation by minimization: Realistic data volumes, represented by sparse data matrices, can be handled; Images can be computed per shot record, or per receiver gather, meaning that only one domain need be properly sampled; By using all multiple scattering (surface and internal) three images are obtained (‘image decomposition’), one from the primary response $$({\rm yielding }\,{{\bf R}}_{\rm a}^ \cup )$$, one from the surface multiple scattering response $$({\rm yielding }\,{{\bf R}}_{\rm b}^ \cup )$$ and one from the internal multiple scattering response $$({\rm yielding }\,{{\bf R}}_{\rm c}^ \cup )$$. Note again the significant advantage: instead of an elaborate removal process, the surface and internal multiple scattering are utilized. From the above we may conclude that image decomposition $$({{\bf R}}_{}^ \cup\,{\rm into}\,{{\bf R}}_{\rm a}^ \cup ,{{\bf R}}_{\rm b}^ \cup ,{{\bf R}}_{\rm c}^ \cup )$$ also allows us to separate primary, surface-related and internal multiple wavefields (see Fig. 4):   $$$${{\bf P}}_{\rm pr}^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = \sum\limits_{\rm m = 1}^{\rm M} {{{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm m}^ - ){{\bf R}}_{\rm a}^ \cup ({\rm z}_{\rm m}^ - ,{\rm z}_{\rm m}^ - ){{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_{\rm s}^{}){{\bf S}}_{}^ + ({{\rm z}_{\rm s}})}\quad \left( {\rm primaries} \right)$$$$ (34a)  $${{\bf M}}_{\rm sfc}^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = \sum\limits_{\rm m = 1}^{\rm M} {{{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm m}^ - ){{\bf R}}_{\rm b}^ \cup ({\rm z}_{\rm m}^ - ,{\rm z}_{\rm m}^ - ){{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_0^ + )\delta {\underline{\bf {P} }}_{}^ + ({\rm z}_0^ + ;{{\rm z}_{\rm s}})}\quad \left( {\rm surface\, multiples} \right)$$ (34b)  $${{\bf M}}_{{\mathop{\rm int}} }^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = \sum\limits_{\rm m = 2}^{\rm M} {{{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm m}^ - ){{\bf R}}_{\rm c}^ \cup ({\rm z}_{\rm m}^ - ,{\rm z}_{\rm m}^ - )\sum\limits_{\rm n = 1}^{{\rm m} - 1} {{{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_{\rm n}^ + )\delta {{\bf P}}_{}^ + ({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}})} } .\quad\left( {\rm internal\, multiples} \right).$$ (34c) These wavefields are computed in each iteration by FWMod (from depth image to time measurements, see Fig. 6). Note also that all decomposed wavefields include the missing data as well (wavefield reconstruction) and superposition leads to the full response $${{\bf P}}_{}^ - = {{\bf P}}_{\rm pr}^ - + {{\bf M}}_{\rm sfc}^ - + {{\bf M}}_{{\mathop{\rm int}} }^ -$$ (compare with expressions 23a and 23b). This is an important conclusion: using the prior information that reflectivity is an angle-dependent scalar at each depth level, ultimate wavefield separation and wavefield reconstruction should not only be part of pre-processing but it also should be an integral part of the migration algorithm. Similar to expressions (24), the combined image can be determined by a data-driven illumination- weighted addition at each subsurface gridpoint k:   $$\vec{\rm R}_{}^ \cup ({\rm k}) = [{\rm A}_{\rm k}^2\vec{\rm R}_{\rm a}^ \cup ({\rm k}) + {\rm B}_{\rm k}^2\vec{R}_{\rm b}^ \cup ({\rm k}) + {\rm C}_{\rm k}^2\vec{\rm R}_{\rm c}^ \cup ({\rm k})]/[{\rm A}_{\rm k}^2 + {\rm B}_{\rm k}^2 + {\rm C}_{\rm k}^2]$$ (35)with $${\rm A}_{\rm k}^2 = \sum\nolimits_{\rm j} {{\rm a}_{\rm kj}^2} ,{\rm B}_{\rm k}^2 = \sum\nolimits_{\rm j} {\rm b_{\rm kj}^2}\, {\rm and}\,{\rm C}_{\rm k}^2 = \sum\nolimits_{\rm j} {{\rm c}_{\rm kj}^2}$$, where akj equals the measured energy of the primary incident wavefield at gridpoint k, $${\rm b_{\rm kj}}$$ equals the measured energy of the secondary incident wavefield at gridpoint k and ckj equals the measured energy of the higher-order incident wavefield at gridpoint k, all being initially generated by primary source j. Using multiple scattering to image the lower side of boundaries as well In the foregoing we have explained how to image the upper side of boundaries by utilizing primary and multiple scattering, yielding the reflectivity operator R∪ at each depth level. If we want to image the lower side of the boundaries as well, the utilization of internal multiple scattering is a must. Having already determined R∪, we have very good knowledge of $${{\bf P}}_{\rm pr}^ -$$ and $${{\bf M}}_{\rm sfc}^ -$$ at $${\rm z}_0^ +$$ (see expressions 34a and 34b). This means that we also have very good knowledge of the internal multiples in the measurements: $${{\bf P}}_{}^ - - {{\bf P}}_{\rm pr}^ - - {{\bf M}}_{\rm sfc}^ -$$ at $${\rm z}_0^ +$$. The expression of $${{\bf M}}_{{\mathop{\rm int}} }^ -$$ is given by (see Fig. 5):   $${{\bf M}}_{{\mathop{\rm int}} }^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = \sum\limits_{\rm m = 2}^{\rm M} {\delta {{\bf M}}_{{\mathop{\rm int}} }^ - ({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}})} = \sum\limits_{\rm m = 2}^{\rm M} {{{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm m}^ - ){{\bf R}}_{\rm c}^ \cup ({\rm z}_{\rm m}^ - ,{\rm z}_{\rm m}^ - ){{\bf M}}_{{\mathop{\rm int}} }^ + ({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}})}$$ (36a)with   $${{\bf M}}_{{\mathop{\rm int}} }^ + ({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}) = \sum\limits_{\rm n = 1}^{{\rm M} - 1} {{{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_{\rm n}^ + ){{\bf R}}_{}^ \cap ({\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ + ){{\bf P}}_{}^ - ({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}})}$$ (36b)or, changing the summation sequence,   $${{\bf M}}_{{\mathop{\rm int}} }^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = \sum\limits_{\rm n = 1}^{{\rm M} - 1} {{{\bf X}}_{\rm pr}^ \cup ({\rm z}_0^ + ,{\rm z}_{\rm n}^ + ){{\bf R}}_{}^ \cap ({\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ + ){{\bf P}}_{}^ - ({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}})}$$ (36c)with   $${{\bf X}}_{\rm pr}^ - ({\rm z}_0^ + ;{\rm z}_{\rm n}^ + ) = \sum\limits_{\rm m = n + 1}^{\rm M} {{{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm m}^ - ){{\bf R}}_{\rm c}^ \cup ({\rm z}_{\rm m}^ - ,{\rm z}_{\rm m}^ - ){{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_{\rm n}^ - )} .$$ (36d) Figure 5. View largeDownload slide Grouping of the internal multiples $$(\delta {{\bf M}}_{{\mathop{\rm int}} }^ - )$$ according to generation by upper-side reflectivity $${{{\bf R}}^ \cup }({\rm z}_{\rm m}^ - ,{\rm z}_{\rm m}^ - )$$, see left-hand side, and by lower-side reflectivity $${{{\bf R}}^ \cap }({\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ + )$$, see right-hand side. Figure 5. View largeDownload slide Grouping of the internal multiples $$(\delta {{\bf M}}_{{\mathop{\rm int}} }^ - )$$ according to generation by upper-side reflectivity $${{{\bf R}}^ \cup }({\rm z}_{\rm m}^ - ,{\rm z}_{\rm m}^ - )$$, see left-hand side, and by lower-side reflectivity $${{{\bf R}}^ \cap }({\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ + )$$, see right-hand side. We will use expressions (36c) and (36d) for the estimation of lower-side reflectivity $${{{\bf R}}^ \cap }({\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ + )$$, n = 1, 2, ………M − 1. If we substitute the three wavefield components of incident wavefield $${{\bf P}}_{}^ -$$ in (36c), we obtain:   $${{\bf M}}_{{\mathop{\rm int}} }^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = \sum\limits_{\rm n = 1}^{{\rm M} - 1} {{{\bf X}}_{\rm pr}^ \cup ({\rm z}_0^ + ,{\rm z}_{\rm n}^ + ){{\bf R}}_{}^ \cap ({\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ + )[{{\bf P}}_{\rm pr}^ - ({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}})} + {{\bf M}}_{\rm sfc}^ - ({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}}) + {{\bf M}}_{{\mathop{\rm int}} }^ - ({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}})]$$ (37a)where wavefields $${{\bf P}}_{\rm pr}^ -$$, $${{\bf M}}_{\rm sfc}^ -$$ and $${{\bf M}}_{{\mathop{\rm int}} }^ -$$ at $${\rm z}_{\rm n}^ +$$ are known and computed by FWMod:   $${{\bf P}}_{\rm pr}^ - ({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}}) = \sum\limits_{\ell = {\rm n + 1}}^{\rm M} {{{\bf W}}({\rm z}_{\rm n}^ + ,{\rm z}_\ell ^ - ){{\bf R}}_{\rm a}^ \cup ({\rm z}_\ell ^ - ,{\rm z}_\ell ^ - ){{\bf W}}({\rm z}_\ell ^ - ,{\rm z}_{\rm s}^{}){{\bf S}}_{}^ + ({{\rm z}_{\rm s}})}$$ (37b)  $${{\bf M}}_{\rm sfc}^ - ({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}}) = \sum\limits_{\ell = {\rm n + 1}}^{\rm M} {{{\bf W}}({\rm z}_{\rm n}^ + ,{\rm z}_\ell ^ - ){{\bf R}}_{\rm b}^ \cup ({\rm z}_\ell ^ - ,{\rm z}_\ell ^ - ){{\bf W}}({\rm z}_\ell ^ - ,{\rm z}_0^{})\delta {{\bf {P} }}_{}^ + ({\rm z}_0^{};{{\rm z}_{\rm s}})}$$ (37c)  $${{\bf M}}_{{\mathop{\rm int}} }^ - ({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}}) = \sum\limits_{\ell = {\rm n + 1}}^{\rm M} {{{\bf W}}({\rm z}_{\rm n}^ + ,{\rm z}_\ell ^ - ){{\bf R}}_{\rm c}^ \cup ({\rm z}_\ell ^ - ,{\rm z}_\ell ^ - )\sum\limits_{\ell = 1}^{{\rm m} - 1} {{{\bf W}}({\rm z}_\ell ^ - ,{\rm z}_{\rm n}^ + )\delta {{\bf P}}_{}^ + ({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}})} } .$$ (37d) Expressions (37a)–(37d) lead to minimization equation (omitting the depth levels):   $$\sum\limits_{\rm j} {\sum\limits_\omega {{{\left\| {\left[ {[{{\bf P}}_{}^ - - {{\bf P}}_{\rm pr}^ - - {{\bf M}}_{\rm sfc}^ - ] - {{\bf X}}_{\rm pr}^ \cup {{\bf R}}_{\rm a}^ \cap {{\bf P}}_{\rm pr}^ - - {{\bf X}}_{\rm pr}^ \cup {{\bf R}}_{\rm b}^ \cap {{\bf M}}_{\rm sfc}^ - - {{\bf X}}_{\rm pr}^ \cup {{\bf R}}_{\rm c}^ \cap {{\bf M}}_{{\mathop{\rm int}} }^ - } \right]{{\vec{{\rm I}}}_{\rm j}}} \right\|}^2}} } = {\rm minimum\, at}\,{\rm z}_{\rm n}^ + .$$ (38) Compare with minimization (33). Note that in minimization (38) the measurements in $${{\bf P}}_{}^ -$$ are known, the reflectivity operators $${{\bf R}}_{\rm a}^ \cap$$, $${{\bf R}}_{\rm b}^ \cap$$ and $${{\bf R}}_{\rm c}^ \cap$$ are unknown and $${{\bf P}}_{\rm pr}^ -$$, $${{\bf M}}_{\rm sfc}^ -$$, $${{\bf M}}_{{\mathop{\rm int}} }^ -$$ are given by expressions (37b)–(37d). It is important to realize that R∩ = −R∪ in areas where wave conversion can be neglected (see Appendix  B). This is the start of iterative minimization (38). Closed-loop, full-wavefield imaging scheme for sparse data Fig. 6 shows the schematic of closed-loop, full wavefield migration (FWM). Bear in mind that in many practical situations already some kind of estimate of the image is available (think not only in terms of previous time lapse images but think also at output already available from mainstream seismic processing); this knowledge should be optionally used as input in the first iteration of FWM. In the feedback loop full wavefield modeling (FWMod) transforms the (initial) image into full wavefield measurements (‘image transformation’). These measurements are presented in terms of primaries, surface multiples and internal multiples according to expressions (34a)–(34c). This (initial) wavefield decomposition is used in the construction of three separate images (image decomposition): first-order image, second-order image and higher-order image. This is input for the next iteration. As a last step, these three images may be combined into one final image according to expression (35). In the downgoing path the reflectivities are computed. In the upgoing path the different images are transformed into different wavefields. These wavefields are compared with the measured wavefields, allowing (i) a scaling correction of the internal wavefields, (ii) filling in the missing data, (iii) estimating the direct and ghost source properties and (iv) updating the surface reflectivity by refining the detector deghosting process (marine data). Roundtrips are continued until simulated wavefields and measured wavefields are sufficiently close (Berkhout 2014). Figure 6. View largeDownload slide In closed-loop, full-wavefield migration (FWM) the seismic image is transformed into simulated measurements (FWMod) and they are compared with the real measurements. After minimization the error measurements are input for the next iteration. Decomposed wavefields and decomposed images are essential constituents in the algorithm. Figure 6. View largeDownload slide In closed-loop, full-wavefield migration (FWM) the seismic image is transformed into simulated measurements (FWMod) and they are compared with the real measurements. After minimization the error measurements are input for the next iteration. Decomposed wavefields and decomposed images are essential constituents in the algorithm. Numerical illustration To illustrate the closed-loop, full wavefield imaging theory we consider a simple subsurface model with three angle-dependent reflectors. The simplicity of the model allows us to identify all primary and multiple-scattering wavefields as they occur in the theory. Fig. 7 shows the model with the angle-dependent reflectivity operators (R∪ at $${\rm z}_1^ - ,{\rm z}_2^ - ,{\rm z}_3^ -$$ and R∩ at $${\rm z}_1^ + ,{\rm z}_2^ +$$). Figure 7. View largeDownload slide Subsurface model with three angle-dependent reflectors. The reflection operators R∪ and R∩ for each reflector are displayed within the angle range of the seismic data, the dashed line representing the imaginary values (which are estimated as well). The seven source locations in this experiment are indicated. Detectors are situated over the full range (9000 m) with a sampling interval of 15 m. Figure 7. View largeDownload slide Subsurface model with three angle-dependent reflectors. The reflection operators R∪ and R∩ for each reflector are displayed within the angle range of the seismic data, the dashed line representing the imaginary values (which are estimated as well). The seven source locations in this experiment are indicated. Detectors are situated over the full range (9000 m) with a sampling interval of 15 m. Figs 8(a)–(c) visualize the wavefields that play a principal role in the today's mainstream seismic processing packages: Fig. 8(a) shows the downgoing primary and secondary source wavefields that leave the surface $${\rm z}_0^ +$$, being one column (j = 4) of   $${{{\bf Q}}^ + } = {{{\bf S}}^ + } + {{{\bf R}}^ \cap }{{{\bf P}}^ - },$$ (39a)where S + represents the impulsive primary source wavefields and R∩P − represents the highly dispersive secondary source wavefields. It is assumed that the ghost wavefields have already been removed. Fig. 8(b) shows the upgoing reflected wavefields (PP) that arrive at the surface $${\rm z}_0^ +$$, being one column (j = 4) of   $${{\bf P}}_{}^ - = {{\bf P}}_0^ - + {{\bf M}}_0^ - ,$$ (39b)where $${{\bf P}}_0^ -$$ represents the surface-free response $$({{\bf X}}_0^ \cup {{{\bf S}}^ + })$$ and $${{\bf M}}_0^ -$$ represents the surface-related multiples $$({{\bf X}}_0^ \cup {{{\bf R}}^ \cap }{{{\bf P}}^ - })$$. The separation is realized by a surface-related multiple removal algorithm. Fig. 8(c) shows the arriving wavefields (PP) at $${\rm z}_0^ +$$ that are generated by S + only, being one column (j = 4) of   $${{\bf P}}_0^ - = {{\bf P}}_{\rm pr}^ - + {\mathop{\bf M}^{\frown }} {}_{{\rm int}}^ - ,$$ (39c)where $${{\bf P}}_{\rm pr}^ -$$ represents the primary response $$({{\bf X}}_{\rm pr}^ \cup {{{\bf S}}^ + })$$ and $$\displaystyle{\mathop{\bf M}^{\frown }} {}_{{\rm int}}^ -$$ represents the internal multiples in $${{\bf P}}_0^ -$$. The separation is realized by an internal multiple removal algorithm. Figure 8. View largeDownload slide Wavefields at $${\rm z}_0^ +$$ (generated by the central source j = 4) that play a principal role in today's mainstream seismic processing packages. $$\displaystyle{\mathop{\bf M}^{\frown }}{}_{{\rm int}}^ - {\vec{\rm I}_{\rm j}}$$ contains the internal multiples being generated by primary source $${{\bf S}}_0^ + {\vec{\rm I}_{\rm j}}$$. The internal multiples being generated by secondary source $${{{\bf R}}^ \cap }{{{\bf P}}^ - }{\vec{\rm I}_{\rm j}}$$ are part of the surface-related multiples $${{\bf M}}_0^ - {\vec{\rm I}_{\rm j}}$$. Figure 8. View largeDownload slide Wavefields at $${\rm z}_0^ +$$ (generated by the central source j = 4) that play a principal role in today's mainstream seismic processing packages. $$\displaystyle{\mathop{\bf M}^{\frown }}{}_{{\rm int}}^ - {\vec{\rm I}_{\rm j}}$$ contains the internal multiples being generated by primary source $${{\bf S}}_0^ + {\vec{\rm I}_{\rm j}}$$. The internal multiples being generated by secondary source $${{{\bf R}}^ \cap }{{{\bf P}}^ - }{\vec{\rm I}_{\rm j}}$$ are part of the surface-related multiples $${{\bf M}}_0^ - {\vec{\rm I}_{\rm j}}$$. The primary wavefield $${{\bf P}}_{\rm pr}^ -$$ is used as input in first-order migration (today's standard). $$\displaystyle{\mathop{\bf M}^{\frown }}{}_{{\rm int}}^ -$$ represents the group of internal multiples that is generated by the primary sources (columns of S + ). The other group of internal multiples $$(\displaystyle{\mathop{\bf M}^{\frown^{\!\!\!\!\!\!\frown} }}{}_{{\mathop{\rm int}} }^ - )$$ can be found in $${{\bf M}}_0^ -$$ and is generated by the secondary sources (columns of R∩P − ) at the surface: $$\displaystyle{\mathop{\bf M}^{\frown^{\!\!\!\!\!\!\frown}}}{}_{{\mathop{\rm int}} }^ - = \displaystyle{\mathop{\bf M}^{\frown}}{}_{{\mathop{\rm int}} }^ - {{{\bf A}}^ \cap }{{\bf P}}_{}^ -$$ with $${{{\bf A}}^ \cap } = {[{{{\bf S}}^ + }]^{ - 1}}{{\bf R}}_{}^ \cap$$. It is interesting to see that $$\displaystyle{\mathop{\bf M}^{\frown}}{}_{{\mathop{\rm int}} }^ -$$ is the prediction operator that determines $$\displaystyle{\mathop{\bf M}^{\frown^{\!\!\!\!\!\!\frown} }}{}_{{\mathop{\rm int}} }^ -$$ from the input data $${{\bf P}}_{}^ -$$. Note here the complete similarity with $${{\bf P}}_0^ -$$ functioning as the prediction operator for the determination of surface-related multiples $${{\bf M}}_0^ -$$: $${{\bf M}}_0^ - = {{\bf P}}_0^ - {{{\bf A}}^ \cap }{{\bf P}}_{}^ -$$. Next, Figs 9(a)–(c) visualize the wavefields that play a principal role in closed-loop, multiple-scattering imaging: Fig. 9(a) shows the downgoing primary and secondary source wavefields that leave the surface $${\rm z}_0^ +$$, being one column (j = 4) of   $${{{\bf Q}}^ + } = {{{\bf W}}^*}{{{\bf S}}^ + } + {{{\bf R}}^ \cap }({{{\bf W}}^ - }{{{\bf S}}^ - } + {{{\bf P}}^ - }),$$ (40a)where the primary sources $${{\bf S}}_{}^ +$$ are situated at $${\rm z}_{\rm s}^{} = 50\,{\rm m}$$, where $${{\bf W}}_{}^*$$ transforms the downgoing wavefield of these sources to $${\rm z}_0^ +$$ and where the ghost sources $${{{\bf R}}^ \cap }({{\bf W}}_{}^ - {{\bf S}}_{}^ - )$$ are included in the surface multiple scattering $$({{{\bf R}}^ \cap }{\underline{\bf {P} }}^ - )$$. Fig. 9(b) shows the upgoing reflected wavefields (PP) that arrive at the surface $${\rm z}_0^ +$$, being one column (j = 4) of   $${{\bf P}}_{}^ - = {{\bf P}}_{\rm pr}^ - + {{\bf M}}_{}^ - ,$$ (40b)where $${{\bf P}}_{\rm pr}^ - = {{\bf X}}_{\rm pr}^ \cup {{{\bf S}}^ + }$$ and $${{\bf M}}_{}^ -$$ contains all multiples (surface +internal). Fig. 9(c) shows the upgoing surface- and internal-multiple scattering wavefields that arrive at $${\rm z}_0^ +$$, being one column (j = 4) of   $${{\bf M}}_{}^ - = {{\bf M}}_{\rm sfc}^ - + {{\bf M}}_{{\mathop{\rm int}} }^ - ,$$ (40c)where $${{\bf M}}_{\rm sfc}^ - = {{\bf X}}_{\rm pr}^ \cup ({{{\bf R}}^ \cap }{{\bf {P} }}_{\!\!\!\!\!\!-}^ - )$$ contains all surface multiples and $${{\bf M}}_{{\mathop{\rm int}} }^ -$$ contains all internal multiples: $$( \displaystyle{\mathop{\bf M}^{\frown}}{}_{{\mathop{\rm int}} }^ - + \displaystyle{\mathop{\bf M}^{\frown}}{}_{{\mathop{\rm int}} }^ - {{{\bf A}}^ \cap }{{\bf {P} }}_{\!\!\!\!\!\!-}^ - )$$. Figure 9. View largeDownload slide Wavefields at $${\rm z}_0^ +$$ (generated by the central source j = 4) that play a principal role in closed-loop, full wavefield imaging. Note that $${{\bf M}}_{}^ - {\vec{\rm I}_{\rm j}}$$ contains all multiples (surface + internal), $${{\bf M}}_{\rm sfc}^ - {\vec{\rm I}_{\rm j}}$$ contains only the surface multiples and $${{\bf M}}_{{\mathop{\rm int}} }^ - {\vec{\rm I}_{\rm j}}$$ contains all internal multiples $$(\displaystyle{\mathop{\bf M}^{\frown}}{}_{{\mathop{\rm int}}}^ - {\vec{\rm I}_{\rm j}} + \displaystyle{\mathop{\bf M}^{\frown^{\!\!\!\!\!\!\frown} }}{}_{{\mathop{\rm int}}}^ - {\vec{\rm I}_{\rm j}})$$. Figure 9. View largeDownload slide Wavefields at $${\rm z}_0^ +$$ (generated by the central source j = 4) that play a principal role in closed-loop, full wavefield imaging. Note that $${{\bf M}}_{}^ - {\vec{\rm I}_{\rm j}}$$ contains all multiples (surface + internal), $${{\bf M}}_{\rm sfc}^ - {\vec{\rm I}_{\rm j}}$$ contains only the surface multiples and $${{\bf M}}_{{\mathop{\rm int}} }^ - {\vec{\rm I}_{\rm j}}$$ contains all internal multiples $$(\displaystyle{\mathop{\bf M}^{\frown}}{}_{{\mathop{\rm int}}}^ - {\vec{\rm I}_{\rm j}} + \displaystyle{\mathop{\bf M}^{\frown^{\!\!\!\!\!\!\frown} }}{}_{{\mathop{\rm int}}}^ - {\vec{\rm I}_{\rm j}})$$. The three response wavefields $${{\bf P}}_{\rm pr}^ - ,{{\bf M}}_{\rm sfc}^ - ,{{\bf M}}_{{\mathop{\rm int}} }^ -$$ are used as input to the closed-loop, multiple-scattering algorithm. Bear in mind that the surface multiples $$({{\bf M}}_{\rm sfc}^ - )$$ do not contain internal multiples, but they do contain the ghost source response ($${{\bf M}}_{\rm sfc}^ - = {{\bf P}}_{\rm pr}^ - {{{\bf A}}^ \cap }{{\bf {P} }}_{\!\!\!\!\!\!-}^ -$$). Next, Figs 10(a)–(c) visualize the decomposed wavefields that enter and leave the upper part of the third reflector, followed by the decomposed images in the ray parameter domain: Fig. 10(a) shows one column (j = 4) of the incident wavefield $${{\bf W}}_{}^ + {{\bf S}}_{}^ +$$ and its response $${{\bf P}}_{\rm pr}^ -$$ at $${\rm z}_3^ -$$. Using all seven sources, first-order image $${{\bf \tilde{R}}}_{\rm a}^ \cup$$ is obtained at $${\rm z}_3^ -$$. Note the acquisition shadow zones in the first-order image due to the coarse sampling of the man-made primary sources (think of a cross section). Fig. 10(b) shows one column (j = 4) of the incident wavefield $${{\bf W}}_{}^ + ({{{\bf R}}^ \cap }{{\bf P}}_{}^ - )$$ and its response $${{\bf M}}_{\rm sfc}^ -$$ at $${\rm z}_3^ -$$. Using all seven sources, second-order image $${{\bf \tilde{R}}}_{\rm b}^ \cup$$ is obtained at $${\rm z}_3^ -$$. Note that in the second-order image the shadow zones have largely disappeared and the illumination area has increased due to the denser and wider illumination by the nature-made secondary sources at the surface. Fig. 10(c) shows one column (j = 4) of the incident wavefield $${{\bf W}}_{}^ + \sum\nolimits_{\rm n = 1}^2 {{{{\bf R}}^ \cap }{{\bf P}}_{}^ - }$$ and its response $${{\bf M}}_{{\mathop{\rm int}} }^ -$$ at $${\rm z}_3^ -$$. Using all seven sources, high-order image $${{\bf \tilde{R}}}_c^ \cup$$ is obtained at $${\rm z}_3^ -$$. Similar to the surface multiples, note the extra contribution of the internal multiples with respect to the primaries. Figure 10. View largeDownload slide Decomposed incident wavefields (left-hand side) with their responses (in the middle) at the upper part of the third reflector, being generated by the central source j = 4. By using all seven sources, the decomposed angle-dependent reflectivities are obtained (right-hand side). Note that the reflectivities are displayed in the ray-parameter domain for each spatial location along the reflector. Figure 10. View largeDownload slide Decomposed incident wavefields (left-hand side) with their responses (in the middle) at the upper part of the third reflector, being generated by the central source j = 4. By using all seven sources, the decomposed angle-dependent reflectivities are obtained (right-hand side). Note that the reflectivities are displayed in the ray-parameter domain for each spatial location along the reflector. In this illustration multiple scattering from dipping boundaries does not occur, minimizing thus the transformation from low to high angles and vice versa. However, we do already see the increase in angle range and, last but not least, we also see that multiples reveal the angle-dependent reflection information. First results from complex examples confirm the expectation that the more complex the geology, the more valuable the contribution from multiple scattering (current research). Finally, Figs 11(a)–(c) visualize the decomposed wavefields that enter and leave the lower part of the second reflector, followed by the decomposed images in the ray parameter domain: Fig. 11(a) shows one column (j = 4) of the incident wavefield $${{\bf P}}_{\rm pr}^ -$$ and its response $$\delta {{\bf P}}_{\rm pr}^ + = {{\bf R}}_{\rm a}^ \cap {{\bf P}}_{\rm pr}^ -$$ at $${\rm z}_2^ +$$. Using all seven sources leads to the angle-dependent image $${{\bf \tilde{R}}}_a^ \cap$$ at $${\rm z}_2^ +$$. Fig. 11(b) shows one column (j = 4) of the incident wavefield $${{\bf M}}_{\rm sfc}^ -$$ and its response $$\delta {{\bf M}}_{\rm sfc}^ + = {{\bf R}}_{\rm b}^ \cap {{\bf M}}_{\rm sfc}^ -$$ at $${\rm z}_2^ +$$. Using all seven sources leads to the angle-dependent image $${{\bf \tilde{R}}}_b^ \cap$$ at $${\rm z}_2^ +$$. Fig. 11(c) shows one column (j = 4) of the incident wavefield $${{\bf M}}_{{\mathop{\rm int}} }^ -$$ and its response $$\delta {{\bf M}}_{{\mathop{\rm int}} }^ + = {{\bf R}}_{\rm c}^ \cap {{\bf M}}_{{\mathop{\rm int}} }^ -$$ at $${\rm z}_2^ +$$. Using all seven sources leads to the angle-dependent image $${{\bf \tilde{R}}}_c^ \cap$$ at $${\rm z}_2^ +$$. Figure 11. View largeDownload slide Decomposed incident wavefields (left-hand side) with their responses (in the middle) at the lower part of the second reflector, being generated by the central source j = 4. By using all seven sources, the decomposed angle-dependent reflectivities are obtained (right-hand side). Note again that the reflectivities are displayed in the ray-parameter domain for each spatial location along the reflector. Figure 11. View largeDownload slide Decomposed incident wavefields (left-hand side) with their responses (in the middle) at the lower part of the second reflector, being generated by the central source j = 4. By using all seven sources, the decomposed angle-dependent reflectivities are obtained (right-hand side). Note again that the reflectivities are displayed in the ray-parameter domain for each spatial location along the reflector. Of course, reflection images of the lower part of reflectors can only be obtained by utilizing internal multiple scattering (see also Fig. 5). Note that in some areas, for instance think of shallow reflectors in the unconsolidated near-surface layer on land, images of the upper part are extremely difficult to acquire in practice. Hence, for these situations new opportunities arise if we aim at images of the lower part ($${{\bf R}}_{}^ \cap$$ instead of $${{\bf R}}_{}^ \cup$$) by using illumination with wavefields that arrive from below: $${{\bf P}}_{}^ - ({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}}) = \sum\limits_{\rm n = m + 1}^{\rm M} {{{\bf W}}({\rm z}_{\rm m}^ + ,{\rm z}_{\rm n}^ - ){{\bf R}}_{\rm c}^ \cup ({\rm z}_{\rm n}^ - ,{\rm z}_{\rm n}^ - ){{\bf P}}_{}^ + ({\rm z}_{\rm n}^ - ;{{\rm z}_{\rm s}})}$$. But also in borehole seismic imaging new opportunities are created. This example illustrates that multiple scattering is not only important to strengthen the response from the upper part of geological boundaries, it is our only hope to image those boundaries which responses from above are not measured at all. CONCLUSIONS In the theoretical situation of perfect data (no noise, no missing data) we have top-quality, fully filled seismic data matrices and primary sources are able to realize complete illumination, even inside geological shadow zones. For such perfect data volumes direct inversion algorithms can be applied to remove surface-related and internal multiple scattering. Using this multiple-free result as input, first-order redatuming and first-order primary migration can provide us with the full image. In practice, however, seismic measurements always contain noise and we never have complete 3-D data volumes at our disposal. We actually deal with sparse data matrices that cause acquisition-generated (non-geological) shadow zones in the subsurface. Without the introduction of assumptions, real data matrices cannot be used in direct inversion and redatuming. Therefore, in practice we badly need to include multiple-scattering wavefields to fill up the first-order illumination gaps that are caused by the sparse measurements—possibly in combination with geological shadow zones—and the background noise. Note that multiple-scattering contributions are ‘free of charge’ and also include the ghost source as extra illumination (when applicable). The author's message is that in the real world multiple scattering is a blessing; it must not be removed but it must be utilized. Instead of applying a complex and elaborate multiple scattering removal process, both surface and internal multiple scattering can be utilized to give an extra contribution to the illuminating wavefield (‘triple illumination’). Internal multiple scattering even allows us to obtain a reflection image from the lower part of reflectors! This is explained by the closed-loop, full-wavefield imaging theory. The full-wavefield extension of the imaging theory also explains that decomposition of the seismic response in first-order (‘primaries’), second-order (‘surface multiples’) and higher-order (‘internal multiples’) scattering should be an integral part of the migration algorithm. In closed-loop FWM the concepts of image decomposition and wavefield decomposition enhance each other (Fig. 6). By comparing the decomposed wavefield components (34a)–(34c) with the measured data, the estimates of the direct and ghost source properties (S + , δS + ) and the surface reflectivity (R∩) can be updated. In addition, the estimates of the missing data can be improved by wavefield reconstruction per decomposed wavefield (part of FWMod). All these capabilities tell us that the closed-loop, full-wavefield migration algorithm should be considered as the nucleus of the future seismic processing system. Ultimately it will position pre-processing as part of an ‘iteration zero’ activity. By preserving the multiple scattering in the seismic measurements, we also preserve a fundamental physical property—the minimum-phase property—in the measurements. This means that an objective way of testing the validity of the estimated wavefields at each depth level can be carried out. It is proposed as an integral part of the closed-loop, full-wavefield imaging process. Epilogue: an image of the future In the history of seismic imaging the presence of strong multiple scattering wavefields has always been perceived as bad news. This explains why a large amount of effort was spent to develop smart algorithms for the attenuation and removal of the ‘damned’ multiples. The author was part of this painstaking endeavour for many years. Today we start to realize that multiple scattering wavefields provide us with extra illumination. There exist billions of different multiples, reaching areas where primaries may have disappeared below the noise (‘shadow zones’) or are not present at all (from the lower part of boundaries). In addition, multiples contain extra information on velocity and absorption, as they travel not only once but many times through the subsurface. The challenge is not to attenuate or remove them, but to find out how to extract subsurface information from these extremely complex phenomena. This article aims at changing the traditional aversion in the seismic industry against multiple scattering by convincing the reader that the presence of multiples means actually good news. This has been done by showing how the conventional theory of second- and higher-order seismic wavefields can be reformulated such that the result leads to a transparent theoretical framework (with simple WRW-modules) and a practical imaging algorithm (with decomposed images). The output consists of angle-dependent reflectivity images that honour wavefield conversion (no acoustic assumption) and include transmission effects (Appendix  B). One important aspect is to organize multiples not according to the well-known sequences of down- and upward traveling (ray)paths, but according to downward-scattered wavefields (R∩P − ) at each depth level (m = 0, 1, 2, …….). This R∩-related organization leads to new theoretical insights that tell us that the very complex multiple scattering generating systems in the Earth have the minimum-phase property. This means that under perfect data conditions multiple-scattering removal systems are causal, allowing the prediction of multiple-scattering from measurements of the past only (in multiple scattering the past determines the future). The minimum-phase property makes clear why algorithms such as direct inversion and Marchenko redatuming can perform well under perfect data conditions (complete sampling, noise-free). Today, seismic data has been recorded in most areas of the world and seismic images are widely available. This situation certainly applies to the many mature exploration and production basins. The multiple-scattering imaging algorithm welcomes such an initial image. It is directly used in the feedback loop of Fig. 6 to make a first estimate of decomposed wavefields down $$({{{\bf W}}^ + }{{\bf S}}_{}^ + ,{{{\bf W}}^ + }\delta {{\bf {P} }}_{\!\!\!\!\!\!-}^ + ,{{\bf M}}_{{\mathop{\rm int}} }^ + )$$ and up $$({{\bf P}}_{\rm pr}^ - ,{{\bf M}}_{\rm sfc}^ - ,{{\bf M}}_{{\mathop{\rm int}} }^ - )$$. Of course, in the special situation that no initial seismic image is available yet, current mainstream processing—with data interpolation, deblending, multiple removal and first-order migration as the main steps—has still the important task to estimate an initial image (‘iteration zero’ in closed-loop full-wavefield imaging). Following this way of thinking, current mainstream processing in the industry would shift to the early exploration phase, thus becoming a pre-imaging tool in virginal areas. According to the author, removal of multiple scattering is more difficult than imaging of multiple scattering. Or, in other words, removal of multiple scattering should not be done before imaging but it should be done as part of imaging. Removal before imaging means prediction and subtraction of multiples in the real-data time domain. The slightest phase error in the prediction result leads already to a large residue, particularly in the higher part of the seismic spectrum. Moreover the data-domain prediction process is very complex, as the prediction operator needs many low-noise shot records that are well sampled in both the source and detector domain. On the other hand, in multiple-scattering imaging the subsurface information in an initial image is utilized to estimate the decomposed wavefields (down as well as up) by image transformation (FWMod). These decomposed wavefields provide ‘triple illumination’ and are again used to construct new decomposed images, etc. (see the closed-loop scheme in Fig. 6). Bear in mind that imaging is a massive double stacking process (focusing in emission + focusing in detection) that is known to be forgiving for initial errors and robust in improving signal-to-noise ratio. Decomposition of wavefields and images may provide new opportunities in the design of (blended) data acquisition geometries. The contribution of all different wavefields to the final image is very nonlinear in reflectivity and may change the outcome of the optimum source and detector distributions. For instance, our current research indicates that the traditional symmetric sampling concept may need be replaced by an asymmetric one where only one domain need be well sampled. This appears even more valid in the situation of complex geology, where limited-angle, first-order wavefields are transformed to high-order wavefields that create incidence from all sides. Finally, I remark about velocities. Traditionally, migration velocities are estimated with the aid of primary wavefields. However, higher-order reflections are more sensitive to velocity errors than primaries. This means that there exist opportunities to extend velocity estimation to multiple scattering wavefields for more accurate velocities. We are working on full-automatic, anisotropic velocity estimation for multiple-scattering imaging. In this ambitious effort, estimation and analysis of propagation operator W plays a key role. Ultimately, we have the ambition to apply seismic data acquisition and multiple-scattering imaging in exactly the same way as self-driving cars collect data and simultaneously image their immediate environment. This approach may revolutionize the way we acquire and image time laps data by considering the sequence of surveys as one big data volume that defines the image of a geological medium with unknown time-variant components. Acknowledgements The author is grateful to the Delphi sponsors for their continuing financial support and for the stimulating discussions on industry needs during the Delphi meetings. In addition, the author would like to thank Dr Eric Verschuur who supported him with the generation of the numerical illustrations. REFERENCES Araujo F.V., Weglein A.B., Carvalho P.M., Stolt R.H., 1994. Inverse scattering series for multiple attenuation: an example with surface and internal multiples, in SEG Expanded Abstracts , pp. 704– 706. Berkhout A.J., 1971. Minimum Phase in Sampled-Signal Theory , TU Delft. Berkhout A.J., 1974. Related properties of minimum-phase and zero-phase time functions, Geophys. Prospect. , 22, 683– 709. https://doi.org/10.1111/j.1365-2478.1974.tb00111.x Google Scholar CrossRef Search ADS   Berkhout A.J., 1982. Seismic Migration: Imaging of Acoustic Energy by Wavefield Extrapolation, A: Theoretical Aspects , 2nd edn, Elsevier. Berkhout A.J., 2014. Review paper: an outlook on the future of seismic imaging, Part II: Full-wavefield migration, Geophys. Prospect. , 62( 5), 931– 949. https://doi.org/10.1111/1365-2478.12154 Google Scholar CrossRef Search ADS   Broggini F., Snieder R., 2012. Connection of scattering principles: a visual and mathematical tour, Eur. J. Phys. , 33, 593– 613. https://doi.org/10.1088/0143-0807/33/3/593 Google Scholar CrossRef Search ADS   Davydenko M., Verschuur D.J., van Groenestijn G.J.A., 2015. Full wavefield migration applied to field data, in 77th Annual International Conference and Exhibition, EAEG, Extended Abstracts , We-N101-10. Verschuur D.J., Berkhout A.J., 1997. Estimation of multiple scattering by iterative inversion, Part II: Practical aspects and examples, Geophysics , 62, 1596– 1611. https://doi.org/10.1190/1.1444262 Google Scholar CrossRef Search ADS   Wapenaar K., Thorbecke J., van der Neut J., Broggini F., Slob E., Snieder R., 2014. Marchenko imaging, Geophysics , 79, WA39– WA57. https://doi.org/10.1190/geo2013-0302.1 Google Scholar CrossRef Search ADS   Weglein A.B., Liu F., Wang Z., Li X., Liang H., 2010. The inverse scattering series depth imaging algorithms: development, tests and progress towards field data application, in SEG, Soc. Expl. Geophys., Expanded abstracts , Denver, pp. 4133– 4138. Weglein A.B., 2014. Multiples: signal or noise?, in SEG Technical Program Expanded Abstracts , pp. 4393– 4399. Weglein A.B., 2015. Primaries: the only events that can be migrated and for which migration has meaning, Leading Edge , 33, 808– 813. https://doi.org/10.1190/tle34070808.1 Google Scholar CrossRef Search ADS   Whitmore N.D., Valenciano A.A., Sollner W., 2010. Imaging of primaries and multiples, using a dual-sensor towed streamer, in 80th Ann. Internat. Mtg, SEG, Expanded Abstracts , pp. 3187– 3192. APPENDIX A: MULTIPLE REMOVAL WITH THE MINIMUM-PHASE CONDITION Using main text expressions (12a) and (12b), we may conclude that at $${\rm z}_0^ +$$:   $$(1)\,\,{{{\bf E}}_0} = [{{\bf I}} - {{\bf \hat{Y}}}_0^{}][{{\bf I}} + {{\bf Y}}] - {{\bf I}} = {{\bf Y}} - {{\bf \hat{Y}}}_0^{} - {{\bf \hat{Y}}}_0^{}{{\bf Y}}$$ (A1a)  $$(2)\,\,\Delta {{\bf Y}}_0^{} = {{{\bf E}}_0}{[{{\bf I}} + {{\bf Y}}]^{ - 1}} \approx {{{\bf E}}_0}[{{\bf I}} - {{{\bf \hat{Y}}}_0}],$$ (A1b)where $${{\bf Y}} = {{\bf P}}_{}^ - {{{\bf A}}^ \cap }\quad {\rm and} \quad{{{\bf Y}}_0} = {{\bf P}}_0^ - {{{\bf A}}^ \cap }$$. Expressions (A1a) and (A1b) lead to the iterative solution at $${\rm z}_0^ +$$:   $$\Delta {{\bf Y}}_0^{(\rm i)} = {{\bf E}}_0^{({\rm i} - 1)}[{{\bf I}} - {{\bf Y}}_0^{({\rm i} - 1)}],$$ (A2a)  $${{\bf \hat{Y}}}_0^{(\rm i)} = {{\bf Y}}_0^{({\rm i} - 1)} + \Delta {{\bf Y}}_0^{(\rm i)},$$ (A2b)  $${{{\bf \hat{E}}}^{(\rm i)}}_0 = {{\bf Y}} - {\alpha ^{(\rm i)}}[{{\bf \hat{Y}}}_0^{(\rm i)} + {({{\bf \hat{Y}}}_0^{}{{\bf Y}})^{(\rm i)}}],$$ (A2c)  $${{\bf Y}}_0^{(\rm i)} = {\alpha ^{(\rm i)}}{{\bf \hat{Y}}}_0^{(\rm i)}\quad {\rm and} \quad{[{{\bf Y}}_0^{}{{\bf Y}}]^{(\rm i)}} = {\alpha ^{(\rm i)}}{[{{\bf \hat{Y}}}_0^{}{{\bf Y}}]^{(\rm i)}}.$$ (A2d) Note that $$[{{\bf Y}}_0^{(\rm i)} + {({{\bf Y}}_0^{}{{\bf Y}})^{(\rm i)}}]$$ is used as interpolated data at the missing data points. If we start in the first iteration (i = 1) with $${{\bf \hat{Y}}}_0^{(0)} = {{\bf Y}}$$, then it follows from expression (A1a) that $${{\bf E}}_0^{(0)} = - {{{\bf Y}}^2}$$ and the first estimates of $$\Delta {{\bf Y}}_0^{(1)}$$ and $${{\bf Y}}_0^{(1)}$$ can be computed with expressions (A2a)–(A2d). The elements in the diagonal matrix α(i) are determined such that the error vectors $$({\vec{\rm E}_{\rm j}} = {{\bf E}}{\vec{\rm I}_{\rm j}})$$ in error matrix E0 have minimum length at the data points:   $$E_0^{(\rm i)} = \sum\limits_{\rm j} {\sum\limits_\omega {{{\left\| {\left[ {{{\bf Y}} - {\alpha ^{(\rm i)}}[{{\bf \hat{Y}}}_0^{(\rm i)} + {{({{\bf \hat{Y}}}_0^{}{{\bf Y}})}^{(\rm i)}}]} \right]{{\vec{{\rm I}}}_{\rm j}}} \right\|}^{2}} = {\rm minimum}} } .$$ (A3) Note that in the last iteration α(i) ≈ I. If we also include errors in the input matrix, $${{\bf Y}} = {{\bf \hat{Y}}} + \Delta {{\bf Y}}$$, expressions (A1a) and (A1b) need be generalized to:   \begin{eqnarray} (1)\,{{{\bf E}}_{\rm tot}} &=& [{{\bf I}} - {{\bf \hat{Y}}}_0^{}][{{\bf I}} + {{\bf \hat{Y}}}] - {{\bf I}} = ({{\bf I}} - {{\bf Y}}_0^{} + \Delta {{\bf Y}}_0^{})({{\bf I}} + {{\bf Y}} - \Delta {{\bf Y}}) - {{\bf I}}\nonumber\\ &=& \Delta {{\bf Y}}_0^{}[{{\bf I}} + {{\bf Y}}] - [{{\bf I}} - {{\bf Y}}_0^{}]\Delta {{\bf Y}}\nonumber\\ & =& {{{\bf E}}_0} + {{\bf E}}. \end{eqnarray} (A4a)  $$(\rm 2a)\, \Delta {{\bf Y}}_0^{} = {{{\bf E}}_0}{[{{\bf I}} + {{\bf \hat{Y}}}]^{ - 1}} \approx {{{\bf E}}_0}[{{\bf I}} - {{{\bf \hat{Y}}}_0}],$$ (A4b)  $$(2{\rm b})\,\Delta {{\bf Y}} = - {[{{\bf I}} - {{{\bf \hat{Y}}}_0}]^{ - 1}}{{\bf E}} \approx - [{{\bf I}} + {{\bf \hat{Y}}}]{{\bf E}}.$$ (A4c) Now, $${{\bf Y}}_0^{}$$ and Y are alternately updated until E0 and E are sufficiently small. It is interesting that in this dual updating scheme both multiple removal and data interpolation are driven by the minimum-phase condition:   $$$$[{{\bf I}} - {{\bf Y}}_0^{}][{{\bf I}} + {{\bf Y}}] = {{\bf I}}.$$$$ (A5) APPENDIX B: ACOUSTIC VERSUS WEAK ELASTIC VERSUS FULL ELASTIC If T represents the transmission operator and we write T = I + δT then δT represents the forward scattering operator. If we explicitly include the transmission operator in our wavefield expressions, the backscattered wavefields $$\delta \vec{\rm P}_{\rm j}^ - = {{{\bf R}}^ \cup }\vec{\rm P}_{\rm j}^ +$$ and $$\delta \vec{\rm P}_{\rm j}^ + = {{{\bf R}}^ \cap }\vec{\rm P}_{\rm j}^ -$$ need be extended to secondary sources (Berkhout 2014):   $$\delta \vec{{\rm S}}_{\rm j}^ - \left( {{\rm z}_{\rm n}^ - ;{{\rm z}_0}} \right)\ = {{{\bf R}}^ \cup }\left( {{\rm z}_{\rm n}^ - ,{\rm z}_{\rm n}^ - } \right){\vec{\rm P}_{\rm j}}^ + \left( {{\rm z}_{\rm n}^ - ;{{\rm z}_0}} \right)\ + \delta {{{\bf T}}^ - }\left( {{\rm z}_{\rm n}^ - ,{\rm z}_{\rm n}^ + } \right){\vec{\rm P}_{\rm j}}^{ - }\left( {{\rm z}_{\rm n}^ + ;{{\rm z}_0}} \right)$$ (B1a)  $$\delta \vec{{\rm S}}_{\rm j}^ + \left( {{\rm z}_{\rm n}^ + ;{\rm z}_0^{}} \right) = {{{\bf R}}^ \cap }\left( {{\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ + } \right)\vec{\rm P}_{\rm j}^ - \left( {{\rm z}_{\rm n}^ + ;{\rm z}_0^{}} \right)\ + \delta {{{\bf T}}^ + }\left( {{\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ - } \right)\vec{\rm P}_{\rm j}^ + \left( {{\rm z}_{\rm n}^ - ;{\rm z}_0^{}} \right),$$ (B1b)representing the elastic interaction equations for P-waves. Note that R and δT play the same role in the scattering process. Eqs (B1a) and (B1b) can be seen as the extension of Huygens’ principle to inhomogeneous, elastic media. Note also that $$\delta \vec{{\rm S}}_{\rm j}^ - = \vec{\rm Q}_{\rm j}^ - - \vec{\rm P}_{\rm j}^ -$$ and $$\delta \vec{{\rm S}}_{\rm j}^ + = \vec{\rm Q}_{\rm j}^ + - \vec{\rm P}_{\rm j}^ +$$, see Fig. B1. Bear in mind that wavefields $${\rm Q}_{\rm j}^ \pm$$ are leaving a depth level and wavefields $${\rm P}_{\rm j}^ \pm$$ are entering a depth level. From expressions (B1a) and (B1b) it follows that for the elastic situation $$\delta \vec{{\rm S}}_{\rm j}^ -$$ and $$\delta \vec{{\rm S}}_{\rm j}^ +$$ are different. This difference can be written as:   $$\delta \vec{{\rm S}}_{\rm j}^ - \left( {{\rm z}_{\rm n}^ - ;{\rm z}_{\rm 0}^ - } \right) - \delta \vec{{\rm S}}_{\rm j}^ + \left( {{\rm z}_{\rm n}^ + ;{\rm z}_{\rm 0}^ + } \right) = \left[ {{{{\bf R}}^ \cup }\left( {{\rm z}_{\rm n}^ - ,{\rm z}_{\rm n}^ - } \right) - \delta {{{\bf T}}^ + }\left( {{\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ - } \right)} \right]{\vec{\rm P}_{\rm j}}^ + \left( {{\rm z}_{\rm n}^ - ;{\rm z}_0^{}} \right)\ - \left[ {{{{\bf R}}^ \cap }\left( {{\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ + } \right) - \delta {{{\bf T}}^ - }\left( {{\rm z}_{\rm n}^ - ,{\rm z}_{\rm n}^ + } \right)} \right]{\vec{\rm P}_{\rm j}}^ - \left( {{\rm z}_{\rm n}^ + ;{\rm z}_0^{}} \right).$$ (B2) Figure B1. View largeDownload slide The secondary sources $$(\delta \vec{{\rm S}}_{\rm j}^ \pm )$$ at a subsurface gridpoint represent the difference between the leaving $$(\vec{\rm Q}_{\rm j}^ \pm )$$ and the incident $$(\vec{\rm P}_{\rm j}^ \pm )$$ wavefields at that gridpoint. Hence, the secondary sources quantify the complex elastic wavefield interactions at each gridpoint. Figure B1. View largeDownload slide The secondary sources $$(\delta \vec{{\rm S}}_{\rm j}^ \pm )$$ at a subsurface gridpoint represent the difference between the leaving $$(\vec{\rm Q}_{\rm j}^ \pm )$$ and the incident $$(\vec{\rm P}_{\rm j}^ \pm )$$ wavefields at that gridpoint. Hence, the secondary sources quantify the complex elastic wavefield interactions at each gridpoint. It is easy to see from expression (B2) that $$\delta \vec{{\rm S}}_{\rm j}^ - = \delta \vec{{\rm S}}_{\rm j}^ +$$ if we assume δT + = R∪ and δT − = R∩ (acoustic situation). The acoustic assumption means a significant simplification of any imaging algorithm. The author thinks that this simplification is generally not justified in consolidated geological media and, therefore, should not be used in seismic imaging. Fig. B2 illustrates that if the shear contrast cannot be ignored, the acoustic simplification (δT = R) does not hold anymore. However, Fig. B2 also shows that a better approximation is obtained if we make the weak-elastic assumption: R∩ = −R∪ and δT − = −δT + or (R∪ − δT + ) = −(R∩ − δT − ). This assumption leads to:   $$\delta \vec{{\rm S}}_{\rm j}^ - \left( {{\rm z}_{\rm n}^ - ;{\rm z}_{\rm 0}^ - } \right) - \delta \vec{{\rm S}}_{\rm j}^ + \left( {{\rm z}_{\rm n}^ + ;{\rm z}_{\rm 0}^ + } \right) \approx \left[ {{{{\bf R}}^ \cup }\left( {{\rm z}_{\rm n}^ - ,{\rm z}_{\rm n}^ - } \right) - \delta {{{\bf T}}^ + }\left( {{\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ - } \right)} \right]\left[ {{{\vec{\rm P}}_{\rm j}}^ + \left( {{\rm z}_{\rm n}^ - ;{\rm z}_0^{}} \right)\ + {{\vec{\rm P}}_{\rm j}}^ - \left( {{\rm z}_{\rm n}^ + ;{\rm z}_0^{}} \right)} \right]$$ (B3a)or leaving the depth level out of the notation:   $$\delta \vec{{\rm S}}_{\rm j}^ - - \delta \vec{{\rm S}}_{\rm j}^ + \approx {{{\bf R}}^ - }({\vec{\rm P}_{\rm j}}^ + + {\vec{\rm P}_{\rm j}}^ - )\quad {\rm at} \quad {{\rm z}_{\rm n}},$$ (B3b)where R − = (R∪ − δT + ) = −(R∩ − δT − ) equals the conversion operator and where $$\delta \vec{{\rm S}}_{\rm j}^ - - \delta \vec{{\rm S}}_{\rm j}^ + = (\vec{\rm Q}_{\rm j}^ - + {\vec{\rm P}_{\rm j}}^ + ) - (\vec{\rm Q}_{\rm j}^ + + {\vec{\rm P}_{\rm j}}^ - )$$, being the difference between the total field above and below zn. Figure B2. View largeDownload slide Scattering is angle-dependent (upper part shows the forward and backward scattering operators; lower part shows the hybrid scattering operators). The angle-dependent wave conversion operator R − = (R∪ − δT + ) = −(R∩ − δT − ) is a measure for the wave conversion in a subsurface gridpoint. In subsurface areas with a significant R − , secondary S-sources are a logic extension to the secondary P-sources to generate the S-waves during the wavefield conversion process. Figure B2. View largeDownload slide Scattering is angle-dependent (upper part shows the forward and backward scattering operators; lower part shows the hybrid scattering operators). The angle-dependent wave conversion operator R − = (R∪ − δT + ) = −(R∩ − δT − ) is a measure for the wave conversion in a subsurface gridpoint. In subsurface areas with a significant R − , secondary S-sources are a logic extension to the secondary P-sources to generate the S-waves during the wavefield conversion process. The larger the difference between $$\delta \vec{{\rm S}}_{\rm j}^ -$$ and $$\delta \vec{{\rm S}}_{\rm j}^ +$$, the larger R − and, therefore, the larger the wave conversion. Note that the weak-elastic assumption also means:   $$\delta \vec{{\rm S}}_{\rm j}^ - + \delta \vec{{\rm S}}_{\rm j}^ + \approx {{{\bf R}}^ + }({\vec{\rm P}_{\rm j}}^ + - {\vec{\rm P}_{\rm j}}^ - )\quad {\rm at} \quad {{\rm z}_{\rm n}},$$ (B3c)where R + = (R∪ + δT + ) = −(R∩ + δT − ) represents the combined backward and forward scattering operator. Note that $$\delta \vec{{\rm S}}_{\rm j}^ - + \delta \vec{{\rm S}}_{\rm j}^ + = (\vec{\rm Q}_{\rm j}^ - + \vec{\rm Q}_{\rm j}^ + ) - ({\vec{\rm P}_{\rm j}}^ - + {\vec{\rm P}_{\rm j}}^ + ),$$ being the difference between the total leaving and the total incoming wavefields. Computation of R + and R − is an important option in our closed-loop, full-wavefield imaging algorithm. It leads to two new images that reveal the subsurface areas in which large wave conversion occurs. In these areas R and δT need be separately imaged. It is part of our current research on extending elastic P-wave imaging to S-wave imaging. APPENDIX C: MINIMIZATION ALGORITHM Let us consider a set of physical experiments, quantified by the following vector-matrix equations:   $$\begin{array}{*{20}{c}} {\vec{Y}_j^{} = {{\bf G}}\vec{X}_j^{}}&\quad{{\rm{for}}}&\quad {j{\rm{ }} = {\rm{ }}1,{\rm{ }}2, \ldots .,} \end{array}$$ (C1)where the elements of the system matrix G, indicated by Gkℓ, are the unknowns. These unknown elements can be estimated by a weighted least-squares minimization of the difference vector:   $$\Delta \vec{E}_j^\dagger {\Lambda _{jj}}\Delta {\vec{E}}\!\!\!\!\!\underline{\,\,\,}_j^{}\ = \ \vec{I}_j^\dagger \left[ {{{\bf Y}} - {{\bf \tilde{G}}}{{\bf \tilde{X}}}} \right]\ \Lambda {\left[ {{{\bf Y}} - {{\bf \tilde{G}}}{{\bf \tilde{X}}}} \right]^H}\vec{I}_j^{}\quad{\rm is\, minimum\,for\,all}\,j,$$ (C2)superscript † meaning a row vector, Λ being a user-specified weight matrix and H denoting Hermitian. In addition, $${{\bf \tilde{G}}} = {{\bf G}}{{{\bf L}}^{ - 1}}$$ and $${{\bf \tilde{X}}} = {{\bf LX}}$$, where operators L and L − 1 transform the internal coordinates of G and X to a suitable domain. In seismic imaging applications this transformation is often a τ − p transform, making the connection of G and X a connection per single plane wave component that occurs at each subsurface gridpoint. By setting the derivative of this weighted difference to $${\tilde{G}_{k\ell }}$$ and $${\tilde{{G} }\!\!\!\!\underline{\,\,\,}_{k\ell }}$$ equals zero for all k  and  ℓ, where $${{\bf \tilde{{G}}\!\!\!\!\underline{\,\,\,}}} = {{{\bf \tilde{G}}}^H}$$, then the solution is given by the normal equations (formulated in terms of matrices):   $${{\bf \tilde{G}}}\left[ {{{\bf \tilde{X}}}\Lambda {{{{\bf \tilde{X}}}}^H}} \right] = {{\bf Y}}\Lambda {{{\bf \tilde{X}}}^H},$$ (C3)where matrix Y contains all column vectors $${\vec{Y}_j}$$ and matrix $${{\bf \tilde{X}}}$$ contains all column vectors $${{\bf L}}{\vec{X}_j}$$. In its simplest version Λ is a diagonal matrix that contains the weights Λjj. To avoid inversion of correlation matrix $$[ {{{\bf \tilde{X}}}\Lambda {{{{\bf \tilde{X}}}}^H}} ]$$, we will use an iterative process. If we choose the elements of diagonal matrix Λ such that matrix $$[ {{{\bf \tilde{X}}}\Lambda {{{{\bf \tilde{X}}}}^H}} ]$$ has unit diagonal elements, then the first iteration is given by:   $${{{\bf \tilde{G}}}^{(1)}} = {{\bf Y}}\Lambda {{{\bf \tilde{X}}}^H},$$ (C4a)meaning that the first estimate of each element of $${{\bf \tilde{G}}}$$ is given by:   $${{{\bf \tilde{G}}}_{k\ell }}^{(1)} = \vec{I}_k^\dagger \left[ {{{\bf Y}}\Lambda {{{{\bf \tilde{X}}}}^H}} \right]\vec{I}_\ell ^{}.$$ (C4b) Of course, if a prior estimate of $${{\bf \tilde{G}}}$$ is already available, this prior estimate replaces expression (C4b). Using this first-order estimate, the error vector is given by   $$\Delta \vec{E}_j^{(1)} = \vec{Y}_j^{} - {{{\bf \tilde{G}}}^{(1)}}{{\bf L}}\vec{X}_j^{}.$$ (C4c) The first-order error vector is zero if the system would be orthogonal, meaning that $${{\bf \tilde{X}}}{{{\bf \tilde{X}}}^H}$$ equals a diagonal matrix. In practice such systems are rare and, therefore, let us look at the error vector of the second iteration:   \begin{eqnarray} \Delta \vec{E}_j^{(2)} &=& \vec{Y}_j^{} - {{{{\bf \tilde{G}}}}^{(2)}}{{\bf L}}\vec{X}_j^{}\nonumber\\ &=& \left[ {\vec{Y}_j^{} - {{{{\bf \tilde{G}}}}^{(1)}}{{\bf L}}\vec{X}_j^{}} \right] - \left[ {{{{{\bf \tilde{G}}}}^{(2)}} - {{{{\bf \tilde{G}}}}^{(1)}}} \right]{{\bf L}}\vec{X}_j^{}\nonumber\\ &=& \Delta \vec{E}_j^{(1)} - \left[ {{{{{\bf \tilde{G}}}}^{(2)}} - {{{{\bf \tilde{G}}}}^{(1)}}} \right]{{\bf L}}\vec{X}_j^{}. \end{eqnarray} (C5) If we substitute (C5) in eq. (C2), then we obtain the second-order minimization equation:   $$\Delta \vec{E}_j^{(2)}{\Lambda _{jj}}{\left[ {\Delta \vec{E}_j^{(2)}} \right]^H}\ = \ \left[ {\Delta \vec{E}_j^{(1)} - \Delta {{{{\bf \tilde{G}}}}^{(2)}}{{\bf L}}\vec{X}_j^{}} \right]\ {\Lambda _{jj}}{\left[ {\Delta \vec{E}_j^{(1)} - \Delta {{{{\bf \tilde{G}}}}^{(2)}}{{\bf L}}\vec{X}_j^{}} \right]^H}{\rm is\,minimum\,for\,all}\,j,$$ (C6)where $$\Delta {{{\bf \tilde{G}}}^{(2)}} = {{{\bf \tilde{G}}}^{(2)}} - {{{\bf \tilde{G}}}^{(1)}}$$. Minimization of (C6) leads to the second-order solution:   $$\Delta {{{\bf \tilde{G}}}^{(2)}} = \Delta {{{\bf E}}^{(1)}}\Lambda {{{\bf \tilde{X}}}^H},$$ (C7a)or   $${{{\bf \tilde{G}}}^{(2)}} = {{{\bf \tilde{G}}}^{(1)}} + \Delta {{{\bf E}}^{(1)}}\Lambda {{{\bf \tilde{X}}}^H}.$$ (C7b) Second-order solution (C7b) can be generalized to:   $${{{\bf \tilde{G}}}^{(i)}} = {{{\bf \tilde{G}}}^{({i} - 1)}} + \Delta {{{\bf E}}^{({i} - 1)}}\Lambda {{{\bf \tilde{X}}}^H}$$ (C8)with  \begin{equation*}{{{\bf \tilde{G}}}^{(0)}} = {{\bf 0}},\quad\Delta {{{\bf E}}^{(0)}} = {{\bf Y}}\quad{\rm{and}}\quad\Delta {{{\bf E}}^{(1)}} = {{\bf Y}} - {{{\bf \tilde{G}}}^{(1)}}{{\bf \tilde{X}}}.\end{equation*} Solution (C8) is used in all iterative estimation algorithms for the wavefield operators. Note that if we deal with the related equations:   $$\begin{array}{*{20}{c}} {\vec{Y}_k^{} = {{\bf X}}\vec{G}_k^{}}&{{\rm{for}}}&\quad{k = 1,{\rm{ }}2, \ldots .\,,} \end{array}$$ (C9a)the iterative solution is given by:   $${{{\bf \tilde{G}}}^{(i)}} = {{{\bf \tilde{G}}}^{({i} - 1)}} + {{{\bf \tilde{X}}}^H}\Lambda \Delta {{{\bf E}}^{({i} - 1)}},$$ (C9b)where $$\Delta {{\bf E}} = {{\bf Y}} - {{\bf \tilde{X}\tilde{G}}}.$$ Note also that in the special situation of matrix inversion,   $${{\bf I}} = {{\bf FX}}\quad{\rm or}\quad{{\bf I}} = {{\bf XF}},$$ (C10a)the iterative solution is given by:   $${{{\bf \tilde{F}}}^{(i)}} = {{{\bf \tilde{F}}}^{({i} - 1)}} + \Delta {{{\bf E}}^{({i} - 1)}}\Lambda {{{\bf \tilde{X}}}^H}\quad{\rm or}\quad{{{\bf \tilde{F}}}^{(i)}} = {{{\bf \tilde{F}}}^{({i} - 1)}} + {{{\bf \tilde{X}}}^H}\Lambda \Delta {{{\bf E}}^{({i} - 1)}},$$ (C10b)where   \begin{equation*}\Delta {{\bf E}} = {{\bf I}} - {{\bf \tilde{F}\tilde{X}}}\quad{\rm or}\quad\Delta {{\bf E}} = {{\bf I}} - {{\bf \tilde{X}\tilde{F}}}\end{equation*} APPENDIX D: SUMMARY OF CONCEPTS AND NOTATION If the time domain is transformed to the frequency domain, a seismic data volume consists of orthogonal slices. Each frequency slice is represented by the data matrix. A measured physical wavefield is represented by one column of the data matrix, ordered according to the x–y observation coordinates (detector positions). A row represents a common detector gather, now ordered by the x-y emission coordinates (source positions). Each column and row may have its own source and detector depth level. Primary sources are man-made and often positioned (close) at the surface. Wavefields consist of downgoing and upgoing constituents, indicated by the uppercases + and − respectively. The definition down and up depend on the orientation of the coordinate system used. Downgoing wavefields illuminate the upper side of reflectors, upgoing wavefields illuminate the lower-side of reflectors. This article focuses on P-waves, but wave conversion during the PP-scattering process is taken into account. We used the following notation: $${{{\bf P}}^ + }({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}})$$: downgoing wavefields that enter the upper part of depth level $${\rm z}_{\rm m}^{}$$ (realizing illumination from above). $${{{\bf P}}^ - }({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})$$: upgoing wavefields that enter the lower part of depth level $${\rm z}_{\rm m}^{}$$ (realizing illumination from below). $${{{\bf Q}}^ + }({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})$$: downgoing wavefields that leave the lower part of depth level $${\rm z}_{\rm m}^{}$$; note that $${{{\bf P}}^ + }({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}) = {{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_{{\rm m} - 1}^ + ){{{\bf Q}}^ + }({\rm z}_{{\rm m} - 1}^ + ;{{\rm z}_{\rm s}})$$, where each column of W equals a local downward wavefield propagation operator (generally being anisotropic). $${{{\bf Q}}^ - }({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}})$$: upgoing wavefields that leave the upper part of depth level $${\rm z}_{\rm m}^{}$$; note that $${{{\bf P}}^ - }({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}}) = {{\bf W}}({\rm z}_{\rm m}^ + ,{\rm z}_{\rm m + 1}^ - ){{{\bf Q}}^ - }({\rm z}_{\rm m + 1}^ - ;{{\rm z}_{\rm s}})$$, where each column of W equals the local upward wavefield propagation operator (generally being anisotropic). Decomposition of the downgoing wavefields for triple illumination from above:   \begin{equation*} {{{\bf P}}^ + }({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}) = {{\bf W}}({\rm z}_{\rm m}^ - ,{{\rm z}_{\rm s}}){{{\bf S}}^ + }({{\rm z}_{\rm s}}) + {{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_0^ + )[{{{\bf R}}^ \cap }({\rm z}_0^ + ,{\rm z}_0^ + ){{{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}})] + {{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_{\rm n}^ + )\sum\limits_{\rm n = 1}^{{\rm m} - 1} {{{{\bf R}}^ \cap }({\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ + ){{{\bf P}}^ - }({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}})}, \end{equation*} where the first term equals the downgoing wavefield generated by the primary sources, the second term equals the downgoing wavefield generated by the secondary sources at the surface (surface backscattering) and the third term equals the down going wavefield generated by the secondary sources in the subsurface (internal backscattering). Decomposition of the upgoing wavefields for triple illumination from below:   \begin{equation*} \quad{{{\bf P}}^ - }({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}) = {{\bf P}}_{\rm pr}^ - ({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}) + {{\bf M}}_{\rm sfc}^ - ({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}) + {{\bf M}}_{{\mathop{\rm int}} }^ - ({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}), \end{equation*} where the first term equals the first-order response from the man-made sources (‘primaries’), $${{\bf P}}_{\rm pr}^ - = {{\bf X}}_{\rm pr}^ \cup {{{\bf W}}^ + }{{{\bf S}}^ + }$$ with $${{\bf X}}_{\rm pr}^ \cup = \sum\nolimits_{\rm n > {\rm m}} {{{{\bf W}}^ - }{{\bf R}}_{\rm a}^ \cup {{{\bf W}}^ + }}$$, the second term equals the response from the surface backscattering (‘surface multiples’) $${{\bf M}}_{\rm sfc}^ - = {{\bf X}}_{\rm pr}^ \cup {{{\bf W}}^ + }\delta {{\underline{\bf P}}^ + }$$ with $$\delta {{\underline{\bf P}}^ + } = {{\bf R}}_{}^ \cap {{\underline{\bf P}}^ - }$$ at $${\rm z}_0^ +$$ and $${{\bf X}}_{\rm pr}^ \cup = \sum\nolimits_{\rm n > {\rm m}} {{{{\bf W}}^ - }{{\bf R}}_{\rm b}^ \cup {{{\bf W}}^ + }}$$, and the third term equals the response from the internal backscattering (‘internal multiples’), $${{\bf M}}_{{\mathop{\rm int}} }^ - = {{\bf X}}_{\rm pr}^ \cup \sum\nolimits_{\ell < {\rm m}} {{{{\bf W}}^ + }\delta {{{\bf P}}^ + }}$$ with $$\delta {{{\bf P}}^ + } = {{\bf R}}_{}^ \cap {{{\bf P}}^ - }$$ at $${\rm z}_\ell ^ +$$ with ℓ > 0 and $${{\bf X}}_{\rm pr}^ \cup = \sum\nolimits_{\rm n > {\rm m}} {{{{\bf W}}^ - }{{\bf R}}_{\rm c}^ \cup {{{\bf W}}^ + }}$$. One column of back scattering operators R∪ (up) and R∩ (down) defines the angle-dependent scattering properties in one gridpoint. These properties are determined by the faster changes of the density (ρ) and the velocities (cp, cs) in that gridpoint. Operators R∪ and R∩ are output of the closed-loop, full-wavefield imaging algorithm. By definition they include wave conversion. One column of propagation operators W − (up) and W + (down) defines the anisotropic propagation properties between gridpoints. These properties are determined by the slower changes of the velocity (cp) between these gridpoints. Operators W − and W + are output of the closed-loop, full-wavefield imaging algorithm. In this article they propagate P-waves only. Wavefield extrapolation may be carried out on the down-and upgoing constituents (P ± , Q ± ) or on the combined versions (Q + + P − , Q + − P − ) and (Q − + P + , Q − − P + ). Estimation of operators R∪ and R∩ always occurs on the separate constituents (P ± , Q ± ), preferably in the linear Radon (τ-p) domain. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Closed-loop multiple-scattering imaging with sparse seismic measurements

, Volume 212 (3) – Mar 1, 2018
25 pages

/lp/ou_press/closed-loop-multiple-scattering-imaging-with-sparse-seismic-IWJ71nxa8l
Publisher
Oxford University Press
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx476
Publisher site
See Article on Publisher Site

### Abstract

Summary In the theoretical situation of noise-free, complete data volumes (‘perfect data’), seismic data matrices are fully filled and multiple-scattering operators have the minimum-phase property. Perfect data allow direct inversion methods to be successful in removing surface and internal multiple scattering. Moreover, under these perfect data conditions direct source wavefields realize complete illumination (no irrecoverable shadow zones) and, therefore, primary reflections (first-order response) can provide us with the complete seismic image. However, in practice seismic measurements always contain noise and we never have complete data volumes at our disposal. We actually deal with sparse data matrices that cannot be directly inverted. The message of this paper is that in practice multiple scattering (including source ghosting) must not be removed but must be utilized. It is explained that in the real world we badly need multiple scattering to fill the illumination gaps in the subsurface. It is also explained that the proposed multiple-scattering imaging algorithm gives us the opportunity to decompose both the image and the wavefields into order-based constituents, making the multiple scattering extension easy to apply. Last but not least, the algorithm allows us to use the minimum-phase property to validate and improve images in an objective way. Seismic interferometry, Wave propagation, Wave scattering and diffraction INTRODUCTION Multiple scattering events (‘multiples’) have always played a vital role in seismic imaging. Traditionally, they are considered as unwanted events (being referred to as ‘shot-generated noise’) that interfere with—or may even overshadow—the primary reflections (the latter being considered as the ‘shot-generated signal’). A representative example of scientific arguments to consider multiple scattering as shot-generated noise can be found in Weglein (2014, 2015). In the past 50 yr enormous investments have been made to remove multiples, making the seismic reflection response linear in reflectivity. However, recently the seismic community started to realize that lack of understanding was the reason of this linearization strategy. There is increasing theoretical and practical evidence that multiples contain subsurface information that cannot be extracted from primary reflections, see for instance Whitmore et al. (2010) and Davidenko et al. (2015). Hence, multiples should not be removed as noise but they should be utilized as signal. The only noise left in seismic data is background noise (Berkhout 2014). The author considers utilization of multiples the next big step forward in the practice of seismic imaging (often referred to as ‘seismic migration’). We will see that this also applies to the utilization of source ghosts (the interaction of near-surface sources with the surface). Feedback model at the surface $${(z_{0}^{+})}$$ In marine data—and in hard-surface land data—multiples that are generated by the surface reflector are strong and must be addressed. Fig. 1(a) shows the feedback model of the surface-related multiple generation process (Berkhout 1982). Figure 1. View largeDownload slide (a) Feedback model at the surface, showing the generation of the surface-related multiples $$({{\bf M}}_0^ - = {{\bf P}}_0^ - {{{\bf A}}^ \cap }{{\bf P}}_{}^ - )$$. Note the important role of multiple generation operator [I + A∩P − ] in the downgoing (Q + ) and upgoing (P − ) wavefields at $${\rm z}_0^ +$$, where A∩ = (S + ) − 1R∩. (b) Feedback model at the surface, showing the generation of both the source ghost response and the surface-related multiples in the data (‘surface multiple scattering’). Note the extended multiple generation operator $$[{{\bf I}} + {{{\bf A}}^ \cap }{{{\underline{\bf P}}}^ - }]$$ in the downgoing (Q + ) and upgoing (P − ) wavefields. Extended multiple removal has a double function, i.e., simultaneously removal of the source ghost response and the surface-related multiples $$({{\bf M}}_0^ - )$$. Figure 1. View largeDownload slide (a) Feedback model at the surface, showing the generation of the surface-related multiples $$({{\bf M}}_0^ - = {{\bf P}}_0^ - {{{\bf A}}^ \cap }{{\bf P}}_{}^ - )$$. Note the important role of multiple generation operator [I + A∩P − ] in the downgoing (Q + ) and upgoing (P − ) wavefields at $${\rm z}_0^ +$$, where A∩ = (S + ) − 1R∩. (b) Feedback model at the surface, showing the generation of both the source ghost response and the surface-related multiples in the data (‘surface multiple scattering’). Note the extended multiple generation operator $$[{{\bf I}} + {{{\bf A}}^ \cap }{{{\underline{\bf P}}}^ - }]$$ in the downgoing (Q + ) and upgoing (P − ) wavefields. Extended multiple removal has a double function, i.e., simultaneously removal of the source ghost response and the surface-related multiples $$({{\bf M}}_0^ - )$$. Using the matrix operator notation (see also Appendix  D), the expression of the upgoing wavefields $$({{\bf P}}_{}^ - )$$ arriving at the surface (z0) can be written as:   $${{\bf P}}_{}^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = {{\bf X}}_{}^ \cup ({\rm z}_0^ + ,{\rm z}_0^ + ){{{\bf S}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}})$$ (1a)or   $${{{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = {{\bf X}}_0^ \cup ({\rm z}_0^ + ,{\rm z}_0^ + ){{{\bf Q}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}})$$ (1b)with $${\rm z}_0^ +$$ being the lower part of the surface and   $${{{\bf Q}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = {{{\bf S}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) + {{{\bf R}}^ \cap }({\rm z}_0^ + ,{\rm z}_0^ + ){{{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}}),$$ (1c)where $${{{\bf S}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = {{{\bf W}}^*}({\rm z}_0^ + ,{\rm z}_{\rm s}^{}){{{\bf S}}^ + }({{\rm z}_{\rm s}})$$ with W being the scattering-free propagation matrix and the complex conjugate version W* being its inverse. Note that source depth level zs may be generalized to zs(x, y), meaning that each column of matrix S + (zs) may have its own depth level. In the above expressions each column of wavefield matrix Q + represents the total downgoing wavefield that leaves the surface, each column of the source matrix S + represents the downgoing source (array) wavefield that excludes the source ghost (will be included later), one column of the reflectivity matrix $${{\bf R}}_{}^ \cap$$ represents the angle-dependent reflection operator for upgoing wavefields at each gridpoint of the surface and one column of the wavefield matrix $${{\bf P}}_{}^ -$$ represents the upgoing wavefields in one shot record. Finally, the surface-included earth's impulse response matrix $${{\bf X}}_{}^ \cup$$ transforms the downgoing source wavefields $${{\bf S}}_{}^ +$$ into the upgoing reflected wavefields $${{\bf P}}_{}^ -$$ and the surface-excluded earth's impulse response matrix $${{\bf X}}_0^ \cup$$ transforms the downgoing total wavefields $${{\bf Q}}_{}^ +$$ into the upgoing reflected wavefields $${{\bf P}}_{}^ -$$. In the ideal situation S + (zs) equals the unity matrix (S + = I) with a unit source element at each gridpoint of depth level zs. Note also that from the physics point of view, $${{\bf X}}_0^ \cup$$ represents one roundtrip (from surface into the subsurface and back to the surface) and $${{\bf X}}_{}^ \cup$$ represents many roundtrips. This can be easily seen from the feedback model in Fig. 1(a). For the response of a single, possibly blended source wavefield—blending meaning the creation of a decomposable areal source wavefield—matrix S + is reduced to column vector $$\vec{{\rm S}}_{\rm j}^ + = {{\bf S}}_{}^ + {\vec{\rm I}_{\rm j}}$$ and expression (1c) becomes a vector-matrix equation:   $$\vec{\rm Q}_{\rm j}^ + ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = \vec{{\rm S}}_{\rm j}^ + ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) + {{{\bf R}}^ \cap }({\rm z}_0^ + ,{\rm z}_0^ + )\vec{\rm P}_{\rm j}^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}),$$ (1d)showing that each gridpoint at the surface functions as one secondary source element of the naturally blended source array $${{{\bf R}}^ \cap }\vec{\rm P}_{\rm j}^ -$$. It is the response of this naturally blended source array that results in the surface-related multiples: $$({{\bf M}}_0^ - {\vec{\rm I}_{\rm j}} = {{\bf X}}_0^ \cup [{{{\bf R}}^ \cap }{{\bf P}}_{}^ - ]{\vec{\rm I}_{\rm j}} = {{\bf X}}_0^ \cup \delta {{\bf P}}_{}^ + {\vec{\rm I}_{\rm j}} = {{\bf X}}_0^ \cup \delta \vec{\rm P}_{\rm j}^ + )$$. If $${{\bf D}}_{}^ - ({{\rm z}_{\rm d}};{\rm z}_0^ + )$$ represents the detector matrix that selects the upgoing wavefields at detector level zd, where D − includes the detector (array) directivity as well, then the wavefield matrix $${{\bf P}}_{}^ -$$ is given by:   \begin{eqnarray} {{{\bf P}}^ - }({{\rm z}_{\rm d}};{{\rm z}_{\rm s}}) = {{{\bf D}}^ - }({{\rm z}_{\rm d}};{\rm z}_0^ + ){{\bf P}}_{}^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = {{{\bf D}}^ - }({{\rm z}_{\rm d}};{\rm z}_0^ + ){{\bf X}}_{}^ \cup ({\rm z}_0^ + ,{\rm z}_0^ + ){{{\bf S}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}})\end{eqnarray} (2a)or   $${{{\bf P}}^ - }({{\rm z}_{\rm d}};{{\rm z}_{\rm s}}) = {{{\bf D}}^ - }({{\rm z}_{\rm d}};{\rm z}_0^ + ){{\bf X}}_0^ \cup ({\rm z}_0^ + ,{\rm z}_0^ + ){{{\bf Q}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}})$$ (2b)with   $${{\bf D}}_{}^ - ({{\rm z}_{\rm d}};{\rm z}_0^ + ) = {{\bf D}}_{}^ - ({{\rm z}_{\rm d}}){{\bf W}}_{}^*({{\rm z}_{\rm d}},{\rm z}_0^ + ).$$ (2c)Similar to source level zs, detector level zd may be generalized to zd(x, y), meaning that each column of matrix D − (zd) may have its own depth level. Expressions (2a)–(2c) quantify the influence of acquisition matrices, $${{\bf S}}_{}^ + ({{\rm z}_{\rm s}})$$ and D − (zd) at the start and the end of the wavefield life cycle. In the 3-D practice both matrices are far from complete and, therefore, their influence is significant on the data space. This influence occurs in the common source gathers (columns of $${{\bf P}}_{}^ -$$) and common detector gathers (rows of $${{\bf P}}_{}^ -$$) as spatial aliasing and spatial band limitation. In industrial acquisition design, the elements of $${{\bf S}}_{}^ +$$ and D − are chosen such that aliasing and band limitation are minimized under pre-specified economical and logistical constraints. These constraints unavoidably lead to sparse 3-D acquisition matrices that cannot be directly inverted. We will see that this practical reality has an important influence on the design of effective multiple removal algorithms. The same is true for full wavefield redatuming algorithms and full wavefield imaging algorithms. Surface scattering removal for perfect data volumes Using expressions (1b) and (1c) and (2b) and (2c), we may write (see Fig. 1a):   \begin{eqnarray} {{{\bf P}}^ - }({{\rm z}_{\rm d}};{{\rm z}_{\rm s}}) &=& {{{\bf D}}^ - }({{\rm z}_{\rm d}};{\rm z}_0^ + ){{{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}})\nonumber\\ &=& {{{\bf D}}^ - }({{\rm z}_{\rm d}};{\rm z}_0^ + ){{\bf X}}_0^ \cup ({\rm z}_0^ + ,{\rm z}_0^ + ){{{\bf Q}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) \nonumber\\ &=& {{{\bf D}}^ - }({{\rm z}_{\rm d}};{\rm z}_0^ + ){{\bf X}}_0^ \cup ({\rm z}_0^ + ,{\rm z}_0^ + )[{{{\bf S}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) + {{{\bf R}}^ \cap }({\rm z}_0^ + ,{\rm z}_0^ + ){{{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}})] \nonumber\\ &=& {{\bf P}}_0^ - ({{\rm z}_{\rm d}};{\rm z}_{\rm s}^{})[{{\bf I}} + {{{\bf A}}^ \cap }({\rm z}_{\rm s}^{},{\rm z}_{\rm d}^{}){{{\bf P}}^ - }({\rm z}_{\rm d}^{};{{\rm z}_{\rm s}})], \end{eqnarray} (3a)where   $${{\bf P}}_0^ - ({{\rm z}_{\rm d}};{{\rm z}_{\rm s}}) = {{{\bf D}}^ - }({{\rm z}_{\rm d}};{\rm z}_0^ + ){{\bf X}}_0^ \cup ({\rm z}_0^ + ,{\rm z}_0^ + ){{{\bf S}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}}),$$ (3b)  $${{\bf A}}_{}^ \cap ({{\rm z}_{\rm s}},{{\rm z}_{\rm d}}) = {[{{{\bf S}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}})]^{ - 1}}{{{\bf R}}^ \cap }({\rm z}_0^ + ,{\rm z}_0^ + ){[{{\bf D}}_{}^ - ({{\rm z}_{\rm d}};{\rm z}_0^ + )]^{ - 1}}.$$ (3c)Note also that:   $${{{\bf Q}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = {{\bf S}}_{}^ + ({\rm z}_0^ + ;{{\rm z}_{\rm s}})[{{\bf I}} + {{{\bf A}}^ \cap }({{\rm z}_{\rm s}},{{\rm z}_{\rm d}}){{{\bf P}}^ - }({\rm z}_{\rm d}^{};{{\rm z}_{\rm s}})].$$ (3d)If we leave the source and detector level out of the notation, we may simplify expressions (3d) and (3a):   $${{{\bf Q}}^ + } = {{\bf S}}_{}^ + [{{\bf I}} + {{{\bf A}}^ \cap }{{{\bf P}}^ - }] \quad{\rm at}\quad {\rm z}_0^ +$$ (4a)  $${{{\bf P}}^ - } = {{\bf P}}_0^ - [{{\bf I}} + {{{\bf A}}^ \cap }{{{\bf P}}^ - }]\quad{\rm at}\quad {\rm z}_{\rm d}^{},$$ (4b)showing that [I + A∩P − ] is the surface multiple generation operator, transforming $${{\bf S}}_{}^ +$$ into $${{\bf Q}}_{}^ +$$ and $${{\bf P}}_0^ -$$ into P − . From expressions (4a) and (4b) it directly follows that   $${{\bf S}}_{}^ + = {{\bf Q}}_{}^ + {[{{\bf I}} + {{{\bf A}}^ \cap }{{{\bf P}}^ - }]^{ - 1}} \quad{\rm at}\quad {\rm z}_0^ + ,$$ (5a)  $${{\bf P}}_0^ - = {{{\bf P}}^ - }{[{{\bf I}} + {{{\bf A}}^ \cap }{{{\bf P}}^ - }]^{ - 1}} \quad{\rm at}\quad {\rm z}_{\rm d}^{},$$ (5b)showing that [I + A∩P − ] − 1 equals the surface multiple removal operator, transforming $${{\bf Q}}_{}^ +$$ into $${{\bf S}}_{}^ +$$ and P − into $${{\bf P}}_0^ -$$. Making use of the series expansion of the removal operator, direct inversion expression (5b) may also be written as:   \begin{eqnarray} {{\bf P}}_0^ - &=& {{{\bf P}}^ - } - {{{\bf P}}^ - }[{{{\bf A}}^ \cap }{{{\bf P}}^ - }] + {{{\bf P}}^ - }{[{{{\bf A}}^ \cap }{{{\bf P}}^ - }]^2} - {{{\bf P}}^ - }{[{{{\bf A}}^ \cap }{{{\bf P}}^ - }]^3} + \cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\nonumber\\ &=& {{{\bf P}}^ - } - [{{{\bf P}}^ - }{{{\bf A}}^ \cap }]{{{\bf P}}^ - } + {[{{{\bf P}}^ - }{{{\bf A}}^ \cap }]^2}{{{\bf P}}^ - } - {[{{{\bf P}}^ - }{{{\bf A}}^ \cap }]^3}{{{\bf P}}^ - } + \cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\quad {\rm at}\quad {\rm z}_{\rm d}^{}. \end{eqnarray} (5c)Expression (5c) represents the well-known inverse scattering series at the surface (Berkhout 1982; Araujo et al.1994; Verschuur & Berkhout 1997). If all elements of surface operator A∩ exist and are known (‘perfect data’), we have two ways of surface-multiple removal by applying direct inversion: Eq. (5b) is used, meaning that multiple-generating matrix operator [I + A∩P − ] need be inverted for; Inverse scattering series (5c) is used, meaning that each previous term need be multiplied with [A∩P − ] at the source side or [P − A∩] at the detector side. Now, let us return to expressions (4a) and (4b). If we pose that the determination of impulse response (IR) matrix $${{\bf X}}_0^ \cup$$ is the ultimate aim of surface-related pre-processing (leading to a subsurface response that has been made independent of the conditions at the surface), then we may write:   $${{\bf X}}_0^ \cup = {({{{\bf D}}^ - })^{ - 1}}{{{\bf P}}^ - }{({{{\bf Q}}^ + })^{ - 1}} = {({{{\bf D}}^ - })^{ - 1}}\frac{{{{\bf P}}_0^ - [{{\bf I}} + {{{\bf A}}^ \cap }{{{\bf P}}^ - }]}}{{{{\bf S}}_{}^ + [{{\bf I}} + {{{\bf A}}^ \cap }{{{\bf P}}^ - }]}} = {({{{\bf D}}^ - })^{ - 1}}{{\bf P}}_0^ - {({{\bf S}}_{}^ + )^{ - 1}},$$ (6)assuming again that the elements of A∩ exist at every gridpoint of the source and detector depth levels. Expression (6) reveals that in the situation of a perfect data volume (no noise, no missing data), surface multiples do not provide extra information for the determination of subsurface IR matrix $${{\bf X}}_0^ \cup$$! Expression (6) also explains that for the determination of $${{\bf X}}_0^ \cup$$ the complex division with double illumination matrix $$({{\bf Q}}_{}^ + )$$ can be avoided by (i) surface multiple removal using the inverse scattering series, followed by (ii) source-detector deconvolution:   $$({\rm i})\,\,{{\bf P}}_0^ - = {{{\bf P}}^ - }{({{\bf I}} + {{{\bf A}}^ \cap }{{{\bf P}}^ - })^{ - 1}} = {{{\bf P}}^ - } - {{{\bf P}}^ - }({{{\bf A}}^ \cap }{{{\bf P}}^ - }) + {{{\bf P}}^ - }{({{{\bf A}}^ \cap }{{{\bf P}}^ - })^2} - {{{\bf P}}^ - }{({{{\bf A}}^ \cap }{{{\bf P}}^ - })^3} + \cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot$$ (7a)or   $${{\bf P}}_0^ - = {({{\bf I}} + {{{\bf P}}^ - }{{{\bf A}}^ \cap })^{ - 1}}{{{\bf P}}^ - } = {{{\bf P}}^ - } - ({{{\bf P}}^ - }{{{\bf A}}^ \cap }){{{\bf P}}^ - } + {({{{\bf P}}^ - }{{{\bf A}}^ \cap })^2}{{{\bf P}}^ - } - {({{{\bf P}}^ - }{{{\bf A}}^ \cap })^3}{{{\bf P}}^ - } + \cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot$$ (7b)and   $$({\rm ii})\,\,{{\bf X}}_0^ \cup = {({{{\bf D}}^ - })^{ - 1}}{{\bf P}}_0^ - {({{{\bf S}}^ + })^{ - 1}}$$ (7c)at each gridpoint on detector level zd and source level zs respectively. Surface scattering removal for sparse data volumes Due to economical and logical constraints in practice, D − and S + are sparse matrices, meaning that [D − ] − 1 and [S + ] − 1 do not exist and therefore operator A∩ does not exist (from theoretical to real data). Moreover, $${{\bf P}}_{}^ -$$ contains always noise and full knowledge of true source and detector behaviour is not available (D − and S + are only approximately known). This means that direct inversion as formulated by expression (6) causes significant problems in practice. In the following, an alternative way of thinking is introduced with respect handling multiples in sparse data volumes. First, let us reveal a fundamental property of the multiple-generating system. From expression (4b) it follows that:   $$\begin{array}{@{}l@{}} {{\bf P}}_0^ - = {{{\bf P}}^ - } - {{\bf P}}_0^ - {{{\bf A}}^ \cap }{{{\bf P}}^ - } = [{{\bf I}} - {{\bf P}}_0^ - {{{\bf A}}^ \cap }]{{{\bf P}}^ - }\quad {\rm at} \quad {\rm z}_{\rm d}^{}, \end{array}$$ (8a)showing that $${{\bf P}}_0^ - {{{\bf A}}^ \cap }$$ functions as a prediction operator. Using expression (7b), we may also write:   $${{\bf P}}_0^ - = {[{{\bf I}} + {{{\bf P}}^ - }{{{\bf A}}^ \cap }]^{ - 1}}{{{\bf P}}^ - }\quad {\rm at} \quad {\rm z}_{\rm d}^{}.$$ (8b)If we compare expression (8a) with expression (8b) it follows that:   $${[{{\bf I}} + {{{\bf P}}^ - }{{{\bf A}}^ \cap }]^{ - 1}} = [{{\bf I}} - {{\bf P}}_0^ - {{{\bf A}}^ \cap }]\quad {\rm at} \quad {\rm z}_{\rm d}^{}$$ (9a)or   $${[{{\bf I}} - {{\bf P}}_0^ - {{{\bf A}}^ \cap }]^{ - 1}} = [{{\bf I}} + {{{\bf P}}^ - }{{{\bf A}}^ \cap }]\quad {\rm at} \quad {\rm z}_{\rm d}^{}.$$ (9b)Expressions (9a) and (9b) tell us that if surface operator A∩ exists, meaning that at every gridpoint wavefields are available, we may conclude that: The inverse of causal operator [I + P − A∩] exists The inversion result $$[{{\bf I}} - {{\bf P}}_0^ - {{{\bf A}}^ \cap }]$$ is again causal. These two conditions are exactly the two necessary and sufficient conditions for an operator to have the minimum-phase property (Berkhout 1971). Hence, expression (9b) extends the classical minimum-phase property of a time-series to the 3-D seismic surface multiple-generating system. Note that the minimum-phase condition for a wavefield is significantly more than the causality condition. In all areas of science (think particularly at today's climate sciences) the practical relevance of minimum phase is that the future can be fully predicted by knowledge of the past! Later in this paper, it will be shown that this fundamental property also applies to the Earth's internal multiple-generating system. Minimum-phase test for multiple removal solutions From expression (9b) it follows that solution $${{\bf P}}_0^ -$$ must obey the condition:   $$[{{\bf I}} - {{\bf P}}_0^ - {{{\bf A}}^ \cap }][{{\bf I}} + {{{\bf P}}^ - }{{{\bf A}}^ \cap }] = {{\bf I}} \qquad \hbox{(minimum-phase condition)}$$ (10a)or, simplifying the notation by writing $${{{\bf Y}}_0} = {{\bf P}}_0^ - {{{\bf A}}^ \cap }$$ and Y = P − A∩,   $$$$[{{\bf I}} - {{\bf Y}}_0^{}][{{\bf I}} + {{\bf Y}}] = {{\bf I}}.$$$$ (10b)If $${{{\bf \hat{Y}}}_0}$$ represents some initial estimate of Y0, then error matrix E0 is given by:   $${{{\bf E}}_0} = [{{\bf I}} - {{\bf \hat{Y}}}_0^{}][{{\bf I}} + {{\bf Y}}] - {{\bf I}}.$$ (11a)Note that E0 is a measure for the violation of the minimum-phase property due to a wrong prediction result. It is also a measure for the difference $$[{{{\bf Y}}_0} - {{{\bf \hat{Y}}}_0}]$$. If we write $${{{\bf Y}}_0} = {{{\bf \hat{Y}}}_0} + \Delta {{{\bf Y}}_0}$$, then it follows from expression (11a) that   $$[{{\bf I}} - {{\bf \hat{Y}}}_0^{}][{{\bf I}} + {{\bf Y}}] = {{\bf I}} + \Delta {{\bf Y}}_0^{}[{{\bf I}} + {{\bf Y}}] = {{\bf I}} + {{{\bf E}}_0}.$$ (11b)Using expressions (11a) and (11b) we may conclude that at $${\rm z}_0^ +$$:   $${\rm (1)}\,\,{{{\bf E}}_0} = [{{\bf I}} - {{\bf \hat{Y}}}_0^{}][{{\bf I}} + {{\bf Y}}] - {{\bf I}} = {{\bf Y}} - {{\bf \hat{Y}}}_0^{} - {{\bf \hat{Y}}}_0^{}{{\bf Y}}$$ (12a)  $${\rm (2)}\,\,\Delta {{\bf Y}}_0^{} = {{{\bf E}}_0}{[{{\bf I}} + {{\bf Y}}]^{ - 1}} \approx {{{\bf E}}_0}[{{\bf I}} - {{{\bf \hat{Y}}}_0}].$$ (12b)In Appendix  A the iterative solution scheme is shown, allowing also an error in input matrix Y. In the following, another application of the minimum-phase property is proposed. Today, many surface multiple removal algorithms are in use. Most of them belong to the data-driven SRME family, based on the feedback model (Fig. 1a). Although different in-house versions show different performances, at the end of the day all multiple scattering solutions must obey the minimum-phase condition (10a). This does not only apply to the output data $${{\bf P}}_0^ -$$ but also to the interpolated data at the missing data points of input data P − (Appendix  A). Let us assume that the performance of a specific multiple removal algorithm need be objectively tested. Then input-output product $$[{{\bf I}} - {{\bf \hat{Y}}}_0^{}][{{\bf I}} + {{\bf \hat{Y}}}]$$ can be easily computed (product of two data matrices). If we exclude the local data around the origin (‘offset = 0, time = 0’), each column of this product represents an error shot record $$({{{\bf E}}_0}{\vec{\rm I}_{\rm j}})$$. The smaller the samples in the error records, the better the performance of the algorithm. Expression (12b) shows that the error records can be used to improve the solution (see Appendix  A):   $${({{\bf P}}_0^ - {{{\bf A}}^ \cap })^{\rm (i)}} = {({{\bf P}}_0^ - {{{\bf A}}^ \cap })^{({\rm i} - 1)}} + {{\bf E}}_0^{({\rm i} - 1)}[{{\bf I}} - {({{\bf P}}_0^ - {{{\bf A}}^ \cap })^{({\rm i} - 1)}}]$$ (13a)  $${({{{\bf P}}^ - }{{{\bf A}}^ \cap })^{(\rm i)}} = {({{{\bf P}}^ - }{{{\bf A}}^ \cap })^{({\rm i} - 1)}} + [{{\bf I}} + {({{\bf P}}_{}^ - {{{\bf A}}^ \cap })^{({\rm i} - 1)}}{{{\bf E}}^{({\rm i} - 1)}}],$$ (13b)where i = 1 represents the first improvement step. For algorithms with a very good performance, one step may be already sufficient. Note that $${{{\bf P}}_0}{{{\bf A}}^ \cap } = {{\bf X}}_0^ \cup {{{\bf R}}^ \cap }$$. Hence, by using knowledge of the surface reflectivity the IR matrix $${{\bf X}}_0^ \cup$$ is known (of course within the constraint of the limited seismic bandwidth). If we apply the Fourier transform with respect to the source and detector coordinates (meaning multiplication by Fourier matrix F and its Hermitian version FH), minimum-phase condition (10a) becomes a condition for a plane wave response.   $$[{{\bf I}} - ({{\bf F}}{{\bf P}}_0^ - {{{\bf F}}^{\rm H}})({{\bf F}}{{{\bf A}}^ \cap }{{{\bf F}}^{\rm H}})][{{\bf I}} + ({{\bf F}}{{{\bf P}}^ - }{{{\bf F}}^{\rm H}})({{\bf F}}{{{\bf A}}^ \cap }{{{\bf F}}^{\rm H}})] = {{\bf I}}$$ (14a)or   $$[{{\bf I}} - {{\bf \tilde{P}}}_0^ - {{{\bf \tilde{A}}}^ \cap }][{{\bf I}} + {{{\bf \tilde{P}}}^ - }{{{\bf \tilde{A}}}^ \cap }] = {{\bf I}}.$$ (14b)For a 1-D subsurface the data matrices $${{\bf \tilde{P}}}_{}^ - {\rm and}\, {{\bf \tilde{P}}}_0^ -$$ become diagonal matrices, each diagonal element representing one frequency component of a plane wave response, and minimum phase condition (14b) can be simplified to the minimum-phase condition of a time-series for each plane wave (Berkhout 1971, 1974). Including the source ghost in surface multiple scattering In marine data we deal with prominent ghosts at the source and the detector side:   $${{\bf S}}({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = {{{\bf W}}^*}({\rm z}_0^ + ,{\rm z}_{\rm s}^{}){{{\bf S}}^ + }({{\rm z}_{\rm s}}) + {{{\bf R}}^ \cap }({\rm z}_0^ + ,{\rm z}_0^ + ){{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm s}^{}){{{\bf S}}^ - }({{\rm z}_{\rm s}})$$ (15a)  $${{\bf D}}({{\rm z}_{\rm d}};{\rm z}_0^ + ) = {{\bf D}}_{}^ - ({{\rm z}_{\rm d}}){{\bf W}}_{}^*({{\rm z}_{\rm d}},{\rm z}_0^ + ) + {{\bf D}}_{}^ + ({{\rm z}_{\rm d}}){{\bf W}}({{\rm z}_{\rm d}},{\rm z}_0^ + ){{{\bf R}}^ \cap }({\rm z}_0^ + ,{\rm z}_0^ + ),$$ (15b)leading to generalization of expression (2a): $${{\bf P}}({{\rm z}_{\rm d}};{{\rm z}_{\rm s}}) = {{{\bf D}}^ - }({{\rm z}_{\rm d}};{\rm z}_0^ + ){{\bf X}}_{}^ \cup ({\rm z}_0^ + ,{\rm z}_0^ + ){{\bf S}}({\rm z}_0^ + ;{{\rm z}_{\rm s}})$$. It is common practice to apply detector deghosting to the data. Looking at expressions (15b), this implies ‘division’ by the operator $${{\bf D}}({{\rm z}_{\rm d}};{\rm z}_0^ + )$$. Although it is for the theoretical framework not required, let us assume that this detector deghosting process has been successfully applied. For source deghosting, the situation may be much more complex due to the nonlinear interference of the source with the surface during firing time. It leads to different values of operator R∩ in (15a) and (15b), meaning that reciprocity does not hold anymore. An alternative solution is proposed by utilizing the source ghost (R∩W − S − ) for extra illumination (extended double illumination), by including it in the Q-matrix (see Fig. 1b):   \begin{eqnarray} {{{\bf Q}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) &=& {{{\bf W}}^*}({\rm z}_0^ + ,{\rm z}_{\rm s}^{}){{{\bf S}}^ + }({{\rm z}_{\rm s}}) + {{{\bf R}}^ \cap }({\rm z}_0^ + ,{\rm z}_0^ + )[{{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm s}^{}){{{\bf S}}^ - }({{\rm z}_{\rm s}}) + {{{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}})] \nonumber\\ &=& {{{\bf W}}^*}({\rm z}_0^ + ,{\rm z}_{\rm s}^{}){{{\bf S}}^ + }({{\rm z}_{\rm s}}) + {{{\bf R}}^ \cap }({\rm z}_0^ + ,{\rm z}_0^ + ){{{\underline{\bf P}}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) \nonumber \\ & =& {{{\bf S}}^ + }({\rm z}_0^ + ;{\rm z}_{\rm s}^{}) + \delta {{{\underline{\bf P}}}^ + }({\rm z}_0^ + ;{\rm z}_{\rm s}^{}) \end{eqnarray} (16)with $${{\underline{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = {{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm s}^{}){{{\bf S}}^ - }({{\rm z}_{\rm s}}) + {{{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}})$$ and $$\delta {{\underline{\bf P}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = {{{\bf R}}^ \cap }({\rm z}_0^ + ,{\rm z}_0^ + ){{\underline{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}})$$. The rationale for including this option is that the ghost source and the multiple source are both coded by the surface reflectivity operator R∩. It is also interesting to realize that the expression of the data with surface-related multiples and source ghost response (‘surface multiple scattering’) is given by (see Fig. 1b):   $${{\bf P}}_{}^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = {{\bf P}}_0^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) + {{\bf P}}_0^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}){{{\bf A}}^ \cap }({{\rm z}_{\rm s}};{\rm z}_0^ + ){\underline{\bf P }}_{}^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}})$$ (17a)with   $${{\underline{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = {{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm s}^{}){{{\bf S}}^ - }({{\rm z}_{\rm s}}) + {{{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}}).$$ (17b)Note that the extended multiple removal process removes simultaneously the multiples of both the direct source and the ghost source: output $${{\bf P}}_0^ -$$ is multiple-free and ghost-free. Imaging with surface multiple scattering The expression of the upgoing wavefields that arrive at the surface $$({\rm z}_0^ + )$$ is given by:   \begin{eqnarray}{{\underline{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = {{\bf X}}_0^ \cup ({\rm z}_0^ + ,{\rm z}_0^ + ){{{\bf Q}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) + {{{\bf W}}^ - }({\rm z}_0^ + ,{{\rm z}_{\rm s}}){{{\bf S}}^ - }({{\rm z}_{\rm s}}),\end{eqnarray} (18a)where   \begin{eqnarray} {{{\bf Q}}^ + }({\rm z}_0^ + ;{\rm z}_{\rm s}^{}) &=& {{{\bf S}}^ + }({\rm z}_0^ + ;{\rm z}_{\rm s}^{}) + {{{\bf R}}^ \cap }({\rm z}_0^ + ,{\rm z}_0^ + ){{{\underline{\bf P}}}^ - }({\rm z}_0^ + ;{\rm z}_{\rm s}^{})\nonumber\\ &=& {{{\bf S}}^ + }({\rm z}_0^ + ;{\rm z}_{\rm s}^{})[{{\bf I}} + {{{\bf A}}^ \cap }({\rm z}_{\rm s}^{},{\rm z}_0^ + ){{{\underline{\bf P}}}^ - }({\rm z}_0^ + ;{\rm z}_{\rm s}^{})] \end{eqnarray} (18b)represents extended double illumination. If we neglect the internal multiples for the moment, we may write for the Earth's first-order transfer matrix (WRW-model):   $${{\bf X}}_0^ \cup ({\rm z}_0^ + ,{\rm z}_0^ + ) = \sum\limits_{\rm n = 1}^{{\rm m} - 1} {[{{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm n}^ - )} {{{\bf R}}^ \cup }({\rm z}_{\rm n}^ - ,{\rm z}_{\rm n}^ - ){{\bf W}}({\rm z}_{\rm n}^ - ,{\rm z}_0^ + )] + {{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm m}^ - ){\underline{\bf {X}}}_{}^ \cup ({\rm z}_{\rm m}^ + ,{\rm z}_{\rm m}^ + ){{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_0^ + ),$$ (18c)where matrix W equals the wavefield propagation matrix and where one column of reflectivity matrix R∪ equals the angle-dependent reflection operator for the downgoing wavefields at one gridpoint of depth level zm. Let us look at perfect data volumes first. Imaging with surface multiple scattering for perfect data In the seismic migration algorithm the downward travelling wavefields (Q + ) that leave $${\rm z}_0^ +$$ are numerically transported into the subsurface by forward extrapolation. It yields the full illuminating wavefields (P + ) that arrive at the upper part of depth level $${\rm z}_{\rm m}^{}$$, being indicated by $${\rm z}_{\rm m}^ -$$ (m = 1, 2, ……). Again, if we neglect the internal multiples for the moment (see Fig. 2):   \begin{eqnarray} {{{\bf P}}^ +}({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}) &=& {{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_0^ + ){{{\bf Q}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) \nonumber \\ &=& {{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_0^ + ){{{\bf S}}^ + }({\rm z}_0^ + ;{{\rm z}_{\rm s}})[{\bf I} + {{{\bf A}}^ \cap }({\rm z}_{\rm s}^{},{\rm z}_0^ + ){{{\underline{\bf P}}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}})] \nonumber \\ &=& {{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_{\rm s}^{}){{{\bf S}}^ + }({{\rm z}_{\rm s}})[{\bf I} + {{{\bf A}}^ \cap }({\rm z}_{\rm s}^{},{\rm z}_0^ + ){{{\underline{\bf P}}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}})]\nonumber \\ &=& {{\bf P}}_0^ + ({\rm z}_{\rm m}^ - ;{\rm z}_{\rm s}^ + )[{\bf I} + {{{\bf A}}^ \cap }({\rm z}_{\rm s}^{},{\rm z}_0^ + ){{{\underline{\bf P}}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}})], \end{eqnarray} (19a)where $${{\bf P}}_0^ + = {{\bf W}}_{}^ + {{\bf S}}_{}^ +$$ represents the illuminating wavefield without surface ghosts and without multiples. Similarly, the subsurface response—being represented by the upward travelling wavefields P − arriving at $${\rm z}_0^ +$$—are numerically transported into the subsurface by inverse extrapolation, yielding the response wavefields Q − that leave $${\rm z}_{\rm m}^ -$$ (m = 1, 2, ……):   \begin{eqnarray} {{{\bf Q}}^ - }({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}) &=& {{{\bf W}}^*}({\rm z}_{\rm m}^ - ,{\rm z}_0^ + )\left[ {{{{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) - \sum\limits_{\rm n = 1}^{{\rm m} - 1} {{{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm n}^ - )} {{\bf R}}_{}^ \cup ({\rm z}_{\rm n}^ - ,{\rm z}_{\rm n}^ - ){{{\bf P}}^ + }({\rm z}_{\rm n}^ - ;{{\rm z}_{\rm s}})} \right] \nonumber\\ &=& {{{\bf W}}^*}({\rm z}_{\rm m}^ - ,{\rm z}_0^ + )\left[ {{{\bf P}}_0^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) - \sum\limits_{\rm n = 1}^{{\rm m} - 1} {{{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm n}^ - )} {{\bf R}}_{}^ \cup ({\rm z}_{\rm n}^ - ,{\rm z}_{\rm n}^ - ){{\bf P}}_0^ + ({\rm z}_{\rm n}^ - ;{{\rm z}_{\rm s}})} \right]\left[ {{{\bf I}} + {{\bf A}}_{}^ \cap ({\rm z}_{\rm s}^{},{\rm z}_0^ + ){{{\underline{\bf P}}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}})} \right] \nonumber\\ &=& {{\bf Q}}_0^ - ({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})\left[ {{{\bf I}} + {{\bf A}}_{}^ \cap ({\rm z}_{\rm s}^{},{\rm z}_0^ + ){{{\underline{\bf P}}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}})} \right], \end{eqnarray} (19b)where $${{\bf Q}}_0^ -$$ represents the response wavefields at $${\rm z}_{\rm m}^ -$$without surface ghosts and without multiples, the term ∑W − R∪P + represents the response wavefields from shallower depth levels $${\rm z}_{\rm n}^ -$$ (zn < zm) that need be removed and $${{\bf P}}_0^ + = {{{\bf W}}^ + }{{{\bf S}}^ + }$$ being the direct source wavefields that arrive at $${\rm z}_{\rm n}^ -$$. Note that subtraction term ∑W − R∪P + is needed to keep Q − causal. At level $${\rm z}_{\rm m}^ -$$ full response wavefields Q − and full incident wavefields P + are related according to (see Fig. 2):   $${{{\bf Q}}^ - }({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}) = {{{\underline{\bf X}} }}_{}^ \cup ({\rm z}_{\rm m}^ - ,{\rm z}_{\rm m}^ - ){{{\bf P}}^ + }({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}).$$ (20a) Figure 2. View largeDownload slide Feedback model for the wavefields at depth level $${\rm z}_{\rm m}^ -$$ for the simplified situation that the internal multiples in the overburden can be neglected. Note that in practice, surface operator R∩ may be a time-variant operator due to nonlinearity at the source and a changing sea state. Figure 2. View largeDownload slide Feedback model for the wavefields at depth level $${\rm z}_{\rm m}^ -$$ for the simplified situation that the internal multiples in the overburden can be neglected. Note that in practice, surface operator R∩ may be a time-variant operator due to nonlinearity at the source and a changing sea state. Applying direct inversion and using expressions (19a) and (19b), transfer operator $${\underline{\bf {X}}}_{}^ \cup$$ is given by:   \begin{eqnarray} {\underline{\bf {X}}}_{}^ \cup ({\rm z}_{\rm m}^ - ,{\rm z}_{\rm m}^ - ) &=& {{\bf Q}}_{}^ - ({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}){[{{\bf P}}_{}^ + ({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}})]^{ - 1}}\nonumber\\ &=& {{\bf Q}}_0^ - ({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}){[{{\bf P}}_0^ + ({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}})]^{ - 1}}. \end{eqnarray} (20b)Expression (20b) reveals that if we deal with perfect data (no noise, no missing data), direct inversion can be applied and the result is insensitive to the surface multiple scattering (source ghost response + the surface-related multiples). In other words, with perfect data the surface multiple scattering does not give an extra contribution to the estimate of transfer operator $${\underline{\bf {X}}}_{}^ \cup$$. Note that $${\underline{\bf {X}}}_{}^ \cup$$ contains the reflectivity at $${\rm z}_{\rm m}^ -$$, i.e., $${\underline{\bf {X}}}_{}^ \cup = {{\bf R}}_{}^ \cup + {{\bf X}}_{}^ \cup$$. This means that for perfect data R∪ can be completely determined by surface-free response $${{\bf P}}_0^ -$$. If we apply Fourier transformation with respect to the detector coordinates on depth level $${\rm z}_{\rm m}^ -$$ (multiplication by Fourier matrix F to expression (20b) on the detector side), reflection coefficient matrix $${{{\bf \tilde{R}}}^ \cup }$$ can be found at zero intercept time:   $${{\bf \tilde{R}}}_{}^ \cup ({\rm z}_{\rm m}^ - ,{\rm z}_{\rm m}^ - ) = {{\bf F}}({\rm z}_{\rm m}^ - ){{\bf Q}}_{}^ - ({\rm z}_{\rm m}^ - ;{{\rm z}_0}){[{{\bf P}}_{}^ + ({\rm z}_{\rm m}^ - ;{{\rm z}_0})]^{ - 1}} = {{\bf F}}({\rm z}_{\rm m}^ - ){{\bf Q}}_0^ - ({\rm z}_{\rm m}^ - ;{{\rm z}_0}){[{{\bf P}}_0^ + ({\rm z}_{\rm m}^ - ;{{\rm z}_0})]^{ - 1}} \quad {\rm at} \quad \tau = 0,$$ (20c)where $${{\bf \tilde{R}}}_{}^ \cup = {{\bf F}}{{\bf R}}_{}^ \cup$$. Bear in mind that the columns of $${{{\bf \tilde{R}}}^ \cup }$$ represent the angle-dependent plane-wave reflection coefficients at each gridpoint of depth level $${\rm z}_{\rm m}^ -$$. Note also that expression (20c) is consistent with expression (20b) and confirms that in the situation of perfect data surface multiple scattering does not contribute to the redatuming result and to the image. Next, we will focus on the realistic situation of sparse data volumes. Imaging with surface multiple scattering for sparse data In the situation of sparse data, direct inversion combined with applying the imaging principle is replaced by ‘reflectivity estimation by minimization’ (m = 1, 2, …….):   $${\rm E}({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}) = \sum\limits_j {\sum\limits_\omega {{{\left\| {[{{\bf Q}}_{}^ - ({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}) - {{\bf R}}_{}^ \cup ({\rm z}_{\rm m}^ - ,{\rm z}_{\rm m}^ - ){{\bf P}}_{}^ + ({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}})]{{\vec{{\rm I}}}_{\rm j}}} \right\|}^2}} } = {\rm minimum},$$ (21a)where j indicates the shot number, $${\vec{\rm I}_{\rm j}}$$ selects the right column,   $${{{\bf Q}}^ - }({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}) = {{{\bf W}}^*}({\rm z}_{\rm m}^ - ,{\rm z}_0^ + )\left[ {{{{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) - \sum\limits_{\rm n = 1}^{{\rm m} - 1} {{{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm n}^ - )} {{\bf R}}_{}^ \cup ({\rm z}_{\rm n}^ - ,{\rm z}_{\rm n}^ - ){{{\bf P}}^ + }({\rm z}_{\rm n}^ - ;{{\rm z}_{\rm s}})} \right]$$ (21b)represents the causal upgoing wavefields that leave depth level $${\rm z}_{\rm m}^ -$$ and   $${{\bf P}}_{}^ + ({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}) = {{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_{\rm s}^{}){{\bf S}}_{}^ + ({{\rm z}_{\rm s}}) + {{\bf W}}({\rm z}_{\rm m}^ - ;{\rm z}_0^ + )[{{\bf R}}_{}^ \cap ({\rm z}_0^ + ,{\rm z}_0^ + ){\underline{\bf P}}_{}^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}})]$$ (21c)represents the extended double illuminating wavefields that arrive at depth level $${\rm z}_{\rm m}^ -$$. Note that (21a) physically means ‘optimized double focusing’ (Berkhout 1982). Using expression (21c), minimization by (21a) can be refined to (notation is simplified by omitting the depth levels):   $$E = \sum\limits_{\rm j} {\sum\limits_\omega {{{\left\| {[{{\bf Q}}_{}^ - - {{\bf R}}_{\rm a}^ \cup {{{\bf W}}^ + }{{\bf S}}_{}^ + - {{\bf R}}_{\rm b}^ \cup {{{\bf W}}^ + }\delta {\underline{\bf P}}_{}^ + ]{{\vec{{\rm I}}}_{\rm j}}} \right\|}^2}} } = {\rm minimum\,at}\ {\rm z}_{\rm m}^ - ,$$ (22)where $$\delta {{\underline{\bf {P} }}^ + } = {{{\bf R}}^ \cap }{{\underline{\bf {P} }}^ - }$$ at $${\rm z}_0^ +$$ and $${{\bf R}}_{\rm a}^ \cup {{{\bf W}}^ + }{{\bf S}}_{}^ + + {{\bf R}}_{\rm b}^ \cup {{{\bf W}}^ + }\delta {\underline{\bf {P} }}_{}^ + = {{\bf R}}_{}^ \cup {{\bf P}}_{}^ +$$ represents extended wide-angle reflection process at depth level $${\rm z}_{\rm m}^ -$$. From minimization expression (22) we see the benefit of reflectivity estimation by utilizing the surface multiple scattering: Realistic data volumes, represented by parse data matrices can be handled (from direct inversion to iterative minimization); Images can be computed per shot record, or per receiver gather, meaning that only one domain need be properly sampled (impossible with data-driven multiple removal); By utilizing surface multiple scattering two images are obtained (‘image decomposition’), one from the primary source response $$({\rm yielding}\,\,{{\bf R}}_{\rm a}^ \cup )$$ and one from the secondary source response $$({\rm yielding}\,\,{{\bf R}}_{\rm b}^ \cup )$$. Note the significant advantage: instead of applying an elaborate removal process, the ghost source response and the surface multiples are utilized. From the above we may conclude that surface multiple scattering (source ghost response + surface multiples) should not be removed but it should be utilized. In addition, we may conclude that image decomposition$$({{\bf R}}_{}^ \cup\, {\rm into}\,{{\bf R}}_{\rm a}^ \cup ,{{\bf R}}_{\rm b}^ \cup )$$ also allows us to separate primaries and surface multiples by using the WRW-model (‘wavefield decomposition’):   $$$${{\bf P}}_{\rm pr}^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = \sum\limits_{\rm m = 1}^{\rm M} {{{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm m}^ - ){{\bf R}}_{\rm a}^ \cup ({\rm z}_{\rm m}^ - ,{\rm z}_{\rm m}^ - ){{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_{\rm s}^{}){{\bf S}}_{}^ + ({{\rm z}_{\rm s}})}\quad \left( {\rm primaries} \right)$$$$ (23a)  $${{\bf M}}_{\rm sfc}^ - ({\rm z}_0^ + ;{{\rm z}_0}) = \sum\limits_{\rm m = 1}^{\rm M} {{{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm m}^ - ){{\bf R}}_{\rm b}^ \cup ({\rm z}_{\rm m}^ - ,{\rm z}_{\rm m}^ - ){{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_0^ + )\delta {\underline{\bf {P} }}_{}^ + ({\rm z}_0^ + ;{{\rm z}_0})}\quad \left( {\rm surface\, multiples} \right)$$ (23b)with $$\delta {{\underline{\bf {P} }}^ + } = {{{\bf R}}^ \cap }({{{\bf W}}^ - }{{{\bf S}}^ - } + {{{\bf P}}^ - })$$. Note also that all wavefield components include the missing data as well (wavefield reconstruction) and superposition leads to the full response: $${{\bf P}}_{}^ - = {{\bf P}}_{\rm pr}^ - + {{\bf M}}_{\rm sfc}^ -$$ at the surface, where $${{\bf M}}_{\rm sfc}^ -$$ includes the primaries from the ghost source. Of course, both images $$({{\bf R}}_{\rm a}^ \cup ,{{\bf R}}_{\rm b}^ \cup )$$ can be combined into one image $$({{\bf R}}_{}^ \cup )$$ at each subsurface gridpoint k by a data-driven illumination-weighted addition at each subsurface gridpoint k:   $$\vec{\rm R}_{}^ \cup ({\rm k}) = [{\rm A}_{\rm k}^2\vec{\rm R}_{\rm a}^ \cup ({\rm k}) + {\rm B}_{\rm k}^2\vec{R}_{\rm b}^ \cup ({\rm k})]/[{\rm A}_{\rm k}^2 + {\rm B}_{\rm k}^2]$$ (24)with $${\rm A}_{\rm k}^2 = \sum\nolimits_{\rm j} {{\rm a}_{\rm kj}^2}\, {\rm and}\, {\rm B}_{\rm k}^2 = \sum\nolimits_{\rm j} {\rm b_{\rm kj}^2}$$, where akj equals the measured energy of the primary incident wavefield at gridpoint k and $${\rm b_{\rm kj}}$$ equals the measured energy of the secondary incident wavefield at gridpoint k, both being initially generated by primary source j. Using the prior information that reflectivity is an angle-dependent scalar at each depth level (R does not contain traveltimes), wavefield separation and wavefield reconstruction should be an integral part of the closed-loop migration algorithm. In the following we will see that the same conclusion holds for the internal multiple generating system. Feedback model below the surface $$({\rm z}_1^ + ,{\rm z}_2^ + ,{\rm z}_3^ + ,........)$$ Fig. 3(a) visualizes the feedback model at $${\rm z}_{\rm m}^ +$$. Note that IR matrix $${{\underline{\bf {X}} }}_{}^ \cap$$ contains all multiples in the overburden $$(z < {\rm z}_{\rm m}^ + )$$. If we would neglect the internal multiples, we may write:   $$\delta {\underline{\bf {X}}}_{}^ \cap ({\rm z}_{\rm m}^ + ,{\rm z}_{\rm m}^ + ) = {{\bf W}}({\rm z}_{\rm m}^ + ,{\rm z}_0^ + ){{\bf R}}_{}^ \cap ({\rm z}_0^ + ,{\rm z}_0^ + ){{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm m}^ + ),$$ (25a)leading to Fig. 2. An interesting refinement of expression (25a) is given by:   $${\underline{\bf {X}}}_{\rm pr}^ \cap ({\rm z}_{\rm m}^ + ,{\rm z}_{\rm m}^ + ) = \sum\limits_{\rm n = 0}^{\rm m} {[{{\bf W}}({\rm z}_{\rm m}^ + ,{\rm z}_{\rm n}^ + ){{\bf R}}_{}^ \cap ({\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ + ){{\bf W}}({\rm z}_{\rm n}^ + ,{\rm z}_{\rm m}^ + )]} ,$$ (25b)being the WRW-model of $${{\underline{\bf {X}}}^ \cap }$$ (linearity in R∩). It generates the third-order multiples in the data. In the following we will use all orders of multiples. At $${\rm z}_{\rm m}^ +$$ the feedback equations can be written as:   \begin{eqnarray} {{{\bf P}}^ - }({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}}) &=& {{\bf X}}_0^ \cup ({\rm z}_{\rm m}^ + ,{\rm z}_{\rm m}^ + ){{{\bf Q}}^ + }({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}}) \nonumber\\ &=& {{\bf X}}_0^ \cup ({\rm z}_{\rm m}^ + ,{\rm z}_{\rm m}^ + )[{{\bf W}}({\rm z}_{\rm m}^ + ,{{\rm z}_{\rm s}}){{{\bf S}}^ + }({{\rm z}_{\rm s}}) + {{{\underline{\bf {X} }}}^ \cap }({\rm z}_{\rm m}^ + ,{\rm z}_{\rm m}^ + ){{{\bf P}}^ - }({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})] \nonumber\\ &=& {{\bf X}}_0^ \cup ({\rm z}_{\rm m}^ + ,{\rm z}_{\rm m}^ + ){{\bf Q}}_0^ + ({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})[{{\bf I}} + {{{\bf A}}^ \cap }({{\rm z}_{\rm s}},{\rm z}_{\rm m}^ + ){{{\bf P}}^ - }({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})] \nonumber\\ &=& {{\bf P}}_0^ - ({{\rm z}_0};{{\rm z}_{\rm s}})[{{\bf I}} + {{{\bf A}}^ \cap }({{\rm z}_{\rm s}},{\rm z}_{\rm m}^ + ){{{\bf P}}^ - }({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})], \end{eqnarray} (26a)where (m = 1, 2, ………):   $${{\bf Q}}_0^ + ({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}}) = {{\bf W}}({\rm z}_{\rm m}^ + ,{{\rm z}_{\rm s}}){{\bf S}}_{}^ + ({{\rm z}_{\rm s}}),$$ (26b)  $${{\bf A}}_{}^ \cap ({{\rm z}_{\rm s}},{\rm z}_{\rm m}^ + ) = {[{{\bf W}}({\rm z}_{\rm m}^ + ,{{\rm z}_{\rm s}}){{{\bf S}}^ + }({{\rm z}_{\rm s}})]^{ - 1}}{{\underline{\bf {X} }}^ \cap }({\rm z}_{\rm m}^ + ,{\rm z}_{\rm m}^ + )$$ (26c)and where $${{\underline{\bf {X} }}^ \cap }({\rm z}_{\rm m}^ + ,{\rm z}_{\rm m}^ + )$$ equals the full up-down IR matrix of the overburden $$(z \le {\rm z}_{\rm m}^ + )$$ that includes all internal multiples. Similar to the expression of the down-up IR matrix $${{\underline{\bf {X} }}^ \cup } = {{{\bf R}}^ \cup } + {{{\bf X}}^ \cup }$$ at $${\rm z}_{\rm m}^ -$$, we may write for the expression of the up-down IR matrix at $${\rm z}_{\rm m}^ +$$: $${{\underline{\bf {X} }}^ \cap } = {{{\bf R}}^ \cap } + {{{\bf X}}^ \cap }$$. We will use this property to image R∩. From expression (26a) it follows that:   $${{{\bf Q}}^ + }({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}}) = {{\bf Q}}_0^ + ({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})[{{\bf I}} + {{{\bf A}}^ \cap }({{\rm z}_{\rm s}},{\rm z}_{\rm m}^ + ){{{\bf P}}^ - }({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})].$$ (26d) Figure 3. View largeDownload slide (a) Feedback model in the subsurface (m > 0), showing the generation of the internal multiples in the data. Note again the important role of multiple generation operator [I + A∩P − ] at $${\rm z}_{\rm m}^ +$$. Compare with Fig. 1(a). (b) Extension of (a), showing also the propagation operators in layer $$({\rm z}_{\rm m}^ + ,{\rm z}_{{\rm m} + 1}^ - )$$. Compare with Fig. 2. Figure 3. View largeDownload slide (a) Feedback model in the subsurface (m > 0), showing the generation of the internal multiples in the data. Note again the important role of multiple generation operator [I + A∩P − ] at $${\rm z}_{\rm m}^ +$$. Compare with Fig. 1(a). (b) Extension of (a), showing also the propagation operators in layer $$({\rm z}_{\rm m}^ + ,{\rm z}_{{\rm m} + 1}^ - )$$. Compare with Fig. 2. Note the resemblance between eqs (3a)–(3d) and (26a)–(26d). Similar to what we saw at the surface $$({\rm z}_0^ + )$$, expression (26a) shows that $$[{{\bf I}} + {{{\bf A}}^ \cap }{{\bf P}}_{}^ - ]$$ is the overburden multiple generation operator at $${\rm z}_{\rm m}^ +$$, transforming $${{\bf P}}_0^ -$$ into $${{\bf P}}_{}^ -$$ and $${{\bf Q}}_0^ +$$ into $${{\bf Q}}_{}^ +$$. From expression (26a) it also directly follows that (omitting the depth level indication in the notation):   $${{\bf P}}_0^ - = {{{\bf P}}^ - }{[{{\bf I}} + {{{\bf A}}^ \cap }{{{\bf P}}^ - }]^{ - 1}}\quad {\rm at} \quad{\rm z}_{\rm m}^ +$$ (27a)  $${{\bf Q}}_0^ + = {{{\bf Q}}^ + }{[{{\bf I}} + {{{\bf A}}^ \cap }{{{\bf P}}^ - }]^{ - 1}} \quad {\rm at} \quad {\rm z}_{\rm m}^ + ,$$ (27b)showing that [I + A∩P − ] − 1 equals the overburden multiple removal operator, transforming $${{\bf P}}_{}^ -$$ into $${{\bf P}}_0^ -$$ and $${{\bf Q}}_{}^ +$$ into $${{\bf Q}}_0^ +$$ at $${\rm z}_{\rm m}^ +$$. Considering the large similarity between the feedback model at the surface $$({\rm z}_0^ + )$$ and in the subsurface $$({\rm z}_{\rm m}^ + )$$, we can derive in exactly the same way ($${\underline{\bf {X} }}_{}^ \cap$$ plays the same role as R∩ and W + S + replaces S + ):   $${[{{\bf I}} + {{{\bf P}}^ - }{{{\bf A}}^ \cap }]^{ - 1}} = [{{\bf I}} - {{\bf P}}_0^ - {{{\bf A}}^ \cap }]\quad {\rm at} \quad {\rm z}_{\rm m}^ +$$ (28a)or   $${[{{\bf I}} - {{\bf P}}_0^ - {{{\bf A}}^ \cap }]^{ - 1}}[{{\bf I}} + {{{\bf P}}^ - }{{{\bf A}}^ \cap }] = {{\bf I}}\quad {\rm at} \quad {\rm z}_{\rm m}^ +$$ (28b)where   {{\bf A}}_{}^ \cap = $${{\bf A}}_{}^ \cap ({{\rm z}_{\rm s}},{\rm z}_{\rm m}^ + ) = {[{{\bf W}}({\rm z}_{\rm m}^ + ,{{\rm z}_{\rm s}}){{{\bf S}}^ + }({{\rm z}_{\rm s}})]^{ - 1}}{{\underline{\bf {X} }}^ \cap }({\rm z}_{\rm m}^ + ,{\rm z}_{\rm m}^ + )$$. (28c) Eq. (28a) shows that not only at the surface (m = 0) but also in the subsurface (m = 1, 2, …. .) the multiple-generating operator [I + P − A∩] has the minimum-phase property. This means that at each depth level the wavefields can be validated by testing whether they fulfil the minimum-phase condition (28b) and, optionally, the error shot records can be visualized at each depth level as a QC for the full wavefield imaging system. We consider the minimum-phase based validation and updating capability a new opportunity to improve our migration algorithms on field data in an objective way. Imaging with all multiple scattering events for perfect data If we make use of expressions (27a) and (27b) and assume perfect data (no noise, no missing data), then direct inversion leads to (m = 1, 2, …….):   $${{\bf X}}_0^ \cup ({\rm z}_{\rm m}^ + ,{\rm z}_{\rm m}^ + ) = {{{\bf P}}^ - }({\rm z}_{\rm m}^ + ;{\rm z}_{\rm s}^{}){[{{{\bf Q}}^ + }({\rm z}_{\rm m}^ + ;{\rm z}_{\rm s}^{})]^{ - 1}} = {{\bf P}}_0^ - ({\rm z}_{\rm m}^ + ;{\rm z}_{\rm s}^{}){[{{\bf Q}}_0^ + ({\rm z}_{\rm m}^ + ;{\rm z}_{\rm s}^{})]^{ - 1}}.$$ (29) Hence, we see that overburden operator [I + A∩P − ] drops out and, therefore, overburden multiple scattering (surface-related + internal) does not contribute to redatuming result $${{\bf X}}_0^ \cup ({\rm z}_{\rm m}^ + ,{\rm z}_{\rm m}^ + )$$. This means that under the ideal condition of perfect data the combination of removing all multiples from the overburden (surface and internal) followed by traditional redatuming (‘first-order redatuming’) will reveal the desired transfer operator $${{\bf X}}_0^ \cup$$. It also explains why the concept of Marchenko redatuming (Broggini & Snieder 2012) is valid in the situation of perfect data (Wapenaar et al.2014). Using expression (29) and moving inside layer $$({\rm z}_{\rm m}^ + ,{\rm z}_{{\rm m} + 1}^ - )$$ from $${\rm z}_{\rm m}^ +$$ to $${\rm z}_{{\rm m} + 1}^ -$$ by wavefield extrapolation, then we obtain (see Fig. 3b):   \begin{eqnarray} {{\bf Q}}_{}^ - ({\rm z}_{{\rm m} + 1}^ - ;{{\rm z}_{\rm s}}) &=& {{{\bf W}}^*}({\rm z}_{{\rm m} + 1}^ - ,{\rm z}_{\rm m}^ + ){{\bf P}}_{}^ - ({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}}) \nonumber\\ &=& {{{\bf W}}^*}({\rm z}_{{\rm m} + 1}^ - ,{\rm z}_{\rm m}^ + ){{\bf P}}_0^ - ({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})[{{\bf I}} + {{{\bf A}}^ \cap }({\rm z}_{\rm s}^{},{\rm z}_{\rm m}^ + ){{\bf P}}_{}^ - ({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})] \nonumber\\ &=& {{\bf Q}}_0^ - ({\rm z}_{{\rm m} + 1}^ + ;{{\rm z}_{\rm s}})[{{\bf I}} + {{{\bf A}}^ \cap }({\rm z}_{\rm s}^{},{\rm z}_{\rm m}^ + ){{\bf P}}_{}^ - ({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})], \end{eqnarray} (30a)  \begin{eqnarray} {{\bf P}}_{}^ + ({\rm z}_{{\rm m} + 1}^ - ;{{\rm z}_{\rm s}}) &=& {{\bf W}}({\rm z}_{{\rm m} + 1}^ - ,{\rm z}_{\rm m}^ + ){{\bf Q}}_{}^ + ({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}}) \nonumber\\ &=& {{\bf W}}({\rm z}_{{\rm m} + 1}^ - ,{\rm z}_{\rm m}^ + ){{\bf Q}}_0^ + ({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})[{{\bf I}} + {{{\bf A}}^ \cap }({\rm z}_{\rm s}^{},{\rm z}_{\rm m}^ + ){{\bf P}}_{}^ - ({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})] \nonumber\\ &=& {{\bf P}}_0^ + ({\rm z}_{{\rm m} + 1}^ + ;{{\rm z}_{\rm s}})[{{\bf I}} + {{{\bf A}}^ \cap }({\rm z}_{\rm s}^{},{\rm z}_{\rm m}^ + ){{\bf P}}_{}^ - ({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})]. \end{eqnarray} (30b) Using (30a) and (30b), expression (29) can be extended to:   $${{\bf {X} }}_{}^ \cup ({\rm z}_{{\rm m} + 1}^ - ,{\rm z}_{{\rm m} + 1}^ - ) = {{\bf Q}}_{}^ - ({\rm z}_{{\rm m} + 1}^ - ;{\rm z}_{\rm s}^{}){[{{\bf P}}_{}^ + ({\rm z}_{{\rm m} + 1}^ - ;{\rm z}_{\rm s}^{})]^{ - 1}} = {{\bf Q}}_0^ - ({\rm z}^-_{\rm m+1}{\rm z}_{\rm s}^{}){[{{\bf P}}_0^ + ({\rm z}_{\rm m+1}^-{\rm z}_{\rm s}^{})]^{ - 1}}.$$ (31a) Taking into account that $${{\underline{\bf {X} }}^ \cup } = {{{\bf X}}^ \cup } + {{{\bf R}}^ \cup }$$, we may write (compare with expression 20c):   $${{\bf \tilde{R}}}_{}^ \cup ({\rm z}_{{\rm m} + 1}^ - ,{\rm z}_{{\rm m} + 1}^ - ) = {{\bf F}}({\rm z}_{{\rm m} + 1}^ - ){{\bf Q}}_0^ - ({\rm z}_{{\rm m} + 1}^ - ;{\rm z}_{\rm s}^{}){[{{\bf P}}_0^ + ({\rm z}_{{\rm m} + 1}^ - ;{\rm z}_{\rm s}^{})]^{ - 1}}\quad {\rm at} \quad \tau = 0,$$ (31b)yielding the angle-dependent reflection coefficients at each gridpoint of depth level $${\rm z}_{{\rm m} + 1}^ -$$. Expression (31b) is consistent with redatuming output (29): ‘in the ideal situation of perfect data multiples drop out and they do not contribute to the image’. This means that in the situation of perfect data volumes, it is justified to follow the traditional approach and apply multiple removal prior to imaging. It also explains why the concept of direct inversion works well in the theoretical situation of perfect data (Weglein et al. 2010). Imaging with all multiple scattering events for sparse data As was already mentioned at the start, in practical situations the theoretical situation of perfect data does not exist. Similar to what we saw when dealing with surface multiple scattering, in practical situations (i) multiple scattering should not be removed but utilized and (ii) direct inversion must be replaced by minimization. Now, if we look again at expressions (21a)–(21d) and we include internal multiples, then these expressions can be extended to (see Fig. 4):   $$E({\rm z}_{{\rm m} + 1}^ - ;{{\rm z}_{\rm s}}) = \sum\limits_{\rm j} {\sum\limits_\omega {{{\left\| {[{{\bf Q}}_{}^ - ({\rm z}_{{\rm m} + 1}^ - ;{{\rm z}_{\rm s}}) - {{\bf R}}_{}^ \cup ({\rm z}_{{\rm m} + 1}^ - ,{\rm z}_{{\rm m} + 1}^ - ){{\bf P}}_{}^ + ({\rm z}_{{\rm m} + 1}^ - ;{{\rm z}_{\rm s}})]{{\vec{{\rm I}}}_{\rm j}}} \right\|}^{2}}} } = {\rm minimum},$$ (32a)where   $${{{\bf Q}}^ - }({\rm z}_{{\rm m} + 1}^ - ;{{\rm z}_{\rm s}}) = {{{\bf W}}^*}({\rm z}_{{\rm m} + 1}^ - ,{\rm z}_0^ + )\left[ {{{{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}}) - \sum\limits_{\ell = 1}^{\rm m} {{{\bf W}}({\rm z}_0^ + ,{\rm z}_\ell ^ - )} {{\bf R}}_{}^ \cup ({\rm z}_\ell ^ - ,{\rm z}_\ell ^ - ){{{\bf P}}^ + }({\rm z}_\ell ^ - ;{{\rm z}_{\rm s}})} \right]$$ (32b)represents the upgoing wavefields that leave depth level $${\rm z}_{{\rm m} + 1}^ -$$,   $${{\bf P}}_{}^ + ({\rm z}_{{\rm m} + 1}^ - ;{{\rm z}_{\rm s}}) = {{\bf W}}({\rm z}_{{\rm m} + 1}^ - ,{\rm z}_{\rm s}^{}){{\bf S}}_{}^ + ({{\rm z}_{\rm s}}) + {{\bf W}}({\rm z}_{{\rm m} + 1}^ - ,{\rm z}_0^ + )\delta {\underline{\bf {P} }}_{}^ + ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) + \sum\limits_{\rm n = 1}^{\rm m} {{{\bf W}}({\rm z}_{{\rm m} + 1}^ - ,{\rm z}_{\rm n}^ + )\delta {{\bf P}}_{}^ + ({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}})}$$ (32c)represents the triple downgoing wavefields that arrive at depth level $${\rm z}_{{\rm m} + 1}^ -$$ (triple illumination) and where $$\delta {\underline{\bf P}}_{}^ + ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = {{\bf R}}_{}^ \cap ({\rm z}_0^ + ,{\rm z}_0^ + ){\underline{\bf {P} }}_{}^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}})\, {\rm and}\, \delta {{\bf P}}_{\!\!\!\!\!\!-}^ + ({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}}) = {{\bf R}}_{}^ \cap ({\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ + ){{\bf P}}_{}^ - ({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}})$$ represent the secondary illuminating wavefields that leave the surface $${\rm z}_0^ +$$ and depth level $${\rm z}_{\rm n}^ +$$ (n = 1, 2, …. . m) respectively. Simplifying the notation by omitting the depth levels, minimization expression (32a) may also be written as (m = 1, 2, ……):   $${{\rm E}_{{\rm m + 1}}} = \sum\limits_{\rm j} {\sum\limits_\omega {{{\left\| {\left[ {{{\bf Q}}_{}^ - - {{\bf R}}_{\rm a}^ \cup {{{\bf W}}^ + }{{\bf S}}_{}^ + - {{\bf R}}_{\rm b}^ \cup {{{\bf W}}^ + }\delta {\underline{\bf {P} }}_{}^ + - {{\bf R}}_{\rm c}^ \cup \sum\limits_{\rm n = 1}^{\rm m} {{{{\bf W}}^ + }\delta {{\bf P}}_{}^ + } } \right]{{\vec{{\rm I}}}_{\rm j}}} \right\|}^{2}}} } = {\rm minimum\, at}\ {\rm z}_{\rm m}^ - .$$ (33) Figure 4. View largeDownload slide Wavefield propagation in the subsurface, showing the generation of multiple scattering (surface + internal) in the downgoing and upgoing wavefields during propagation. It explains how full wavefield extrapolation need be applied. In our full wavefield migration algorithm (FWM) this is implemented in a closed-loop. Bear in mind that if the source ghost is included in the surface multiple scattering, we need to take $${{{\bf R}}^ \cap }{{\underline{\bf {P} }}^ - }$$ for n = 0, where $${{\underline{\bf {P} }}^ - } = {{{\bf W}}^ - }{{{\bf S}}^ - } + {{{\bf P}}^ - }$$ at the surface. Figure 4. View largeDownload slide Wavefield propagation in the subsurface, showing the generation of multiple scattering (surface + internal) in the downgoing and upgoing wavefields during propagation. It explains how full wavefield extrapolation need be applied. In our full wavefield migration algorithm (FWM) this is implemented in a closed-loop. Bear in mind that if the source ghost is included in the surface multiple scattering, we need to take $${{{\bf R}}^ \cap }{{\underline{\bf {P} }}^ - }$$ for n = 0, where $${{\underline{\bf {P} }}^ - } = {{{\bf W}}^ - }{{{\bf S}}^ - } + {{{\bf P}}^ - }$$ at the surface. From expression (33) we see the benefit of triple reflectivity estimation by minimization: Realistic data volumes, represented by sparse data matrices, can be handled; Images can be computed per shot record, or per receiver gather, meaning that only one domain need be properly sampled; By using all multiple scattering (surface and internal) three images are obtained (‘image decomposition’), one from the primary response $$({\rm yielding }\,{{\bf R}}_{\rm a}^ \cup )$$, one from the surface multiple scattering response $$({\rm yielding }\,{{\bf R}}_{\rm b}^ \cup )$$ and one from the internal multiple scattering response $$({\rm yielding }\,{{\bf R}}_{\rm c}^ \cup )$$. Note again the significant advantage: instead of an elaborate removal process, the surface and internal multiple scattering are utilized. From the above we may conclude that image decomposition $$({{\bf R}}_{}^ \cup\,{\rm into}\,{{\bf R}}_{\rm a}^ \cup ,{{\bf R}}_{\rm b}^ \cup ,{{\bf R}}_{\rm c}^ \cup )$$ also allows us to separate primary, surface-related and internal multiple wavefields (see Fig. 4):   $$$${{\bf P}}_{\rm pr}^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = \sum\limits_{\rm m = 1}^{\rm M} {{{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm m}^ - ){{\bf R}}_{\rm a}^ \cup ({\rm z}_{\rm m}^ - ,{\rm z}_{\rm m}^ - ){{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_{\rm s}^{}){{\bf S}}_{}^ + ({{\rm z}_{\rm s}})}\quad \left( {\rm primaries} \right)$$$$ (34a)  $${{\bf M}}_{\rm sfc}^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = \sum\limits_{\rm m = 1}^{\rm M} {{{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm m}^ - ){{\bf R}}_{\rm b}^ \cup ({\rm z}_{\rm m}^ - ,{\rm z}_{\rm m}^ - ){{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_0^ + )\delta {\underline{\bf {P} }}_{}^ + ({\rm z}_0^ + ;{{\rm z}_{\rm s}})}\quad \left( {\rm surface\, multiples} \right)$$ (34b)  $${{\bf M}}_{{\mathop{\rm int}} }^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = \sum\limits_{\rm m = 2}^{\rm M} {{{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm m}^ - ){{\bf R}}_{\rm c}^ \cup ({\rm z}_{\rm m}^ - ,{\rm z}_{\rm m}^ - )\sum\limits_{\rm n = 1}^{{\rm m} - 1} {{{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_{\rm n}^ + )\delta {{\bf P}}_{}^ + ({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}})} } .\quad\left( {\rm internal\, multiples} \right).$$ (34c) These wavefields are computed in each iteration by FWMod (from depth image to time measurements, see Fig. 6). Note also that all decomposed wavefields include the missing data as well (wavefield reconstruction) and superposition leads to the full response $${{\bf P}}_{}^ - = {{\bf P}}_{\rm pr}^ - + {{\bf M}}_{\rm sfc}^ - + {{\bf M}}_{{\mathop{\rm int}} }^ -$$ (compare with expressions 23a and 23b). This is an important conclusion: using the prior information that reflectivity is an angle-dependent scalar at each depth level, ultimate wavefield separation and wavefield reconstruction should not only be part of pre-processing but it also should be an integral part of the migration algorithm. Similar to expressions (24), the combined image can be determined by a data-driven illumination- weighted addition at each subsurface gridpoint k:   $$\vec{\rm R}_{}^ \cup ({\rm k}) = [{\rm A}_{\rm k}^2\vec{\rm R}_{\rm a}^ \cup ({\rm k}) + {\rm B}_{\rm k}^2\vec{R}_{\rm b}^ \cup ({\rm k}) + {\rm C}_{\rm k}^2\vec{\rm R}_{\rm c}^ \cup ({\rm k})]/[{\rm A}_{\rm k}^2 + {\rm B}_{\rm k}^2 + {\rm C}_{\rm k}^2]$$ (35)with $${\rm A}_{\rm k}^2 = \sum\nolimits_{\rm j} {{\rm a}_{\rm kj}^2} ,{\rm B}_{\rm k}^2 = \sum\nolimits_{\rm j} {\rm b_{\rm kj}^2}\, {\rm and}\,{\rm C}_{\rm k}^2 = \sum\nolimits_{\rm j} {{\rm c}_{\rm kj}^2}$$, where akj equals the measured energy of the primary incident wavefield at gridpoint k, $${\rm b_{\rm kj}}$$ equals the measured energy of the secondary incident wavefield at gridpoint k and ckj equals the measured energy of the higher-order incident wavefield at gridpoint k, all being initially generated by primary source j. Using multiple scattering to image the lower side of boundaries as well In the foregoing we have explained how to image the upper side of boundaries by utilizing primary and multiple scattering, yielding the reflectivity operator R∪ at each depth level. If we want to image the lower side of the boundaries as well, the utilization of internal multiple scattering is a must. Having already determined R∪, we have very good knowledge of $${{\bf P}}_{\rm pr}^ -$$ and $${{\bf M}}_{\rm sfc}^ -$$ at $${\rm z}_0^ +$$ (see expressions 34a and 34b). This means that we also have very good knowledge of the internal multiples in the measurements: $${{\bf P}}_{}^ - - {{\bf P}}_{\rm pr}^ - - {{\bf M}}_{\rm sfc}^ -$$ at $${\rm z}_0^ +$$. The expression of $${{\bf M}}_{{\mathop{\rm int}} }^ -$$ is given by (see Fig. 5):   $${{\bf M}}_{{\mathop{\rm int}} }^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = \sum\limits_{\rm m = 2}^{\rm M} {\delta {{\bf M}}_{{\mathop{\rm int}} }^ - ({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}})} = \sum\limits_{\rm m = 2}^{\rm M} {{{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm m}^ - ){{\bf R}}_{\rm c}^ \cup ({\rm z}_{\rm m}^ - ,{\rm z}_{\rm m}^ - ){{\bf M}}_{{\mathop{\rm int}} }^ + ({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}})}$$ (36a)with   $${{\bf M}}_{{\mathop{\rm int}} }^ + ({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}) = \sum\limits_{\rm n = 1}^{{\rm M} - 1} {{{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_{\rm n}^ + ){{\bf R}}_{}^ \cap ({\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ + ){{\bf P}}_{}^ - ({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}})}$$ (36b)or, changing the summation sequence,   $${{\bf M}}_{{\mathop{\rm int}} }^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = \sum\limits_{\rm n = 1}^{{\rm M} - 1} {{{\bf X}}_{\rm pr}^ \cup ({\rm z}_0^ + ,{\rm z}_{\rm n}^ + ){{\bf R}}_{}^ \cap ({\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ + ){{\bf P}}_{}^ - ({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}})}$$ (36c)with   $${{\bf X}}_{\rm pr}^ - ({\rm z}_0^ + ;{\rm z}_{\rm n}^ + ) = \sum\limits_{\rm m = n + 1}^{\rm M} {{{\bf W}}({\rm z}_0^ + ,{\rm z}_{\rm m}^ - ){{\bf R}}_{\rm c}^ \cup ({\rm z}_{\rm m}^ - ,{\rm z}_{\rm m}^ - ){{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_{\rm n}^ - )} .$$ (36d) Figure 5. View largeDownload slide Grouping of the internal multiples $$(\delta {{\bf M}}_{{\mathop{\rm int}} }^ - )$$ according to generation by upper-side reflectivity $${{{\bf R}}^ \cup }({\rm z}_{\rm m}^ - ,{\rm z}_{\rm m}^ - )$$, see left-hand side, and by lower-side reflectivity $${{{\bf R}}^ \cap }({\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ + )$$, see right-hand side. Figure 5. View largeDownload slide Grouping of the internal multiples $$(\delta {{\bf M}}_{{\mathop{\rm int}} }^ - )$$ according to generation by upper-side reflectivity $${{{\bf R}}^ \cup }({\rm z}_{\rm m}^ - ,{\rm z}_{\rm m}^ - )$$, see left-hand side, and by lower-side reflectivity $${{{\bf R}}^ \cap }({\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ + )$$, see right-hand side. We will use expressions (36c) and (36d) for the estimation of lower-side reflectivity $${{{\bf R}}^ \cap }({\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ + )$$, n = 1, 2, ………M − 1. If we substitute the three wavefield components of incident wavefield $${{\bf P}}_{}^ -$$ in (36c), we obtain:   $${{\bf M}}_{{\mathop{\rm int}} }^ - ({\rm z}_0^ + ;{{\rm z}_{\rm s}}) = \sum\limits_{\rm n = 1}^{{\rm M} - 1} {{{\bf X}}_{\rm pr}^ \cup ({\rm z}_0^ + ,{\rm z}_{\rm n}^ + ){{\bf R}}_{}^ \cap ({\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ + )[{{\bf P}}_{\rm pr}^ - ({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}})} + {{\bf M}}_{\rm sfc}^ - ({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}}) + {{\bf M}}_{{\mathop{\rm int}} }^ - ({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}})]$$ (37a)where wavefields $${{\bf P}}_{\rm pr}^ -$$, $${{\bf M}}_{\rm sfc}^ -$$ and $${{\bf M}}_{{\mathop{\rm int}} }^ -$$ at $${\rm z}_{\rm n}^ +$$ are known and computed by FWMod:   $${{\bf P}}_{\rm pr}^ - ({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}}) = \sum\limits_{\ell = {\rm n + 1}}^{\rm M} {{{\bf W}}({\rm z}_{\rm n}^ + ,{\rm z}_\ell ^ - ){{\bf R}}_{\rm a}^ \cup ({\rm z}_\ell ^ - ,{\rm z}_\ell ^ - ){{\bf W}}({\rm z}_\ell ^ - ,{\rm z}_{\rm s}^{}){{\bf S}}_{}^ + ({{\rm z}_{\rm s}})}$$ (37b)  $${{\bf M}}_{\rm sfc}^ - ({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}}) = \sum\limits_{\ell = {\rm n + 1}}^{\rm M} {{{\bf W}}({\rm z}_{\rm n}^ + ,{\rm z}_\ell ^ - ){{\bf R}}_{\rm b}^ \cup ({\rm z}_\ell ^ - ,{\rm z}_\ell ^ - ){{\bf W}}({\rm z}_\ell ^ - ,{\rm z}_0^{})\delta {{\bf {P} }}_{}^ + ({\rm z}_0^{};{{\rm z}_{\rm s}})}$$ (37c)  $${{\bf M}}_{{\mathop{\rm int}} }^ - ({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}}) = \sum\limits_{\ell = {\rm n + 1}}^{\rm M} {{{\bf W}}({\rm z}_{\rm n}^ + ,{\rm z}_\ell ^ - ){{\bf R}}_{\rm c}^ \cup ({\rm z}_\ell ^ - ,{\rm z}_\ell ^ - )\sum\limits_{\ell = 1}^{{\rm m} - 1} {{{\bf W}}({\rm z}_\ell ^ - ,{\rm z}_{\rm n}^ + )\delta {{\bf P}}_{}^ + ({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}})} } .$$ (37d) Expressions (37a)–(37d) lead to minimization equation (omitting the depth levels):   $$\sum\limits_{\rm j} {\sum\limits_\omega {{{\left\| {\left[ {[{{\bf P}}_{}^ - - {{\bf P}}_{\rm pr}^ - - {{\bf M}}_{\rm sfc}^ - ] - {{\bf X}}_{\rm pr}^ \cup {{\bf R}}_{\rm a}^ \cap {{\bf P}}_{\rm pr}^ - - {{\bf X}}_{\rm pr}^ \cup {{\bf R}}_{\rm b}^ \cap {{\bf M}}_{\rm sfc}^ - - {{\bf X}}_{\rm pr}^ \cup {{\bf R}}_{\rm c}^ \cap {{\bf M}}_{{\mathop{\rm int}} }^ - } \right]{{\vec{{\rm I}}}_{\rm j}}} \right\|}^2}} } = {\rm minimum\, at}\,{\rm z}_{\rm n}^ + .$$ (38) Compare with minimization (33). Note that in minimization (38) the measurements in $${{\bf P}}_{}^ -$$ are known, the reflectivity operators $${{\bf R}}_{\rm a}^ \cap$$, $${{\bf R}}_{\rm b}^ \cap$$ and $${{\bf R}}_{\rm c}^ \cap$$ are unknown and $${{\bf P}}_{\rm pr}^ -$$, $${{\bf M}}_{\rm sfc}^ -$$, $${{\bf M}}_{{\mathop{\rm int}} }^ -$$ are given by expressions (37b)–(37d). It is important to realize that R∩ = −R∪ in areas where wave conversion can be neglected (see Appendix  B). This is the start of iterative minimization (38). Closed-loop, full-wavefield imaging scheme for sparse data Fig. 6 shows the schematic of closed-loop, full wavefield migration (FWM). Bear in mind that in many practical situations already some kind of estimate of the image is available (think not only in terms of previous time lapse images but think also at output already available from mainstream seismic processing); this knowledge should be optionally used as input in the first iteration of FWM. In the feedback loop full wavefield modeling (FWMod) transforms the (initial) image into full wavefield measurements (‘image transformation’). These measurements are presented in terms of primaries, surface multiples and internal multiples according to expressions (34a)–(34c). This (initial) wavefield decomposition is used in the construction of three separate images (image decomposition): first-order image, second-order image and higher-order image. This is input for the next iteration. As a last step, these three images may be combined into one final image according to expression (35). In the downgoing path the reflectivities are computed. In the upgoing path the different images are transformed into different wavefields. These wavefields are compared with the measured wavefields, allowing (i) a scaling correction of the internal wavefields, (ii) filling in the missing data, (iii) estimating the direct and ghost source properties and (iv) updating the surface reflectivity by refining the detector deghosting process (marine data). Roundtrips are continued until simulated wavefields and measured wavefields are sufficiently close (Berkhout 2014). Figure 6. View largeDownload slide In closed-loop, full-wavefield migration (FWM) the seismic image is transformed into simulated measurements (FWMod) and they are compared with the real measurements. After minimization the error measurements are input for the next iteration. Decomposed wavefields and decomposed images are essential constituents in the algorithm. Figure 6. View largeDownload slide In closed-loop, full-wavefield migration (FWM) the seismic image is transformed into simulated measurements (FWMod) and they are compared with the real measurements. After minimization the error measurements are input for the next iteration. Decomposed wavefields and decomposed images are essential constituents in the algorithm. Numerical illustration To illustrate the closed-loop, full wavefield imaging theory we consider a simple subsurface model with three angle-dependent reflectors. The simplicity of the model allows us to identify all primary and multiple-scattering wavefields as they occur in the theory. Fig. 7 shows the model with the angle-dependent reflectivity operators (R∪ at $${\rm z}_1^ - ,{\rm z}_2^ - ,{\rm z}_3^ -$$ and R∩ at $${\rm z}_1^ + ,{\rm z}_2^ +$$). Figure 7. View largeDownload slide Subsurface model with three angle-dependent reflectors. The reflection operators R∪ and R∩ for each reflector are displayed within the angle range of the seismic data, the dashed line representing the imaginary values (which are estimated as well). The seven source locations in this experiment are indicated. Detectors are situated over the full range (9000 m) with a sampling interval of 15 m. Figure 7. View largeDownload slide Subsurface model with three angle-dependent reflectors. The reflection operators R∪ and R∩ for each reflector are displayed within the angle range of the seismic data, the dashed line representing the imaginary values (which are estimated as well). The seven source locations in this experiment are indicated. Detectors are situated over the full range (9000 m) with a sampling interval of 15 m. Figs 8(a)–(c) visualize the wavefields that play a principal role in the today's mainstream seismic processing packages: Fig. 8(a) shows the downgoing primary and secondary source wavefields that leave the surface $${\rm z}_0^ +$$, being one column (j = 4) of   $${{{\bf Q}}^ + } = {{{\bf S}}^ + } + {{{\bf R}}^ \cap }{{{\bf P}}^ - },$$ (39a)where S + represents the impulsive primary source wavefields and R∩P − represents the highly dispersive secondary source wavefields. It is assumed that the ghost wavefields have already been removed. Fig. 8(b) shows the upgoing reflected wavefields (PP) that arrive at the surface $${\rm z}_0^ +$$, being one column (j = 4) of   $${{\bf P}}_{}^ - = {{\bf P}}_0^ - + {{\bf M}}_0^ - ,$$ (39b)where $${{\bf P}}_0^ -$$ represents the surface-free response $$({{\bf X}}_0^ \cup {{{\bf S}}^ + })$$ and $${{\bf M}}_0^ -$$ represents the surface-related multiples $$({{\bf X}}_0^ \cup {{{\bf R}}^ \cap }{{{\bf P}}^ - })$$. The separation is realized by a surface-related multiple removal algorithm. Fig. 8(c) shows the arriving wavefields (PP) at $${\rm z}_0^ +$$ that are generated by S + only, being one column (j = 4) of   $${{\bf P}}_0^ - = {{\bf P}}_{\rm pr}^ - + {\mathop{\bf M}^{\frown }} {}_{{\rm int}}^ - ,$$ (39c)where $${{\bf P}}_{\rm pr}^ -$$ represents the primary response $$({{\bf X}}_{\rm pr}^ \cup {{{\bf S}}^ + })$$ and $$\displaystyle{\mathop{\bf M}^{\frown }} {}_{{\rm int}}^ -$$ represents the internal multiples in $${{\bf P}}_0^ -$$. The separation is realized by an internal multiple removal algorithm. Figure 8. View largeDownload slide Wavefields at $${\rm z}_0^ +$$ (generated by the central source j = 4) that play a principal role in today's mainstream seismic processing packages. $$\displaystyle{\mathop{\bf M}^{\frown }}{}_{{\rm int}}^ - {\vec{\rm I}_{\rm j}}$$ contains the internal multiples being generated by primary source $${{\bf S}}_0^ + {\vec{\rm I}_{\rm j}}$$. The internal multiples being generated by secondary source $${{{\bf R}}^ \cap }{{{\bf P}}^ - }{\vec{\rm I}_{\rm j}}$$ are part of the surface-related multiples $${{\bf M}}_0^ - {\vec{\rm I}_{\rm j}}$$. Figure 8. View largeDownload slide Wavefields at $${\rm z}_0^ +$$ (generated by the central source j = 4) that play a principal role in today's mainstream seismic processing packages. $$\displaystyle{\mathop{\bf M}^{\frown }}{}_{{\rm int}}^ - {\vec{\rm I}_{\rm j}}$$ contains the internal multiples being generated by primary source $${{\bf S}}_0^ + {\vec{\rm I}_{\rm j}}$$. The internal multiples being generated by secondary source $${{{\bf R}}^ \cap }{{{\bf P}}^ - }{\vec{\rm I}_{\rm j}}$$ are part of the surface-related multiples $${{\bf M}}_0^ - {\vec{\rm I}_{\rm j}}$$. The primary wavefield $${{\bf P}}_{\rm pr}^ -$$ is used as input in first-order migration (today's standard). $$\displaystyle{\mathop{\bf M}^{\frown }}{}_{{\rm int}}^ -$$ represents the group of internal multiples that is generated by the primary sources (columns of S + ). The other group of internal multiples $$(\displaystyle{\mathop{\bf M}^{\frown^{\!\!\!\!\!\!\frown} }}{}_{{\mathop{\rm int}} }^ - )$$ can be found in $${{\bf M}}_0^ -$$ and is generated by the secondary sources (columns of R∩P − ) at the surface: $$\displaystyle{\mathop{\bf M}^{\frown^{\!\!\!\!\!\!\frown}}}{}_{{\mathop{\rm int}} }^ - = \displaystyle{\mathop{\bf M}^{\frown}}{}_{{\mathop{\rm int}} }^ - {{{\bf A}}^ \cap }{{\bf P}}_{}^ -$$ with $${{{\bf A}}^ \cap } = {[{{{\bf S}}^ + }]^{ - 1}}{{\bf R}}_{}^ \cap$$. It is interesting to see that $$\displaystyle{\mathop{\bf M}^{\frown}}{}_{{\mathop{\rm int}} }^ -$$ is the prediction operator that determines $$\displaystyle{\mathop{\bf M}^{\frown^{\!\!\!\!\!\!\frown} }}{}_{{\mathop{\rm int}} }^ -$$ from the input data $${{\bf P}}_{}^ -$$. Note here the complete similarity with $${{\bf P}}_0^ -$$ functioning as the prediction operator for the determination of surface-related multiples $${{\bf M}}_0^ -$$: $${{\bf M}}_0^ - = {{\bf P}}_0^ - {{{\bf A}}^ \cap }{{\bf P}}_{}^ -$$. Next, Figs 9(a)–(c) visualize the wavefields that play a principal role in closed-loop, multiple-scattering imaging: Fig. 9(a) shows the downgoing primary and secondary source wavefields that leave the surface $${\rm z}_0^ +$$, being one column (j = 4) of   $${{{\bf Q}}^ + } = {{{\bf W}}^*}{{{\bf S}}^ + } + {{{\bf R}}^ \cap }({{{\bf W}}^ - }{{{\bf S}}^ - } + {{{\bf P}}^ - }),$$ (40a)where the primary sources $${{\bf S}}_{}^ +$$ are situated at $${\rm z}_{\rm s}^{} = 50\,{\rm m}$$, where $${{\bf W}}_{}^*$$ transforms the downgoing wavefield of these sources to $${\rm z}_0^ +$$ and where the ghost sources $${{{\bf R}}^ \cap }({{\bf W}}_{}^ - {{\bf S}}_{}^ - )$$ are included in the surface multiple scattering $$({{{\bf R}}^ \cap }{\underline{\bf {P} }}^ - )$$. Fig. 9(b) shows the upgoing reflected wavefields (PP) that arrive at the surface $${\rm z}_0^ +$$, being one column (j = 4) of   $${{\bf P}}_{}^ - = {{\bf P}}_{\rm pr}^ - + {{\bf M}}_{}^ - ,$$ (40b)where $${{\bf P}}_{\rm pr}^ - = {{\bf X}}_{\rm pr}^ \cup {{{\bf S}}^ + }$$ and $${{\bf M}}_{}^ -$$ contains all multiples (surface +internal). Fig. 9(c) shows the upgoing surface- and internal-multiple scattering wavefields that arrive at $${\rm z}_0^ +$$, being one column (j = 4) of   $${{\bf M}}_{}^ - = {{\bf M}}_{\rm sfc}^ - + {{\bf M}}_{{\mathop{\rm int}} }^ - ,$$ (40c)where $${{\bf M}}_{\rm sfc}^ - = {{\bf X}}_{\rm pr}^ \cup ({{{\bf R}}^ \cap }{{\bf {P} }}_{\!\!\!\!\!\!-}^ - )$$ contains all surface multiples and $${{\bf M}}_{{\mathop{\rm int}} }^ -$$ contains all internal multiples: $$( \displaystyle{\mathop{\bf M}^{\frown}}{}_{{\mathop{\rm int}} }^ - + \displaystyle{\mathop{\bf M}^{\frown}}{}_{{\mathop{\rm int}} }^ - {{{\bf A}}^ \cap }{{\bf {P} }}_{\!\!\!\!\!\!-}^ - )$$. Figure 9. View largeDownload slide Wavefields at $${\rm z}_0^ +$$ (generated by the central source j = 4) that play a principal role in closed-loop, full wavefield imaging. Note that $${{\bf M}}_{}^ - {\vec{\rm I}_{\rm j}}$$ contains all multiples (surface + internal), $${{\bf M}}_{\rm sfc}^ - {\vec{\rm I}_{\rm j}}$$ contains only the surface multiples and $${{\bf M}}_{{\mathop{\rm int}} }^ - {\vec{\rm I}_{\rm j}}$$ contains all internal multiples $$(\displaystyle{\mathop{\bf M}^{\frown}}{}_{{\mathop{\rm int}}}^ - {\vec{\rm I}_{\rm j}} + \displaystyle{\mathop{\bf M}^{\frown^{\!\!\!\!\!\!\frown} }}{}_{{\mathop{\rm int}}}^ - {\vec{\rm I}_{\rm j}})$$. Figure 9. View largeDownload slide Wavefields at $${\rm z}_0^ +$$ (generated by the central source j = 4) that play a principal role in closed-loop, full wavefield imaging. Note that $${{\bf M}}_{}^ - {\vec{\rm I}_{\rm j}}$$ contains all multiples (surface + internal), $${{\bf M}}_{\rm sfc}^ - {\vec{\rm I}_{\rm j}}$$ contains only the surface multiples and $${{\bf M}}_{{\mathop{\rm int}} }^ - {\vec{\rm I}_{\rm j}}$$ contains all internal multiples $$(\displaystyle{\mathop{\bf M}^{\frown}}{}_{{\mathop{\rm int}}}^ - {\vec{\rm I}_{\rm j}} + \displaystyle{\mathop{\bf M}^{\frown^{\!\!\!\!\!\!\frown} }}{}_{{\mathop{\rm int}}}^ - {\vec{\rm I}_{\rm j}})$$. The three response wavefields $${{\bf P}}_{\rm pr}^ - ,{{\bf M}}_{\rm sfc}^ - ,{{\bf M}}_{{\mathop{\rm int}} }^ -$$ are used as input to the closed-loop, multiple-scattering algorithm. Bear in mind that the surface multiples $$({{\bf M}}_{\rm sfc}^ - )$$ do not contain internal multiples, but they do contain the ghost source response ($${{\bf M}}_{\rm sfc}^ - = {{\bf P}}_{\rm pr}^ - {{{\bf A}}^ \cap }{{\bf {P} }}_{\!\!\!\!\!\!-}^ -$$). Next, Figs 10(a)–(c) visualize the decomposed wavefields that enter and leave the upper part of the third reflector, followed by the decomposed images in the ray parameter domain: Fig. 10(a) shows one column (j = 4) of the incident wavefield $${{\bf W}}_{}^ + {{\bf S}}_{}^ +$$ and its response $${{\bf P}}_{\rm pr}^ -$$ at $${\rm z}_3^ -$$. Using all seven sources, first-order image $${{\bf \tilde{R}}}_{\rm a}^ \cup$$ is obtained at $${\rm z}_3^ -$$. Note the acquisition shadow zones in the first-order image due to the coarse sampling of the man-made primary sources (think of a cross section). Fig. 10(b) shows one column (j = 4) of the incident wavefield $${{\bf W}}_{}^ + ({{{\bf R}}^ \cap }{{\bf P}}_{}^ - )$$ and its response $${{\bf M}}_{\rm sfc}^ -$$ at $${\rm z}_3^ -$$. Using all seven sources, second-order image $${{\bf \tilde{R}}}_{\rm b}^ \cup$$ is obtained at $${\rm z}_3^ -$$. Note that in the second-order image the shadow zones have largely disappeared and the illumination area has increased due to the denser and wider illumination by the nature-made secondary sources at the surface. Fig. 10(c) shows one column (j = 4) of the incident wavefield $${{\bf W}}_{}^ + \sum\nolimits_{\rm n = 1}^2 {{{{\bf R}}^ \cap }{{\bf P}}_{}^ - }$$ and its response $${{\bf M}}_{{\mathop{\rm int}} }^ -$$ at $${\rm z}_3^ -$$. Using all seven sources, high-order image $${{\bf \tilde{R}}}_c^ \cup$$ is obtained at $${\rm z}_3^ -$$. Similar to the surface multiples, note the extra contribution of the internal multiples with respect to the primaries. Figure 10. View largeDownload slide Decomposed incident wavefields (left-hand side) with their responses (in the middle) at the upper part of the third reflector, being generated by the central source j = 4. By using all seven sources, the decomposed angle-dependent reflectivities are obtained (right-hand side). Note that the reflectivities are displayed in the ray-parameter domain for each spatial location along the reflector. Figure 10. View largeDownload slide Decomposed incident wavefields (left-hand side) with their responses (in the middle) at the upper part of the third reflector, being generated by the central source j = 4. By using all seven sources, the decomposed angle-dependent reflectivities are obtained (right-hand side). Note that the reflectivities are displayed in the ray-parameter domain for each spatial location along the reflector. In this illustration multiple scattering from dipping boundaries does not occur, minimizing thus the transformation from low to high angles and vice versa. However, we do already see the increase in angle range and, last but not least, we also see that multiples reveal the angle-dependent reflection information. First results from complex examples confirm the expectation that the more complex the geology, the more valuable the contribution from multiple scattering (current research). Finally, Figs 11(a)–(c) visualize the decomposed wavefields that enter and leave the lower part of the second reflector, followed by the decomposed images in the ray parameter domain: Fig. 11(a) shows one column (j = 4) of the incident wavefield $${{\bf P}}_{\rm pr}^ -$$ and its response $$\delta {{\bf P}}_{\rm pr}^ + = {{\bf R}}_{\rm a}^ \cap {{\bf P}}_{\rm pr}^ -$$ at $${\rm z}_2^ +$$. Using all seven sources leads to the angle-dependent image $${{\bf \tilde{R}}}_a^ \cap$$ at $${\rm z}_2^ +$$. Fig. 11(b) shows one column (j = 4) of the incident wavefield $${{\bf M}}_{\rm sfc}^ -$$ and its response $$\delta {{\bf M}}_{\rm sfc}^ + = {{\bf R}}_{\rm b}^ \cap {{\bf M}}_{\rm sfc}^ -$$ at $${\rm z}_2^ +$$. Using all seven sources leads to the angle-dependent image $${{\bf \tilde{R}}}_b^ \cap$$ at $${\rm z}_2^ +$$. Fig. 11(c) shows one column (j = 4) of the incident wavefield $${{\bf M}}_{{\mathop{\rm int}} }^ -$$ and its response $$\delta {{\bf M}}_{{\mathop{\rm int}} }^ + = {{\bf R}}_{\rm c}^ \cap {{\bf M}}_{{\mathop{\rm int}} }^ -$$ at $${\rm z}_2^ +$$. Using all seven sources leads to the angle-dependent image $${{\bf \tilde{R}}}_c^ \cap$$ at $${\rm z}_2^ +$$. Figure 11. View largeDownload slide Decomposed incident wavefields (left-hand side) with their responses (in the middle) at the lower part of the second reflector, being generated by the central source j = 4. By using all seven sources, the decomposed angle-dependent reflectivities are obtained (right-hand side). Note again that the reflectivities are displayed in the ray-parameter domain for each spatial location along the reflector. Figure 11. View largeDownload slide Decomposed incident wavefields (left-hand side) with their responses (in the middle) at the lower part of the second reflector, being generated by the central source j = 4. By using all seven sources, the decomposed angle-dependent reflectivities are obtained (right-hand side). Note again that the reflectivities are displayed in the ray-parameter domain for each spatial location along the reflector. Of course, reflection images of the lower part of reflectors can only be obtained by utilizing internal multiple scattering (see also Fig. 5). Note that in some areas, for instance think of shallow reflectors in the unconsolidated near-surface layer on land, images of the upper part are extremely difficult to acquire in practice. Hence, for these situations new opportunities arise if we aim at images of the lower part ($${{\bf R}}_{}^ \cap$$ instead of $${{\bf R}}_{}^ \cup$$) by using illumination with wavefields that arrive from below: $${{\bf P}}_{}^ - ({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}}) = \sum\limits_{\rm n = m + 1}^{\rm M} {{{\bf W}}({\rm z}_{\rm m}^ + ,{\rm z}_{\rm n}^ - ){{\bf R}}_{\rm c}^ \cup ({\rm z}_{\rm n}^ - ,{\rm z}_{\rm n}^ - ){{\bf P}}_{}^ + ({\rm z}_{\rm n}^ - ;{{\rm z}_{\rm s}})}$$. But also in borehole seismic imaging new opportunities are created. This example illustrates that multiple scattering is not only important to strengthen the response from the upper part of geological boundaries, it is our only hope to image those boundaries which responses from above are not measured at all. CONCLUSIONS In the theoretical situation of perfect data (no noise, no missing data) we have top-quality, fully filled seismic data matrices and primary sources are able to realize complete illumination, even inside geological shadow zones. For such perfect data volumes direct inversion algorithms can be applied to remove surface-related and internal multiple scattering. Using this multiple-free result as input, first-order redatuming and first-order primary migration can provide us with the full image. In practice, however, seismic measurements always contain noise and we never have complete 3-D data volumes at our disposal. We actually deal with sparse data matrices that cause acquisition-generated (non-geological) shadow zones in the subsurface. Without the introduction of assumptions, real data matrices cannot be used in direct inversion and redatuming. Therefore, in practice we badly need to include multiple-scattering wavefields to fill up the first-order illumination gaps that are caused by the sparse measurements—possibly in combination with geological shadow zones—and the background noise. Note that multiple-scattering contributions are ‘free of charge’ and also include the ghost source as extra illumination (when applicable). The author's message is that in the real world multiple scattering is a blessing; it must not be removed but it must be utilized. Instead of applying a complex and elaborate multiple scattering removal process, both surface and internal multiple scattering can be utilized to give an extra contribution to the illuminating wavefield (‘triple illumination’). Internal multiple scattering even allows us to obtain a reflection image from the lower part of reflectors! This is explained by the closed-loop, full-wavefield imaging theory. The full-wavefield extension of the imaging theory also explains that decomposition of the seismic response in first-order (‘primaries’), second-order (‘surface multiples’) and higher-order (‘internal multiples’) scattering should be an integral part of the migration algorithm. In closed-loop FWM the concepts of image decomposition and wavefield decomposition enhance each other (Fig. 6). By comparing the decomposed wavefield components (34a)–(34c) with the measured data, the estimates of the direct and ghost source properties (S + , δS + ) and the surface reflectivity (R∩) can be updated. In addition, the estimates of the missing data can be improved by wavefield reconstruction per decomposed wavefield (part of FWMod). All these capabilities tell us that the closed-loop, full-wavefield migration algorithm should be considered as the nucleus of the future seismic processing system. Ultimately it will position pre-processing as part of an ‘iteration zero’ activity. By preserving the multiple scattering in the seismic measurements, we also preserve a fundamental physical property—the minimum-phase property—in the measurements. This means that an objective way of testing the validity of the estimated wavefields at each depth level can be carried out. It is proposed as an integral part of the closed-loop, full-wavefield imaging process. Epilogue: an image of the future In the history of seismic imaging the presence of strong multiple scattering wavefields has always been perceived as bad news. This explains why a large amount of effort was spent to develop smart algorithms for the attenuation and removal of the ‘damned’ multiples. The author was part of this painstaking endeavour for many years. Today we start to realize that multiple scattering wavefields provide us with extra illumination. There exist billions of different multiples, reaching areas where primaries may have disappeared below the noise (‘shadow zones’) or are not present at all (from the lower part of boundaries). In addition, multiples contain extra information on velocity and absorption, as they travel not only once but many times through the subsurface. The challenge is not to attenuate or remove them, but to find out how to extract subsurface information from these extremely complex phenomena. This article aims at changing the traditional aversion in the seismic industry against multiple scattering by convincing the reader that the presence of multiples means actually good news. This has been done by showing how the conventional theory of second- and higher-order seismic wavefields can be reformulated such that the result leads to a transparent theoretical framework (with simple WRW-modules) and a practical imaging algorithm (with decomposed images). The output consists of angle-dependent reflectivity images that honour wavefield conversion (no acoustic assumption) and include transmission effects (Appendix  B). One important aspect is to organize multiples not according to the well-known sequences of down- and upward traveling (ray)paths, but according to downward-scattered wavefields (R∩P − ) at each depth level (m = 0, 1, 2, …….). This R∩-related organization leads to new theoretical insights that tell us that the very complex multiple scattering generating systems in the Earth have the minimum-phase property. This means that under perfect data conditions multiple-scattering removal systems are causal, allowing the prediction of multiple-scattering from measurements of the past only (in multiple scattering the past determines the future). The minimum-phase property makes clear why algorithms such as direct inversion and Marchenko redatuming can perform well under perfect data conditions (complete sampling, noise-free). Today, seismic data has been recorded in most areas of the world and seismic images are widely available. This situation certainly applies to the many mature exploration and production basins. The multiple-scattering imaging algorithm welcomes such an initial image. It is directly used in the feedback loop of Fig. 6 to make a first estimate of decomposed wavefields down $$({{{\bf W}}^ + }{{\bf S}}_{}^ + ,{{{\bf W}}^ + }\delta {{\bf {P} }}_{\!\!\!\!\!\!-}^ + ,{{\bf M}}_{{\mathop{\rm int}} }^ + )$$ and up $$({{\bf P}}_{\rm pr}^ - ,{{\bf M}}_{\rm sfc}^ - ,{{\bf M}}_{{\mathop{\rm int}} }^ - )$$. Of course, in the special situation that no initial seismic image is available yet, current mainstream processing—with data interpolation, deblending, multiple removal and first-order migration as the main steps—has still the important task to estimate an initial image (‘iteration zero’ in closed-loop full-wavefield imaging). Following this way of thinking, current mainstream processing in the industry would shift to the early exploration phase, thus becoming a pre-imaging tool in virginal areas. According to the author, removal of multiple scattering is more difficult than imaging of multiple scattering. Or, in other words, removal of multiple scattering should not be done before imaging but it should be done as part of imaging. Removal before imaging means prediction and subtraction of multiples in the real-data time domain. The slightest phase error in the prediction result leads already to a large residue, particularly in the higher part of the seismic spectrum. Moreover the data-domain prediction process is very complex, as the prediction operator needs many low-noise shot records that are well sampled in both the source and detector domain. On the other hand, in multiple-scattering imaging the subsurface information in an initial image is utilized to estimate the decomposed wavefields (down as well as up) by image transformation (FWMod). These decomposed wavefields provide ‘triple illumination’ and are again used to construct new decomposed images, etc. (see the closed-loop scheme in Fig. 6). Bear in mind that imaging is a massive double stacking process (focusing in emission + focusing in detection) that is known to be forgiving for initial errors and robust in improving signal-to-noise ratio. Decomposition of wavefields and images may provide new opportunities in the design of (blended) data acquisition geometries. The contribution of all different wavefields to the final image is very nonlinear in reflectivity and may change the outcome of the optimum source and detector distributions. For instance, our current research indicates that the traditional symmetric sampling concept may need be replaced by an asymmetric one where only one domain need be well sampled. This appears even more valid in the situation of complex geology, where limited-angle, first-order wavefields are transformed to high-order wavefields that create incidence from all sides. Finally, I remark about velocities. Traditionally, migration velocities are estimated with the aid of primary wavefields. However, higher-order reflections are more sensitive to velocity errors than primaries. This means that there exist opportunities to extend velocity estimation to multiple scattering wavefields for more accurate velocities. We are working on full-automatic, anisotropic velocity estimation for multiple-scattering imaging. In this ambitious effort, estimation and analysis of propagation operator W plays a key role. Ultimately, we have the ambition to apply seismic data acquisition and multiple-scattering imaging in exactly the same way as self-driving cars collect data and simultaneously image their immediate environment. This approach may revolutionize the way we acquire and image time laps data by considering the sequence of surveys as one big data volume that defines the image of a geological medium with unknown time-variant components. Acknowledgements The author is grateful to the Delphi sponsors for their continuing financial support and for the stimulating discussions on industry needs during the Delphi meetings. In addition, the author would like to thank Dr Eric Verschuur who supported him with the generation of the numerical illustrations. REFERENCES Araujo F.V., Weglein A.B., Carvalho P.M., Stolt R.H., 1994. Inverse scattering series for multiple attenuation: an example with surface and internal multiples, in SEG Expanded Abstracts , pp. 704– 706. Berkhout A.J., 1971. Minimum Phase in Sampled-Signal Theory , TU Delft. Berkhout A.J., 1974. Related properties of minimum-phase and zero-phase time functions, Geophys. Prospect. , 22, 683– 709. https://doi.org/10.1111/j.1365-2478.1974.tb00111.x Google Scholar CrossRef Search ADS   Berkhout A.J., 1982. Seismic Migration: Imaging of Acoustic Energy by Wavefield Extrapolation, A: Theoretical Aspects , 2nd edn, Elsevier. Berkhout A.J., 2014. Review paper: an outlook on the future of seismic imaging, Part II: Full-wavefield migration, Geophys. Prospect. , 62( 5), 931– 949. https://doi.org/10.1111/1365-2478.12154 Google Scholar CrossRef Search ADS   Broggini F., Snieder R., 2012. Connection of scattering principles: a visual and mathematical tour, Eur. J. Phys. , 33, 593– 613. https://doi.org/10.1088/0143-0807/33/3/593 Google Scholar CrossRef Search ADS   Davydenko M., Verschuur D.J., van Groenestijn G.J.A., 2015. Full wavefield migration applied to field data, in 77th Annual International Conference and Exhibition, EAEG, Extended Abstracts , We-N101-10. Verschuur D.J., Berkhout A.J., 1997. Estimation of multiple scattering by iterative inversion, Part II: Practical aspects and examples, Geophysics , 62, 1596– 1611. https://doi.org/10.1190/1.1444262 Google Scholar CrossRef Search ADS   Wapenaar K., Thorbecke J., van der Neut J., Broggini F., Slob E., Snieder R., 2014. Marchenko imaging, Geophysics , 79, WA39– WA57. https://doi.org/10.1190/geo2013-0302.1 Google Scholar CrossRef Search ADS   Weglein A.B., Liu F., Wang Z., Li X., Liang H., 2010. The inverse scattering series depth imaging algorithms: development, tests and progress towards field data application, in SEG, Soc. Expl. Geophys., Expanded abstracts , Denver, pp. 4133– 4138. Weglein A.B., 2014. Multiples: signal or noise?, in SEG Technical Program Expanded Abstracts , pp. 4393– 4399. Weglein A.B., 2015. Primaries: the only events that can be migrated and for which migration has meaning, Leading Edge , 33, 808– 813. https://doi.org/10.1190/tle34070808.1 Google Scholar CrossRef Search ADS   Whitmore N.D., Valenciano A.A., Sollner W., 2010. Imaging of primaries and multiples, using a dual-sensor towed streamer, in 80th Ann. Internat. Mtg, SEG, Expanded Abstracts , pp. 3187– 3192. APPENDIX A: MULTIPLE REMOVAL WITH THE MINIMUM-PHASE CONDITION Using main text expressions (12a) and (12b), we may conclude that at $${\rm z}_0^ +$$:   $$(1)\,\,{{{\bf E}}_0} = [{{\bf I}} - {{\bf \hat{Y}}}_0^{}][{{\bf I}} + {{\bf Y}}] - {{\bf I}} = {{\bf Y}} - {{\bf \hat{Y}}}_0^{} - {{\bf \hat{Y}}}_0^{}{{\bf Y}}$$ (A1a)  $$(2)\,\,\Delta {{\bf Y}}_0^{} = {{{\bf E}}_0}{[{{\bf I}} + {{\bf Y}}]^{ - 1}} \approx {{{\bf E}}_0}[{{\bf I}} - {{{\bf \hat{Y}}}_0}],$$ (A1b)where $${{\bf Y}} = {{\bf P}}_{}^ - {{{\bf A}}^ \cap }\quad {\rm and} \quad{{{\bf Y}}_0} = {{\bf P}}_0^ - {{{\bf A}}^ \cap }$$. Expressions (A1a) and (A1b) lead to the iterative solution at $${\rm z}_0^ +$$:   $$\Delta {{\bf Y}}_0^{(\rm i)} = {{\bf E}}_0^{({\rm i} - 1)}[{{\bf I}} - {{\bf Y}}_0^{({\rm i} - 1)}],$$ (A2a)  $${{\bf \hat{Y}}}_0^{(\rm i)} = {{\bf Y}}_0^{({\rm i} - 1)} + \Delta {{\bf Y}}_0^{(\rm i)},$$ (A2b)  $${{{\bf \hat{E}}}^{(\rm i)}}_0 = {{\bf Y}} - {\alpha ^{(\rm i)}}[{{\bf \hat{Y}}}_0^{(\rm i)} + {({{\bf \hat{Y}}}_0^{}{{\bf Y}})^{(\rm i)}}],$$ (A2c)  $${{\bf Y}}_0^{(\rm i)} = {\alpha ^{(\rm i)}}{{\bf \hat{Y}}}_0^{(\rm i)}\quad {\rm and} \quad{[{{\bf Y}}_0^{}{{\bf Y}}]^{(\rm i)}} = {\alpha ^{(\rm i)}}{[{{\bf \hat{Y}}}_0^{}{{\bf Y}}]^{(\rm i)}}.$$ (A2d) Note that $$[{{\bf Y}}_0^{(\rm i)} + {({{\bf Y}}_0^{}{{\bf Y}})^{(\rm i)}}]$$ is used as interpolated data at the missing data points. If we start in the first iteration (i = 1) with $${{\bf \hat{Y}}}_0^{(0)} = {{\bf Y}}$$, then it follows from expression (A1a) that $${{\bf E}}_0^{(0)} = - {{{\bf Y}}^2}$$ and the first estimates of $$\Delta {{\bf Y}}_0^{(1)}$$ and $${{\bf Y}}_0^{(1)}$$ can be computed with expressions (A2a)–(A2d). The elements in the diagonal matrix α(i) are determined such that the error vectors $$({\vec{\rm E}_{\rm j}} = {{\bf E}}{\vec{\rm I}_{\rm j}})$$ in error matrix E0 have minimum length at the data points:   $$E_0^{(\rm i)} = \sum\limits_{\rm j} {\sum\limits_\omega {{{\left\| {\left[ {{{\bf Y}} - {\alpha ^{(\rm i)}}[{{\bf \hat{Y}}}_0^{(\rm i)} + {{({{\bf \hat{Y}}}_0^{}{{\bf Y}})}^{(\rm i)}}]} \right]{{\vec{{\rm I}}}_{\rm j}}} \right\|}^{2}} = {\rm minimum}} } .$$ (A3) Note that in the last iteration α(i) ≈ I. If we also include errors in the input matrix, $${{\bf Y}} = {{\bf \hat{Y}}} + \Delta {{\bf Y}}$$, expressions (A1a) and (A1b) need be generalized to:   \begin{eqnarray} (1)\,{{{\bf E}}_{\rm tot}} &=& [{{\bf I}} - {{\bf \hat{Y}}}_0^{}][{{\bf I}} + {{\bf \hat{Y}}}] - {{\bf I}} = ({{\bf I}} - {{\bf Y}}_0^{} + \Delta {{\bf Y}}_0^{})({{\bf I}} + {{\bf Y}} - \Delta {{\bf Y}}) - {{\bf I}}\nonumber\\ &=& \Delta {{\bf Y}}_0^{}[{{\bf I}} + {{\bf Y}}] - [{{\bf I}} - {{\bf Y}}_0^{}]\Delta {{\bf Y}}\nonumber\\ & =& {{{\bf E}}_0} + {{\bf E}}. \end{eqnarray} (A4a)  $$(\rm 2a)\, \Delta {{\bf Y}}_0^{} = {{{\bf E}}_0}{[{{\bf I}} + {{\bf \hat{Y}}}]^{ - 1}} \approx {{{\bf E}}_0}[{{\bf I}} - {{{\bf \hat{Y}}}_0}],$$ (A4b)  $$(2{\rm b})\,\Delta {{\bf Y}} = - {[{{\bf I}} - {{{\bf \hat{Y}}}_0}]^{ - 1}}{{\bf E}} \approx - [{{\bf I}} + {{\bf \hat{Y}}}]{{\bf E}}.$$ (A4c) Now, $${{\bf Y}}_0^{}$$ and Y are alternately updated until E0 and E are sufficiently small. It is interesting that in this dual updating scheme both multiple removal and data interpolation are driven by the minimum-phase condition:   $$$$[{{\bf I}} - {{\bf Y}}_0^{}][{{\bf I}} + {{\bf Y}}] = {{\bf I}}.$$$$ (A5) APPENDIX B: ACOUSTIC VERSUS WEAK ELASTIC VERSUS FULL ELASTIC If T represents the transmission operator and we write T = I + δT then δT represents the forward scattering operator. If we explicitly include the transmission operator in our wavefield expressions, the backscattered wavefields $$\delta \vec{\rm P}_{\rm j}^ - = {{{\bf R}}^ \cup }\vec{\rm P}_{\rm j}^ +$$ and $$\delta \vec{\rm P}_{\rm j}^ + = {{{\bf R}}^ \cap }\vec{\rm P}_{\rm j}^ -$$ need be extended to secondary sources (Berkhout 2014):   $$\delta \vec{{\rm S}}_{\rm j}^ - \left( {{\rm z}_{\rm n}^ - ;{{\rm z}_0}} \right)\ = {{{\bf R}}^ \cup }\left( {{\rm z}_{\rm n}^ - ,{\rm z}_{\rm n}^ - } \right){\vec{\rm P}_{\rm j}}^ + \left( {{\rm z}_{\rm n}^ - ;{{\rm z}_0}} \right)\ + \delta {{{\bf T}}^ - }\left( {{\rm z}_{\rm n}^ - ,{\rm z}_{\rm n}^ + } \right){\vec{\rm P}_{\rm j}}^{ - }\left( {{\rm z}_{\rm n}^ + ;{{\rm z}_0}} \right)$$ (B1a)  $$\delta \vec{{\rm S}}_{\rm j}^ + \left( {{\rm z}_{\rm n}^ + ;{\rm z}_0^{}} \right) = {{{\bf R}}^ \cap }\left( {{\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ + } \right)\vec{\rm P}_{\rm j}^ - \left( {{\rm z}_{\rm n}^ + ;{\rm z}_0^{}} \right)\ + \delta {{{\bf T}}^ + }\left( {{\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ - } \right)\vec{\rm P}_{\rm j}^ + \left( {{\rm z}_{\rm n}^ - ;{\rm z}_0^{}} \right),$$ (B1b)representing the elastic interaction equations for P-waves. Note that R and δT play the same role in the scattering process. Eqs (B1a) and (B1b) can be seen as the extension of Huygens’ principle to inhomogeneous, elastic media. Note also that $$\delta \vec{{\rm S}}_{\rm j}^ - = \vec{\rm Q}_{\rm j}^ - - \vec{\rm P}_{\rm j}^ -$$ and $$\delta \vec{{\rm S}}_{\rm j}^ + = \vec{\rm Q}_{\rm j}^ + - \vec{\rm P}_{\rm j}^ +$$, see Fig. B1. Bear in mind that wavefields $${\rm Q}_{\rm j}^ \pm$$ are leaving a depth level and wavefields $${\rm P}_{\rm j}^ \pm$$ are entering a depth level. From expressions (B1a) and (B1b) it follows that for the elastic situation $$\delta \vec{{\rm S}}_{\rm j}^ -$$ and $$\delta \vec{{\rm S}}_{\rm j}^ +$$ are different. This difference can be written as:   $$\delta \vec{{\rm S}}_{\rm j}^ - \left( {{\rm z}_{\rm n}^ - ;{\rm z}_{\rm 0}^ - } \right) - \delta \vec{{\rm S}}_{\rm j}^ + \left( {{\rm z}_{\rm n}^ + ;{\rm z}_{\rm 0}^ + } \right) = \left[ {{{{\bf R}}^ \cup }\left( {{\rm z}_{\rm n}^ - ,{\rm z}_{\rm n}^ - } \right) - \delta {{{\bf T}}^ + }\left( {{\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ - } \right)} \right]{\vec{\rm P}_{\rm j}}^ + \left( {{\rm z}_{\rm n}^ - ;{\rm z}_0^{}} \right)\ - \left[ {{{{\bf R}}^ \cap }\left( {{\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ + } \right) - \delta {{{\bf T}}^ - }\left( {{\rm z}_{\rm n}^ - ,{\rm z}_{\rm n}^ + } \right)} \right]{\vec{\rm P}_{\rm j}}^ - \left( {{\rm z}_{\rm n}^ + ;{\rm z}_0^{}} \right).$$ (B2) Figure B1. View largeDownload slide The secondary sources $$(\delta \vec{{\rm S}}_{\rm j}^ \pm )$$ at a subsurface gridpoint represent the difference between the leaving $$(\vec{\rm Q}_{\rm j}^ \pm )$$ and the incident $$(\vec{\rm P}_{\rm j}^ \pm )$$ wavefields at that gridpoint. Hence, the secondary sources quantify the complex elastic wavefield interactions at each gridpoint. Figure B1. View largeDownload slide The secondary sources $$(\delta \vec{{\rm S}}_{\rm j}^ \pm )$$ at a subsurface gridpoint represent the difference between the leaving $$(\vec{\rm Q}_{\rm j}^ \pm )$$ and the incident $$(\vec{\rm P}_{\rm j}^ \pm )$$ wavefields at that gridpoint. Hence, the secondary sources quantify the complex elastic wavefield interactions at each gridpoint. It is easy to see from expression (B2) that $$\delta \vec{{\rm S}}_{\rm j}^ - = \delta \vec{{\rm S}}_{\rm j}^ +$$ if we assume δT + = R∪ and δT − = R∩ (acoustic situation). The acoustic assumption means a significant simplification of any imaging algorithm. The author thinks that this simplification is generally not justified in consolidated geological media and, therefore, should not be used in seismic imaging. Fig. B2 illustrates that if the shear contrast cannot be ignored, the acoustic simplification (δT = R) does not hold anymore. However, Fig. B2 also shows that a better approximation is obtained if we make the weak-elastic assumption: R∩ = −R∪ and δT − = −δT + or (R∪ − δT + ) = −(R∩ − δT − ). This assumption leads to:   $$\delta \vec{{\rm S}}_{\rm j}^ - \left( {{\rm z}_{\rm n}^ - ;{\rm z}_{\rm 0}^ - } \right) - \delta \vec{{\rm S}}_{\rm j}^ + \left( {{\rm z}_{\rm n}^ + ;{\rm z}_{\rm 0}^ + } \right) \approx \left[ {{{{\bf R}}^ \cup }\left( {{\rm z}_{\rm n}^ - ,{\rm z}_{\rm n}^ - } \right) - \delta {{{\bf T}}^ + }\left( {{\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ - } \right)} \right]\left[ {{{\vec{\rm P}}_{\rm j}}^ + \left( {{\rm z}_{\rm n}^ - ;{\rm z}_0^{}} \right)\ + {{\vec{\rm P}}_{\rm j}}^ - \left( {{\rm z}_{\rm n}^ + ;{\rm z}_0^{}} \right)} \right]$$ (B3a)or leaving the depth level out of the notation:   $$\delta \vec{{\rm S}}_{\rm j}^ - - \delta \vec{{\rm S}}_{\rm j}^ + \approx {{{\bf R}}^ - }({\vec{\rm P}_{\rm j}}^ + + {\vec{\rm P}_{\rm j}}^ - )\quad {\rm at} \quad {{\rm z}_{\rm n}},$$ (B3b)where R − = (R∪ − δT + ) = −(R∩ − δT − ) equals the conversion operator and where $$\delta \vec{{\rm S}}_{\rm j}^ - - \delta \vec{{\rm S}}_{\rm j}^ + = (\vec{\rm Q}_{\rm j}^ - + {\vec{\rm P}_{\rm j}}^ + ) - (\vec{\rm Q}_{\rm j}^ + + {\vec{\rm P}_{\rm j}}^ - )$$, being the difference between the total field above and below zn. Figure B2. View largeDownload slide Scattering is angle-dependent (upper part shows the forward and backward scattering operators; lower part shows the hybrid scattering operators). The angle-dependent wave conversion operator R − = (R∪ − δT + ) = −(R∩ − δT − ) is a measure for the wave conversion in a subsurface gridpoint. In subsurface areas with a significant R − , secondary S-sources are a logic extension to the secondary P-sources to generate the S-waves during the wavefield conversion process. Figure B2. View largeDownload slide Scattering is angle-dependent (upper part shows the forward and backward scattering operators; lower part shows the hybrid scattering operators). The angle-dependent wave conversion operator R − = (R∪ − δT + ) = −(R∩ − δT − ) is a measure for the wave conversion in a subsurface gridpoint. In subsurface areas with a significant R − , secondary S-sources are a logic extension to the secondary P-sources to generate the S-waves during the wavefield conversion process. The larger the difference between $$\delta \vec{{\rm S}}_{\rm j}^ -$$ and $$\delta \vec{{\rm S}}_{\rm j}^ +$$, the larger R − and, therefore, the larger the wave conversion. Note that the weak-elastic assumption also means:   $$\delta \vec{{\rm S}}_{\rm j}^ - + \delta \vec{{\rm S}}_{\rm j}^ + \approx {{{\bf R}}^ + }({\vec{\rm P}_{\rm j}}^ + - {\vec{\rm P}_{\rm j}}^ - )\quad {\rm at} \quad {{\rm z}_{\rm n}},$$ (B3c)where R + = (R∪ + δT + ) = −(R∩ + δT − ) represents the combined backward and forward scattering operator. Note that $$\delta \vec{{\rm S}}_{\rm j}^ - + \delta \vec{{\rm S}}_{\rm j}^ + = (\vec{\rm Q}_{\rm j}^ - + \vec{\rm Q}_{\rm j}^ + ) - ({\vec{\rm P}_{\rm j}}^ - + {\vec{\rm P}_{\rm j}}^ + ),$$ being the difference between the total leaving and the total incoming wavefields. Computation of R + and R − is an important option in our closed-loop, full-wavefield imaging algorithm. It leads to two new images that reveal the subsurface areas in which large wave conversion occurs. In these areas R and δT need be separately imaged. It is part of our current research on extending elastic P-wave imaging to S-wave imaging. APPENDIX C: MINIMIZATION ALGORITHM Let us consider a set of physical experiments, quantified by the following vector-matrix equations:   $$\begin{array}{*{20}{c}} {\vec{Y}_j^{} = {{\bf G}}\vec{X}_j^{}}&\quad{{\rm{for}}}&\quad {j{\rm{ }} = {\rm{ }}1,{\rm{ }}2, \ldots .,} \end{array}$$ (C1)where the elements of the system matrix G, indicated by Gkℓ, are the unknowns. These unknown elements can be estimated by a weighted least-squares minimization of the difference vector:   $$\Delta \vec{E}_j^\dagger {\Lambda _{jj}}\Delta {\vec{E}}\!\!\!\!\!\underline{\,\,\,}_j^{}\ = \ \vec{I}_j^\dagger \left[ {{{\bf Y}} - {{\bf \tilde{G}}}{{\bf \tilde{X}}}} \right]\ \Lambda {\left[ {{{\bf Y}} - {{\bf \tilde{G}}}{{\bf \tilde{X}}}} \right]^H}\vec{I}_j^{}\quad{\rm is\, minimum\,for\,all}\,j,$$ (C2)superscript † meaning a row vector, Λ being a user-specified weight matrix and H denoting Hermitian. In addition, $${{\bf \tilde{G}}} = {{\bf G}}{{{\bf L}}^{ - 1}}$$ and $${{\bf \tilde{X}}} = {{\bf LX}}$$, where operators L and L − 1 transform the internal coordinates of G and X to a suitable domain. In seismic imaging applications this transformation is often a τ − p transform, making the connection of G and X a connection per single plane wave component that occurs at each subsurface gridpoint. By setting the derivative of this weighted difference to $${\tilde{G}_{k\ell }}$$ and $${\tilde{{G} }\!\!\!\!\underline{\,\,\,}_{k\ell }}$$ equals zero for all k  and  ℓ, where $${{\bf \tilde{{G}}\!\!\!\!\underline{\,\,\,}}} = {{{\bf \tilde{G}}}^H}$$, then the solution is given by the normal equations (formulated in terms of matrices):   $${{\bf \tilde{G}}}\left[ {{{\bf \tilde{X}}}\Lambda {{{{\bf \tilde{X}}}}^H}} \right] = {{\bf Y}}\Lambda {{{\bf \tilde{X}}}^H},$$ (C3)where matrix Y contains all column vectors $${\vec{Y}_j}$$ and matrix $${{\bf \tilde{X}}}$$ contains all column vectors $${{\bf L}}{\vec{X}_j}$$. In its simplest version Λ is a diagonal matrix that contains the weights Λjj. To avoid inversion of correlation matrix $$[ {{{\bf \tilde{X}}}\Lambda {{{{\bf \tilde{X}}}}^H}} ]$$, we will use an iterative process. If we choose the elements of diagonal matrix Λ such that matrix $$[ {{{\bf \tilde{X}}}\Lambda {{{{\bf \tilde{X}}}}^H}} ]$$ has unit diagonal elements, then the first iteration is given by:   $${{{\bf \tilde{G}}}^{(1)}} = {{\bf Y}}\Lambda {{{\bf \tilde{X}}}^H},$$ (C4a)meaning that the first estimate of each element of $${{\bf \tilde{G}}}$$ is given by:   $${{{\bf \tilde{G}}}_{k\ell }}^{(1)} = \vec{I}_k^\dagger \left[ {{{\bf Y}}\Lambda {{{{\bf \tilde{X}}}}^H}} \right]\vec{I}_\ell ^{}.$$ (C4b) Of course, if a prior estimate of $${{\bf \tilde{G}}}$$ is already available, this prior estimate replaces expression (C4b). Using this first-order estimate, the error vector is given by   $$\Delta \vec{E}_j^{(1)} = \vec{Y}_j^{} - {{{\bf \tilde{G}}}^{(1)}}{{\bf L}}\vec{X}_j^{}.$$ (C4c) The first-order error vector is zero if the system would be orthogonal, meaning that $${{\bf \tilde{X}}}{{{\bf \tilde{X}}}^H}$$ equals a diagonal matrix. In practice such systems are rare and, therefore, let us look at the error vector of the second iteration:   \begin{eqnarray} \Delta \vec{E}_j^{(2)} &=& \vec{Y}_j^{} - {{{{\bf \tilde{G}}}}^{(2)}}{{\bf L}}\vec{X}_j^{}\nonumber\\ &=& \left[ {\vec{Y}_j^{} - {{{{\bf \tilde{G}}}}^{(1)}}{{\bf L}}\vec{X}_j^{}} \right] - \left[ {{{{{\bf \tilde{G}}}}^{(2)}} - {{{{\bf \tilde{G}}}}^{(1)}}} \right]{{\bf L}}\vec{X}_j^{}\nonumber\\ &=& \Delta \vec{E}_j^{(1)} - \left[ {{{{{\bf \tilde{G}}}}^{(2)}} - {{{{\bf \tilde{G}}}}^{(1)}}} \right]{{\bf L}}\vec{X}_j^{}. \end{eqnarray} (C5) If we substitute (C5) in eq. (C2), then we obtain the second-order minimization equation:   $$\Delta \vec{E}_j^{(2)}{\Lambda _{jj}}{\left[ {\Delta \vec{E}_j^{(2)}} \right]^H}\ = \ \left[ {\Delta \vec{E}_j^{(1)} - \Delta {{{{\bf \tilde{G}}}}^{(2)}}{{\bf L}}\vec{X}_j^{}} \right]\ {\Lambda _{jj}}{\left[ {\Delta \vec{E}_j^{(1)} - \Delta {{{{\bf \tilde{G}}}}^{(2)}}{{\bf L}}\vec{X}_j^{}} \right]^H}{\rm is\,minimum\,for\,all}\,j,$$ (C6)where $$\Delta {{{\bf \tilde{G}}}^{(2)}} = {{{\bf \tilde{G}}}^{(2)}} - {{{\bf \tilde{G}}}^{(1)}}$$. Minimization of (C6) leads to the second-order solution:   $$\Delta {{{\bf \tilde{G}}}^{(2)}} = \Delta {{{\bf E}}^{(1)}}\Lambda {{{\bf \tilde{X}}}^H},$$ (C7a)or   $${{{\bf \tilde{G}}}^{(2)}} = {{{\bf \tilde{G}}}^{(1)}} + \Delta {{{\bf E}}^{(1)}}\Lambda {{{\bf \tilde{X}}}^H}.$$ (C7b) Second-order solution (C7b) can be generalized to:   $${{{\bf \tilde{G}}}^{(i)}} = {{{\bf \tilde{G}}}^{({i} - 1)}} + \Delta {{{\bf E}}^{({i} - 1)}}\Lambda {{{\bf \tilde{X}}}^H}$$ (C8)with  \begin{equation*}{{{\bf \tilde{G}}}^{(0)}} = {{\bf 0}},\quad\Delta {{{\bf E}}^{(0)}} = {{\bf Y}}\quad{\rm{and}}\quad\Delta {{{\bf E}}^{(1)}} = {{\bf Y}} - {{{\bf \tilde{G}}}^{(1)}}{{\bf \tilde{X}}}.\end{equation*} Solution (C8) is used in all iterative estimation algorithms for the wavefield operators. Note that if we deal with the related equations:   $$\begin{array}{*{20}{c}} {\vec{Y}_k^{} = {{\bf X}}\vec{G}_k^{}}&{{\rm{for}}}&\quad{k = 1,{\rm{ }}2, \ldots .\,,} \end{array}$$ (C9a)the iterative solution is given by:   $${{{\bf \tilde{G}}}^{(i)}} = {{{\bf \tilde{G}}}^{({i} - 1)}} + {{{\bf \tilde{X}}}^H}\Lambda \Delta {{{\bf E}}^{({i} - 1)}},$$ (C9b)where $$\Delta {{\bf E}} = {{\bf Y}} - {{\bf \tilde{X}\tilde{G}}}.$$ Note also that in the special situation of matrix inversion,   $${{\bf I}} = {{\bf FX}}\quad{\rm or}\quad{{\bf I}} = {{\bf XF}},$$ (C10a)the iterative solution is given by:   $${{{\bf \tilde{F}}}^{(i)}} = {{{\bf \tilde{F}}}^{({i} - 1)}} + \Delta {{{\bf E}}^{({i} - 1)}}\Lambda {{{\bf \tilde{X}}}^H}\quad{\rm or}\quad{{{\bf \tilde{F}}}^{(i)}} = {{{\bf \tilde{F}}}^{({i} - 1)}} + {{{\bf \tilde{X}}}^H}\Lambda \Delta {{{\bf E}}^{({i} - 1)}},$$ (C10b)where   \begin{equation*}\Delta {{\bf E}} = {{\bf I}} - {{\bf \tilde{F}\tilde{X}}}\quad{\rm or}\quad\Delta {{\bf E}} = {{\bf I}} - {{\bf \tilde{X}\tilde{F}}}\end{equation*} APPENDIX D: SUMMARY OF CONCEPTS AND NOTATION If the time domain is transformed to the frequency domain, a seismic data volume consists of orthogonal slices. Each frequency slice is represented by the data matrix. A measured physical wavefield is represented by one column of the data matrix, ordered according to the x–y observation coordinates (detector positions). A row represents a common detector gather, now ordered by the x-y emission coordinates (source positions). Each column and row may have its own source and detector depth level. Primary sources are man-made and often positioned (close) at the surface. Wavefields consist of downgoing and upgoing constituents, indicated by the uppercases + and − respectively. The definition down and up depend on the orientation of the coordinate system used. Downgoing wavefields illuminate the upper side of reflectors, upgoing wavefields illuminate the lower-side of reflectors. This article focuses on P-waves, but wave conversion during the PP-scattering process is taken into account. We used the following notation: $${{{\bf P}}^ + }({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}})$$: downgoing wavefields that enter the upper part of depth level $${\rm z}_{\rm m}^{}$$ (realizing illumination from above). $${{{\bf P}}^ - }({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})$$: upgoing wavefields that enter the lower part of depth level $${\rm z}_{\rm m}^{}$$ (realizing illumination from below). $${{{\bf Q}}^ + }({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}})$$: downgoing wavefields that leave the lower part of depth level $${\rm z}_{\rm m}^{}$$; note that $${{{\bf P}}^ + }({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}) = {{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_{{\rm m} - 1}^ + ){{{\bf Q}}^ + }({\rm z}_{{\rm m} - 1}^ + ;{{\rm z}_{\rm s}})$$, where each column of W equals a local downward wavefield propagation operator (generally being anisotropic). $${{{\bf Q}}^ - }({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}})$$: upgoing wavefields that leave the upper part of depth level $${\rm z}_{\rm m}^{}$$; note that $${{{\bf P}}^ - }({\rm z}_{\rm m}^ + ;{{\rm z}_{\rm s}}) = {{\bf W}}({\rm z}_{\rm m}^ + ,{\rm z}_{\rm m + 1}^ - ){{{\bf Q}}^ - }({\rm z}_{\rm m + 1}^ - ;{{\rm z}_{\rm s}})$$, where each column of W equals the local upward wavefield propagation operator (generally being anisotropic). Decomposition of the downgoing wavefields for triple illumination from above:   \begin{equation*} {{{\bf P}}^ + }({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}) = {{\bf W}}({\rm z}_{\rm m}^ - ,{{\rm z}_{\rm s}}){{{\bf S}}^ + }({{\rm z}_{\rm s}}) + {{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_0^ + )[{{{\bf R}}^ \cap }({\rm z}_0^ + ,{\rm z}_0^ + ){{{\bf P}}^ - }({\rm z}_0^ + ;{{\rm z}_{\rm s}})] + {{\bf W}}({\rm z}_{\rm m}^ - ,{\rm z}_{\rm n}^ + )\sum\limits_{\rm n = 1}^{{\rm m} - 1} {{{{\bf R}}^ \cap }({\rm z}_{\rm n}^ + ,{\rm z}_{\rm n}^ + ){{{\bf P}}^ - }({\rm z}_{\rm n}^ + ;{{\rm z}_{\rm s}})}, \end{equation*} where the first term equals the downgoing wavefield generated by the primary sources, the second term equals the downgoing wavefield generated by the secondary sources at the surface (surface backscattering) and the third term equals the down going wavefield generated by the secondary sources in the subsurface (internal backscattering). Decomposition of the upgoing wavefields for triple illumination from below:   \begin{equation*} \quad{{{\bf P}}^ - }({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}) = {{\bf P}}_{\rm pr}^ - ({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}) + {{\bf M}}_{\rm sfc}^ - ({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}) + {{\bf M}}_{{\mathop{\rm int}} }^ - ({\rm z}_{\rm m}^ - ;{{\rm z}_{\rm s}}), \end{equation*} where the first term equals the first-order response from the man-made sources (‘primaries’), $${{\bf P}}_{\rm pr}^ - = {{\bf X}}_{\rm pr}^ \cup {{{\bf W}}^ + }{{{\bf S}}^ + }$$ with $${{\bf X}}_{\rm pr}^ \cup = \sum\nolimits_{\rm n > {\rm m}} {{{{\bf W}}^ - }{{\bf R}}_{\rm a}^ \cup {{{\bf W}}^ + }}$$, the second term equals the response from the surface backscattering (‘surface multiples’) $${{\bf M}}_{\rm sfc}^ - = {{\bf X}}_{\rm pr}^ \cup {{{\bf W}}^ + }\delta {{\underline{\bf P}}^ + }$$ with $$\delta {{\underline{\bf P}}^ + } = {{\bf R}}_{}^ \cap {{\underline{\bf P}}^ - }$$ at $${\rm z}_0^ +$$ and $${{\bf X}}_{\rm pr}^ \cup = \sum\nolimits_{\rm n > {\rm m}} {{{{\bf W}}^ - }{{\bf R}}_{\rm b}^ \cup {{{\bf W}}^ + }}$$, and the third term equals the response from the internal backscattering (‘internal multiples’), $${{\bf M}}_{{\mathop{\rm int}} }^ - = {{\bf X}}_{\rm pr}^ \cup \sum\nolimits_{\ell < {\rm m}} {{{{\bf W}}^ + }\delta {{{\bf P}}^ + }}$$ with $$\delta {{{\bf P}}^ + } = {{\bf R}}_{}^ \cap {{{\bf P}}^ - }$$ at $${\rm z}_\ell ^ +$$ with ℓ > 0 and $${{\bf X}}_{\rm pr}^ \cup = \sum\nolimits_{\rm n > {\rm m}} {{{{\bf W}}^ - }{{\bf R}}_{\rm c}^ \cup {{{\bf W}}^ + }}$$. One column of back scattering operators R∪ (up) and R∩ (down) defines the angle-dependent scattering properties in one gridpoint. These properties are determined by the faster changes of the density (ρ) and the velocities (cp, cs) in that gridpoint. Operators R∪ and R∩ are output of the closed-loop, full-wavefield imaging algorithm. By definition they include wave conversion. One column of propagation operators W − (up) and W + (down) defines the anisotropic propagation properties between gridpoints. These properties are determined by the slower changes of the velocity (cp) between these gridpoints. Operators W − and W + are output of the closed-loop, full-wavefield imaging algorithm. In this article they propagate P-waves only. Wavefield extrapolation may be carried out on the down-and upgoing constituents (P ± , Q ± ) or on the combined versions (Q + + P − , Q + − P − ) and (Q − + P + , Q − − P + ). Estimation of operators R∪ and R∩ always occurs on the separate constituents (P ± , Q ± ), preferably in the linear Radon (τ-p) domain. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Mar 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations