Real-time monitoring method and system for fracturing cracks based on time-varying electric field dynamic data

Abstract

The application provides a fracturing fracture real-time monitoring method and system based on time-varying electric field dynamic data, and specifically, time-varying electric field data sequences collected by electric field sensors arranged in a fracturing working area are acquired; for each monitoring time step k in the electric field data sequences, the following sequential Bayesian inversion is performed: a. based on samples of a fracture state parameter posterior distribution of a time step k-1, a state parameter prior distribution of a time step k is constructed by using kernel density estimation; b. a target posterior probability density function composed of a prior distribution, a likelihood function and a physical regularization term is established; c. in Markov chain Monte Carlo iteration, the target posterior probability density function is sampled; d. after iteration, a sample set output by the Markov chain is taken as a fracture state parameter posterior distribution of the time step k, and fracture three-dimensional geometric morphology and distribution information is extracted 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

[0002] 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 transformation and the ultimate productivity of the oil and gas well. Real-time monitoring of fracture expansion during fracturing operations can optimize construction parameters, control fracture growth, predict production capacity, and assess environmental risks. Microseismic monitoring, the most widely used fracture monitoring technology, uses seismic waves generated by rock fracturing to delineate the fracture's impact area.

[0003] However, microseismic events only represent shear failure points in the rock and cannot directly reflect the main fracture morphology of effective proppant placement. Electromagnetic fracture monitoring technology utilizes electrical anomalies formed by injecting conductive fracturing fluid into high-resistivity reservoir rock. Electrode arrays deployed on the surface or in wells monitor changes in the electric or magnetic field to invert fracture geometry. Compared to microseismic techniques, 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 potential. Traditional inversion methods often lack effective physical constraints, resulting in inversion results that may not conform to rock mechanics principles. Furthermore, the computational complexity of three-dimensional electromagnetic forward models is enormous. When combined with inversion algorithms such as Monte Carlo, which require tens of thousands or even millions of iterations, real-time monitoring becomes nearly impossible. Bayesian inversion methods based on Markov Chain Monte Carlo (MCMC) can assess uncertainty, but their standard sampling strategies are inefficient and converge slowly, making them difficult to adapt to the real-time requirements of dynamic fracture growth during fracturing. Summary of the Invention

[0004] 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:

[0005] Obtaining a time-varying electric field data sequence collected by electric field sensors deployed in the fracturing area;

[0006] For each monitoring time step k in the electric field data series, the following sequential Bayesian inversion is performed:

[0007] 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;

[0008] b. a physical regularization term is obtained based on the physical constraint of the regional geostress field and the stress intensity factor at the crack tip, and a likelihood function is obtained according to the matching degree of the forward electric field data under the given state parameter and the measured electric field data; a target posterior probability density function is established by the prior distribution, the likelihood function and the physical regularization term;

[0009] c. in the Markov chain Monte Carlo iteration, sampling is performed on the target posterior probability density function;

[0010] d. after the iteration is completed, a sample set output by the Markov chain is taken as a crack state parameter posterior distribution at a time step k, and crack three-dimensional geometric morphology and distribution information is extracted therefrom.

[0011] Optionally, the sampling on the target posterior probability density function is specifically:

[0012] c1. according to the crack propagation rate estimated at a previous time step k-1, adjusting the covariance related to crack propagation in the multivariate normal proposal distribution, and generating a first-stage candidate sample therefrom;

[0013] c2. a multi-fidelity strategy is used to evaluate the likelihood function of the first-stage candidate sample: an initial likelihood value is calculated by using a low-fidelity model, when a preset switching criterion is met, an accurate likelihood value is calculated by starting a high-fidelity model, and the results of the two models are combined for correction;

[0014] c3. when the first-stage candidate sample is rejected, the gradient of the target posterior probability density at the rejected sample is calculated, and a second-stage candidate sample is generated by using a shrunk proposal distribution along the gradient direction.

[0015] Optionally, the kernel density estimation is used to construct the state parameter prior distribution at the time step k, and specifically:

[0016] The kernel density estimation method is used to construct the prior distribution at the time step k which can reflect the uncertainty of the crack evolution over time according to the crack state parameter posterior distribution at the time step k-1.

[0017] Optionally, the physical regularization term is obtained based on the physical constraint of the regional geostress field and the stress intensity factor at the crack tip, and specifically:

[0018] The stress intensity factor at the crack tip of the candidate crack geometric morphology is compared with the preset rock fracture toughness, and when the stress intensity factor exceeds the rock fracture toughness, a preset low-probability penalty is applied to the candidate crack geometric morphology.

[0019] Optionally, the likelihood function is obtained according to the matching degree of the forward electric field data under the given state parameter and the measured electric field data, and specifically:

[0020] For any given state parameter sample, calculate theoretical electric field data by forward modeling;

[0021] Calculate the residual of the theoretical electric field data and the measured electric field data at all sensor positions;

[0022] Construct the exponential term of the Gaussian likelihood function based on the weighted quadratic norm of the residual, and take the Gaussian likelihood function as the likelihood function; wherein the weight is determined by the covariance of the measurement noise.

[0023] Optionally, the step c1 is specifically:

[0024] Estimate the propagation rate of the fracture based on the posterior distribution of time steps k-1 and k-2, and increase the components corresponding to the length and opening of the fracture in the covariance matrix of the multivariate normal proposal distribution according to the propagation rate.

[0025] Optionally, the step c2 is specifically:

[0026] A low-fidelity forward model with low computational cost is used to preliminarily evaluate the likelihood of the first-stage candidate sample;

[0027] 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 accurate likelihood value, and the accurate likelihood value is corrected in combination with the calculation results of the low-fidelity and high-fidelity models.

[0028] Optionally, the step c3 is specifically:

[0029] When the first-stage candidate sample is rejected, the gradient of the target posterior probability density of the rejected sample position is calculated; and a second-stage candidate sample pointing to a region with higher probability is generated from a variance-shrunk proposal distribution by using the gradient information.

[0030] In addition, the present application also provides a real-time monitoring system for fracturing fractures based on time-varying electric field dynamic data, comprising the following modules:

[0031] A data acquisition module is used to acquire the time-varying electric field data sequence collected by the electric field sensors arranged in the fracturing work area;

[0032] An inversion monitoring module is used to perform the following sequential Bayesian inversion for each monitoring time step k in the electric field data sequence:

[0033] a. Based on the sample of the fracture state parameter posterior distribution of time step k-1, a kernel density estimation is used to construct the prior distribution of the state parameter of time step k; the state parameter includes the length, opening, azimuth angle and inclination angle of the fracture;

[0034] b. a physical regularization term is obtained based on a physical constraint of a regional geostress field and a stress intensity factor at a crack tip, and a likelihood function is obtained according to matching degrees of forward electric field data and measured electric field data under given state parameters; a target posterior probability density function is established by the prior distribution, the likelihood function and the physical regularization term;

[0035] c. sampling is performed on the target posterior probability density function in Markov chain Monte Carlo iteration;

[0036] d. after iteration ends, a sample set output by the Markov chain is taken as a crack state parameter posterior distribution at a time step k, and crack three-dimensional geometric morphology and distribution information is extracted from the sample set.

[0037] Optionally, the sampling on the target posterior probability density function is specifically:

[0038] c1. according to a crack propagation rate estimated at a previous time step k-1, adjusting a covariance related to crack propagation in a multivariate normal proposal distribution, and generating a first-stage candidate sample from the multivariate normal proposal distribution;

[0039] c2. a multi-fidelity strategy is used to evaluate a likelihood function of the first-stage candidate sample: an initial likelihood value is calculated by using a low-fidelity model, when a preset switching criterion is met, an accurate likelihood value is calculated by using a high-fidelity model, and a result of the two models is combined for correction;

[0040] c3. when the first-stage candidate sample is rejected, a target posterior probability density gradient at the rejected sample is calculated, and a second-stage candidate sample is generated by using a shrunk proposal distribution in a direction of the gradient.

[0041] Optionally, the kernel density estimation is used to construct the state parameter prior distribution at the time step k, and specifically:

[0042] a kernel density estimation method is used to construct a prior distribution at the time step k, which can reflect uncertainty of crack evolution over time, according to a crack state parameter posterior distribution at a previous time step k-1.

[0043] Optionally, the physical regularization term is obtained based on a physical constraint of a regional geostress field and a stress intensity factor at a crack tip, and specifically:

[0044] a stress intensity factor at a crack tip of a candidate crack geometric morphology 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 applied to the candidate crack geometric morphology.

[0045] Optionally, the likelihood function is obtained according to matching degrees of the forward electric field data and the measured electric field data under a given state parameter, and specifically is:

[0046] For any given state parameter sample, theoretical electric field data are calculated through a forward model;

