Fracturing crack real-time monitoring method and system based on time-varying electric field dynamic data

Abstract

The invention provides a fracturing crack real-time monitoring method and system based on time-varying electric field dynamic data, and the method specifically comprises the steps: obtaining an electric field data sequence which changes with time and is collected by an electric field sensor disposed in a fracturing work area; for each monitoring time step k in the electric field data sequence, executing the following sequential Bayesian inversion: a, constructing state parameter prior distribution of the time step k by adopting kernel density estimation based on a sample of crack state parameter posteriori distribution of the time step k-1; b, establishing a target posterior probability density function composed of prior distribution, a likelihood function and a physical regularization item; c, in Markov chain Monte Carlo iteration, sampling is carried out on the target posterior probability density function; and d, after iteration is finished, taking a sample set output by the Markov chain as crack state parameter posteriori distribution of the time step k, and extracting crack three-dimensional geometrical morphology and distribution information from the sample set.

Description

Technical Field

[0001] The present invention belongs to the field of real-time monitoring and relates to a method and system for real-time monitoring of fracturing cracks based on time-varying electric field dynamic data. Background Art Hydraulic fracturing is key to developing unconventional energy sources such as shale oil and gas. It creates a network of artificial fractures with high conductivity within underground reservoirs. The geometry of the fractures, including their length, height, aperture, and spatial orientation, determines the effectiveness of reservoir stimulation and the ultimate productivity of oil and gas wells. Real-time monitoring of fracture expansion during fracturing operations can optimize construction parameters, control fracture growth, predict productivity, and assess environmental risks. Microseismic monitoring is the most widely used fracture monitoring technology. It uses seismic waves generated by rock fractures to delineate the fracture's impact area. However, microseismic events only represent the shear failure point of the rock and cannot directly reflect the main fracture morphology of effective proppant placement. Electromagnetic fracture monitoring utilizes the electrical anomalies formed by injecting conductive fracturing fluid into high-resistivity reservoir rock. Electrode arrays deployed at the surface or in wells monitor changes in the electric or magnetic field to invert the fracture geometry. Compared to microseismic technology, electromagnetic methods can more directly track the fracture body filled with conductive fluid. However, electromagnetic inversion is a typical nonlinear, highly ill-conditioned geophysical problem with severe multi-solution problems. Traditional inversion methods often lack effective physical constraints, resulting in inversion results that may not conform to the principles of rock mechanics. In addition, the three-dimensional electromagnetic forward model is computationally intensive. When combined with inversion algorithms such as Monte Carlo that require tens of thousands or even millions of iterations, real-time monitoring becomes almost impossible. Although the Bayesian inversion method based on Markov Chain Monte Carlo (MCMC) can assess uncertainty, its standard sampling strategy is inefficient and converges slowly, making it difficult to adapt to the real-time requirements of dynamic crack expansion during fracturing. Summary of the Invention

[0002] To address the above problems, the present invention proposes a real-time monitoring method for hydraulic fractures based on time-varying electric field dynamic data, comprising the following steps: Obtaining a time-varying electric field data sequence collected by electric field sensors deployed in the fracturing area; For each monitoring time step k in the electric field data series, the following sequential Bayesian inversion is performed: a. Based on the sample of the posterior distribution of the fracture state parameters at time step k-1, kernel density estimation is used to construct the prior distribution of the state parameters at time step k; the state parameters include the length, aperture, azimuth and inclination of the fracture; b. Determine a physical regularization term based on physical constraints calculated from the regional geostress field and the fracture tip stress intensity factor, and determine a likelihood function based on the degree of match between the forward modeled electric field data and the measured electric field data under given state parameters; and establish a target posterior probability density function consisting of the prior distribution, the likelihood function, and the physical regularization term. c. sampling the target posterior probability density function in a Markov chain Monte Carlo iteration; d. After the iteration, the sample set output by the Markov chain is used as the posterior distribution of the crack state parameters at time step k, and the three-dimensional geometric morphology and distribution information of the crack is extracted from it.

[0003] Optionally, sampling the target posterior probability density function is specifically: c1. Based on the crack propagation rate estimated at the previous time step k-1, adjust the covariance related to crack propagation in the multivariate normal proposal distribution and generate the first-stage candidate samples from it; c2. Using a multi-fidelity strategy to evaluate the likelihood function of the candidate samples in the first stage: using the low-fidelity model to calculate the initial likelihood value, when the preset switching criteria are met, the high-fidelity model is started to calculate the exact likelihood value, and the results of the two models are combined for correction; c3. When the first-stage candidate sample is rejected, the target posterior probability density gradient at the rejected sample is calculated, and along the gradient direction, the shrunken proposal distribution is used to generate the second-stage candidate sample.

[0004] Optionally, the kernel density estimation is used to construct the prior distribution of the state parameters at time step k, specifically: The kernel density estimation method is used to construct the prior distribution of the crack state parameters at time step k-1, which can reflect the uncertainty of the crack evolution over time.

