Search

Accurately injecting and extrapolating multicomponent seismic data is a major step in elastic reverse time migration (RTM) and full-waveform inversion. The theoretical basis of this process is the seismic representation theorem. We propose a numerical approach based on the collocated-grid MacCormack finite-difference (FD) scheme to implement the injection and extrapolation. The MacCormack FD scheme defines wavefield parameters on the same grid. When dealing with the boundary condition in the representation theorem, it is more advantageous than the commonly used staggered-grid FD scheme which defines the stress and velocity on offset grids; thus, interpolations are required to link the boundary values to nearby nodes. In contrast, with the MacCormack scheme, the traction and particle velocity data can be directly injected into the boundary nodes as equivalent concentrated forces, force dipoles, and double couples. The accuracy of this approach is validated by comparing the regenerated wavefield with the original wavefield. Our primary purpose is to develop an injection method for elastic RTM, but this approach also can be used to numerically investigate other interesting issues in the elastic wave injection, including errors and artifacts caused by injection from unclosed boundaries, incomplete boundary values, and uncertain velocity models. In particular, seismic data are usually acquired on the ground surface where the traction-free condition is naturally maintained. We simulate a land acquisition scenario by generating a synthetic multicomponent data set on a traction-free surface. The reconstructed wavefields from velocity data are compared with the original forward-propagated waves to evaluate the accuracy of the boundary treatments. Two multicomponent elastic RTM examples are presented to demonstrate the applications of the proposed injection method. By combining the regenerated receiver-side wavefield and the vectorized imaging condition, the PP, PS, SP, and SS images are calculated.

SUMMERYDue to the far-field effect of the India–Eurasia collision, the Tien Shan orogenic belt has been undergoing reactivation and modification. Two end-member models of the geodynamic mechanisms are (1) surface uplift due to crustal shortening caused by lithospheric compression and (2) mountain formation resulting from thermal upwelling of asthenospheric mantle materials generated by lithospheric subduction. However, the topography along the Tien Shan orogenic belt changes significantly, and the deep structure and dynamic process are quite different beneath the Tien Shan orogenic belt from both geological and geophysical observations. Therefore, the reactivation and modification of the Tien Shan orogenic belt are likely influenced by both geodynamic mechanisms, which also generate various thermal anomalies in the crust. Seismic Lg-wave attenuation is very sensitive to crustal material composition and status and can point to the presence of partial melting within the crust resulting from mantle upwelling. In this study, we develop a high-resolution Lg-wave attenuation model between 0.05 and 10.0 Hz in Northwest China and use lateral attenuation variations to infer thermal structures in the crust. The central Tien Shan is characterized by prominent low-QLg anomalies, whereas relatively high-QLg distributions are imaged beneath the eastern and western Tien Shan. The surface uplift and crustal deformation are mostly related to the upwelling of hot mantle materials in the central Tien Shan and are likely induced by lithospheric compression in the eastern and western Tien Shan. However, low-Q anomalies are observed in the junction between the Pamir Plateau and western Tien Shan, indicating that the uplift in the south of the western Tien Shan is related to thermal subduction-induced upwelling and intrusion into the crust due to the collision between the Indian and Eurasian plates. The Kazakh Shield, characterized by pronounced high-QLg values, is likely a cold and hard terrane, and hence blocks the far-field effect of the India–Eurasia collision.

Background. Myasthenia gravis (MG) is an acquired autoimmune disease. The main clinical features of MG are skeletal muscle fatigue and pathological fatigue, which worsen at night or after fatigue, such as dyspnea, dysphagia, and systemic weakness. Plasma exchange (PE) is often used in patients with acute exacerbation of MG. Intravenous immunoglobulin (IVIG) is a collection of immunoglobulins from thousands of donors. IVIG can replace a variety of immunosuppressants or PE. However, the effect of PE or IVIG on patients’ consciousness, immune function, and prognosis is not clear. Objective. A prospective randomized test of the effects of PE combined with immunoglobulin on consciousness, immune function, and prognosis in patients with myasthenia gravis crisis (MGC). Methods. Sixty patients with MGC treated from February 2019 to April 2021 were enrolled in our hospital. The cases who received PE were set as the PE group, and those who received PE combined with immunoglobulin were set as the PE+immunoglobulin group. The efficacy, clinical score, state of consciousness, immune function, acetylcholine receptor antibody (AChR-Ab), lymphocyte (LYM), albumin (ALB) levels, and the incidence of adverse reactions were compared. Results. The improvement rate was 100.005% in the treatment group and 83.33% in the PE group. After treatment, the clinical score of the PE+immunoglobulin group was lower than that of the PE group, and the clinical relative score of the PE+immunoglobulin group was higher than that of the PE group ( P < 0.05 ). The number of conscious people in the PE+immunoglobulin group was more than that in the PE group ( P < 0.05 ). Immunoglobulin A, immunoglobulin M, immunoglobulin G, and immunoglobulin G in the PE+immunoglobulin group were higher than those in the PE group ( P < 0.05 ). The levels of AChR-Ab and ALB in the PE+immunoglobulin group were higher than those in the PE group, while the level of LYM in the PE+immunoglobulin group was lower than that in the PE group. The incidence of skin system, gastrointestinal system, nervous system, and systemic damage in the PE+immunoglobulin group was lower than that in the PE group ( P < 0.05 ). Conclusion. The treatment of MGC with PE combined with immunoglobulin can not only effectively enhance the consciousness and immune function of patients but also effectively promote the prognosis, and the safety of treatment can be guaranteed.

