Monitoring time-lapse changes inside the subsurface is of great significance to many geotechnical applications, such as storage of gasses in underground geological formations. Minute differences in the seismic wavefield between an initial baseline and a subsequent monitor survey have to be detected in order to observe fluid flow inside subsurface reservoirs. This problem becomes even more challenging when the reservoir is situated underneath a series of complex, highly reflective layers. Such an overburden will generate strong multiple reflections that will interfere with the reflections of the target zone. Ideally, a methodology is designed in order to remove these internal multiples to allow a clear view of the reservoir response for time-lapse analysis. The Marchenko method can redatumthe seismic wavefield to arbitrary depth levels or points in the subsurface, while accounting for all orders of internal multiple reflections. This method, therefore, has great potential to solve some of the time-lapse issues, as it is able to closely examine specific zones of interest in the subsurface without distortions from surrounding layers. Time-lapse studies are often hampered by irregular or imperfect sampling, whereas the Marchenko method relies on densely sampled, co-located sources and receivers. It is, therefore, important that the Marchenko method is able to handle more complex acquisition geometries. This can either be achieved by interpolating the reflection data as a pre-processing step or by correcting for errors inside the Marchenko scheme. Here, point-spread functions are introduced that describe the imperfections in the reflection data. These imperfections distort the focusing and Green’s functions retrieved from the Marchenko method. Next, each iteration of theMarchenko scheme is extended to deblur the imperfect focusing and Green’s functions by multidimensional deconvolution with these point-spread functions. Additionally, a slight modification is required to ensure stability of the new scheme. This new iterative Marchenko scheme is computationally more expensive, but removes all sampling artifacts. Finally, the migrated images of the target zone show significant improvements, when using either the new scheme or interpolation as pre-processing step...

Reflection seismology aims to estimate the Earth's subsurface elastic parameters for further investigation by geologists and engineers. This involves generating elastic waves using seismic sources and recording the Earth's response with receivers. The subsurface model is typically considered a combination of a background model and a short-wavelength reflectivity model. There are two main paths to estimate these parameters: non-linear waveform inversion to directly compute the elastic parameters or depth migration to estimate a structural image or reflectivity of the subsurface.

Reverse-Time Migration (RTM) is a common depth migration technique that migrates recorded wavefields from the space-time domain to the space-depth domain. It utilizes the Born approximation and the adjoint of the Born operator to produce an RTM image. However, RTM can suffer from errors, such as noise, temporal and spatial limitations, and multiple reflections.

Least-Squares Reverse-Time Migration (LSRTM) is used to overcome some of these errors. LSRTM involves resolving the reflectivity model by least-squares inversion, which is computationally expensive. Gradient-based optimization algorithms are often employed to reduce the computational burden, but they still require solving the wave equation and its adjoint for a large model in multiple iterations. One way to reduce the computational cost is by limiting the computational domain to a target region of interest.

Target-oriented LSRTM, known as TOLSRTM, focuses on the wavefield just above the target by bypassing the overburden. This approach proves beneficial when the overburden generates strong internal multiple reflections that obscure the reflections from the target area. However, a redatuming method is required to predict all orders of multiples. Marchenko redatuming is a data-driven technique that predicts the Green's functions at the boundary of the target region, incorporating all orders of internal multiples. It allows for double-sided redatuming, considering both the source and receiver perspectives. By combining the LSRTM algorithm and Marchenko double-focusing, a target-oriented LSRTM algorithm is devised that can predict interactions between the target and overburden and remove the effects of the overburden in the image. Predicting these interactions results in an artifact-free image, a better convergence rate, and a high-resolution image of the target.

Target-oriented migration algorithms typically consider only the upper horizontal boundary of the region of interest (ROI), neglecting wavefields entering the ROI from the medium beneath the lower boundary. To address this, a target-enclosed LSRTM algorithm is proposed, including both the ROI's upper and lower boundaries. Including the lower boundary provides transmission information and can improve inversion convergence. In addition, this algorithm is adopted for virtual receivers created by Marchenko redatuming. In the case of physical receivers at the boundaries of the target zone, the target-enclosed algorithm can incorporate the transmission information emanating from the lower boundary to the upper one. Consequently, when the initial model is far from the actual model, the resulting image partly recovers the long wavelength part of the model in agreement with the Born approximation criteria. Moreover, when an initial model closer to the actual model is used, the algorithm can partially recover the vertical interfaces of the perturbation. In the case of virtual receivers at the boundaries of the target zone, since the Marchenko redatuming is performed in the initial background model, the redatumed wavefields at the lower boundary suffer from kinematic errors. Therefore, the algorithm can not recover the long wavelength part of the model.

