一种改进的基于网格搜索的微地震震源定位方法

震源定位是微地震监测技术要解决的主要问题.目前,井下微地震监测多采用走时拟合法计算震源位置.常规方法受到环境噪声、初至拾取误差、速度模型误差等因素的影响,定位结果存在一定误差.为了提高定位精度,本文提出了一种改进的基于网格搜索的微地震震源定位方法.本文方法根据P波的偏振特征参数计算概率密度函数求取震源方位角,并采用改进的目标函数和搜索算法计算震源的径向距离和深度.模型数据和实际资料的处理结果表明,本文方法具有较强的抗噪性,计算得到的震源方位角更加接近真实值;与常规目标函数相比,本文方法采用的目标函数具有更好的收敛性,其定位结果受初至拾取误差和速度模型误差的影响更小;本文提出的搜索算法能够消除由于错误拾取造成的观测到时中的异常值对定位结果的影响.     

Measuring changes in fracture properties from temporal variations in anisotropic attenuation of microseismic waveforms

We investigate fracture‐induced attenuation anisotropy in a cluster of events from a microseismic dataset acquired during hydraulic fracture stimulation. The dataset contains 888 events of magnitude ?3.0 to 0.0. We use a log‐spectral‐amplitude‐ratio method to estimate change in over a half‐hour time period where fluid is being injected and an increase in fracturing from S‐wave splitting analysis has been previously inferred. A Pearson's correlation analysis is used to assess whether or not changes in attenuation with time are statistically significant. P‐waves show no systematic change in during this time. In contrast, S‐waves polarised perpendicular to the fractures show a clear and statistically significant increase with time, whereas S‐waves polarised parallel to the fractures show a weak negative trend. We also compare between the two S‐waves, finding an increase in with time. A poroelastic rock physics model of fracture‐induced attenuation anisotropy is used to interpret the results. This model suggests that the observed changes in t* are related to an increase in fracture density of up to 0.04. This is much higher than previous estimates of 0.025 ± 0.002 based on S‐wave velocity anisotropy, but there is considerably more scatter in the attenuation measurements. This could be due to the added sensitivity of attenuation measurement to non‐aligned fractures, fracture shape, and fluid properties. Nevertheless, this pilot study shows that attenuation measurements are sensitive to fracture properties such as fracture density and aspect ratio.     

A multi-channel approach for automatic microseismic event association using RANSAC-based Arrival Time Event Clustering (RATEC)

In the presence of background noise, arrival times picked from a surface microseismic data set usually include a number of false picks that can lead to uncertainty in location estimation. To eliminate false picks and improve the accuracy of location estimates, we develop an association algorithm termed RANSAC-based Arrival Time Event Clustering (RATEC) that clusters picked arrival times into event groups based on random sampling and fitting moveout curves that approximate hyperbolas. Arrival times far from the fitted hyperbolas are classified as false picks and removed from the data set prior to location estimation. Simulations of synthetic data for a 1-D linear array show that RATEC is robust under different noise conditions and generally applicable to various types of subsurface structures. By generalizing the underlying moveout model, RATEC is extended to the case of a 2-D surface monitoring array. The effectiveness of event location for the 2-D case is demonstrated using a data set collected by the 5200-element dense Long Beach array. The obtained results suggest that RATEC is effective in removing false picks and hence can be used for phase association before location estimates.     

Improvements in mining induced microseismic source locations at the Lucky Friday mine using an automated whole-waveform analysis system

