Misfit Function Research Articles

An Effective Inverse Procedure for Identifying DEMParameters of Rock-Like Materials

This study aims to identify discrete element model parameters of rock-like materials. An inverse procedure is developed to determine the discrete element model parameters from experimental measurements. This involves the solution of an inverse problem through minimizing the misfit function which describes the error between numerical computation and experiment by an optimization procedure. In this procedure, the discrete element method is adopted as the numerical calculation method of the forward problem. The orthogonal experimental design is used for parameter sensitivity analysis. Besides, the approximation model with radial basis function is adopted instead of the actual calculation model to reduce the time of forward calculation. The ant-colony optimization algorithm is employed as the inverse operator. Therefore, the parameters of the discrete element model are optimized by this procedure. The three-point bending experiment with discrete element simulation is provided to verify the validity and accuracy of the inversion results. The results indicate that it can rapidly obtain the available and reliable model parameters just through a few sets of experimental data. As a result, this inverse procedure can be applied more widely to parameter identification of the discrete element model for brittle materials.

Amplitude-based misfit functions in viscoelastic full-waveform inversion applied to walk-away vertical seismic profile data

Viscoelastic full-waveform inversion is applied to walk-away vertical seismic profile data acquired at a producing heavy-oil field in Western Canada for the determination of subsurface velocity models (P-wave velocity [Formula: see text] and S-wave velocity [Formula: see text]) and attenuation models (P-wave quality factor [Formula: see text] and S-wave quality factor [Formula: see text]). To mitigate strong velocity-attenuation trade-offs, a two-stage approach is adopted. In Stage I, [Formula: see text] and [Formula: see text] models are first inverted using a standard waveform-difference (WD) misfit function. Following this, in Stage II, different amplitude-based misfit functions are used to estimate the [Formula: see text] and [Formula: see text] models. Compared to the traditional WD misfit function, the amplitude-based misfit functions exhibit stronger sensitivity to attenuation anomalies and appear to be able to invert [Formula: see text] and [Formula: see text] models more reliably in the presence of velocity errors. Overall, the root-mean-square amplitude-ratio and spectral amplitude-ratio misfit functions outperform other misfit function choices. In the final outputs of our inversion, significant drops in the [Formula: see text] to [Formula: see text] ratio (~1.6) and Poisson’s ratio (~0.23) are apparent within the Clearwater Formation (depth ~0.45–0.50 km) of the Mannville Group in the Western Canada Sedimentary Basin. Strong [Formula: see text] (~20) and [Formula: see text] (~15) anomalies are also evident in this zone. These observations provide information to help identify the target attenuative reservoir saturated with heavy-oil resources.

Wave-equation dispersion inversion of Love waves

We present a theory for wave-equation inversion of Love-wave dispersion curves, in which the misfit function is the sum of the squared differences between the wavenumbers along the predicted and observed dispersion curves. Similar to inversion of Rayleigh-wave dispersion curves, the complicated Love-wave arrivals in traces are skeletonized as simpler data, namely, the picked dispersion curves in the [Formula: see text] domain. Numerical solutions to the SH-wave equation and an iterative optimization method are then used to invert these dispersion curves for the S-wave velocity model. This procedure, denoted as wave-equation dispersion inversion of Love waves (LWD), does not require the assumption of a layered model or smooth velocity variations, and it is less prone to the cycle-skipping problems of full-waveform inversion. We demonstrate with synthetic and field data examples that LWD can accurately reconstruct the S-wave velocity distribution in a laterally heterogeneous medium. Compared with Rayleigh waves, inversion of the Love-wave dispersion curves empirically exhibits better convergence properties because they are completely insensitive to the P-velocity variations. In addition, Love-wave dispersion curves for our examples are simpler than those for Rayleigh waves, and they are easier to pick in our field data with a low signal-to-noise ratio.

Least-squares reverse time migration via linearized waveform inversion using a Wasserstein metric