[0005] Optionally, the physical regularization term is obtained based on the physical constraints calculated based on the regional in-situ stress field and the stress intensity factor at the crack tip, specifically: The fracture tip stress intensity factor of the candidate fracture geometry is compared with a preset rock fracture toughness, and when the stress intensity factor exceeds the rock fracture toughness, a preset low probability penalty is imposed on the candidate fracture geometry.

[0006] Optionally, the likelihood function is obtained according to the matching degree between the forward modeled electric field data and the measured electric field data under given state parameters, specifically: For any given state parameter sample, the theoretical electric field data is calculated through the forward model; Calculating the residuals between the theoretical electric field data and the measured electric field data at all sensor positions; An exponential term of a Gaussian likelihood function is constructed based on a weighted quadratic norm of the residual, and the Gaussian likelihood function is used as the likelihood function; wherein the weight is determined by the covariance of the measurement noise.

[0007] Optionally, the step c1 is specifically: The crack propagation rate is estimated based on the posterior distribution of time steps k-1 and k-2, and the components corresponding to the crack length and aperture in the covariance matrix of the multivariate normal proposal distribution are increased according to the propagation rate.

[0008] Optionally, the step c2 is specifically: A computationally inexpensive low-fidelity forward model is used to perform a preliminary likelihood assessment on the first-stage candidate samples; When the preliminary likelihood evaluation result meets a preset switching criterion, a high-fidelity forward model with high computational cost and high accuracy is started to calculate the precise likelihood value, and the precise likelihood value is corrected in combination with the calculation results of the low-fidelity and high-fidelity models.

[0009] Optionally, the step c3 is specifically: When a candidate sample in the first stage is rejected, the gradient of the target posterior probability density of the rejected sample position is calculated; and using the gradient information, a second-stage candidate sample pointing to a higher probability area is generated from a variance-shrinking proposal distribution.

[0010] In addition, the present invention also provides a real-time monitoring system for hydraulic fractures based on time-varying electric field dynamic data, comprising the following modules: A data acquisition module is used to obtain a time-varying electric field data sequence collected by electric field sensors deployed in the fracturing area; The inversion monitoring module is configured to perform the following sequential Bayesian inversion for each monitoring time step k in the electric field data sequence: a. Based on the sample of the posterior distribution of the fracture state parameters at time step k-1, kernel density estimation is used to construct the prior distribution of the state parameters at time step k; the state parameters include the length, aperture, azimuth and inclination of the fracture; b. Determine a physical regularization term based on physical constraints calculated from the regional geostress field and the fracture tip stress intensity factor, and determine a likelihood function based on the degree of match between the forward modeled electric field data and the measured electric field data under given state parameters; and establish a target posterior probability density function consisting of the prior distribution, the likelihood function, and the physical regularization term. c. sampling the target posterior probability density function in a Markov chain Monte Carlo iteration; d. After the iteration, the sample set output by the Markov chain is used as the posterior distribution of the crack state parameters at time step k, and the three-dimensional geometric morphology and distribution information of the crack is extracted from it.

[0011] Optionally, sampling the target posterior probability density function is specifically: c1. Based on the crack propagation rate estimated at the previous time step k-1, adjust the covariance related to crack propagation in the multivariate normal proposal distribution and generate the first-stage candidate samples from it; c2. Using a multi-fidelity strategy to evaluate the likelihood function of the candidate samples in the first stage: using the low-fidelity model to calculate the initial likelihood value, when the preset switching criteria are met, the high-fidelity model is started to calculate the exact likelihood value, and the results of the two models are combined for correction; c3. When the first-stage candidate sample is rejected, the target posterior probability density gradient at the rejected sample is calculated, and along the gradient direction, the shrunken proposal distribution is used to generate the second-stage candidate sample.

[0012] Optionally, the kernel density estimation is used to construct the prior distribution of the state parameters at time step k, specifically: The kernel density estimation method is used to construct the prior distribution of the crack state parameters at time step k-1, which can reflect the uncertainty of the crack evolution over time.

[0013] Optionally, the physical regularization term is obtained based on the physical constraints calculated based on the regional in-situ stress field and the stress intensity factor at the crack tip, specifically: The fracture tip stress intensity factor of the candidate fracture geometry is compared with a preset rock fracture toughness, and when the stress intensity factor exceeds the rock fracture toughness, a preset low probability penalty is imposed on the candidate fracture geometry.

[0014] Optionally, the likelihood function is obtained according to the matching degree between the forward modeled electric field data and the measured electric field data under given state parameters, specifically: For any given state parameter sample, the theoretical electric field data is calculated through the forward model; Calculating the residuals between the theoretical electric field data and the measured electric field data at all sensor positions; An exponential term of a Gaussian likelihood function is constructed based on a weighted quadratic norm of the residual, and the Gaussian likelihood function is used as the likelihood function; wherein the weight is determined by the covariance of the measurement noise.

