Lp-NORM DECONVOLUTION1

The purpose of deconvolution is to retrieve the reflectivity from seismic data. To do this requires an estimate of the seismic wavelet, which in some techniques is estimated simultaneously with the reflectivity, and in others is assumed known. The most popular deconvolution technique is inverse filtering. It has the property that the deconvolved reflectivity is band-limited. Band-limitation implies that reflectors are not sharply resolved, which can lead to serious interpretation problems in detailed delineation. To overcome the adverse effects of band-limitation, various alternatives for inverse filtering have been proposed. One class of alternatives is L_p-norm deconvolution, L₁norm deconvolution being the best-known of this class. We show that for an exact convolutional forward model and statistically independent reflectivity and additive noise, the maximum likelihood estimate of the reflectivity can be obtained by L_p-norm deconvolution for a range of multivariate probability density functions of the reflectivity and the noise. The L_∞-norm corresponds to a uniform distribution, the L₂-norm to a Gaussian distribution, the L₁-norm to an exponential distribution and the L₀-norm to a variable that is sparsely distributed. For instance, if we assume sparse and spiky reflectivity and Gaussian noise with zero mean, the L_p-norm deconvolution problem is solved best by minimizing the L₀-norm of the reflectivity and the L₂-norm of the noise. However, the L₀-norm is difficult to implement in an algorithm. From a practical point of view, the frequency-domain mixed-norm method that minimizes the L₁norm of the reflectivity and the L₂-norm of the noise is the best alternative. L_p-norm deconvolution can be stated in both time and frequency-domain. We show that both approaches are only equivalent for the case when the noise is minimized with the L₂-norm. Finally, some L_p-norm deconvolution methods are compared on synthetic and field data. For the practical examples, the wide range of possible L_p-norm deconvolution methods is narrowed down to three methods with p= 1 and/or 2. Given the assumptions of sparsely distributed reflectivity and Gaussian noise, we conclude that the mixed L₁norm (reflectivity) L₂-norm (noise) performs best. However, the problems inherent to single-trace deconvolution techniques, for example the problem of generating spurious events, remain. For practical application, a greater problem is that only the main, well-separated events are properly resolved.     

Determination of source-time function of intermediate and small earthquakes from empirical Green’s functions

The regularization method was used to invert the source-time function of four larger aftershocks of the Luquan, Yun-nan,M _s = 6.1 earthquake of April 18, 1984. Near-field digital accelerograms were used in the inversion. Instead of calculating synthetic seismograms theoretically, an accelerogram of smaller earthquake was used as an empirical Green’s function. The results obtained showed that the source-time function of smaller aftershock was usually a simple pulse, while that of larger events was more or less complicated. The rupture velocities and the average particle velocities of the aftershocks were determined. All the results were in good agreement with the result estimated from the present seismic source theory. The Chinese version of this paper appeared in the Chinese edition ofActa Seismologica Sinica,15, 22–31, 1993. This study is supported by the Western Yunnan Experimental Site for Earthquake Prediction and the Chinese Joint Seismological Science Foundation.     

Source Parameters for 1999 North Anatolian Fault Zone Aftershocks

We develop a data set of aftershock recordings of the 1999, M = 7.4 Izmit and M = 7.2 Duzce (Turkey) earthquakes to study their source parameters. We combined seismograms from 44 stations maintained by several sources (organizations) to obtain a unified data set of events (2.1 ≤ M_w ≤ 5.5). We calculate source parameters of these small earthquakes by two methods that use different techniques to address the difficulty in obtaining source spectra for small earthquakes subject to interference from site response. One method (program NetMoment (NM), Hutchings, 2004) uses spectra of direct S waves in a simultaneous inversion of local high-frequency network data to estimate seismic moment, source corner frequency (f_c), site attenuation (k) and whole-path Q. This approach takes advantage of the source commonality in all recordings for a particular earthquake by fitting a common Brune source spectrum to the data with a and individual k. The second approach (Mayeda et al., 2003) uses the coda method (CM) to obtain “nonmodel-based” source spectra and moment estimates from selected broadband recording sites. We found that both methods do well for events that allow the comparison with seismic moment estimates derived from waveform modeling. Also, source spectra obtained from the two methods are very closely matched for most of the events they have in common. We use an F test to examine the trade-off between k and f_c picks identified by the direct S-wave method. About half of the events could be constrained to have less than a 50% average uncertainty in f_c and k. We used these source spectra solutions to calculate energy and apparent stress and compare these to estimates from the selected “good quality” source spectra from CM. Both studies have values mutually consistent and show a similar increase in apparent stress with increasing moment. This result has added merit due to the independent approaches to calculate apparent stress. We conclude that both methods are at least partially validated by our study, and they both have usefulness for different circumstances of recording local small earthquakes. CM would work well in studies for which there is a broad magnitude range of events and NM works well for local events recorded by band-limited recorders.     

