Seismic imaging is a method for generating images of the subsurface of the Earth without the need to drill down to observe it directly. Using a seismic source, high-amplitude acoustic waves are generated, which propagate through the Earth. When these waves encounter interfaces between different layers, part of their energy is reflected, traveling back upwards. At the surface, these scattered waves are recorded using an array of sensors. These recorded wavefields are then used as input for a seismic imaging algorithm, which attempts to reconstruct an accurate image of the subsurface based on these measurements.
While many different seismic imaging methods exist, this thesis focuses specifically on Full-Wavefield Migration (FWM). FWM is a least-squares migration method, based on iteratively updating the image in order to minimize the misfit between the recorded data and the forward modelled wavefield at the surface. Using this iterative approach, FWM is able to incorporate multiple scattering effects into the imaging process, increasing the accuracy of the resulting images. Also, by using explicit convolutional operators based on the one-way wave equation, the computational cost of FWMis reduced compared to alternative iterative imaging methods based on finite-difference modelling.
In this thesis, we describe a number of recent advancements to the FWM method. Our main focus is the extension of the FWM method to include converted waves. In order to take these effects into account, we need accurate reflection and transmission operators. However, the true, elastic reflection and transmission operators are notoriously non-linear, making them difficult to work with. Therefore, we introduce a novel set of approximations of these operators, which we name the extended Shuey approximations. To benchmark these approximations, we apply them in a simple, 1.5D scenario. This test shows that the extended Shuey approximations yield improved results for forward modelling and imaging, compared to the conventional Shuey approximation.
We then use these extended Shuey approximations to derive accurate reflection and transmission operators for the 2D case. Combining these operators with the existing theoretical framework of FWM we develop a robust imaging algorithm which takes converted waves into account. We then apply this algorithm, which we name elastic FWM, to two synthetic models. We show that, for these models, the elastic FWM method out-performs the conventional, acoustic FWM method. We also demonstrate that, using this method, additional information regarding the elastic medium parameters of the subsurface can be recovered.
Finally, we examine the known issue of slow convergence for the conventional, acoustic FWM algorithm. We introduce a novel preconditioner, based on approximating the pseudo-inverse using Proper Orthogonal Decomposition (POD). Using this preconditioner, we demonstrate improved convergence for the synthetic Marmousi model and a field data set from the Vøring basin in Norway.

The introduction and adoption of seismic full waveform inversion (FWI) revolutionized Earth's subsurface imaging practices. FWI uses all the information in the seismic data (amplitudes and phases) to reconstruct a detailed Earth's subsurface model. However, it does come with limitations. Beyond the reach of refractions and diving waves, FWI cannot effectively reconstruct subsurface layers. This led to the development of reflection waveform inversion (RWI), which exclusively uses the pair of transmission-after-reflection wavepaths to sample deeper compared to FWI. RWI reconstructs the background velocity model of the subsurface by alternating between a migration loop and a tomography loop.

Despite the theoretical appeal, RWI has its own share of limitations. This dissertation investigates the barriers limiting the optimal performance of reflection waveform inversion in the context of one-way RWI (ORWI), a variation of reflection waveform inversion that adopts one-way wavefield propagators to forward model seismic reflection data. After exploring the barriers, the dissertation offers solutions to improve the reliability, accuracy, and convergence of ORWI.

The dissertation acknowledges several barriers that limit the optimal performance of conventional/standard ORWI. First, ORWI relies on accurate subsurface images. However, limited-resolution images with unpreserved amplitudes, resulting from the migration loop, lead to suboptimal background velocity updates. This issue also extends to the tomography loop, where limited-resolution tomographic wavepaths impede optimal background velocity updates. Second, ORWI overlooks the immediate impact that updating the velocity model has on the reflectivity model, as the reflectors' positions in depth remain fixed while the background velocity is updated. This oversight leads to inconsistent reflectivity and velocity models in ORWI, introducing full-wave inconsistencies in the short-offset residual waveforms for tomography. Third, similar to other seismic waveform inversion techniques, ORWI suffers from the detrimental effect of including cycle-skipped data from long offsets.

To mitigate the barriers, the dissertation proposes a range of solutions. Initially, the dissertation introduces a computationally efficient high-resolution migration algorithm called preconditioned least-squares wave-equation migration (PLS-WEM) through depth-dependent gradient preconditioning. PLS-WEM reconstructs high-resolution, amplitude-preserved seismic images in fewer iterations.