[0015] Optionally, the step c1 is specifically: The crack propagation rate is estimated based on the posterior distribution of time steps k-1 and k-2, and the components corresponding to the crack length and aperture in the covariance matrix of the multivariate normal proposal distribution are increased according to the propagation rate.

[0016] Optionally, the step c2 is specifically: A computationally inexpensive low-fidelity forward model is used to perform a preliminary likelihood assessment on the first-stage candidate samples; When the preliminary likelihood evaluation result meets a preset switching criterion, a high-fidelity forward model with high computational cost and high accuracy is started to calculate the precise likelihood value, and the precise likelihood value is corrected in combination with the calculation results of the low-fidelity and high-fidelity models.

[0017] Optionally, the step c3 is specifically: When a candidate sample in the first stage is rejected, the gradient of the target posterior probability density of the rejected sample position is calculated; and using the gradient information, a second-stage candidate sample pointing to a higher probability area is generated from a variance-shrinking proposal distribution.

[0018] Compared with the existing technology, the present invention uses a sequential Bayesian inversion framework to use the inversion results of the previous time step as the prior information of the current time step, ensuring the temporal continuity of crack extension. At the same time, physical constraints based on the ground stress field and the stress intensity factor at the crack tip are integrated into the objective function to make the inversion results conform to the laws of rock mechanics. The crack extension information at the previous moment is used to guide the generation of proposed samples. A multi-fidelity model strategy is adopted to evaluate the likelihood function by combining a high-efficiency low-precision model with a high-precision model. A gradient-based remedial sampling mechanism is introduced to improve sampling efficiency and convergence speed, reducing computational cost while ensuring accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] Figure 1 Schematic diagram of the k-1 moment posterior sample and k moment prior distribution; Figure 2 Schematic diagram of the probability penalty function based on the physical regularization term; Figure 3 A comparison chart of the standard proposal distribution and the adaptive proposal distribution; Figure 4 Schematic diagram of gradient-based delayed rejection sampling. DETAILED DESCRIPTION

[0020] To make the purpose, technical solutions, and advantages of the embodiments of this application more clear, the technical solutions in the embodiments of this application will be clearly and completely described below in conjunction with the drawings in the embodiments of this application. Obviously, the described embodiments are part of the embodiments of this application, not all of the embodiments. Based on the embodiments in this application, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of this application.

[0021] In this application, the term "plurality" refers to two or more. In addition, it should be understood that in the description of this application, the terms "first" and "second" are used only for the purpose of distinguishing descriptions and should not be understood as indicating or implying relative importance or order.

[0022] In a first embodiment of the present invention, a method for real-time monitoring of hydraulic fractures based on time-varying electric field dynamic data is proposed, comprising the following steps: Obtaining a time-varying electric field data sequence collected by electric field sensors deployed in the fracturing area; For example, dipole-dipole electrode pairs are deployed as sensors on the surface around the fracturing well or in adjacent wells. A data acquisition system records the time-varying potential difference data during fracturing fluid injection, forming an electric field data sequence containing multiple time steps. For example, data is recorded once every minute, forming an observation data set with time steps k = 1, 2, 3, etc.

[0023] For each monitoring time step k in the electric field data series, the following sequential Bayesian inversion is performed: a. Based on the sample of the posterior distribution of the fracture state parameters at time step k-1, kernel density estimation is used to construct the prior distribution of the state parameters at time step k; the state parameters include the length, aperture, azimuth and inclination of the fracture; Using the posterior sample set obtained at the previous time step k-1, a continuous, nonparametric probability density function is fitted using the Gaussian kernel density estimation method. This probability density function serves as the joint prior distribution of the four state parameters of the crack length, aperture, azimuth, and inclination at the current time step k. For the initial time step k=1, an uninformative prior, such as a uniform distribution, is used.

[0024] b. Determine a physical regularization term based on physical constraints calculated from the regional geostress field and the fracture tip stress intensity factor, and determine a likelihood function based on the degree of match between the forward modeled electric field data and the measured electric field data under given state parameters; and establish a target posterior probability density function consisting of the prior distribution, the likelihood function, and the physical regularization term. The regional geostress field refers to the naturally occurring stress state in the rock mass within a specific geological region, unaffected by engineering activities such as drilling and fracturing. It is usually described by the magnitude and direction of three mutually perpendicular principal stresses (maximum horizontal principal stress, minimum horizontal principal stress, and vertical stress). The geostress field is the most critical factor controlling the initiation and propagation direction of hydraulic fracturing cracks. Cracks generally propagate in the direction of the maximum principal stress and open perpendicularly to the direction of the minimum principal stress. The crack tip stress intensity factor is a parameter in fracture mechanics that quantifies the concentration of the stress field near the crack tip. In one embodiment, the scalar posterior probability density function is proportional to the product of the prior distribution, the likelihood function, and the physical regularization term. The likelihood function adopts a Gaussian form, and its exponential term is the L2 norm residual between the electric field data calculated using the three-dimensional finite element forward model and the measured electric field data. The physical regularization term is implemented through a penalty function. When the stress intensity factor at the crack tip is less than the fracture toughness of the formation, or the crack propagation direction deviates too much from the direction of the maximum principal stress, the function assigns an extremely low probability value, thereby penalizing the fracture morphology that does not conform to the principles of rock mechanics in the inversion.

