Developed Inversion Algorithm Research Articles

Three-Dimensional Inversion of Long-Offset Transient Electromagnetic Method over Topography

The long-offset transient electromagnetic method (LOTEM) is widely employed in geophysical exploration, including environmental investigation, mineral exploration, and geothermal resource exploration. However, most interpretations of LOTEM data assume a flat Earth, and the commonly used one-dimensional (1D) interpretation encounters significant challenges in achieving reliable geological interpretations when topography is ignored. To address these challenges, this study presents an effective three-dimensional (3D) LOTEM inversion method. In this study, we discretize the simulation domain using unstructured tetrahedra to accurately simulate complex geological structures. The finite-element time-domain (FETD) method is utilized to calculate the LOTEM forward responses, and the limited-memory BFGS (L-BFGS) optimization method is employed for 3D LOTEM inversion. To avoid explicit calculation of sensitivity, we obtain the product of the transposed sensitivity matrix and the vector through adjoint forward modeling. Several synthetic models are used to verify the developed program, and the influence of topography on LOTEM inversion is examined. The numerical results demonstrate that topography can significantly impact the inversion result, potentially leading to incorrect geological interpretations. Finally, the developed inversion algorithm is applied to a realistic ore model from Voisey’s Bay, Labrador, Canada. The 3D inversion successfully reconstructs the spatial distribution of the ore body, further confirming the effectiveness of the developed algorithm.

3D Bayesian Inversion of Potential Fields: The Quebec Oka Carbonatite Complex Case Study

Potential fields methods, based on the exploitation of gravity and magnetic fields, are among the most important methods to recover fundamental information on the Earth crust structure at global, regional and local scales. The bottleneck for this kind of geophysical methods is often represented by the development of ad-hoc techniques to fully exploit the available data. In fact, each different technique can observe the effect of a single property of the subsurface and when we want to estimate this property from the observed field (the so-called inverse problem), several problems such as non-uniqueness and instability arise. A possible solution to these problems consists in jointly inverting, in a consistent way, different observed fields, possibly also incorporating all the available geological constraints. In the current work, we present an innovative Bayesian algorithm aimed at performing a full 3D joint inversion of gravity and magnetic fields constrained by geological a-priori qualitative information. The algorithm is tested on a real-case scenario, namely, a local study to estimate a complete 3D model of the Oka carbonatite complex. This complex is a composite pluton in Quebec (Canada), important for mining operations related to critical raw material such as Niobium and other rare earth. This example shows the reliability of the developed inversion algorithm and gives hints on the fundamental role that potential fields can play in mining activities.

Effective 3D-transient electromagnetic inversion using finite-element method with a parallel direct solver

The correct interpretation of time-domain or transient electromagnetic (TEM) data relies on effective inversion techniques. Herein, we have developed an effective 3D-TEM inversion method based on the vector finite-element method. The domain was discretized with unstructured tetrahedra to simulate complex geologic structures. We develop a second-order backward-Euler scheme, rigorously with nonuniform time-step sizes, to approximate the time derivative of the electric field. The waveform was directly simulated and considered in the forward modeling and inversion. In a 3D inversion, the most time-consuming element of the workflow is solving a large system of linear equations. To alleviate this problem, we adopt a parallel direct solver for solving a large system of equations involved in forward modeling and sensitivity calculations. We use the Gauss-Newton approach to solve the inverse problem and formulate the least-squares problem to obtain model updates. We avoided the explicit sensitivity calculation by solving the matrix-vector multiplication between the sensitivity matrix (or its transpose) and vector. The developed inversion algorithm is validated using several realistic synthetic models. Finally, we apply the method to practical TEM data with complex topography for coal mine goaf detection. We thus obtain a detailed 3D resistivity model that revealed the layering structure and some localized anomalies related to the coal mine goaf.

Gradient-Based Electromagnetic Inversion for Metasurface Design Using Circuit Models

A gradient-based optimization algorithm that is capable of directly designing a metasurface at the circuit parameter level for a desired field (amplitude and phase) pattern or a desired power (phaseless) pattern on some region of interest (ROI) external to the metasurface boundary is presented. Specifically, the inversion algorithm designs the microwave admittance profile of each subwavelength element of the metasurface when a three-layer admittance model is assumed. To this end, a forward model that maps the admittance profile of each layer of the metasurface to the desired field on the ROI is developed. Then, for the inverse design problem, a data misfit cost functional is defined and minimized over the unknown admittance profile using analytically derived gradients and step lengths. The developed inversion algorithm is then utilized to design metasurfaces capable of beam forming in the far-field zone.