Following that, by incorporating PLS-WEM into standard ORWI, the dissertation enhances ORWI, achieving improved reflectivity imaging and thereby reconstructing tomograms that are more representative of the true subsurface layers. The dissertation also proposes the following data solutions: (a) Muting short-offset residual waveforms in the tomography data to reduce the adverse imprint of inconsistencies between the reflectivity and velocity models on the tomographic gradient of ORWI. (b) Building on (a), extending the migration offset to the maximum effective migration offset (MEMO) to enhance both the signal-to-noise ratio and the illumination of the reflectivity model. (c) Introducing a data selection algorithm to minimize the impact of cycle-skipped long-offset data.

The dissertation then presents high-resolution ORWI (HR-ORWI) technology, which leverages depth-dependent gradient preconditioning in both migration and tomography loops to reconstruct optimal tomograms in fewer cycles.

The dissertation next evaluates three approaches to depth-dependent gradient preconditioning: conventional, source-interference-free, and source-interference-inclusive. Numerical results show the superiority of the source-interference-inclusive approach, offering enhanced resolution, reduced computational demands, and resilience to source interference.

Lastly, the dissertation develops a mathematical framework that integrates early-arrival waveform inversion with ORWI through the subspace gradient method, combining the strengths of transmission and transmission-after-reflection wavepaths to enhance tomogram reconstruction.

In conclusion, this dissertation offers a comprehensive set of solutions to overcome the limitations of ORWI, facilitating its broader adoption and application in seismic exploration and velocity model building.

Seismic data reconstruction addresses the challenge of accurately restoring incomplete or damaged seismic datasets, which is crucial for subsurface imaging and exploration. The missing data is often due to equipment failure, signal loss, environmental noise, or poor geological conditions. It may also result from limitations in data acquisition, such as sparse sampling and noise interference. This thesis proposes a Clean Convergent alternating Projections onto Convex Sets (CCP) method for data reconstruction. By incorporating an initial value tweaking step into the alternating projections process in the Convergent alternating Projections onto Convex Sets (CP) method, this method aims to reduce ringing noise in the reconstruction results by the CP method. Initially, the traditional CP method and its application in seismic data reconstruction are introduced, highlighting its shortcomings in addressing ringing noise. To overcome this issue, a CCP method based on a non-local means algorithm for initial value tweaking is proposed and applied within the outer loop of the CP method. The theoretical foundation and implementation steps of the CCP method are discussed in detail, and its effectiveness is validated through a series of experiments. The experimental results demonstrate that, compared to the traditional CP method, the CCP method significantly reduces ringing noise and improves the quality of the reconstructed data. Various datasets, including images, 2D seismic section data from the SEAM II Arid model, and 3D seismic model cubes, are used to showcase the reconstruction capabilities of the CCP method in different scenarios. Finally, the thesis provides an in-depth discussion of the CCP method, including intermediate results in the outer loop, the necessity of data preconditioning, and parameter testing, offering a practical set of control parameter settings. The study shows that the CCP method has broad application prospects in data reconstruction, providing valuable references for future research.

Imaging and inversion with seismic data recorded with sources and receivers at the surface are powerful tools to infer knowledge about the subsurface. However, creating an image with seismic data is unfortunately not as easy as taking a picture with a smartphone. The estimated subsurface models in many situations are far from ideal due to the low quality nature of the data. One of the reasons can be weathering of the near-surface geology that generates unconsolidated material characterized by slow velocity with rapidly varying, heterogeneous and season-dependent nature. Acquiring seismic data on such near-surface leads to complex wave propagation, posing challenges to imaging and inversion. In this dissertation, we tackle the weathering effects during seismic data processing, imaging and inversion with low-rank-based methods.

One approach to tackle the weathering effects on seismic data is removing them during seismic data processing. To do so for 2D data, we propose a model-independent low-rank-based near-surface estimation and correction in the midpoint-offset-frequency domain. In this domain, ideal data exhibit low rank structures, which get destroyed due to the influence of the weathering layers. Accordingly, the method makes use of the redundant nature of seismic data that allows for accurate approximation by low-rank matrices. To estimate the time shifts that compensate for the weathering effects, we cross-correlate a data set influenced by the near-surface weathering layers with its low-rank approximated version. Since we estimate time shifts (commonly referred to as statics) and no longer the directly low-rank approximated data, we avoid losses of the amplitude information. To improve the estimated statics and to alleviate the need for accurate rank selection for low-rank approximation, we implement the method in an iterative and multi-scale fashion. Since the low-rank approximation deteriorates at high frequencies, we utilize its better performance at low frequencies and exploit the common statics amongst different frequency bands. Using synthetic and field data, we demonstrate the performance of the proposed proposed, which requires no knowledge of the subsurface model, demands minimal data pre-processing, and provides accurate solutions with high computational efficiency compared to existing techniques.