[0025] c. In the Markov chain Monte Carlo iteration, the target posterior probability density function is sampled: c1. Based on the crack propagation rate estimated at the previous time step k-1, adjust the covariance related to crack propagation in the multivariate normal proposal distribution and generate the first-stage candidate samples from it; The average growth rates of crack length and aperture are calculated from the posterior samples at time step k-1. This growth rate information is used to update the covariance matrix of a four-dimensional multivariate normal proposal distribution. This increases the variances corresponding to the length and aperture dimensions, allowing it to explore a larger parameter space and reflect the growth trend of the crack. The variances corresponding to other parameters, such as azimuth and inclination, are kept small. Samples are then drawn from this adjusted proposal distribution to generate candidate state parameter samples.

[0026] c2. Using a multi-fidelity strategy to evaluate the likelihood function of the candidate samples in the first stage: using the low-fidelity model to calculate the initial likelihood value, when the preset switching criteria are met, the high-fidelity model is started to calculate the exact likelihood value, and the results of the two models are combined for correction; The low-fidelity model preferably adopts a three-dimensional finite element electrical method forward model with a coarse mesh, which has a fast calculation speed but low accuracy. The high-fidelity model preferably adopts a three-dimensional finite element electrical method forward model with a fine mesh, which has a slow calculation speed but high accuracy. The preset switching criterion is that the likelihood value calculated by the low-fidelity model is greater than a preset threshold, for example, greater than 0.1. When the criterion is met, the high-fidelity model calculation is started, and the Gaussian process regression model is used to fit the difference between the low-fidelity and high-fidelity model calculation results. The corrected likelihood value is the difference between the low-fidelity likelihood value and the Gaussian process regression prediction.

[0027] c3. When the first-stage candidate sample is rejected, the target posterior probability density gradient at the rejected sample is calculated, and along the gradient direction, the shrunken proposal distribution is used to generate the second-stage candidate sample.

[0028] After rejecting a first-stage candidate sample based on the Metropolis-Hastings criterion, the adjoint state method is used to calculate the gradient of the target posterior probability density function with respect to the state parameter at the rejected sample point. This gradient direction is then used as the mean, and a shrinking multivariate normal distribution with a smaller covariance than the first-stage proposal distribution is used as the new proposal distribution to generate a remedial second-stage candidate sample.

[0029] d. After the iteration, the sample set output by the Markov chain is used as the posterior distribution of the crack state parameters at time step k, and the three-dimensional geometric morphology and distribution information of the crack is extracted from it.

[0030] After reaching a preset number of iterations or when the Markov chain converges, all accepted samples are combined to form a posterior probability distribution for the four parameters of crack length, aperture, azimuth, and inclination at time step k. The best estimate of each parameter is obtained by calculating the mean or mode of this sample set, and the variance or confidence interval is calculated to quantify the uncertainty. Finally, based on these best estimates of parameters, 3D visualization software is used to plot the 3D geometry and distribution of the cracks in the underground space.

[0031] More specifically, the kernel density estimation is used to construct the prior distribution of the state parameters at time step k, specifically: The kernel density estimation method is used to construct the prior distribution of the crack state parameters at time step k-1, which can reflect the uncertainty of the crack evolution over time.

[0032] After the Bayesian inversion for time step k-1 is complete, a posterior probability distribution is obtained. This distribution is typically represented by a large number of weighted samples or particles, for example, 10,000 particles. Each particle represents a possible fracture state, with state parameters such as length, aperture, and azimuth. Kernel density estimation utilizes these discrete particles to construct a continuous and smooth probability density function. This method places a kernel function, preferably a Gaussian kernel, at the location of each particle. All kernel functions are then superimposed to form an estimate of the probability density across the entire parameter space.

[0033] The constructed continuous distribution is used as the prior distribution of time step k. It inherits all the information of the previous moment, and its form is completely driven by data and does not depend on any preset distribution form. Therefore, it can well obtain the complex distribution characteristics such as multi-peak or asymmetric that may be produced by crack evolution, such as Figure 1As shown in Figure 2 . Furthermore, uncertainty is added by adjusting the kernel function's bandwidth parameter. For example, if a fracture is expected to grow rapidly, the kernel function bandwidth corresponding to the fracture length parameter can be increased. This makes the prior distribution flatter in the length dimension, allowing the inversion process to explore a wider range of fracture lengths at time step k, reflecting the uncertainty in the fracture's temporal evolution.

