Joint migration inversion: features and challenges

Joint migration inversion is a recently proposed technology, accommodating velocity model building and seismic migration in one integrated process. Different from the widely accepted full waveform inversion technology, it uses imaging parameters, i.e. velocities and reflectivities of the subsurface, to parameterize its solution space. The unique feature of this new technology is its explicit capability to exploit multiples in its inversion scheme, which are treated as noise by most current technologies. In this paper, we comprehensively evaluate the state-of-the-art joint migration inversion technology from various angles: we first benchmark its performance, on both velocity model building and seismic imaging, against that of the well-accepted workflow comprising full waveform inversion and reverse-time migration using a fully controlled 2D realistic synthetic dataset. Next, we demonstrate its application on a 2D field dataset. Last, we use another 2D synthetic dataset to clearly illustrate the challenges the current joint migration inversion technology is facing. With this paper, we transparently reveal the pros of cons of the current joint migration inversion, and we will also point out the imminent research directions joint migration inversion technology should focus on in the next phase for it to be more widely accepted by the geophysics community.


Introduction
Joint migration inversion ( JMI) is a recently proposed technology accommodating velocity model building and seismic migration in one integrated process (Berkhout 2014c), and it possesses flavours of both full waveform inversion (FWI) and imaging-based tomography methods (Verschuur et al. 2016). Compared to other velocity model building technologies and seismic imaging technologies, JMI is unique in two aspects: (1) multiples are treated as signals instead of noise and exploited explicitly for both velocity model building and seismic migration, and (2) velocity model building is addressed along with seismic migration in a closed-loop manner. Full wavefield modelling (FWMod), which is an acoustic wave propagation engine based on one-way operators (Sun et al. 2018), is a key component in JMI as it is capable of using imaging parameters (Berkhout 2014a), i.e. a velocity model and a reflectivity model, to carry out forward modelling. As a matter of fact, it is this unique modelling engine that makes explicit use of both surface and internal multiples possible in JMI. The migration scheme used in JMI is also based on FWMod, and is named full wavefield migration (FWM) (Berkhout 2014b). In terms of model updating, there are both similarities and differences between JMI and FWI. JMI and FWI are similar in their optimization schemes as they are both based on gradients of the models (Sun et al. 2019). Contrary to FWI, the forward modelling in JMI generates the wavefields with all orders of multiples included, while the backward propagation of the surface data via FWMod only considers the last bounce of any reflection event. In a nutshell, JMI presents a self-consistent logic to handle the problem of velocity model building and seismic imaging, and it aims to minimize the mismatch between the simulated data and the measurements via updating both a velocity model and a reflectivity model in a closed-loop manner. Similarities between JMI and some established seismic imaging and inversion technologies do exist, for instance, the velocity update in JMI is similar to the migration-based traveltime tomography (Clément et al. 2001) and the reflectivity update in JMI is similar to the machinery of least-squares migration (Nemeth et al. 1999).
Although in physics the reflectivity model should depend on concrete angles of the incident wavefield, in reality this will create an unacceptably large solution space for JMI, producing an over-parameterization of the inverse problem, which in turn renders the subsequent inversion process intractable. As a result, currently JMI only adopts an angleindependent reflectivity model in its implementation (Sun et al. 2019). Due to this approximation, amplitude versus offset (AVO) effects in the measured data cannot be correctly dealt with, and as a trade-off solution we have to apply JMI to the limited offset data.
Full waveform inversion (FWI) has been attracting attention in the field of oil and gas exploration for more than one decade, and it aims to use the original seismic data as much as possible to extract physical parameters of the subsurface (Virieux et al. 2014). One crucial application of FWI is to derive a high-quality velocity model and many successful applications have been reported so far, for instance Liu et al. (2014). Reverse-time migration (RTM) has been another well-established seismic imaging technology (Baysal et al. 1983) for years, but even now it is still one of the topics being actively researched (Fei et al. 2015). As RTM is based on two-way wave equations, it is capable of handling very complex imaging situations such as dealing with turning waves and imaging complex subsurface structures including salt bodies. In production, now it is common practice to combine FWI with RTM, i.e. the FWI-RTM workflow (Virieux & Operto 2009), to obtain seismic images with more accuracy and higher resolution.
For JMI to be better understood by the seismic community, in this paper we aim to reveal the pros and cons of the state-of-the-art JMI technology in depth, especially in comparison with the well-accepted FWI-RTM workflow. We will first use a realistic 2D deep-water model to carry out a comparison study between JMI and FWI-RTM, and via this study an application niche for the current JMI technology will be pointed out. With insights gained from this comparison study, we will apply the current JMI technology to a 2D deep-water field dataset to both demonstrate its application value and test its robustness in the real world. Finally, we will use another realistic 2D shallow water model to illustrate challenges and shortcomings the current JMI technology is facing, and this will also

