Iterative Damped Least-squares Research Articles

Shallow crustal model of the DehDasht in Zagros, Iran, using Rayleigh wave tomography

Surface wave tomography is considered an efficient method to reveal the physical properties of subsurface geological structures, particularly beneath geologically complex areas where conventional active seismic methods are problematic. The salt-rich DehDasht structural basin, in the SW of Iran, is an example of this area, which is characterized by high seismicity. Despite the complex geology and the great thickness of evaporitic deposits in this region, surface wave tomography can retrieve high-resolution images by applying a dense network of stations. Therefore, >11,000 micro-earthquakes (M ≤ 3.0) recorded by 116 broadband stations in 9 months were processed to detect geological structures and a crustal model beneath the DehDasht region. After selecting strict data selection criteria, the observed dispersion curves of Rayleigh waves were then measured using the multiple filter analysis method, and the 2D tomography was done for periods of 0.3 to 5.0 s. A non-linear iterative damped least-squares inversion procedure was then applied to each local dispersion curve to obtain the 1D VS model of the uppermost 5 km, and the results were inserted into the original cells to construct a quasi-3D VS model. Although our resulted VS model has high correlations with known geological units and tectonic features of the study area, we validated it by inverting the observed gravity data independently and consequently recovering the location of the anomalies at the basin and/or around faults, anticlines, and synclines. Generally, the observed low-and high-velocity anomalies are associated with the sedimentary materials.

Analysis of the shallow crustal structure of the Bam Fault Zone, Iran, using the surface wave and first-arrival P-wave tomography approaches

The crustal structure around the Bam fault located southeast of Iran was investigated by applying surface and P-wave tomography to locally recorded aftershocks of the destructive Bam earthquake (Mw 6.6) on December 26, 2003. The aftershock dataset was recorded for 1 month after the main shock by 23 temporary seismic stations operated by the International Institute of Earthquake Engineering and Seismology (IIEES). After applying the selection criteria, 331 well-located events, including more than 1700 waveforms, remained. By calculating the surface wave dispersion measurements in the period range of 0.6–5 s, a tomographic procedure was then performed, and the local dispersion curves were extracted for each geographic grid point (1.8 km × 1.8 km). We applied an iterative damped least-squares inversion procedure to obtain the VSV and VSH models at depths of up to 6.5 km. The discrepancy between these VSV and VSH models also resulted in the radial anisotropy model. Moreover, the P-wave first-arrival times were inverted by the computer program SIMULPS to calculate the VP model for a depth range of 6.5–14 km, while the calculated VP model does not have sufficient resolution at depths of up to 6 km. The tomographic results indicate that a ~ 2 km thick upper sedimentary layer with a shear wave velocity of less than 2.4 km/s and relatively negative radial anisotropy covers the study area (except Bam City, with a thin 500 m positive radial anisotropy layer), which overlay on relatively positive radial anisotropy. Our results indicate that the main shock occurred on the Bam fault escarpment; then the rupturing was shifted to the co-seismic surface rupture (approximately 5 km west of the escarpment) at a depth of ~8–10 km. This dividing of energy on the fault plane (an oblique Y shape) leads to damping most of the energy as a deformation in the lower sedimentary layer.

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.

An iterative damped least-squares algorithm for simultaneously monitoring the development of hemorrhagic and secondary ischemic lesions in brain injuries.

Electrical impedance tomography (EIT) is a non-invasive and real-time imaging method that has the potential to be used for monitoring intracerebral hemorrhage (ICH). Recent studies have proposed that ischemia secondary to ICH occurs simultaneously in the brain. Real-time monitoring of the development of hemorrhage and risk of secondary ischemia is crucial for clinical intervention. However, few studies have explored the performance of EIT monitoring in cases where hemorrhage and secondary ischemia exist. When these lesions get close to each other, or their conductivity and volume changes differ greatly, it becomes challenging for dynamic EIT algorithms to simultaneously reconstruct subtle injuries. To address this, an iterative damped least-squares (IDLS) algorithm is proposed in this study. The quality of the IDLS algorithm was assessed using blur radius and temporal response during computer simulation and a phantom 3D head-shaped model where bidirectional disturbance targets were simulated. The results showed that the IDLS algorithm enhanced contrast and concurrently reconstructed bidirectional disturbance targets in images. Moreover, it showed superior performance in decreasing the blur radius and was time cost-effective. With further improvement, the IDLS algorithm has the potential to be used for monitoring the development of hemorrhage and risk of ischemia secondary to ICH. Graphical abstract (a) and (b) are simulation images of bidirectional disturbance targets with different change ratios of volume (Vr) and conductivity (σr) based on the damped least-squares (DLS) algorithm and iterative damped least-squared (IDLS) algorithm, respectively. (c) shows the performance metrics of blur radius and temporal response with different volume ratio (corresponding to Vr). (d) shows the performance metrics of blur radius and temporal response with different conductivity change percentage (corresponding to σr).

