A hydraulic fracturing monitoring system and method based on high-frequency pressure water hammer wave and cepstrum comprehensive analysis

Abstract

The application discloses a hydraulic fracturing monitoring system and method based on high-frequency pressure water hammer wave and cepstrum comprehensive analysis, and relates to the technical field of oil and gas field hydraulic fracturing monitoring. The technical points are as follows: the method comprises distributed multi-modal intelligent sensor front-end data acquisition and preprocessing, cepstrum signal enhancement and feature extraction based on physical constraints, hierarchical ground stress driven mechanical model structure deviation detection and self-adaptive switching inversion, dynamic basic probability distribution correction and evidence conflict processing of self-diagnosis quality factors, cross-module expert system comprehensive reasoning and result output; relying on the integrated process of multi-source data fusion and multi-model linkage analysis, the wave velocity accurate measurement, cepstrum feature extraction, fracture parameter iteration, bridge plug sealing evaluation and comprehensive reasoning are synchronously completed, the two problems of traditional technology wave velocity calculation distortion and large sealing determination error are solved from the root, and the overall precision and integrity of the whole process monitoring of hydraulic fracturing are greatly improved.

Description

Technical Field

[0001] This invention relates to the field of hydraulic fracturing monitoring technology in oil and gas fields, specifically to a hydraulic fracturing monitoring system and method based on high-frequency pressure water shock wave and cepstrum comprehensive analysis. Background Technology

[0002] Hydraulic fracturing is a key technology for the development of unconventional reservoirs such as shale oil and gas and tight oil and gas. The geometric morphology of the fracture network formed during fracturing, such as the number, location, height, and fluid injection distribution of fractures, directly determines the production enhancement effect. Accurately obtaining these parameters is crucial for optimizing fracturing design and evaluating well completion quality.

[0003] Currently, fracturing diagnostic technologies mainly include microseismic monitoring, logging, tracers, and bottom hole pressure analysis. Microseismic monitoring is costly and susceptible to noise interference. Logging and tracer methods can only be implemented after fracturing and cannot provide real-time feedback. Traditional pressure analysis mainly relies on seepage theory and is difficult to directly obtain the number of fractures and fluid inflow distribution.

[0004] In deep shale gas multi-stage fracturing operations, the fluid inside the wellbore is often accompanied by pressure drop and degassing, and the formation contains complex stress shielding layers. In this scenario, existing technologies suffer from the following problems: First, in pressure drop and degassing scenarios, the volume fraction of air bubbles is non-uniformly distributed along the wellbore depth, causing local abrupt changes in water hammer wave velocity with depth. Existing technologies use constant wave velocity for cepstral time-frequency conversion and phase alignment, which leads to severe phase distortion and energy dissipation of the reflected signal in the cepstral domain. This causes the weak fracture reflection characteristics to be submerged by background noise, making it impossible to accurately extract the true spatial location and opening state of the fracture. Second, when the fracture propagates through a shielding layer with a significant stress gradient, the fracture morphology changes from radial propagation to restricted propagation. The original single model structure assumptions fundamentally deviate from the actual physical process. Existing technologies often use fixed single mechanical models, lacking real-time detection and adaptive switching mechanisms for model structure deviations. This causes iterative calculations to fall into pseudo-convergence, resulting in a systematic deviation between the final output fracture height and total volume, thus misleading the comprehensive diagnosis of the overall fracturing effect. Summary of the Invention

[0005] To achieve the above objectives, the present invention provides the following technical solution: A hydraulic fracturing monitoring method based on high-frequency pressure water hammer wave and cepstrum analysis includes the following steps: acquiring high-frequency pressure, in-situ stress, temperature, and vibration data from the wellbore; performing adaptive baseline drift correction and multi-source fusion estimation on the high-frequency pressure to obtain the reference wave velocity, data quality factor, pure pressure sequence, and water hammer wave two-way propagation time; based on the reference wave velocity and water hammer wave two-way propagation time, performing cepstrum transformation, dynamic time warping alignment, and feature extraction on the pure pressure sequence to obtain fracture reflection feature locations and multiple profiles, including a spatial wave velocity field profile, a bubble volume fraction profile, a fracture opening width profile, and a fluid inflow intensity profile; and combining in-situ stress and multiple profiles to perform structural deviation detection and adaptive analysis on the mechanical model for hydraulic fracturing monitoring. The inversion process is switched to obtain the crack height, average crack width, and total crack volume. Multiple types of bridge plug seal assessment evidence are defined and assigned initial basic probabilities. The basic probabilities are dynamically corrected using a data quality factor. The conflict coefficient between evidence is calculated, and conflict is redistributed when it exceeds the conflict threshold. The posterior probabilities of good seal, seal failure, and uncertainty are fused to output the bridge plug seal status classification. Based on the pump stop pressure drop curve, G-function closure pressure analysis is performed to obtain the closure pressure and filtration coefficient. The crack reflection characteristics, inflow intensity profile, crack height, bridge plug seal status classification, closure pressure, and filtration coefficient are transformed into fact objects with confidence levels. A pre-built inference rule base is loaded for forward chain inference to output a comprehensive diagnostic conclusion.

[0006] Furthermore, the pressure cross-correlation wave velocity estimator, vibration propagation wave velocity estimator, and theoretical wave velocity estimator are run in parallel to obtain preliminary wave velocity estimates. The preliminary wave velocity estimates are then substituted into the discretized mass and momentum conservation equations for physical verification. The residuals between the predicted and measured pressure values ​​are calculated, and the estimated value that passes the verification and has the smallest residual is selected as the reference wave velocity. If all preliminary wave velocity estimates fail the verification, the wave velocity prediction model based on a long short-term memory network is invoked. The time series data within a fixed time window before pump shutdown is used as input, and the predicted wave velocity value is output as the reference wave velocity.

[0007] Furthermore, the dynamic time warping alignment and feature extraction process includes: using the cepstrum profile of the first periodic window as an initial reference, performing forward alignment with the cepstrum profile of subsequent periodic windows through a dynamic time warping algorithm to obtain a forward cumulative warping mapping function; using the periodic window with the highest signal-to-noise ratio as a reference, performing reverse alignment through a dynamic time warping algorithm to obtain a reverse cumulative warping mapping function; extracting the local mapping ratio from the cumulative warping mapping function, calculating the local wave velocity varying along the well depth based on the reference wave velocity and the local mapping ratio to form a spatial wave velocity field profile; comparing the forward cumulative warping mapping function and the reverse cumulative warping mapping function, if the difference in displacement between the two at the same water hammer wave two-way propagation time position exceeds a preset difference threshold, then using the wave velocity value in the spatial wave velocity field profile as a constraint for weighted correction to obtain the final water hammer wave two-way propagation time axis aligned profile.

[0008] Furthermore, the mechanical model includes the PKN mechanical model, the KGD mechanical model, and the P3D mechanical model; the process of structural deviation detection and adaptive switching inversion processing of the mechanical model for hydraulic fracturing monitoring includes: constructing a stress profile based on the geostress data collected along the wellbore, calculating the stress gradient to locate the stress shielding layer; starting from a preset initial mechanical model and an initial fracture height value, using an ensemble Kalman filter algorithm to assimilate the pressure transient features extracted from the high-frequency pressure data and dynamically update the mechanical model parameters; after each iteration, detecting whether the fracture height boundary predicted by the current mechanical model continuously exceeds the located stress shielding layer boundary, and whether the relative error between the predicted fracture volume and the actual target volume exceeds the relative error threshold; if both exceed the relative error threshold, it is determined that the current mechanical model has deviated structurally; calling the spatial wave velocity field profile for wave velocity and stress dual constraint verification; if the verification confirms structural deviation, then automatically switching to different fracture propagation mechanical models to continue iteration based on the average stress gradient and thickness of the stress shielding layer, to obtain the fracture height, average fracture width, and total fracture volume.

[0009] Furthermore, it also includes a physical constraint fracture propagation early warning based on bubble volume fraction profile and fracture opening width profile. If, in a region outside the fracture height boundary of the current iteration, the inverted bubble volume fraction exceeds a first preset threshold and the fracture opening width exceeds a second preset threshold, then the region is recorded as an abnormal region, and the depth location of the abnormal region is further determined: if the abnormal region is located inside the reservoir, then in the next iteration, the search range of the fracture height is expanded towards the direction of the abnormal region; if the abnormal region exceeds the reservoir interlayer, then it is determined to be a risk of excessive fracture propagation, the search range is no longer expanded, and an early warning is output; if, within the fracture height range of the current iteration, the fracture opening width changes non-monotonically along the depth direction, then wave velocity stress dual constraint verification or pseudo-convergence detection is triggered.

[0010] Furthermore, the process of dynamically correcting the basic probability using data quality factors includes: presetting high-factor thresholds, medium-factor thresholds, and low-factor thresholds; if the data quality factor > high-factor threshold, increasing the basic probability allocation value of the most probable state proposition in the evidence body; if the low-factor threshold ≤ data quality factor < medium-factor threshold, decreasing the basic probability allocation value of all propositions in the evidence body according to a linear relationship; if the medium-factor threshold ≤ data quality factor ≤ high-factor threshold, keeping the initial basic probability allocation unchanged; if the data quality factor < low-factor threshold, determining that the evidence body is unusable and transferring all its basic probability allocations to uncertain propositions; the conflict redistribution process includes: the conflict coefficient is the sum of the probability products of mutually exclusive propositions between evidence bodies, uncertain propositions do not participate in the conflict coefficient calculation; if the conflict coefficient is greater than the preset conflict threshold, using a weighted correction method to preprocess the evidence probability before performing evidence synthesis, while eliminating evidence bodies with a data quality factor lower than the medium-factor threshold.

[0011] Furthermore, the process of performing G-function closure pressure analysis based on the pump shutdown pressure drop curve includes: converting the actual time after pump shutdown into a dimensionless G-function; plotting the curve of the derivative of pressure with respect to the G-function as a function of change, identifying a stable linear decreasing segment within the range where the G-function value is greater than a preset value; performing linear regression on the linear decreasing segment, and extrapolating the regression line to intersect the axis where the pressure derivative is equal to zero, obtaining the G-function value corresponding to the intersection point; and finding the pressure value corresponding to the G-function value at the intersection point on the curve of pressure changing with the G-function, which is taken as the closure pressure.

[0012] Furthermore, the forward chain reasoning process includes: sequentially scanning the rules in the reasoning rule base, matching whether the condition part of the rule matches the fact objects in the working memory; if the match is successful, triggering the conclusion part of the rule, adding the new fact and its propagated confidence to the working memory, wherein the confidence of the new fact is equal to the rule confidence multiplied by the minimum value of the confidence of each fact object in the condition; repeating the scanning and matching until no new facts are generated, and outputting a comprehensive diagnostic conclusion including fracture effectiveness evaluation, bridge plug sealing comprehensive risk level and fracturing effect comprehensive score.