When seismic data acquired on complex near-surface are additionally subsampled for economical reasons, such as monitoring of sequestrated carbon dioxide and hydrogen, the problem is further exacerbated. Both the weathering layers and randomized subsampling render coherent energy incoherent. Therefore, they both contribute to destruction of the low-rank structure commonly associated with statics-free densely-sampled data. Frugal data acquisition in complex near-surface regimes makes separation of the distinct sampling and weathering effects on the rank structure difficult, which as a result lead to poor reconstruction. To overcome that, we propose to reconstruct the data with joint rank-reduction-based near-surface correction and interpolation. The method simultaneously accounts for the weathering and subsampling effects to provide accurate reconstruction. Since low-rank approximation is used for near-surface correction, we also utilize it in rank-minimization interpolation as a cost-free initial solution to the optimization problem. As both near-surface correction and interpolation operate in the midpoint-offset domain, we avoid the cost of transformations back and forth from the source-receiver to midpoint-offset transform domain. Consequently, the proposed reconstruction, which shows its potential on synthetic and field data, additionally increases the computational efficiency.

While the aforementioned near-surface correction deals with 2D data, the Earth is a 3D object that requires acquisition of 5D data for proper subsurface model estimation. For 5D data, the limitations and challenges of conventional near-surface correction methods are magnified. To avoid them, we propose a 5D model-independent low-rank-based near-surface correction. To compute the singular value decomposition of 5D data volumes with 1 temporal and 4 spatial dimensions, which is necessary for low-rank approximation, we need to perform matricization of the 5D data, i.e. organization of the 5D data into matrices. At the same time, it is essential that the chosen organization domain reveals the underlying low-rank structure. Therefore, we first analyze different matricization domains that can be used to organize the 5D data. Similar to the 2D case, we show that --- in the potential domain --- the near-surface weathering layers render coherent energy incoherent, which results in slowly decaying singular values compared to the statics-free data that are of low-rank nature. The proposed method, which we show on synthetic and field data, enjoys the same benefits of the proposed method for 2D data, in addition to being able to capture the 3D nature of the Earth.

Due to the complex nature of the near-surface and due to its impact on the subsurface model, the near-surface model gets treated separately from the subsurface model. However, the optimal goal is not to remove the near-surface effects with data processing, but to accurately estimate near- and sub-surface models simultaneously. To do so, we use the inherent scale separation of joint migration inversion that estimates a low-wavenumber velocity and high-wavenumber reflectivity. Since rapid variations in surface elevation and near-surface model result in high wavenumber effects, they end up affecting the reflectivity model. At the same time, the estimated reflectivity influences velocity estimation. Consequently, JMI provides erroneous subsurface models in the presence of complex weathering layers. To mitigate that, we use multi-scale low-rank updates in the reflectivity domain. The proposed method reduces the near-surface effects at the initial iterations, but it allows more details of the near-surface model to enter the solution at later iterations. In the end, we estimate accurate near- and sub-surface models simultaneously without the need to bypass the weathering layers.

Seismic survey design deals with determining the acquisition parameters that lead to the best possible imaging and characterization of the subsurface. The design of the survey is constrained by health, safety and environmental considerations and the available budget, seeking for a balance between quality and cost. Because seismic exploration is a widely used geophysical method for revealing underground resources, information about the subsurface is available in many areas. Therefore, it can potentially be used for purposes supplementary to exploration such as the monitoring of producing fields and fluids injection. However, the available budget for these purposes is usually lower than for exploration, and it becomes a priority to maximize the benefits derived from a potentially cheaper acquisition. In this thesis, we propose new methods for the analysis and design of seismic surveys that are based on previous knowledge from existing subsurface models and that aimto maximize image quality with the lowest acquisition efforts.