[0034] More specifically, the physical constraints calculated based on the regional in-situ stress field and the stress intensity factor at the crack tip yield a physical regularization term, specifically: The fracture tip stress intensity factor of the candidate fracture geometry is compared with a preset rock fracture toughness, and when the stress intensity factor exceeds the rock fracture toughness, a preset low probability penalty is imposed on the candidate fracture geometry.

[0035] This embodiment introduces physical constraints based on the principles of fracture mechanics. The fracture toughness of rock, denoted by KIC, is an inherent property of a material that measures its ability to resist crack propagation. For certain rock types, such as shale, the fracture toughness may be preset to 0.7 MPa times the square root of a meter. For each candidate fracture geometry generated during the inversion process, a mechanical forward model is used to calculate the stress intensity factor (KI) at the fracture tip under the current in-situ stress field.

[0036] According to linear elastic fracture mechanics, cracks will only grow unstably when the stress intensity factor (KI) reaches or exceeds the fracture toughness (KIC). Therefore, if the calculated results of a candidate crack morphology show that its KI value is much greater than the KIC value, for example, the calculated KI is 2.5 and the KIC is only 0.7, then this state is physically unstable and unrealistic because it should have expanded to another morphology long ago. Figure 2 This is a schematic diagram of a probability penalty function based on a physical regularization term. Those skilled in the art should know that the probability penalty function is not limited to Figure 2 To incorporate this information into the inversion, an extremely low probability penalty is imposed on such physically invalid candidate fractures, for example, by multiplying their posterior probability by a very small value, such as 10 to the power of -10. This eliminates such solutions during the sampling process, guiding the inversion results toward more physically reasonable fracture morphologies.

[0037] More specifically, the likelihood function is obtained based on the matching degree between the forward modeled electric field data and the measured electric field data under given state parameters, specifically: For any given state parameter sample, the theoretical electric field data is calculated through the forward model; Calculating the residuals between the theoretical electric field data and the measured electric field data at all sensor positions; An exponential term of a Gaussian likelihood function is constructed based on a weighted quadratic norm of the residual, and the Gaussian likelihood function is used as the likelihood function; wherein the weight is determined by the covariance of the measurement noise.

[0038] The likelihood function measures the degree of match between a given candidate fracture model and actual observed data. In this embodiment, for a candidate fracture state, such as a 10-meter-long, 5-mm-wide fracture, an electromagnetic field forward model, such as a finite element model, is first used to calculate the theoretical electric field or voltage response that the fracture would generate at all sensor locations, such as sixteen electrode pairs arranged on the surface, if it were present underground. This data vector containing sixteen theoretical voltage values ​​is then compared one by one with the sixteen voltage values ​​obtained from actual field measurements, and the difference between them, the residual vector, is calculated. The Gaussian likelihood function is constructed based on the assumption that the observation error follows a Gaussian distribution. Its core component is an exponential term, which is the weighted quadratic norm of the residual. The weight is determined by the inverse of the covariance matrix of the measurement noise. For example, if the standard deviation of the measurement noise of a sensor is 0.1 millivolt, its data is highly reliable and its weight in the calculation is high; conversely, if the standard deviation of the noise of another sensor is 5 millivolts, its weight is low.

[0039] More specifically, the step c1 is as follows: The crack propagation rate is estimated based on the posterior distribution of time steps k-1 and k-2, and the components corresponding to the crack length and aperture in the covariance matrix of the multivariate normal proposal distribution are increased according to the propagation rate.

[0040] To improve the efficiency of Markov Chain Monte Carlo sampling, especially when dealing with dynamic evolution problems, this embodiment employs an adaptive proposal distribution strategy. This proposal distribution is used to generate new candidate samples in the parameter space. Historical information is used to predict the dynamic behavior of the crack. Specifically, the rate of change of the key crack parameters is estimated by comparing the mean or mode of the posterior distribution at time steps k-2 and k-1. For example, if the average crack length at time step k-2 is 3.2 meters and increases to 3.8 meters at time step k-1, the crack growth rate can be estimated to be approximately 0.6 per time step.

[0041] The estimated propagation rate is used to adjust the proposed distribution for time step k. The proposed distribution is usually a multivariate normal distribution, whose shape is controlled by the covariance matrix. Based on the estimated propagation rate of 0.6 meters, the variance term corresponding to the crack length parameter in the covariance matrix will be increased accordingly, for example, the value of the diagonal element will be increased from 0.05 to 0.3, as shown in the following example: Figure 3As shown in Figure 2, this will cause the proposed distribution to jump with larger steps in the crack length dimension, allowing for more efficient exploration following the actual crack extension trend. The adaptive mechanism enables the sampling process to focus computing resources on the regions in the parameter space most likely to change, thereby accelerating convergence.

