P-wave receiver functions (RFs), which utilize P-wave conversions to probe subsurface structures, face significant challenges in sedimentary environments. Specifically, strong reverberations generated by ultra-low-velocity sedimentary layers distort RF waveforms and obscure crustal signals, posing challenges for robust shallow crustal imaging. We develop a novel Bayesian joint inversion framework that simultaneously utilizes three complementary data sets-reverberant receiver functions, dereverberated receiver functions, and surface wave dispersion-to address this challenge. Our approach employs Unscented Kalman Inversion, a derivative-free method that efficiently handles nonlinear joint inversion problems. Synthetic tests demonstrate that our joint inversion recovers sediment thickness and Moho depth with uncertainties of +/- 0.50 km and +/- 1.0 km, respectively. Application to real data from the Songliao Basin verifies the approach, successfully reconstructing sediment thickness and Moho depth beneath sedimentary cover. This methodology demonstrates potential for advancing crustal investigations in complex sedimentary settings, such as continental rift basins and oceanic margins, where sedimentary sequences of variable thickness often obscure deeper structures.

Abstract The Jiaodong gold deposits are the largest gold mining base in China, with a rare mineralization scale globally. Their formation differs from existing metallogenic models and may represent a unique type of gold deposit that has not yet been fully and systematically defined. Ore-controlling structures are key factors controlling fluid migration and precipitation, and their study is crucial for understanding metallogenic law and guiding deep exploration. Based on data from a linear short-period dense array spanning the major tectonic units and ore clusters in Jiaodong, this study employs double beamforming tomography to obtain a high-resolution S-wave velocity structure within the upper 8 km of the crust in the region, systematically revealing the major ore-controlling structures and the shallow architecture of the boundary between the North China craton (NCC) and the Yangtze craton (YC). The results show that the Jiaoxibei ore cluster is primarily controlled by a series of low-angle detachment faults, such as the Sanshandao, Jiaojia, Zhaoping, and Qixia faults, which exhibit distinct listric low-velocity anomalies. In contrast, the Muping ore cluster is controlled by steeply dipping faults, such as the Muru fault, which show near-vertical low-velocity anomalies. Significant differences exist in the crustal structure across the Wulian–Yantai fault zone. The eastern YC shows signs of subduction beneath the western NCC, suggesting that this fault zone may be the northern boundary of the collisional orogenic belt between the NCC and YC. Double beamforming tomography effectively enhances weak, but coherent signals and clearly delineates the key structures within the major ore clusters, providing important seismological evidence for understanding the metallogenic process and mechanism of the Jiaodong gold deposits.

We develop a Fortran package with high programming optimization and parallel computing for simulating high-frequency (>1 Hz) teleseismic wavefields using a hybrid numerical method that couples the finite-difference (FD) and frequency-wavenumber (FK) methods. This method can simulate the interactions of incoming teleseismic wavefields with local heterogeneities but reduce computational region to a much smaller localized domain, which can significantly reduce the computing cost of the high-frequency teleseismic wavefields. The local heterogeneities are allowed to vary arbitrarily in a localized heterogeneous domain. In this package, the geographical locations of earthquakes are permitted, which can consider the real azimuthal effect of the source. Numerical benchmark tests first demonstrate the effectiveness of the developed method for P- and S-wave receiver functions (RFs). The consistent travel times of synthetic and theoretical RFs phases demonstrate its high accuracy. Application on a dense array generally obtains consistent RFs profiles with observed ones and successfully reproduces the observed common-converted-point (CCP) stacking image, which further verifies the effectiveness of the presented method. In addition, statistics of the timeconsuming of typical models illustrate the high efficiency of this package, which needs very little computing resources even to be feasible on a laptop.

Delineating subsurface interfaces is a crucial step in site selection and characterization for various subsurface applications, such as the geologic carbon sequestration, the radioactive waste disposal and hydrocarbon exploration and production. 3D seismic surveys are widely used for identifying subsurface interfaces and geologic features. Due to the large volumes and the complexity of seismic data, manual interpretation of subsurface interfaces is extremely time-consuming, and the interpretation results can be greatly affected by the subjectivity of a particular interpreter. With the latest advances in deep neural networks (DNNs), automatic seismic interpretation methods based on DNNs emerged. Most of the DNN-based seismic interpretation methods are supervised learning methods, which require large amount of labeled data for network training. We have developed an unsupervised learning method with deep fully convolutional networks (FCNs) for rapid subsurface interface identification based on self-learning algorithms, which does not require manual data labeling and specific training datasets. The characteristics of subsurface interfaces are represented as numerical constraints added to the specially designed loss function for constructing the FCN model. Physical constraints are further applied in postprocessing the outputs of the FCN model to improve the resolution of the identified subsurface interfaces. Application of the unsupervised learning method on a real seismic dataset collected at a potential CO2 storage site demonstrates that the proposed method yields rapid and accurate identification of subsurface interfaces with relatively strong acoustic impedance contrast in seismic images. The proposed approach may assist in automatic subsurface interface identification in real time and facilitate building geological models for subsurface applications.