Inversion of long-period body-wave data for the source process of the Gonghe, Qinghai, China earthquake

Long period body waves data recorded by the China Digital Seismograph Network (CDSN) are inverted for the seismic moment tensors of the April 26, 1990, Gonghe, QinghaiM _S=6.9 earthquake and itsM _S=5.0 after-shock occurred on May 7, 1990. In the inversion, the generalized reflection-transmission coefficient matrix method is used to generate Green’s function. From the inversion it is obtained that the rupture process of theM _S=5.0 aftershock is relatively simple, and that of the main shock is rather complex. There are at least two events during main shock rupture process with an interval about 35 seconds. The focal mechanisms of two events are roughly the same as that of the aftershock, all of them were mainly reverse dip-slipping faulting with minor left-lateral strike-slip motion. These results indicate that the Gonghe earthquake was the result of the farther extension of one NWW-SEE striking buried fault on the southern margin of Gonghe basin from shallower depth to deeper depth and from NW to SE under the action of a nearly horizontal NE direction compressive stress. Contribution No. 95A0111, Institute of Geophysics, SSB, China.     

Iterated extended Kalman filter method for time-lapse seismic full-waveform inversion

Time-lapse seismic data is useful for identifying fluid movement and pressure and saturation changes in a petroleum reservoir and for monitoring of CO₂ injection. The focus of this paper is estimation of time-lapse changes with uncertainty quantification using full-waveform inversion. The purpose of also estimating the uncertainty in the inverted parameters is to be able to use the inverted seismic data quantitatively for updating reservoir models with ensemble-based methods. We perform Bayesian inversion of seismic waveform data in the frequency domain by combining an iterated extended Kalman filter with an explicit representation of the sensitivity matrix in terms of Green functions (acoustic approximation). Using this method, we test different strategies for inversion of the time-lapse seismic data with uncertainty. We compare the results from a sequential strategy (making a prior from the monitor survey using the inverted baseline survey) with a double difference strategy (inverting the difference between the monitor and baseline data). We apply the methods to a subset of the Marmousi2 P-velocity model. Both strategies performed well and relatively good estimates of the monitor velocities and the time-lapse differences were obtained. For the estimated time-lapse differences, the double difference strategy gave the lowest errors.     

CO2 saturation estimates at Sleipner (North Sea) from seismic tomography and rock physics inversion

CO₂ saturations are estimated at Sleipner using a two-step imaging workflow. The workflow combines seismic tomography (full-waveform inversion) and rock physics inversion and is applied to a two-dimensional seismic line located near the injection point at Sleipner. We use baseline data (1994 vintage, before CO₂ injection) and monitor data that was acquired after 12 years of CO₂ injection (2008 vintage). P-wave velocity models are generated using the Full waveform inversion technology and then, we invert selected rock physics parameters using an rock physics inversion methodology. Full waveform inversion provides high-resolution P-wave velocity models both for baseline and monitor data. The physical relations between rock physics properties and acoustic wave velocities in the Utsira unconsolidated sandstone (reservoir formation) are defined using a dynamic rock physics model based on well-known Biot–Gassmann theories. For data prior to injection, rock frame properties (porosity, bulk and shear dry moduli) are estimated using rock physics inversion that allows deriving physically consistent properties with related uncertainty. We show that the uncertainty related to limited input data (only P-wave velocity) is not an issue because the mean values of parameters are correct. These rock frame properties are then used as a priori constraint in the monitor case. For monitor data, the Full waveform inversion results show nicely resolved thin layers of CO₂–brine saturated sandstones under intra-reservoir shale layers. The CO₂ saturation estimation is carried out by plugging an effective fluid phase in the rock physics model. Calculating the effective fluid bulk modulus of the brine–CO₂ mixture (using Brie equation in our study) is shown to be the key factor to link P-wave velocity to CO₂ saturation. The inversion tests are done with several values of Brie/patchiness exponent and show that the CO₂ saturation estimates are varying between 0.30 and 0.90 depending on the rock physics model and the location in the reservoir. The uncertainty in CO₂ saturation estimation is usually lower than 0.20. When the patchiness exponent is considered as unknown, the inversion is less constrained and we end up with values of exponent varying between 5 and 20 and up to 33 in specific reservoir areas. These estimations tend to show that the CO₂–brine mixing is between uniform and patchy mixing and variable throughout the reservoir.     

