Simultaneous Imaging of Primaries and Multiples Based on Wavelet Estimation and Stereographic Imaging Condition

Multiple reflections contain abundant structural information about the subsurface because of their smaller reflection angles and wider coverage. Multiple imaging has increasingly attracted attention. However, the imaging of multiples first needs to separate the multiples, which is a time-consuming and high-cost job. Therefore, the simultaneous imaging of primaries and multiples, which can make good use of multiple reflections with no need for multiple separation, has increasingly attracted attention. However, there are still some challenges to the conventional method, including wavelet selection and image artifacts suppression. In this paper, an improved method of the simultaneous imaging of primaries and multiples is proposed by introducing two creative strategies. First, considering that the wavelet, which is usually estimated by signal analysis under certain assumptions, is crucial to the simultaneous imaging, we propose to estimate the wavelet by iterative SRME (surface-related multiple elimination) using the original data containing multiples based on wave equation theory instead of signal analysis. Second, the stereographic imaging condition is introduced to suppress the crosstalk artifacts in the image. According to the numerical examples and field data test, the feasibility and effectiveness of our approach is verified.


I. INTRODUCTION
When a seismic wave propagates in the subsurface and encounters formation interfaces, reflection or transmission occurs. Since the structures are complex with large impedance contrasts between some layers, multiple reflections, which are reflected more than once, will generate, which make the seismic data complex and difficult to deal with.
An accurate image of the subsurface is of vital importance to oil and gas exploration and natural earthquake research. Traditionally, only primary reflections (primaries) are considered as signals in the typical imaging strategy, The associate editor coordinating the review of this manuscript and approving it for publication was Madhu S. Nair .
whereas other events, including multiple reflections (multiples), are regarded as noise that needs to be suppressed as much as possible before imaging. Generally, multiples are suppressed based on the time difference between the primary and multiple reflections, the predictability of multiple reflections, or the periodicity of multiple reflections [1]- [6].
However, multiples are also real responses reflected from formation interfaces under the surface. Multiples usually propagate with longer wave paths, smaller reflection angles and wider illumination areas than primaries [7], [8]; therefore, they theoretically contain more information about the subsurface structures. If appropriately imaged, multiples can provide wider illumination, higher fold, and better imaging results for the subsurface than primaries [9]- [13]. As a result, multiple imaging instead of multiple suppression has drawn VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ increasing attention in recent years. Owing to the strong energy of surface-related multiples, most studies focus on the imaging of surface-related multiples. As in these studies, the word ''multiple'' refers to surface-related multiple in our research, unless otherwise noted. Currently, most commonly used multiple imaging techniques are performed by altering the boundary conditions, replacing the source wavelet with the recorded data and using the multiples as the input receiver data instead of primary reflections first. Then, with Kirchhoff migration [14], one-way wave equation migration [9], [15] or reverse time migration (RTM) [11], [13], [16], multiples can be imaged to the right position. However, the accuracy of Kirchhoff migration may not be guaranteed in some complex cases because of its high-frequency approximation, and the computational cost of RTM is considerable. Considering both accuracy and efficiency, one-way wave equation migration is chosen to perform the imaging in our approach.
The implementation of the multiple imaging techniques discussed above needs to separate multiples from the recorded seismic data first, which is a time-consuming and high-cost job. In addition, these techniques cannot effectively deal with the combination of the separated primary image and multiple image. To avoid the complex process of multiple separation, the simultaneous imaging of primary and multiple reflections has been developed by replacing the source wavelet through adding the recorded data on the source wavelet and using the recorded data as receiver data [17], [18].
However, an arbitrary wavelet usually cannot appropriately and measurably tackle the combination of primary and multiple images. If the energy of the arbitrary wavelet is too small, the image result of primary imaging will be affected, whereas when the energy of the selected wavelet is too large, the imaging information of multiple waves will be covered by the primary image. Moreover, an improper wavelet may cause low resolution and inaccurate image location. Therefore, an accurate wavelet is extremely vital to the image quality of the simultaneous imaging of primaries and multiples. In addition, the cross-correlation of unrelated wavefields by using the conventional imaging condition results in serious crosstalk artifacts, which is quite a tough problem that is similar to the conventional multiple imaging. Modifying the imaging condition is proposed to eliminate the crosstalk artifacts in multiple imaging [19], [20]. However, few studies have been done on crosstalk suppression for simultaneous imaging.
Since an incorrect wavelet may lead to a terrible image, an accurate estimation of the real wavelet should be done for the simultaneous imaging of primaries and multiples. The current methods for wavelet estimation are mostly realized based on signal theory of the convolution model under certain assumptions [21]. In addition, the methods are all carried out without the effect of multiples. How to deal with the wavelet estimation for the data with multiples is worthy of pondering. Inspired by multiple estimation, a method based on wave equation theory may be a better choice compared to a method based on signal analysis theory. In the early stage, wavelet estimation is extremely important to the surface-related multiple elimination (SRME) method [5], [22] and the extended method [23], which are proposed based on wave equation theory. With the development of SRME by iterative inversion [6], [24], the wavelet is no longer necessary. However, we do reverse thinking in this paper: the related theories can provide us a way to estimate the wavelet based on wave equation theory.
For the suppression of crosstalk artifacts in the image, Jiang et al. [25] suggest that one solution is to separate primaries and different-order multiples before imaging [26], [27], and the other solution is least-square migration [13], [28], [29], whereas they are all at the cost of additional large computational costs. Muijs et al. [19], Ravasi [30] propose to suppress the crosstalk artifacts by introducing the 2D deconvolution imaging condition, but residual artifacts are caused by the interference of upgoing and downgoing waves not associated with the same reflector. Nevertheless, the method should be of value in the development of a new way to suppress crosstalk by modifying the imaging condition [20], [31]. Li et al. propose to apply the stereographic imaging condition to suppress the crosstalk artifacts for multiple imaging [20]. Examples verify that the crosstalk artifacts can be well suppressed. Thus, we extend the method used in multiple imaging to deal with the crosstalk artifacts in the simultaneous imaging of primaries and multiples.
In this paper, we propose an improved method of the simultaneous imaging of primaries and multiples mainly by introducing two creative strategies to solve the problems discussed above. First, to better describe and balance the contribution of the primary image and multiple image, we suggest to estimate the wavelet through iterative SRME based on wave equation theory to effectively guarantee the accuracy of wavelet estimation for data with multiples, and then use the estimated wavelet together with the recorded data as the source wavefield, which is quite different from selecting an arbitrary wavelet in conventional simultaneous imaging. Second, the stereographic imaging condition [20], [32] is introduced to suppress the crosstalk artifacts for simultaneous imaging. Numerical and field data examples are tested to demonstrate the feasibility and validity of our method.