The thesis concludes with a discussion of the results obtained from applying the algorithms to marine datasets. The images resulting from the Marchenko double-focusing based target-oriented LSRTM algorithm show improvements in both resolution and artifact reduction by suppressing the overburden generated internal multiple effects. Moreover, the double-focusing enables the user to reduce the computational costs of the LSRTM algorithm and choose finer spatial sampling for the image.

An appendix proposes a formulation for integrating the target-oriented algorithms with non-linear inversion like Full Waveform Inversion (FWI). The results of this proposed algorithm show its effectiveness by reducing the internal multiple related artifacts and increasing resolution and faster convergence.

Motivated by environmental concern, the industry has been developing an alternative marine seismic source, in particular a marine vibrator. By spreading the emitted energy out over time, vibrator sources are perceived to be less intrusive to marine mammals. It is also believed that vibrators have greater control of the emitted source wavelet than can be achieved with traditional airguns. With the added control, it is possible to only emit portions of the frequency spectrum, which in turn allows for many applications such as deblending and the ability to avoid masking mammal communications. To effectively implement these, two methodologies are proposed to interpolate the frequency data that are not emitted. The first is a deep learning approach utilizing a U-Net architecture, with a custom frequency loss function. The second is a sparse optimization method that approximates the reflectivity series of the subsurface using known frequency content. By assuming that the signal can be represented sparsely and that all frequencies interact with the subsurface interfaces similarly at all frequencies, the frequency spectrum can be reconstructed. Both of the presented methods are tasked to interpolate the missing frequency band(s) in North Sea shot data. It is found that both methods are able to interpolate narrow 2.5 Hz bands, but are unable to accurately reconstruct wider (ex. 10 Hz), frequency bands. Overall, the U-Net shows better results than the sparse optimization method when the frequency gaps are positioned closely.

Damages within asphalt have been interesting phenomena in asphalt engineering, the detection of which is significant for maintenance of road sections. This project focuses on cracks and delaminations. An attempt was made to filter radar image data with a method based on a VNA-antenna-multilayered system model as well as the data from two specific measurements, aiming at better visualizing cracks as well as other features in radar image, and the results were checked and analysed. This part of work has provided an application of the aforementioned method of radar image data filtering as well as the points worth noticing and avoiding when making this application.
For delaminations, machine learning algorithms, first the EM algorithm and then the YOLO v3 algorithm, were used as an attempt to highlight and detect them. Though the results still need improving, it is still valuable that the workload for human intervention can be alleviated with the help of these algorithms and that better performance can be expected based on current work, with the increasing amount of data with high quality achieved in the future.

In this thesis I conducted ground penetrating radar (GPR) and ground conductivity meter (GCM) surveys to detect the presence of simulated clandestine burials at the Amsterdam Research Initiative for Subsurface Taphonomy and Anthropology (ARISTA) test facility, and determine their characteristic response in this environment; providing valuable insights and recommendations for forensic investigations. I performed four days of GPR and GCM surveys over three simulated clandestine burials at ARISTA. I collected common-offset GPR data to investigate changes to burial detectability due to different central antenna frequencies (250 MHz and 500 MHz), different GPR instruments (NOGGIN or pulseEKKO), changes to survey grid orientation relative to burials, and increased soil moisture content in the survey area. Additionally, I acquired common-source GPR data to examine the efficacy of electromagnetic interferometry (EI) and adaptive subtraction (AS) methods in improving burial detectability. I conducted GCM surveys with two coil configurations, (horizontal co-planar (HCP) and vertical co-planar (VCP)), three intercoil spacings (0.32 m, 0.71 m, 1.18 m), two different line spacings (0.5 m, 0.25 m) and in the presence of variable soil moisture content. I also performed low induction number (LIN) correction and elevation correction procedures on GCM data to determine the extent to which these influence the detectability of clandestine burials in this environment.
In common-offset radargrams characteristic burial anomalies take on many forms, appearing as disruptions to existing features (direct-wave arrivals and soil horizons) and as isolated reflection events (hyperbolic events and burial length horizontal anomalies). In timeslices, burials are characterized by high or low amplitude rectangular anomalies. When used in conjunction, radargrams and timeslices produced characteristic responses regardless of survey grid orientation, consistent with the locations of the burials. Increased soil moisture at the site improved the detectability of burials and the 250 MHz antenna was found to be superior to the 500 MHz antenna in obtaining a characteristic burial response, though both were successful to a large extent. EI and AS processing techniques were successful in removing direct-wave contributions in radargrams, though detectability was not significantly improved when compared to raw data. Overall, the three burials were detected using GPR to various extents, and in future work thorough historical data in addition to zero-measurements should be obtained for all burials in order to investigate the source of these differences. GCM surveys conducted in this work were largely unsuccessful in detecting simulated clandestine burials due to significant conductive noise sources (metal fence, sensors, etc.) and the limited conductivity contrast in the soil. Low conductivity zones were detected over some burials using HCP at an intercoil spacing of 1.18 m, however, confidence in the validity of these responses is low due to the dominating noise sources.