We investigate the 4 August 2020 Beirut accidental chemical explosion based on regional seismic data recorded at both on- and off-shore stations. The Lg-wave body-wave magnitude is mb(Lg)=3.30±0.46 for the Beirut explosion. The explosive yield obtained using an empirical magnitude-yield relation based on a fully buried explosion source model is only 0.112 kt. Alternatively, the yield estimated using an empirical relation between the yield and crater size is 1.22 kt, with the uncertainties between 0.48 and 2.3 kt. The latter is closer to reality. The P/S spectral amplitude ratios, including Pg/Lg, Pn/Lg, and Pn/Sn, are calculated for the Beirut explosion and nearby natural earthquakes. We find the P/S spectral ratios are effective in discriminating the explosion from earthquakes in the Northwestern Arabia plate. By comparing the spectral ratios of large open-pit explosions, including the Beirut, Xiangshui, and Tianjin explosions, with those from historical nuclear explosions, buried small chemical explosions, and natural earthquakes, we further investigate the detailed differences of network-averaged P/S spectral ratios between different source types.

During hydraulic fracturing (HF) stimulation for unconventional reservoir development, seismic attenuation has a significant influence on high-frequency microseismic data. Attenuation also provides important information for characterizing reservoir structure and changes to it due to HF injections. However, the attenuation effect is typically not considered in microseismic analysis. We have adopted the spectral ratio and centroid-frequency shift methods to estimate the subsurface attenuation (the factor [Formula: see text]) in a tight oil reservoir in the Ordos Basin, China. The P- and S-wave attenuations are calculated using the 3C waveform data recorded by a single-well downhole geophone array during a 12-stage HF stimulation. Both methods provide similar results (with differences in [Formula: see text] of absolute values less than 0.010 for P- and S-waves). For individual events, their median [Formula: see text] values calculated from different geophones are selected to represent the average attenuation. Spatiotemporal variations in attenuation are obtained by investigating [Formula: see text] values along propagating rays linking different source–receiver pairs. The [Formula: see text] values derived at different HF stages reveal significant attenuation in the targeted tight sandstone layer (0.030–0.062 for [Formula: see text] and 0.026−0.058 for [Formula: see text]), and the attenuation is apparently increased by fluid injection activities. We explain the sudden decrease in attenuation near the geophone array as a result of high shale content using log data from a horizontal treatment well. The consistency between the [Formula: see text] values and horizontal well-log data, as well as the HF process, indicates the reliability and robustness of the attenuation results. By studying spatiotemporal variations in attenuation, the changes in subsurface structures may be quantitatively characterized, thereby creating a reliable basis for microseismic modeling and data processing and providing additional information on monitoring the HF process.

SUMMARY On 2019 March 21, an explosion accidentally occurred at a chemical plant in Xiangshui, Yancheng City, Jiangsu Province, China. Using broad-band digital seismic data from East China, South Korea and Japan, we investigate properties of the Xiangshui explosion as well as two nearby chemical explosions and four nearby natural earthquakes in Jiangsu Province, East China. From Lg and Rayleigh waves recorded by regional networks, both body wave magnitude mb (Lg) and surface wave magnitude Ms (Rayleigh) are calculated for these events. The magnitudes of the Xiangshui explosion are mb (Lg) = 3.39 ± 0.24 and Ms = 1.95 ± 0.27, respectively. Both the empirical magnitude–yield relation for buried explosion and empirical yield–crater dimension relation for open-pit explosion are adopted for investigating the explosive yield. The result from the yield–crater dimension relation is approximately 492 ton, which is consistent with the ground truth and considerably larger than that from the buried source model. This also reveals that, for Xiangshui explosion, the explosion to seismic energy conversion rate is approximately one-third compared to a similar sized fully confined explosion. By comparing the body wave and surface wave magnitudes from explosions and nearby earthquakes, we find that the mb:Ms discriminant calculated at regional distances cannot properly distinguish explosions from natural earthquakes. However, the P/S spectral ratios Pg/Lg, Pn/Lg and Pn/Sn from the same data set can be good discriminants for identifying explosions from earthquakes.