Estimation of fracture weaknesses and integrated attenuation factors from azimuthal variations in seismic amplitudes

Seismic wave propagation in fractured reservoirs exhibits anisotropy and attenuation, which are in turn related to fracture properties (e.g., fracture density) and fluid parameters (e.g., moduli and viscosity). Based on the linear slip theory, stiffness parameters can be determined for fractured and dissipative rocks, from which integrated attenuation factors involving host-rock intrinsic attenuation and fracture-induced attenuation emerge. Using a simplified mathematical form for these stiffness parameters, a linearized mathematical relationship directly relating the reflection coefficient to fracture weaknesses and integrated attenuation factors is available. A two-step inversion approach, involving (1) an iterative damped least-squares algorithm to predict P- and S-wave moduli using seismic angle gathers along the fracture orientation azimuth and (2) an iterative inversion method to estimate fracture weaknesses and integrated attenuation factors from azimuthal amplitude differences, is examined. The objective function for the second step is constructed based on a Bayesian framework. Synthetic testing confirms that fracture weaknesses and integrated attenuation factors are stably determined from seismic amplitudes exhibiting a moderate signal-to-noise ratio. The approach is applied to a field data set from a fractured carbonate reservoir. We observe that geologically reasonable results of fracture weaknesses and integrated attenuation factors are obtained. We conclude that this estimation procedure provides a reliable tool in fracture prediction and inverted attenuation factors appear as additional proofs to identify fluid type in fractured reservoirs.

Rayleigh waves from correlation of seismic noise in Great Island of Tierra del Fuego, Argentina: Constraints on upper crustal structure

In this study, the ambient seismic noise cross-correlation technique is applied to estimate the upper structure of the crust beneath Great Island of Tierra del Fuego (TdF), Argentina, by the analysis of short-period Rayleigh wave group velocities. The island, situated in the southernmost South America, is a key area of investigation among the interaction between the South American and Scotia plates and is considered as a very seismically active one. Through cross-correlating the vertical components of ambient seismic noise registered at four broadband stations in TdF, we were able to extract Rayleigh waves which were used to estimate group velocities in the period band of 2.5–16 s using a time-frequency analysis. Although ambient noise sources are distributed inhomogeneously, robust empirical Green's functions could be recovered from the cross-correlation of 12 months of ambient noise. The observed group velocities were inverted considering a non-linear iterative damped least-squares inversion procedure and several 1-D shear wave velocity models of the upper crust were obtained. According to the inversion results, the S-wave velocity ranges between 1.75 and 3.7 km/s in the first 10 km of crust, depending on the pair of stations considered. These results are in agreement to the major known surface and sub-surface geological and tectonic features known in the area. This study represents the first ambient seismic noise analysis in TdF in order to constraint the upper crust beneath this region. It can also be considered as a successful feasibility study for future analyses with a denser station deployment for a more detailed imaging of structure.

Optimal fourth-order staggered-grid finite-difference scheme for 3D frequency-domain viscoelastic wave modeling