[0047] Residuals of the theoretical electric field data and the measured electric field data at all sensor positions are calculated;

[0048] An exponential term of a Gaussian likelihood function is constructed based on a weighted quadratic norm of the residuals, and the Gaussian likelihood function is taken as the likelihood function; wherein the weight is determined by a covariance of measurement noise.

[0049] Optionally, the step c1 is specifically:

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

[0051] Optionally, the step c2 is specifically:

[0052] A low-fidelity forward model with low calculation cost is used to preliminarily evaluate the likelihood of the first-stage candidate sample;

[0053] When the preliminary likelihood evaluation result meets a preset switching criterion, a high-fidelity forward model with high calculation cost and high accuracy is started to calculate an accurate likelihood value, and the accurate likelihood value is corrected in combination with calculation results of the low-fidelity and high-fidelity models.

[0054] Optionally, the step c3 is specifically:

[0055] When the first-stage candidate sample is rejected, the gradient of the target posterior probability density at the position of the rejected sample is calculated; and a second-stage candidate sample pointing to a region with higher probability is generated from a proposal distribution with variance shrinkage by using the gradient information.

[0056] Compared with the prior art, the sequential Bayesian inversion framework is used to ensure the continuity of the fracture propagation in time by taking the inversion result of the previous time step as the prior information of the current time step; meanwhile, the physical constraint based on the in-situ stress field and the stress intensity factor of the fracture tip is integrated into the target function, so that the inversion result conforms to the rock mechanics law; the fracture propagation information of the previous moment is used to guide the generation of the proposal sample, the multi-fidelity model strategy is used to efficiently evaluate the likelihood function by combining the low-precision model with the high-precision model, and the gradient-based remedial sampling mechanism is introduced, so that the sampling efficiency and the convergence speed are improved, and the calculation cost is reduced under the premise of ensuring the accuracy. BRIEF DESCRIPTION OF DRAWINGS

[0057] Figure 1 A schematic diagram of posterior samples at time k-1 and prior distribution at time k;

[0058] Figure 2 A schematic diagram of probability penalty function based on physical regularization term;

[0059] Figure 3 A comparison diagram of standard proposal distribution and adaptive proposal distribution;

[0060] Figure 4 A schematic diagram of gradient-based delayed rejection sampling. DETAILED DESCRIPTION

[0061] In order to make the purposes, technical solutions and advantages of the embodiments of the present application clearer, the technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only some of the embodiments of the present application, rather than all the embodiments of the present application. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative work belong to the scope of protection of the present application.

[0062] The plurality in the present application refers to two or more. In addition, it should be understood that in the description of the present application, the words "first", "second", etc. are only used for the purpose of distinguishing the description, and cannot be understood as indicating or implying relative importance, nor can it be understood as indicating or implying order.

[0063] In a first embodiment of the present application, a real-time monitoring method for fracturing fractures based on time-varying electric field dynamic data is proposed, comprising the following steps:

[0064] Obtaining the time-varying electric field data sequence collected by the electric field sensor arranged in the fracturing work area;

[0065] Arranging, for example, a dipole-dipole electrode pair as a sensor around the surface of the fracturing well or in the adjacent well, recording the time-varying potential difference data in the fracturing fluid injection process through a data acquisition system, forming an electric field data sequence containing multiple time steps, for example, recording data once every minute, constituting an observation data set of time step k=1, 2, 3...

[0066] For each monitoring time step k in the electric field data sequence, the following sequential Bayesian inversion is performed:

[0067] a. Based on the samples of the fracture state parameter posterior distribution at time step k-1, a kernel density estimation is used to construct the prior distribution of the state parameter at time step k; the state parameter includes the length, opening, azimuth angle and inclination angle of the fracture;

[0068] Using the posterior sample set obtained at the previous time step k-1, a continuous non-parametric probability density function is fitted by Gaussian kernel density estimation method, which serves as the joint prior distribution of the four state parameters of fracture length, aperture, azimuth and dip at the current time step k. For the initial time step k = 1, an uninformative prior is adopted, for example, a uniform distribution.

[0069] b. A physical regularization term is obtained based on the regional in-situ stress field and the physical constraint of fracture tip stress intensity factor, and a likelihood function is obtained according to the matching degree of the forward electric field data and the measured electric field data under the given state parameters; an objective posterior probability density function is established by the prior distribution, the likelihood function and the physical regularization term;