JMI versus FWI-RTM on a 2D deep-water synthetic dataset
We use a realistic 2D deep-water synthetic dataset as the test bed for the comparison study between JMI and FWI-RTM. The velocity and density models used in this example are shown in figure 1. This velocity model is built via our known geological information, while the density model is derived from the velocity model using Gardner's equation (Gardner et al. 1974) except for the water layer where 1000 kg m −3 is used as the water density. Dimensions of our models are 10 km long by 6 km deep, and the water depth of this model is ∼500 m. A time domain finite difference method, based on first-order acoustic wave equations, is used for forward modelling. The source wavelet adopted in our simulation is a Ricker source wavelet with a dominant frequency of 15 Hz. For the data acquisition plan, we adopt the fixed-spread scheme: all receivers are located at 10 m below the water surface and cover the entire 10 km range of the model with an interval of 25 m; all sources are located at 5 m below the water surface with their x coordinates varying between 1 and 9 km with an interval of 100 m. The time sampling rate is 4 ms, and the recording length is 5 s. Due to the fact that RTM cannot easily handle surface-related multiples in imaging, in our forward simulation an absorbing surface condition is used in modelling so only internal multiples exist in our synthetic dataset. Our primary goal in this comparison study via this 2D realistic synthetic dataset is to benchmark the performance of the current JMI against the FWI-RTM workflow for both velocity model building and seismic imaging. The FWI used in this study exploits both refraction and reflection signals to derive the velocity model, and the frequency range used is between 3 and 10 Hz. The global correlation norm (Choi & Alkhalifah 2012) is used as the objective function in our FWI, and the source wavelet estimation is carried out in the frequency domain at every iteration (Kim et al. 2011). The zero time-lag cross-correlation-based imaging condition (Claerbout 1985) is used in our RTM, and the frequency range exploited by our RTM is between 3 and 40 Hz. In our FWI and RTM, the modelling engine is a time-domain finite difference method based on the second-order acoustic wave equation to avoid the inverse crime situation as our synthetic dataset is simulated via the first-order acoustic wave equations.
As pointed out, the current JMI cannot easily handle the AVO effects in the input data due to an angle-independent reflectivity model adopted in its parameterization (Sun et al. 2019), and hence in this study we limit the maximum sourcereceiver offset to 1.5 km for JMI. In other words, our current JMI actually only works on a limited-offset subset instead of the complete dataset. The frequency range used by our JMI is between 3 and 40 Hz. In addition, we treat the correct source wavelet as available a priori information for JMI. The initial velocity model provided to both JMI and FWI-RTM is shown in figure 2a, and compared to the ground truth, i.e. figure 1a, it is a very smooth starting model. Figure 2 parts b and c show velocity models derived by FWI and JMI, respectively. Clearly, FWI has a very strong capability to recover a velocity model close to the true velocity model. The velocity model derived by JMI is still very smooth, and compared to the initial velocity model, only the shallow part of the model is updated, but overall there is not much improvement gained from JMI on the initial velocity model. The limited-offset data used by JMI should be an important reason leading to this effect as it is usually the far offset data that can provide enough information to update the deeper part of the velocity model. Figure 3 shows migration images obtained by FWI-RTM and JMI. As shown in figures 2b and 2c, there exist huge differences in the velocity models for this migration step. As the FWI derived velocity model is much closer to the ground truth, the FWI-RTM image shown in figure 3a gives far more accurate subsurface structures than the JMI image, especially for the deeper areas. Figure 3c illustrates a comparison among the ground truth, the FWI-RTM image and the JMI image on a zoomed-in area, which is in the green box on  figure 3a and b. Clearly, it can be observed that the JMI result is quite off in depth. Furthermore, the shallow part of the JMI image looks wobbly and more smeared. As a matter of fact, all these unsatisfying features of the JMI image are due to the inter-dependency between the velocity model and the reflectivity model in JMI. As now the velocity model is far from the ground truth, a 'more wrong' image, which in fact is the reflectivity model in JMI, does help JMI to have a lower data mismatch between the measurements and the simulated data. In other words, the inversion result of JMI is trapped in a local minimum. Based on this comparison study, we can conclude that the current JMI using only limited-offset data is inferior to the well-established FWI-RTM workflow. Because the AVO effects cannot be modelled correctly in the current JMI, to mitigate this simulation inaccuracy in our inversion, we then have to use a limited-offset subset of the complete dataset, but this in turn comprises JMI's capability for velocity model building. To better use the power of the current JMI technology in handling multiples, we now adopt the mentality of bootstrapping by directly using a high quality velocity model, for instance the FWI derived velocity model shown in figure 2b, as the initial velocity model for our JMI. Surprisingly, with such a high quality initial velocity model in place, our current JMI, although still only based on the limitedoffset data, shows the capability to further update fine details in the velocity model, as shown in figure 2d. For example, in our true velocity model shown in figure 1a, there exist some lateral blocky features as circled out by the two green circles. These lateral blocky structures are not recovered by FWI as shown in figure 2b, but now are clearly visible in the FWI-JMI velocity model shown in figure 2d. In addition, as now the FWI-JMI velocity model is much closer to the ground truth than the previous JMI velocity model shown in figure 2c, the migration image from FWI-JMI, shown in figure 4a, also looks much improved: all subsurface structures are imaged accurately with a better vertical resolution than the FWI-RTM image shown in figure 3a. We further use the FWI-JMI velocity model as the migration velocity model for RTM, and the corresponding image is shown in figure 4b. Comparing the image from FWI-JMI to both FWI-RTM image and FWI-JMI-RTM image, we can find that the FWI-JMI image shows fewer artefacts and higher resolution than the other two images, suggesting that internal multiples are better accommodated by JMI to improve the final image quality in this limited-offset situation. For instance, events below the deepest authentic reflector, as circled on figure 4b, are purely cross-talk noises due to the improper handling of internal multiples by RTM. As shown in figure 4a, all these erroneous events are excellently suppressed by JMI. In addition, figure 4a also bears more accurate estimations for depth information of reflectors than 5 figure 3a, and this is due to the uplift brought by JMI to the velocity model.
Although we believe that the improvement of image resolution from our FWI-JMI workflow is the result of the better handling of internal multiples, there still exists another possible explanation that this resolution improvement is from the fact that only the limited-offset data is used in JMI. After all, an inaccurate migration velocity leads to poorly aligned common image point (CIP) gathers, and hence stack of the full-offset CIP gathers in turn smears out events compromising the image resolution. To rule out this potential cause and justify the value of exploitation of multiples for image resolution enhancement in JMI, even only using the limitedoffset data, we decide to carry out a dedicated comparison study for imaging in the angle domain between RTM and FWM (the imaging part of JMI) with the complete dataset used and the migration velocity model fixed at the FWI-JMI velocity model shown in figure 2d. Our angle gather reverse time migration (AGRTM) calculates the angle gathers via the optical flow idea. To address the AVO effect in the complete dataset in FWM, an angle-dependent reflectivity model in the linear Radon domain, i.e. the − p domain, is adopted (Davydenko & Verschuur 2017), and in this paper we refer to this imaging scheme as angle-dependent full-wavefield migration (ADFWM). Since ADFWM uses the ray parameter p to describe the wavefield propagation direction at each subsurface layer in migration, the wavefield propagation direction thus defined becomes a function of both the ray parameter p and the medium velocity at the point of interest in the subsurface: where p is the ray parameter, v is the medium velocity and is the 'apparent' propagation direction (or 'apparent' propagation angle), which is with respect to the vertical direction. In Davydenko & Verschuur (2017), they directly used ray parameter gathers as equivalent angle gathers, but in this paper because we need to compare our results against those of AGRTM, we have to convert ray parameter gathers to 'apparent' propagation angle gathers via = sin −1 (pv) .
As our AGRTM uses the optical flow idea to calculate angle gathers, the angle defined in that case is actually only between 0 and 90º as that angle is a specular reflection angle per se. To better match those angle gathers in AGRTM, in ADFWM the angle defined by equation (2) is treated by its absolute value during the stacking process, for instance, the ADFWM image covering the angle range of [20, 30] degrees stacks all the angle gathers between −30 and −20º and between 20 and 30º. We would like to point out that small differences exist between AGRTM images and ADFWM im-ages due to the inherent difference in the definition of angle gathers in AGRTM and ADFWM. Figure 5 parts a, c and e shows the migration image of AGRTM in the angle ranges of [0, 10], [30,40] and [60,70] degrees, and figure 5 parts b, d and f shows the corresponding angle-domain image from ADFWM. From these comprehensive comparisons, several effects can be clearly observed: results of ADFWM do look very similar to those from AGRTM, and this shows that using a ray parameter-based reflectivity model is a feasible way to carry out migration in the angle domain; ADFWM images have better consistencies at high angles, and this is the result of its proper handling of both primaries and multiples; ADFWM images look sharper than the corresponding AGRTM images, which is also the result of its better handling of multiples so that energies are better collapsed to correct events rather than to smear them out during the imaging process. In addition, as the inversion engine in FWM/JMI contains a feedback loop, which is like the one used in the least-squares migration process, we believe this also helps to improve the image resolution. Figure 5 should be convincing evidence to support our conclusion that resolution improvement of the FWI-JMI image in figure 4a over the corresponding RTM image in figure 4b is from JMI's better handling of multiples.
Our systematic benchmarking tests between JMI and FWI-RTM on this 2D deep-water synthetic dataset demonstrate that the current JMI technology, even though it can only deal with limited-offset data, is still a powerful technology for both velocity model building and seismic imaging in the deep-water scenario as long as a high-quality initial velocity model is available. In reality, this high-quality initial velocity model could be derived by FWI, migration velocity analysis (MVA) or some other available means.