SUMMARYThe crustal attenuation structure can effectively reveal the rheology and thermal properties of different geological blocks, and can provide seismological constraints on regional tectonic evolution. Based on 11 306 vertical-component Lg-wave seismograms recorded at 111 broad-band stations from 891 crustal earthquakes that occurred between 1994 and 2020, a broad-band Lg-wave attenuation model is obtained for Southeast Asia. This study demonstrates the capability of applying crustal Lg-wave attenuation inversion in a complex region mixed with continents, islands and marginal seas. The resolution approaches 2° in most parts in the study region. Lg blockages are observed at places with sharp Moho depth changes. The resultant Q models are consistent with regional geologic structures provided by previous studies. Prominent low-attenuation anomalies are located in the Sundaland Core containing stable ancient crust, including Indochina, Malay Peninsula, East Sumatra, Sunda Shelf and Borneo Core. Regions with strong attenuation are associated with complex tectonic conditions, such as the Indo–Australian subduction zone, sutures in Sarawak and Sabah. The observed Lg-wave attenuation characteristics provide constraints on the tectonic affinities and evolutions of the geological blocks. The results show that the Borneo Core remained stable since its accretion with the Sundaland Core. Ancient blocks are characterized by weak Lg attenuation, whereas geologically younger terranes are often characterized by strong Lg attenuation, which can be exploited to better understand the separation and convergence of plates during the tectonic evolution.

We have developed an optimal method to determine expansion parameters for flexible stencils in 2D scalar-wave finite-difference frequency-domain (FDFD) simulation. Our stencil only requires the involved grid points to be paired and rotationally symmetric around the central point. We apply this method to the transition zone in discontinuous-grid modeling, in which the key issue is designing particular FDFD stencils to correctly propagate the wavefield passing through the discontinuous interface. Our method can work in an FDFD discontinuous grid with arbitrary integer coarse- to fine-grid spacing ratios. Numerical examples are developed to determine how to apply this optimal method to discontinuous-grid FDFD schemes with spacing ratios of 3 and 5. The synthetic wavefields are highly consistent to those calculated using the conventional dense uniform grid, and the memory requirement and computational costs are greatly reduced. For velocity models with large contrasts, our discontinuous-grid FDFD method can significantly improve the computational efficiency in forward modeling, imaging, and full-waveform inversion.

Through the study of microseismic focal mechanisms, information such as fracture orientation, event magnitude, and in situ stress status can be quantitatively obtained, thus providing a reliable basis for unconventional oil and gas exploration. Most source inversion methods assume that the medium is isotropic. However, hydraulic fracturing is usually conducted in sedimentary rocks, which often exhibit strong anisotropy. Neglecting this anisotropy may cause errors in focal mechanism inversion results. We have developed a microseismic focal mechanism inversion method that considers velocity anisotropy in a vertically transverse isotropic medium. To generate synthetic data, we adopt the moment-tensor model to represent microearthquake sources. We use a staggered-grid finite-difference method to calculate synthetic seismograms in anisotropic media. We perform seismic moment-tensor (SMT) inversion with only P-waves by matching the synthetic and observed waveforms. The synthetic and field data sets are used to test the inversion method. For the field data set, we investigate the inversion stability using randomly selected partial data sets in the calculation. We pay special attention to analyze the sensitivity of the inversion. We test and evaluate the impact of noise in the data and errors in the model parameters ([Formula: see text], [Formula: see text], and [Formula: see text]) on the SMT inversion using synthetic data sets. The results indicate that, for a surface acquisition system, our method can tolerate moderate noise in the data and deviations in the anisotropy parameters can cause errors in the SMT inversion, especially for dip-slip events and the inverted percentages of non-double-couple (DC) components. According to our study, including anisotropy in the model is important to obtain reliable non-DC components of moment tensors for hydraulic fracturing-induced microearthquakes.

The applicability of the empirical magnitude–yield relations developed for northeast China and Korean Peninsula explosions was investigated for data from northwest China. We collected regional broadband digital seismic data from 13 chemical explosions (CEx) detonated between 6 September and 10 October 2018, on the eastern margin of the Junggar basin, northwest China, five nuclear tests at the Semipalatinsk nuclear test site, and eight natural earthquakes. Both Lg and Rayleigh-wave magnitudes (mb(Lg) and Ms, respectively) were estimated for these events. Similar to the North Korean test site, the mb(Lg)–Ms discriminant did not properly distinguish explosions from natural earthquakes at the Semipalatinsk test site. However, network-averaged P/S spectral ratios (Pg/Lg, Pn/Lg, and Pn/Sn) did successfully discriminate explosions from earthquakes at both the North Korean and the Semipalatinsk test sites at frequencies above 2.0 Hz. Based on 13 known-yield CEx, we selected an empirical magnitude–yield relation to constrain the explosive yields of five historical nuclear tests at the Semipalatinsk test site. The resulting yields are lower than those previously obtained from teleseismic observations.