Each violin generates a unique sound with a distinct set of acoustic properties. In particular, the Stradivarius violins are said to have a very unique sound. The ultimate goal of this thesis is to aid in recreating the unique sound of the Stradivarius. Though, the field of research surrounding the violin and its acoustic properties is very broad and many aspects of this research, such as subtle structural differences between violins, fall beyond the scope of this thesis. Rather, the goal of this thesis is to provide the tools to find the differences between different violins. This thesis is split into two parts. In the first part, a literature review, describing the different parts of the violin, is presented along with an overview of different experiments and finite element models, describing the plate dynamics of the violin. In the second part, a model of the violin is devised, with modelling cycles based on consecutive simplification, and implemented in Python using finite difference equations and the method of lines. This model is made with the goal to provide a framework for creating a realistic model of the violin where differences in timbre between violins can shown. In the literature review, four subjects are addressed. First, an explanation on timbre is given with an illustration of its perception. Additionally the timbre of historical and contemporary violins is compared. Second, the radiated sound of the violin is discussed. In specific, the problem with directivity of radiation is brought to attention. Third, an overview of the bridge is given with a strong emphasis on the influence of the bridge admittance on the radiated sound. Last, the corpus is dissected into different parts which are investigated with optical sensors and modelled with the software packages Comsol and Abaqus. The model in the second part of the thesis is based on principle of consecutive simplification. To set the groundwork for further iterations of the model, the most important assumptions are kept track off. The method of lines is then used to simulate the vibrations of a wooden square plate using the flexural wave equation as the governing equation of motion. In this simulation, the forward Euler method is used for the numerical time integration. The steady state solutions of the numerical model are coherent with the analytical solutions to the flexural wave equation. However, the time dependent results diverge from the results shown in the literature study. This is either due numerical time integration method being unstable or due to the physical properties of the wood not being translated by the stiffness matrix in the model.

Weather and climate models require knowledge about the surface temperature and surface heat flux to make predictions for climatology, farming, wind power prediction, etc. Most current models make use of Monin-Obukhov Similarity Theory (MOST). This theory expresses multiple turbulent fluxes in terms of their mean gradients with height. MOST results in an accurate universal form for the non-dimensional wind, temperature and humidity profiles from just above the surface and up to the top of the Atmospheric Surface Layer. The description of the near-surface temperature profile above a rough surface is, however, non-robust and ill-defined in the lowest layer of the atmosphere, the Roughness Sublayer. This thesis presents a new flux-gradient framework based on surface geometry. We apply it to the near-surface temperature profile of high-resolution experimental data from fibre optic measurements over a grass surface. The framework was developed to universally compare 3 different types of models at the surface: a Direct Numerical Simulation in a smooth pipe validated model; a constant surface gradient over grass model and a tall canopy model. Where the temperature is scaled with the already in use turbulent heat flux scale, the new length scale scales the height based on a finite surface gradient which is observed in experimental data. The experimental data can in this framework collapse onto one universal profile with the use of a reference temperature taken at the top of the grass. The framework thus extends MOST with a universal non-dimensional temperature profile in the roughness sublayer.

Seismic reflection imaging aims to generate a representation of the subsurface of the earth using acoustic or elastic waves recorded in the form of seismic data. During the processing of the data for imaging, we pick a part of the data as signal and discard the rest as noise. Since so-called primary wavefields carry the single-scattering reflection response of the subsurface, they are assumed to be sufficient to carry out reflection imaging. This also means that significant effort is devoted towards removing not only the noise but also the multiple reflection events to avoid imaging artefacts, also known as cross-talk. Surface-related multiples are the multiples generated during the marine acquisition by at least one downward reflection at the water-air boundary, and tend to be the strongest in amplitude compared to the other multiples. Over the past several decades, many novel techniques have been developed to remove the surface-related multiples effectively. While we are getting better at primary-multiple separation, it is still a very challenging and expensive problem. In recent years, multiples have gained recognition as valuable signals, not just noise. Multiples also contain the reflection responses of the subsurface and since they travel different paths they often contain additional information about the subsurface compared to the primaries-only wavefields. Imaging with primaries and multiples without separation is the way forward as it avoids the expensive multiple removal steps along with providing (potentially) additional illumination from the multiples...

The goal of seismic imaging is to create a model of the subsurface from samples of a transmitted wavefield that is reflected in soil layers.
The created model can be used to locate storage possibilities for CO2 or H2 or to find oil and gas reservoirs or other natural resources.
Seismic imaging relies on the propagation of wavefields, which can be computed by evaluating the Rayleigh integral, a process that often requires a lot of time and resources.
In this thesis we develop methods to reduce error in the computation of the Rayleigh integral whilst remaining time-efficient.
To achieve this we first give the derivation of the Rayleigh integral and provide a few methods to evaluate integrals numerically using interpolating polynomials.
Afterwards, we adapt the integration methods to the Rayleigh integral and arrive at the main algorithms for its evaluation, at the end we compare these algorithms.

