A reservoir wave velocity variation imaging method suitable for non-uniform scattering media

Abstract

The application discloses a reservoir wave velocity variation imaging method suitable for non-uniform scattering medium and relates to the technical field of geophysical monitoring. Firstly, continuous seismic waveforms are collected, tail wave interference is used to extract time shift and construct an observation data set; based on radiation transmission theory, a zero collision differential Monte Carlo is used to calculate a global sensitive kernel and construct a sensitive kernel matrix; a plurality of Voronoi dissection models are adaptively generated according to the sensitive kernel distribution to realize parameter dimension reduction; high-dimensional sensitive kernels are mapped to a low-dimensional space through a projection matrix, and a plurality of medium models are obtained through damped inversion; imaging results are obtained by averaging the plurality of models, and the uncertainty of the results is evaluated by standard deviation. The application is adapted to the non-uniform scattering characteristics of the reservoir, weakens the ill-conditioned and grid effect of traditional grid inversion, improves the imaging resolution and stability, and is suitable for reservoir dynamic monitoring such as hydraulic fracturing, geothermal energy and carbon storage.

Description

Technical Field

[0001] This invention relates to the field of geophysical monitoring technology, and more specifically to a reservoir wave velocity variation imaging method applicable to non-uniform scattering media. Background Technology

[0002] With the continuous development of unconventional oil and gas resources, hydraulic fracturing has become an important technical means for enhancing the production of unconventional reservoirs such as coalbed methane and shale gas. During fracturing, a large amount of high-pressure fluid is injected into the underground reservoir, causing a series of complex geophysical responses, including fracture propagation, stress redistribution, and pore fluid migration. Therefore, real-time and accurate monitoring of changes in the properties of the underground medium during fracturing is of great significance for evaluating the fracturing effect, optimizing construction parameters, and ensuring the safe development of the reservoir.

[0003] Existing fracturing monitoring technologies mainly include microseismic location, seismic tomography, and repeated seismic monitoring. Microseismic monitoring can reflect the fracture propagation process, but its results are highly dependent on the accuracy of event location and the characteristics of the seismic source distribution. While seismic tomography can reconstruct the subsurface medium structure, it typically requires high-density X-ray coverage and tends to produce low-resolution imaging results under limited observation conditions at the fracturing site. In contrast, codawave interferometry (CWI) utilizes the information from multiple scattering of seismic waves in complex media. By measuring the changes in wake time before and after disturbance, it estimates the changes in subsurface medium wave velocity. It exhibits extremely high sensitivity to minute media disturbances and has been widely applied in recent years in fields such as volcanic activity monitoring, fault zone evolution research, groundwater change monitoring, and geothermal reservoir dynamic monitoring.

[0004] However, traditional wake wave interferometry typically assumes a uniform distribution of medium changes within the study area. The results obtained only reflect the average effect of wave velocity changes within the observation area, making it difficult to accurately describe the spatial distribution characteristics of local fracture development zones and fluid-affected zones during fracturing. To achieve spatial variation imaging of subsurface media, researchers have proposed using wake wave sensitive kernels to establish the relationship between observed measurements and subsurface medium disturbances, and to recover the spatial distribution of wave velocity changes through inversion methods. Currently, the most widely used sensitive kernels include diffusion approximation sensitive kernels and radiative transfer sensitive kernels. While diffusion approximation sensitive kernels are suitable for strongly scattering homogeneous media conditions, in actual fracturing reservoirs, the development of fractures, bedding, and structures leads to significant spatial non-uniformity in the scattering medium, making it difficult for diffusion approximation models to accurately describe the wake wave propagation process. In contrast, sensitive kernels based on radiative transfer theory can more realistically characterize the energy propagation characteristics in complex scattering media, exhibiting higher physical accuracy.

[0005] However, the calculation of radiative transfer sensitive nuclei involves a large number of random scattering path integrals, resulting in high computational costs. When the monitoring area is large or the inversion grid is dense, traditional Monte Carlo methods often require a large number of particle simulations to obtain stable results, leading to low computational efficiency and making it difficult to meet the needs of rapid imaging in actual fracturing monitoring. In addition, existing wake imaging methods usually use fixed regular grid parameterized subsurface models. Due to the uneven distribution of actual observation stations, limited observation coverage, and uneven spatial distribution of sensitive nuclei, fixed grid inversion is prone to problems such as excessively high parameter dimensionality, enhanced ill-conditionedness, and unstable imaging results, thus affecting the accuracy of restoring the spatial location and variation amplitude of subsurface anomalies.

