Least-squares Inversion Method Research Articles

An efficient approach for inverting rock exhumation from thermochronologic age–elevation relationship

Abstract. This study implements the least-squares inversion method for solving the exhumation history from the thermochronologic age–elevation relationship (AER) based on the linear equation among exhumation rate, age and total exhumation from the closure depth to the Earth surface. Modeling experiments suggest significant and systematic influence of initial geothermal model, the a priori exhumation rate and the time interval length on the a posteriori exhumation history. Lessons learned from the experiments include that (i) the modern geothermal gradient can be used for constraining the initial geothermal model, (ii) a relatively high a priori exhumation rate would lead to systematically lower a posteriori exhumation and vice versa, (iii) the variance of the a priori exhumation rate controls the variation in the inverted exhumation history, and (iv) the choice of time interval length should be optimized for resolving the potential temporal changes in exhumation. To mitigate the dependence of inverted erosion history on these initial parameters, we implemented a new stepwise inverse modeling method for optimizing the model parameters by comparing the observed and predicted thermochronologic data and modern geothermal gradients. Finally, method demonstration was performed using four synthetic datasets and three natural examples of different exhumation rates and histories. It is shown that the inverted rock exhumation histories from the synthetic datasets match the whole picture of the “truth”, although the temporal changes in the magnitude of exhumation are underestimated. Modeling of the datasets from natural samples produces geologically reasonable exhumation histories. The code and data used in this work are available on Zenodo (https://doi.org/10.5281/zenodo.10839275).

Seismic Imaging of Lithospheric Structure Beneath Central-East Java Region, Indonesia: Relation to Recent Earthquakes

The Central and East Java region, which is part of the Sunda Arc, has an important role in producing destructive earthquakes and volcanic complexes in Indonesia as a result of the convergence between the Indo-Australian plate that subducts under the Eurasian plate. In this study, the 3-D seismic velocity structure (Vp, Vs, and Vp/Vs) of the crust and upper mantle was determined to reveal the presence of the subducted slab, volcanic sources, and seismogenic features in the Central and East Java region. We have manually re-picked P- and S-arrivals of 1,488 events from January 2009 to September 2017 recorded at 27 stations of the BMKG network. An iterative damped least-squares inversion method was applied to simultaneously calculate both hypocenter relocations and velocity structure beneath this high-risk region to a depth of 200 km. We then compare the tomographic results and seismicity to interpret structural features in the seismic zones. The subducted slab is dipping toward the north, imaged by the high-velocity regions with low Vp/Vs at depths of about 50–100 km. Low-velocity anomalies with high Vp/Vs above the slab at a depth of ∼100 km, imply the possible location of partial melting from slab dehydration. Fluids and melts are ascending to feed the volcanoes i. e., Merapi-Merbabu, Wilis, Pandan, Semeru, Bromo, and Ijen that also have similar low-velocity anomalies at 10–30 km depths, suggesting the presence of the sedimentary basin or magma reservoir. We also have redetermined the hypocenter location of the 2021 (Mw 6.1) Malang earthquake at 8.94oS, 112.45oE, with a depth of 59.7 km. The location error in the x, y, and z directions are 3.08, 6.39, and 11.91 km, respectively. This intraslab event with a thrusting mechanism is located in the high-velocity region and close to the intermediate-depth seismic clusters, which indicates the geometry of the oceanic slab. In the region of 1994 (Mw 7.8) Banyuwangi earthquake, we found a low-velocity anomaly at ∼ 50 km depth that might be associated with the presence of subducting seamount that is more hydrated than the surrounding slab. The slip over this subducting seamount caused the tsunamigenic earthquake.

Application of the Electrical Resistivity Method for Site Investigation in University of Anbar, Ar-Ramadi City, Western Iraq

The 2D resistivity imaging technique was applied in an engineering study for the investigation of subsurface weakness zones within University of Anbar, western Iraq. The survey was carried out using Dipole-dipole array with an n-factor of 6 and a-spacing values of 2 m and 5 m. The inverse models of the 2D electrical imaging clearly show the resistivity contrast between the anomalous parts of the weakness zones and the background resistivity distribution. The thickness and shape of the subsurface weakness zones were well defined from the 2D imaging using Dipole-dipole array of 2 m a-spacing. The thickness of the weakness zone ranges between 9.5 m to 11.5 m. Whereas the Dipole-dipole array with a-spacing of 5 m and n-factor of 6 allocated the geoelectrical stratigraphic layers sequence in low-accuracy of weakness zones, but deeper than the inverse model of 2 m a-spacing. This survey was made to explain the correlation between the weakness zone and the deeper layers in the study area. It points out that the deeper layers were not affected in the weakness zones. The inverse model was produced using the Standard Least-Squares Inversion Method and the Robust Inversion Model Constraints Method. The first method had a gradational boundary of the weakness zones and the second had sharper and straighter boundaries of fractures and voids within the weakness zones.

A Multi-Observation Least-Squares Inversion for GNSS-Acoustic Seafloor Positioning

Monitoring deformation on the seafloor is a major challenge for modern geodesy and a key to better understanding tectonic processes and assess related hazards. The extension of the geodetic networks offshore can be achieved by combining satellite positioning (GNSS) of a surface platform with acoustic ranging to seafloor transponders. This approach is called GNSS-Acoustic (GNSS-A). The scope of this work is to provide a tool to identify and quantify key points in the error budget of such experiment. For this purpose, we present a least-squares inversion method to determine the absolute position of a seafloor transponder array. Assuming the surface platform is accurately positioned by GNSS, the main observables are the two-way travel time in water between the transponders on the seafloor and the surface platform acoustic head. To better constrain transponder positions, we also consider the baseline lengths and the relative depth-differences between different pairs of them. We illustrate the usefulness of our forward modeling approach and least-square inversion by simulating different experimental protocols (i.e., platform trajectories, with or without information on the distance and depth between transponders). We find that the overall accuracy of a GNSS-A experiment is significantly improved with additional information about the relative depths of the instruments. Baseline lengths also improve the accuracy, but only when combined with depth differences. The codes in Python3 used in this article are freely available online.

Measurement of a very large negative reactivity inserted by rapid withdrawal of a partial fuel loading in Kyoto University Critical Assembly

ABSTRACT Solid-moderated cores of the Kyoto University Critical Assembly (KUCA) have an emergency shutdown system to avoid any reactivity-initiated accidents. When a reactor protection system of KUCA detects some abnormal conditions, the emergency system withdraws downward several fuel and reflector assemblies assigned as a central fuel loading from the core in 10 and several seconds. The withdrawal of the fuel loading inserts a very large negative reactivity into the core. The negative reactivity has been measured by the integral count technique; however, the application of the integral count technique to this measurement must have been impossible in principle. The integral count technique assumes an instantaneous insertion of a large negative reactivity but the assumption must fail for this measurement. An accurate measurement of the negative reactivity is of extreme importance from the viewpoint of reactor safety. In this study, the least-squares inverse kinetics method was employed to determine the negative reactivity from time-sequence counts data of a neutron counter. The negative reactivity determined by this method was in good agreement with that calculated by the continuous-energy Monte Carlo code MVP (version 3) with the nuclear data library JENDL-4.0.

A least-squares correlation-based full traveltime inversion for shallow subsurface velocity reconstruction

Accurate estimation of shallow subsurface velocity models with complex topography is crucial for statics corrections and imaging deep structures. The correlation-based wave-equation traveltime inversion (CWTI) method is suitable for estimating such shallow subsurface velocity structures. However, the CWTI objective function suffers an inherent resolution-loss problem because the traveltime weighted crosscorrelation misfit does not fall to zero even when the model is perfectly matched. Furthermore, the Born-approximation-based CWTI gradient cannot provide an effective model update during each iteration. To overcome these problems, we have developed a least-squares correlation-based full traveltime inversion (LCFTI) method, in which the least-squares correlation-based objective function was designed to minimize the traveltime weighted difference between the autocorrelation and the crosscorrelation. By incorporating the autocorrelation, LCFTI indicates better convergence and higher resolution than CWTI. The LCFTI model updates are derived using the Rytov approximation to avoid incorrect model updates by emphasizing phase matching. Furthermore, to accurately simulate wave propagation, we use the spectral-element method as the modeling engine, in which the mesh of the complex topography is flexibly represented. Synthetic and field data examples are performed to demonstrate the effectiveness of the proposed method in shallow subsurface velocity reconstruction.

Warping least-squares inversion for 4D velocity change: Gulf of Mexico case study

We have developed a new trace-based, warping least-squares inversion method to quantify 4D velocity changes. There are two steps to solve for these velocity changes: (1) dynamic warping with phase constraints to align the baseline and monitor traces and (2) least-squares inversion for 4D velocity changes incorporating the time shifts and 4D amplitude differences (computed after trace alignment by warping). We have demonstrated this new inversion workflow using simple synthetic layered models. For the noise-free case, phase-constrained warping is superior to standard, amplitude-based warping by improving trace alignment, resulting in more accurate inverted velocity changes (less than 1% error). For synthetic data with 6% rms noise, inverted velocity changes are reasonably accurate (less than 10% error). Additional inversion tests with migrated finite-difference data shot over a realistic anticline model result in less than 10% error. The inverted velocity changes on a 4D field data set from the Gulf of Mexico are more interpretable and consistent with the dynamic reservoir model than those estimated from the conventional time-strain method.

Origin and causes of complex magnetic anomalies of Dokdo and neighboring guyots

To identify the causes of complex magnetic anomalies on the three guyots located in the northeastern part of the Ulleung Basin, East Sea, a linear least-squares inversion method was implemented. Modeling under the assumption that the guyots are uniformly magnetized did not yield reliable results. However, an inhomogeneously-magnetized model shows relatively high goodness-of-fit parameters and low residual magnetic anomalies, indicating that these guyots are considerably nonuniformly magnetized. Our modeling results show that all three guyots show similar magnetic patterns. Magnetic anomalies north of guyots have different magnetic polarities at similar depths in the south, indicating that these magnetic anomalies were formed at different times. These results suggest that relatively short wavelengths and complicated magnetic anomalies in the study area are the compound result of different magnetic polarities. Paleomagnetic analysis on trachyandesite lavas and pyroclastic breccia indicates a thermo-chemical remagnetization activity after Reunion Subchron (2.138–2.122 Ma). This result also suggests that remagnetization contributes to the creation of complex magnetic anomalies to some extent as well. In addition, the upward continuation results up to 1500 m using a FFT routine show high positive magnetic anomalies with a linearly segmented shape at the base of guyots. These linear and positive magnetic anomalies are interpreted as volcanic eruptions along fissures because mafic lava should have much higher magnetic anomalies compared to continental crust, where guyots are located. Therefore, we interpret that the complex magnetic anomalies distributed in the study area are due to polarity reversals occurring at least twice, and were formed during the formation of the current guyots as volcanic eruptions derived from fissures.

Characterization of Pan-African Aquifer Layers by the Least Squares Inversion Method Applied on Geoelectric Data

Geoelectric data obtained from forty (40) vertical electrical soundings collected with a Schlumberger device in the Adamawa plateau region, also known as the Cameroon water tower, have been treated by the least-squares inversion method. In order to study the nature and thickness of the aquifer and the necessary geoelectric parameters, quantitative and qualitative interpretations of the data were made. The results obtained showed that: about four to five geoelectric layers have been delimited in the study area with a dominance of the KH curve, which can be used as a reference for future studies. The first two layers constitute an association of clay and laterite with resistivity values ranging from 58 to 9122 Ω&middotm and whose thickness is between 0.6 and 13.4 m. The third layer is a potentially aquiferous laterite composed of clay, laterite and especially clay sand and cracked/good granite, with a dominance of sandy alteration whose resistivity values are between 81 and 960 Ω&middotm and its thickness between 12.2 and 26.8 m. The fourth and fifth layers are made up of cracked/good granite with a resistivity ranging from 12-10705 Ω&middotm with an average value of 1817 Ω&middotm. This study also shows that the North-East, South-West and South sectors could be the groundwater convergence zones and that the average depth of the basement aquifer roof is about 28.3 m. The geoelectric sections of certain demarcated vertical electrical sounding stations are consistent with the geologic description of the area.

Seismic resolution enhancement in shale-oil reservoirs

We developed a case study of seismic resolution enhancement for shale-oil reservoirs in the Q Depression, China, featured by rhythmic bedding. We proposed an innovative method for resolution enhancement, called the full-band extension method. We implemented this method in three consecutive steps: wavelet extraction, filter construction, and data filtering. First, we extracted a constant-phase wavelet from the entire seismic data set. Then, we constructed the full-band extension filter in the frequency domain using the least-squares inversion method. Finally, we applied the band extension filter to the entire seismic data set. We determined that this full-band extension method, with a stretched frequency band from 7–70 to 2–90 Hz, may significantly enhance 3D seismic resolution and distinguish reflection events of rhythmite groups in shale-oil reservoirs.

A new prestack three-parameter amplitude variation with offset inversion method

Prestack seismic inversion is an important method in reservoir prediction and fluid identification. In this paper, we develop a new strategy to carry out a prestack three-parameter AVO inversion to obtain P-wave velocity, S-wave velocity, and density. Using the Aki and Richards equation to calculate the PP reflection coefficients, which is the approximation to the Zoeppritz PP equation, we re-derive the relationship between reflection coefficients and their partial derivatives with respect to P-wave velocity, S-wave velocity, and density. We transform the objective function into the form of matrix multiplication and solve it using a damped least-squares method and Tikhonov regularization method. Then we directly obtain the inversion results, including P-wave velocity, S-wave velocity and density. To demonstrate the validity of the proposed prestack inversion method, we set up a 1D forward simulation model and add uniform noise into the synthetic seismic records. Then we test the proposed method on the overthrust model and field data. In order to demonstrate its high accuracy and the efficiency of the Tikhonov regularization method, we carry out a damped least-squares inversion method using the same models as a comparison. We find that the proposed prestack seismic data AVO inversion method has higher anti-noise capability, accuracy, and efficiency, both in synthetic records and field data using the Tikhonov regularization method.

Performance Analysis of Low-Flux Least-Squares Single-Pixel Imaging

A single-pixel camera is able to computationally form spatially resolved images using one photodetector and a spatial light modulator. The images it produces in low-light-level operation are imperfect, even when the number of measurements exceeds the number of pixels, because its photodetection measurements are corrupted by Poisson noise. Conventional performance analysis for single-pixel imaging generates estimates of mean-square error (MSE) from Monte Carlo simulations, which require long computational times. In this letter, we use random matrix theory to develop a closed-form approximation to the MSE of the widely used least-squares inversion method for Poisson noise-limited single-pixel imaging. We present numerical experiments that validate our approximation and a motivating example showing how our framework can be used to answer practical optical design questions for a single-pixel camera.

Free-surface multiple attenuation for blended data

Blended data sets are now being acquired because of improved efficiency and reduction in cost compared with conventional seismic data acquisition. We have developed two methods for blended data free-surface multiple attenuation. The first method is based on an extension of surface-related multiple elimination (SRME) theory, in which free-surface multiples of the blended data can be predicted by a multidimensional convolution of the seismic data with the inverse of the blending operator. A least-squares inversion method is used, which indicates that crosstalk noise existed in the prediction result due to the approximate inversion. An adaptive subtraction procedure similar to that used in conventional SRME is then applied to obtain the blended primary — this can damage the energy of primaries. The second method is based on inverse data processing (IDP) theory adapted to blended data. We derived a formula similar to that used in conventional IDP, and we attenuated free-surface multiples by simple muting of the focused points in the inverse data space (IDS). The location of the focused points in the IDS for blended data, which can be calculated, is also related to the blending operator. We chose a singular value decomposition-based inversion algorithm to stabilize the inversion in the IDP method. The advantage of IDP compared with SRME is that, it does not have crosstalk noise and is able to better preserve the primary energy. The outputs of our methods are all blended primaries, and they can be further processed using blended data-based algorithms. Synthetic data examples show that the SRME and IDP algorithms for blended data are successful in attenuating free-surface multiples.

Application of the least-squares inversion method: Fourier series versus waveform inversion

We describe an implicit link between waveform inversion and Fourier series based on inversion methods such as gradient, Gauss–Newton, and full Newton methods. Fourier series have been widely used as a basic concept in studies on seismic data interpretation, and their coefficients are obtained in the classical Fourier analysis. We show that Fourier coefficients can also be obtained by inversion algorithms, and compare the method to seismic waveform inversion algorithms. In that case, Fourier coefficients correspond to model parameters (velocities, density or elastic constants), whereas cosine and sine functions correspond to components of the Jacobian matrix, that is, partial derivative wavefields in seismic inversion. In the classical Fourier analysis, optimal coefficients are determined by the sensitivity of a given function to sine and cosine functions. In the inversion method for Fourier series, Fourier coefficients are obtained by measuring the sensitivity of residuals between given functions and test functions (defined as the sum of weighted cosine and sine functions) to cosine and sine functions. The orthogonal property of cosine and sine functions makes the full or approximate Hessian matrix become a diagonal matrix in the inversion for Fourier series. In seismic waveform inversion, the Hessian matrix may or may not be a diagonal matrix, because partial derivative wavefields correlate with each other to some extent, making them semi-orthogonal. At the high-frequency limits, however, the Hessian matrix can be approximated by either a diagonal matrix or a diagonally-dominant matrix. Since we usually deal with relatively low frequencies in seismic waveform inversion, it is not diagonally dominant and thus it is prohibitively expensive to compute the full or approximate Hessian matrix. By interpreting Fourier series with the inversion algorithms, we note that the Fourier series can be computed at an iteration step using any inversion algorithms such as the gradient, full-Newton, and Gauss–Newton methods similar to waveform inversion.

Impedance inversion based on L1 norm regularization

Impedance is an important parameter in geophysics which is greatly relevant to petroleum reservoir. Traditional iterative inversion methods, such as damped least-squares inversion method, often lead to unreliable inversion results when there is strong noise in the seismic records. In this paper, we propose a new impedance inversion method called L1 norm regularization inversion method which has high anti-noise ability. The accuracy of impedances derived from the proposed method is much higher than those derived from damped least-squares inversion method using synthetic seismic record and field data. We set up a forward simulation model and add different intensities of uniform and Gaussian noise into the synthetic seismic records respectively. We find that when there is strong noise in the synthetic record, L1 norm regularization inversion method can obtain the inverted impedance with higher accuracy compared with damped least-squares inversion method. At last, it is found that impedance derived from L1 norm regularization method has a higher accuracy than impedance derived from damped-squares inversion method when there is strong noise both in synthetic record and field data. We conclude that the proposed inversion method shows higher anti-noise ability and can provide impedances with higher accuracy.

Crustal structure of Damavand volcano, Iran, from ambient noise and earthquake tomography

We investigated the shear wave velocity structure beneath the Damavand volcano, Iran, using analysis of Rayleigh wave empirical Green’s functions obtained from the cross-correlation of 1-year seismic ambient noise and 26 local earthquakes. After the observed phase, velocity dispersion curves were inverted to obtain phase velocity tomographic maps; the synthetic phase velocity dispersion curves were then estimated for each evenly spaced geographic grid point. Estimated synthetic phase velocities were inverted using the nonlinear damped least-squares inversion method to obtain a quasi-3D shear wave velocity for the study area. Analysis of obtained quasi-3D shear wave velocity model reveals the presence of three distinct low or high anomalies in the upper 5 km. A low-velocity layer (VS ~ 2.8 km/s) was observed in the upper ~2 km of the crust that includes the sedimentary sequence of carbonate, siliciclastic, and volcanic rocks. According to our results, there is an indication of a high-velocity layer (VS ≥ 3.0 km/s) at the depth range of 2 to 5 km, which may indicate the presence of dense, cooled magma. There is a clear indication of intrusion of the low-velocity anomaly/body into the higher velocity layer, in the depth range of approximately 3.0 to 4.5 km. The corresponding low-velocity body is interpreted as a hot magma chamber associated with the young eruption of Damavand. The chamber is located southwest of the crater, which stretches from the west to beneath the Damavand cone.

Measurement of large negative reactivity of an accelerator-driven system in the Kyoto University Critical Assembly

Large negative reactivity of a subcritical system driven by a pulsed 14 MeV neutron source has been measured in the Kyoto University Critical Assembly. The subcriticality of the accelerator-driven system (ADS) ranged in effective multiplication factor roughly from 0.98 to 0.92, which corresponded to an operational range of an actual ADS proposed by Japan Atomic Energy Agency. As the measurement technique, pulsed neutron method, power spectral analysis for pulsed neutron source, accelerator-beam trip method were employed. From neutron count decay data obtained by the pulsed neutron experiment, not only the prompt-neutron decay constant of fundamental mode but also a higher spatial mode could be derived. The subcriticality was also determined from the fundamental decay constant. The measured cross-power spectral density consisted of a familiar correlated reactor-noise component and many uncorrelated delta-function-like peaks at the integral multiple of pulse repetition frequency. The fundamental prompt-neutron decay constant, i.e., the subcriticality determined from the latter uncorrelated peaks was consistent with that obtained by the above pulsed neutron experiment. However, the magnitude of the former correlated component was reduced with an increase in the subcriticality and eventually this component became almost white at deeply subcritical state ranging in the multiplication factor under 0.95. Consequently, the determination of the decay constant from the correlated component was impossible under such a subcritical state. As data analysis method for the beam trip experiment, both the conventional integral count method and the least-squares inverse kinetics method (LSIKM) were employed. The LSIKM analysis led to the consistent subcriticality with that obtained by the pulsed neutron experiment, while the integral count method significantly underestimated the subcriticality. This underestimation originated from a residual background count, which was maintained after the beam trip. The LSIKM was mostly not influenced by such a slight count rate.

Determination of Subcritical Reactivity of a Thermal Accelerator-Driven System from Beam Trip and Restart Experiment

An experimental technique based on an accelerator-beam trip or restart operation is proposed to determine the subcritical reactivity of an accelerator-driven system (ADS). Applying the least-squares inverse kinetics method to the data analysis, the subcriticality can be inferred from time-sequence neutron count data after these operations. A series ofbeam trip and restart experiments with 14 MeV neutrons were carried out in a thermal ADS of Kyoto University Critical Assembly (KUCA), to demonstrate the applicability of the proposed technique. The subcriticalities evaluated using neutroncounters far from the DT target were consistent with those obtained in a previous pulsed neutron experiment. However, a counter placed close to the target significantly overestimated the subcriticality. The present technique is expected to be available for subcriritcality measurement at startup and shutdown of various ADSs.

An Extended Rod Drop Method Applicable to Subcritical Reactor System Driven by Neutron Source

Control rod drop method based on the integral count technique and the extrapolation one was extended to apply to a subcritical reactor system driven by a neutron source. Similarly to the conventional method, the extended method is capable of rapid estimation of rod reactivity worth and is essential from the educational viewpoint. For the demonstration of the applicability, a series of rod drop experiments was carried out in subcritical and critical states of a reactor system, and rod reactivity worth and effective neutron source strength were determined from the counting data of a neutron detector. In addition, the time-sequence count data were acquired and the least-squares inverse kinetics analysis was made to justify the present extended method. The rod worth and source strength obtained by the present method were consistent with those obtained by the least-squares inverse kinetics method.

An Extended Rod Drop Method Applicable to Subcritical Reactor System Driven by Neutron Source

Control rod drop method based on the integral count technique and the extrapolation one was extended to apply to a subcritical reactor system driven by a neutron source. Similarly to the conventional method, the extended method is capable of rapid estimation of rod reactivity worth and is essential from the educational viewpoint. For the demonstration of the applicability, a series of rod drop experiments was carried out in subcritical and critical states of a reactor system, and rod reactivity worth and effective neutron source strength were determined from the counting data of a neutron detector. In addition, the time-sequence count data were acquired and the least-squares inverse kinetics analysis was made to justify the present extended method. The rod worth and source strength obtained by the present method were consistent with those obtained by the least-squares inverse kinetics method.