[0013] A hydraulic fracturing monitoring system integrating high-frequency pressure water hammer and cepstrum analysis includes: a data acquisition and wave velocity fusion verification module, which acquires high-frequency pressure, geostress, temperature, and vibration data from the wellbore, performs adaptive baseline drift correction, and obtains the reference wave velocity, data quality factor, pure pressure sequence, and water hammer two-way propagation time through multi-source fusion estimation and physical verification; a physically constrained cepstrum enhancement and feature inversion module, which performs cepstrum transformation and dynamic time warping alignment based on the pure pressure sequence and reference wave velocity, inverts the spatial wave velocity field profile, and extracts fracture reflection features, bubble volume fraction profile, fracture opening width profile, and fluid injection intensity profile by combining neural network classification and acoustic attenuation verification; and a mechanical model deviation detection and adaptive inversion module, which locates stress by combining geostress data and inversion profile. The shielding layer uses a filtering algorithm to assimilate parameters. When a deviation from the mechanical model structure is detected, it performs wave velocity and stress dual-constraint verification and adaptively switches the mechanical model, outputting crack height, average crack width, and total crack volume. The evidence dynamic correction and conflict fusion module defines multiple types of bridge plug seal evaluation evidence bodies, dynamically corrects the basic probability allocation using data quality factors, calculates the conflict coefficient, and performs conflict redistribution when the threshold is exceeded. It fuses and outputs the bridge plug seal state classification and performs closure pressure analysis based on the pump stop pressure drop curve to obtain the closure pressure and filtration coefficient. The cross-level expert rule reasoning and comprehensive diagnosis module transforms the inversion profile, crack geometric parameters, seal state classification, and closure pressure parameters into fact objects with confidence, loads the rule base, performs forward chain reasoning, and outputs a comprehensive diagnostic conclusion.

[0014] Furthermore, in the data acquisition and wave velocity fusion verification module, pressure cross-correlation, vibration propagation, and theoretical wave velocity estimators are run in parallel. The residuals are calculated by substituting them into the discretized mass and momentum conservation equations, and the minimum residual value is selected as the reference wave velocity. If all values ​​exceed the limit, the long short-term memory network is called to output the predicted wave velocity value. In the physical constraint cepstral enhancement and feature inversion module, forward and reverse dynamic time warping alignment chains are established. The local mapping ratio is extracted to calculate the corrected local wave velocity to form a spatial wave velocity field profile. A lightweight convolutional neural network is used to separate features, and the profile is generated through acoustic attenuation model verification and multi-sensor spatial collaborative constraints. In the mechanical model deviation detection and adaptive inversion module, a state vector assimilation high-frequency pressure transient feature is constructed. When the crack height boundary continuously exceeds the stress shielding layer boundary and the volume error exceeds the limit, The system determines structural deviations and automatically switches the mechanical model based on stress shielding layer parameters to continue iterating until the dual convergence constraints are met. In the evidence dynamic correction and conflict fusion module, probability enhancement, maintenance, decay, or transfer corrections are performed based on the data quality factor interval. The sum of the probabilities of mutually exclusive propositions between different evidence bodies is calculated as the conflict coefficient. If the threshold is exceeded, a weighted correction method is used for preprocessing before evidence synthesis. The pressure and derivative curves of the dimensionless G-function are plotted and extrapolated to obtain the closing pressure. In the cross-level expert rule reasoning and comprehensive diagnosis module, the fact objects in the rule matching working memory are scanned sequentially. If a match is successful, the new fact and its propagated confidence are added to the working memory. The confidence propagation is based on the product of the rule confidence and the minimum confidence of the fact in the condition. The scan is repeated until no new facts are generated and a diagnostic report is output.

[0015] This invention provides a hydraulic fracturing monitoring system and method based on high-frequency pressure water hammer wave and cepstral analysis, which has the following advantages: 1. This invention obtains dynamic reference wave velocity through multi-source fusion estimation and combines it with dynamic time warping alignment to effectively eliminate phase distortion caused by the non-uniform distribution of water hammer wave velocity due to fluid degassing, and accurately extracts crack reflection characteristics, spatial wave velocity field profile and bubble volume fraction profile. At the same time, it introduces a mechanical model structure deviation detection and adaptive switching inversion mechanism, which automatically switches the model when the crack penetrates the stress shielding layer, avoids false convergence, eliminates systematic deviations in the inversion of crack geometric parameters (height, width, volume), and improves inversion reliability.

[0016] 2. This invention dynamically corrects the basic probability allocation of the DS evidence theory based on data quality factors and performs conflict redistribution. It fuses multiple types of bridge plug sealing assessment evidence, effectively suppressing the interference of low-quality data or abnormal sensors on diagnostic results. Furthermore, it combines G-function closure pressure analysis to obtain the true closure pressure and filtration coefficient, providing accurate quantitative basis for the actual filtration behavior of fracturing fluid in the fracture. This upgrades the assessment of bridge plug sealing status from a single probability inference to a reliable graded early warning supported by multi-dimensional evidence.

[0017] 3. This invention constructs a cross-module expert-level comprehensive reasoning and decision support capability. Through forward chain reasoning, it transforms inversion parameters, evaluation results, and G-function analysis conclusions into unified fact objects, loads a cross-module rule base covering crack effectiveness, construction risks, and effect evaluation for automatic reasoning, and outputs a comprehensive score for fracturing effect, early warning of micro-leakage and excessive crack propagation, and construction optimization suggestions to form a comprehensive diagnostic conclusion that can directly guide on-site operations. Attached Figure Description

[0018] Figure 1 This is an overall flowchart of the method in Embodiment 1 of the present invention.

[0019] Figure 2 This is a flowchart of the mechanical model inversion and adaptive switching of the present invention.

[0020] Figure 3 This is a flowchart of the multi-evidence dynamic fusion and comprehensive evaluation process of the present invention.

[0021] Figure 4 This is a system framework diagram of Embodiment 2 of the present invention. Detailed Implementation

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

[0023] Example 1: Please see Figures 1 to 3 This embodiment provides a hydraulic fracturing monitoring method based on high-frequency pressure water hammer and cepstral comprehensive analysis. The system mentioned in this method refers to the computing and processing system used to execute this hydraulic fracturing monitoring method, which includes an edge computing unit with distributed intelligent sensor nodes and a central processing server. This method includes steps 1 to 5.

[0024] Step 1: Data acquisition and preprocessing of distributed multimodal intelligent sensors at the front end: Distributed multimodal intelligent sensor arrays are deployed at key locations such as the wellhead of oil and gas wells, fracturing sections, and above and below bridge plugs, with all sensor nodes starting up synchronously.

[0025] Step 1.1: Synchronous acquisition of raw data from multiple sources.

[0026] Each sensor node integrates a high-frequency pressure sensing unit, a micro-in-situ stress sensing unit, a temperature sensing unit, a triaxial vibration sensing unit, and an edge computing unit. The high-frequency pressure sensing unit continuously acquires high-frequency pressure time-series data within the wellbore at a sampling rate of 1000 times / s. The micro-in-situ stress sensing unit collects vertical in-situ stress, maximum horizontal principal stress, minimum horizontal principal stress, and Young's modulus data at various depths in segments. The temperature sensing unit collects real-time data on the wellbore fluid temperature and the wellbore ambient temperature. The triaxial vibration sensing unit collects triaxial vibration acceleration data at a sampling rate of 2000 times / s.

[0027] Step 1.2: Multi-source fusion estimation and physical verification of water hammer wave propagation velocity.

[0028] The edge computing units run three preliminary wave velocity estimators in parallel: Pressure cross-correlation wave velocity estimator: Calculates the cross-correlation function of the pressure signal sequences collected by two adjacent nodes, finds the delay time corresponding to the peak value of the cross-correlation function, and divides the known distance between the two nodes by this delay time to obtain the estimated cross-correlation wave velocity. Vibration propagation wave velocity estimator: Based on the arrival time difference of the peak value of fluid disturbance in the vibration acceleration signal; records the time when the same vibration peak value arrives at two vibration sensors at different depths, and divides the depth difference by the time difference to obtain the estimated vibration propagation wave velocity. Theoretical wave velocity estimator: Based on the real-time temperature, pressure, and gas content of the current well fluid, and using fluid property models such as the Wood equation, calculates the theoretical water hammer wave velocity to obtain the estimated theoretical wave velocity.

[0029] Subsequently, a physical verification model is invoked to filter the above estimates. This physical verification model consists of a set of discretized mass and momentum conservation equations. The pressure gradient at the current moment, the fluid density converted from temperature and pressure, and the three preliminary wave velocity estimates are substituted into the physical verification model to calculate the pressure prediction value after the next small time interval. The pressure prediction value corresponding to each preliminary wave velocity estimate is compared with the actual collected pressure value to calculate the residual. The residual threshold is set to 2% of the pressure measurement range. This is an empirical value based on the accuracy requirements of conventional pressure sensors and on-site measurement statistics, used to determine whether the estimated value conforms to physical laws: if the residual is ≤2%, it means that the preliminary wave velocity estimate meets physical consistency and the estimation result is valid; if the residual is >2%, it means that the deviation is too large. The estimated value with the smallest residual is determined to be the one that best conforms to physical laws and is used as the global reference wave velocity output at the current moment, denoted as the reference wave velocity.

[0030] If the residuals calculated after substituting the three preliminary wave velocity estimates into the physical verification model are all greater than 2% of the pressure measurement range, the current operating condition is determined to be abnormal. At this time, the system automatically switches to standby mode and calls the wave velocity prediction model built based on a long short-term memory network. This model uses the time series data of pressure, temperature, and displacement within a fixed 5-second time window before pump shutdown as the input sequence, extracts time-dependent features through long short-term memory units, and outputs the predicted wave velocity value as a supplementary value to the baseline wave velocity. The training data of this model comes from the measured wave velocity data of historical fractured wells, and the training objective is to minimize the mean square error between the predicted value and the measured value. If there is at least one residual of the estimated value ≤ 2% of the pressure measurement range, the preliminary wave velocity estimate with the smallest residual is directly selected as the baseline wave velocity, and no abnormality processing is triggered.

[0031] Step 1.3: Adaptive baseline drift correction and data quality factor calculation.

[0032] The edge computing unit performs adaptive baseline drift correction on the high-frequency pressure time series data; a five-point window mid-range filtering algorithm is used to remove spike noise generated by sand impact, etc.: five consecutive sampling points are selected, their amplitudes are sorted from smallest to largest, and the median value is used to replace the amplitude of the current sampling point; a local weighted regression method is used to fit the low-frequency trend term of the pressure signal, and the bandwidth parameter of the local weighted regression is taken as 1 / 3 of the window width; the original pressure sequence is subtracted from the trend term to obtain the pure pressure sequence.

[0033] Calculate the signal-to-noise ratio (SNR): Take the root mean square (RMS) of the amplitude of the pure pressure sequence during the water hammer oscillation period after pump shutdown, minus the mean value of that period, as the effective value of the signal. The mean value of that period is the arithmetic mean of the pure pressure sequence. Take the RMS of the pure pressure sequence during the static period before fracturing operations begin as the effective value of the noise. Square the effective value of the signal and the effective value of the noise, calculate the ratio, take the common logarithm, and multiply by 10 to obtain the signal-to-noise ratio in decibels.

[0034] Calculate the signal distortion index: Statistically measure the amplitude of the pure pressure sequence within a preset time window, with a window length of 10 seconds; calculate the mean and standard deviation of the pressure amplitude within the window; count the number of sampling points whose amplitude exceeds the mean plus three times the standard deviation; divide this number by the total number of sampling points in the window to obtain the signal distortion index.

[0035] The signal-to-noise ratio (SNR) and signal distortion index are combined to form the data quality factor. If the signal distortion index > 20%, the mandatory rule is applied first, and the data quality factor is directly forced to 0.15. When the signal distortion index ≤ 20%, it is calculated according to the formula: Data quality factor = Normalized signal-to-noise ratio × 0.6 + (1 - Signal distortion index) × 0.4. Wherein, the normalized signal-to-noise ratio is equal to the signal-to-noise ratio value divided by 30dB and truncated to the interval between 0 and 1. 30dB is the maximum SNR statistical value of the field-measured pressure signal under ideal working conditions. Values ​​exceeding 30dB are counted as 1. 0.6 and 0.4 are weighting coefficients, determined based on sensitivity analysis of historical data from 30 field-tested fractured wells. The data quality factor ranges from 0 to 1.