[0006] Therefore, how to propose a reservoir wave velocity change imaging method suitable for non-uniform scattering media, which can improve the computational efficiency of the sensitive core, achieve adaptive optimization of the subsurface model parameterization, and improve the stability and resolution of subsurface medium wave velocity change imaging during fracturing, thus providing a more reliable technical means for dynamic monitoring of hydraulic fracturing, is a problem that urgently needs to be solved by those skilled in the art. Summary of the Invention

[0007] In view of this, the present invention provides a reservoir wave velocity change imaging method applicable to non-uniform scattering media, which makes full use of the wake multiple scattering information to improve the stability and spatial resolution of underground medium change imaging under non-uniform scattering media conditions.

[0008] To achieve the above objectives, the present invention adopts the following technical solution: An imaging method for reservoir wave velocity variations applicable to non-uniform scattering media includes the following steps: S1. Collect continuous waveform data recorded by multiple receiving points within the monitoring area, preprocess the raw data, and construct a reference waveform and the current waveform; S2. Calculate the wake time shift information between the reference waveform and the current waveform using the wake interferometry method to obtain observation data; S3. The zero-collision differential Monte Carlo method is used to simulate the propagation process of the wake wave in the non-uniform scattering medium, calculate the radiation transfer sensitive nuclei corresponding to each location in the study area, and perform statistical analysis on the radiation transfer sensitive nuclei to obtain the overall sensitive nuclei distribution; S4. Combine the radiation transfer sensitive kernels obtained under multiple receiving points, multiple time windows, and multiple virtual seismic sources to construct a sensitive kernel matrix; S5. Determine the spatial distribution probability of Voronoi kernel points based on the overall sensitive kernel distribution, and generate multiple random Voronoi subdivision models within the study area; S6. Using the sensitive kernel matrix and observation data, inversion solutions are performed under multiple Voronoi subdivision models to obtain multiple sets of subsurface medium change models; S7. Statistically average multiple subsurface medium change models and calculate the corresponding standard deviations to obtain the final subsurface medium change imaging results and its stability evaluation results.

[0009] Preferably, S2 includes: S21. Extract multiple non-overlapping wake time windows, and perform sliding cross-correlation calculations on the reference waveform and the current waveform within the same time window. The cross-correlation formula is: ; In the formula, For cross-correlation coefficients, For reference waveform, The current waveform after time shift The waveform after sliding, This is the center moment of the wake window; S22. Take the time shift corresponding to the maximum value of the cross-correlation coefficient. The wake time shift of the current time window Then, the relative change in wave velocity is calculated. ; S23. Summarize the wake time shifts corresponding to all virtual sources, receivers, and wake time windows. Relative wave velocity change This forms an observation dataset.

[0010] Preferably, S3 includes: S31. Establish a non-uniform scattering medium model, and given the spatially varying scattering coefficient and absorption coefficient, calculate the extinction coefficient; S32. Using the extinction coefficient, the zero-collision differential Monte Carlo method is employed to simulate the propagation process of the wake wave in a non-uniform scattering medium model, and the contribution of the random walk path of the particles to the observed quantities is statistically analyzed to obtain the radiative transfer sensitive kernel. S33. Map all radiation transfer sensitive nuclei onto the spatial grid of the study area according to spatial coordinates, and accumulate and statistically analyze the sensitivity values ​​in each grid cell to obtain the overall sensitive nuclei distribution.

[0011] Preferably, in S31, the spatial position Extinction coefficient at The calculation formula is: ; In the formula, The scattering coefficient is... is the absorption coefficient.

[0012] Preferably, S32 includes: S321. Extract the upper bound of the global extinction coefficient of the study area. The distance of free propagation of particles is obtained by random number sampling. ; S322. The particle starts from the virtual hypocenter and travels a free propagation distance along a random propagation direction. Upon reaching the target mesh, read the extinction coefficient of the target mesh and compare it with the upper bound of the global extinction coefficient. The collision type is determined by comparison. If the true collision probability condition is met, the particle scattering propagation direction is updated. If it is determined to be a zero collision event, the particle continues to move forward with its original propagation direction. S323. Repeatedly execute the free path sampling, particle propagation, and collision detection process until the total particle propagation time reaches the preset wake time window duration; S324. Track all grid cells that the particle passes through throughout its entire path, and based on the radiative transfer response relationship between changes in the underground medium and the time shift of the wake wave, accumulate the particle path contribution weight corresponding to each passing grid cell; S325. Filter the effective particles that finally arrive at the receiving point, repeatedly release a large number of particles to complete the global path statistics, normalize the cumulative contribution weight of each grid unit, and output the radiation transfer sensitive kernel under the current virtual source, receiving point, and wake time window matching.