[0070] The regional in-situ stress field refers to the stress state naturally existing in the rock mass in a specific geological region without being affected by drilling, fracturing and other engineering activities. It is usually described by the magnitude and direction of three mutually perpendicular principal stresses (the maximum horizontal principal stress, the minimum horizontal principal stress and the vertical stress). The in-situ stress field is the most critical factor controlling the initiation and propagation direction of hydraulic fracture. The fracture generally propagates along the direction of the maximum principal stress and opens perpendicular to the direction of the minimum principal stress. The fracture tip stress intensity factor is a parameter in fracture mechanics, which quantifies the stress concentration near the tip of the fracture. In an embodiment, the 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 by the three-dimensional finite element forward model and the measured electric field data. The physical regularization term is realized by a penalty function, which gives a very low probability value when the stress intensity factor at the tip of the fracture is less than the fracture toughness of the formation, or the propagation direction of the fracture deviates too much from the direction of the maximum principal stress, so as to punish the fracture morphology that does not comply with the principles of rock mechanics in the inversion.

[0071] c. In the Markov chain Monte Carlo iteration, the objective posterior probability density function is sampled:

[0072] c1. According to the estimated fracture propagation rate at the previous time step k-1, the covariance related to fracture propagation in the multivariate normal proposal distribution is adjusted, and a first-stage candidate sample is generated therefrom;

[0073] The average growth rates of the fracture length and aperture are calculated from the posterior samples at the time step k-1, and these growth rate information is used to update the covariance matrix of a four-dimensional multivariate normal proposal distribution, so as to increase the variances corresponding to the length and aperture dimensions, so that it can explore a larger parameter space and reflect the growth trend of the fracture. The variances corresponding to the remaining parameters such as azimuth and dip are kept small. Sampling is performed from the adjusted proposal distribution to generate a candidate state parameter sample.

[0074] c2. Evaluate the likelihood function of the first-stage candidate sample using a multi-fidelity strategy: use a low-fidelity model to calculate an initial likelihood value, when a preset switching criterion is met, start a high-fidelity model to calculate an accurate likelihood value, and combine the results of the two models to correct them;

[0075] The low-fidelity model preferably uses a three-dimensional finite element electrical method forward model with coarse grid subdivision, which has fast calculation speed but low accuracy. The high-fidelity model preferably uses a three-dimensional finite element electrical method forward model with fine grid subdivision, which has 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 is started to calculate, and a Gaussian process regression model is used to fit the difference between the results of the low-fidelity and high-fidelity models, and the corrected likelihood value is the low-fidelity likelihood value plus the difference predicted by the Gaussian process regression.

[0076] c3. When the first-stage candidate sample is rejected, calculate the target posterior probability density gradient at the rejected sample, and generate a second-stage candidate sample using a shrunk proposal distribution along the gradient direction.

[0077] When the first-stage candidate sample is rejected according to the Metropolis-Hastings criterion, the gradient of the target posterior probability density function with respect to the state parameter at the rejected sample point is calculated using the adjoint state method. Then, taking the gradient direction as the mean, a shrunk 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.

[0078] d. After the iteration ends, the sample set output by the Markov chain is taken as the posterior distribution of the crack state parameter at time step k, and the three-dimensional geometric shape and distribution information of the crack is extracted from it.

[0079] After reaching the preset number of iterations or the convergence of the Markov chain, all accepted samples are combined to form the posterior probability distribution of the crack length, opening, azimuth angle and inclination angle at time step k. The best estimate value of each parameter is obtained by calculating the mean or mode of the sample set, and the variance or confidence interval is calculated to quantify the uncertainty. Finally, based on the best estimate value of these parameters, a three-dimensional visualization software is used to draw the three-dimensional geometric shape and distribution position of the crack in the underground space.

[0080] More specifically, the kernel density estimation is used to construct the prior distribution of the state parameter at time step k, specifically:

[0081] Using the kernel density estimation method, the prior distribution at time step k that can reflect the uncertainty of the crack evolution over time is constructed according to the posterior distribution of the crack state parameter at time step k-1.

[0082] 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.

[0083] 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 1 As 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.

[0084] 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:

[0085] 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.

[0086] 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.

[0087] 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 2For a probabilistic penalty function based on a physical regularization term, one skilled in the art will appreciate that the probabilistic penalty function is not limited to Figure 2 the probabilistic penalty function shown. To incorporate this information into the inversion, a very low probability penalty is imposed on such physically implausible candidates, e.g. multiplying their posterior probability by a very small number, say ten to the minus ten, so that such solutions are excluded in the sampling process, guiding the inversion result towards physically more plausible fracture shapes.