Joint inversion, such as the combination of receiver function and surface wave dispersion, can significantly improve subsurface imaging by exploiting their complementary sensitivities. Bayesian methods have been demonstrated to be effective in this field. However, there are practical challenges associated with this approach. Notably, most Bayesian methods, such as the Markov Chain Monte Carlo method, are computationally intensive. Additionally, accurately determining the data noise across different data sets to ensure effective inversion is often a complex task. This study explores the unscented Kalman inversion (UKI) as a potential alternative. Through a data-driven approach to adjust estimated noise levels, we can achieve a balance between actual noise and the weights assigned to different data sets, enhancing the effectiveness of the inversion process. Synthetic tests of joint inversion of receiver function and surface wave dispersions indicate that the UKI can provide robust solutions across a range of data noise levels. Furthermore, we apply the UKI to real data from seismic arrays in Pamir and evaluate the accuracy of the joint inversion through posterior Gaussian distribution. Our results demonstrate that the UKI presents a promising supplement to conventional Bayesian methods in the joint inversion of geophysical data sets with superior computational efficiency.

SUMMARY Higher frequency receiver function (RF) analysis based on dense nodal arrays has been widely used for imaging crustal structures. However, the scattered Rayleigh waves generated by the steep topography including mountain ranges and basin-range junction zones, have become a significant interference that can lead to false structures in RF images. In this study, we propose a novel method to remove scattered Rayleigh waves from RF profiles by using a high-resolution linear Radon transform. Based on the difference in the apparent velocity of Rayleigh and converted waves at interfaces, we construct a scheme to design an optimal filter mask. Synthetic and observed data show that this method can be an effective tool to remove high-amplitude Rayleigh waves and preserve low-amplitude converted waves almost harmlessly. Modelling tests also show that it is suitable for non-uniform station spacing, white noise and models that include dipping interfaces.

First-arrival traveltime tomography has been widely used for upper crustal velocity modeling, but it usually suffers from the problem of complex surface topography. To overcome this problem, we have developed a new topography-dependent eikonal tomography scheme that combines a developed accurate and efficient traveltime modeling method and introduces a flexible and robust adjoint inversion scheme in the presence of irregular topography. A surface-flattening scheme is used to handle the irregular surface, where the real model is discretized by curvilinear grids and the irregular free surface is mathematically flattened through the transformation from Cartesian to curvilinear coordinates. Based on this parameterization, the forward traveltime modeling is conducted by a monotone fast-sweeping method that discretizes the factored topography-dependent eikonal equation with a point-source condition. This algorithm can circumvent the source-singularity problem and decrease the numerical error in the vicinity of a point source in the curvilinear system. Then, the gradient-based inversion is used to minimize the misfit function, which is achieved by a matrix-free adjoint-state method without cumbersome ray tracing and explicit estimation of the Fréchet derivative matrix in the curvilinear coordinate system. The new tomographic scheme is evaluated through numerical examples with different seismic structures with complex topography, and then applied to a wide-angle profile acquired in the northeastern Tibetan Plateau. The results validate the effectiveness and efficiency of our tomography scheme in constructing shallow crustal velocity models with irregular topography.

In the geophysical joint inversion, the gradient and Bayesian Markov Chain Monte Carlo (MCMC) sampling-based methods are widely used owing to their fast convergences or global optimality. However, these methods either require the computation of gradients and easily fall into local optimal solutions, or cost much time to carry out the millions of forward calculations in a huge sampling space. Different from these two methods, taking advantage of the recently developed unscented Kalman method in computational mathematics, we extend an iterative gradient-free Bayesian joint inversion framework, i.e., Multi-task unscented Kalman inversion (MUKI). In this new framework, information from various observations is incorporated, the model is iteratively updated in a derivative-free way, and a Gaussian approximation to the posterior distribution of the model parameters is obtained. We apply the MUKI to the joint inversion of receiver functions and surface wave dispersion, which is well-established and widely used to construct the crustal and upper mantle structure of the earth. Based on synthesized and real data, the tests demonstrate that MUKI can recover the model more efficiently than the gradient-based method and the Markov Chain Monte Carlo method, and it would be a promising approach to resolve the geophysical joint inversion problems.