[0036] Data quality factors > 0.8 are marked as high confidence, data quality factors 0.4 ≤ data quality factor ≤ 0.8 are marked as medium confidence, and data quality factors < 0.4 are marked as low confidence. When the signal distortion index exceeds 20%, regardless of the calculated value of the data quality factor, the data for that period is forcibly marked as low confidence, and an excessive distortion note is added to the quality label. The confidence thresholds of 0.4 and 0.8 are based on the following: field statistics show that exceeding this ratio usually means that the sensor has been subjected to severe sand impact or fluid impact interference.

[0037] The data quality factor is packaged as a quality label and uploaded alongside the original measurement value.

[0038] Step 1.4: Measurement and verification of the two-way propagation time of water hammer waves.

[0039] The node measures the two-way propagation time of the water hammer wave in real time. The measurement method is as follows: at the instant the pump stops, the moment when the sensor detects the pressure rise corresponding to the main reflecting surface of the fracture at the bottom of the well is recorded as the start time; the moment when the sensor subsequently detects the first pressure rise reflected back from the bottom of the well to the wellhead is recorded as the return time; the time difference between the two is the two-way propagation time of the water hammer wave. The measured time is compared with the theoretical two-way propagation time of the water hammer wave, which is equal to twice the wellbore length divided by the reference wave velocity. If the relative deviation between the measured value and the theoretical value exceeds 5%, the measured time is discarded, and the wave velocity estimation process in step 1.2 is repeated, without using the measured time as a constraint condition. If the deviation still exceeds 5% after the second estimation, the wellbore condition is judged to be abnormal, such as degassing or leakage, a warning is marked, and the forced iteration is terminated. The deviation threshold of 5% is set according to the normal fluctuation range statistically observed on site. Exceeding this range indicates that there may be a large error in the wave velocity estimation.

[0040] Each node packages and uploads the baseline wave velocity, data quality factor, clean pressure sequence, and other multi-source feature data.

[0041] Step 2: Cepstral signal enhancement and feature extraction based on physical constraints: Using the pure pressure sequence, reference wave velocity, water hammer wave two-way propagation time and data quality factor output from Step 1 as input, complete the cepstral signal processing.

[0042] Step 2.1: Dynamic periodic window division based on reference wave velocity.

[0043] The duration of a complete water hammer wave cycle is equal to the two-way propagation time of the water hammer wave. The start and end times of each cycle window are calculated in real time based on the reference wave velocity: the start of the first window is the pump stop time, and the end is the pump stop time plus the two-way propagation time of the water hammer wave; the start of the second window is the end of the first window, and the end is the start plus the two-way propagation time of the water hammer wave; and so on, dividing the pressure signal after the pump stops into multiple continuous cycle windows.

[0044] Step 2.2: Cepstral transformation for each periodic window.

[0045] A cepstral transform is performed on the pure pressure signal within each periodic window to obtain the cepstral profile of that window. The cepstral transform is calculated as follows: First, a fast Fourier transform is performed on the pressure signal within the window to obtain the spectrum; the square of the spectrum amplitude is taken to obtain the power spectrum; the natural logarithm of each value of the power spectrum is taken to obtain the logarithmic power spectrum; then, an inverse fast Fourier transform is performed on the logarithmic power spectrum to obtain the cepstral sequence. The independent variable of the cepstral sequence is called the inverted frequency. When there is a reflected wave with a two-way propagation time of τ in the pressure signal, a peak will appear at τ in the cepstral. The horizontal axis corresponding to the cepstral peak value is equal to the two-way propagation time of the reflected signal, and the vertical axis reflects the relative intensity of the reflected signal.

[0046] Step 2.3: Inversion and phase alignment of the spatial non-uniform wave velocity field of water hammer waves based on dynamic time warping.

[0047] While completing the cepstral profile phase alignment, the spatial distribution of water hammer wave velocity along the well depth is inverted, and the wave velocity anomaly is quantitatively converted into fracture opening width and bubble volume fraction.

[0048] Step 2.3.1: Establish forward and reverse dynamic time warping alignment chains.

[0049] The cepstrum profile of the first period window is used as the initial reference. The cepstrum profile of the second period window is dynamically time-warped and aligned with that of the first period window. The dynamic time-warping algorithm finds a curved path between these two cepstrum profile sequences, minimizing the cumulative Euclidean distance between corresponding points on the two sequences along the two-way propagation time axis of the water hammer wave. This curved path consists of a set of coordinate pairs, each representing a two-way propagation time position of the water hammer wave in the first sequence that matches a two-way propagation time position of the water hammer wave in the second sequence. The warping mapping function of the second period relative to the first period is derived from this curved path. The cepstrum profile of the second period is stretched or compressed according to this mapping function and used as a new reference. The cepstrum profile of the third period window is dynamically time-warped and aligned with the aligned second period to obtain the warping mapping function of the third period. This process is repeated to complete the forward alignment of all periods. The cumulative result of the warping mapping function of each period during the forward alignment process is called the forward cumulative warping mapping function.

[0050] Simultaneously, a reverse serial alignment chain is established: the period with the highest signal-to-noise ratio is selected from all periods as the reference period, denoted as the stable reference period; the adjacent periods after this period are dynamically time-warped and aligned with it to obtain the reverse warping mapping function, and the reverse alignment of all periods is completed in sequence; the cumulative result of the reverse warping mapping function of each period during the reverse alignment process is called the reverse cumulative warping mapping function.

[0051] Step 2.3.2: Invert the spatial non-uniform wave velocity field from the mapping function.

[0052] The cumulative mapping function relative to the reference window for each period window is extracted. This function describes the new position of each water hammer wave's two-way propagation time position on the reference cepstral profile after stretching or compression, mapped onto the current period-aligned profile. The first-order difference of this mapping function is calculated point by point to obtain the local mapping ratio. If the local mapping ratio corresponding to a certain water hammer wave two-way propagation time interval is <1, it indicates that the interval is compressed in the current period, meaning that the actual propagation time of the water hammer wave on that propagation path is shorter, and therefore the actual wave speed is higher than the reference wave speed. If the local mapping ratio is >1, it indicates that the wave speed is stretched, meaning that the actual wave speed is lower than the reference wave speed. If the local mapping ratio is =1, it indicates that the wave speed of that segment is consistent with the reference wave speed, and no correction is required.

[0053] The iterative correction method is used to invert the non-uniform spatial wave velocity field, as follows: Initial depth estimation: First, assume uniform wave velocity and use the reference wave velocity as a reference value to convert the two-way propagation time of the water hammer wave into the downhole depth; the conversion relationship is: the downhole depth equals the reference wave velocity multiplied by the one-way propagation time, and the one-way propagation time is half of the two-way time; this depth value is the initial depth; Local wave velocity calculation: According to the local mapping ratio, calculate the corrected local wave velocity point by point; the correction formula is: the local wave velocity equals the reference wave velocity divided by the local mapping ratio corresponding to the depth; Depth re-integration: Re-integrate the calculated local wave velocity along the propagation time to obtain a new depth value; the integration method is: the depth equals the integral of the local wave velocity with respect to the one-way propagation time; Iterative convergence: Repeat the local wave velocity calculation and depth re-integration until the relative change of the root mean square between the depth values ​​calculated by two adjacent iterations is less than 1%; at this point, the depth and local wave velocity are considered to be consistent; the local wave velocity obtained after convergence is plotted as a continuous profile curve along the depth direction, called the spatial wave velocity field profile; this profile also outputs the corresponding converged depth coordinates.

[0054] Step 2.3.3: Wave velocity mutation detection and physical interpretation.

[0055] The spatial wave velocity field profile is scanned along the depth direction to detect local abrupt changes in wave velocity. For each pair of adjacent depth points, the relative rate of change of wave velocity is calculated, which is the absolute value of the difference between the wave velocity of the later point and the wave velocity of the earlier point divided by the wave velocity of the earlier point. If the relative rate of change of wave velocity exceeds 5%, the location is marked as a wave velocity abrupt change zone. If it does not exceed 5%, it is considered as normal background fluctuation and is not marked. The scan continues to the next depth point. This 5% change threshold is set based on the superposition of the wave velocity change range caused by the typical fracturing fluid temperature gradient and the background noise level. Exceeding this value is considered a significant anomaly.

[0056] Analyze the characteristics of abrupt wave velocity changes: If the wave velocity decreases by more than 10% relative to the initial wave velocity of the abrupt change zone within three consecutive depth points, and this downward trend does not immediately rebound at subsequent depth points, it is determined that there is a sudden drop in wave velocity at this location. This phenomenon indicates the presence of a large amount of bubble precipitation, which causes a significant decrease in sound velocity, indicating that the pressure drop induced by the opening of the fracture has induced degassing. If the wave velocity exhibits a non-monotonic local disturbance, such as decreasing first and then increasing or increasing first and then decreasing, it is determined that there is a discontinuity in the medium at this location, i.e., the opening or closing boundary of the fracture. If the wave velocity shows a unidirectional increasing or decreasing trend along the entire wellbore, it is determined to be a background drift caused by changes in the temperature field or the fluid compressibility coefficient. After fitting this trend, it is subtracted from the spatial wave velocity field profile to obtain the residual wave velocity anomaly. The detected wave velocity abrupt change zone location is output together with the depth coordinates.

[0057] Step 2.3.4: Consistency verification and correction of forward and reverse alignment results based on wave velocity physical constraints.

[0058] The forward cumulative regularization mapping function and the reverse cumulative regularization mapping function are compared point by point. At each position of the two-way propagation time of the water hammer wave, the displacements of the forward and reverse cumulative regularization mapping functions are the forward and reverse displacements, respectively. If the difference between the forward and reverse displacements exceeds a preset difference threshold, the wave velocity value in the spatial wave velocity field profile is used as an external constraint to correct the position. The difference threshold is half of the water hammer wave two-way propagation time sampling interval of 0.005s, i.e., 0.0025s. The correction rule is that the final displacement is equal to the weighted average of the forward and reverse displacements. The weight is proportional to the displacement deviation. The larger the deviation, the higher the weight of the expected displacement calculated based on the local wave velocity. The corrected alignment result is the final water hammer wave two-way propagation time axis alignment profile.

[0059] Step 2.4: Adaptive construction of short-time cepstrum.

[0060] The water hammer wave period dynamically changes with the reference wave velocity and well length, therefore the short-time window length is also an adaptively adjusted value. Based on the long-axis aligned profile of the final two-way propagation time of the water hammer wave, it is transformed into a two-dimensional time-frequency distribution map, called a short-time cepstrum, through a short-time sliding window transformation. During the transformation, the cepstrum sequence in the long-axis aligned profile of the final two-way propagation time of the water hammer wave is first segmented along the time axis according to the short-time window length. The default value of the short-time window length is 1 / 2 of the water hammer wave period. A fast Fourier transform is performed on the cepstrum segment in each short-time window to obtain the distribution of the cepstrum amplitude with frequency in that window. The spectra of all windows are arranged in chronological order to form a two-dimensional matrix, with time on the horizontal axis and frequency on the vertical axis. The matrix element values ​​are the cepstrum amplitude intensities at that time and frequency position.