[0088] More specifically, the likelihood function is derived from the match between the forward modeled electric field data and the measured electric field data for a given state parameter, and is specifically:

[0089] For any given state parameter sample, the theoretical electric field data is calculated by a forward modeling;

[0090] The residuals between the theoretical electric field data and the measured electric field data at all sensor locations are calculated;

[0091] A Gaussian likelihood function is constructed based on the weighted quadratic norm of the residuals, and the Gaussian likelihood function is taken as the likelihood function; wherein the weights are determined by the covariance of the measurement noise.

[0092] The likelihood function measures the match between a given candidate fracture model and the actual observation data. In this embodiment, for a candidate fracture state, e.g. a ten-meter long fracture with a five-millimeter aperture, the theoretical electric field or voltage response at all sensor locations, e.g. sixteen electrode pairs arranged on the ground, is calculated by a forward electromagnetic field model, e.g. a finite element model. This data vector containing sixteen theoretical voltage values is compared with the sixteen voltage values actually measured in the field, and the difference between them, i.e. the residual vector, is calculated. The construction of the Gaussian likelihood function is based on the assumption that the observation error follows a Gaussian distribution. Its core part is an exponential term, which is the weighted quadratic norm of the residuals. The weights are determined by the inverse of the covariance matrix of the measurement noise. For example, if the measurement noise standard deviation of a certain sensor is 0.1 millivolt, indicating that its data is very reliable, it will have a high weight in the calculation; on the contrary, if the noise standard deviation of another sensor is 5 millivolts, its weight will be low.

[0093] More specifically, the step c1 is specifically:

[0094] The expansion rate of the fracture is estimated based on the posterior distribution at time steps k-1 and k-2, and according to the expansion rate, the components corresponding to the fracture length and aperture in the covariance matrix of the multivariate normal proposal distribution are increased.

[0095] To improve the efficiency of MCMC sampling, especially when dealing with dynamic evolution problems, an adaptive proposal distribution strategy is adopted in this embodiment. The proposal distribution is used to generate new candidate samples in the parameter space. Historical information is utilized to predict the dynamic behavior of the fracture, i.e. by comparing the mean or mode of the posterior distribution at time step k-2 and k-1 to estimate the rate of change of the main fracture parameters. For example, if the average fracture length at time step k-2 is 3.2 meters, and it grows to 3.8 meters at time step k-1, the rate of fracture propagation can be estimated to be about 0.6 per time step.

[0096] The estimated propagation rate is used to adjust the proposal distribution at time step k. The proposal distribution is usually a multivariate normal distribution, whose shape is controlled by the covariance matrix. According to the estimated propagation rate of 0.6 meters, the variance term corresponding to the fracture 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 Figure 3 This will make the proposal distribution jump a larger step in the fracture length dimension, so that it can more efficiently follow the actual expansion trend of the fracture to explore. The adaptive mechanism enables the sampling process to concentrate computational resources in the area of the parameter space where changes are most likely to occur, thereby speeding up the convergence.

[0097] More specifically, the step c2 is specifically:

[0098] A low-fidelity forward model with low computational cost is used to perform preliminary likelihood evaluation on the first-stage candidate samples;

[0099] 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 exact likelihood value, and the exact likelihood value is corrected in combination with the calculation results of the low-fidelity and high-fidelity models.

[0100] To solve the problem of long computation time of high-precision forward modeling, two kinds of forward models with different precision are used. A low-precision model, such as a finite element model using coarse grid or a simplified analytical solution, is very fast, which may only take 0.5 second. A high-precision model is a finite element model using fine grid, which is more accurate but time-consuming, which may take 5 minutes. For each candidate sample generated by the Monte Carlo sampler, a preliminary likelihood value is calculated by calling the fast low-precision model. A switching criterion is set, such as a likelihood threshold. If the preliminary likelihood value calculated by the low-precision model is lower than the threshold, it means that the candidate sample is extremely mismatched with the observed data, and the sample is directly rejected without further expensive calculation. For example, if 90% of the candidate samples are quickly filtered out by the low-precision model, a lot of calculation time can be saved. Only when the preliminary likelihood value is higher than the threshold, indicating that the sample has a certain rationality, the high-precision model is started for accurate calculation. In order to further improve the accuracy, a small number of samples with both high-precision and low-precision model results are used to establish a correction model to correct the calculation results of the high-precision model, so as to obtain a fast and accurate likelihood evaluation, thereby improving the overall calculation efficiency by tens of times on the premise of ensuring accuracy.