For years, severe rockburst problems at the Lucky Friday mine in northern Idaho have been a persistent safety hazard and an impediment to production. An MP250 based microseismic monitoring system, which uses simple voltage threshold picking of first arrivals, has been used in this mine since 1973 to provide source locations and energy estimates of seismic events. Recently, interest has been expressed in developing a whole waveform microseismic monitoring system for the mine to provide more accurate source locations and information about source characteristics. For this study, we have developed a prototype whole-waveform microseismic monitoring system based on a 80386 computer equipped with a 50 kHz analog-digital convertor board. The software developed includes a data collection program, a data analysis program, and an event detection program. Whole-waveform data collected and analyzed using this system during a three-day test have been employed to investigate sources of error in the hypocenter location process and to develop an automatic phase picker appropriate for microseismic events.Comparison of hypocenter estimates produced by the MP250 system to those produced by the whole-waveform system shows that significant timing errors are common in the MP250 system and that these errors caused a large part of the scatter evident in the daily activity plots produced at the mine. Simulations and analysis of blast data show that analytical control over the solutions is strongly influenced by the array geometry. Within the geophone array, large errors in the velocity model or moderate timing errors may result in small changes in the solution, but outside the array, the solution is very sensitive to small changes in the data.Our whole-waveform detection program picks event onset times and determines event durations by analysis of a segmented envelope function (SEF) derived from the microseismic signal. The detection program has been tested by comparing its arrival time picks to those generated by human analysis of the data set. The program picked 87% of the channels that were picked by hand with a standard error of 0.75 milliseconds. Source locations calculated using times provided by our entire waveform detection program were similar to those calculated using hand-picked arrival times. In particular, they show far less scatter than source locations calculated using arrival times based on simple voltage threshold picking of first arrivals.     

A fast algorithm for 3D azimuthally anisotropic velocity scan

The conventional velocity scan can be computationally expensive for large‐scale seismic data sets, particularly when the presence of anisotropy requires multiparameter scanning. We introduce a fast algorithm for 3D azimuthally anisotropic velocity scan by generalizing the previously proposed 2D butterfly algorithm for hyperbolic Radon transforms. To compute semblance in a two‐parameter residual moveout domain, the numerical complexity of our algorithm is roughly as opposed to of the straightforward velocity scan, with N being the representative of the number of points in a particular dimension of either data space or parameter space. Synthetic and field data examples demonstrate the superior efficiency of the proposed algorithm.     

Parsimonious wave‐equation travel‐time inversion for refraction waves

We present a parsimonious wave‐equation travel‐time inversion technique for refraction waves. A dense virtual refraction dataset can be generated from just two reciprocal shot gathers for the sources at the endpoints of the survey line, with N geophones evenly deployed along the line. These two reciprocal shots contain approximately 2N refraction travel times, which can be spawned into refraction travel times by an interferometric transformation. Then, these virtual refraction travel times are used with a source wavelet to create N virtual refraction shot gathers, which are the input data for wave‐equation travel‐time inversion. Numerical results show that the parsimonious wave‐equation travel‐time tomogram has about the same accuracy as the tomogram computed by standard wave‐equation travel‐time inversion. The most significant benefit is that a reciprocal survey is far less time consuming than the standard refraction survey where a source is excited at each geophone location.     

Time‐lapse pre‐stack seismic data registration and inversion for CO2 sequestration study at Cranfield

Pre‐stack seismic data are indicative of subsurface elastic properties within the amplitude versus offset characteristic and can be used to detect elastic rock property changes caused by injection. We perform time‐lapse pre‐stack 3‐D seismic data analysis for monitoring sequestration at Cranfield. The time‐lapse amplitude differences of Cranfield datasets are found entangled with time‐shifts. To disentangle these two characters, we apply a local‐correlation‐based warping method to register the time‐lapse pre‐stack datasets, which can effectively separate the time‐shift from the time‐lapse seismic amplitude difference without changing the original amplitudes. We demonstrate the effectiveness of our registration method by evaluating the inverted elastic properties. These inverted time‐lapse elastic properties can be reliably used for monitoring plumes.     

Surface‐wave inversion for a P‐velocity profile with a constant depth gradient of the squared slowness