SUMMARY Wave-equation-based traveltime tomography has been extensively applied in both global tomography and seismic exploration. Typically, the traveltime Fréchet derivative is obtained using the first-order Born approximation, which is only satisfied for weak velocity perturbations and small phase shifts (i.e. the weak-scattering assumption). Although the small phase-shift restriction can be handled with the Rytov approximation, the weak velocity-perturbation assumption is still a major limitation. The recently developed generalized Rytov approximation (GRA) method can achieve an improved phase accuracy of the forward-scattered wavefield, in the presence of large-scale and strong velocity perturbations. In this paper, we combine GRA with the classical finite-frequency theory and propose a GRA-based traveltime sensitivity kernel (GRA-TSK), which overcomes the weak-scattering limitation of the conventional finite-frequency methods. Numerical examples demonstrate that the accumulated time delay of forward-scattered waves caused by large-scale smooth perturbations can be correctly handled by the GRA-TSK, regardless of the magnitude of the velocity perturbations. Then, we apply the new sensitivity kernel to solve the traveltime inverse problem, and we propose a matrix-free Gauss–Newton method that has a faster convergence rate compared with the gradient-based method. Numerical tests show that, compared with the conventional adjoint traveltime tomography, the proposed GRA-based traveltime tomography can obtain a more accurate model with a faster convergence rate, making it more suited for recovering the large-intermediate scale of the velocity model, even for strong-perturbation and complex subsurface structures.

Many full-waveform inversion schemes are based on the iterative perturbation theory to fit the observed waveforms. When the observed waveforms lack low frequencies, those schemes may encounter convergence problems due to cycle skipping when the initial velocity model is far from the true model. To mitigate this difficulty, we have developed a new objective function that fits the seismic-waveform intensity, so the dependence of the starting model can be reduced. The waveform intensity is proportional to the square of its amplitude. Forming the intensity using the waveform is a nonlinear operation, which separates the original waveform spectrum into an ultra-low-frequency part and a higher frequency part, even for data that originally do not have low-frequency contents. Therefore, conducting multiscale inversions starting from ultra-low-frequency intensity data can largely avoid the cycle-skipping problem. We formulate the intensity objective function, the minimization process, and the gradient. Using numerical examples, we determine that the proposed method was very promising and could invert for the model using data lacking low-frequency information.

The discontinuous-grid method can greatly reduce the storage requirement and computational cost in finite-difference modeling, especially for models with large velocity contrasts. However, this technique is mostly applied to time-domain methods. We have developed a discontinuous-grid finite-difference scheme for frequency-domain 2D scalar wave modeling. Special frequency-domain finite-difference stencils are designed in the fine-coarse grid transition zone. The coarse-to-fine-grid spacing ratio is restricted to [Formula: see text], where [Formula: see text] is a positive integer. Optimization equations are formulated to obtain expansion coefficients for irregular stencils in the transition zone. The proposed method works well when teamed with commonly used 9- and 25-point schemes. Compared with the conventional frequency-domain finite-difference method, the proposed discontinuous-grid method can largely reduce the size of the impedance matrix and number of nonzero elements. Numerical experiments demonstrated that the proposed discontinuous-grid scheme can significantly reduce memory and computational costs, while still yielding almost identical results compared with those from conventional uniform-grid simulations. When tested for a very long elapsed time, the frequency-domain discontinuous-grid method does not show instability problems as its time-domain counterpart usually does.

The seismological characteristics of the 15 February 2013 Chelyabinsk bolide explosion are investigated based on seismograms recorded at 50 stations with epicentral distances ranging from 229 to 4324 km. By using 8–25 s vertical-component Rayleigh waveforms, we obtain a surface-wave magnitude of 4.17±0.31 for this event. According to the relationship among the Rayleigh-wave magnitude, burst height and explosive yield, the explosion yield is estimated to be 686 kt. Using a single-force source to fit the observed Rayleigh waveforms, we obtain a single force of 1.03×1012 N, which is equivalent to the impact from the shock wave generated by the bolide explosion.

We propose a method based on the Poynting vector that combines angle-domain imaging and image amplitude correction to overcome the shortcomings of reverse-time migration that cannot handle different angles during wave propagation. First, the local image matrix (LIM) and local illumination matrix are constructed, and the wavefield propagation directions are decomposed. The angle-domain imaging conditions are established in the local imaging matrix to remove low-wavenumber artifacts. Next, the angle-domain common image gathers are extracted and the dip angle is calculated, and the amplitude-corrected factors in the dip angle domain are calculated. The partial images are corrected by factors corresponding to the different angles and then are superimposed to perform the amplitude correction of the final image. Angle-domain imaging based on the Poynting vector improves the computation efficiency compared with local plane-wave decomposition. Finally, numerical simulations based on the SEG/EAGE velocity model are used to validate the proposed method.