Least-squares reverse time migration (LSRTM), an effective tool for imaging the structures of the earth from seismograms, can be characterized as a linearized waveform inversion problem. We have investigated the performance of three minimization functionals as the [Formula: see text] norm, the hybrid [Formula: see text] norm, and the Wasserstein metric ([Formula: see text] metric) for LSRTM. The [Formula: see text] metric used in this study is based on the dynamic formulation of transport problems, and a primal-dual hybrid gradient algorithm is introduced to efficiently compute the [Formula: see text] metric between two seismograms. One-dimensional signal analysis has demonstrated that the [Formula: see text] metric behaves like the [Formula: see text] norm for two amplitude-varied signals. Unlike the [Formula: see text] norm, the [Formula: see text] metric does not suffer from the differentiability issue for null residuals. Numerical examples of the application of three misfit functions to LSRTM on synthetic data have demonstrated that, compared to the [Formula: see text] norm, the hybrid [Formula: see text] norm and [Formula: see text] metric can accelerate LSRTM and are less sensitive to non-Gaussian noise. For the field data application, the [Formula: see text] metric produces the most reliable imaging results. The hybrid [Formula: see text] norm requires tedious trial-and-error tests for the judicious threshold parameter selection. Hence, the more automatic [Formula: see text] metric is recommended as a robust alternative to the customary [Formula: see text] norm for time-domain LSRTM.

Enhanced performance of ISC focal mechanism computations as a result of automatic first-motion polarity picking optimization

The International Seismological Centre (ISC) routinely calculates and makes available automatic earthquake focal mechanisms by combining reported parametric data (first-motion polarities) available in the reviewed ISC Bulletin and auto-picked first-motion polarities obtained from waveform data using a broadband automatic picker. In order to further enhance the robustness of the auto-picked polarities, we set up an optimization strategy which is carried out using the neighbourhood algorithm on a 24-processor mini computer cluster. The aim is to minimize an objective misfit function which takes into account the data uncertainties and compares the first P-wave arrival times and polarities of a large dataset of nearly 18,000 manual picks and the associated auto-picked waveform phase arrivals. The optimization yielded an overall increase of matching auto-picks from 15 to 30% in comparison with the default setup of the automatic picker. We then applied the optimized automatic picker to a set of earthquakes from the reviewed ISC Bulletin where we could not obtain well-constrained mechanism solutions using its default setup. As a result of using the optimized picker, we obtained well-trusted mechanism solutions for 28% of these cases by increasing the number of first motion auto-picked polarities, and hence minimizing the station azimuthal gap in some cases, and/or correcting some of the erroneous auto-picked polarities where possible.

Kantorovich-Rubinstein misfit for inverting gravity-gradient data by the level-set method

We have developed a novel Kantorovich-Rubinstein (KR) norm-based misfit function to measure the mismatch between gravity-gradient data for the inverse gradiometry problem. Under the assumption that an anomalous mass body has an unknown compact support with a prescribed constant value of density contrast, we implicitly parameterize the unknown mass body by a level-set function. Because the geometry of an underlying anomalous mass body may experience various changes during inversion in terms of level-set evolution, the classic least-squares ([Formula: see text]-norm-based) and the [Formula: see text]-norm-based misfit functions for governing the level-set evolution may potentially induce local minima if an initial guess of the level-set function is far from that of the target model. The KR norm from the optimal transport theory computes the data misfit by comparing the modeled data and the measured data in a global manner, leading to better resolution of the differences between the inverted model and the target model. Combining the KR norm with the level-set method yields a new effective methodology that is not only able to mitigate local minima but is also robust against random noise for the inverse gradiometry problem. Numerical experiments further demonstrate that the new KR norm-based misfit function is able to recover deep dipping flanks of SEG/EAGE salt models even at extremely low signal-to-noise ratios. The new methodology can be readily applied to gravity and magnetic data as well.

Notes on the meaning of adjoint wave equation derived via Lagrange multiplier approach

Full-wave adjoint inversion method can be formulated via the Lagrange multiplier approach, taking the forward seismic wave equation as an additional independent constraint on synthetic wavefield. The adjoint wave equation is explained as the Lagrange multiplier where the local minimum existed, in case the model parameters are unperturbed. Such an explanation is unsatisfactory, due to the fact that the variation of synthetic wavefield equals to zero in case the model parameters are unperturbed. To avoid such a contradiction, we give a more clear explanation about the meaning of the adjoint wave equation, that it serves for converting the gradient of misfit function with respect to synthetic wavefield to the gradient of misfit function with respect to the model parameter.