Based on the recently developed theory of Unscented Kalman Inversion in computational mathematics, we proposed a Bayesian joint inversion framework, i.e., Multi-task Unscented Kalman Inversion (MTUKI), and apply it to the joint inversion of receiver function (RF) and surface wave dispersion (SWD). This method can share information between different observations in a derivative-free way and provide an efficient Gaussian approximation to the posterior distribution of model parameters (thickness and S-wave velocity in each layer of media). The theory and experiments show that our proposed framework demonstrates superior performance in terms of robustness, accuracy, and high efficiency.

The conventional elastic least-squares reverse time migration (LSRTM) generally inverts the parameter perturbation of the model rather than the reflectivity of reflected P- and S-modes, which leads to difficulty in directly interpreting the physical properties of the subsurface media. However, an accurate velocity model that is needed by the separation of seismic records of conventional LSRTM is usually unavailable in real data, which limits its application. In this study, we introduce a new practical correlative LSRTM (CLSRTM) scheme based on wave mode decomposition without amplitude and phase distortion, which frees from separation of seismic records. In this study, we deduced the migration and the de-migration operators using the decoupled P- and S-wave equations in heterogeneous media, which needs no extra wavefield decomposition in simulated data. To accelerate the convergence and improve the efficiency of the inversion, we adopted an analytical step-length formula that can be incidentally computed during the necessary de-migration process and the L-BFGS algorithm. Two numerical examples demonstrate that the proposed method can compensate the energy of deep structures, and generate clear images with balanced amplitudes and enhanced resolution even for the fault structures beneath the salt dome.

Seismic wavefield numerical simulation is an important base for crust-mantle structure imaging and deep exploration. The classical teleseismic wavefield simulation methods, including analytical method, semi-analytical method and numerical method, are mainly based on one-dimensional Earth model. These algorithms can efficiently calculate the synthetic seismogram, but the lateral heterogeneity of the medium is not being incorporated. With the improvement of computer performance, numerical simulation methods for three-dimensional seismic have been developed rapidly, and have been widely used in both local and regional seismic wave simulations. But due to the limitation of computing resources, the implementation of high frequency seismic wavefield numerical simulation based on global scale is still a great challenge. In recent years, hybrid numerical simulation methods of teleseismic wavefield have been developed, in which the target simulation region is decomposed into two scales (global scale and local scale). In global scale, high-frequency synthetic seismogram is calculated through fast algorithm based on one-dimensional earth model assumption. With injection method, three-dimensional numerical methods (spectral element method, finite difference method, etc.) are used to simulate the propagation of seismic wave in three-dimensional heterogeneous medium in local target scale, to achieve the balance between efficiency and accuracy. With the development of dense array observation, scientific research puts forward higher requirements for the resolution of underground structure imaging. Accurate and efficient hybrid simulation method of seismic wavefield will play an important role in the field of high-resolution seismic imaging. In this paper, we systematically summarize one-dimensional simulation methods of teleseismic wavefield numerical simulations, as well as the principle and application of hybrid numerical simulation methods of teleseismic wavefield.

High-efficient seismic wave forward modeling is important for the investigation of seismic wave phenomenon in complex media and the imaging of subsurface structures of the Earth. In the framework of spectral-element method (SEM), we improve the computation efficiency of the previous element-by-element method and propose a new element-by-element parallel spectral-element method (EBE-SEM) for solving 3D seismic wave equation to obtain seismic wavefield. The essential idea of EBE-SEM is that the products of stiffness matrix and solution vector are operated on element level, and the products are equally distributed to each CPU processor. This strategy can greatly increase the parallelization of SEM. Because Gauss-Legendre-Lobatto integration points coincide with the interpolation points in SEM, the storage of dense element stiffness matrix can be transformed to the storage of the determinant and inversion of element Jacobian matrix, and therefore the efficiency of EBE-SEM can be further increased. To eliminate the reflected wave from truncated boundary, we solve the second-order perfectly matched layer (PML) absorbing boundary condition that is constructed at the boundary of the computational domain. The validity and efficiency of the EBE-SEM are demonstrated by two numerical examples.