Petrophysical inversion of AVA data

We investigate the interactions between the elastic parameters, V_P, V_S and density, estimated by non-linear inversion of AVA data, and the petrophysical parameters, depth (pressure), porosity, clay content and fluid saturation, of an actual gas-bearing reservoir. In particular, we study how the ambiguous solutions derived from the non-uniqueness of the seismic inversion affect the estimates of relevant rock properties. It results that the physically admissible values of the rock properties greatly reduce the range of possible seismic solutions and this range contains the actual values given by the well. By means of a statistical inversion, we analyse how approximate a priori knowledge of the petrophysical properties and of their relationships with the seismic parameters can be of help in reducing the ambiguity of the inversion solutions and eventually in estimating the petrophysical properties of the specific target reservoir. This statistical inversion allows the determination of the most likely values of the sought rock properties along with their uncertainty ranges. The results show that the porosity is the best-resolved rock property, with its most likely value closely approaching the actual value found by the well, even when we insert somewhat erroneous a priori information. The hydrocarbon saturation is the second best-resolved parameter, but its most likely value does not match the well data. The depth of the target interface is the least-resolved parameter and its most likely value is strongly dependent on a priori information. Although no general conclusions can be drawn from the results of this exercise, we envisage that the proposed AVA–petrophysical inversion and its possible extensions may be of use in reservoir characterization.     

Simultaneous inversion of velocity structures and hypocentral locations: Application to the Friuli seismic area NE Italy

A layeredP- andS-wave velocity model is obtained for the Friuli seismic area using the arrival time data ofP- andS-waves from local earthquakes. A damped least-squares method is applied in the inversion.The data used are 994P-wave arrival times for 177 events which have epicenters in the region covered by the Friuli seismic network operated by Osservatorio Geofisico sperimentale (OGS) di Trieste, which are jointly inverted for the earthquake hypocenters andP-wave velocity model. TheS-wave velocity model is estimated on the basis of 978S-wave arrival times and the hypocenters obtained from theP-wave arrival time inversion. We also applied an approach thatP- andS-wave arrival time data are jointly used in the inversion (Roecker, 1982). The results show thatS-wave velocity structures obtained from the two methods are quite consistent, butP-wave velocity structures have obvious differences. This is apparent becauseP-waves are more sensitive to the hypocentral location thanS-waves, and the reading errors ofS-wave arrival times, which are much larger than those ofP-waves, bring large location errors in the joint inversion ofP- andS-wave arrival time. The synthetic data tests indicated that when the reading errors ofS-wave arrivals are larger than four times that ofP-wave arrivals, the method proposed in this paper seems more valid thanP- andS-wave data joint inversion. Most of the relocated events occurred in the depth range between 7 and 11 km, just above the biggest jump in velocity. This jump might be related to the detachment line hypothesized byCarulli et al. (1982). From the invertedP- andS-wave velocities, we obtain an average value 1.82 forV _p /V _s in the first 16 km depth.     

2014年金寨M_L3.9震群序列遗漏地震检测及发震构造分析

针对2014年8月—2015年1月安徽金寨发生的M_L3.9震群,利用匹配滤波技术补充台网目录遗漏的地震事件,再利用波形互相关震相检测技术标定P波和S波到时,进一步采用双差定位方法对震群进行重定位,结合震源机制解等分析此次震群活动可能的发震构造。计算结果显示,通过互相关扫描检测到1376个地震台网常规分析遗漏的地震,数量约为台网目录给出的585个事件的2.35倍。检测到的遗漏地震震级估算为M_L0～2.3,通过震级-频次统计分析,加入遗漏地震后地震目录的完整性在M_L0～1.5范围内有较明显的改善。重定位后地震走时残差更小,水平位置更集中,沿NNE向断裂F和NW向青山-晓天断裂呈现近直立的条带状分布。结合地质构造、震源机制解和水库因素,推测2014年金寨M_L3.9震群可能是由周边水库水下渗引起NW向青山-晓天断裂与NNE向断裂F慢滑动而触发的。     