Misfit function for full waveform inversion based on the Wasserstein metric with dynamic formulation

Conventional full waveform inversion (FWI) using least square distance (L2 norm) between the observed and predicted seismograms suffers from local minima. Recently, the Wasserstein metric (W1 metric) has been introduced to FWI to compute the misfit between two seismograms. Instead of comparisons bin by bin, the W1 metric allows to compare signal intensities across different coordinates. This measure has great potential to account for time and space shifts of events within seismograms. However, there are two main challenges in application of the W1 metric to FWI. The first one is that the compared signals need to satisfy nonnegativity and mass conservation assumptions. The second one is that the computation of W1 metric between two seismograms is a computationally expensive problem. In this paper, a strategy is used to satisfy the two assumptions via decomposition and recombination of original seismic data. In addition, the computation of the W1 metric based on dynamic formulation is formulated as a convex optimization problem. A primal-dual hybrid gradient method with linesearch has been developed to solve this large-scale optimization problem on GPU device. The advantages of the new method are that it is easy to implement and has high computational efficiency. Compared to the L2 norm based FWI, the computation time of the proposed method will approximately increase by 11% in our case studies. A 1D time-shift signals case study has indicated that the W1 metric is more effective in capturing time shift and makes the misfit function more convex. Two applications to synthetic data using transmissive and reflective recording geometries have demonstrated the effectiveness of the W1 metric in mitigating cycle-skipping issues. We have also applied the proposed method to SEG 2014 benchmark data, which has further demonstrated that the W1 metric can mitigate local minima and provide reliable velocity estimations without using low frequency information in the recorded data.

Linking continental erosion to marine sediment transport and deposition: A new implicit and O(N) method for inverse analysis

The marine sedimentary record contains unique information about the history of erosion, uplift and climate of the adjacent continent. Inverting this record has been the purpose of many numerical studies. However, limited attention has been given to linking continental erosion to marine sediment transport and deposition in large-scale surface process evolution models. Here we present a new numerical method for marine sediment transport and deposition that is directly coupled to a landscape evolution algorithm solving for the continental fluvial and hillslope erosion equations using implicit and O(N) algorithms. The new method takes into account the sorting of grain sizes (e.g., silt and sand) in the marine domain using a non-linear multiple grain-size diffusion equation and assumes that the sediment flux exported from the continental domain is proportional to the bathymetric slope. Specific transport coefficients and compaction factors are assumed for the two different grain sizes to simulate the stratigraphic architecture. The resulting set of equations is solved using an efficient (O(N) and implicit) algorithm. It can thus be used to invert stratigraphic geometries using a Bayesian approach that requires a large number of simulations. This new method is used to invert the sedimentary geometry of a natural example, the Ogooué Delta (Gabon), over the last ∼5 Myr. The objective is to unravel the set of erosional histories of the adjacent continental domain compatible with the observed geometry of the offshore delta. For this, we use a Bayesian inversion scheme in which the misfit function is constructed by comparing four geometrical parameters between the natural and the simulated delta: the volume of sediments stored in the delta, the surface slope, the initial and the final shelf lengths. We find that the best-fit values of the transport coefficients for silt in the marine domain are in the range of 300−500 m2/yr, in agreement with previous studies on offshore diffusion. We also show that, in order to fit the sedimentary geometry, erosion rate on the continental domain must have increased by a factor of 6 to 8 since 5.3 Ma.

A graph space optimal transport distance as a generalization of Lp distances: application to a seismic imaging inverse problem