3-D Magnetotelluric Inversion and Application Using the Edge-Based Finite Element With Hexahedral Mesh

Three-dimensional (3-D) inversion technique has become an important and practical approach for magnetotelluric (MT) data interpretation. In this article, we developed a 3-D parallelized MT inversion scheme using the edge-based finite element method and applied the developed method to the newly collected MT data in the Xinjiang Luntai area. The distorted hexahedral element is adopted to incorporate topography into the forward modeling and inversion for complicated scenarios. We use the Gauss–Newton optimization method to minimize the objective functional for MT inversion. The developed algorithm is parallelized using MPI over frequencies and parallel direct solvers when solving the forward and adjoint problems for each frequency. We compare the performance of the least-square QR (LSQR) factorization and preconditioned conjugate gradient (PCG) solvers for the model update within each Gauss–Newton iteration and found that the LSQR solver is more stable. The developed inversion algorithm is validated using several synthetic models. Finally, we applied the inversion algorithm to the subsurface resistivity imaging in the Luntai area. The recovered geoelectric model from full 3-D inversion fits well with the known geological and geophysical information. The recovered model shows a low resistivity layer which may be caused by the salt strata. Besides, the inversion results reveal the movement tectonic in this survey area within a depth of 9 km.

Satellite-based view of the aerosol spatial and temporal variability in the Córdoba region (Argentina) using over ten years of high-resolution data

Space-based observations offer a unique opportunity to investigate the atmosphere and its changes over decadal time scales, particularly in regions lacking in situ and/or ground based observations. In this study, we investigate temporal and spatial variability of atmospheric particulate matter (aerosol) over the urban area of Córdoba (central Argentina) using over ten years (2003–2015) of high-resolution (1 km) satellite-based retrievals of aerosol optical depth (AOD). This fine resolution is achieved exploiting the capabilities of a recently developed inversion algorithm (Multiangle implementation of atmospheric correction, MAIAC) applied to the MODIS sensor datasets of the NASA-Terra and -Aqua platforms. Results of this investigation show a clear seasonality of AOD over the investigated area. This is found to be shaped by an intricate superposition of aerosol sources, acting over different spatial scales and affecting the region with different yearly cycles. During late winter and spring (August-October), local as well as near- and long-range transported biomass burning (BB) aerosols enhance the Córdoba aerosol load, and AOD levels reach their maximum values (>0.35 at 0.47 µm). The fine AOD spatial resolution allowed to disclose that, in this period, AOD maxima are found in the rural/agricultural area around the city, reaching up to the city boundaries pinpointing that fires of local and near-range origin play a major role in the AOD enhancement. A reverse spatial AOD gradient is found from December to March, the urban area showing AODs 40–80% higher than in the city surroundings. In fact, during summer, the columnar aerosol load over the Córdoba region is dominated by local (urban and industrial) sources, likely coupled to secondary processes driven by enhanced radiation and mixing effects within a deeper planetary boundary layer (PBL). With the support of modelled AOD data from the Modern-Era Retrospective Analysis for Research and Application (MERRA), we further investigated into the chemical nature of AOD. The results suggest that mineral dust is also an important aerosol component in Córdoba, with maximum impact from November to February. The use of a long-term dataset finally allowed a preliminary assessment of AOD trends over the Córdoba region. For those months in which local sources and secondary processes were found to dominate the AOD (December to March), we found a positive AOD trend in the Córdoba outskirts, mainly in the areas with maximum urbanization/population growth over the investigated decade. Conversely, a negative AOD trend (up to −0.1 per decade) is observed all over the rural area of Córdoba during the BB season, this being attributed to a decrease of fires both at the local and the continental scale.

In silico simulation and in vitro evaluation of an elastomeric scaffold using ultrasonic shear wave imaging