Surface waves are often used to estimate a near‐surface shear‐velocity profile. The inverse problem is solved for the locally one‐dimensional problem of a set of homogeneous horizontal elastic layers. The result is a set of shear velocities, one for each layer. To obtain a P‐wave velocity profile, the P‐guided waves should be included in the inversion scheme. As an alternative to a multi‐layered model, we consider a simple smooth acoustic constant‐density velocity model, which has a negative constant vertical depth gradient of the squared P‐wave slowness and is bounded by a free surface at the top and a homogeneous half‐space at the bottom. The exact solution involves Airy functions and provides an analytical expression for the dispersion equation. If the ratio is sufficiently small, the dispersion curves can be picked from the seismic data and inverted for the continuous P‐wave velocity profile. The potential advantages of our model are its low computational cost and the fact that the result can serve as a smooth starting model for full‐waveform inversion. For the latter, a smooth initial model is often preferred over a rough one. We test the inversion approach on synthetic elastic data computed for a single‐layer P‐wave model and on field data, both with a small ratio. We find that a single‐layer model can recover either the shallow or deeper part of the profile but not both, when compared with the result of a multi‐layer inversion that we use as a reference. An extension of our analytic model to two layers above a homogeneous half‐space, each with a constant vertical gradient of the squared P‐wave slowness and connected in a continuous manner, improves the fit of the picked dispersion curves. The resulting profile resembles a smooth approximation of the multi‐layered one but contains, of course, less detail. As it turns out, our method does not degrade as gracefully as, for instance, diving‐wave tomography, and we can only hope to fit a subset of the dispersion curves. Therefore, the applicability of the method is limited to cases where the ratio is small and the profile is sufficiently simple. A further extension of the two‐layer model to more layers, each with a constant depth gradient of the squared slowness, might improve the fit of the modal structure but at an increased cost.     

Effects of variations in fluid properties and fracture geometry on dispersion,anisotropy, and reflection in media with planar fractures

Our objective was to discover the effect of variations in fluid properties and fracture geometry on the velocity of seismic wave propagation in fluid‐saturated media with parallel planar fractures. We used numerical models calculated by analytical solutions to examine the behaviour of P‐wave phase velocity dispersion in the normal direction to layering, in non‐porous and porous media with planar fractures. We also examined the anisotropy of low frequency phase and group velocities of fast and slow P‐waves and angular‐dependent reflection coefficients in media with planar fractures, under conditions of saturation by fluids with varying bulk moduli, densities, and fracture apertures. We defined several parametre , , and characterising dispersion, characterising anisotropy, characterising the difference between fast and slow modes, and R₀ and characterising reflection. Our results show that the behaviour of dispersion shows wider stopbands in the case of gas saturation. Concavity indicator of dispersion for gas saturation was greater than that for liquid saturation and is usually greater than one. Anisotropy is more sensitive to bulk modulus contrast than to density contrast between the solid and the fluid, and is more sensitive to density contrast than to bulk modulus contrast. The case of gas saturation usually had a greater negative R₀ and a greater value of compared with those of brine and heavy and light oil saturations. Our results are helpful in distinguishing fluid types saturating geophysical fractures and estimating the aperture and spacing of planar fractures. In seismic exploration, bulk modulus and fluid density can provide useful information in distinguishing among brine, oil, and gas; fracture geometry is important to estimate the permeability of reservoirs.     

Homotopy scattering series for seismic forward modelling with variable density and velocity

We have derived a convergent scattering series solution for the frequency-domain wave equation in acoustic media with variable density and velocity. The convergent scattering series solution is based on the homotopy analysis of a vectorial integral equation of the Lippmann–Schwinger type. By using the Green's function and partial integration, we have derived the vectorial integral equation of the Lippmann–Schwinger type that involves the pressure gradient field as well as the pressure field from the wave equation. The vectorial Lippmann–Schwinger equation can in principle be solved via matrix inversion, but the computational cost of matrix inversion scales like $N^3$, where $N$ is the number of grid blocks. The computational cost can be significantly reduced if one solves the vectorial Lippmann–Schwinger equation iteratively. A simple iterative solution is the Born series, but it is only convergent when the scattering potential is sufficiently small. In this study, we have used the so-called homotopy analysis method to derive an iterative solution for the vectorial Lippmann–Schwinger equation which can be made convergent even in strongly scattering media. The computational cost of our convergent scattering series scales as $N^2$. Our algorithm, which is based on the homotopy analysis method, involves a convergence control operator that we select using hierarchical matrices. We use a three-layer model and a resampled version of the SEG/EAGE salt model to show the performance of the developed convergent scattering series.     