At 11:30 on September 3, 2017 (Beijing time), a strong earthquake occurred in Democratic People′s Republic of Korea (DPRK). International seismic monitoring agencies, e.g., the China Earthquake Network Center (CENC) and the United States Geological Survey, suspected that it is an explosion. Based on seismic data from the China National Digital Seismic Network (CNDSN) and Global Seismic Network (GSN), we investigated characteristics of this event and an aftershock 8 min after the main event. The P- and S-wave excitation functions of explosion and earthquake sources are scaled differently. Therefore, the P/S-type spectral ratios can be an effective discriminant for separating explosions from earthquakes. Using the P/S spectral ratios Pg/Lg, Pn/Lg and Pn/Sn as discriminants, we confirmed the 3 September 2017 event was an explosion. For the aftershock occurred after the main event, we identified it is a collapse, likely caused by the failure of the explosion generated cavity. Using a pre-calibrated regional seismic network in Northeast China and the Korean peninsula, and the regional Lg-wave attenuation model developed previously, we obtained the Lg wave body wave magnitudes for the 3 September 2017 main event and its aftershock to be m b(Lg)=5.6±0.2, and 3.95±0.04. We used a group of historical events to calibrate the regional network for calculating Rayleigh wave magnitude. After correcting for the site responses, the network averaged surface wave magnitudes for the main shock and the aftershock were obtained to be M s=5.1±0.2 and 3.95±0.08. To test body-wave versus surface-wave magnitude as a potential discriminant, we compared the M s (Rayleigh) and m b (Lg) for all 6 DPRK nuclear explosions and a group of earthquakes in Northeast China and the Korean peninsula. The explosion and earthquake populations were largely overlapped with each other. The above results show that the P/S ratio method is a more effective discriminant than the m b(Lg)- M s criterion in Northeast China and the Korean Peninsula. The seismic yield of an underground nuclear explosion can be estimated from its magnitude using a calibrated empirical magnitude-yield relation. However, the DPRK test site (DPRKTS) is an uncalibrated test site. Considering that the DPRKTS is located at a granite site in a stable geology platform, we adopted the fully-coupled hard-rock site equation used at the Novaya Zemlya test site to calculate the yield at the DPRKTS. The estimated yield for the 3 September 2017 explosion was 56 kt using this relationship and assuming a normally scaled burial depth. Transferring the measurement error of ±0.2 magnitude unit to the yield calculation introduced uncertainties between 30 and 100 kt. However, if the explosion was over buried at depths between 1000 and 2400 m, the yield could be increased to 100–200 kt.

For elastic-wave propagation in strong-contrast media, we developed the theory and method of a new one-way operator, which can overcome the weak perturbation assumption of a traditional elastic thin-slab propagator. In the framework of the new propagator, the elastic model is separated into several domains along the strong-contrast boundaries. Different background velocities are chosen for different domains so that each domain can be considered as a weakly heterogeneous media with a sharp boundary. The new propagator adopts a hybrid approach to extrapolate the wavefields in each domain: a traditional thin-slab propagator for weak heterogeneities combined with a reflection/transmission operator for wave exchanging at the sharp boundary. It inherits the dual domain (wavenumber-space domain), depth-marching implementation of an elastic thin-slab propagator. The boundary is discretized into a collection of boundary segments based on the marching depth, and the reflection/transmission operators are estimated locally at the boundary segment by making a tangent plane approximation. By neglecting the reverberation between the boundary segments, the reflected/transmitted wave from the boundary can be gradually incorporated into wave propagation in a one-way fashion. We validated the accuracy and efficiency of the hybrid propagator by showing numerical simulation results on two simple salt models. Finally, we applied the hybrid propagator to subsalt imaging using a synthetic data set generated from a 2D subsalt model.

We characterize the seismic events that occurred in North Korea on 9 September 2016 and South Korea on 12 September 2016. The 9 September 2016 event was identified as an explosion, and the two 12 September 2016 events were identified as natural earthquakes using the P/S (P-and S-wave) spectral ratios, Pg/Lg, Pn/Lg, and Pn/Sn as discriminants. The explosive event was relocated within the North Korean nuclear test site using a relative location method and the 2006 North Korea underground nuclear test as the master event, and the epicenter was identified at 41.2976 degrees N latitude and 129.0804 degrees E longitude. From the regional Lg and Rayleigh waves, the body-and surface-wave magnitudes for the 9 September 2016 event were calculated as m(b) (Lg) = 4: 8 +/- 0: 2 and M-s = 4: 2 +/- 0: 1. By adopting an empirical magnitude-yield relation for the body-wave magnitude, and assuming that the explosion was fully coupled and detonated at a normally scaled depth, we estimated that the seismic yield was similar to 6 kt, and the uncertainty range was between 3 and 11 kt. If an overburied depth range between 780 and 1200 m was applied, then the yield would be increased to 16-22 kt.