地震偏移反演成像的迭代正则化方法研究

利用伴随算子L^*,直接的偏移方法通常导致一个低分辨率或模糊的地震成像.线性化偏移反演方法需求解一个最小二乘问题.但直接的最小二乘方法的数值不稳定,为目视解译带来困难.本文建立约束正则化数学模型,研究了地震偏移反演成像问题的迭代正则化求解方法.首先对最小二乘问题施加正则化约束,接着利用梯度迭代法求解反演成像问题,特别是提出了共轭梯度方法的混合实现技巧.为了表征该方法的可实际利用性,分别对一维,二维和三维地震模型进行了数值模拟.结果表明该正则偏移反演成像方法是有效的,对于实际的地震成像问题有着良好的应用前景.     

Precise determination of focal parameters for July 20, 1995 M_L=4.1 earthquake sequence in the Huailai basin

IntroductionSince the late 1970s, the quickly developed global digital seismograph network has been providing high quality recordings of large earthquakes in global scale, based on which digital seismology has made great progress. Compared with large earthquakes, moderate and small sized shocks have more frequent occurrence, and comprise clues to geological tectonics and tectonic stress field in a region. Preceding and following a large earthquake, usually occur numbers of small events that im…     

An insightful parametrization for the flatlander's interpretation of time‐lapsed seismic data

An approximation is developed that allows mapped 4D seismic amplitudes and time‐shifts to be related directly to the weighted linear sum of pore pressure and saturation changes. The weights in this relation are identified as key groups of parameters from a petroelastic model and include the reservoir porosity. This dependence on groups of parameters explains the inherent non‐uniqueness of this problem experienced by previous researchers. The proposed relation is of use in 4D seismic data feasibility studies and inversion and interpretation of the 4D seismic response in terms of pore pressure and water saturation changes. A further result is drawn from analysis of data from the North Sea and West Africa, which reveals that the relative interplay between the effects of pore pressure and saturation changes on the seismic data can be simplified to the control of a single, spatially variant parameter C_S/C_P. Combining these results with those from published literature, we find that C_S/C_P = 8 appears to be a generality across a range of clastic reservoirs with a similar mean porosity. Using this C_S/C_P value, an in situ seismic‐scale constraint for the rock stress sensitivity component of the petroelastic model is constructed considering this component carries the largest uncertainty.     

Determination of Stress State from Focal Mechanisms of Microseismic Events Induced During Hydraulic Injection at the Hijiori Hot Dry Rock Site

—?The stress state at the Hijiori hot dry rock site was estimated based on the inversion from focal mechanisms of microseismic events induced during hydraulic injection experiments. The best fit stress model obtained by inverting 58 focal mechanisms of seismic events simultaneously indicates that the maximum principal stress σ₁ is vertical, while the minimum principal stress σ₃ is horizontal and trends north-south. The average misfit between the stress model and all the data is 6.8°. The inversion results show that the average misfit is small enough to satisfy the assumption of homogeneity in the focal mechanism data and that the 95% confidence regions of σ₁ and σ₃ are well constrained, i.e., they do not overlap, suggesting that the inversion results are acceptable. The stress estimates obtained by the focal mechanism inversion essentially agree with other stress estimates previously obtained. It is therefore concluded that the focal mechanism inversion method provides a useful tool for estimating the stress state. The hypocentral distributions of microseismic events associated with the hydraulic fracturing experiments are distributed around the plane that spreads to almost east–west from the injection wells and declines to the north at a high angle. The vertical orientation and east–west strike of the seismic events are essentially coplanar with the caldera ring-fault structure in the southern portion of the Hijiori Caldera. This indicates that tensile fractures of intact rock were not being created, but pre-existing fractures were being re-opened and developed in the direction of the maximum horizontal principal stress, although microseismic events were caused by shear failures.     

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

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

Double parameterized regularization inversion method for migration velocity analysis in transversely isotropic media with a vertical symmetry axis