We investigate an optimal fourth-order staggered-grid finite-difference scheme for 3D frequency-domain viscoelastic wave modeling. An anti-lumped mass strategy is incorporated to minimize the numerical dispersion. The optimal finite-difference coefficients and the mass weighting coefficients are obtained by minimizing the misfit between the normalized phase velocities and the unity. An iterative damped least-squares method, the Levenberg–Marquardt algorithm, is utilized for the optimization. Dispersion analysis shows that the optimal fourth-order scheme presents less grid dispersion and anisotropy than the conventional fourth-order scheme with respect to different Poisson's ratios. Moreover, only 3.7 grid-points per minimum shear wavelength are required to keep the error of the group velocities below 1%. The memory cost is then greatly reduced due to a coarser sampling. A parallel iterative method named CARP-CG is used to solve the large ill-conditioned linear system for the frequency-domain modeling. Validations are conducted with respect to both the analytic viscoacoustic and viscoelastic solutions. Compared with the conventional fourth-order scheme, the optimal scheme generates wavefields having smaller error under the same discretization setups. Profiles of the wavefields are presented to confirm better agreement between the optimal results and the analytic solutions.

Source-Type-Specific Inversion of Moment Tensors

The mapping of the fit of seismic moment tensor (MT) solutions in source-type space helps to characterize uncertainty and solution uniqueness. Current practice relies on the forward testing of a distribution of randomly generated MTs in source-type space, which is slow and does not necessarily recover the true maximum fit surface. We design an iterative damped least-squares inversion scheme to invert waveforms and/or P -wave first motions for best-fitting MT solutions for specific source types. An event associated with the sinkhole at the Napoleonville salt dome, Louisiana, an industrial quarry explosion, and an earthquake at The Geysers geothermal field, northern California, are presented as examples. We find that the inversion method is more accurate and successful than the random-search approach in recovering the region of best-fitting MT solutions or source types and is substantially faster. The approach also enables the determination of the best-fitting MT for specified source types such as pure double couples, tensile cracks, or explosions, as well as compound mechanisms in a single numerical framework. Online Material: Description and figures comparing inversion performance using spherical and Cartesian eigenvector parameterizations; figures showing the revised full moment tensor solution of event TE1 and additional network sensitivity solutions of event TE3.

Estimation of shear wave velocity profiles by the inversion of spatial autocorrelation coefficients

The subsurface shear-wave velocity (Vs) is considered to be a key parameter for site characterization and assessment of earthquake hazard because of its great influence on local ground-motion amplification. Array microtremor measurements are widely used for the estimation of shear-wave velocities. Compared to other methods such as frequency-wavenumber (f-k) methods, the spatial autocorrelation (SPAC) method requires fewer sensors and thus is relatively easier to implement and gives robust estimations of shear-wave velocity profiles for depths down to a few hundred meters. The quantity derived from observed data is the SPAC coefficient, which is a function of correlation distance, frequency and phase velocity. Generally, estimation of Vs profiles is a two stage process: Estimation of the dispersion data from the SPAC coefficients and inversion of the dispersion data for shear-wave velocity structure. In this study, instead of inverting dispersion curves, a more practical approach is used; that is, observed SPAC coefficients are directly inverted for the S-wave velocities. A synthetic case and a field data application are presented to test the potential of the inversion algorithm. We obtain an iterative damped least-squares solution with differential smoothing. The differential smoothing approach constrains the change in shear-wave velocities of the adjacent layers and thus stabilizes the inversion.

Spectral Shear-Wave Ground-Motion Scaling in Switzerland

Three-component recordings of local and regional earthquakes and explosions are used to assess the spectral characteristics of attenuation, excitation, and duration of ground motion in Switzerland. The data set consists of 292 events in the magnitude range from 2.0 to 5.2, with a total of 2958 waveforms recorded in Switzerland and the German border region from 1984 to 2000. Distance ranges from 5 to 350 km. Stations are located on National Earthquake Hazard Reduction Program 1994 site classes A and B outcrops. Empirical excitation, site, and attenuation terms are derived for the Fourier spectra and peak ground velocities of the ground motion in the frequency range 1-15 Hz by applying an iterative damped least-squares regression. These results are used to calibrate effective theoretical attenuation and excitation models. A grid search through the parameter space is then applied to obtain the quality factor Q ( f ) and a piecewise linear geometrical spreading function G ( r ), allowing complex behavior of attenuation. Optimum results are obtained for Q ( f ) = 270 f 0.50, G ( r ) = r -1.1 for 0-50 km, r -0.6 for 50-70 km, r +0.2 for 70-100 km, and r -0.5 for distances greater than 100 km. The increased amplitudes in the distance range 70-100 km can be explained by the reflection of shear waves at the Moho. This reflected wave energy also leads to an increased duration of ground motion for the same range. The inverted excitation terms are modeled, based on Brune's source spectrum, applying an effective stress parameter of 5-10 bars and a high-frequency roll off (e- πκ fk) κ = 0.015. The small values of κ agree with the mean site condition. A regional site difference between the Alps and the Alpine Foreland site is obtained with a factor of two higher amplitudes in the latter. Manuscript received 25 August 2001.