II. SIMULTANEOUS IMAGING OF PRIMARIES AND MULTIPLES A. CONVENTIONAL SIMULTANEOUS IMAGING BASED ON ONE-WAY WAVE EQUATION MIGRATION
In conventional prestack depth imaging, the imaging techniques based on one-way wave equation theory are widely used to realize the structural imaging because the one-way wave equation takes both the accuracy and efficiency into account. These imaging techniques can image structures with relatively high accuracy and less computational cost. Therefore, we use the one-way wave equation migration method to realize simultaneous imaging of primaries and multiples in this paper.
Assume that sources on the surface transmit waves into the earth. After removing the source and receiver effects, the primaries can be approximately expressed in the frequency domain by matrices as [11], [22], [26], [33] P 0 (z 0 ) = X 0 (z 0 , z 0 )S(z 0 ) where P 0 denotes the matrix of the primaries, X 0 denotes the matrix of the subsurface response irrespective of multiples, S denotes the matrix of the source wavelet signature, and z 0 is the surface notation. Equation (1) reveals the basic idea of primary imaging. That is, structural images can be obtained by extrapolating the source wavefield and the receiver wavefield into the earth interior on the basis of the one-way wave equation migration first, then the imaging condition [34] can be applied to obtain the related reflectivity information from the forward-extrapolated source wavefield and backwardextrapolated receiver wavefield at each depth-step [20]. This process is shown in Fig. 1a.
If we use P to denote the recorded data containing both the primary and the multiple reflections with different orders, the relationship between them can be expressed as where M is the multiples, and M 1 , M 2 and M 3 denote the first-, second-, third-order multiples, respectively. According to the generation mechanism of multiple reflections, the multiples can be further expressed as Here, the surface reflectivity is assumed to be -1. Equation (3) describes the basic idea of multiple imaging. Fig. 1b shows this concept by taking the imaging of the first-order multiples as an example, and the same goes for the higher-order multiples. The recorded data can be treated as a virtual source. Then, the virtual source and the multiples separated from the recorded data are forward-extrapolated and backwardextrapolated, respectively, at each depth-step, and the multiple image at location R' in Fig. 1b can be obtained by applying the conventional imaging condition.
Substituting equations (1) and (3), we can obtain the following expression P 0 (z 0 ) + M(z 0 ) = X 0 (z 0 , z 0 )(S(z 0 ) − P(z 0 )) = P(z 0 ) (4) Similar to equations (1) and (3), we can achieve the basic idea of the simultaneous imaging of primaries and multiples as shown in Fig. 1c. Simultaneous imaging can be performed by replacing the wavelet function with the estimated source wavelet together with the recorded data as the source wavefield and using the recorded data as the receiver wavefield first. Then, extrapolating the source and receiver wavefields at each depth-step using one-way wave equation migration, the simultaneous image can finally be given by the conventional imaging condition as where I is the simultaneous image, the subscripts F and B represent the forward and backward extrapolation of the source and receiver data, respectively. The operator '' * '' stands for the crosscorrelation. As a result of the recorded data being known, the image quality depends only on the unknown wavelet. Therefore, an appropriate wavelet is extremely important to simultaneous imaging. The image can be further written by expanding equation (5) as Here, the first part contains the true image of the primaries using the estimated wavelet and the true image of different order multiples with matched wavefields, whereas the second part leads to undesired crosstalk artifacts in the image. If a proper imaging condition can be chosen to suppress the crosstalk artifacts, an image with higher quality can be obtained. From the theoretical derivation and analysis above, we faced two important challenges during the conventional simultaneous imaging of primaries and multiples: one is the inaccurate wavelet and the other is the artifacts, which should VOLUME 8, 2020 be avoided by a proper imaging condition. Therefore, to accurately implement the simultaneous imaging of primaries and multiples, an accurate wavelet should be extracted from the original data. In addition, a proper imaging condition should be selected.