[0061] Using the measured two-way propagation time of water hammer waves in the well as a physical benchmark, the expected occurrence range of the fracture reflection signal in the short-time cepstrum is estimated. The two-way propagation time of water hammer waves corresponding to fracture reflection should be equal to the two-way time of the sound wave propagating from the wellhead to the fracture location and back to the wellhead. This time value is between 0 and the two-way propagation time of water hammer waves. Based on this expected range and the current data quality factor, the optimal short-time Fourier transform parameters are adaptively selected.

[0062] The intensity of the reflected signal is determined by the peak amplitude A of the cepstral amplitudes of the previous few cycles: the adjacent region is a rectangular area enclosed by one water hammer cycle before and after the peak time position and five frequency sampling points above and below the frequency position, excluding the peak point itself. The arithmetic mean of all cepstral amplitudes within this region is the average amplitude A1 of the adjacent region; A > 3A1 indicates strong reflection, 1.5A1 ≤ A ≤ 3A1 indicates moderate reflection, and A < 1.5A1 indicates weak reflection. If the peak is at the beginning or end of the sequence and lacks a complete cycle on one side, only the existing one-sided region is used to calculate the average amplitude; if there are no complete cycles on either side, and the total number of cycles is less than 2, the peak is judged to have weak reflection by default.

[0063] When the reflection is strong and the current data quality factor is >0.8, a longer short time window length is used to obtain higher frequency resolution, with the window length being 1 / 4 of the water hammer wave period; when the data quality factor is <0.5 or the reflection is weak, the window length is automatically shortened to 1 / 8 of the water hammer wave period, and the overlap rate of adjacent windows is increased to 75% to balance time resolution and frequency resolution; for other cases, the default window length is used, which is 1 / 2 of the water hammer wave period, and the overlap rate is kept at 50%.

[0064] A 50% overlap rate is the default setting for standard short-time Fourier transform, ensuring a balance between time resolution and computational efficiency. A 75% overlap rate is used to increase the smoothness between windows when dealing with low-quality data; empirically, every 25% increase in overlap rate can reduce time-domain jitter by approximately 15%. The data quality factor thresholds of 0.8 and 0.5 are based on the following: 0.8 corresponds to the empirical boundary between high signal-to-noise ratio and low distortion; field statistics show that above this value, signal distortion is low and it can be used for high-resolution analysis; 0.5 corresponds to a moderate level of signal-to-noise ratio and distortion. Below this value, the short time window needs to be shortened to ensure time resolution.

[0065] The constructed short-time cepstral maps were uniformly scaled to 300 pixels on the time axis and 500 pixels on the frequency axis before being input into the classification model. The scaling was performed using bilinear interpolation to preserve texture features.

[0066] Step 2.5: Feature separation of classification model.

[0067] The constructed short-time cepstral image is input into a lightweight convolutional neural network classification model. The network structure of this classification model is as follows: the input layer receives a two-dimensional image with a time axis of 300 points and a frequency axis of 500 points; the first convolutional layer uses a 3×3 convolutional kernel and outputs 16 feature maps; the second convolutional layer uses a 5×5 convolutional kernel and outputs 32 feature maps; the first pooling layer uses 2×2 max pooling; the third convolutional layer uses a 3×3 convolutional kernel and outputs 64 feature maps; the second pooling layer is halved again; these are flattened into feature vectors; the first fully connected layer outputs 128 nodes; the second fully connected layer outputs 2 nodes, corresponding to crack reflections and harmonic interference, respectively, which are normalized by the Softmax activation function and then output as classification probabilities.

[0068] The convolutional neural network was trained on a large dataset of cepstral time-frequency maps of fracturing conditions. In the training dataset, each two-dimensional time-frequency map was manually labeled to indicate the continuous blocky texture regions belonging to fracture reflection and the discrete spot-like texture regions belonging to perforation cluster harmonic interference. During training, cross-entropy was used as the loss function, and an adaptive moment estimation optimizer was used to update the network weights.

[0069] After the convolutional neural network is trained, during forward inference, a short-time cepstral graph is input, and the crack reflection probability, harmonic interference probability, and the corresponding feature's position coordinates in the graph are output. If the crack reflection probability is greater than 0.7, the feature is determined to be a true crack reflection. If the crack reflection probability is less than or equal to 0.7, the feature is determined not to meet the confidence requirement for true crack reflection and is marked as suspected interference or weak reflection. It is not included in the subsequent crack feature set, but it can still participate in trend verification in the multi-sensor spatial collaborative constraint in step 2.7. Finally, it is not output as a valid crack. The reflection probability threshold of 0.7 is determined based on the balance point of the accuracy-recall curve on the training set. Setting it to 0.7 can obtain better overall performance.

[0070] Step 2.6: Posterior verification of the physical model of acoustic wave attenuation.

[0071] The output of the classification model needs to be verified by the physical model of acoustic attenuation. When acoustic waves propagate in downhole fluids, their energy decreases exponentially with increasing propagation distance. For typical water-based fracturing fluids, the acoustic attenuation coefficient is 0.005 dB / m, for oil-based fracturing fluids it is 0.003 dB / m, and for oil-water mixed fracturing fluids it is 0.004 dB / m. All coefficients are obtained from laboratory tests and field calibration.

[0072] Based on the center position of the crack reflection characteristics identified by the classification model on the two-way propagation time axis of the water hammer wave, the theoretical propagation distance of the signal is deduced. The theoretical propagation distance equals the two-way propagation time of the water hammer wave multiplied by the reference wave velocity and then divided by 2. Substituting the theoretical propagation distance into the acoustic wave attenuation physical model, the theoretical energy attenuation is calculated. Theoretical energy attenuation = initial energy × 10 (﹣声波衰减系数×理论传播距离 / 10) The initial energy is taken as the actual amplitude of this feature in the short-time cepstrum.

[0073] Compare the theoretical energy attenuation with the actual signal energy detected by the network: Relative deviation = |Actual detected signal energy - Theoretical energy attenuation| ÷ Theoretical energy attenuation × 100%; If the relative deviation > 20%, the network output confidence is determined to be too low, and the confidence of this feature is automatically multiplied by an attenuation factor of 0.5; at the same time, secondary verification is triggered: cross-correlation analysis is performed between the cepstrum profile of this period and the cepstrum profile of the adjacent period, and the peak value of the cross-correlation function is calculated; if the cross-correlation peak value < 0.6, the confidence of this feature is further multiplied by a secondary attenuation factor of 0.3; if the cross-correlation peak value ≥ 0.6, the secondary attenuation factor is not multiplied; the threshold of 0.6 comes from the statistical distribution of the cross-correlation peak value between the actual fracture reflection and the profile of the adjacent period in the measured data of 50 wells, and values ​​below 0.6 are considered outliers; if the relative deviation ≤ 20%, the network output confidence is determined to be acceptable, and the original confidence is maintained without attenuation; the above threshold and attenuation factor are determined based on the statistical results of the error distribution of the measured data of multiple wells.

[0074] Step 2.7: Time consistency verification and multi-sensor spatial coordination constraints.

[0075] After the posterior verification, a time consistency verification is further performed; for three consecutive periodic windows, the number of times the classification model outputs crack reflections is counted; if a certain reflection feature appears only in one period and does not appear in the two adjacent periods, then the feature is determined to be a periodic occasional interference, and its posterior probability is automatically multiplied by a time consistency decay factor of 0.3.

[0076] For sensor nodes deployed at different depths along the wellbore, the spatial propagation law of water hammer waves from the same source is used for multi-sensor spatial collaborative constraints. In the cepstral profiles of the same water hammer wave collected by sensor nodes at different depths, the two-way propagation time of the water hammer wave reflected from the actual fracture should follow these rules: the closer the sensor node is to the fracture location, the smaller the two-way propagation time of the water hammer wave at the fracture reflection peak in its cepstral profile, because the propagation path is shorter and the round-trip time is shorter; the farther the sensor node is, the larger the two-way propagation time of the water hammer wave. The identification results of multiple sensor nodes are sorted according to their depth positions to examine the two-way propagation time of the water hammer wave at the fracture reflection peak. Whether the length value exhibits a V-shaped distribution pattern with the crack depth as the minimum value; if all sensors are located on the same side of the crack, it should be a monotonic change proportional to the distance; if the identification result of a certain node deviates significantly from this pattern while other nodes are consistent, the system automatically reduces the confidence level of the identification result of that node, with a reduction coefficient of 0.5, which is an empirical value determined based on the linear relationship between the degree of deviation and the confidence level reduction in Monte Carlo simulation; the water hammer wave two-way propagation time axis of the node is linearly interpolated and corrected according to the water hammer wave two-way propagation time values ​​of two adjacent nodes: the corrected water hammer wave two-way propagation time value is equal to the average of the water hammer wave two-way propagation time values ​​of adjacent nodes.

[0077] The deviation is defined as the difference between the sorted position of the two-way propagation time of the water hammer wave and the sorted position of the depth at that node is greater than 2. This threshold of 2 is based on Monte Carlo simulation: the probability that random noise causes the sorted position of the two-way propagation time of the water hammer wave to deviate from the sorted position of the depth by more than 2 positions is less than 5%, so it is used as an effective discrimination threshold.

[0078] Step 2.8: In-depth mining of wave velocity anomaly information, inversion of crack opening width and bubble volume fraction.

[0079] The output spatial wave velocity field profile is further analyzed to quantitatively convert wave velocity anomalies into key engineering parameters of underground fractures and fluid conditions.

[0080] Step 2.8.1: Non-background processing and anomaly separation.

[0081] First, the background wave velocity drift caused by temperature changes and hydrostatic pressure gradients is eliminated. Based on real-time data from downhole temperature sensing units and pressure sensors, the theoretical background wave velocity under pure water conditions at this depth is dynamically calculated according to fluid property models such as the Wood equation or empirical formulas. The residual wave velocity anomaly is the measured local wave velocity minus the theoretical background wave velocity. When the absolute value of the residual wave velocity anomaly exceeds 1% of the reference wave velocity, it is determined that there is a physical anomaly at this depth, triggering subsequent quantitative inversion. The 1% threshold comes from the superposition of measurement errors of typical sensors and natural fluctuations.

[0082] Step 2.8.2: Inversion of bubble volume fraction.

[0083] A physical model of fluid bubble volume fraction is used, which is based on the Wood equation in classical acoustic wave theory and establishes a quantitative relationship between bubble volume fraction and water hammer wave velocity. The Wood equation is a nonlinear equation, requiring the bubble volume fraction to be solved based on measured wave velocities. Parameters such as measured local wave velocity, pure fluid wave velocity, gas wave velocity, fluid density, and gas density are substituted into the Wood equation, and numerical methods, such as the bisection method or Newton's iteration method, are used to solve for the bubble volume fraction. Initially, the Wood equation is used to obtain an initial value of the bubble volume fraction through nonlinear solution. If this initial value is less than 2%, a linear approximation formula can be used subsequently. Rapid calculation: If the initial value is ≥2%, nonlinear solution results must be used, and linear approximation is not allowed; the bubble volume fraction is approximately equal to the wave velocity decrease percentage divided by the conversion factor; wave velocity decrease percentage = (reference wave velocity - measured local wave velocity) / reference wave velocity; the conversion factor is differentiated according to the fracturing fluid type: 25~30 for water-based fracturing fluid, 20~24 for oil-based fracturing fluid, and 22~26 for oil-water mixture; when outputting the results, the inversion method used and the error range that may be introduced by the linear approximation must be indicated, and the error range shall not exceed ±10%; output the bubble volume fraction distribution profile along the wellbore depth.