Biodegradable elastomeric scaffolds for soft tissue repair represent a growing area of biomaterials research. Mechanical strength is one of the key factors to consider in the evaluation of candidate materials and the designs for tissue scaffolds. It is desirable to develop non-invasive evaluation methods of the mechanical property of scaffolds which would provide options for monitoring temporal mechanical property changes in situ. In this paper, we conduct in silico simulation and in vitro evaluation of an elastomeric scaffold using a novel ultrasonic shear wave imaging (USWI). The scaffold is fabricated from a biodegradable elastomer, poly(carbonate urethane) urea using salt leaching method. A numerical simulation is performed to test the robustness of the developed inversion algorithm for the elasticity map reconstruction which will be implemented in the phantom experiment. The generation and propagation of shear waves in a homogeneous tissue-mimicking medium with a circular scaffold inclusion is simulated and the elasticity map is well reconstructed. A PVA phantom experiment is performed to test the ability of USWI combined with the inversion algorithm to non-invasively characterize the mechanical property of a porous, biodegradable elastomeric scaffold. The elastic properties of the tested scaffold can be easily differentiated from the surrounding medium in the reconstructed image. The ability of the developed method to identify the edge of the scaffold and characterize the elasticity distribution is demonstrated. Preliminary results in this pilot study support the idea of applying the USWI based method for non-invasive elasticity characterization of tissue scaffolds.

Approximate solutions of acoustic 3D integral equation and their application to seismic modeling and full-waveform inversion

Over the recent decades, a number of fast approximate solutions of Lippmann–Schwinger equation, which are more accurate than classic Born and Rytov approximations, were proposed in the field of electromagnetic modeling. Those developments could be naturally extended to acoustic and elastic fields; however, until recently, they were almost unknown in seismology. This paper presents several solutions of this kind applied to acoustic modeling for both lossy and lossless media. We evaluated the numerical merits of those methods and provide an estimation of their numerical complexity. In our numerical realization we use the matrix-free implementation of the corresponding integral operator. We study the accuracy of those approximate solutions and demonstrate, that the quasi-analytical approximation is more accurate, than the Born approximation. Further, we apply the quasi-analytical approximation to the solution of the inverse problem. It is demonstrated that, this approach improves the estimation of the data gradient, comparing to the Born approximation. The developed inversion algorithm is based on the conjugate-gradient type optimization. Numerical model study demonstrates that the quasi-analytical solution significantly reduces computation time of the seismic full-waveform inversion. We also show how the quasi-analytical approximation can be extended to the case of elastic wavefield.

Fast 1D inversion of logging-while-drilling resistivity measurements for improved estimation of formation resistivity in high-angle and horizontal wells

We have developed an efficient inversion method to estimate layer-by-layer electric resistivity from logging-while-drilling electromagnetic induction measurements. The method assumes a 1D model based on planarly layered transversely isotropic formations with known bed boundaries, penetrated by arbitrary well trajectories. Forward simulations are based on a 1D reduced model in which borehole and mandrel effects are assumed to be negligible. The stopping criteria, regularization term, regularization parameter, the Jacobian matrix, and inversion variables are automatically estimated by the inversion algorithm, thereby minimizing the required user interaction and expertise. Numerical inversion results of challenging synthetic and field data confirmed the high stability and superior approximation properties of the developed inversion algorithm. We evaluated results indicating that triaxial induction measurements provide significantly more stable inversion results than conventional coaxial measurements, especially in the presence of anisotropic formations.

Application of Cauchy-type integrals in developing effective methods for depth-to-basement inversion of gravity and gravity gradiometry data

One of the most important applications of gravity surveys in regional geophysical studies is determining the depth to basement. Conventional methods of solving this problem are based on the spectrum and/or Euler deconvolution analysis of the gravity field and on parameterization of the earth’s subsurface into prismatic cells. We have developed a new method of solving this problem based on 3D Cauchy-type integral representation of the potential fields. Traditionally, potential fields have been calculated using volume integrals over the domains occupied by anomalous masses subdivided into prismatic cells. This discretization can be computationally expensive, especially in a 3D case. The technique of Cauchy-type integrals made it possible to represent the gravity field and its gradients as surface integrals. In this approach, only the density contrast surface between sediment and basement needed to be discretized for the calculation of gravity field. This was especially significant in the modeling and inversion of gravity data for determining the depth to the basement. Another important result was developing a novel method of inversion of gravity data to recover the depth to basement, based on the 3D Cauchy-type integral representation. Our numerical studies determined that the new method is much faster than conventional volume discretization method to compute the gravity response. Our synthetic model studies also showed that the developed inversion algorithm based on Cauchy-type integral is capable of recovering the geometry and depth of the sedimentary basin effectively with a complex density profile in the vertical direction.