B. WAVELET ESTIMATION USING MULTIPLES BASED ON WAVE EQUATION THEORY
Most of the current wavelet estimation methods are based on signal analysis theory under certain assumptions, and the accuracy and reliability of the estimated wavelet need to be further improved. Moreover, these methods are all adaptive to the multiple-free seismic data, but failure to deal with the seismic data with multiples. It means that multiple suppression should be done first, which is time-consuming and easily introduces errors during the suppression of multiples. Therefore, it is important to develop an accurate wavelet estimation method that can work for seismic data with multiples. To solve this problem, a wavelet estimation method extended from SRME is presented in this paper based on propagation theory of the seismic wave and the relationship of primary and multiple reflections. Compared to the traditional methods, the proposed method not only can estimate the wavelet with higher accuracy based on wave equation theory but also can be applied to data with multiples.
Similar to equation (1), after removing the source and receiver effects and considering that the sources and receivers do not show any variations during the seismic survey, the recorded data including multiples can be expressed by matrices as where X denotes the matrix of the subsurface response including multiples. Fig. 2 shows the forward model of seismic data without and with multiples, it shows the comparison between equations (1) and (7), where R − (z 0 ) defines the reflection operator for upgoing incident wavefields at the surface level. According to Fig. 2, the relationship between X 0 and X can be given as The resulting relationship between P 0 and P can be written as Further derivation can lead to the following expression where The matrix A(z 0 ) is defined as the surface operator. If we may assume the reflection operator at the surface R − (z 0 ) = −I, then the estimation of A means an estimation of [S(z 0 )] −1 .
In addition, since we assume that the directional source effects have been corrected in advance and that the sources do not show any variations, then A can be represented by a scaled unity matrix. The scaling factor defines one Fourier component of the inverse of the source wavelet as Let A(ω) = S −1 (ω), which represents the inverse source wavelet in the frequency domain. Based on this simplification and equation (10), the iterative procedure of SRME can be expressed as where n denotes the times of iteration. The surface function A(ω) can be optimized for each iteration by minimizing where E{} means averaging over all frequencies within the seismic bandwidth. The estimation of A(ω) can be realized by a number of linear optimization processes [6], [22], [24]. Through a stabilized inversion of the surface function, the source signature S(ω) can be estimated. Then, the estimated wavelet in the time domain can be obtained through the inverse Fourier transform where s(t) denotes the estimated wavelet. This iterative procedure can not only improve the accuracy of multiple estimation but also provide a way to obtain an estimated wavelet through the inversion approach. The wavelet is estimated based on wave equation theory, which is more advanced than signal theory. Therefore, for the seismic data with strong multiples, the iterative SRME provides an optional and accurate method for wavelet estimation. With the accurately estimated wavelet, the simultaneous imaging of primaries and multiples can be done by applying a proper imaging condition.