[0013] Preferably, the radiative transfer response relationship between changes in the underground medium and wake time shift is as follows: ; In the formula, In order to provide a virtual earthquake source Receiving point Total propagation time of the wake wave Under certain conditions, spatial location Radiative transfer sensitive nuclei at the location; For particles from virtual seismic sources Duration Along direction propagation to spatial location The probability of; For particles from spatial position Remaining time Along the reverse direction Transmit to receiving point The probability, For particles from virtual seismic sources Total propagation time of the wake Successfully transmitted to the receiving point The total propagation probability; For spatial dimensions The unit spherical surface area in the middle, For the intermediate propagation time integral variable, The unit vector is the direction of propagation.

[0014] Preferably, S4 includes: The radiation transfer sensitive kernels corresponding to each virtual source, each receiving point, and each wake time window are expanded into one-dimensional vectors in grid order, and all one-dimensional vectors are spliced ​​together to form a sensitive kernel matrix.

[0015] Preferably, S5 includes: A probability density function is constructed based on the overall sensitive kernel distribution, and Voronoi kernel points are randomly generated within the study area according to the probability density function; the generation probability of Voronoi kernel points is related to the overall sensitive kernel intensity. With each Voronoi core point as the center, Voronoi cells are divided according to the nearest neighbor principle, and multiple grid cells within each Voronoi cell share the same model parameters. Voronoi kernel points were generated independently and randomly multiple times, and spatial partitioning was performed to obtain multiple sets of stochastic Voronoi partitioning models with structural differences.

[0016] Preferably, S6 includes: Establish a projection mapping relationship between the regular mesh model and the Voronoi model, and based on the projection mapping relationship, apply the sensitive kernel matrix... Mapped to the parameter space corresponding to each Voronoi partitioning model; Based on the mapped sensitive kernel matrix and the observed data, damped constraint inversion solutions are performed under each Voronoi partitioning model to obtain multiple sets of subsurface medium change models.

[0017] As can be seen from the above technical solution, compared with the prior art, the present invention discloses a reservoir wave velocity variation imaging method suitable for non-uniform scattering media. By introducing the zero-collision differential Monte Carlo sensitive kernel calculation method and the sensitive kernel-constrained Voronoi adaptive parameterized inversion strategy, it can effectively solve the problems of difficult calculation of tail wave imaging sensitivity, excessive inversion dimension and insufficient imaging stability in complex non-uniform scattering media, and improve the reliability and spatial resolution of underground medium anomaly recovery. Attached Figure Description

[0018] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on the provided drawings without creative effort.

[0019] Figure 1 A flowchart of the method provided by the present invention; Figure 2 A comparison diagram of waveforms before and after disturbance in a multi-scattering medium; Figure 3This is a diagram of the overall structure of the sensitive kernel matrix; Figure 4 This is a map showing the overall spatial distribution of sensitive nuclei. Figure 5 Voronoi subdivision and single inversion diagram; Figure 6 This is an image showing the dt imaging results of the tail wave travel time change during fracturing based on the reference waveform before fracturing; Figure 7 The image shows the relative travel time variation dt / t during fracturing based on the reference waveform before fracturing. Figure 8 The image shows the dt imaging results of the tail wave travel time change after fracturing based on the reference waveform before fracturing. Figure 9 The results are three-stage Voronoi imaging based on a full-time reference waveform; Figure 10 This is a diagram showing the results of a three-stage conventional inversion. Detailed Implementation

[0020] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0021] refer to Figure 1 This embodiment discloses a reservoir wave velocity variation imaging method applicable to non-uniform scattering media, including the following steps: S1. Collect continuous waveform data recorded by multiple receiving points within the monitoring area, preprocess the raw data, and construct a reference waveform and the current waveform.

[0022] Preprocessing includes waveform mean removal, detrending, bandpass filtering, and time synchronization correction.

[0023] S2. Calculate the wake time shift information between the reference waveform and the current waveform using the wake interferometry method to obtain observation data, including: S21. Extract multiple non-overlapping wake time windows, and perform sliding cross-correlation calculations on the reference waveform and the current waveform within the same time window. The cross-correlation formula is: ; In the formula, For cross-correlation coefficients, For reference waveform, The current waveform after time shift The waveform after sliding, This is the center moment of the wake window; S22. Take the time shift corresponding to the maximum value of the cross-correlation coefficient. The wake time shift of the current time window ,Right now Then according to the formula The relative change in wave velocity was calculated. ; S23. Summarize the wake time shifts corresponding to all virtual sources, receivers, and wake time windows. Relative wave velocity change This forms an observation dataset.