To mitigate the local minima problem in full waveform inversion, the iteration process can be carried out in a multi-scale way, with the inversion starting from large-scale perturbations and gradually changing to small-scale structures. However, this technique is often limited due to the lack of low-frequency information in the observed data. In a scattering process, the scale length of the velocity perturbation is not only related to the frequency, but also associated with the scattering angles. We introduce a scattering angle filter in frequency-domain full waveform inversion. Choosing proper scattering angles in the iteration process can control the wavenumber components entering into the velocity updating, and accomplish multi-scale waveform inversion. Using the slant-stacking method, we decompose the source- and receiver-side waves into local plane waves, from which the scattering angles between the incident and scattering directions are calculated and used as the constraint to determine the wavenumber components to be retrieved. Small scattering angles are related to large-scale model perturbations. Therefore, by filtering out large scattering angles, we can invert the large-scale background perturbations first. Then, by gradually increasing the band pass of scattering angles, we can retrieve fine-scale structures. Numerical examples demonstrate that when the initial model has large velocity errors and the low-frequency information is unavailable in the data, the multi-scale strategy based on the angle-domain wavenumber filter can largely improve the convergence in the initial iteration stage and make the searching towards the global minimum.

We have developed a general optimal method for 2D frequency-domain finite-difference simulation of the scalar wave equation. For a given finite-difference stencil, this method can generate the dispersion equation and optimize the expansion coefficients. Many commonly used frequency-domain finite-difference schemes (e.g., grids with different numbers of points, rotated grids, and grid spaces with different aspect ratios) can be derived as special cases under this framework. The possibility of expanding this method to 3D does exist. Based on the 2D scalar wave equation, the optimized coefficients of 25-point, 9-point, 17-point, and 15-point schemes have been worked out. The dispersion analysis indicates that our 25-point scheme has much higher accuracy than the average-derivative method 25-point scheme. The number of grid points per the smallest wavelength is reduced from 2.78 to 2.13 for a maximum phase velocity errors of 1%. The synthetic seismograms and the wavefield snapshots calculated using our optimal 25-point finite-different scheme give smaller dispersions than other finite-difference schemes.

Summary Conventional migrations are apt to structure imaging, but often incapable of generating true-amplitude subsurface image. We propose a true-amplitude Gaussian beam migration (GBM) method under the framework of seismic illumination compensation. A novel scheme based on the GBM is developed to estimate the point spread functions, with which the illumination compensation can be efficiently implemented in the local wave-number domain. The total computational cost of the proposed true-amplitude imaging includes one Born modelling process and two conventional GBM processes, which is more efficient than the true-amplitude imaging using least squares migration which requires multiple iterations. Numerical examples using synthetic data demonstrated the effectiveness and efficiency of the proposed method.

The Pn phase for crustal earthquakes refracts through the uppermost mantle to become the first seismic‐wave arrival at distances of ∼200 to ∼1500 km. The amplitude of Pn is particularly sensitive to the mantle lid radial velocity gradient, and the elastic geometric spreading is frequency dependent. Lateral variations in the radial velocity gradient and in the crustal thickness encountered by Pn phases that traverse boundaries of geological provinces further complicate the frequency‐dependent geometric spreading behavior. We use 2D finite‐difference (FD) calculations to explore the effects of smoothly varying lateral structures in the lid and crust on Pn phases in the 0.4–10 Hz frequency range. Pn geometric spreading for 1D structures with mantle lid gradients is calculated first to ensure that the FD calculations correctly capture the expected frequency dependence. Then, 2D models with smooth laterally varying radial lid gradients are computed to explore how the phase samples the structures to give the overall geometric spreading. The direction of propagation through the laterally varying structure affects the frequency dependence and the transition from source‐side to receiver‐side dominance in the spreading. Laterally varying crustal thickness is explored in the context of the North Korean nuclear test site, for which many observations of Pn are made by stations with paths either through continental structure or along paths with variation from continental to oceanic to arc structure (in Japan). Frequency‐dependent amplitude effects are observed on the paths through oceanic structure. 2D FD calculations for structures with a central region of thin oceanic crust reproduce the observed first‐order frequency dependence of Pn along corresponding paths, including rapid amplitude decay for frequencies above 6 Hz. These analyses demonstrate that even though Pn propagation is very complex, modeling approaches can guide interpretations of the behavior when constraints on the structure are available. Electronic Supplement: Figures of Pn ‐wave snapshots, first Pn arrival times, and frequency dependence of the Pn amplitudes.