C. STEREOGRAPHIC IMAGING CONDITION FOR SIMULTANEOUS IMAGING
The imaging condition is of great importance to the image quality. The crosscorrelation imaging condition is widely used in one-way wave equation migration since it is robust and fast, allowing structural images with high accuracy to be obtained under complex cases. It can be carried out by crosscorrelating the extrapolated source and receiver wavefields with zero lag in space and time. If extrapolated accurately, we can obtain images of high accuracy at zero crosscorrelation lag. However, this conclusion is not applicable to the simultaneous imaging of primaries and multiples. Equation (6) shows that the crosscorrelation of those unrelated wavefields results in tiresome crosstalk artifacts with the time constraint only, causing a poor result of simultaneous imaging. Suppressing the crosstalk artifacts to obtain an accurate and reliable image of the simultaneous imaging presents a considerable challenge.
Considering that the unrelated wavefield events will surely have different local slopes at the same time because of their different propagating directions and the stereographic imaging condition matches the forward-extrapolated source and backward-extrapolated receiver wavefields using both propagating times and local slopes of the extrapolated wavefields [20], [32], the unrelated wavefields with different propagating directions at the same time can be distinguished by the extra local slope constraint, and the crosstalk artifacts in the simultaneous image will be well suppressed. Consequently, in this paper, the stereographic imaging condition is suggested to suppress the crosstalk energy in simultaneous imaging of primaries and multiples.
An extra parameter, local slope, is considered in the stereographic imaging condition. First, the extrapolated source and receiver wavefields are decomposed by local slant stacks at every position x =(x, z) and time t. The source and receiver wavefields can be expressed as [32] U S (x, t) = p W S (x, p, t)dp (16) where U S (x,t) denotes the forward-propagated source wavefield,Û R (x,t) denotes the back-propagated receiver wavefield, p is the local slope function of position x and time t, W S andŴ R are the decomposed forward-extrapolated source and backward-extrapolated receiver wavefields, respectively, which are obtained by the local slant stacks over a small number of neighboring x locations through a window (19) where w is a windowing function, i.e., in the simplest case a rectangle function, x 1 and x 2 represent the range of the window, x 0 fixes the center of the window, and τ is the time intercept at x = x 0 [35], [36]. Then, the stereographic imaging condition is implemented by crosscorrelating the decomposed source and receiver wavefields, W S andŴ R , and superimposing over the decomposed variables In equation (20), the decomposed forward-extrapolated source wavefields, W S , and the backward-extrapolated receiver wavefields,Ŵ R , propagating with different directions at the same time, can be distinguished by taking the local spatial coherence into account. In contrast, in the crosscorrelation imaging condition, (21) the decomposed wavefields with different local slopes are all summed together. Therefore, the wavefields propagating with different directions at the same time cannot be distinguished. Thus, crosstalk artifacts are introduced by the constraint of propagating times only. However, the stereographic imaging condition can provide a reasonable way to suppress the crosstalk artifacts with the extra local slope constraint for simultaneous imaging. Fig. 3 shows the simultaneous imaging procedure using the stereographic imaging condition. Assuming the bold line represents the extrapolated source wavefield and the thin line represents the extrapolated receiver wavefield, Fig. 3a shows the imaging process of unrelated wavefields. Taking the forward-propagated source wavelet and the backward-propagated first-order multiple as an example, the two events are unrelated with different propagating directions. As seen in the figure, the two events have different local slopes at the same time (p 1 = p 2 ) when they correlate in positions x 1 and x 2 . If the crosscorrelation imaging condition is applied, artifacts will be introduced in positions x 1 and x 2 with the constraint of the propagating time only. However, the stereographic imaging condition will avoid the crosstalk artifacts by the extra local slope constraint. The forward-extrapolated source wavefield and backward-extrapolated receiver wavefield propagating in different directions will not generate crosstalk any more in positions x 1 and x 2 (as shown in Fig. 3a). The case is FIGURE 3. Illustration of the simultaneous imaging procedure using stereographic imaging condition: (a) imaging process of unrelated wavefields and (b) imaging process of related wavefields. The red and green lines are the tangent lines of the extrapolated source and receiver wavefields, respectively, at a certain position, and represent the local slopes. VOLUME 8, 2020 the same for the other unrelated wavefields. Consequently, the crosstalk artifacts introduced in the simultaneous imaging of primaries and multiples can be eliminated with the application of the stereographic imaging condition. Fig. 3b shows the imaging process of the related wavefields. Taking the the forward-propagated source wavelet and the backwardpropagated primary as an example, the local slope is equivalent at the same time, and the real image will generate in positions x 1 and x 2 .
The local slant stacks are implemented at every position for both the source and receiver wavefields, which greatly raise the requirement of computational cost and memory. The simplification on the window is an effective way to decrease the cost [35]. In addition, the constraint of estimated reflector slopes is introduced to further improve the computational effectiveness in our method.
In equation (21), the 3D vector p is the local slope function of the positions (x, z) and time (t). The local slope can be uniquely determined by two components, the local slope component in the x−t domain, p t , and the local slope component in the x−z domain, p z , of the extrapolated wavefields. The local slope component in the x−z domain is determined by the local slope of the imaged reflectors. Therefore, the correspondence between the local slopes p of the decomposed extrapolated-source and receiver wavefields occurs only in planes dipping with the slope of the reflectors at every position x = (x, z) in space [32]. If there is no knowledge of the reflector slopes, we need to loop over a wide range of possible slopes in the vicinity of the expected reflector slopes, which will greatly increase the computation costs.
Therefore, we introduce plane-wave destructors and apply them in approximating the local slopes of the imaged reflectors. Then, the approximated reflector slopes can be used to constrain the crosscorrelation of the decomposed forward-extrapolated source wavefield and backward-extrapolated receiver wavefield by choosing the local slopes around the approximated slopes during the loop of the local slopes in the x−z domain [20]. The loop over large range of all slopes in the x−z domain can be avoided, and the efficiency can be greatly improved. In addition, relatively accurate local slopes in the x−z domain can effectively ensure the accuracy of the simultaneous imaging of primaries and multiples.
Once the reflector slope is estimated, the local slope component in the x − z domain can be determined, and the image will form in places where the decomposed forward-extrapolated source wavefield and backwardextrapolated receiver wavefield have the same local slope in the x − t domain [20], [32]. To obtain an image of high accuracy, small sampling intervals and a plentiful number of local slopes in the x − t domain are preferred to choose.
The simultaneous imaging process of primaries and multiples based on one-way wave equation migration, which is proposed in this paper, can be performed as shown in Fig. 4. The first important step is the wavelet estimation based on iterative SRME using multiples. After obtaining an estimated wavelet, the recorded data together with the wavelet and the recorded data can be regarded as the source and receiver wavefields, respectively, and extrapolated based on one-way wave equation theory. Then, we can perform wavefield decomposition of the extrapolated wavefields basing on local slant stacks and use plane-wave destructors to estimate the reflector slopes. Finally, stereographic imaging can be completed by crosscorrelating the wavefields in the decomposed domain followed by a summation over the decomposition variables with the constraint of the reflector slopes, and the simultaneous image of primaries and multiples can eventually be obtained.
In this section, we have discussed the conventional simultaneous imaging of primaries and multiples, pointed out the main challenges, and presented two creative strategies to address these problems: 1) a wavelet estimation method using multiples is presented based on wave equation theory to ensure the accuracy of the estimated wavelet for the data with multiples, and 2) the stereographic imaging condition is suggested to address the tiresome crosstalk artifacts to improve the image quality. Next, numerical and field data examples will be used to further validate the effectiveness of the proposed method.