[0084] Step 2.8.3: Crack opening width inversion.

[0085] A quantitative mapping model between fracture opening width and wave velocity was established. The core principle is that fracture opening leads to discontinuity in the propagation medium, causing wavefront diffraction and energy attenuation when elastic waves pass through, manifested as a local decrease in wave velocity. Through theoretical derivation and numerical simulation calibration, such as based on a displacement discontinuity model, the mapping relationship was established: the opening width equals the absolute value of the residual wave velocity anomaly multiplied by a calibration coefficient. The calibration coefficient depends on the acoustic properties of the rock and the fracture morphology. Rock acoustic properties include P-wave velocity and density. The calibration coefficient is derived based on the theory of elastic wave propagation in fractured media. The fracture opening width calibration coefficient is 0.01~0.03 mm / (m / s) for hard strata and 0.04~0.1 mm / (m / s) for soft strata. mm / (m / s); Simultaneously, using the crack reflection feature location output in step 2.5, the cepstral amplitude intensity corresponding to that location is extracted from the short-time cepstral plot, and the crack reflection coefficient is calculated in reverse using the acoustic attenuation model. Then, another opening width estimate is obtained through the theoretical relationship between the reflection coefficient and the opening width, which is similar to the relationship between crack impedance and opening size in elastic wave theory. Finally, the two estimates are fused using a weighted average method to obtain a more reliable final value for the opening width. Among them, the weight of the residual wave velocity anomaly estimate is taken as the data quality factor of spatial interpolation at that depth, and the weight of the reflection coefficient estimate is taken as the confidence level of the crack reflection feature output in step 2.5.

[0086] Step 2.8.4: Multiphysics consistency verification and risk warning.

[0087] The inverted bubble volume fraction and fracture opening width are spatially aligned within the same depth profile for consistency inference and verification. If both an increase in bubble volume fraction and an increase in crack opening width occur at a certain depth, the physical mechanism of pressure drop degassing caused by crack propagation is verified, increasing the confidence level of both. If the two show opposite trends in two or more consecutive adjacent depth points, a dual verification of wave velocity and stress is triggered. If the verification passes, i.e., the wave velocity mutation zone coincides with the stress shielding layer and the pressure transient characteristics are normal, only the physical mechanism is marked as pending confirmation, and no warning is triggered. If the verification fails, it is determined to be an abnormal working condition, such as abnormal fracturing fluid performance or fluid channeling, and a first-level warning is output, and a comprehensive diagnosis is performed in conjunction with the evidence conflict handling mechanism in step 4.3. When the fracture opening width exceeds the set safety threshold, it is determined that there is a risk of excessive fracture propagation at this location, providing a quantitative basis for subsequent construction optimization. The safety threshold is 2.5 mm, and its determination is based on the linear elastic fracture mechanics model in rock mechanics. The critical condition for excessive fracture propagation is that the net pressure inside the fracture exceeds the difference between the minimum horizontal principal stress of the interlayer and the minimum horizontal principal stress of the reservoir. By substituting this stress difference, interlayer thickness, rock Young's modulus, and Poisson's ratio into the stress difference-opening width relationship, the critical opening width when the fracture penetrates the interlayer can be estimated.

[0088] It should be noted that all the conversion coefficients, calibration coefficients and thresholds involved in the above inversion are exemplary preferred values, and those skilled in the art can adjust them according to specific circumstances.

[0089] Step 2.9: Multi-cycle adaptive weighted accumulation and generation of liquid inflow intensity profile.

[0090] After completing all the above checks and constraints, weighting coefficients are dynamically assigned to the cepstral amplitudes of each period in the final two-way propagation time-axis aligned profile of the water hammer wave, based on the data quality factor of each period window and the signal-to-noise ratio within that period. The weighting coefficient is calculated as follows: the weighting coefficient for that period is equal to the data quality factor of that period multiplied by the normalized signal-to-noise ratio of that period, and then divided by the sum of the weighting coefficients of all periods. Wherein, the normalized signal-to-noise ratio of that period is equal to the signal-to-noise ratio of that period divided by the maximum value of the signal-to-noise ratio of all periods. The cepstral amplitudes in the final two-way propagation time-axis aligned profile of the water hammer wave for each period are weighted and summed according to the above weighting coefficients to obtain the inflow intensity profile. The value corresponding to each depth position in the inflow intensity profile represents the relative inflow ratio at that position, and the larger the value, the greater the inflow.

[0091] Meanwhile, the spatial wave velocity field profile, bubble volume fraction profile, and crack opening width profile obtained in step 2.3 are used as additional outputs for cross-validation with the stress shielding layer positioning results in the subsequent expert inference layer.

[0092] Step 3: Layered stress-driven mechanical model structural deviation detection and adaptive switching inversion system: The mechanical model includes the PKN (Perkins-Kern-Norden) model, the KGD (Khristianovic-Geertsma-deKlerk) model, and the P3D (pseudo-3D fracture propagation) model; the fracture height is iteratively calculated using the fluid inflow intensity profile, fracture reflection characteristic location, spatial wave velocity field profile, bubble volume fraction profile, and fracture opening width profile output in Step 2, as well as the layered stress data and high-frequency pressure transient characteristics collected in Step 1 as inputs.

[0093] Step 3.1: Constructing a layered stress sensor network to define stress profiles and locate stress shielding layers.

[0094] The micro-stress sensing unit deploys a measuring point every 10 meters along the wellbore depth. Each measuring point contains three mutually orthogonal micro-stress sensing elements, which measure the vertical stress, maximum horizontal principal stress, and minimum horizontal principal stress, respectively. A longitudinal stress profile curve is constructed using these data, and the measuring points are connected using cubic spline interpolation. The stress profile curve is traversed along the depth direction, and the stress gradient between adjacent depth points is calculated in real time. When the stress gradient of three or more consecutive measuring points within a certain depth interval is >0.1 MPa / m, it is determined that a stress-shielding layer or a weak stress interlayer exists in that interval, and its upper and lower boundaries and average stress gradient are recorded. If the stress gradient is ≤0.1 MPa / m, the interval is determined to be a normal stress continuity layer and is not considered a shielding layer. This gradient threshold of 0.1 MPa / m is an empirical value set based on the statistical values ​​of stress differences in typical shale reservoirs.

[0095] Step 3.2: Iterative initialization and parameter assimilation of crack height.

[0096] The initial value of the fracture height is half the reservoir thickness, and the PKN model is selected as the default initial model. In each iteration, the ensemble Kalman filter algorithm is used for parameter assimilation. First, the state vector of the mechanical model is established, which includes the rock Young's modulus, fracturing fluid filtration efficiency, net pressure in the fracture, net pressure distribution in the fracture, and filtration coefficient distribution. The fracture geometric parameters are used as the post-processed output after filtering and do not participate in the assimilation loop. A set of 100 samples is generated for each parameter, and the initial value of each sample is randomly perturbed within a reasonable range: the perturbation range of the rock Young's modulus is ±10%, and the perturbation range of the fracturing fluid filtration efficiency is ±20%. The sample size of 100 and the perturbation range are commonly used settings in engineering practice, which can ensure statistical stability without excessively increasing the computational burden.

[0097] The transient pressure features at the moment of microfracture opening were extracted from high-frequency pressure time-series data, including: pressure rise steepness, energy surge increment in the 300Hz to 500Hz frequency band, and acoustic emission event count rate detected by vibration sensors; the pressure rise steepness is the pressure rise amplitude per millisecond, the energy surge increment is obtained by calculating the signal energy after bandpass filtering of this frequency band, and the acoustic emission event count rate is the number of events per second with vibration peak values ​​exceeding 0.05g, where g is the gravitational acceleration. This vibration threshold of 0.05g is an industry-standard way to distinguish acoustic emission from formation fractures from wellbore environmental interference; these features were used as observation vectors.

[0098] The calculation steps of ensemble Kalman filtering are as follows: calculate the ensemble mean and ensemble covariance of the predicted fracture volume for each sample; calculate the ensemble mean and ensemble covariance of the observation vectors; calculate the Kalman gain matrix, which is equal to the state prediction covariance multiplied by the inverse of the observation prediction covariance; multiply the Kalman gain by the deviation between the observed and predicted values ​​to obtain the state correction; add the correction to the state vector of each sample to obtain the updated parameter set; take the mean of the updated set as the parameter estimate for this iteration; through the above process, the numerical parameters such as the rock Young's modulus and fracturing fluid filtration efficiency in the mechanical model are updated online and dynamically.

[0099] Step 3.3: Model structure deviation detection and wave velocity stress dual constraint verification.

[0100] After parameter assimilation in each iteration, a model structure deviation detection is performed. The upper and lower boundaries of the fracture height calculated in the current iteration are compared with the stress shielding layer boundary located in step 3.1. If the following two conditions are met simultaneously, it is determined that the structural assumptions of the current mechanical model have fundamentally deviated from the actual fracture propagation morphology: the calculated upper boundary of the fracture continuously exceeds the measured upper boundary depth of the stress shielding layer in three consecutive iterations, with an excess of >0.5m; the relative error between the fracture volume predicted by the current mechanical model and the target volume obtained by multiplying the actual injected fluid volume by the fluid efficiency is still >5%. The excess threshold of 0.5m and the relative error threshold of 5% are both set based on the statistical distribution of fitting errors from multiple wells in the field.

[0101] At this point, the spatial wave velocity field profile output in step 2.3 is further used for double verification: depth overlap is defined as an absolute difference between the depth of the center of the wave velocity abrupt change zone and the depth of the stress shielding layer boundary being less than 2m; significant offset is defined as an absolute difference greater than 5m; and partial overlap is defined as being between 2m and 5m, with the following processing rules: If the location of the wave velocity abrupt change zone coincides with the depth of the stress shielding layer boundary, the confidence level of the model deviation is increased; if the location of the wave velocity abrupt change zone is significantly offset from the boundary of the stress shielding layer, the wave velocity re-estimation in step 1.2 is triggered, and the location of the stress shielding layer may have an error; when the wave velocity re-estimation in step 1.2 is triggered, the depth range of the current wave velocity anomaly zone is used as prior information, and a local constraint of this depth range is added to the physical verification model in step 1.2, that is, the wave velocity estimate within this range is forced to converge to the current inversion value, and the wave velocity estimation process is re-executed; if there is a partial overlap, the original confidence level is maintained, no additional operation is triggered, and the iteration continues.

[0102] The error threshold remains unchanged after model switching: After switching to KGD or P3D model, the volume error threshold remains 5%, and the height excess threshold remains 0.5 meters; if the crack height exceeds the stress shielding layer boundary but the volume error meets the threshold, i.e., pseudo-convergence, then switch the mechanical model again first, trying in the order of PKN→KGD→P3D; if all models have been tried and pseudo-convergence is still found, then force a re-iteration and relax the volume error threshold to 10%; if convergence is still not achieved after reaching the maximum number of iterations of 50, then output the current optimal solution and mark the convergence anomaly, and link it with the bridge plug seal evaluation in step 4: if the bridge plug seal status is likely to leak or the data is insufficient, then attribute the anomaly to the data distortion caused by seal failure and indicate it in the diagnostic report.

[0103] Step 3.4: Stress boundary deviation drives adaptive switching of the model.