Geometric spreading of the Pn phase for shallow sources, which refracts through the uppermost mantle and is the first seismic‐wave arrival at distances of ∼200 to ∼1500 km, is expected to be frequency dependent for most velocity structures. The smoothly varying distance and frequency dependence of Pn geometric spreading predicted for simple 1D spherical structures are dependent upon precise interference of multiple underside reflections from the Moho boundary. This interference is strongly impacted by lateral variations in the velocity structure. Prior work has shown that the presence of 2D random heterogeneities in the mantle lid or topography of the Moho boundary modifies Pn geometric spreading for ∼1 Hz signals, suppressing the distinctive character of the distance dependence for a 1D structure and approaching a power‐law behavior as the root mean square (rms) strength of the heterogeneity increases. Here, 2D finite‐difference calculations of the effects of random heterogeneities in the mantle lid on the Pn geometric spreading are extended to frequencies of up to 10 Hz, quantifying the frequency dependence. Although the shape of the distance dependence is still systematically modified, the strength of the frequency dependence is actually increased from that for homogeneous models for many suites of models with varying rms perturbations and correlation lengths of exponential random heterogeneity distributions superimposed on constant and linear gradient background mantle lid velocity structures. This indicates that assumption of frequency‐independent power‐law spreading for Pn directly introduces artificial frequency dependence to inferred attenuation coefficients. Similar behavior is expected for Sn phases. Efforts to characterize the mantle lid heterogeneity spectrum are required to overcome this trade‐off with anelasticity representations. Electronic Supplement: Figures showing calculated Pn ‐wave spectral amplitudes versus distance.

Illumination and resolution analysis provides vital information regarding the response of an imaging system to subsurface structures. However, generating the resolution function is often computationally intensive, which prevents it from being widely used in practice. This problem is particularly severe for the time-domain migration method, such as reverse time migration (RTM). To solve this problem, we have developed a fast full-wave-based illumination and resolution analysis method. The source- and receiver-side waves are extrapolated to the subsurface for simulating the imaging process. To create the relations linking the incident and scattering directions to the target wavenumber components, we have adopted an efficient Poynting-vector-based method for wavefield angle decomposition. By taking the approximation that the time-domain wavefield preserves the source spectrum during propagation, massive input/output and trace-by-trace Fourier transform can be avoided and the finite-frequency calculated broadband signal can be directly converted to the wavenumber domain for calculating the point spreading functions (PSFs). Combining these approaches, the resulting method avoids intensive calculations, massive input/output, and huge storage requirements commonly involved in generating the illumination and resolution functions. The method is highly efficient and particularly suitable to team with RTM for resolution analysis. Numerical examples are used to validate this method. We have determined how to calculate the acquisition dip response and PSFs, based on which, the quality of the depth image can be significantly improved.

Seismic migration moves reflections to their true subsurface positions and yields seismic images of subsurface areas. However, due to limited acquisition aperture, complex overburden structure and target dipping angle, the migration often generates a distorted image of the actual subsurface structure. Seismic illumination and resolution analyses provide a quantitative description of how the above-mentioned factors distort the image. The point spread function (PSF) gives the resolution of the depth image and carries full information about the factors affecting the quality of the image. The staining algorithm establishes a correspondence between a certain structure and its relevant wavefield and reflected data. In this paper, we use the staining algorithm to calculate the PSFs, then use these PSFs for extracting the acquisition dip response and correcting the original depth image by deconvolution. We present relevant results of the SEG salt model. The staining algorithm provides an efficient tool for calculating the PSF and for conducting broadband seismic illumination and resolution analyses.

Seismology plays an important role in characterizing potential underground nuclear tests. Using broad-band digital seismic data from Northeast China, South Korea and Japan, we investigated the properties of the recent seismic event occurred in North Korea on 2016 January 6. Using a relative location method and choosing the previous 2006 explosion as the master event, the 2016 event was located within the North Korean nuclear test site, with its epicentre at latitude 41.3003°N and longitude 129.0678°E, approximately 900 m north and 500 m west of the previous event on 2013 February 12. Based on the error ellipse, the relocation uncertainty was approximately 70 m. Using the P/S spectral ratios, including Pg/Lg, Pn/Lg and Pn/Sn, as the discriminants, we identify the 2016 event as an explosion rather than an earthquake. The body-wave magnitude calculated from regional wave Lg is mb(Lg) equal to 4.7 ± 0.2. Adopting an empirical magnitude–yield relation, and assuming that the explosion is fully coupled and detonated at a normally scaled depth, we find that the seismic yield is about 4 kt, with the uncertainties allowing a range from 2 to 8 kt.