As an alternative to the industry standard processing of total field data, a Basin-hopping inversion can be used. Another alternative is the use of vector field magnetometers, which can be processed with a Gauss-Newton inversion. Both the Basin-hopping and Guass-Newton inversion consist of a two step inversion. First, an inversion to estimate the centre location and magnetic dipole moment of a point source is carried out, the so called point-source inversion. Second, the retrieved magnetic dipole moment can be used to reconstruct the magnetic moment with a prolate spheroid model, which is called the prolate-spheroid inversion.

First, both the TMI and vector data inversion are tested on 3000 noise-free synthetic datasets located in the North Sea where a prolate spheroid with varying shape, orientation and magnetic susceptibility creates an anomaly field. The largest global root mean square (GNRMS) error for the point-source inversion of the Gauss-Newton vector data inversion has a value in the order 10⁻⁹ whereas the Basin-hopping inversion for the TMI data returns for only 61% of the cases a GNRMS error in the order 10⁻⁷ or lower.
The result of the point-source inversion serves as a plateau for the prolate spheroid inversion; when the point-source inversion does not reach the set level of accuracy, it cannot be determined if the prolate spheroid model describes the true object well. For both the TMI and vector data, the UXO was correctly classified in 72% of the 3000 datasets. 
Second, comparing the current industry standard to the Basin-hopping point-source inversion only gave inconclusive results. The estimated centre location of both methods are within a 0.15 m radius but the disagreement in magnetic moment leaves an advice on which technique would work better indecipherable.

Failure of subsea wind turbine cables are the main cause for wind farm downtimes. Furthermore, 80% of insurance payouts to wind farms come from cable repairs and maintenance. Cable tracking is part of the wind farm cables maintenance scheme. Surveys are required for localisation and determining the depth of burial of the cables. Conventional surveying methods have shortcomings when it comes to localising offshore buried power cables, which can be resolved by using electromagnetic methods. From the measured magnetic field created by the cable, its position and burial depth can be inverted. Forward modelling the magnetic response of cables improves the understanding of their electromagnetic behaviour in one dimensional complex submarine environments. Maxwell's equations are the cornerstones of understanding the behaviour of electromagnetic fields, forming the basis of the models. Implementing numerical modelling packages like Empymod accommodate for boundary conditions in one-dimensional layered media to more accurately investigate the magnetic response of electric cables. Land measurements show that Empymod's numerical models can model cables in the accuracy range of 2-15% depending on the distance. New numerical model analyses show the influence of water depth, underlying layers, depth of burial, geometrical complexities and adjacent cables on the magnetic field. Surrounding resistivity contrasts together with changes in water depth change the measured magnitude of the magnetic field by up to 10%. Geometrical complexities have large effects on the magnetic field, sometimes completely changing the expected response. Therefore, twisted and irregularly shaped cables can pose problems when inverting for their position during surveys.

The magnetic method is one of the geophysical techniques that prevails in the detection and identification of Uneploaded Ordnance (UXO). Inversion of the magnetic data, allows us to recover the position and the magnetic dipole moment of the object. The typical use of a model is a sphere, though it does not
provide information about the shape of the object. This thesis focuses on modeling a prolate spheroid, that not only recovers the position and the magnetic moment moment, but also the orientation and the dimensions of body. The results were fairly similar, though the residuals in the case of the sphere were less than the ones of the spheroid. In addition, measurements of real UXOS were used to verify the code. The inversion yielded reasonable models in four cases, but did not manage to recover an accurate model for one them. Lastly, the octupole moment was modeled as well, though its contribution to the results was deemed negligible.