[0101] More specifically, the step c3 is specifically:

[0102] When the first-stage candidate sample is rejected, the gradient of the target posterior probability density at the position of the rejected sample is calculated, and a second-stage candidate sample pointing to a region with higher probability is generated from a variance-shrunk proposal distribution by using the gradient information.

[0103] 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 the embodiment, when a first-stage candidate sample is rejected, it is not immediately abandoned. Instead, the gradient of the target posterior probability density at the position of the first-stage candidate sample is calculated. The gradient is a vector pointing to the direction in which the probability density grows fastest in the parameter space, such as Figure 4The second proposal is guided by the gradient information. A second stage of candidate samples 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 in the direction of the just computed gradient from the current position. For example, if the gradient indicates that increasing the fracture opening will most rapidly increase the posterior probability, then the second proposal will have a larger opening value than the first proposal. At the same time, the covariance matrix of the second proposal distribution is shrunk, for example, its variance is one fourth of the first proposal variance, indicating a smaller local search. The second proposal has a higher probability of being accepted because it is intentionally pushed towards a region of higher probability, thus turning a failed attempt into an effective move, improving the overall efficiency of the sampler.

[0104] In a first embodiment of the present application, a real-time monitoring system for hydraulic fracture is provided, which is based on time-varying electric field dynamic data, comprising the following modules:

[0105] A data acquisition module is configured to acquire time-varying electric field data sequences collected by electric field sensors arranged in a hydraulic fracturing area;

[0106] An inversion monitoring module is configured to, for each monitoring time step k in the electric field data sequences, perform the following sequential Bayesian inversion:

[0107] a. based on the samples of the fracture state parameter posterior distribution at time step k-1, a kernel density estimation is used to construct the state parameter prior distribution at time step k; the state parameters include the length, opening, azimuth angle and inclination angle of the fracture;

[0108] b. a physical regularization term is obtained based on the regional stress field and the fracture tip stress intensity factor, and a likelihood function is obtained according to the matching degree between the forward electric field data and the measured electric field data under the given state parameters; a target posterior probability density function is established by the prior distribution, the likelihood function and the physical regularization term;

[0109] c. in the Markov chain Monte Carlo iteration, the target posterior probability density function is sampled;

[0110] d. after the iteration is completed, the sample set output by the Markov chain is taken as the fracture state parameter posterior distribution at time step k, and the three-dimensional geometric shape and distribution information of the fracture is extracted therefrom.

[0111] More specifically, the sampling of the target posterior probability density function is specifically:

[0112] c1. according to the estimated fracture propagation rate at the previous time step k-1, the covariance related to the fracture propagation in the multivariate normal proposal distribution is adjusted, and the first stage candidate sample is generated therefrom;

[0113] c2. evaluating a likelihood function of the first-stage candidate sample using a multi-fidelity strategy: calculating an initial likelihood value using a low-fidelity model, starting a high-fidelity model to calculate an accurate likelihood value when a preset switching criterion is met, and correcting the results of the two models in combination;

[0114] c3. when the first-stage candidate sample is rejected, calculating a target posterior probability density gradient at the rejected sample, and generating a second-stage candidate sample using a shrunk proposal distribution along the gradient direction.

[0115] More specifically, the kernel density estimation is used to construct the state parameter prior distribution at time step k, specifically:

[0116] The kernel density estimation method is used to construct the prior distribution at time step k that can reflect the uncertainty of crack evolution over time according to the crack state parameter posterior distribution at time step k-1.

[0117] More specifically, the physical regularization term is obtained based on the regional in-situ stress field and the crack tip stress intensity factor calculation, specifically:

[0118] The crack tip stress intensity factor of the candidate crack geometry is compared with the preset rock fracture toughness, and when the stress intensity factor exceeds the rock fracture toughness, a preset low probability penalty is applied to the candidate crack geometry.

[0119] More specifically, the likelihood function is obtained according to the matching degree of the forward electric field data and the measured electric field data under a given state parameter, specifically:

[0120] For any given state parameter sample, theoretical electric field data is calculated through a forward model;

[0121] The residuals of the theoretical electric field data and the measured electric field data at all sensor locations are calculated;

[0122] Based on the weighted quadratic norm of the residuals, the exponential term of the Gaussian likelihood function is constructed, and the Gaussian likelihood function is taken as the likelihood function; wherein the weight is determined by the covariance of the measurement noise.