JMI on a 2D field dataset
With the experience gained from our previous 2D synthetic deep-water model, in this section we apply our current JMI to a 2D deep-water field dataset. Data preconditioning is first carried out to prepare this dataset for JMI, and relevant processing manipulations include noise attenuation, removal of surface-related multiples, redatuming the data to the surface level, source wavelet estimation, reconstruction of the missing near-offset data and creating split-spread records via reciprocity. The x dimension of this survey is 10 km long. After data processing, receivers cover the complete x range with an interval of 20 m. Sources are placed between 1 and 9 km in the x direction with an interval of 100 m. Time sampling rate is 6 ms, and the record length is 4.5 s. The valid frequency range of this dataset is between 3 and 40 Hz. As in this dataset the near source-receiver offset part (from −150 to 150 m) is reconstructed, which unavoidably brings uncertainties and errors into the input wavefields for JMI, we use a slightly longer offset of 2.5 km as our offset limit to select the limited-offset subset. Figure 6a shows several shot gathers with their corresponding sources around x = 5 km, and the complexity of this data can be clearly observed. In this application of JMI, we use a velocity model derived by MVA as the initial velocity model where a big salt body is in existence, as shown in figure 6b. Our JMI updated velocity model is shown in figure 6c, and it can be seen that JMI not only derives a new velocity layer below the top of salt, which is located at ∼1 km deep on figure 6c and pointed at by a green arrow, but also seems to recover more velocity details inside the salt body.
The JMI image is shown in figure 6d. As a comparison, we carry out RTM using the JMI updated velocity model shown in figure 6c for migration; the image is shown in figure 6e. From these two images, two effects can be observed. First, our JMI image is more consistent with our JMI updated velocity model, and at the new velocity layer around 1 km deep (pointed out by a green arrow on figure 6d) we see a corresponding structure. In the RTM image, no such structure exists, which may indicate that some signals are not handled correctly by RTM. Second, for the event around ∼4 km deep, the JMI image shows much better event continuity than that of the RTM image (pointed out by red arrows on figures 6d and 6e). We believe these two differences reflect the contributions of multiples and transmission and reflection effects, which are properly handled by JMI but neglected by RTM.
To justify the correctness of our JMI results of this 2D field dataset, we refer to a well log from the nearest available production well, which is ∼12 km away from the area where this 2D field dataset is acquired. The relevant part of this well log is shown on figure 6f. Furthermore, we squeeze it proportionally to the depth scale of figure 6c and put it at the relative location on that picture, and readers should be able to see a relatively good match between our inverted velocity model and the well log. According to this well log, we indeed successfully recover the high velocity layer below the top of the salt body-in our JMI results, both the velocity model and the migration image successfully reveal the existence of this layer. Overall, we believe the application of JMI on this 2D field dataset gives very encouraging results, and it also shows the consistency of JMI in both velocity model building and seismic imaging.