[0024] like Figure 2 As shown, when the properties of the underground medium change, the direct wave is only affected along a limited propagation path, while the wake wave, due to undergoing a large number of scattering processes, has a wider propagation path coverage and is therefore more sensitive to small changes in the medium. Figure 2 The left side of the image shows the wake waveform obtained during the reference period, while the right side shows the wake waveform obtained during the monitoring period. When the subsurface medium changes, the wake waveform will exhibit a slight time shift. The time shift within different time windows is calculated using a sliding time window cross-correlation method, thus obtaining the wake time shift observation data.

[0025] S3. The zero-collision differential Monte Carlo method is used to simulate the propagation process of the wake wave in a non-uniform scattering medium. The radiation transfer sensitive nuclei at each location within the study area are calculated, and statistical analysis is performed on the radiation transfer sensitive nuclei to obtain the overall sensitive nuclei distribution, including: S31. Establish a model of a non-uniform scattering medium, given a spatially varying scattering coefficient. With absorption coefficient Calculate the extinction coefficient.

[0026] Spatial location Extinction coefficient at The calculation formula is: .

[0027] S32. Using the extinction coefficient, the zero-collision differential Monte Carlo method is employed to simulate the propagation process of the wake wave in a non-uniform scattering medium model. The contribution of the random walk path of particles to the observed quantities is statistically analyzed to obtain the radiative transfer sensitive kernel.

[0028] To avoid directly solving complex radiative transfer equations in non-homogeneous media, this embodiment employs a zero-collision method for propagation simulation, including: S321. Extract the upper bound of the global extinction coefficient of the study area. The distance of free propagation of particles is obtained by random number sampling. .

[0029] ;in, is a uniform random number in the interval (0, 1).

[0030] S322. The particle starts from the virtual hypocenter and travels a free propagation distance along a random propagation direction. Upon reaching the target mesh, read the extinction coefficient of the target mesh and compare it with the upper bound of the global extinction coefficient. The collision type is determined by comparison (specifically, the probability of accepting a real collision is calculated based on the ratio between the local scattering extinction term at the current position and the global maximum extinction coefficient, and whether the candidate collision is converted into a real collision event is determined by random sampling). If the real collision probability condition is met, the particle scattering propagation direction is updated according to the local scattering law; when no real collision occurs, it is recorded as a zero collision event, and the particle continues to move forward with its original propagation direction.

[0031] In the zero-collision differential Monte Carlo simulation, to avoid the problem of a large number of path resampling caused by the change of scattering coefficient in the non-uniform scattering medium in space, the traditional Monte Carlo method adopts a zero-collision sampling strategy based on the upper bound extinction coefficient.

[0032] First, based on the scattering coefficient, absorption coefficient, and additional attenuation term at each spatial location within the study area, the maximum extinction coefficient over the entire computational domain is determined:

[0033] in: The absorption coefficient at spatial location (x, y) is represented. The scattering coefficient at spatial location (x, y) represents the scattering coefficient. This represents the virtual collision coefficient used to satisfy the upper bound condition. It is the maximum extinction coefficient that satisfies the extinction coefficient constraints at all locations within the entire model region.

[0034] During particle propagation, instead of directly sampling the free path based on the actual extinction coefficient at the current position, a unified upper bound extinction coefficient is used to sample candidate collision positions. For each candidate collision position... First, calculate the true extinction coefficient based on the medium parameters at that location:

[0035] Then, the probability of actual collision acceptance is calculated:

[0036] Generate uniformly distributed random numbers between 0 and 1:

[0037] Compare the random number with the collision probability: when < When a real collision event occurs at the current location, the particle is considered to enter a real scattering or absorption process; when ≥ If the current position is considered to be a zero-collision event, the particle continues to propagate along its original direction without changing its direction of propagation; only the particle weight is updated according to the weight correction formula.

[0038] Furthermore, after a real collision event occurs, the collision type is determined based on the scattering ratio:

[0039] When ξs < If ξs is a new random number, it is determined to be a scattering event, and the propagation direction is regenerated according to the scattering phase function; otherwise, it is determined to be a non-scattering attenuation event, and the particle energy weight is updated according to the absorption coefficient.

[0040] Therefore, the zero-collision differential Monte Carlo method in this embodiment is not a simple random judgment, but rather a probabilistic simulation of the radiation transmission path in a spatially non-uniform scattering medium: (1) Establish candidate collision processes using the maximum extinction coefficient across the entire region; (2) Determine the true collision acceptance probability by using the ratio of the local true extinction coefficient to the maximum extinction coefficient; (3) Further determine the scattering behavior by using the ratio of the scattering coefficient to the true extinction coefficient.