[0123] More specifically, the step c1 is specifically:

[0124] The crack propagation rate is estimated based on the posterior distributions at time steps k-1 and k-2, and according to the propagation rate, the components corresponding to the crack length and opening in the covariance matrix of the multivariate normal proposal distribution are increased.

[0125] More specifically, the step c2 is specifically:

[0126] a low-cost low-fidelity forward model is used to perform a preliminary likelihood evaluation on the first-stage candidate samples;

[0127] when the preliminary likelihood evaluation result meets a preset switching criterion, a high-cost high-fidelity forward model is started to calculate an accurate likelihood value, and the accurate likelihood value is corrected in combination with the calculation results of the low-fidelity and high-fidelity models.

[0128] More specifically, the step c3 specifically includes:

[0129] when the first-stage candidate sample is rejected, a gradient of a target posterior probability density of the rejected sample position is calculated, and a second-stage candidate sample pointing to a region with higher probability is generated from a variance-shrunk proposal distribution by using the gradient information.

[0130] From the above description of the embodiments, those skilled in the art can clearly understand that the present application can be implemented by means of software and the necessary general hardware platform. Based on such an understanding, the technical solutions of the present application can be embodied in the form of a software product, which can be stored in a storage medium, such as a ROM / RAM, a magnetic disk, an optical disk, etc., and includes a number of instructions to make a computer device (which can be a personal computer, a server, or a network device, etc.) execute the methods described in the various embodiments or some parts of the embodiments.

[0131] Each of the embodiments in the specification is described in a progressive manner, and the same or similar parts between the embodiments can be referred to each other. Each embodiment mainly describes the differences from other embodiments. In particular, for the system or system embodiments, since they are basically similar to the method embodiments, the description is relatively simple, and the relevant parts can be referred to the part of the method embodiments. The above-described system and system embodiments are only illustrative, and the units described as separate components can be or can not be physically separated, and the components displayed as units can be or can not be physical units, i.e., they can be located in one place or distributed on multiple network units. According to the actual needs, some or all of the modules can be selected to achieve the purpose of the embodiment. Those skilled in the art can understand and implement without creative labor.

[0132] The above provides a method for providing commodity object information and an electronic device, and the principles and implementation modes of the present application are described by applying specific examples in this paper. The above example is only used to help understand the method and its core idea; at the same time, for those skilled in the art, according to the idea of the present application, the specific implementation mode and application range will be changed. In summary, the content of the specification should not be understood as a limitation of the present application.

Claims

1. A method for real-time monitoring of a hydraulic fracture based on time-varying electric field dynamic data, characterized in that, The method comprises the following steps: obtaining a time-varying electric field data sequence collected by an electric field sensor arranged in a fracturing work area; for each monitoring time step k in the electric field data sequence, performing the following sequential Bayesian inversion: a. based on the samples of the fracture state parameter posterior distribution at time step k-1, constructing the state parameter prior distribution at time step k using kernel density estimation; the state parameters include the length, opening, azimuth angle and inclination angle of the fracture; b. obtaining a physical regularization term based on the physical constraints calculated from the regional stress field and the stress intensity factor at the fracture tip, and obtaining a likelihood function according to the matching degree of the forward electric field data under a given state parameter and the measured electric field data; establishing a target posterior probability density function composed of the prior distribution, the likelihood function and the physical regularization term; c. sampling the target posterior probability density function in Markov chain Monte Carlo iteration; d. after iteration, taking the sample set output by the Markov chain as the fracture state parameter posterior distribution at time step k, and extracting the fracture three-dimensional geometric shape and distribution information therefrom; sampling the target posterior probability density function, specifically: c1. adjusting the covariance related to fracture propagation in the multivariate normal proposal distribution according to the estimated fracture propagation rate at the previous time step k-1, and generating a first-stage candidate sample therefrom; c2. evaluating the likelihood function of the first-stage candidate sample using a multi-fidelity strategy: using a low-fidelity model to calculate an initial likelihood value, when a preset switching criterion is met, starting a high-fidelity model to calculate an accurate likelihood value, and combining the results of the two models for correction; c3. when the first-stage candidate sample is rejected, calculating the gradient of the target posterior probability density at the rejected sample, and generating a second-stage candidate sample using a shrunk proposal distribution along the gradient direction; constructing the state parameter prior distribution at time step k using kernel density estimation, specifically: using kernel density estimation method, constructing the prior distribution at time step k which can reflect the uncertainty of fracture evolution over time according to the fracture state parameter posterior distribution at time step k-1; obtaining a physical regularization term based on the physical constraints calculated from the regional stress field and the stress intensity factor at the fracture tip, specifically: comparing the stress intensity factor of the candidate fracture geometric shape with the preset rock fracture toughness, and applying a preset low-probability penalty to the candidate fracture geometric shape when the stress intensity factor exceeds the rock fracture toughness.