Robust f‐x projection filtering for simultaneous random and erratic seismic noise attenuation

Linear prediction filters are an effective tool for reducing random noise from seismic records. Unfortunately, the ability of prediction filters to enhance seismic records deteriorates when the data are contaminated by erratic noise. Erratic noise in this article designates non‐Gaussian noise that consists of large isolated events with known or unknown distribution. We propose a robust f‐x projection filtering scheme for simultaneous erratic noise and Gaussian random noise attenuation. Instead of adopting the ?₂‐norm, as commonly used in the conventional design of f‐x filters, we utilize the hybrid ‐norm to penalize the energy of the additive noise. The estimation of the prediction error filter and the additive noise sequence are performed in an alternating fashion. First, the additive noise sequence is fixed, and the prediction error filter is estimated via the least‐squares solution of a system of linear equations. Then, the prediction error filter is fixed, and the additive noise sequence is estimated through a cost function containing a hybrid ‐norm that prevents erratic noise to influence the final solution. In other words, we proposed and designed a robust M‐estimate of a special autoregressive moving‐average model in the f‐x domain. Synthetic and field data examples are used to evaluate the performance of the proposed algorithm.     

Characterizing the degree of amplitude‐variation‐with‐offset nonlinearity in seismic physical modelling reflection data

The nonlinearity of the seismic amplitude‐variation‐with‐offset response is investigated with physical modelling data. Nonlinearity in amplitude‐variation‐with‐offset becomes important in the presence of large relative changes in acoustic and elastic medium properties. A procedure for pre‐processing physical modelling reflection data is enacted on the reflection from a water‐plexiglas boundary. The resulting picked and processed amplitudes are compared with the exact solutions of the plane‐wave Zoeppritz equations, as well as approximations that are first, second, and third order in , , and . In the low angle range of 0°–20°, the third‐order plane‐wave approximation is sufficient to capture the nonlinearity of the amplitude‐variation‐with‐offset response of a liquid‐solid boundary with , , and ρ contrasts of 1485–2745 m/s, 0–1380 m/s, and 1.00–1.19 gm/cc respectively, to an accuracy value of roughly 1%. This is in contrast to the linear Aki–Richards approximation, which is in error by as much as 25% in the same angle range. Even‐order nonlinear corrective terms are observed to be primarily involved in correcting the angle dependence of , whereas the odd‐order nonlinear terms are involved in determining the absolute amplitude‐variation‐with‐offset magnitudes.     

Research Note: Experimental measurements of Q‐contrast reflections

While seismic reflection amplitudes are generally determined by real acoustical impedance contrasts, there has been recent interest in reflections due to contrasts in seismic‐Q. Herein we compare theoretical and modelled seismic reflection amplitudes for two different cases of material contrasts. In case A, we examine reflections from material interfaces that have a large contrast in real‐valued impedance () with virtually no contrast in seismic‐Q. In case B, we examine reflections from material interfaces that have virtually no contrast in but that have very large seismic‐Q contrasts. The complex‐valued reflection coefficient formula predicts non‐zero seismic reflection amplitudes for both cases. We choose physical materials that typify the physics of both case A and case B. Physical modelling experiments show significantly large reflections for both cases – with the reflections in the two cases being phase shifted with respect to each other, as predicted theoretically. While these modelling experiments show the existence of reflections that are predicted by theory, there are still intriguing questions regarding the size of the Q‐contrast reflections, the existence of large Q‐contrast reflections in reservoir rocks and the possible application of Q‐reflection analysis to viscosity estimation in heavy oilfields.     