[0104] After detecting structural deviations, the mechanical model is automatically switched based on the parameters of the stress shielding layer. Decision rule 1: If the average stress gradient of the stress shielding layer is >0.2MPa / m and the thickness of the shielding layer is >5m, the system automatically switches to the KGD mechanical model. Decision rule 2: If the calculated crack height simultaneously crosses two or more stress shielding layers of different depth ranges, the system automatically switches to the P3D mechanical model. After the switch is completed, the system takes the crack geometry and hydrodynamic state at the time of the switch as the starting point for the new iteration and continues to iterate under the new mechanical model framework. The crack geometry includes the current crack half-length and average crack width, and the hydrodynamic state includes the net pressure distribution within the crack.

[0105] After the model adaptively switches, the volume error threshold and height excess threshold used for convergence determination remain unchanged.

[0106] Step 3.5: Double convergence constraint and iteration termination.

[0107] A dual convergence constraint is introduced; the iteration termination condition requires that the following conditions be met simultaneously: Condition 1: Volume error < dynamic threshold; the volume error is equal to the absolute value of the difference between the fracture volume predicted by the mechanical model and the target volume obtained by multiplying the actual injected fluid volume by the fluid efficiency, divided by the target volume; the dynamic threshold is taken as 5% of the current fracture volume and 0.5m. 3 The larger of the two, 0.5m 3 The absolute lower limit of volume error is derived from 1% of the typical value of single-stage injection volume in small-scale fracturing tests in the field, avoiding overly strict relative error thresholds for small-volume fractures; Condition 2: Consistency of top and bottom stress boundaries; that is, the upper and lower limits of the currently calculated fracture height do not exceed the boundary range of the stress shielding layer measured by the intelligent stress sensor; if the stress shielding layer does not exist, this condition is automatically met; if the calculated height exceeds the boundary, the excess must monotonically decrease in five consecutive iterations and eventually converge to less than 0.1m.

[0108] If the calculated crack height exceeds the position of the stress shielding layer measured by the sensor, but the volume error still meets condition one, it is judged as pseudo-convergence, and a forced re-iteration or model switching is triggered again; if none of the above conditions are met, the iteration continues until the maximum number of iterations is 50. If convergence is not achieved after reaching the maximum number of iterations, the current optimal solution is output and marked as not strictly converged.

[0109] Step 3.6: Early warning of crack propagation based on physical constraints of bubble volume fraction and crack opening width.

[0110] Read the bubble volume fraction distribution profile and crack opening width distribution profile obtained from step 2.8.

[0111] If, in the region outside the upper and lower boundaries of the fracture height in the current iteration, the inverted bubble volume fraction exceeds the first preset threshold (e.g., 0.1%), and the fracture opening width exceeds the second preset threshold (e.g., opening width > 1 mm), it is first determined whether the currently calculated fracture height anomaly region is still located inside the reservoir, i.e., between the upper and lower interlayers. If the anomaly region is located inside the reservoir, the fracture height search range is expanded in the next iteration, extending the search range towards the anomaly region by 50% of the difference between the depth of the anomaly region and the boundary. If the anomaly region exceeds the reservoir interlayer range, it is directly determined that there is a risk of excessive fracture propagation and an early warning is output, without further expanding the search range.

[0112] If, within the fracture height range of this iteration, the fracture opening width exhibits a non-monotonic change along the wellbore depth direction, i.e., it is not the typical pattern where the fracture opening is largest at the center and decreases towards both sides, then it is determined that there is a false convergence situation where the fracture height is too large but the volume fitting still matches, triggering the wave velocity and stress dual constraint verification in step 3.3 or the false convergence detection in step 3.5.

[0113] The above thresholds are all based on the statistical distribution of opening width of typical hydraulic fractures and the critical value of fracturing fluid bubble precipitation. In actual application, they can be adjusted according to specific reservoir conditions. This step enables the mechanical model inversion to obtain an independent source of constraint information that does not depend on the pressure response, effectively suppressing the ambiguity in mechanical inversion.

[0114] After the iteration terminates, the final parameters such as crack height, average crack width, and total crack volume are output.

[0115] Step 4: Dynamic basic probability allocation correction and evidence conflict handling based on the self-diagnostic quality factor of distributed multimodal intelligent sensors: Integrate the three types of analysis results of pressure waveform analysis, cepstral peak characteristics, and multi-cycle signal stability, and combine them with DS evidence theory to complete the bridge plug sealing status assessment.

[0116] Step 4.1: Definition of the evidence body for bridge plug seal assessment and initial basic probability assignment.

[0117] Four types of bridge plug seal assessment evidence bodies are defined: pressure waveform analysis evidence body, cepstral peak feature evidence body, multi-cycle stability evidence body, and temperature and vibration sensor data evidence body. The initial basic probability allocation for each evidence body is assigned probability mass to three sealing state assumptions: good seal, seal failure, and uncertainty, with a sum of 1. The preset values ​​are as follows: pressure waveform analysis evidence body: good seal 0.6, seal failure 0.3, uncertainty 0.1; cepstral peak feature evidence body: good seal 0.7, seal failure 0.2, uncertainty 0.1; multi-cycle stability evidence body: good seal 0.5, seal failure 0.3, uncertainty 0.2; temperature and vibration sensor data evidence body: good seal 0.6, seal failure 0.3, uncertainty 0.1. These initial values ​​are derived from the statistical distribution of bridge plug seal verification results after fracturing in multiple wells.

[0118] Step 4.2: Dynamic basic probability allocation correction based on data quality factor.

[0119] The basic probability allocation of each piece of evidence is corrected online using the data quality factor transmitted in real time by distributed multimodal intelligent sensors. If the data quality factor > 0.8, the probability value of the proposition with the largest initial basic probability allocation in that piece of evidence (i.e., the state most favored by that piece of evidence) is increased by an increase of (data quality factor - 0.8) × increase coefficient 0.2. The probability values ​​of other propositions are reduced proportionally to ensure that the sum of the total probabilities of that piece of evidence remains 1, and then normalization is performed. The increase coefficient 0.2 is determined based on grid search optimization using a training set of 30 wells. If 0.4 ≤ data quality factor ≤ 0.8, no correction is made, and the initial basic probability allocation is maintained. If 0.2 ≤ data quality factor < 0.4, the probability is reduced linearly, and the initial probability of all propositions is adjusted accordingly. The initial basic probability allocation is multiplied by the data quality factor and divided by 0.4, then renormalized. If the data quality factor is less than 0.2, the evidence is automatically deemed unusable, and all its basic probability allocations are transferred to an uncertain proposition. The high factor threshold of 0.8, the medium factor threshold of 0.4, and the low factor threshold of 0.2 are determined based on the cross-validation results of the data quality factor and diagnostic accuracy. 0.8 is the high confidence threshold, 0.4 is the low confidence threshold, and 0.2 is the unusable threshold. Below 0.2, the signal-to-noise ratio is extremely low or the distortion is severe, failing to provide effective information. If the signal distortion index is greater than 20%, regardless of the calculated value of the data quality factor, it is forcibly set to the low confidence lower limit of 0.2, and then the correction is performed according to the above rules.

[0120] Step 4.3: Conflict of evidence detection and redistribution.

[0121] Calculate the conflict coefficient in the DS evidence synthesis rule; for two pieces of evidence, the conflict coefficient is defined as the sum of the probability products of mutually exclusive propositions between the pieces of evidence, including well sealed, unsealed, and uncertain; uncertain propositions are not included in the calculation of the conflict coefficient; the conflict coefficient ranges from 0 to 1, and the larger the conflict coefficient, the more significant the conflict between the pieces of evidence.

[0122] If the conflict coefficient is >0.5, it is determined that there is a significant conflict among multiple pieces of evidence. In this case, a weighted correction method can be used to preprocess the evidence before synthesis, or the confidence level of the conclusion can be adjusted according to the degree of conflict after synthesis according to the standard DS rules. At the same time, the rule is added: evidence with a data quality factor lower than 0.4 is directly removed and no longer participates in the synthesis. If the conflict coefficient is ≤0.5, it is determined that the consistency between the evidence is good, and it is directly fused according to the DS synthesis rules. The conflict threshold of 0.5 comes from the commonly used definition of high conflict in DS theory, that is, when the conflict coefficient exceeds 0.5, the contradiction between the evidence is significant and intervention is required.

[0123] Step 4.4: Tiered early warning and uncertainty visualization.

[0124] After dynamic correction and conflict resolution, the evidence is fused according to the DS evidence theory synthesis rules to obtain the posterior probabilities of three propositions: well sealed, poorly sealed, and uncertain. The following four-level decision conclusions are output according to the criteria below: Well sealed: Probability of well sealed ≥ 0.75 and probability of uncertainty < 0.1; Suspected minor leakage: 0.5 < probability of well sealed < 0.75, or 0.1 ≤ probability of uncertainty ≤ 0.3; High probability of leakage: 0.25 ≤ probability of well sealed ≤ 0.5, and probability of poorly sealed > 0.3. If data is insufficient, retesting is required: the probability of good sealing is <0.25, or there is evidence of a data quality factor <0.2, or the probability of uncertainty is >0.3; all conclusions are accompanied by uncertainty boundaries, indicated as a percentage.

[0125] The threshold values ​​mentioned above are based on the following criteria: 0.75 is derived from the high confidence threshold in typical engineering decision-making, corresponding to a reliability of over 85%; 0.5 is the probability midpoint, representing the critical value of uncertainty; 0.25 is the low confidence threshold, below which evidence is considered insufficient; 0.3 is the maximum acceptable uncertainty ratio, derived from the empirical value that the false alarm rate should not exceed 30% in multi-well statistics; and 0.2 is the data unavailability threshold, consistent with step 4.2.

[0126] Step 4.5: G-function closure pressure analysis.

[0127] Based on the pressure drop curve after pump shutdown, the G-function graphical method is used to calculate the closure pressure and filtration coefficient. Specifically, the method involves: taking the pump shutdown moment as the starting point and extracting a data sequence where the pressure tends to stabilize; converting the actual time into a dimensionless G-function, where the G-function value depends on the ratio of the pump shutdown time to the total pumping time; plotting the pressure versus G-function curve and the pressure derivative versus G-function curve; analyzing the pressure derivative versus G-function curve within the G-function range where the G-function value is greater than 1.5 to identify a stable linear decreasing segment; performing linear regression on this segment and extrapolating the regression line to intersect the axis where the pressure derivative equals 0, denoted as Gc at this intersection point; finally, finding the pressure value corresponding to Gc on the pressure versus G-function curve, which is the fracture closure pressure; simultaneously, calculating the filtration coefficient based on the slope of this linear segment, construction parameters, and the fracture height output in step 3; and outputting the closure pressure and filtration coefficient; construction parameters include displacement, fracturing fluid viscosity, and fluid compressibility.

[0128] Step 5: Cross-module expert system integrated reasoning and result output: The crack location results, liquid ingress intensity profile, spatial wave velocity field profile, bubble volume fraction profile, and crack opening width profile output in Step 2, the crack height and total crack volume output in Step 3, the bridge plug sealing probability classification output in Step 4.4, and the closing pressure value and filtration coefficient output in Step 4.5 are all converted into fact objects with confidence.