We present the method that is currently used for the computation of the Rayleigh integral and a modification that extends the method to samples on a rectilinear grid.
The adjusted version of this method did not perform consistently in our results, although it achieved better results than the original method on average.
Secondly we present a combined interpolation method (that is, it interpolates the full integrand without taking advantage that one of the two functions in the integrand is known) named Simpson's rule, which only performed better in the high sample points per wavelength range with minimal noise.
An alteration to the method for semi-equidistant grids also did not perform consistently and is not worth the extra operations.
Next was the separate interpolation algorithm (this time taking advantage of the known part of the integrand), this method is difficult to implement and had only slightly better performance than the combined method.
The final method was an extension of the separate interpolation algorithm to non-equidistant grids, this method performed undoubtedly the best of all presented methods (often better by a factor of 15 compared to the currently used method), this method is therefore recommended for evaluating the Rayleigh integral.
The amount of operations needed for all presented methods is approximately the same and the execution time needed for the calculation depends on the implementation of the problem.

In conclusion, we found that the separate interpolation method for non-equidistant grids delivered the most accurate approximations of the Rayleigh integral. This method can be used in seismic imaging to increase the accuracy of models of the subsurface.

In this thesis we have analysed the behaviour of a physics informed neural network and it’s competence in predicting a wave in a non-homogeneous medium. During this project we have used a fully connected network with labelled input data of a 2D acoustic wave. On top of this we used a special loss function that calculated whether the output of the network satisfies the wave equation. Our experiment consisted of the tuning of the hyper parameters, analysing the optimal choice of activation function and the optimisation of the input data and improving the loss function. During this project the unpredictable nature of machine learning has become very clear. We have experimented with several activation functions and have found that the optimal choice of activation function depends on how long you are willing to train the network, as the development of the loss function differs immensely between activation functions. When we looked at the optimal scaling of the input values we find that a non-trivial scaling seems to work better than for example, normalisation of these values. Furthermore we have tried to improve the sampling of the points we use to calculate whether the prediction of the neural network satisfies the wave equation and got interesting results. When we implement all op- timisation techniques, we find that the neural network is extremely capable of predicting the wave’s behaviour in a high contrast media within the time frames of the input data. Prediction outside of this time frame does work but the results do deteriorate especially in the positive time direction. Predicting in the negative time direction yield slightly better results.

Conventional imaging techniques from medical ultrasound and seismic imaging are not suitable for multilayered domains where there is no separation between received reflection events. This happens in layered structures where the thickness of the individual layer is small with respect to the induced wavelengths. An example of these structures is the wall of a typical blast furnace. Imaging these domains require a method based on the interference pattern in the frequency domain. This method is already widely used for imaging concrete structures. However, these consists of only one or two layers. This study introduces a FullWaveform Inversion scheme which extends the usability of themethod to multilayered domains. Itminimizes the error between measurements and a forward model. To reduce the computation load, a simplified forward model based on a plane wave response in a vertically layered medium is adopted. Pre-processing steps are necessary to remove the effect of the source wavelet and to enhance the peaks and troughs positions. A sensitivity analysis provides information about the usability of the method. This shows only layers with standing wave frequencies in the range of the induced source wavelet can be imaged. Furthermore, results show it is impossible tomeasure the wave velocity and thickness of a single layer simultaneously given our simplified forward model. Multichannel measurements over an aperture of 500 mm are required to transform a spherical wave response to its plane wave equivalent in order to fit the forward model. Tapering the sides of the aperture is essential to reduce simulation artefacts. The inversion scheme is successfully tested on synthetic data of a blast furnace wall in the acoustic and elastic domain using realistic simulations of the measured data.

During the past decade, time-lapse seismic technology has been widely applied in hydrocarbon reservoir management. It is a very powerful method to obtain information on reservoir changes in the inter-well regions. This information helps to identify bypassed hydrocarbons and extend the economic life of a field. In a typical scenario, one baseline survey and subsequent monitoring surveys are acquired over time. The survey geometry is usually exactly repeated and well-sampled to mitigate acquisition effects on the next steps in the process. By processing and comparing all the datasets, some physical changes, e.g. reflection amplitude and travel-time changes, can be estimated. These time-lapse changes are then used to calculate interpretable parameter changes in dynamic reservoir rock and fluid properties, e.g. pore pressure and fluid saturation.