Optimal transport distance is an appealing tool to measure the discrepancy between datasets in the frame of inverse problems, for its ability to perform global comparisons and its convexity with respect to shifted patterns in the compared quantities. However, solving inverse problems might require to compare signed quantities, while the optimal transport theory has been developed for the comparison of probability measures. In this study we propose to circumvent this difficulty by applying optimal transport to the comparison of the graph of the data rather than the data itself. We investigate this approach in the frame of seismic imaging, where each channel of the oscillatory data is interpreted as a discrete point cloud. We demonstrate that the corresponding misfit function can be computed through the solution of series of linear sum assignment problem (LSAP), while, based on the adjoint state technique, its gradient can be computed from the assignment solution of these LSAP. We illustrate how this approach yields a convex misfit function in the frame of seismic imaging using the full waveform. We show how an efficient strategy, based on a specific LSAP solver, the auction algorithm, can be designed. We illustrate the interest of the approach on a realistic 2D visco-acoustic seismic imaging problem. The proposed strategy relaxes the constraint on the accuracy of the initial model, outperforming the conventional least-squares approach and a previously proposed optimal transport based approach. The computational time increases by only a few percents compared with the least-squares approach, opening the way to applications to 3D field data in the near future.

A multiscale approach to full-waveform inversion using a sequence of time-domain misfit functions

Most of the approaches designed to avoid cycle skipping in full-waveform inversion (FWI) involve calculating a sequence of inversions in a multiscale fashion. We have adopted an alternative strategy, which is inverting a sequence of different misfit functions in the time domain. This is an implicit multiscale approach in the sense that the used misfit functions are sensitive to different wavelengths, but all of the inversion steps use the same modeling algorithm and the same model grid. In the first and third inversion steps, the transmitted (early arrivals) and reflected (late arrivals) components of the wavefield envelopes are respectively fitted. The second step promotes a smooth transition between the first and third steps, by using the envelope of the complete waveform. Because fitting just the envelope of the reflected waves has a minor effect on the misfit function of the whole data set, the phases of the reflected waves are mostly fitted in the fourth step, which is based on the waveform misfit function preserving only the late arrivals. The third and fourth steps are of crucial importance to fit the reflected events. We test the sequential inversion approach with the Marmousi model using data sets with different frequencies, obtaining better estimates of the velocity field than those obtained with the classic FWI. The solutions obtained with classic FWI and sequential inversion approach degrade with a progressively higher peak frequency data set, but the classic FWI solution degrades more rapidly.

Adaptive traveltime inversion

We have developed a method to obtain a misfit function for robust waveform inversion. In this method, called adaptive traveltime inversion (ATI), a matching filter that matches predicted data to measured data is computed. If the velocity model is relatively accurate, the resulting matching filter is close to a Dirac delta function. Its traveltime shift, which characterizes the defocusing of the matching filter, is computed by minimization of the crosscorrelation between a penalty function such as [Formula: see text] and the matching filter. ATI is constructed by minimization of the least-squares errors of the calculated traveltime shift. Further analysis indicates that the resulting traveltime shift corresponds to a first-order moment, the mean value of the resulting matching filter distribution. We extend ATI to a more general misfit function formula by computing different order moment of the resulting matching filter distribution. Choosing the penalty function in adaptive waveform inversion (AWI) as [Formula: see text], the misfit function of AWI is the second-order moment, the variance of the resulting matching filter distribution with zero mean. Because our ATI method is based on a global comparison using deconvolution, such as AWI, it can resolve the “cycle skipping” issue. We evaluate our ATI misfit function and compare it with state-of-the-art options such as least-squares inversion (L2 norm), wave-equation traveltime inversion, and AWI using schematic examples before moving to more complex examples, such as the Marmousi model. For the Marmousi model, starting with a 1D [Formula: see text] model, with data without low frequencies (no energy below 3 Hz), a meaningful estimation of the P-wave velocity model is recovered. Our ATI misfit function (first-order moment) indicates comparable performance with the AWI misfit function (the second-order moment). We also include a real data example from the Gulf of Mexico to demonstrate the effectiveness of our method.

An Optimized Deep Network Representation of Multimutation Differential Evolution and its Application in Seismic Inversion

Seismic inversion problems are well-known to be nonlinear and their misfit functions often involve many local minima. Global optimization methods are capable of converging to the global minimum of a misfit function, thus, they are promising in seismic inversion. As a global optimization method, multimutation differential evolution (MMDE) has been proven to be effective in solving high-dimensional seismic inversion problems. However, it is challenging to choose the optimal parameters for MMDE to achieve the best performance in seismic inversion. In this paper, we propose a new deep network based on MMDE and name it as MMDE-Net, which enables us to learn the optimal parameters by using a network training procedure rather than empirically choosing them. Benefiting from the learned parameters, MMDE-Net has advantages over MMDE in applications. Numerical examples based on synthetic and field data set clearly indicate that MMDE-Net can provide faster convergence speed and better inversion result than conventional methods in seismic inversion.