분해능 모델 제한자를 사용하는 자력탐사자료의 2차원 역산

새로운 2차원 자력탐사자료 역산 방법을 개발하였다. 중,자력탐사와 같은 포텐셜 자료의 역산에서 가장 문제가 되는 점은 비유일해 문제이다. 일반적으로 자력탐사 자료의 역산은 모델변수의 수가 자료의 수보다 훨씬 많은 불충분 문제이며, 이는 비유일해 문제를 더욱 심화시키게 된다. 일반적인 최소자승법을 자력탐사자료의 역산에 적용하게 되면, 이 상대가 지표면에 집중되는 결과를 초래한다. 본 연구에서는 이러한 비유일해 문제를 극복하기 위하여 모델분해능에 근거한 새로운 모델제한자를 제안하였다. 이 모델제한자는 분해능이 높은 모델변수에는 큰 제한을 가하고, 작은 모델변수에는 약한 제한을 가하게 된다. 따라서 분해능이 낮은 심부의 모델변수도 효과적으로 추정할 수 있다. 개발된 역산 알고리듬을 이용하여, 전형적인 모델에 대한 이론자료의 역산에 적용하였다. 또한 옥천대에서 얻어진 항공자력탐사자료 역산에 적용하였다. We developed a method for inverting magnetic data to image 2D susceptibility models. The major difficulty in the inversion of the potential data is the nonuniqueness. Furthermore, generally the number of inversion blocks are greater than the number of the magnetic data available, and thus the magnetic inversion leads to under-determined problem, which aggravates the nonuniqueness. When the magnetic data were inverted by the general least-squares method, the anomalous susceptibility would be concentrated near the surface in the inverted section. To overcome this nonuniqueness problem, we propose a new resolution model constraint that is calculated from the parameter resolution. The model constraint imposes large penalty on the model parameter with good resolution, on the other hand small penalty on the model parameter with poor resolution. Thus, the deep-seated model parameter, generally having poor resolution, can be effectively resolved. The developed inversion algorithm is applied to the inversion of the synthetic data for typical models of magnetic anomalies and is tested on real airborne data obtained at the Okcheon belt of Korea.

자력이상 3차원 모델링 및 역산

자력탐사자료의 3차원 역산법을 개발하였다. 자력탐사자료의 역산에서 가장 문제가 되는 점은 비유일해 문제와 방대한 계산시간이다. 일반적으로 자력탐사자료의 역산은 모델변수의 수가 자료의 수보다 훨씬 많아 비유일해 문제를 더욱 심화시키게 된다. 또한 자력탐사자료는 심도 분해능이 매우 낮다. 비유일해 문제를 극복하기 위하여 분해능이 높은 모델변수에는 큰 제한을 가하고, 작은 모델변수에는 약한 제한을 가하는 분해능 모델제한자를 제안하고, 이를 적용하여 분해능이 낮은 모델변수도 효과적으로 추정할 수 있었다. 또한 대형 행렬식을 웨이블릿 변환을 통하여 희소행렬로 변환하고, 역행렬의 계산에 병렬계산 방식을 적용하여 계산시간을 획기적으로 절감하였다. 수치실험을 통하여 개발된 3차원 역산알고리듬의 타당성을 검토하였다. 또한 금산 지역에서 얻어진 항공자력탐사자료의 역산에 적용하였다. We developed a method for inverting magnetic data to recover the 3D susceptibility models. The major difficulty in the inversion of the potential data is the non-uniqueness and the vast computing time. The insufficient number of data compared with that of inversion blocks intensifies the non-uniqueness problem. Furthermore, there is poor depth resolution inherent in magnetic data. To overcome this non-uniqueness problem, we propose a resolution model constraint that imposes large penalty on the model parameter with good resolution; on the other hand, small penalty on the model parameter with poor resolution. Using this model constraint, the model parameter with a poor resolution can be effectively resolved. Moreover, the wavelet transform and parallel solving were introduced to save the computing time. Through the wavelet transform, a large system matrix was transformed to a sparse matrix and solved by a parallel linear equation solver. This procedure is able to enormously save the computing time for the 3D inversion of magnetic data. The developed inversion algorithm is applied to the inversion of the synthetic data for typical models of magnetic anomalies and real airborne data obtained at the Geumsan area of Korea.