[0042] More specifically, the step c2 is as follows: A computationally inexpensive low-fidelity forward model is used to perform a preliminary likelihood assessment on the first-stage candidate samples; When the preliminary likelihood evaluation result meets a preset switching criterion, a high-fidelity forward model with high computational cost and high accuracy is started to calculate the precise likelihood value, and the precise likelihood value is corrected in combination with the calculation results of the low-fidelity and high-fidelity models.

[0043] To address the issue of excessive computational time associated with high-precision forward models, two forward models with different accuracies are used simultaneously. A low-fidelity model, such as a finite element model using a coarse mesh or a simplified analytical solution, is very fast, perhaps taking only 0.5 seconds to compute. A high-fidelity model, a finite element model using a fine mesh, is more accurate but computationally expensive, potentially taking up to 5 minutes. For each candidate sample generated by the Monte Carlo sampler, the fast low-fidelity model is invoked to calculate a preliminary likelihood. A switching criterion, such as a likelihood threshold, is set. If the preliminary likelihood calculated by the low-fidelity model falls below the threshold, indicating that the candidate sample is extremely mismatched with the observed data, the sample is rejected without further expensive computation. For example, if 90% of candidate samples are quickly filtered out by the low-fidelity model, significant computational time can be saved. Only when the preliminary likelihood value exceeds the threshold, indicating that the sample has a certain degree of plausibility, is the high-fidelity model activated for precise computation. To further improve accuracy, a small number of samples with both high- and low-fidelity model results will be used to establish a correction model to correct the calculation results of the high-fidelity model, obtaining a likelihood assessment that is both fast and accurate, thereby increasing the overall computing efficiency by dozens of times while ensuring accuracy.

[0044] More specifically, step c3 is as follows: When a candidate sample in the first stage is rejected, the gradient of the target posterior probability density of the rejected sample position is calculated; and using the gradient information, a second-stage candidate sample pointing to a higher probability area is generated from a variance-shrinking proposal distribution.

[0045] An improved sampling strategy called delayed rejection is adopted to make full use of the information contained in the rejected samples. In the standard Monte Carlo method, a rejected proposal step is usually regarded as a wasted calculation. In this embodiment, when a candidate sample in the first stage is rejected, it is not abandoned immediately. Instead, the gradient of the target posterior probability density at the point of the candidate sample in the first stage is calculated. The gradient is a vector that points to the direction of the fastest growth of the probability density in the parameter space, such as Figure 4 As shown. Gradient information is used to guide the second proposal. A second-stage candidate sample is generated from a new proposal distribution. The center of the new proposal distribution is no longer the position of the current chain, but a small step away from the current position along the gradient direction just calculated. For example, if the gradient shows that increasing the crack opening can improve the posterior probability the fastest, then the second proposed sample will have a larger opening value than the first proposal. At the same time, the covariance matrix of the second proposal distribution will be shrunk, for example, its variance is one-fourth of the variance of the first proposal, indicating a smaller local search. The second proposed sample has a greater chance of being accepted because it is intentionally pushed to an area with higher probability, thereby converting a failed attempt into a valid move and improving the overall efficiency of the sampler.

[0046] In a first embodiment of the present invention, a real-time monitoring system for hydraulic fractures based on time-varying electric field dynamic data is provided, comprising the following modules: A data acquisition module is used to obtain a time-varying electric field data sequence collected by electric field sensors deployed in the fracturing area; The inversion monitoring module is configured to perform the following sequential Bayesian inversion for each monitoring time step k in the electric field data sequence: a. Based on the sample of the posterior distribution of the fracture state parameters at time step k-1, kernel density estimation is used to construct the prior distribution of the state parameters at time step k; the state parameters include the length, aperture, azimuth and inclination of the fracture; b. Determine a physical regularization term based on physical constraints calculated from the regional geostress field and the fracture tip stress intensity factor, and determine a likelihood function based on the degree of match between the forward modeled electric field data and the measured electric field data under given state parameters; and establish a target posterior probability density function consisting of the prior distribution, the likelihood function, and the physical regularization term. c. sampling the target posterior probability density function in a Markov chain Monte Carlo iteration; d. After the iteration, the sample set output by the Markov chain is used as the posterior distribution of the crack state parameters at time step k, and the three-dimensional geometric morphology and distribution information of the crack is extracted from it.

[0047] More specifically, the sampling of the target posterior probability density function is specifically as follows: c1. Based on the crack propagation rate estimated at the previous time step k-1, adjust the covariance related to crack propagation in the multivariate normal proposal distribution and generate the first-stage candidate samples from it; c2. Using a multi-fidelity strategy to evaluate the likelihood function of the candidate samples in the first stage: using the low-fidelity model to calculate the initial likelihood value, when the preset switching criteria are met, the high-fidelity model is started to calculate the exact likelihood value, and the results of the two models are combined for correction; c3. When the first-stage candidate sample is rejected, the target posterior probability density gradient at the rejected sample is calculated, and along the gradient direction, the shrunken proposal distribution is used to generate the second-stage candidate sample.