Joint Multiscale Direct Envelope Inversion of Phase and Amplitude in the Time–Frequency Domain

Time–frequency analysis can reveal local variations and allow for separation of phase and amplitude information of nonstationary seismic waveforms. Seismic signals are used since long as a robust tool for inversion of underground structures, as has been the practice in geophysical exploration. However, the mixing of phase and amplitude in seismic data increases the nonlinearity of seismic inversion. The authors first use Gabor transform to separate the phase and amplitude information of envelope data, and then introduce an adaptive factor into the misfit function to redistribute the weight of phase and amplitude information for direct envelope inversion (DEI) in the time–frequency domain. By adopting this procedure, greater flexibility can be achieved in operating the local phase of envelope and waveform spectra to enhance stability of multiscale phase inversion. For DEI, the direct envelope Frechet derivative is used, and thus, no weak scattering assumption is imposed on the joint multiscale DEI of phase and amplitude (PADEI). Compared with the DEI method, the PADEI can better recover the deeper parts of salt-bottom and subsalt structures by boosting the signal energy and weakening the nonlinearity of the waveform inversion.

Acoustic and geoacoustic inverse problems in randomly perturbed shallow-water environments.

The main goal of this paper is to estimate the regional acoustic and geoacoustic shallow-water environment from data collected by a vertical hydrophone array and transmitted by distant time-harmonic point sources. The aim is to estimate the statistical properties of the random fluctuations of the index of refraction in the water column and the characteristics of the sea bottom. It is explained from first principles how acoustic wave propagation can be expressed as Markovian dynamics for the complex mode amplitudes of the sound pressure. This makes it possible to express the cross moments of the sound pressure in terms of the parameters to be estimated. Then it is shown how the estimation problem can be formulated as a nonlinear inverse problem using this formulation, which can be solved by minimization of a misfit function. The method is applied to experimental data collected by the Acoustic Laboratory for Marine Applications system. A Bayesian analysis quantifies the uncertainty of the estimation.

Frequency‐domain double‐plane‐wave least‐squares reverse time migration

ABSTRACTLeast‐squares reverse time migration is often formulated as an iterative updating process, where estimating the gradient of the misfit function is necessary. Traditional time domain shot‐profile least‐squares reverse time migration is computationally expensive because computing the gradient involves solving the two‐way wave equation several times in every iteration. To reduce the computational cost of least‐squares reverse time migration, we propose a double‐plane‐wave least‐squares reverse time migration method based on a misfit function for frequency‐domain double‐plane‐wave data. In double‐plane‐wave least‐squares reverse time migration, the gradient is computed by multiplying frequency‐domain plane‐wave Green's functions with the corresponding double‐plane‐wave data residual. Because the number of plane‐wave Green's functions used for migration is relatively small, they can be pre‐computed and stored in a computer's discs or memory. We can use the pre‐computed plane‐wave Green's functions to obtain the gradient without solving the two‐way wave equation in each iteration. Therefore, the migration efficiency is significantly improved. In addition, we study the effects of using sparse frequency sampling and sparse plane‐wave sampling on the proposed method. We can achieve images with correct reflector amplitudes and reasonable resolution using relatively sparse frequency sampling and plane‐wave sampling, which are larger than that determined by the Nyquist theorem. The well‐known wrap‐around artefacts and linear artefacts generated due to under‐sampling frequency and plane wave can be suppressed during iterations in cases where the sampling rates are not excessively large. Moreover, implementing the proposed method with sparse frequency sampling and sparse plane‐wave sampling further improves the computational efficiency. We test the proposed double‐plane‐wave least‐squares reverse time migration on synthetic models to show the practicality of the method.

Joint inversion methods with relative density offset correction for muon tomography and gravity data, with application to volcano imaging