III. EXPERIMENTS AND ANALYSIS
Two models and a field data test are used to demonstrate the validity of the proposed method. Note that without the requirement of multiple separation in our method, only the conventional pre-processing is needed to the original seismic data before the imaging, similar to the conventional primary imaging.

A. FLAT-LAYERED MODEL
A flat-layered model is first used to verify the effectiveness. The velocity model is shown in Fig. 5 with a model size 99052 VOLUME 8, 2020 A known modeling model should be the best way to test the effectiveness of the wavelet estimation since the used wavelet is known, and the consistency comparison between the estimated wavelet and the used wavelet can effectively judge the validity of the method. Hence, we use the modeling data to test the effectiveness of the wavelet estimation method presented in this paper. First, we apply the iterative SRME to obtain the surface operator, and then the estimated wavelet can be obtained, which is compared with the real wavelet, as shown in Fig. 6. The dotted line denotes the real wavelet, which is used to simulate the recordings; the solid line represents the estimated wavelet using multiples based on iterative SRME theory. The estimated wavelet is highly consistent with the real wavelet. Therefore, it can be proved that our method is valid in estimating the source wavelet. The presented method is carried out based on wave equation theory, which is more advanced than the signal analysis, avoiding the assumptions in methods based on convolution models and providing a method of wavelet estimation for data containing multiples.
After wavelet estimation, the wavelet and the recorded data can be used as the source wavefield to do the forward extrapolation, and the recorded data can be regarded as the receiver wavefield to do the backward extrapolation based on one-way wave equation migration theory. Because the model is simple, we only take one shot record, which is shot in the middle of the model, as an example to validate the advantages of our method. The source and receiver wavefields at the first time sample point are composed as shown in Fig. 7a and Fig. 7b. As shown in the figure, there are obvious multiple reflections in the recorded data that may contain plentiful information about the subsurface. Conventionally, an arbitrary wavelet is used to obtain the simultaneous image. However, the inappropriate waveform of the chosen wavelet may lead to an unreliable image. In addition, the improper energy of the chosen wavelet cannot effectively deal with the combination of the image contribution between the primary image and the multiple image. To emphasize the advantages of our method, we choose an arbitrary wavelet to do the simultaneous imaging, which is the same as the conventional method. We also use the estimated wavelet to perform the simultaneous imaging. The comparison of the images using both wavelets is shown in Fig. 8. Fig. 8a is the conventional primary image using an arbitrarily selected wavelet. Fig. 8b shows the conventional simultaneous image with the same wavelet. Fig. 8c and Fig. 8d are the simultaneous images using the conventional imaging condition and the stereographic imaging condition, respectively, with the accurate wavelet that is estimated from the shot data by our method. From the comparison of Fig. 8a with Fig. 8c, we can see that the energy of the selected wavelet is so strong that the image information from the multiple reflections is completely covered up in Fig. 8b. In contrast, if the selected wavelet is too weak, the accuracy of the primary image will be seriously impacted by the strong multiples. Moreover,  the image location of the first layer is compared as indicated by the white box in Fig. 8. Note that the real depth of the first layer is located in the middle of the white box. It can be seen from the figure that the erroneously chosen wavelet yields an incorrect image location, as shown in Fig. 8a and Fig. 8b, while the estimated wavelet by our method leads to an accurate location of the images in Fig. 8c and Fig. 8d. In addition, the comparison of Fig. 8c and Fig. 8b also demonstrates that not only can more information be gained with the image supplemented with multiples but also the resolution is greatly improved by using the accurate wavelet. Therefore, the comparison of the image results verifies the importance of wavelet extraction in the simultaneous imaging of primaries and multiples.
However, heavy crosstalk artifacts are introduced because of the crosscorrelation of the unrelated wavefields. To remove the crosstalk, we introduce the stereographic imaging condition to perform the simultaneous imaging by the method   suppressed. The feasibility of this method can be effectively verified.
To further emphasize the contribution of multiple reflections to the image, we compare the illumination of the shot record with different methods in Fig. 9. Fig. 9a is the illumination of the conventional primary imaging with an arbitrary wavelet, which corresponds to the image in Fig. 8a. Fig. 9b shows the illumination of the conventional simultaneous imaging with the same wavelet as in Fig. 9a, which   FIGURE 12. Velocity of the field data.
corresponds to the image in Fig. 8b. Differing from that in Fig. 8a, the image in Fig. 8b contains not only the primary image but also the multiple image. It can be seen that the illumination energy in Fig. 9b is enhanced and that the illumination range is enlarged compared to that in Fig. 9a, which shows that the application of multiple reflections can provide more information about the subsurface. However, the improvement is not significant enough because the wavelet is inappropriate. Fig. 9c shows the illumination of the simultaneous imaging with the accurate wavelet estimated by our method, which corresponds to the image in Fig. 8c and Fig. 8d. We can see that the illumination is greatly improved compared to that of Fig. 9a and Fig. 9b. Therefore, more reliable information should be obtained by applying the proposed simultaneous imaging method using both primary and multiple reflections compared to the that of conventional imaging method using primary reflections only.