Clay has been used for many hundreds of years in building dikes. Their physical properties and thickness are essential for a dike to maintain its function. For the last fifteen years electromagnetic induction instruments (EMI) have been playing a growing role in mapping electromagnetic soil properties. Layered models with sharp boundaries between different soil units represents the soil profile accurately. Therefore, two database inversion procedures are developed and proposed. Synthetic data sets and field measurements have been used to validate and compare our retrieved electrical conductivity/resistivity models. Field measurements of a so-called multi-receiver EMI instruement (CMD-Explorer and CMD-MiniExplorer) have been used to perform the database inversion on, and are being compared to ERT data and hand drill data. The EMI instruments (CMD-Explorer and CMD-MiniExplorer) consisting of three receiving antennas, fixed operating frequency and two measurement configurations will not allow us more information than of a estimate three-layered model. Therefore, the calibration process here proposed is based on finding the best fit of the electromagnetic responses (in the least squares sense) in a 5-D solution space created by calculating the electromagnetic responses using five medium parameters and the layered earth response. These medium parameters are the three layered conductivity’s and thicknesses of the first two layers. Using these five medium parameters, a database is built based on predefined ranges of medium parameters and is being used to perform the inversion procedure with. Two methods of inversion have been proposed: constrained and full database inversion. The synthetic models have been recovered with an average data misfit of 0.0096 and 0.0032 for the constrained data base inversion for respectively model 1 and model 2. For the full database inversion data misfit values of 0.00316 and 0.0021 were accomplished. Comparing ERT and electrical resistivity obtained from EMI measurements, one could see that the trend and first two layer thicknesses were recovered. Verifying EMI measurements using hand drill data, we were able to recover the top layer at few locations.

The hazards of unexploded ordnance threaten the increasing marine construction activities nowadays, which increases the importance of unexploded ordnance detection. Research has shown that transient electromagnetic methods can successfully be used to detect unexploded ordnance on land. New equipment is being developed to make marine unexploded ordnance detection also possible. This study aims to determine which targets can be detected and which not in a marine environment through a numerical feasibility study. Building on an existing geophysical simulation framework, it asks: Under which conditions can we detect a conductor on or below the seafloor using a time domain loop source? Through the three-dimensional modelling of Maxwell’s equations, responses were computed for hollow rectangular targets of different burial depths, sizes, wall thicknesses, andwall conductivities. For the analysis of these responses two quantities were introduced, a net effect and a measurability. Evaluation of these quantities demonstrated the individual impact of the tested parameters on these quantities as well as the relative significance of the influence of these parameters. The results included derived relations for the influence of individual parameters on the net effect, as well as limits on the measurability of targets. Arectangular conductor of 0.1 by 0.1 by 0.4 metres or smaller with a wall thickness of 10 millimetre, buried more than 2 metres under the seafloor is not measurable under the noise assumptions made. The relative significance of the parameters was found to be from most to least significant: burial depth, size, wall conductivity, and wall thickness.

Geophysical monitoring is a popular tool in aquifer characterization and groundwater flow. To address this objective at groundwater extraction site ‘t Klooster, an ERT dataset was analyzed to identify groundwater flow patterns resulting from the injection of warm oxygenated water. Using a petrophysical model, changes in resistivity were converted to estimated temperature changes to visualize the spread of warm oxygenated water. Multi-dimensional analysis of the resistivity response of the subsurface was carried out. This allowed for the division of the subsurface into 4 depth regimes according to their response to well activity. It is shown that wells up to 100m removed from the ERT set-up influenced the temperature distribution. Furthermore, injected oxygenated water highlighted a preferential flow path between the depths of 20 and 35m in a north-west direction. This is in line with global groundwater flow in the area. Groundwater flow effects could not be reliably separated from the effect of well activity, however its effect is recognized both during extraction of groundwater and injection of warm water.

The aim of this study is to look into the effect of water content in clays on the electrical conductivity as obtained from measurements of an frequency domain electromagnetic induction survey and a laboratory resistivity test. The laboratory test results are compared to theoretical conductivity models. The focus is hereby on determining the relationship between decreasing water content in clay and the electrical conductivity. The research has been performed in collaboration with the land geophysics department at Fugro. I would like to thank Veerle Steenhuisen for giving me the opportunity to carry out this research at Fugro and Serkan Elgun for supervising my thesis on behalf of Fugro. I also want to thank Evert Slob and Dominique Ngan-Tillard for the help on academic level and supervising my thesis.