Simultaneous estimation of velocity gradients and anisotropic parameters from seismic reflection data is one of the main challenges in transversely isotropic media with a vertical symmetry axis migration velocity analysis. In migration velocity analysis, we usually construct the objective function using the l₂ norm along with a linear conjugate gradient scheme to solve the inversion problem. Nevertheless, for seismic data this inversion scheme is not stable and may not converge in finite time. In order to ensure the uniform convergence of parameter inversion and improve the efficiency of migration velocity analysis, this paper develops a double parameterized regularization model and gives the corresponding algorithms. The model is based on the combination of the l₂ norm and the non‐smooth l₁ norm. For solving such an inversion problem, the quasi‐Newton method is utilized to make the iterative process stable, which can ensure the positive definiteness of the Hessian matrix. Numerical simulation indicates that this method allows fast convergence to the true model and simultaneously generates inversion results with a higher accuracy. Therefore, our proposed method is very promising for practical migration velocity analysis in anisotropic media.     

BAYESIAN ESTIMATION IN SEISMIC INVERSION. PART I: PRINCIPLES1

This paper gives a review of Bayesian parameter estimation. The Bayesian approach is fundamental and applicable to all kinds of inverse problems. Its basic formulation is probabilistic. Information from data is combined with a priori information on model parameters. The result is called the a posteriori probability density function and it is the solution to the inverse problem. In practice an estimate of the parameters is obtained by taking its maximum. Well-known estimation procedures like least-squares inversion or l₁ norm inversion result, depending on the type of noise and a priori information given. Due to the a priori information the maximum will be unique and the estimation procedures will be stable except (in theory) for the most pathological problems which are very unlikely to occur in practice. The approach of Tarantola and Valette can be derived within classical probability theory. The Bayesian approach allows a full resolution and uncertainty analysis which is discussed in Part II of the paper.     

Source time functions of the Gonghe,China earthquake retrieved from long-period digital waveform data using empirical Green’s function technique

An earthquake ofM _S=6.9 occurred at the Gonghe, Qinghai Province, China on April 26, 1990. Three larger aftershocks took place at the same region,M _S=5.5 on May 7, 1990,M _S=6.0 on Jan. 3, 1994 andM _S=5.7 on Feb. 16, 1994. The long-period recordings of the main shock from China Digital Seismograph Network (CD-SN) are deconvolved for the source time functions by the correspondent recordings of the three aftershocks as empirical Green’s functions (EGFs). No matter which aftershock is taken as EGF, the relative source time functions (RSTFs) obtained are nearly identical. The RSTFs suggest theM _S=6.9 event consists of at least two subevents with approximately equal size whose occurrence times are about 30 s apart, the first one has a duration of 12 s and a rise time of about 5 s, and the second one has a duration of 17 s and a rise time of about 8 s. Comparing the RSTFs obtained from P- and SH-phases respectively, we notice that those from SH-phases are a slightly more complex than those from P-phases, implying other finer subevents exist during the process of the main shock. It is interesting that the results from the EGF deconvolution of long-period wavform data are in good agreement with the results from the moment tensor inversion and from the EGF deconvolution of broadband waveform data. Additionally, the two larger aftershocks are deconvolved for their RSTFs. The deconvolution results show that the processes of theM _S=6.0 event on Jan. 3, 1994 and theM _S=5.7 event on Feb. 16, 1994 are quite simple, both RSTFs are single impulses. The RSTFs of theM _S=6.9 main shock obtained from different stations are noticed to be azimuthally dependent, whose shapes are a slightly different with different stations. However, the RSTFs of the two smaller aftershocks are not azimuthally dependent. The integrations of RSTFs over the processes are quite close to each other, i. e., the scalar seismic moments estimated from different stations are in good agreement. Finally the scalar seismic moments of the three aftershocks are compared. The relative scalar seismic moment of the three aftershocks deduced from the relative scalar seismic moments of theM _S=6.9 main shock are very close to those inverted directly from the EGF deconvolution. The relative scalar seismic moment of theM _S=6.9 main shock calculated using the three aftershocks as EGF are 22 (theM _S=6.0 aftershock being EGF), 26 (theM _S=5.7 aftershock being EGF) and 66 (theM _S=5.5 aftershock being EGF), respectively. Deducing from those results, the relative scalar sesimic moments of theM _S=6.0 to theM _S=5.7 events, theM _S=6.0 to theM _S=5.5 events and theM _S=5.7 to theM _S=5.5 events are 1.18, 3.00 and 2.54, respectively. The correspondent relative scalar seismic moments calculated directly from the waveform recordings are 1.15, 3.43, and 3.05.