Determination of Q β(f) in Different Parts of Kumaon Himalaya from the Inversion of Spectral Acceleration Data

This paper presents the results of a modified two-step inversion algorithm approach to find S wave quality factor Qβ(f) given by Joshi (Bull Seis Soc Am 96:2165–2180, 2006). Seismic moment is calculated from the source displacement spectra of the S wave using both horizontal components. Average value of seismic moment computed from two horizontal components recorded at several stations is used as an input to the first part of inversion together with the spectra of S phase in the acceleration record. Several values of the corner frequency have been selected iteratively and are used as inputs to the inversion algorithm. Solution corresponding to minimum root mean square error (RMSE) is used for obtaining the final estimate of Qβ(f) relation. The estimates of seismic moment, corner frequency and Qβ(f) from the first part of inversion are further used for obtaining the residual of theoretical and observed source spectra which are treated as site amplification terms. The acceleration record corrected for the site amplification term is used for determination of seismic moment from source spectra by using Qβ(f) obtained from first part of inversion. Corrected acceleration record and new estimate of seismic moment are used as inputs to the second part of the inversion scheme which is similar to the first part except for use of input data. The final outcome from this part of inversion is a new Qβ(f) relation together with known values of seismic moment and corner frequency of each input. The process of two-step inversion is repeated for this new estimate of seismic moment and goes on until minimum RMSE is obtained which gives final estimate of Qβ(f) at each station and corner frequency of input events. The Pithoragarh district in the state of Uttarakhand in India lies in the border region of India and Nepal and is part of the seismically active Kumaon Himalaya zone. A network of eight strong motion recorders has been installed in this region since March, 2006. In this study we have analyzed data from 18 local events recorded between March, 2006 and October, 2010 at various stations. These events have been located using HYPO71 and data has been used to obtain frequency-dependent shear-wave attenuation. The Qβ(f) at each station is calculated by using both the north-south (NS) and east-west (EW) components of acceleration records as inputs to the developed inversion algorithm. The average Qβ(f) values obtained from Qβ(f) values at different stations from both NS and EW components have been used to compute a regional average relationship for the Pithoragarh region of Kumaon Himalaya of form Qβ(f) = (29 ± 1.2)f(1.1 ± 0.06).

A 3D parametric inversion algorithm for triaxial induction data

We develop a parametric inversion algorithm to determine simultaneously the horizontal and vertical resistivities of both the formation and invasion zones, invasion radius, bed boundary upper location and thickness, and relative dip angle from electromagnetic triaxial induction logging data. This is a full 3D inverse scattering problem in transversally isotropic media. To acquire sufficient sensitivity to invert for all of these parameters, we collect the data using a multicomponent, multispacing induction array. For each transmitter-receiver spacing this multicomponent tool has sets of three orthogonal transmitter and receiver coils. At each logging point single-frequency data are collected at multiple spacings to obtain information at different depths of investigation. This inversion problem is solved iteratively with a constrained regularized Gauss-Newton minimization scheme. As documented in the literature, the main computational bottleneck when solving this full 3D inverse problem is the CPU time associated with constructing the Jacobian matrix. In this study, to achieve the inversion results within a reasonable computational time, we implement a dual grid approach wherein the Jacobian matrix is computed using a very coarse optimal grid. Furthermore, to regularize the inversion process we use the so-called multiplicative regularization technique. This technique automatically determines the regularization parameter. Synthetic data tests indicate the developed inversion algorithm is robust in extracting formation and invasion anisotropic resistivities, invasion radii, bed boundary locations, relative dip, and azimuth angle from multispacing, multicomponent induction logging data.

Joint inversion of Surface Nuclear Magnetic Resonance and Vertical Electrical Sounding