The position of installed submarine power cables is often not accurately known for several reasons. The precise knowledge of the cable position is important for the maintenance process and necessitates the need for cable tracking systems. Current systems are in many cases imprecise and have a short sensing distance, limited to a few meters. A device with better accuracy and sensing range is in demand for enhancing the cable maintenance process. Three new inversion algorithms are introduced which invert the passive electromagnetic _eld created from an injected signal in the target cable. The cable is treated like a line source. This can be seen as an inverse source location _nding problem, while the source time signature is known. The algorithms are tested on synthetic data, modeling a power cable in homogeneous sea water with noise. Several parameters like sensor array, dip angle of the cable and relative position to the system are analyzed. The inuence of soil and the sea surface are studied on synthetic data created from a numerical three-layer forward model. Additionally, a prototype is developed, and di_erent processing schemes are presented and compared. The system is tested with di_erent inversion algorithms on a _eld cable on land. A solution is found which can determine the cable position accurately in sea water with noise terms for a distance up to 6 m. However, this system is highly sensor array dependent. A second inversion method gives a smaller sensing range of 5m but can be used with more versatility. In a di_erent modeled scenario without noise terms, the cable is buried in marine sediments and air is on top of a layer of sea water. The sensors are in the sea water layer. The soil has a great impact on the accuracy of the inversion. In some cable-system-orientations, the absolute error caused by the layered earth exceeds 1 m. The prototype is successfully tested on a _eld cable for di_erent scenarios. An inverted cable position in a global reference frame can be obtained by including a motion sensor.

The multi-pole PML (MPML) is tested on models that simulate seismic waves traveling through the subsurface. Using a recursive integration technique a stretching function consisting of the sum of multiple stretching functions is implemented in the velocity-stress finite difference time domain wave equations. The MPML is implemented in both the rotated staggered grid (RSG) and the Virieux grid. The performance of the MPML is tested on a square model, rectangular model and a rectangular model with a free-surface and compared to other types of PML’s implemented in these models. The main result is that the MPML can be implemented in the velocity stress wave equations giving stable results similar to other PML types.

Forensic investigations focused on determining clandestine buried weapons, narcotics or even homicide evidence, can be expensive, inefficient and depend greatly on prior information and the tools available. Geophysical tools have potential to improve these investigations, on the ground that they can detect shallow buried objects in a non-invasive way. This report studies the use of electrical resistivity tomography (ERT) in combination with ground penetrating radar (GPR) to detect buried objects at two different sites, both with conditions found here in The Netherlands. The first is a test site at the Technical University of Delft, The Netherlands. Here, two plastic barrels are buried, one is empty and the other is filled with metal rods. The second site is at the ARISTA Facility in Amsterdam, The Netherlands. Here, human cadavers are buried for forensic research. The ERT was used in two ways, the first was to produce 2 dimensional (2D) surveys, which were combined to create a 3 dimensional (3D) model. The second method made measurements that created a 3D model directly. For each the dipole-dipole array was used. At the TU Delft site the results for the 2D method showed clear resistivity anomalies at the locations of the barrels. These anomalies corresponded to clear reflections in the GPR radargrams. The results at the ARISTA facility are inconclusive due to damage in the instrument used. The grid designs made for both ERT methods could however be used in a continuation of this study, and future research should be done within this topic to improve forensic investigations.

Backward erosion piping is a dike failure mechanism. It is the internal erosion process by which sand is eroded away from underneath a dike or levee by seepage flow. This erosion process progresses in the direction opposite to the direction of seepage flow and forms a small pipe directly beneath the dike. As erosion continues, this process can lead to dike failure. During this erosion process, the groundwater flow pattern is subject to continuous change, due to the growth of the pipe. Self-potential (SP) monitoring is sensitive to changes in the groundwater flow pattern, because of the electrokinetic coupling between fluid flow and the streaming potential. The SP field due to flow underneath a test dike was modeled with a FreeCAD -> Gmsh -> pyGIMLi workflow. This workflow can also be used to effectively resolve a wide range of standard and customizable geophysical modeling and inversion tasks. After modeling, field experiments were conducted, on the same test dike, to further assess the possibilities and limitations of SP monitoring to track the progress of backward erosion piping. Given that it is essential to have an accurate resistivity model in order to find the location of the SP source, an integrated electrical resistivity tomography (ERT) and SP monitoring system was designed. The electrodes used in this monitoring system were polarizable stainless steel stakes. The reliability of polarizable electrodes was greatly overestimated, as they turned out to give unstable SP measurements. The reason for the inferior reliability of polarizable compared to non-polarizable electrodes, was found through extensive literature research. The reason being that the largest potential in any electrode originates from the contact between the metal and the electrolytes in solution. The metal of polarizable electrodes is in direct contact with the electrolytes in the soil, which have variable concentrations. Therefore, the potential measured fluctuates together with the concentration of soil electrolytes in contact with the metal. The metal of non-polarizable electrodes, on the other hand, is in contact with a solution of its own salt, which has a constant concentration. Finally, piping is not expected to be measurable with SP monitoring, before a large sand boil is visible in the field. Once a positive SP anomaly develops at the sand boil, changes in the SP field due to the growth of the pipe are expected to be too small relative to the SP anomaly associated with water flow through the sand boil. Even though an integrated ERT and SP monitoring system is known to provide useful information about the hydrology of a dike, such a system is not sensitive enough to be able to monitor the development of a backward erosion pipe.