The separation of P- and S-wavefields is considered to be an effective approach for eliminating wave-mode crosstalk in elastic reverse time migration (RTM). At present, the Helmholtz decomposition method is widely used for isotropic media. However, it tends to change the amplitudes and phases of the separated wavefields compared with the original wavefields. Other methods used to obtain pure P- and S-wavefields include the application of the elastic wave equations of the decoupled wavefields. To achieve high computational accuracy, staggered-grid finite-difference (FD) schemes are usually used to numerically solve the equations by introducing an additional stress variable. However, the computational cost of this method is high because a conventional hybrid wavefield (the P- and S-wavefields are mixed together) simulation must be created before the P- and S-wavefields can be calculated. We have developed the first-order particle velocity equations to reduce the computational cost. The equations can describe four types of particle-velocity wavefields: the vector P-wavefield, the scalar P-wavefield, the vector S-wavefield, and the vector S-wavefield rotated in the direction of the curl factor. Without introducing the stress variable, only the four types of particle velocity variables are used to construct the staggered-grid FD schemes; therefore, the computational cost is reduced. We also develop an algorithm to calculate the P and S propagation vectors using the four particle velocities, which is simpler than the Poynting vector. Finally, we apply the velocity equations and propagation vectors to elastic RTM and angle-domain common-image gather computations. These numerical examples illustrate the efficiency of our methods.

A rock mass often contains joints filled with a viscoelastic medium of which seismic response is significant to geophysical exploration and seismic engineering design. Using the propagator matrix method, an analytical model was established to characterize the seismic response of viscoelastic filled joints. Stress wave propagation through a single joint highly depended on the water content and thickness of the filling as well as the frequency and incident angle of the incident wave. The increase in the water content enhanced the viscosity (depicted by quality factor) of the filled joint, which could promote equivalent joint stiffness and energy dissipation with double effects on stress wave propagation. There existed multiple reflections when the stress wave propagated through a set of filled joints. The dimensionless joint spacing was the main controlling factor in the seismic response of the multiple filled joints. As it increased, the transmission coefficient first increased, then it decreased instead, and at last it basically kept invariant. The effect of multiple reflections was weakened by increasing the water content, which further influenced the variation of the transmission coefficient. The water content of the joint filling should be paid more attention in practical applications.

The 2014 Ludian M(s)6.5 earthquake caused dramatic death and significant economic losses, and received much attention. Many studies have been made on the focal mechanisms of this earthquake sequence. However, due to the limitation of the inversion method and the complex structures of the seismic area, there are still some disputes on this issue, especially concerned with the dip-angle causative fault of the mainshock. In this work, we use regional seismic waveform inversion based on the gCAP method and the full waveform modeling to invert focal mechanism solutions of the mainshock and 9 aftershocks of the Ludian earthquake sequence. Results show that the azimuthal coverage has a significant effect on the reliability and stability of gCAP's inversion. To solve this problem, a two-bandwidth method is proposed in this paper, in which the bandwidths are 0. 01 similar to 0. 05 Hz and 0. 01 similar to 0. 2 Hz, respectively. With this inversion, this work determines that the dip angle and rake of the mainshock 76 degrees similar to 83 degrees and 157 degrees similar to 164 degrees, respectively, implying a high-dip strike-slip rupture. The stabilities of the gCAP method and the full waveform modeling are verified by the good consistence between the derived mechanisms of the aftershocks. Our research suggests that the high-dip character of the Ludian earthquake might lead to the fast and complete energy release, which may be the main reason for the aggravating destruction and the lack of larger aftershocks.

The North China Craton (NCC) is a key region to study the destruction of the ancient craton. Two groups of phases (denoted as “Pw1” and “Pw2”), which are parallel to the PmP phase reflected from the Moho discontinuity and the PLP phase reflected from the Lithosphere and Asthenosphere Boundary (LAB) respectively, are found on the record section of the Rongcheng-Xinzhou-Alxa long-range deep seismic sounding profile. The nature of the two phases is still unclear, although they are clearly observable and reverberant. In this paper, we use travel time inversion and amplitude forward modelling to fit the reflected and refracted phases in the lithosphere. The results show: (1) the Pw1 is a multiple reflected phase which is successively reflected by the crystalline basement, the surface, the Moho and then finally received on the surface; (2) the Pw2 phase is also a multiple reflected phase successively reflected by the crystalline basement, the surface, the LAB interface and then received on the surface. We conclude that the significant velocity difference between the thick sedimentary cover and the crystalline basement in the North China rifted basin may be the main reason for generating the multiple reflections. Furthermore, the two multiple reflections provide potent constraints on the lithospheric velocity model, and constitute seismological evidence for the lithospheric thinning in the eastern NCC.