[0129] Load the pre-built IF-THEN cross-module inference rule base, which contains the following typical rules. The threshold in each rule is an optimal value obtained based on field data statistics or theoretical analysis: Rule 1: If the posterior probability of a good bridge plug seal is >0.8 and the ratio of the G-function closing pressure to the instantaneous pump stop pressure is <0.8, then it is determined that there is a micro-leakage in the bridge plug. The confidence level is equal to 0.7 multiplied by the minimum confidence level of the ratio of the posterior probability of a good bridge plug seal to the G-function ratio. The confidence level of the G-function ratio is the measurement reliability of this ratio, which is determined by the pressure sensor accuracy and data quality factor, and is set to 0.9 by default. The instantaneous pump stop pressure refers to the pressure value recorded at the moment of pump stop, which is the average value of the pressure sensor readings within 0.1 seconds before and after that moment. Rule 2: If the fracture height exceeds 1.5 times the reservoir thickness and the coefficient of variation of the fluid ingress intensity profile is >0.3, then the fracture is determined to have broken through the stress shielding layer, and the confidence level is equal to 0.85 multiplied by the minimum value of the fracture height confidence level and the coefficient of variation confidence level. Rule 3: If the bridge plug seal assessment result indicates insufficient data requiring retesting and the volume error is >10% when the crack height converges iteratively, it is recommended to retest the bridge plug and re-collect pressure data, with a confidence level of 0.9. Rule 4: If the crack reflection probability is greater than 0.8 and the matching error between the position of the main inflow peak and the depth of the perforation cluster in the inflow intensity profile is less than 5m, then the perforation cluster is deemed to be effectively opened, with a confidence level equal to 0.9 multiplied by the crack reflection probability. Rule 5: If the width of the deep crack opening obtained by inversion in step 2.8 is greater than 1 mm, and the crack reflection probability at that location output in step 2.5 is greater than 0.6, then the main crack of the perforation cluster at that location is determined to be effectively opened, and the confidence level is the higher of the two multiplied by 0.9. Rule 6: If the crack reflection probability is greater than 0.7 and the crack opening width obtained from the inversion in step 2.8 is less than 0.2 mm, then it is judged to be high reflection but low opening, suspected to be microcrack development rather than main crack, and the confidence level is the crack reflection probability multiplied by 0.7. Rule 7: If the fracture opening width obtained by inversion in step 2.8 is greater than 1.5 mm and the fracture height exceeds the reservoir thickness when the iteration in step 3.5 ends, it is determined that the fracture height has been excessively extended and there is a risk of perforation of the interlayer. The confidence level is the minimum of the two confidence levels multiplied by 0.85. Rule 8: If the volume fraction of bubbles obtained from the inversion in step 2.8 is greater than 0.001 at three consecutive depth points, it is determined that the fracturing fluid is severely degassed. It is recommended to adjust the fracturing fluid formula or reduce the discharge rate. The confidence level is taken as 0.9 of the average confidence level at each point. Rule 9: If there is a wave velocity abrupt change zone in the spatial wave velocity field profile and the overlap of the abrupt change zone with the boundary depth of the stress shielding layer is >80%, then the wave velocity anomaly and the stress shielding are determined to be of the same origin, and the crack propagation is controlled by the stress interface, with a confidence level of 0.8.

[0130] The forward chain reasoning engine is started; the reasoning process is as follows: all initial facts are stored in the working memory; each rule in the rule base is scanned sequentially, and the condition part of the rule is matched with the facts in the working memory; if the match is successful, the conclusion part of the rule is triggered, and the new facts and their propagated confidence are added to the working memory; the confidence propagation formula is: the confidence of the new fact is equal to the rule confidence multiplied by the arithmetic mean of the conditional fact confidence, and then multiplied by an evidence accumulation enhancement factor; this enhancement factor is not less than 1 and has a maximum of 1.5; at the same time, the lower limit of confidence is set to 0.05, and facts below the lower limit are not included in the reasoning; the above scanning and matching process is repeated until no new facts are generated in a continuous round of scanning.

[0131] After reasoning is completed, a comprehensive diagnostic conclusion is output; the diagnostic conclusion includes: Fracture effectiveness evaluation: number of fractures, fluid ingress ratio of each cluster, and ratio of fracture height to reservoir thickness; Bridge plug seal comprehensive risk level: one of the following: good seal, suspected minor leakage, high probability of leakage, insufficient data; Comprehensive score for fracturing effect: expressed as a percentage, equal to the posterior probability of good sealing multiplied by the normalized value of the ratio of fracture height to reservoir thickness multiplied by 100; normalization rules: the upper limit of normalization is set at 2 times the reservoir thickness; when the ratio of fracture height to reservoir thickness is ≥2, the normalized value is 1.0; when the ratio is <2, the normalized value = actual ratio ÷ 2; comprehensive score for fracturing effect = posterior probability of good sealing × normalized value × 100; Subsequent construction optimization suggestions: such as whether additional perforations are needed, whether secondary setting of bridge plugs is required, and whether the fracturing fluid formula needs to be adjusted; Key physical parameters appendix: provided in text list format, including depth location, spatial wave velocity, residual wave velocity anomaly, inverted bubble volume fraction, inverted fracture opening width, and confidence percentage.

[0132] All conclusions are accompanied by confidence level labels; at the same time, a traceable analysis report is generated, which records the rules and original data sources on which each step of reasoning was based.

[0133] The entire hydraulic fracturing monitoring workflow is completed; key parameters, such as anomalies in the spatial wave velocity field profile and the local wave velocity correction values ​​obtained by inversion, are transmitted back to the front-end edge computing unit in a non-real-time manner for the initialization of preprocessing parameters for the next well or the next fracturing stage, and are not adjusted online adaptively within the same acquisition cycle.

[0134] Example 2: Please see Figure 4 Based on Example 1, this example also provides a hydraulic fracturing monitoring system for high-frequency pressure water hammer wave and cepstral comprehensive analysis, including: a data acquisition and wave velocity fusion verification module, which acquires high-frequency pressure, ground stress, temperature and vibration data of the wellbore, performs adaptive baseline drift correction, and obtains the reference wave velocity, data quality factor, pure pressure sequence and water hammer wave two-way propagation time through multi-source fusion estimation and physical verification; The physical constraint cepstral enhancement and feature inversion module performs cepstral transformation and dynamic time warping alignment based on the pure pressure sequence and reference wave velocity, inverts the spatial wave velocity field profile, and extracts the crack reflection feature location, bubble volume fraction profile, crack opening width profile and liquid ingress intensity profile by combining neural network classification and acoustic attenuation verification. The mechanical model deviation detection and adaptive inversion module combines geostress data and inversion profiles to locate stress shielding layers, uses filtering algorithms to assimilate parameters, performs wave velocity and stress dual-constraint verification when mechanical model deviation is detected, and adaptively switches mechanical models to output crack height, average crack width and total crack volume. The evidence dynamic correction and conflict fusion module defines multiple types of bridge plug seal evaluation evidence bodies, uses data quality factors to dynamically correct the basic probability allocation, calculates the conflict coefficient between evidence and performs conflict redistribution when the threshold is exceeded, fuses and outputs bridge plug seal status classification, and performs closure pressure analysis based on the pump stop pressure drop curve to obtain closure pressure and filtration coefficient. The cross-level expert rule reasoning and comprehensive diagnosis module transforms multi-dimensional inversion profiles, crack geometric parameters, sealing status classifications, and closure pressure parameters into fact objects with confidence levels, loads a preset reasoning rule base to perform forward chain reasoning, and outputs comprehensive diagnostic conclusions.

[0135] Furthermore, the data acquisition and wave velocity fusion verification module is specifically used to perform five-point window mid-range filtering and local weighted regression to remove trend terms, run pressure cross-correlation wave velocity estimator, vibration propagation wave velocity estimator and theoretical wave velocity estimator in parallel, substitute into the discretized mass conservation and momentum conservation equations to calculate the residuals, select the estimated value with the smallest residual as the reference wave velocity, and if all exceed the limit, call the long short-term memory network wave velocity prediction model to output the predicted wave velocity value.

[0136] The physical constraint cepstrum enhancement and feature inversion module is specifically used to establish forward and reverse dynamic time warping alignment chains with the first period window cepstrum profile as the initial reference, extract the local mapping ratio to calculate the corrected local wave velocity to form a spatial wave velocity field profile, use a lightweight convolutional neural network to perform feature separation on the short-time cepstrum, and perform posterior verification through the acoustic attenuation physical model. Combined with multi-sensor spatial collaborative constraints, it generates the liquid inflow intensity profile and the crack opening width and bubble volume fraction profile.

[0137] The mechanical model deviation detection and adaptive inversion module is specifically used to construct stress profiles and locate stress shielding layers based on data from micro-stress sensing units. It establishes a state vector that includes rock Young's modulus, fracturing fluid filtration efficiency, net pressure within the fracture, fracture half-length, and fracture height to assimilate high-frequency pressure transient characteristics. When the fracture height boundary continuously exceeds the stress shielding layer boundary and the volume error exceeds the limit, structural deviation is determined. Based on the stress shielding layer parameters, it automatically switches to a specific mechanical model to continue iterating until the dual convergence constraints are met.

[0138] The evidence dynamic correction and conflict fusion module is specifically used for pressure waveform analysis, cepstral peak characteristics, multi-cycle stability, and temperature and vibration sensing data evidence. Based on the interval where the data quality factor is located, it performs correction operations such as probability enhancement, maintenance, attenuation, or transfer to uncertain propositions. It calculates the sum of the probability products of mutually exclusive propositions between different evidence bodies as the conflict coefficient. If it exceeds the set threshold, it performs evidence synthesis after preprocessing using a weighted correction method. At the same time, it plots the curves of pressure and pressure derivative as a function of dimensionless G to identify the linear descent segment and extrapolate to obtain the closing pressure.

[0139] The cross-level expert rule reasoning and comprehensive diagnosis module is specifically used to build a reasoning rule library containing multiple preset conditions. It sequentially scans the fact objects in the rule matching working memory. If a match is successful, the conclusion section is triggered, and the new fact and its propagated confidence level are added to the working memory. The confidence level propagation is based on the product of the rule confidence level and the minimum confidence level of each fact in the condition. The scanning is repeated until no new facts are generated and a diagnostic report containing the crack effectiveness evaluation and bridge plug sealing risk level is output.

[0140] The system's workflow is as follows: The data acquisition and wave velocity fusion verification module acquires multimodal data, and outputs reference wave velocity, data quality factor, and pure pressure sequence after correction and fusion verification; The physical constraint cepstral enhancement and feature inversion module performs cepstral transformation and dynamic time warping alignment to invert the spatial wave velocity field, bubble volume fraction, fracture opening width, liquid ingress intensity profile, and fracture reflection characteristics; The mechanical model deviation detection and adaptive inversion module combines geostress and inversion profile, and adaptively iteratively calculates fracture height, fracture width, and total volume through parameter assimilation and structural deviation detection; The evidence dynamic correction and conflict fusion module uses the data quality factor to dynamically correct evidence and handle conflicts, outputs bridge plug sealing status classification, and obtains closure pressure and filtration coefficient through G-function analysis; The cross-level expert rule reasoning and comprehensive diagnosis module transforms all results into confidence facts, and generates diagnostic conclusions such as fracture effectiveness, risk level, and comprehensive score through forward chain reasoning, realizing full automation from data acquisition to intelligent diagnosis.

[0141] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented in software, the above embodiments can be implemented, in whole or in part, as a computer program product. Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution.