JMI on a 2D shallow water synthetic dataset-a challenging scenario
So far, we have demonstrated some success of JMI on both a 2D realistic deep-water synthetic dataset and a 2D deepwater field dataset. Nevertheless, as we introduced before, at this stage actually JMI faces a big challenge to handle AVO effects correctly in the data. In this section, we will use a 2D realistic shallow water synthetic dataset to give readers some indications on this challenge. Figure 7a shows the true velocity of this 2D realistic shallow water model and figure 7b shows the corresponding density model, which is built via Gardner's equation (Gardner et al. 1974). Figure 7c shows the shallow part of the true velocity model, where the water bottom topography can be observed and the water depth is ∼30 m. This 2D model is built from true well logs, and it represents a typical geological situation in the Middle East. Dimensions of our models are 10 km long by 6 km deep. The data modelling scheme used in this example is exactly the same as that in our previous 2D deep-water synthetic example. Figure 7d shows the first 3 s of the central shot gather of this dataset. Due to the shallow water situation of this model, strong AVO effects, as pointed out by the green arrows, even in the near-offset part of this shot gather, can be observed easily. As pointed out by the green arrows on figure 7d, the amplitudes of those two events are weak around x = 5 km but get stronger as the source-receiver offset increases. The initial velocity model for JMI is shown in figure 7e, and it is built by smoothing the true velocity model in the z direction. All major structures of the true velocity model have been captured by our initial velocity model, so it actually is a high-quality enough initial velocity model for JMI. For this dataset, we set the limit of the source-receiver offset at 1.5 km for our JMI, and its updated velocity model is shown in figure 7f. As pointed out by the green arrows, quite some erroneous structures are derived by our JMI process, which is mostly due to the AVO challenge our JMI is facing. As this is a shallow water dataset, AVO effects unavoidably exist even in the near-offset part of the dataset, so, different from the deep-water situation shown 12 before, it is not adequate enough for us to just straightforwardly select a limited-offset dataset to feed into our current JMI. Figure 7g shows the JMI image, and due to the erroneous velocity model shown in figure 7f, erroneous migration events are also visible as indicated by green arrows.
The shallow water scenario is a situation where strong AVO effects are almost everywhere in the data, and the easy trick of selecting a limited-offset subset of the complete dataset is no longer a valid workaround for JMI. This controlled 2D synthetic test transparently reveals the challenge faced by our JMI in a clear manner, and we believe it is also a good test bed for our future work to improve the JMI technology.