Three-Dimensional Velocity Structure of the Upper Crust in the Hida Region, Central Honshu, Japan, and Its Relation to Local Seismicity, Quaternary Active Volcanoes and Faults.

The three-dimensional (3-D) velocity structure of the upper crust in the central part of the Hida region, central Honshu, Japan, has been investigated by simultaneous inversion of travel time data for velocity and hypocentral parameters using Thurber's method (1983). The data used for this purpose were 2, 231 P-wave arrival times from 204 local earthquakes observed at 16 high-sensitivity seismograph stations. The iterative damped least-squares inversion used here provided reliable results with the diagonal elements of resolution matrix well exceeding 0.90 and standard errors less than 3% for the central area. The central part of the Hida region covering an extremely low seismicity area has high velocities probably extending down to a mid-crust. The high-velocity area composed of hard metamorphic and granitic complex extends from the western flank of the Hida mountains to the western section of the seismically active Atotsugawa-Ushikubi faults. The axial part of the Hida mountains and its southwestward area, on the other hand, is covered by remarkable low-velocities. The low velocities may be associated with a high thermal state due to active volcanism beneath the mountains involving volcanoes Mts. Tateyama, Yake, Norikura, and Ontake. The southwestern part of the Hida region is also partially covered by low-velocities beneath the Ryohaku mountains, which may also be related to volcanism of Mt. Hakusan. The obtained 3-D velocity structure appears to be qualitatively consistent with the Bouguer gravity anomalies and also with the large-scale 3-D structure so far obtained.

Digital restoration of blurred images in speckle noise by a two-step procedure

A two-step procedure for restoring blurred images in speckle noise is described. The restoration method includes a local-statistics algorithm for smoothing speckle noise and an iterative damped least-squares algorithm for removal of the blur. We present some experimental results that demonstrate the effectiveness of the proposed procedure.

Restoration of bandlimited images by an iterative damped least-squares method with adaptive regularization

A digital method for restoring bandlimited images in the presence of noise is described. The restoration procedure is an iterative damped least-squares algorithm that is based on the principle of damped least squares. This method has the advantage of tolerance to noise by incorporating additional constraints of non-negativity of the object and the adaptive regularization as well as a priori information about the finite extent of the object. We present some results of computer simulations of both 1-D and 2-D images that demonstrate the effectiveness of the proposed method.

Image restoration by an iterative damped least-squares method with non-negativity constraint

In this paper we describe a new digital method for restoring linearly degraded images in the presence of noise. The restoration procedure is an iterative damped least-squares (DLS) algorithm that is based on the principle of damped least squares. This method acquires the advantage of tolerance to noise by incorporating additional constraints of non-negativity of the object and the adaptive technique. We also use a median filter in the proposed algorithm to suppress spiky noise. After discussing the convergence of the iterative DLS algorithm, we present some results of computer simulations that demonstrate the effectiveness of the proposed method.

P-velocity models of the southern Diablo Range, California, from inversion of earthquake and explosion arrival times

Abstract Velocity models are developed for the region in and around the southern Diablo Range, lying just east of seismically active portions of the San Andreas and Calaveras faults in central California. An iterative damped least-squares inversion procedure is used to simultaneously estimate hypocentral parameters, station delays, and velocities in a horizontally layered crustal model. Arrival times from both earthquakes and explosions of known location and origin time are used as data. The procedure used has the potential to detect low-velocity layers if they exist. The models obtained show a rapid increase in velocity to about 5.5 km/sec at 3 km depth, another rapid increase to about 6.0 km/sec at 6 km depth, and a roughly linear increase below 6 km with the velocity reaching approximately 6.7 km/sec at a depth of 15 km.