[0041] S323. Execute the free path sampling, particle propagation, and collision detection process repeatedly until the total particle propagation time reaches the preset wake time window duration.

[0042] S324. Track all grid cells along the entire path of the particle, and based on the radiative transfer response relationship between changes in the underground medium and the time shift of the wake wave, accumulate the particle path contribution weight corresponding to each passing grid.

[0043] The radiative transfer response relationship between changes in the underground medium and wake time shift is as follows: ; In the formula, In order to provide a virtual earthquake source Receiving point Total propagation time of the wake wave Under certain conditions, spatial location Radiative transfer sensitive nuclei at the location; For particles from virtual seismic sources Duration Along direction propagation to spatial location The probability of; For particles from spatial position Remaining time Along the reverse direction Transmit to receiving point The probability, For particles from virtual seismic sources Total propagation time of the wake Successfully transmitted to the receiving point The total propagation probability; For spatial dimensions The unit spherical surface area in the middle, For the intermediate propagation time integral variable, The unit vector is the direction of propagation.

[0044] Sensitive kernel definition In the middle, molecules Indicates the position the particle passed through The probability density, denominator This represents the total probability density of particles from the source to the receiver; the ratio of the two is the weight of the influence of that location on the wake time shift.

[0045] S325. Filter the effective particles that finally arrive at the receiving point, repeatedly release a large number of particles to complete the global path statistics, normalize the cumulative contribution weight of each grid unit, and output the radiation transfer sensitive kernel under the current virtual source, receiving point, and wake time window matching.

[0046] S33. Map all radiation transfer sensitive nuclei onto the spatial grid of the study area according to spatial coordinates, and accumulate and statistically analyze the sensitivity values ​​in each grid cell to obtain the overall sensitive nuclei distribution.

[0047] S4. Combine the radiation transfer sensitive kernels obtained under multiple receiving points, multiple time windows and multiple virtual seismic sources to construct a sensitive kernel matrix.

[0048] Specifically, the radiation transfer sensitive cores corresponding to each virtual seismic source, each receiving point, and each wake time window are expanded into one-dimensional vectors in grid order, and all one-dimensional vectors are spliced ​​together to form a sensitive core matrix.

[0049] Figure 3 The overall structure of the constructed sensitive kernel matrix is ​​presented. As can be seen from the figure, the matrix exhibits a distinct striped and blocky structure, rather than a random distribution. Significant differences in response intensity exist at different receiver points and time windows, indicating that the observation system demonstrates significant spatial selectivity in sampling the subsurface medium.

[0050] To further analyze the spatial constraint capability of the observation system, statistical analysis was performed on all sensitive nuclei. Figure 4The overall spatial distribution of sensitive cores is presented. The black line represents the trajectory of the estimated horizontal wells, and each small dot represents the location of the deployed seismic stations. The same applies to each subsequent figure. Figure 4 The high-value areas are mainly concentrated near the horizontal well trajectory and gradually decrease outwards. This result indicates that wake observations have the strongest constraint on the well perimeter area, while the recovery ability of the outer area is relatively weak due to observation coverage limitations. Therefore, this difference in spatial resolution needs to be fully considered in subsequent parameterization processes.

[0051] S5. Determine the spatial distribution probability of Voronoi kernel points based on the overall sensitive kernel distribution, and generate multiple random Voronoi partitioning models within the study area, including: S51. Construct a probability density function based on the overall sensitive kernel distribution, and randomly generate Voronoi kernel points within the study area according to the probability density function; wherein the generation probability of Voronoi kernel points is related to the overall sensitive kernel intensity.

[0052] In this embodiment, the overall sensitivity kernel is normalized and a probability density function is constructed accordingly. Higher sampling probabilities are assigned to high-sensitivity regions and lower sampling probabilities are assigned to low-sensitivity regions. Then, Voronoi kernel points are randomly generated within the model region.

[0053] S52. Using each Voronoi core point as the center, divide the Voronoi cells according to the nearest neighbor principle. Multiple grid cells within each Voronoi cell share the same model parameters.

[0054] Voronoi units are constructed based on the nearest neighbor principle, satisfying: ; In the formula, For the first One Voronoi core point; For the first One Voronoi core point; For the first One Voronoi unit; It is a measure of spatial distance.

[0055] All mesh elements within each Voronoi cell share the same medium model parameters.