[0142] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A hydraulic fracturing monitoring method based on a comprehensive analysis of high-frequency pressure water hammer waves and cepstral spectrum, characterized in that, Includes the following steps: High-frequency pressure, ground stress, temperature and vibration of the wellbore are collected. Adaptive baseline drift correction and multi-source fusion estimation are performed on the high-frequency pressure to obtain the reference wave velocity, data quality factor, pure pressure sequence and water hammer wave two-way propagation time. Based on the reference wave velocity and the two-way propagation time of the water hammer wave, cepstral transformation, dynamic time warping and alignment and feature extraction are performed on the pure pressure sequence to obtain the crack reflection feature location and multiple profiles, including spatial wave velocity field profile, bubble volume fraction profile, crack opening width profile and liquid ingress intensity profile. By combining in-situ stress and multiple profiles, the mechanical model for hydraulic fracturing monitoring is subjected to structural deviation detection and adaptive switching inversion processing to obtain the fracture height, average fracture width and total fracture volume. Define multiple types of bridge plug seal evaluation evidence bodies and assign initial basic probabilities. Use data quality factors to dynamically correct the basic probabilities, calculate the conflict coefficient between evidence bodies and redistribute conflicts when they exceed the conflict threshold, fuse them to obtain the posterior probabilities of good seal, seal failure and uncertainty, output the bridge plug seal status classification, and perform G-function closure pressure analysis based on the pump stop pressure drop curve to obtain the closure pressure and filtration coefficient. The crack reflection characteristics, influent intensity profile, crack height, bridge plug sealing status classification, closure pressure, and filtration coefficient are transformed into fact objects with confidence levels. A pre-built reasoning rule base is loaded to perform forward chain reasoning and output a comprehensive diagnostic conclusion.

2. The hydraulic fracturing monitoring method based on high-frequency pressure water hammer wave and cepstral comprehensive analysis according to claim 1, characterized in that, The multi-source fusion estimation process includes: The pressure cross-correlation wave velocity estimator, vibration propagation wave velocity estimator, and theoretical wave velocity estimator are run in parallel to obtain preliminary wave velocity estimates. The preliminary wave velocity estimates are substituted into the discretized mass and momentum conservation equations for physical verification. The residuals between the predicted and measured pressure values ​​are calculated, and the estimated value with the smallest residual that passes the verification is selected as the reference wave velocity. If all preliminary wave velocity estimates fail the verification, the wave velocity prediction model based on a long short-term memory network is invoked. The model takes the time series data within a fixed time window before pump shutdown as input and outputs the predicted wave velocity value as the reference wave velocity.

3. The hydraulic fracturing monitoring method based on high-frequency pressure water hammer wave and cepstral analysis according to claim 1, characterized in that, The dynamic time warping alignment and feature extraction process includes: Using the cepstral profile of the first periodic window as the initial reference, the cepstral profile of the subsequent periodic windows is positively aligned using a dynamic time warping algorithm to obtain the positive cumulative warping mapping function. Using the periodic window with the highest signal-to-noise ratio as a reference, the reverse alignment is performed through a dynamic time warping algorithm to obtain the reverse cumulative warping mapping function; The local mapping ratio is extracted from the cumulative regularized mapping function. Based on the reference wave velocity and the local mapping ratio, the local wave velocity varying along the well depth is calculated to form a spatial wave velocity field profile. By comparing the forward cumulative regularization mapping function and the reverse cumulative regularization mapping function, if the difference in displacement between the two functions at the same position of the two-way propagation time of the water hammer wave exceeds the preset difference threshold, the wave velocity value in the spatial wave velocity field profile is used as a constraint for weighted correction to obtain the final two-way propagation time axis aligned profile of the water hammer wave.

4. The hydraulic fracturing monitoring method based on high-frequency pressure water hammer wave and cepstral comprehensive analysis according to claim 1, characterized in that, The mechanical models include the PKN mechanical model, the KGD mechanical model, and the P3D mechanical model; The process of structural deviation detection and adaptive switching inversion processing of the mechanical model for hydraulic fracturing monitoring includes: A stress profile is constructed based on the in-situ stress data collected along the wellbore, and the stress gradient is calculated to locate the stress shielding layer. Starting with a preset initial mechanical model and an initial crack height, the mechanical model parameters are dynamically updated by using an ensemble Kalman filter algorithm to assimilate the transient pressure features extracted from high-frequency pressure data. After each iteration, it is checked whether the crack height boundary predicted by the current mechanical model continues to exceed the boundary of the located stress shielding layer, and whether the relative error between the predicted crack volume and the actual target volume exceeds the relative error threshold. If both exceed the relative error threshold, it is determined that the current mechanical model has deviated from the structure. The space wave velocity field profile was used to perform dual-constraint verification of wave velocity stress. If the structural deviation is verified, the system will automatically switch to different crack propagation mechanical models based on the average stress gradient and thickness of the stress shielding layer to continue iterating and obtain the crack height, average crack width and total crack volume.

5. The hydraulic fracturing monitoring method based on high-frequency pressure water hammer wave and cepstral comprehensive analysis according to claim 1, characterized in that, It also includes physical constraint fracture propagation early warning based on bubble volume fraction profile and fracture opening width profile. If, in the region outside the fracture height boundary of the current iteration, the inverted bubble volume fraction exceeds a first preset threshold and the fracture opening width exceeds a second preset threshold, then the region is recorded as an abnormal region, and the depth location of the abnormal region is further determined: if the abnormal region is located inside the reservoir, then in the next iteration, the search range of the fracture height is expanded towards the abnormal region; if the abnormal region exceeds the reservoir interlayer, then it is determined to be a risk of excessive fracture propagation, the search range is no longer expanded, and an early warning is output. If the crack opening width changes non-monotonically along the depth direction within the crack height range of the current iteration, then the wave velocity stress double constraint verification or pseudo-convergence detection is triggered.

6. The hydraulic fracturing monitoring method based on high-frequency pressure water hammer wave and cepstral analysis according to claim 1, characterized in that, The process of dynamically correcting the basic probability using data quality factors includes: setting high factor threshold, medium factor threshold and low factor threshold; If the data quality factor is greater than the high factor threshold, then increase the basic probability assignment value of the most probable state proposition in the evidence body; If the medium factor threshold ≤ data quality factor ≤ high factor threshold, then the initial basic probability allocation remains unchanged; If the low factor threshold ≤ data quality factor < medium factor threshold, then the basic probability allocation values ​​of all propositions in the evidence body are decayed according to a linear relationship. If the data quality factor is less than the low factor threshold, the evidence is deemed unusable, and all its basic probability allocations are transferred to the uncertain proposition. The conflict redistribution process includes: the conflict coefficient is the sum of the probabilities of mutually exclusive propositions among the evidence bodies; uncertain propositions are not included in the conflict coefficient calculation; if the conflict coefficient is greater than the preset conflict threshold, the evidence probability is preprocessed using a weighted correction method before evidence synthesis is performed, and evidence bodies with data quality factors lower than the medium factor threshold are removed.

7. The hydraulic fracturing monitoring method based on high-frequency pressure water hammer wave and cepstral comprehensive analysis according to claim 1, characterized in that, The process of performing G-function closed-loop pressure analysis based on the pump shutdown pressure drop curve includes: Convert the actual time after the pump stops into a dimensionless G function; Plot the curve of the derivative of pressure with respect to the G function as a function of G, and identify a stable linear decreasing segment in the interval where the value of the G function is greater than the preset value; Perform linear regression on the linearly decreasing segment, and extrapolate the regression line until it intersects the axis where the pressure derivative is zero, to obtain the G function value corresponding to the intersection point; On the curve of pressure changing with the G function, find the pressure value corresponding to the G function value at the intersection point, and use it as the closing pressure.

8. The hydraulic fracturing monitoring method based on high-frequency pressure water hammer wave and cepstral analysis according to claim 1, characterized in that, The forward chain reasoning process includes: sequentially scanning the rules in the reasoning rule base, matching whether the condition part of the rule matches the fact objects in the working memory; if the match is successful, triggering the conclusion part of the rule, adding the new fact and its propagated confidence to the working memory, wherein the confidence of the new fact is equal to the rule confidence multiplied by the minimum value of the confidence of each fact object in the condition; repeating the scanning and matching until no new facts are generated, and outputting a comprehensive diagnostic conclusion including fracture effectiveness evaluation, bridge plug sealing comprehensive risk level and fracturing effect comprehensive score.

9. A hydraulic fracturing monitoring system based on high-frequency pressure water hammer wave and cepstral analysis, characterized in that, include: The data acquisition and wave velocity fusion verification module collects high-frequency pressure, ground stress, temperature and vibration in the wellbore, performs adaptive baseline drift correction, and obtains the reference wave velocity, data quality factor, pure pressure sequence and water hammer wave two-way propagation time through multi-source fusion estimation and physical verification. The physical constraint cepstral enhancement and feature inversion module performs cepstral transformation and dynamic time warping alignment based on the pure pressure sequence and reference wave velocity, inverts the spatial wave velocity field profile, and extracts crack reflection features, bubble volume fraction profile, crack opening width profile and liquid ingress intensity profile by combining neural network classification and acoustic attenuation verification. The mechanical model deviation detection and adaptive inversion module combines geostress data and inversion profiles to locate stress shielding layers, uses filtering algorithms to assimilate parameters, and performs wave velocity and stress dual-constraint verification and adaptively switches mechanical models when a deviation in the mechanical model is detected, outputting crack height, average crack width and total crack volume. The evidence dynamic correction and conflict fusion module defines multiple types of bridge plug seal evaluation evidence bodies, uses data quality factors to dynamically correct the basic probability allocation, calculates the conflict coefficient and performs conflict redistribution when the threshold is exceeded, fuses the output bridge plug seal status classification, and performs closure pressure analysis based on the pump stop pressure drop curve to obtain closure pressure and filtration coefficient. The cross-level expert rule reasoning and comprehensive diagnosis module transforms the inversion profile, crack geometric parameters, sealing status classification, and closure pressure parameters into fact objects with confidence levels, loads the rule base, performs forward chain reasoning, and outputs comprehensive diagnostic conclusions.

10. A hydraulic fracturing monitoring system based on high-frequency pressure water hammer wave and cepstral analysis according to claim 9, characterized in that, In the data acquisition and wave velocity fusion verification module, pressure cross-correlation, vibration propagation and theoretical wave velocity estimators are run in parallel. The residuals are calculated by substituting them into the discretized mass and momentum conservation equations. The minimum residual value is selected as the reference wave velocity. If all values ​​exceed the limit, the long short-term memory network is called to output the predicted wave velocity value. In the physical constraint cepstral enhancement and feature inversion module, a forward and reverse dynamic time warping alignment chain is established, the local mapping ratio is extracted to calculate the corrected local wave velocity to form a spatial wave velocity field profile, a lightweight convolutional neural network is used to separate features, and the profile is generated through acoustic attenuation model verification and multi-sensor spatial collaborative constraints. In the mechanical model deviation detection and adaptive inversion module, a state vector assimilation high-frequency pressure transient feature is constructed. When the crack height boundary continuously exceeds the stress shielding layer boundary and the volume error exceeds the limit, the structure is determined to deviate. Based on the stress shielding layer parameters, the mechanical model is automatically switched to continue iterating until the dual convergence constraints are met. In the evidence dynamic correction and conflict fusion module, the probability of performing the correction is increased, maintained, decreased or transferred according to the data quality factor interval. The sum of the product of the probabilities of mutually exclusive propositions between different evidence bodies is calculated as the conflict coefficient. If the threshold is exceeded, the weighted correction method is used for preprocessing before evidence synthesis is performed. The pressure and derivative curves of the dimensionless G function are plotted and extrapolated to obtain the closure pressure. In the cross-level expert rule reasoning and comprehensive diagnosis module, the fact objects in the rule matching working memory are scanned sequentially. If a match is successful, the new fact and its propagated confidence are added to the working memory. The confidence propagation is based on the product of the rule confidence and the minimum confidence of the fact in the condition. The scanning is repeated until no new facts are generated and a diagnosis report is output.

Images