Mining-induced microseismic event location errors: Accuracy and precision of two location systems

The accuracy and precision of microseismic event locations were measured, analyzed, and compared for two types of location systems: anolog and digital. In the first system, relative times of first arrival were estimated from analog signals using automated hardware circuitry; station positions were estimated from mine map coordinates; and event locations were determined using the BLD (Blake, Leighton, and Duvall) direct solution method. In the second system, arrival times were manually measured during interactive displays of digital waveforms; station coordinates were surveyed; and the SW-GBM (Salamon and Wiebols; Godson, Bridges, and McKavanagh) direct basis function was used to solve for locations. Both systems assume constant isotropic seismic velocity of slightly different signals data sets, calibration blast signals with known source site and origin time, and microseismic event signals, were recorded by each location system employing the same array of high-frequency (5 kHz) accelerometers with 150 m maximum dimension. The calibration blast tests indicated a location precision of ±2 m and accuracy of ±10 m for the analog system. Location precision and accuracy for the digital system measured ±1 m and ±8 m, respectively. Numerical experiments were used to assess the contributions of errors in velocity, arrival times, and station positions on the location accuracy and precision for each system. Measured and estimated errors appropriate to each system for microseismic events were simulated in computing source locations for comparison with exact synthetic event locations. Discrepancy vectors between exact locations and locations calculated with known data errors averaged 7.7 and 1.4 m for the analog and digital systems, respectively. These averages are probably more representative of the location precision of microseismic events, since the calibration blast tests produce impulsive seismic arrivals resulting in smaller arrival-time pick errors in the analog system. For both systems, location accuracy is limited by inadequate modeling of the velocity structure. Consequently, when isotropic velocity models are used in the travel-time inversions, the increased effort expended with the digital location system does not, for the particular systems studied, result in increased accuracy.     

Experimental quantification of the seismoelectric transfer function and its dependence on conductivity and saturation in loose sand

Under certain circumstances, seismic propagation within porous media may be associated to the conversion of mechanical energy to electromagnetic energy, which is known as a seismo‐electromagnetic phenomenon. The propagation of fast compressional P‐waves is more specifically associated to the manifestations of a seismoelectric field linked to the fluid flows within the pores. The analysis of seismoelectric phenomena, which requires the combination of the theory of electrokinetics and Biot's theory of poroelasticity, provides us with transfer function that links the coseismic seismoelectric field E to the seismic acceleration . To measure the transfer function, we have developed an experimental setup enabling seismoelectric laboratory observation in unconsolidated quartz sand within the kilohertz range. The investigation focused on the impact of fluid conductivity and water saturation over the coseismic seismoelectric field. During the experiment, special attention was given to the accuracy of electric field measurements. We concluded that, to obtain a reliable estimate of the electric field amplitude, the dipole from which the potential differences are measured should be of much smaller length than the wavelength of the propagating seismic field. Time‐lapse monitoring of the seismic velocities and seismoelectric transfer functions were performed during imbibition and drainage experiments. In all cases, the quantitative analysis of the seismoelectric transfer function was in good agreement with theoretical predictions. While investigating saturation variations from full to residual water saturation, we showed that the ratio undergoes a switch in polarity at a particular saturation , which also implies a sign change of the filtration, traducing a reversal of the relative fluid displacement with respect to the frame. This sign change at critical saturation stresses a particular behaviour of the poroelastic medium: the dropping of the coseismic electric field to zero traduces the absence of relative pore/fluid displacements representative of a Biot dynamically compatible medium. We concluded from our experimental study in loose sand that the measurements of the coseismic seismoelectric coupling may provide information on fluid distribution within the pores and that the reversal of the seismoelectric field may be used as an indicator of the dynamically compatible state of the medium.     

Improved microseismic event location by inclusion of a priori dip particle motion: a case study from Ekofisk