The method of Surface Nuclear Magnetic Resonance (SNMR) provides a very new technology to directly determine subsurface water distribution. The microscopic magnetization of water molecules is used to derive water content and pore size information from SNMR soundings. The observed similarity and agreement between interpreted aquifer structure from SNMR and resistivity distribution from Vertical Electrical Sounding (VES) has led to our objective to jointly invert both data sets using a generalized petrophysical model based on Archie's Law. To perform inversion of both methods, the Simulated Annealing (SA) technique was applied. Since a very fast numerical solution is available for both geophysical methods, this kind of guided random search algorithm promises better performance than least square methods. The developed inversion algorithm has been applied on a number of different synthetic data to study its properties and prove its reliability. Investigations on well-known test sites where both methods were conducted finally proved the effectiveness of the joint inversion on real data. The interpretation of the subsurface model could be optimized beyond an enhanced spatial resolution to a quantitative interpretation of the ratio of mobile and adhesive water contents, leading to prediction of hydrological parameters from geophysical investigations.

Thin-film elastic-property measurements with laser-ultrasonic SAW spectrometry

We have developed laser-ultrasonic methods to evaluate the elastic properties of thin films. Our method, called SAW spectrometry, used broadband surface acoustic waves (SAWs) generated by a line-focused, pulsed laser. SAWs with frequency components up to ∼400 MHz were detected by a Michelson interferometer. Dispersion relations for the phase velocity vs. frequency were calculated from displacement waveforms acquired over a range of source–detector separations. Quantitative elastic-property information was obtained from the dispersion relations with a newly developed inversion algorithm. The algorithm was based on the delta-function representation of the elastodynamic Green's function for anisotropic layered systems. The methods were tested by applying them to a molybdenum film on a single-crystal silicon wafer. Values obtained for the film's Poisson ratio and Young's modulus were consistent with literature values. We also used SAW spectrometry to investigate several samples containing titanium nitride films on silicon substrates. Assuming literature values for the film density and Poisson ratio, we applied the inversion algorithm to determine the thickness and Young's modulus of each film. Young's modulus in the films was found to increase monotonically with increasing film thickness from 368 to 453 GPa. Inversion results for the film thickness were in very good agreement with destructive measurements of the actual thicknesses.

A new method of retrieving surface bidirectional reflectance from ground measurements: Atmospheric sensitivity study

Commonly, experimental data on surface directional reflectance are not properly corrected for atmospheric effects. In this paper, we describe a new method for rigorous retrieval of the bidirectional reflectance distribution function (BRDF) from ground‐based measurements. It obtains the pointwise BRDF in the view directions and the best fit parameters of the approximating analytical model. Unlike previously developed algorithms, it can perform a retrieval using directional observations made at a single solar angle. This improvement allows one to process experiments with incomplete sets of measurements, such as those contaminated by clouds. The paper describes the performance of the method and its sensitivity to atmospheric parameters. Using rigorous direct modeling of the reflected radiance and developed inversion algorithm, we performed an analysis of BRDF errors, typical for measurements. It shows that the BRDF shape distortions due to diffuse irradiance are considerable in the visible and near‐IR, even for low aerosol content in the atmosphere. This result confirms that an accurate atmospheric correction should be a required processing step in the experimental research of anisotropic surface reflectance.

SAGE II inversion algorithm

This paper provides a detailed description of the current operational SAGE II multichannel data inversion algorithm implemented at NASA Langley Research Center, Hampton, Virginia. This algorithm is compared to an independently developed inversion algorithm from the Laboratory of Atmospheric Optics, University of Lille, Lille, France. Inverted aerosol and ozone profiles from these two algorithms are shown to be similar within their respective uncertainties.

Shaping of arbitrary dose distributions by dynamic multileaf collimation

Traditionally, the shaping of non-uniform dose distributions has been performed by using wedges or compensating filters. The advent of high resolution multileaf collimators may largely eliminate the need for material attenuators for modification of the beam. This is achieved by a new technique for the shaping of arbitrary dose distributions by dynamic motion of the collimator leaves. By employing narrow elementary slit beams that correspond to the smallest possible opening of the multileaf collimator, the optimal density of such slit beams, i.e. opening density, can be determined automatically using a newly developed inversion algorithm. The present method has two major advantages: (1) internal structures in the field can be created, controlled solely by steering the collimator leaves, (2) the opening density determined by the algorithm never gives rise to underdosage: this is important from a radiobiological point of view.