[0048] More specifically, the kernel density estimation is used to construct the prior distribution of the state parameters at time step k, specifically: The kernel density estimation method is used to construct the prior distribution of the crack state parameters at time step k-1, which can reflect the uncertainty of the crack evolution over time.

[0049] More specifically, the physical constraints calculated based on the regional in-situ stress field and the stress intensity factor at the crack tip yield a physical regularization term, specifically: The fracture tip stress intensity factor of the candidate fracture geometry is compared with a preset rock fracture toughness, and when the stress intensity factor exceeds the rock fracture toughness, a preset low probability penalty is imposed on the candidate fracture geometry.

[0050] More specifically, the likelihood function is obtained based on the matching degree between the forward modeled electric field data and the measured electric field data under given state parameters, specifically: For any given state parameter sample, the theoretical electric field data is calculated through the forward model; Calculating the residuals between the theoretical electric field data and the measured electric field data at all sensor positions; An exponential term of a Gaussian likelihood function is constructed based on a weighted quadratic norm of the residual, and the Gaussian likelihood function is used as the likelihood function; wherein the weight is determined by the covariance of the measurement noise.

[0051] More specifically, the step c1 is as follows: The crack propagation rate is estimated based on the posterior distribution of time steps k-1 and k-2, and the components corresponding to the crack length and aperture in the covariance matrix of the multivariate normal proposal distribution are increased according to the propagation rate.

[0052] More specifically, the step c2 is as follows: A computationally inexpensive low-fidelity forward model is used to perform a preliminary likelihood assessment on the first-stage candidate samples; When the preliminary likelihood evaluation result meets a preset switching criterion, a high-fidelity forward model with high computational cost and high accuracy is started to calculate the precise likelihood value, and the precise likelihood value is corrected in combination with the calculation results of the low-fidelity and high-fidelity models.

[0053] More specifically, step c3 is as follows: When a candidate sample in the first stage is rejected, the gradient of the target posterior probability density of the rejected sample position is calculated; and using the gradient information, a second-stage candidate sample pointing to a higher probability area is generated from a variance-shrinking proposal distribution.

[0054] Through the description of the above embodiments, it can be seen that those skilled in the art can clearly understand that the present application can be implemented by means of software plus a necessary general hardware platform. Based on this understanding, the technical solution of the present application, or the part that contributes to the prior art, can be embodied in the form of a software product, which can be stored in a storage medium such as ROM / RAM, a magnetic disk, an optical disk, etc., and includes a number of instructions for enabling a computer device (which can be a personal computer, a server, or a network device, etc.) to execute the methods described in various embodiments of the present application or certain parts of the embodiments.

[0055] Each embodiment in this specification is described in a progressive manner. The same or similar parts between the embodiments can be referred to each other. Each embodiment focuses on the differences from other embodiments. In particular, for system or system embodiments, since they are basically similar to method embodiments, the description is relatively simple. For relevant parts, refer to the partial description of the method embodiment. The system and system embodiments described above are merely schematic, wherein the units described as separate components may or may not be physically separated, and the components displayed as units may or may not be physical units, that is, they may be located in one place, or they may be distributed on multiple network units. Some or all of the modules can be selected according to actual needs to achieve the purpose of the solution of this embodiment. A person of ordinary skill in the art can understand and implement it without expending creative work.

[0056] The above describes in detail the method and electronic device for providing commodity object information provided by this application. Specific examples are used herein to illustrate the principles and implementation methods of this application. The description of the above embodiments is only intended to help understand the method and core concept of this application. At the same time, for those skilled in the art, based on the concept of this application, there may be changes in the specific implementation methods and application scope. In summary, the contents of this specification should not be construed as limiting this application.

Claims

1. A real-time monitoring method for hydraulic fractures based on time-varying electric field dynamic data, characterized in that: The following steps are involved: Obtaining a time-varying electric field data sequence collected by electric field sensors deployed in the fracturing area; For each monitoring time step k in the electric field data series, the following sequential Bayesian inversion is performed: a. Based on the sample of the posterior distribution of the fracture state parameters at time step k-1, kernel density estimation is used to construct the prior distribution of the state parameters at time step k; the state parameters include the length, aperture, azimuth and inclination of the fracture; b. Determine a physical regularization term based on physical constraints calculated from the regional geostress field and the fracture tip stress intensity factor, and determine a likelihood function based on the degree of match between the forward modeled electric field data and the measured electric field data under given state parameters; and establish a target posterior probability density function consisting of the prior distribution, the likelihood function, and the physical regularization term. c. sampling the target posterior probability density function in a Markov chain Monte Carlo iteration; d. After the iteration, the sample set output by the Markov chain is used as the posterior distribution of the crack state parameters at time step k, and the three-dimensional geometric morphology and distribution information of the crack is extracted from it.