Microseismic monitoring in petroleum settings provides insights into induced and naturally occurring stress changes. Such data are commonly acquired using an array of sensors in a borehole, providing measures of arrival times and polarizations. Events are located using 1D velocity models, P‐ and S‐wave arrival times and the azimuths of P‐wave particle motions. However in the case of all the sensors being deployed in a vertical or near‐vertical borehole, such analysis leads to an inherent 180° ambiguity in the source location. Here we present a location procedure that removes this ambiguity by using the dip of the particle motion as an a priori information to constrain the initial source location. The new procedure is demonstrated with a dataset acquired during hydraulic fracture stimulation, where we know which side of the monitoring well the events are located. Using a 5‐step location procedure, we then reinvestigate a microseismic data set acquired in April 1997 at the Ekofisk oilfield in the North Sea. Traveltimes for 2683 candidate events are manually picked. A noise‐weighted analytic‐signal polarization analysis is used to estimate the dip and azimuth of P‐wave particle motions. A modified t‐test is used to statistically assess the reliability of event location. As a result, 1462 events are located but 627 are deemed to be statistically reliable. The application of a hierarchal cluster analysis highlights coherent structures that cluster around wells and inferred faults. Most events cluster at a depth of roughly 3km in the Ekofisk chalk formation but very little seismicity is observed from the underlying Tor chalk formation, which is separated from the Ekofisk formation by an impermeable layer. We see no evidence for seismicity in the overburden but such events may be too distant to detect. The resulting picture of microseismicity at Ekofisk is very different from those presented in previous studies.     

Multi-parameter reflection waveform inversion for acoustic transversely isotropic media with a vertical symmetry axis

Full waveform inversion in transversely isotropic media with a vertical symmetry axis provides an opportunity to better match the data at the near and far offsets. However, multi-parameter full waveform inversion, in general, suffers from serious cycle-skipping and trade-off problems. Reflection waveform inversion can help us recover a background model by projecting the residuals of the reflected wavefield along the reflection wavepath. Thus, we extend reflection waveform inversion to acoustic transversely isotropic media with a vertical symmetry axis utilizing the proper parameterization for reduced parameter trade-off. From a radiation patterns analysis, an acoustic transversely isotropic media with a vertical symmetry axis is better described by a combination of the normal-moveout velocity and the anisotropic parameters η and δ for reflection waveform inversion applications. We design a three-stage inversion strategy to construct the optimal resulting model. In the first stage, we only invert for the background by matching the simulated reflected wavefield from the perturbations of and δ with the observed reflected wavefield. In the second stage, the background and η are optimized simultaneously and the far-offset reflected wavefield mainly contribute to their updates. We perform Born modelling to compute the reflected wavefield for the two stages of reflection waveform inversion. In the third stage, we perform full waveform inversion for the acoustic transversely isotropic media with a vertical symmetry axis to delineate the high-wavenumber structures. For this stage, the medium is described by a combination of the horizontal velocity , η and ε instead of , η and δ. The acoustic multi-parameter full waveform inversion utilizes the diving waves to improve the background as well as utilizes reflection for high-resolution information. Finally, we test our inversion algorithm on the modified Sigsbee 2A model (a salt free part) and a two-dimensional line from a three-dimensional ocean bottom cable dataset. The results demonstrate that the proposed reflection waveform inversion approach can recover the background model for acoustic transversely isotropic media with a vertical symmetry axis starting from an isotropic model. This recovered background model can mitigate the cycle skipping of full waveform inversion and help the inversion recover higher resolution structures.     

Some remarks on Q-compensated sparse deconvolution without knowing the quality factor Q