B. COMPLEX MODEL
To further demonstrate the advantage of our method, especially on the acquisition of high-quality images and more structural information, the proposed method is validated with a more complex model including flat layers, anticline, overlaying structures and a fault. The velocity of the model is shown in Fig. 10. It is discretized into 401 × 251 grid points with a grid interval of 10 m in both the vertical and the horizontal directions. The source is located in the middle of the shot at a depth of 10 m. There are 100 shots with a 20 m shot interval and 101 geophones spaced of 20 m for each shot. Fig. 11 compares the images obtained using different imaging methods. Fig. 11a shows the conventional primary image using the conventional imaging condition. Since the input receiver-side seismic data contains both primaries and multiples, crosstalk caused by the multiples is introduced in the image (seen as the white arrows). Fig. 11b shows the simultaneous image with the estimated accurate wavelet using the conventional imaging condition. Similarly, the crosscorrelation of unrelated wavefields results in the heavy crosstalk artifacts, as shown by the white arrows. Comparing Fig. 11a with Fig. 11b (white line at x = 0.5 km), it can be seen that a wider image range can be gained in Fig. 11b due to VOLUME 8, 2020 the wider illumination and higher fold from multiple reflections. Fig. 11c is the simultaneous image using the method described in this paper. Most of the crosstalk artifacts are well suppressed (seen as the white arrows) with the extra constraint of local slopes in the stereographic imaging condition. The advantages of our method can be effectively validated.
However, there are still some faint residual artifacts, which are generated for some unrelated events that may have similar local slopes at a certain time, so the suppression of these artifacts is still a challenging problem for future study.