Conclusions
In this paper, we comprehensively reviewed the status of the current JMI technology using two 2D realistic synthetic datasets and one 2D field dataset. Due to the assumption that reflection/transmission effects are incident-angleindependent, which helps to avoid the over-parameterization problem, the current JMI technology cannot properly handle the full-offset dataset due to the AVO effects in the data, so our current workaround is for JMI to only work on a limited source-receiver offset subset of the complete dataset. By using a 2D realistic deep-water synthetic dataset, performance of the current JMI is benchmarked against that of the well-established FWI-RTM workflow for both velocity model building and seismic migration, and we revealed one niche for the current JMI technology, which is to use a high-quality initial velocity model for the current JMI to start with in the deep-water scenario. By using the FWI velocity model as the starting point for JMI in this synthetic test, our current JMI not only further updated fine details in the velocity model but also produced a high-quality migration image with much sharper resolution and fewer artefacts from internal multiples. To support our claim that the improved resolution from the JMI image is the result of internal multiples being better handled instead of limited source-receiver offset, we also carried out an imaging comparison study in the angle domain on the full-offset data between ADFWM and AGRTM using the FWI-JMI velocity model as the migration velocity model. We further applied our current JMI technology to a 2D deep-water field dataset, and results are encouraging. Finally, we used a 2D realistic shallow water synthetic dataset to clearly reveal the AVO challenge faced by the current JMI technology: erroneous velocity and image results were derived by our current JMI as the AVO effects in the data cannot be correctly accounted for. We believe with this paper we have convincingly demonstrated the pros and cons of the current JMI technology, and this paves a solid way to further develop this technology.
At this moment, we see several potential directions to resolve the AVO challenge faced by the current JMI technology. One is to adopt a different than least-squares dataresidual based objective function in the current framework of the JMI technology, and some initial efforts show improvement compared to the current JMI technology. Another is to fully honour the physics of wavefield propagation via more accurate propagation and reflection/transmission operators, and some initiatives are also being investigated. In addition, using image-domain criterion for velocity updating might also be a direction to partially mitigate the AVO challenge in JMI. We believe all these directions deserve to be looked at carefully.
Conflict of interest statement. None declared.