In a conventional time-lapse processing workflow, all the multiples are first removed from the data, then independent imaging process is employed to each dataset, given the same propagation velocity model. Later on, to compensate the ignored velocity variations between different surveys, a time-shift map (travel-time differences) is estimated from the calculated images and then applied back to them, yielding the final reflection amplitude differences. However, this conventional processing strategy is usually sensitive to the success of multiple removal and survey repeatability, and also requires well-sampled surveys providing proper illumination. Moreover, artifacts are often generated in addition to the actual time-lapse changes due to the non-repeatable uncertainties during the independent processing steps. Regarding the time-shift-map tool, the relative velocity changes derived from the time-shift map are not the actual velocity changes due to its local 1D subsurface assumption that is embedded.

In order to relax these rigid requirements and have a better velocity change indicator, we propose Simultaneous Joint Migration Inversion (S-JMI) as an effective time-lapse tool for reservoir monitoring, which combines a simultaneous time-lapse data processing strategy with the Joint Migration Inversion (JMI) method. JMI is a full wavefield inversion method that explains the measured reflection data using a parameterization in terms of reflectivities and propagation velocities. JMI is able to make use of multiples and at the same time take velocity variations between surveys into account. The simultaneous strategy, which means fitting all the datasets simultaneously, allows the baseline and monitor parameters to communicate and compensate with each other dynamically during inversion via L2-norm constraints, thus, reducing the non-repeatable uncertainties during the time-lapse processing workflow. As a result, more accurate time-lapse differences can be achieved by S-JMI, compared to inverting each dataset independently. Moreover, in order to get more localized time-lapse velocity differences, we further extend the regular S-JMI to a robust high-resolution S-JMI (HR-S-JMI) process by making a link between the reflectivity/reflectivity-difference and velocity/velocity-difference during inversion. With a complex synthetic example based on the Marmousi model, we demonstrate the performance of the time-shift-map-based method, sequential JMI, the regular S-JMI and HR-S-JMI is improving in this particular order.

Next, we further demonstrate the effectiveness of the proposed method in more real-life cases with a highly realistic synthetic model based on the Grane field, offshore Norway, and a time-lapse field dataset from the Troll Field. Moreover, in order to investigate the feasibility of HR-S-JMI in practice, several numerical experiments based on the realistic Grane model are conducted, regarding the following aspects: noise, including random noise and coherent noise caused by the acoustic assumption; the quality of time-lapse surveys, including sparse surveys, non-repeated surveys, and Ocean Bottom Node (OBN) vs streamer (different types of monitoring surveys); non-repeated sources, including source positioning errors and non-repeated source wavelets; spatial weighting operators in the L2-norm constraints; and sensitivity to weak time-lapse effects. These experiments show that HR-S-JMI is very robust to random noise, coherent noise, survey sparsity, survey non-repeatability, source positioning errors and source wavelet discrepancies. Furthermore, HR-S-JMI remains effective when the spatial weighting operators in the L2-norm constraints are largely relaxed and HR-S-JMI is capable of detecting weak time-lapse changes (e.g. velocity changes down to +/- 35 m/s). These features make it a suitable time-lapse processing solution for cost-effective (semi-)continuous monitoring, termed i4D survey technology, in which inexpensive localized and sparse surveys are employed between the conventional full-field surveys. The simultaneous strategy of S-JMI allows the full-field survey information to compensate the poor illumination of the in-between sparse surveys during process. Furthermore, calender-time constraints are proposed and applied to the parameter differences between the baseline and monitors along the calender-time axis by taking advantage of the feature that time-lapse effects usually develop gradually over time. With a complex synthetic example based on the Marmousi model, we demonstrate that S-JMI is a promising tool to process datasets acquired from (semi-)continuous monitoring, like an i4D survey.

In conclusion, we propose high-resolution simultaneous JMI (HR-S-JMI) as an effective time-lapse processing tool for the following main reasons:
• HR-S-JMI is able to make use of multiples to extend the illumination of the subsurface, instead of removing them;
• HR-S-JMI is an extended imaging process, including automatic velocity updating. Therefore, it takes velocity variations between surveys directly into account;
• HR-S-JMI is a good indicator of velocity changes, it can invert for high-resolution accurate time-lapse velocity changes;
• HR-S-JMI is robust to the uncertainties existing in the monitoring surveys, e.g. noise, sparsity, non-repeatability, source positioning errors, source wavelet discrepancy, etc;
• HR-S-JMI has the ability to detect weak time-lapse changes (velocity changes down to +/- 35 m/s)