Muography is a relatively new geophysical imaging method that uses muons to provide estimates of average densities along particular lines of sight. Muography can only see above the horizontal elevation of the detector and it is therefore attractive to attempt a joint inversion of muography data with gravity data, which is also responsive to density but generally requires combination with another geophysical data set to overcome issues related to non-uniqueness and poor depth resolution. Some previous work has investigated this joint inverse problem and demonstrated the potential improvements to be gained by jointly inverting muography and gravity data. However, there has yet to be a thorough investigation of different numerical approaches for formulating the joint inverse problem. Particularly important is how to account for the fact that even though the two data types are sensitive to the same physical quantity, density, they respond through different response functions. Moreover, the two measurements are affected by different systematic uncertainties that are difficult to model. In this work, we considered an approximation where the two density quantities, inferred from the two data types, can be related by an unknown scalar offset. We considered various existing and new joint inversion methods that might solve this problem and we applied them to a synthetic volcano imaging scenario based on the Puy de Dome volcano in the Central Massif region of France. We used unstructured meshes in our modelling to adequately honour the significant topography in that scenario. Our experiments indicated that the most successful joint inversion method for this type of geological scenario was one in which the data misfit function is reformulated to automatically determine the best-fitting offset following a least-squares minimization argument. However, other approaches showed merit and we suggest several of the investigated methods be applied and compared for any specific joint inversion scenario.

Iterative multiscale gradient computation for heterogeneous subsurface flow

We introduce a semi-analytical iterative multiscale derivative computation methodology that allows for error control and reduction to any desired accuracy, up to fine-scale precision. The model responses are computed by the multiscale forward simulation of flow in heterogeneous porous media. The derivative computation method is based on the augmentation of the model equation and state vectors with the smoothing stage defined by the iterative multiscale method. In the formulation, we avoid additional complexity involved in computing partial derivatives associated to the smoothing step. We account for it as an approximate derivative computation stage. The numerical experiments illustrate how the newly introduced derivative method computes misfit objective function gradients that converge to fine-scale one as the iterative multiscale residual converges. The robustness of the methodology is investigated for test cases with high contrast permeability fields. The iterative multiscale gradient method casts a promising approach, with minimal accuracy-efficiency tradeoff, for large-scale heterogeneous porous media optimization problems.

A robust source-independent misfit function for time domain waveform inversion based on normalized convolved wavefield

It is usually necessary to provide a source wavelet before the waveform inversion, and accuracy of the wavelet determines whether the synthetic data matches the observed record. In addition to directly estimating a wavelet, another method is to construct a source-independent objective function. The source difference disappears in the convolved wavefield waveform inversion through the convolution of one data set and the reference trace selected from the other data set, because two different wavelets are already present in the two convolution terms. However, the convolution wavefield inversion is also seriously affected by the incorrect amplitude information of the seismic waves. The attenuation of seismic waves in underground media is quite common in practice. In contrast, the method based on the global normalized wavefield is less sensitive to the amplitude loss because it compares the phase information of the two data sets. The approach we present in this paper is a combination of two methods, that is to construct the normalized wavefield on the basis of the convolution wavefield. In this way, we alleviate the problem of the traditional convolution wavefield inversion caused by amplitude mismatch, and further extend the applicability of this method. The new normalized convolution wavefield inversion combines the advantages of the two methods, which are independence on the source wavelet and being insensitive to amplitude loss, we give a multi-scale and multi-stage inversion strategy for the new objective function, which can improve the implementation of the multi-scale strategy and effectively deal with missing low frequencies.

Tikhonov-type regularization and the finite element method applied to point source estimation in the atmosphere

This work is concerned with the source identification problem in air pollution modeling. The forward problem is described by an advection-diffusion equation with physically realistic coefficients. It is solved by the combination of adaptive meshes and a stabilized finite element method. The source is estimated by Tikhonov-type regularization with a composite misfit function and an entropic penalty term. The composite misfit is given by the convex combination of ℓ1 and ℓ2 norms of the discrepancy between observed and predicted concentrations. Morozov's discrepancy principle is applied to simultaneously choose the regularization and the misfit combination parameters. Numerical experiments using Copenhagen field data is used to validate numerically the proposed methodology.