C. FIELD DATA
Next, the proposed method is demonstrated with the field data. The estimated velocity model, as shown in Fig. 12, is obtained by velocity analysis. It is composed of 316 (in X) by 600 (in depth) grid points with grid intervals of 25 m (in X) and 8 m (in depth). Fifty shots are extracted from the recorded field data to image the area. One single shot record contains 120 traces with a trace interval of 25 m, and the shot interval is 50 m. The total recording time is 4.8 s, and the sample rate is 8 ms. The source-receiver near offset is 200 m, and the far offset is 3175 m. Fig. 13 compares the images obtained with different imaging methods. First, the conventional primary imaging using the conventional crosscorrelation imaging condition is applied to the original data without consideration of the multiple attenuation and the result is shown in Fig. 13a. Because of the existence of multiples, crosstalk artifacts are introduced in the image under the depth of 3.3 km as shown in the black box. Fig. 13b is the simultaneous image with the estimated wavelet using the conventional crosscorrelation imaging condition. Similarly, crosscorrelation of the unrelated wavefields leads to heavy crosstalk artifacts (see the black boxes). However, more information can obviously be seen on the right side of the model, as shown in the white box, due to the wider illumination and higher fold from multiple reflections. Fig. 13c shows the simultaneous image using the proposed method described in this paper. From this figure, we can see that most of the crosstalk artifacts are well suppressed (see the black boxes) with the constraint of local slopes in the stereographic imaging condition compared with those of the image in Fig. 13b. The effectiveness of our method can be well demonstrated.