2. The method according to claim 1, characterized in that The sampling of the target posterior probability density function is specifically as follows: c1. Based on the crack propagation rate estimated at the previous time step k-1, adjust the covariance related to crack propagation in the multivariate normal proposal distribution and generate the first-stage candidate samples from it; c2. Using a multi-fidelity strategy to evaluate the likelihood function of the candidate samples in the first stage: using the low-fidelity model to calculate the initial likelihood value, when the preset switching criteria are met, the high-fidelity model is started to calculate the exact likelihood value, and the results of the two models are combined for correction; c3. When the first-stage candidate sample is rejected, the target posterior probability density gradient at the rejected sample is calculated, and along the gradient direction, the shrunken proposal distribution is used to generate the second-stage candidate sample.

3. The method according to claim 1, characterized in that The kernel density estimation is used to construct the prior distribution of the state parameters at time step k, specifically: The kernel density estimation method is used to construct the prior distribution of the crack state parameters at time step k-1, which can reflect the uncertainty of the crack evolution over time.

4. The method according to claim 1, wherein The physical regularization term is obtained based on the physical constraints calculated based on the regional in-situ stress field and the stress intensity factor at the crack tip, specifically: The fracture tip stress intensity factor of the candidate fracture geometry is compared with a preset rock fracture toughness, and when the stress intensity factor exceeds the rock fracture toughness, a preset low probability penalty is imposed on the candidate fracture geometry.

5. The method according to claim 1, wherein The likelihood function is obtained based on the matching degree between the forward modeled electric field data and the measured electric field data under given state parameters, specifically: For any given state parameter sample, the theoretical electric field data is calculated through the forward model; Calculating the residuals between the theoretical electric field data and the measured electric field data at all sensor positions; An exponential term of a Gaussian likelihood function is constructed based on a weighted quadratic norm of the residual, and the Gaussian likelihood function is used as the likelihood function; wherein the weight is determined by the covariance of the measurement noise.

6. The method according to claim 2, characterized in that The step c1 is specifically as follows: The crack propagation rate is estimated based on the posterior distribution of time steps k-1 and k-2, and the components corresponding to the crack length and aperture in the covariance matrix of the multivariate normal proposal distribution are increased according to the propagation rate.

7. The method according to claim 2, characterized in that The step c2 is specifically as follows: A computationally inexpensive low-fidelity forward model is used to perform a preliminary likelihood assessment on the first-stage candidate samples; When the preliminary likelihood evaluation result meets a preset switching criterion, a high-fidelity forward model with high computational cost and high accuracy is started to calculate the precise likelihood value, and the precise likelihood value is corrected in combination with the calculation results of the low-fidelity and high-fidelity models.

8. The method according to claim 2, characterized in that The step c3 is specifically as follows: When a candidate sample in the first stage is rejected, the gradient of the target posterior probability density of the rejected sample position is calculated; and using the gradient information, a second-stage candidate sample pointing to a higher probability area is generated from a variance-shrinking proposal distribution.

9. A real-time monitoring system for fracturing cracks based on time-varying electric field dynamic data, characterized in that: Includes the following modules: A data acquisition module is used to obtain a time-varying electric field data sequence collected by electric field sensors deployed in the fracturing area; The inversion monitoring module is configured to perform the following sequential Bayesian inversion for each monitoring time step k in the electric field data sequence: a. Based on the sample of the posterior distribution of the fracture state parameters at time step k-1, kernel density estimation is used to construct the prior distribution of the state parameters at time step k; the state parameters include the length, aperture, azimuth and inclination of the fracture; b. Determine a physical regularization term based on physical constraints calculated from the regional geostress field and the fracture tip stress intensity factor, and determine a likelihood function based on the degree of match between the forward modeled electric field data and the measured electric field data under given state parameters; and establish a target posterior probability density function consisting of the prior distribution, the likelihood function, and the physical regularization term. c. sampling the target posterior probability density function in a Markov chain Monte Carlo iteration; d. After the iteration, the sample set output by the Markov chain is used as the posterior distribution of the crack state parameters at time step k, and the three-dimensional geometric morphology and distribution information of the crack is extracted from it.

10. The system according to claim 9, characterized in that The sampling of the target posterior probability density function is specifically as follows: c1. Based on the crack propagation rate estimated at the previous time step k-1, adjust the covariance related to crack propagation in the multivariate normal proposal distribution and generate the first-stage candidate samples from it; c2. Using a multi-fidelity strategy to evaluate the likelihood function of the candidate samples in the first stage: using the low-fidelity model to calculate the initial likelihood value, when the preset switching criteria are met, the high-fidelity model is started to calculate the exact likelihood value, and the results of the two models are combined for correction; c3. When the first-stage candidate sample is rejected, the target posterior probability density gradient at the rejected sample is calculated, and along the gradient direction, the shrunken proposal distribution is used to generate the second-stage candidate sample.

Images