[0056] S53. Voronoi kernel points are generated independently and randomly multiple times, and spatial partitioning is completed to obtain multiple sets of stochastic Voronoi partitioning models with structural differences.

[0057] To reduce the uncertainty of the results caused by random partitioning, this embodiment sets 20 independent random Voronoi partitions for subsequent parallel inversion.

[0058] S6. Performing inversion calculation using the sensitivity kernel matrix and observation data under multiple Voronoi tessellation models respectively to obtain multiple sets of underground medium variation models, including: S61. Establishing the projection mapping relationship between the regular grid model and the Voronoi model, and based on the projection mapping relationship, the sensitivity kernel matrix is mapped to the parameter space corresponding to each Voronoi tessellation model.

[0059] S62. Converting the high-dimensional sensitivity matrix into a low-dimensional design matrix based on the Voronoi mapping matrix, establishing a linear relationship between the observation data and the low-dimensional model parameters, constructing an objective function by adding a damping constraint term, and solving the wave velocity variation parameters corresponding to each Voronoi model using a regularized least squares method.

[0060] After completing the calculation of the sensitivity kernel matrix, the sensitivity matrix G under the high-resolution spatial grid is first obtained, in which it represents the sensitivity degree of the j-th spatial grid unit corresponding to the i-th observation time window to the wave velocity variation. Since the high-resolution grid will lead to an excessive number of unknown parameters, this embodiment further adopts a Voronoi adaptive parameterization method to reduce the dimension of model parameters.

[0061] For any Voronoi tessellation model, the spatial grid is divided into N Voronoi units: V k (k=1,2,...,N), it is assumed that the wave velocity variation parameter remains consistent within the same Voronoi unit: m(x,y)=m k . Therefore, the original high-dimensional model parameter is mapped to a low-dimensional parameter , where N<M. This is realized through the Voronoi mapping matrix P:

[0062] Meanwhile, the sensitivity matrix is converted into a low-dimensional design matrix:

[0063] Subsequently, a linear inversion model is established according to the time variation observation quantity d obtained by coda wave interferometry:

[0064] Considering the non-uniqueness of reservoir monitoring inversion and the influence of noise, a damping constraint term is introduced to construct an objective function:

[0065] Where: the first term ensures that the inversion results match the actual observation data; the second term limits excessive changes in model parameters and improves inversion stability; and λ is the damping coefficient. The wave velocity variation distribution corresponding to the current Voronoi model is obtained by solving for the first-order minimum of the objective function.

[0066] Since different initial Voronoi node positions will produce different spatial parameterization results, this invention randomly generates multiple Voronoi subdivision models P1, P2, ..., P n Calculate the corresponding inversion results m1, m2, ..., m respectively. n .

[0067] Figure 5 Typical Voronoi subdivision results and corresponding single-pass inversion results are presented. As can be seen from the figures, the Voronoi cells are significantly denser near the well trajectory, while becoming sparser in the outer regions. This parameterization method allows the model's degrees of freedom to match the observational constraints.

[0068] S7. Statistically average multiple subsurface medium change models and calculate the corresponding standard deviations to obtain the final subsurface medium change imaging results and its stability evaluation results.

[0069] This embodiment yielded a total of 20 independent inversion models. The arithmetic mean of all models is used to obtain the ensemble average model, which serves as the final imaging result of subsurface media changes.

[0070] This embodiment employs a multi-model averaging strategy to reduce spatial bias caused by a single parameterization method, thereby improving the stability of anomaly region location in complex fracturing reservoirs.

[0071] In the formula, The total number of groups in a random Voronoi partition; For the first The inversion model corresponding to the group subdivision; For the integrated average model.

[0072] The dispersion of each group of inversion results is quantified by the sample standard deviation, which is used to evaluate the uncertainty and stability of the imaging results. ; In the formula, This represents the standard deviation field corresponding to each inversion model group; the larger the standard deviation value, the higher the uncertainty of the imaging result at that location.

[0073] Figure 6The average inversion results for the fracturing stage are presented, using the pre-fracturing superimposed waveform as the reference waveform. As can be seen from the figure, the main anomaly areas are concentrated near the well trajectory and exhibit a continuous banded distribution. The locally high anomaly areas correspond to the regions with the strongest fracturing action, reflecting the wake response caused by fracture propagation and changes in medium properties during the fracturing process.

[0074] To verify the stability of the results, a control inversion was further performed using the normalized observation dt / t. Figure 7 The corresponding dt / t inversion results are given. Figure 6 In comparison, the anomaly locations and overall spatial patterns are basically consistent, with only a slight decrease in local amplitude. This result indicates that the main anomaly regions recovered by this invention do not depend on the selection of the observation scale and have good stability and robustness.