We investigate the geometric spreading and attenuation of seismic Pn waves in Northeast China and the Korean Peninsula. A high-quality broadband Pn wave data set generated by North Korean nuclear tests is used to constrain the parameters of a frequency-dependent log-quadratic geometric spreading function and a power law Pn Q model. The geometric spreading function and apparent Pn wave Q are obtained for Northeast China and the Korean Peninsula between 2.0 and 10.0 Hz. Using the two-station amplitude ratios of the Pn spectra and correcting them with the known spreading function, we remove the contributions of the source and crust from the apparent Pn Q and retrieve the P wave attenuation information along the pure upper mantle path. We then use both Pn amplitudes and amplitude ratios in a tomographic approach to obtain the upper mantle P wave attenuation in the studied area. The Pn wave spectra observed in China are compared with those recorded in Japan, and the result reveals that the high-frequency Pn signal across the oceanic path attenuated faster compared with those through the continental path.

Internal multiple reflections have different propagation paths than primary reflections and surface-related multiples; thus, they can complement illumination where other waves are unavailable. Consequently, migration of internal multiples may provide vital images in the shadow zone of primary reflections. We have developed an approach to image structures with internal multiples. With this method, internal multiples can be separated and used as data for reverse time migration. Instead of applying the imaging condition directly to the source and receiver wavefields, we decomposed the extrapolated source and receiver waves into up- and downgoing waves. Then we applied the imaging condition to upgoing source and receiver waves to form an up-up (UU) image and also to downgoing source and receiver waves to form a down-down (DD) image. The image of the true reflector and the artifacts behave differently in the DD and UU images. A similarity analysis is conducted to the UU and DD images to separate true images from artifacts. Numerical examples with simple velocity models are used to demonstrate how to construct UU and DD imaging conditions and eliminate migration artifacts. Finally, we test the 2D SMAART JV model with salt structures for migration using internal multiples.

The first-order (or linear) Rytov or Born approximation is the foundation for formulation of wave-equation tomography and waveform inversion, so the validity of the Rytov/Born approximation can substantially affect the applicability of these theories. However, discussions and research reported in literature on this topic are insufficient or limited. In this paper we introduce five variables in scattering theory to help us discuss conditions under which the Rytov approximation, in the form of the finite frequency sensitivity kernels (RFFSK), the basis of waveform inversion and tomography, is valid. The five variables are propagation length L, heterogeneity scale a, wavenumber k, anisotropy ratio ξ, and perturbation strength ɛ. Combined with theoretical analysis and numerical experiments, we conclude that varying the conditions used to establish the Rytov approximation can lead to uninterpretable or undesired results. This conclusion has two consequences. First, one cannot rigorously apply the linear Rytov approximation to all theoretical or practical cases without discussing its validity. Second, the nonlinear Rytov approximation is essential if the linear Rytov approximation is not valid. Different from previous literature, only phase (or travel time) terms for the whole wavefield are discussed. The time shifts of two specific events between the background and observed wavefields measured by cross-correlation will serve as a reference for evaluation of whether the time shifts predicted by the FFSKs are reasonably acceptable. Significantly, the reference “cross-correlation” should be regarded as reliable only if the condition “two specific similar signals” is satisfied. We cannot expect it to provide a reasonable result if this condition is not met. This paper reports its reliability and experimental limitations. Using cross-correlation (CC) samples as the X axis and sensitivity kernel (SK) or ray tracing (RT) samples as the Y axis, a chart of cross validation among the three is established. According to experimental analysis for a random and deterministic medium, the RFFSK works well for diffraction and geometric optics enclosed by the boundaries ka > 1 and (Λ 1, Φ

Due to incomplete aperture coverage and complex overburden structures, the migration process cannot provide a true-amplitude image even though a true-amplitude propagator is used. Amplitude compensation based on source-side illumination ignores the aperture effects on the receiver side, and it may fail to recover the true-reflection/scattering strength of a geologic structure from the image. The structural dip largely controls if the wave incident on the structure can be reflected back and received by the acquisition aperture, so it should be taken into account in removing the acquisition effects from the migration image. We derived a dip-angle domain amplitude correctionfrom the resolution theory. The stacked migration image created by reverse time migration was decomposed into common dip images, which were compensated individually by the corresponding amplitude correction factor. Then, we summed up the corrected images to form a final image. To construct the amplitude correction factor, we generated a monofrequency Green’s function at the shot/geophone location and further decomposed it into incident/scattered plane waves. They were combined based on Snell’s law to construct correction factors for different dips. The final amplitude correction factor was formed by visiting all the shot-geophone pairs in the observation system. We devised efficient algorithms to make the amplitude correction more practical. We evaluated two numerical experiments, a five-layer model and the SEG/EAGE salt model, in which the amplitude correction led to a scalar/pressure image with an amplitude better matching the true impedance contrasts of subsurface structures, especially in areas with steep dips and in subsalt regions.