The subsurface media are not perfectly elastic, thus anelastic absorption, attenuation and dispersion (aka Q filtering) effects occur during wave propagation, diminishing seismic resolution. Compensating for anelastic effects is imperative for resolution enhancement. Q values are required for most of conventional Q-compensation methods, and the source wavelet is additionally required for some of them. Based on the previous work of non-stationary sparse reflectivity inversion, we evaluate a series of methods for Q-compensation with/without knowing Q and with/without knowing wavelet. We demonstrate that if Q-compensation takes the wavelet into account, it generates better results for the severely attenuated components, benefiting from the sparsity promotion. We then evaluate a two-phase Q-compensation method in the frequency domain to eliminate Q requirement. In phase 1, the observed seismogram is disintegrated into the least number of Q-filtered wavelets chosen from a dictionary by optimizing a basis pursuit denoising problem, where the dictionary is composed of the known wavelet with different propagation times, each filtered with a range of possible values. The elements of the dictionary are weighted by the infinity norm of the corresponding column and further preconditioned to provide wavelets of different values and different propagation times equal probability to entry into the solution space. In phase 2, we derive analytic solutions for estimates of reflectivity and Q and solve an over-determined equation to obtain the final reflectivity series and Q values, where both the amplitude and phase information are utilized to estimate the Q values. The evaluated inversion-based Q estimation method handles the wave-interference effects better than conventional spectral-ratio-based methods. For Q-compensation, we investigate why sparsity promoting does matter. Numerical and field data experiments indicate the feasibility of the evaluated method of Q-compensation without knowing Q but with wavelet given.     

Constrained tomography of realistic velocity models in microseismic monitoring using calibration shots

The knowledge of the velocity model in microseismic jobs is critical to achieving statistically reliable microseismic event locations. The design of microseismic networks and the limited sources for calibration do not allow for a full tomographic inversion. We propose optimizing a priori velocity models using a few active shots and a non‐linear inversion, suitable to poorly constrained systems. The considered models can be described by several layers with different P‐ and S‐wave velocities. The velocities may be constant or have 3D gradients; the layer interfaces may be simple dipping planes or more complex 3D surfaces. In this process the P‐ and S‐ wave arrival times and polarizations measured on the seismograms constitute the observed data set. They are used to estimate two misfit functions: i) one based on the measurement residuals and ii) one based on the inaccuracy of the source relocation. These two functions are minimized thanks to a simulated annealing scheme, which decreases the risk of converging to a local solution within the velocity model. The case study used to illustrate this methodology highlights the ability of this technique to constrain a velocity model with dipping layers. This was performed by jointly using sixteen perforation shots recorded during a multi‐stage fracturing operation from a single string of 3C‐receivers. This decreased the location inaccuracies and the residuals by a factor of six. In addition, the retrieved layer dip was consistent with the pseudo‐horizontal trajectories of the wells and the background information provided by the customer. Finally, the theoretical position of each calibration shot was contained in the uncertainty domain of the relocation of each shot. In contrast, single‐stage inversions provided different velocity models that were neither consistent between each other nor with the well trajectories. This example showed that it is essential to perform a multi‐stage inversion to derive a better updated velocity model.     

The rotated Cartesian coordinate method to remove the axial singularity of cylindrical coordinates in finite‐difference schemes for elastic and viscoelastic waves

When modelling the propagation of 3D non‐axisymmetric elastic and viscoelastic waves in cylindrical coordinates using the finite‐difference time‐domain method, a mathematical singularity occurs due to the presence of terms in the elastic and viscoelastic wave equations. For many years, this issue has been impeding the accurate numerical solution near the axis. In this work, we propose a simple but effective method for the treatment of this numerical singularity problem. By rotating the Cartesian coordinate system around the z‐axis in cylindrical coordinates, the numerical singularity problems in both 2D and 3D cylindrical coordinates can be removed. This algorithm has three advantages over the conventional treatment techniques: (i) the excitation source can be directly loaded at , (ii) the central difference scheme with second‐order accuracy is maintained, and (iii) the stability condition at the axis is consistent with the finite‐difference time‐domain in Cartesian coordinates. This method is verified by several 3D numerical examples. Results show that the rotating the Cartesian coordinate method is accurate and stable at the singularity axis. The improved finite‐difference time‐domain algorithm is also applied to sonic logging simulations in non‐axisymmetric formations and sources.