[0075] The same method was then used to invert the monitoring data after the fracturing was completed. Figure 8 Imaging results are presented for the post-fracturing stage. Compared to the fracturing stage, the overall anomalous amplitude is reduced, and the spatial distribution is smoother, but the main anomalous areas still maintain continuity. This indicates that the effects of fractures and pore pressure disturbances do not immediately disappear after fracturing, but rather form a persistent response in the wellbore periphery.

[0076] To further analyze the temporal evolution of subsurface media changes, a full-time superimposed waveform was used as a unified reference waveform. Imaging was performed at three stages: before fracturing, during fracturing, and after fracturing. It is important to note that the physical meaning of these results differs from that of wake imaging based on the pre-fracturing reference waveform. The former reflects the deviation of each stage from the long-term average state, while the latter reflects the changes in each stage relative to the pre-fracturing background state.

[0077] Figure 9 The results of three-stage Voronoi adaptive imaging are presented. The overall anomaly is weak in the pre-fracturing stage, close to the background state; the anomaly is strongest in the mid-fracturing stage, mainly concentrated near the well trajectory; the anomaly weakens in the post-fracturing stage but still remains continuously distributed, indicating that the subsurface medium has undergone an evolution process from the background state to the strongly disturbed state and then to the residual response state.

[0078] To verify the superiority of this invention, a comparative inversion was performed using a traditional regular mesh parameterization method. Figure 10 The corresponding conventional inversion results are given. (Compared to...) Figure 9 In contrast, conventional inversion exhibits a significant grid effect, with coarse anomaly boundaries and discrete artifacts appearing in local areas. However, this invention employs sensitive kernel-constrained Voronoi adaptive parameterization, resulting in more continuous anomaly boundaries, a more stable spatial structure, and a more accurate recovery of the subsurface medium variation characteristics in the fracturing-affected zone.

[0079] The above embodiments demonstrate that the present invention can make full use of wake wave multiple scattering information and combine non-uniform scattering sensitive kernel with Voronoi adaptive parameterization method to achieve stable imaging of spatial changes in underground media, which is particularly suitable for scenarios such as hydraulic fracturing monitoring, geothermal reservoir monitoring and underground engineering safety monitoring.

[0080] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on its differences from other embodiments. Similar or identical parts between embodiments can be referred to interchangeably. For the apparatus disclosed in the embodiments, since they correspond to the methods disclosed in the embodiments, the description is relatively simple; relevant parts can be referred to the method section.

[0081] The above description of the disclosed embodiments enables those skilled in the art to make or use the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A reservoir wave velocity variation imaging method suitable for non-uniform scattering media, characterized in that, Includes the following steps: S1. Collect continuous waveform data recorded by multiple receiving points within the monitoring area, preprocess the raw data, and construct a reference waveform and the current waveform; S2. Calculate the wake time shift information between the reference waveform and the current waveform using the wake interferometry method to obtain observation data; S3. The zero-collision differential Monte Carlo method is used to simulate the propagation process of the wake wave in the non-uniform scattering medium, calculate the radiation transfer sensitive nuclei corresponding to each location in the study area, and perform statistical analysis on the radiation transfer sensitive nuclei to obtain the overall sensitive nuclei distribution; S4. Combine the radiation transfer sensitive kernels obtained under multiple receiving points, multiple time windows, and multiple virtual seismic sources to construct a sensitive kernel matrix; S5. Determine the spatial distribution probability of Voronoi kernel points based on the overall sensitive kernel distribution, and generate multiple random Voronoi subdivision models within the study area; S6. Using the sensitive kernel matrix and observation data, inversion solutions are performed under multiple Voronoi subdivision models to obtain multiple sets of subsurface medium change models; S7. Statistically average multiple subsurface medium change models and calculate the corresponding standard deviations to obtain the final subsurface medium change imaging results and its stability evaluation results.

2. The reservoir wave velocity variation imaging method applicable to non-uniform scattering media according to claim 1, characterized in that, S2 include: S21. Extract multiple non-overlapping wake time windows, and perform sliding cross-correlation calculations on the reference waveform and the current waveform within the same time window. The cross-correlation formula is: ; In the formula, For cross-correlation coefficients, For reference waveform, The current waveform after time shift The waveform after sliding, This is the center moment of the wake window; S22. Take the time shift corresponding to the maximum value of the cross-correlation coefficient. The wake time shift of the current time window Then, the relative change in wave velocity is calculated. ; S23. Summarize the wake time shifts corresponding to all virtual sources, receivers, and wake time windows. Relative wave velocity change This forms an observation dataset.