近十年来,短周期密集台阵被动源探测技术日益成长为国内外深部结构探测领域的一项重要手段.该技术相较于传统宽频带地震探测具有高分辨、省时省钱、绿色环保等优点.尽管低频信号不足,对岩石圈地幔以深结构探测能力有限,但密集的台站间距使得地壳精细结构成像成为可能;台阵下方不同方位射线形成密集的交叉覆盖,从而可通过反演和叠加偏移等手段获取稳定的壳幔结构图像.因此,短周期密集台阵探测技术已广泛应用于深部速度和界面结构成像,以及矿产资源勘查、火山活动监测、微震精定位、发震断层的几何学和运动学特征研究等多种不同领域.本文系统的总结了短周期密集台阵在地壳结构研究和微震定位检测等方面的研究进展.展望未来,以科学问题为导向,利用天然地震及背景噪声观测、重力测量、InSAR数据、GPS测量等多种地球物理数据,开展联合反演和成像,并提取研究对象的多属性特征正日益成为减少解的非唯一性、揭示地质体真实赋存状态的有效途径;该领域方法技术的迅速提升有望大大促进地震学及地球动力学研究,在深部结构成像、矿产资源勘查等领域具有广阔应用前景.

SUMMARY In this paper, we propose a spectral element method (SEM) based on unstructured tetrahedral grids for direct current (dc) resistivity modelling. Unlike the tensor-product of 1-D Gauss–Lobatto–Legendre (GLL) quadrature in conventional SEM, we use Proriol–Koornwinder–Dubiner (PKD) polynomials to form the high-order basis polynomials on tetrahedral grids. The final basis functions are established by using Vandermonde matrix. Compared to traditional SEM, our method truly takes into account the high precision of spectral method and the flexibility of finite element method with unstructured grids for modelling the complex underground structures. After addressing the theory on the construction of basis functions and interpolation and integration nodes, we validate our algorithm using the analytical solutions for a layered earth model and the results from other methods for multiple geoelectrical models. We further investigate a dual-track scheme for improving the accuracy of our SEM by increasing the order of interpolation polynomials or by refining the grids.

青藏高原是全球海拔最高、规模最大、时代最新的陆-陆碰撞造山带.几十个百万年以来高原隆升、喜马拉雅山系崛起是地球演化史上最为壮观的构造事件之一.青藏高原壳幔结构和深部过程备受国内外地球科学家关注.近七十年来的地球物理研究与探索表明:(1)青藏高原地壳巨厚,岩石圈较薄;(2)壳内存在软弱层,但厚度和联通性有限;(3)高原下地壳及Moho面广泛发育叠瓦状反射特征,存在明显的脆性变形;(4)喜马拉雅和拉萨块体南部存在双Moho现象/迹象;(5)印度大陆岩石圈向高原下方俯冲的形态存在显著的东西向差异;(6)高原主体上地幔各向异性以NEE向为主;(7)青藏高原布格重力异常四周高、中间低,异常场边界与地形梯度变化密切相关;(8)高原内部磁异常较弱,周边地区较强,其分界与区域构造边界基本一致;(9)青藏高原水热活动强烈,大地热流值高,主要来自加厚地壳的贡献.但是,有关青藏高原深部过程,诸如是否存在中/下地壳流、印度与欧亚大陆岩石圈的俯冲模式等重大科学问题目前仍存在争议.青藏高原地球物理和动力学研究是一个复杂的系统工程,以科学问题为导向,结合国家重大需求,在关键区域组织实施综合地球物理探测,可望在地学领域取得创新与突破.

The characteristics and formation of the Yadong-Dongqiao-Huluhu Belt (YDHB) are important to the understanding of crustal deformation, material movement and deep dynamics of the Tibetan Plateau. Analyses on multiple elements suggests that: (1) The YDHB is an intracontinental rift under EW-directed extensional force as evidenced by spatial distribution of crust and mantle structures; (2) High seismicity, distribution of abnormal terrestrial heat flow and the stress field by mantle convection indicate that it is an active continental rift; (3) The formation and evolution of the YDHB is the product of intense exchange of material and energy in the Earth's interiors.