One of the most crucial estimates retrieved from measured seismic reflection data is the subsurface image. The image provides detailed information of the subsurface of the Earth. Seismic reflection data consists of so-called primary and multiple reflections. Primary reflections are events that have been reflected a single time, while multiple reflections have been reflected multiple times before they are recorded by the receivers. Most current migration algorithms assume all reflections in the data are primary reflections. Hence, in order to retrieve an accurate image of the subsurface, multiple reflections need to be eliminated before migration. Keeping the multiple reflections in the measured seismic reflection data will lead to a sub-optimal image of the subsurface, because the multiple reflections will be imaged as if they were primary reflections. Such artefacts in the image can cause erroneous interpretation...

Marchenko inversion is a new way to invert seismic or electromagnetic data recorded during geophysical surveys. The inversion method uses Marchenko theory. This is a recent development which enables the retrieval of Green's functions at any place in the subsurface. A non-recursive Marchenko inversion method has already been introduced but in this thesis a recursive Marchenko inversion method is implemented and analysed. A recursive scheme lies at the center of this new method. In this thesis, the new method is implemented and tested on a 1D subsurface model. The recursive scheme is first validated. This is done by computing a reflection response with it and comparing it with a reflection response resulting from forward modeling. After this, the accuracy of retrieved local reflection coefficients from the recursive inversion method is determined. This is done by comparison with exact reflection coefficients of the subsurface model. After this, several different parameters of the used subsurface model, data computation and the recursive inversion method itself are investigated for their influence on the accuracy of the inversion method. In particular interest is the effect of interval time errors because these result in errors that can build up rapidly through the recursion. However, the method has a big advantage. It is shown that the recursive Marchenko inversion method has a way to retrieve the magnitude of made interval time errors and correct for these when interval times are overestimated. In this way the error build up is stopped. In the end, it is shown that the new method delivers high accuracy results and has an advantage in computational expense compared to the existing recursive Marchenko inversion method. It is concluded that the new method shows promising prospects and that it is worthwhile to investigate the method further.

The Marchenko method is being developed to compute redatuming operators from the reflection response based on a new normalization option at the data level. These operators can be used to retrieve the two-way focusing functions in terms of pressure and velocity, respectively, without effects from the ghost, and free-surface and internal multiples. To compute these operators an estimate of the first arrival of the redatuming operator is required. This can be done using the same model as is used for traditional migration. We then propose to delay the need for such a model and remove free-surface and internal multiples in one step while retaining the primary reflections at the original two-way traveltime. We modify the amplitude by removing the transmission effects at both sides of the equations. The resulting data becomes full focusing functions in terms of pressure and velocity at two-way traveltime, respectively. New datasets can then be generated by picking the value at each focusing level and storing it which are free of the ghost effects, free surface and internal multiples.

The finite-element method can easily handle complicated geometries because of the application of unstructured meshes. Unlike the Cartesian grid used in the finite-difference method, the unstructured mesh can follow the sharp interfaces that separate two layers of different properties. Therefore, the finite-element method can provide more accurate solutions for the simulation of seismic wave propagation. Meshes of good quality are required for the finite-element simulation. However, it is not trivial to set up an appropriate mesh. First
of all, the mesh should contain elements of good shapes and sizes. In addition, the sharp interfaces should coincide with the edges of the elements instead of intersecting with them. These requirements are formulated as an optimization problem with three terms, measuring the difference between the actual and prescribed scaling field, shape quality, and the area between prescribed curves and the nearest triangle edges. The solution of the optimization problem should provide the desired mesh. The mesh generator MESH2D was applied to obtain an initial mesh. The Matlab function minFunc was used to search for the minimum of the constructed objective function. Three weights balance the three terms in the objective function. When it comes to complicated models, these weights have to be chosen carefully to produce a reasonable mesh.