2. The method of claim 1, wherein, the likelihood function is obtained according to the matching degree of the forward electric field data under a given state parameter and the measured electric field data, specifically: for any given state parameter sample, calculating the theoretical electric field data through the forward model; calculating the residual of the theoretical electric field data and the measured electric field data at all sensor positions; constructing the exponential term of the Gaussian likelihood function based on the weighted quadratic norm of the residual, and taking the Gaussian likelihood function as the likelihood function; wherein the weight is determined by the covariance of the measurement noise.

3. The method of claim 1, wherein, the step c1, specifically: The propagation rate of the fracture is estimated based on the posterior distribution at time steps k-1 and k-2, and according to the propagation rate, components in the covariance matrix of the multivariate normal proposal distribution corresponding to the length and opening of the fracture are increased.

4. The method of claim 1, wherein, The step c2 is specifically: A low-fidelity forward model with low calculation cost is used to preliminarily evaluate the likelihood of the first-stage candidate sample; When the preliminary evaluation result meets a preset switching criterion, a high-fidelity forward model with high calculation cost and high accuracy is started to calculate an accurate likelihood value, and the accurate likelihood value is corrected in combination with the calculation results of the low-fidelity and high-fidelity models.

5. The method of claim 1, wherein, The step c3 is specifically: When the first-stage candidate sample is rejected, the gradient of the target posterior probability density at the position of the rejected sample is calculated, and a second-stage candidate sample is generated from a proposal distribution with variance shrinkage and in a direction with higher probability by using the gradient information.

6. A real-time monitoring system for fracturing fractures based on time-varying electric field dynamic data, characterized in that, The method comprises the following modules: A data acquisition module is configured to acquire a time-varying electric field data sequence collected by an electric field sensor arranged in a fracturing working area; An inversion monitoring module is configured to, for each monitoring time step k in the electric field data sequence, perform the following sequential Bayesian inversion: a. Based on the samples of the fracture state parameter posterior distribution at time step k-1, a kernel density estimation is used to construct the state parameter prior distribution at time step k; the state parameter includes the length, opening, azimuth angle and inclination angle of the fracture; b. A physical regularization term is obtained based on a regional in-situ stress field and a fracture tip stress intensity factor, and a likelihood function is obtained according to the matching degree of the forward electric field data under a given state parameter and the measured electric field data; a target posterior probability density function is established by the prior distribution, the likelihood function and the physical regularization term; c. In Markov chain Monte Carlo iteration, the target posterior probability density function is sampled; d. After the iteration ends, the sample set output by the Markov chain is taken as the fracture state parameter posterior distribution at time step k, and the three-dimensional geometric shape and distribution information of the fracture is extracted from the sample set; The sampling of the target posterior probability density function is specifically: c1. According to the estimated fracture propagation rate at the previous time step k-1, the covariance of the multivariate normal proposal distribution related to the fracture propagation is adjusted, and a first-stage candidate sample is generated therefrom; c2. A multi-fidelity strategy is used to evaluate the likelihood function of the first-stage candidate sample: an initial likelihood value is calculated by using a low-fidelity model, when a preset switching criterion is met, an accurate likelihood value is calculated by starting a high-fidelity model, and the accurate likelihood value is corrected in combination with the results of the two models; c3. When the first-stage candidate sample is rejected, the gradient of the target posterior probability density at the position of the rejected sample is calculated, and a second-stage candidate sample is generated from a proposal distribution with variance shrinkage and in a direction with higher probability by using the gradient information; The kernel density estimation is used to construct the state parameter prior distribution at time step k, and the construction is specifically: The kernel density estimation method is used to construct the prior distribution at time step k, which can reflect the uncertainty of the fracture evolution over time, according to the fracture state parameter posterior distribution at time step k-1. The physical regularization term is obtained based on the regional geostress field and the physical constraint of the crack tip stress intensity factor calculation, and is specifically: The crack tip stress intensity factor of the candidate crack 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 applied to the candidate crack geometry.

Images