3. The reservoir wave velocity variation imaging method applicable to non-uniform scattering media according to claim 1, characterized in that, S3 include: S31. Establish a non-uniform scattering medium model, and given the spatially varying scattering coefficient and absorption coefficient, calculate the extinction coefficient; S32. Using the extinction coefficient, the zero-collision differential Monte Carlo method is employed to simulate the propagation process of the wake wave in a non-uniform scattering medium model, and the contribution of the random walk path of the particles to the observed quantities is statistically analyzed to obtain the radiative transfer sensitive kernel. S33. Map all radiation transfer sensitive nuclei onto the spatial grid of the study area according to spatial coordinates, and accumulate and statistically analyze the sensitivity values ​​in each grid cell to obtain the overall sensitive nuclei distribution.

4. The reservoir wave velocity variation imaging method applicable to non-uniform scattering media according to claim 3, characterized in that, In S31, spatial location Extinction coefficient at The calculation formula is: ； In the formula, The scattering coefficient is... is the absorption coefficient.

5. The reservoir wave velocity variation imaging method applicable to non-uniform scattering media according to claim 4, characterized in that, S32 includes: S321. Extract the upper bound of the global extinction coefficient of the study area. The distance of free propagation of particles is obtained by random number sampling. ; S322. The particle starts from the virtual hypocenter and travels a free propagation distance along a random propagation direction. Upon reaching the target mesh, read the extinction coefficient of the target mesh and compare it with the upper bound of the global extinction coefficient. The collision type is determined by comparison. If the true collision probability condition is met, the particle scattering propagation direction is updated. If it is determined to be a zero collision event, the particle continues to move forward with its original propagation direction. S323. Repeatedly execute the free path sampling, particle propagation, and collision detection process until the total particle propagation time reaches the preset wake time window duration; S324. Track all grid cells that the particle passes through throughout its entire path, and based on the radiative transfer response relationship between changes in the underground medium and the time shift of the wake wave, accumulate the particle path contribution weight corresponding to each passing grid cell; S325. Filter the effective particles that finally arrive at the receiving point, repeatedly release a large number of particles to complete the global path statistics, normalize the cumulative contribution weight of each grid unit, and output the radiation transfer sensitive kernel under the current virtual source, receiving point, and wake time window matching.

6. The reservoir wave velocity variation imaging method applicable to non-uniform scattering media according to claim 5, characterized in that, The radiative transfer response relationship between changes in the underground medium and wake time shift is as follows: ； In the formula, In order to provide a virtual earthquake source Receiving point Total propagation time of the wake wave Under certain conditions, spatial location Radiative transfer sensitive nuclei at the location; For particles from virtual seismic sources Duration Along direction propagation to spatial location The probability of; For particles from spatial position Remaining time Along the reverse direction Transmit to receiving point The probability, For particles from virtual seismic sources Total propagation time of the wake Successfully transmitted to the receiving point The total propagation probability; For spatial dimensions The unit spherical surface area in the text For the intermediate propagation time integral variable, The unit vector is the direction of propagation.

7. The reservoir wave velocity variation imaging method applicable to non-uniform scattering media according to claim 1, characterized in that, S4 include: The radiation transfer sensitive kernels corresponding to each virtual source, each receiving point, and each wake time window are expanded into one-dimensional vectors in grid order, and all one-dimensional vectors are spliced ​​together to form a sensitive kernel matrix.

8. The reservoir wave velocity variation imaging method applicable to non-uniform scattering media according to claim 1, characterized in that, S5 include: Construct a probability density function based on the overall sensitive kernel distribution, and randomly generate Voronoi kernel points within the study area according to the probability density function; The generation probability of Voronoi kernel points is related to the overall sensitive kernel strength; With each Voronoi core point as the center, Voronoi cells are divided according to the nearest neighbor principle, and multiple grid cells within each Voronoi cell share the same model parameters. Voronoi kernel points were generated independently and randomly multiple times, and spatial partitioning was performed to obtain multiple sets of stochastic Voronoi partitioning models with structural differences.

9. The reservoir wave velocity variation imaging method applicable to non-uniform scattering media according to claim 1, characterized in that, S6 include: Establish a projection mapping relationship between the regular mesh model and the Voronoi model, and based on the projection mapping relationship, apply the sensitive kernel matrix... Mapped to the parameter space corresponding to each Voronoi partitioning model; Based on the mapped sensitive kernel matrix and the observed data, damped constraint inversion solutions are performed under each Voronoi partitioning model to obtain multiple sets of subsurface medium change models.

Images