IV. CONCLUSION
In this paper, we first analyze the main challenges in the conventional simultaneous imaging of primaries and multiples through theoretical derivation, and it is concluded that one of the most important issues is the inaccurate wavelet, and the another issue is the artifacts that should be avoided by a proper imaging condition. Then, two creative strategies are presented to deal with these issues. First, to better describe and balance the contribution of the primary image and the multiple image, an estimated wavelet rather than an arbitrary wavelet together with the recorded data is used as the source wavefield. However, most of the current methods for wavelet estimation are realized based on signal theory under certain assumptions, and they cannot directly deal with the wavelet estimation for the data containing multiples. Therefore, we propose to estimate the wavelet through iterative SRME based on wave equation theory, which not only can avoid the assumptions in the convolution model based on signal theory but also is adaptive to the data containing multiples. Second, the stereographic imaging condition is introduced to deal with the heavy crosstalk artifacts generated during the simultaneous imaging. The conventional imaging condition is realized only by the time constraint, whereas the stereographic imaging condition is performed by multi-parameter constraint, with both local slope and time constraints of the forward-extrapolated and backward-extrapolated wavefields. According to the local slope difference of unrelated events, the crosstalk can be effectively suppressed. By using these two strategies, multiples can directly supplement the primary image, and an accurate imaging result with more information obtained from the multiples and with fewer artifacts can be obtained. According to the experiments and analysis of the numerical and field data examples, the feasibility and effectiveness of the proposed method is verified. As useful signals, multiples could be used not only for imaging but also for updating velocity. The waveform inversion of both primaries and multiples should be a promising study direction, and the proposed strategies, namely, wavelet estimation and artifact suppression, could be applicable as a reference.