CN120159397B - Method for determining position of liquid inlet point by utilizing hydraulic fracturing water hammer signal - Google Patents
Method for determining position of liquid inlet point by utilizing hydraulic fracturing water hammer signalInfo
- Publication number
- CN120159397B CN120159397B CN202510642276.8A CN202510642276A CN120159397B CN 120159397 B CN120159397 B CN 120159397B CN 202510642276 A CN202510642276 A CN 202510642276A CN 120159397 B CN120159397 B CN 120159397B
- Authority
- CN
- China
- Prior art keywords
- cepstrum
- signal
- water hammer
- liquid inlet
- determining
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Classifications
-
- E—FIXED CONSTRUCTIONS
- E21—EARTH OR ROCK DRILLING; MINING
- E21B—EARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
- E21B47/00—Survey of boreholes or wells
-
- E—FIXED CONSTRUCTIONS
- E21—EARTH OR ROCK DRILLING; MINING
- E21B—EARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
- E21B47/00—Survey of boreholes or wells
- E21B47/09—Locating or determining the position of objects in boreholes or wells, e.g. the position of an extending arm; Identifying the free or blocked portions of pipes
Landscapes
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Geology (AREA)
- Mining & Mineral Resources (AREA)
- Geophysics (AREA)
- Environmental & Geological Engineering (AREA)
- Fluid Mechanics (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Geochemistry & Mineralogy (AREA)
- Geophysics And Detection Of Objects (AREA)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
Abstract
本发明公开一种利用水力压裂水击信号确定进液点位置的方法，该方法包括如下步骤：通过快速傅里叶变换分析频谱主峰自适应确定高斯带通滤波器上下截止频率；使用alpha‑trimmed滤波器去除尖峰噪声；采用自相关算法提取水击信号的周期T；对时域滤波后的信号进行去均值处理；根据水击信号周期T自适应确定短时傅里叶变换相关参数，执行短时傅里叶变换操作得到频域二维矩阵；对二维矩阵的每一列进行高斯带通滤波与倒谱变换构建倒谱矩阵；基于倒谱值的符号特征区分桥塞与进液点，结合时深转换确定各进液点实际位置，绘制倒谱云图。与现有方法相比，本发明在自适应滤波、周期提取准确性及倒谱云图去伪影方面实现改进，具备更高的定位精度与更优的实时性能。
The present invention discloses a method for determining the position of the liquid inlet point by using a hydraulic fracturing water hammer signal, the method comprising the following steps: adaptively determining the upper and lower cutoff frequencies of a Gaussian bandpass filter by analyzing the main peak of a spectrum through a fast Fourier transform; using an alpha‑trimmed filter to remove spike noise; using an autocorrelation algorithm to extract the period T of the water hammer signal; performing a mean removal process on the signal after time domain filtering; adaptively determining short-time Fourier transform related parameters according to the period T of the water hammer signal, performing a short-time Fourier transform operation to obtain a two-dimensional matrix in the frequency domain; performing a Gaussian bandpass filter and a cepstrum transform on each column of the two-dimensional matrix to construct a cepstrum matrix; distinguishing a bridge plug from a liquid inlet point based on the sign characteristics of the cepstrum value, determining the actual position of each liquid inlet point in combination with time-depth conversion, and drawing a cepstrum cloud map. Compared with the existing methods, the present invention achieves improvements in adaptive filtering, period extraction accuracy, and cepstrum cloud map artifact removal, and has higher positioning accuracy and better real-time performance.
Description
Technical Field
The invention relates to the technical field of oil and gas reservoir development, in particular to a method for determining the position of a liquid inlet point by utilizing a water hammer signal generated by wellhead pump stopping in a hydraulic fracturing process.
Background
In the development process of oil and gas reservoirs, the hydraulic fracturing technology is widely applied as an important means for improving the single well productivity of low-permeability and compact oil and gas reservoirs. In order to improve the fracturing effect, the actual position of the fracturing fluid entering the stratum, namely the position of a fluid inlet point, is accurately mastered, and is a key problem in fracturing monitoring and effect evaluation. Currently, the fracturing monitoring methods commonly used in the industry include microseism monitoring, perforation acoustic emission monitoring, distributed optical fiber monitoring, ground pressure curve analysis and the like. The microseism monitoring can reflect crack propagation tracks, but has the advantages of high equipment cost, complex construction and limited positioning precision, the distributed optical fiber monitoring needs to be embedded with optical fibers, the construction difficulty is high, the data interpretation is complex, and the ground pressure curve analysis is simple and convenient, but only can qualitatively infer crack behaviors, and the accurate positioning of liquid inlet points is difficult.
However, the existing method for determining the position of the liquid inlet point by using the water hammer signal generated by wellhead pump stopping in the hydraulic fracturing process is mostly dependent on disposable global frequency domain analysis and fixed bandwidth filtering, lacks effective suppression and dynamic bandwidth selection of time peak noise, is difficult to consider the time frequency resolution and the calculation efficiency, and is difficult to realize high-precision and real-time liquid inlet point positioning in a complex noise environment. At present, the prior art for positioning liquid inlet points in the hydraulic fracturing process of a hydrocarbon reservoir mainly comprises two typical methods, namely an edge calculation data processing method based on high-frequency pressure fracture monitoring. According to the method, a high-frequency pressure gauge is arranged at a wellhead four-way valve, and pressure data are transmitted and analyzed in real time through edge computing equipment and a cloud platform so as to detect perforation quality, bridge plug leakage, crack initiation and the like. However, the method needs to arrange a large number of sensors and network equipment, has complex system deployment and high construction and debugging cost, and is difficult to realize accurate positioning and quick response of the liquid injection point. The other type is a fracturing fracture main liquid inlet point depth calculation method based on wellhead water hammer signals. According to the method, through a water hammer pressure wave generated when a pump is stopped, direct current and Gaussian or mask filtering is firstly carried out on an original time domain signal, then window functions such as a Hamming window and a Hanning window are adopted to intercept a water hammer time domain waveform, and finally the depth of a liquid injection point is calculated through global Fourier transform and a frequency difference method. Although the technology avoids the high cost of large-scale microseism monitoring, the technology only depends on one-time DC component removal and unified Gaussian/mask filtering, lacks effective suppression of time peak noise, does not dynamically select filter bandwidth parameters according to the frequency spectrum characteristics of field signals, and has limitations in time-frequency resolution and calculation efficiency in a pure frequency domain global analysis mode, so that high-precision positioning is difficult to obtain while real-time performance is ensured.
Disclosure of Invention
In order to solve the problems of poor real-time performance and low positioning accuracy in a complex noise environment in the prior art, the invention provides a method for determining the position of a liquid inlet point by utilizing a hydraulic fracturing water hammer signal.
Specifically, the invention provides a method for determining the position of a liquid inlet point by utilizing a hydraulic fracturing water hammer signal, which comprises the following steps:
1. Carrying out fast Fourier transform on the noisy time domain water hammer signal acquired by the wellhead of the hydraulic fracturing site, analyzing the main peak of the frequency spectrum of the signal, and adaptively determining the upper and lower cut-off frequencies of the Gaussian band-pass filter according to the main peak;
2. processing the original water hammer signal with noise in a time domain by adopting an alpha-trimmed filter to remove abnormal peak noise, and reserving an effective water hammer signal;
3. extracting the period T of the water hammer signal from the time domain filtered signal by adopting an autocorrelation algorithm;
4. subtracting the mean value of the signals filtered in the second step to eliminate the direct current component;
5. determining relevant parameters of short-time Fourier transform according to the water hammer signal period T obtained in the step three, and performing short-time Fourier transform operation on the signal windows to obtain a frequency domain two-dimensional matrix S;
6. Calculating magnitude spectrums for each column of the matrix S, filtering the magnitude spectrums by using Gaussian band-pass filters constructed by the upper and lower cut-off frequencies determined in the step one, taking logarithm of the filtered result, performing inverse Fourier transform, taking the real part to obtain cepstrum of each window, and splicing the cepstrum of all windows according to columns to form a cepstrum matrix;
7. and D, summing the cepstrum matrix in the step six according to rows, selecting the row with the largest sum as the bridge plug position, calculating the time interval between the bridge plug and the liquid inlet points, multiplying the time interval by the wave velocity to finish time-depth conversion, thereby determining the actual position of each liquid inlet point, and drawing a cepstrum cloud picture.
The method has the advantages that the method fuses the frequency domain self-adaptive Gaussian band-pass filtering and nonlinear time domain filtering technologies, can dynamically adjust the cut-off frequency of the Gaussian band-pass filter according to different well conditions to adapt to noise changes, can efficiently remove abnormal peaks and retain key characteristics of water hammer signals, adopts an autocorrelation algorithm to determine the period T of the water hammer signals, combines the positive and negative characteristic differences of a cepstrum to distinguish bridge plugs from liquid inlet points, and effectively improves the accuracy and instantaneity of liquid inlet point positioning by utilizing the bridge plug position to calculate the liquid inlet point position, and has good engineering applicability and popularization value.
Drawings
In order to more clearly illustrate the practice or prior art of the present invention, a brief description of the drawings will be provided below for the purposes of illustrating the embodiments and the drawings used in the description of the prior art.
FIG. 1 is a schematic flow chart of a method for determining a position of a fluid inlet point by utilizing a hydraulic fracturing water hammer signal according to an embodiment of the invention;
FIG. 2 is an autocorrelation function diagram and a signal period dividing schematic diagram of an autocorrelation algorithm of a method for determining a position of a liquid inlet point by using a hydraulic fracturing water hammer signal according to an embodiment of the present invention;
FIG. 3 is a cepstrum cloud chart of a bridge plug signal of a method for determining the position of a fluid inlet point by utilizing a hydraulic fracturing water hammer signal according to the embodiment of the invention;
FIG. 4 is a cepstrum cloud chart of two fluid inlet points and bridge plugs of a method for determining the position of the fluid inlet points by utilizing hydraulic fracturing water hammer signals according to the embodiment of the invention;
FIG. 5 is a raw waveform of the hydraulic fracturing water hammer signal collected in example 1 of the present invention;
FIG. 6 is a spectrum of the water hammer signal of example 1 of the present invention;
FIG. 7 is a waveform of the strong energy segment of the water hammer signal taken in example 1 of the present invention;
FIG. 8 is a graph comparing the original signal intercepted by example 1 of the present invention with the time-domain filtered signal;
FIG. 9 is a schematic diagram of the autocorrelation algorithm extraction signal period T of example 1 of the present invention;
FIG. 10 is a plot of the frequency domain Gaussian bandpass filter response used in example 1 of the invention;
Fig. 11 is a cepstrum cloud image of a water hammer signal obtained based on cepstrum analysis in example 1 of the present invention.
Detailed Description
The present invention will be described in further detail with reference to the following embodiments and the accompanying drawings, in order to make the objects, technical solutions and advantages of the present invention more apparent. The exemplary embodiments of the present invention and the descriptions thereof are used herein to explain the present invention, but are not intended to limit the invention.
As shown in fig. 1 to 11, the embodiment of the invention discloses a method for determining the position of a liquid inlet point by using a hydraulic fracturing water hammer signal, which comprises the following steps:
Firstly, performing fast Fourier transform on an original noisy time domain water hammer signal acquired by a hydraulic fracturing field wellhead, analyzing the frequency f 0 of a main peak of a frequency spectrum of the original noisy time domain water hammer signal, and adaptively setting the upper and lower cut-off frequencies f low and f high of a Gaussian band-pass filter by taking bandwidths adaptively selected from two sides of the main peak as references so as to match signal characteristics under different well conditions and improve subsequent analysis robustness;
filtering the time signal by adopting an alpha-trimmed filter, removing isolated peak noise, and keeping the effective high-frequency and low-frequency components of the water hammer signal at the same time, so as to provide a pure input signal for the next period extraction;
Extracting a period T of the water hammer signal from the time domain filtered signal by using an autocorrelation algorithm, extracting a sample point difference delta k between adjacent significant peaks by using a normalized autocorrelation function and combining a dynamic threshold value and a minimum peak distance, and calculating the period T=delta k/Fs of the water hammer signal according to the sample point difference delta k, and guiding subsequent operation parameter selection by using the period T;
step four, subtracting the average value from the whole output signal of the step two, removing the direct current component, and avoiding the influence of low-frequency leakage on the frequency domain and the cepstrum analysis precision;
step five, based on the period T calculated in the step three, the window length of short-time Fourier transform is adaptively determined, and the window function type (preferably kaiser windows), the overlapping length and the number of Fourier transform points are adaptively determined, so that each analysis window at least covers 2 water-hammer signal periods and retains the important period characteristics of the water-hammer signals;
Step six, calculating magnitude spectrum |S j | for each row S j of the matrix S, carrying out frequency domain filtering on the magnitude spectrum by adopting a Gaussian band-pass filter H (f) with upper and lower cut-off frequencies set in the step one, taking logarithm of the magnitude spectrum after filtering, carrying out inverse Fourier transform and taking the real part of the magnitude spectrum to obtain cepstrum vectors of each window, and splicing the cepstrum vectors of all windows in rows to form a two-dimensional cepstrum matrix C (i, j);
Step seven, summing the cepstrum matrix C (i, j) obtained in the step six according to the rows, selecting the row with the largest energy E i as the position (corresponding to the larger cepstrum positive peak) of the bridge plug signal, then multiplying the difference value by the sampling time interval delta tau and wave velocity v corresponding to the adjacent row index in the cepstrum by the difference value of the corresponding row index between the bridge plug and the liquid inlet point (corresponding to the smaller cepstrum negative peak), completing time-depth conversion, thereby determining the actual depth position of each liquid inlet point, and finally drawing a cepstrum cloud picture, so as to realize the visual display of the bridge plug and the liquid inlet point.
The method comprises the steps of firstly, performing fast Fourier transform on an original noisy time domain water hammer signal, normalizing and extracting a single-side amplitude spectrum, positioning a main peak frequency f 0 in the spectrum, and then respectively obtaining adaptive bandwidths at two sides of the main peak to determine upper and lower cut-off frequencies f low and f high of a Gaussian band-pass filter, so as to complete the determination of parameters of the adaptive Gaussian band-pass filter and prepare for the adaptive frequency domain filtering of a target frequency band in the subsequent step.
And step two, the alpha-trimmed filter is a nonlinear time domain filter, and the nonlinear filter is adopted to remove isolated peaks on the original noise-carrying signal in the time domain without losing the effective frequency component of the water hammer signal. and in each sliding window, the alpha-trimmed filter firstly eliminates the maximum value and the minimum value of the preset proportion (alpha), and then averages the rest samples to be used as output. By adjusting the alpha value, the noise removal depth and the signal fidelity can be flexibly balanced, the noise is degraded into an average filter when alpha is 0, and the noise is degraded into a median filter when alpha reaches the half window length, so that the peak can be removed, the effective high-frequency and low-frequency components of the signal can be well reserved, and clean time domain signal input is provided for the subsequent period extraction.
The autocorrelation algorithm in the third step is a time domain analysis method based on similarity calculation of different delay versions of signals, and the principle is as follows:
First, a normalized autocorrelation function is calculated for the time-domain filtered signal x [ n ]:
Where k is the number of delayed samples, n is the time index, R [ k ] is the normalized autocorrelation function value when delaying k samples, x [ n+k ] is the signal after delaying k samples, and R [0] is the autocorrelation value when delaying zero.
The non-negative delay portion of the autocorrelation function is then truncated to form vector R + [ k ], and a peak search function (e.g., the findpeaks function of MATLAB) is invoked to locate the first two significant peaks k 1 and k 2 on R + [ k ] with 0.5max (R +) as the dynamic threshold and the number of samples corresponding to the signal sampling frequency Fs as the minimum peak distance. Sample point difference delta k=k 2−k1 corresponding to the two peaks is divided by Fs to obtain period T:
wherein T is a signal period, Δk is a sample point difference corresponding to two peaks, and Fs is a signal sampling frequency.
Fig. 2 shows an autocorrelation function diagram and a signal period division schematic diagram of the autocorrelation algorithm of the present invention. The upper graph is an autocorrelation function graph (the horizontal axis is the number of delay samples, the vertical axis is the autocorrelation value), and the lower graph is a signal period division diagram divided according to an autocorrelation algorithm (the horizontal axis is the number of data points, and the vertical axis is the signal amplitude). The autocorrelation extraction method has stronger robustness to random noise and baseline drift, can accurately identify repeated water-hammer pulse periods in a complex noise environment, and the accurately extracted period T is an important basis for short-time Fourier transform parameter self-adaptive setting in the subsequent step five.
And step four, further executing an average value removing operation on the signal subjected to the time domain filtering in the step two. This operation eliminates the dc component or baseline drift by subtracting its average value from the signal sequence. Eliminating the direct current component helps to avoid strong low frequency leakage and false peaks in the subsequent short-time Fourier transform analysis, thereby improving the resolution and accuracy of the frequency domain analysis.
And step five, based on the water hammer signal period T extracted in the step three, adaptively setting all parameters of short-time Fourier transform. Specifically, firstly, determining a window length peak_distance according to an extracted water hammer signal period T, and adaptively setting three parameters of a window function win type, a window overlap length overlap and a Fourier transform point number nfft. The window length peak_distance is selected such that each window covers at least 2 periods and more. The window function win adopts kaiser windows, so that flexible trade-off between main lobe width and side lobe suppression can be realized, and the adjustable shape parameter beta enables side lobe leakage to be obviously reduced while the frequency resolution is highest, and the method is particularly effective to eliminate spectrum aliasing and spectrum leakage. The window overlap length overlap refers to the number of samples that overlap when two adjacent windows slide. Proper overlapping can ensure smooth transition between adjacent windows and avoid losing mutation information due to window segmentation. Fourier transform points nfft represent the transform points used in performing the discrete fourier transform on each window, i.e., performing the zero padding operation on the window. It may be equal to the window length, or may be an integer power greater than the next 2 of the window length or an integer power +1 greater than the next 2 of the window length, with zero padding enabling higher resolution frequency interpolation. Increasing nfft can increase the frequency domain sampling density, making the spectrum smoother and finer, but the corresponding computation increases. In the present invention, according to the experimental results, the integer power of the next 2 greater than the window length is generally taken. The parameter selection ensures enough frequency resolution and anti-aliasing capability, and also considers time domain continuity and calculation efficiency, thereby laying a solid foundation for the subsequent cepstrum matrix construction. In a specific implementation, parameters of short-time Fourier transform mainly comprise input signal data, signal sampling frequency Fs, window function win, window overlap length overlap and Fourier transform point number nfft, and three parameters, namely a frequency domain two-dimensional matrix S, a frequency axis f and a time axis t, are obtained by using the parameters to call a stft function of MATLAB for short-time Fourier transform, so that a foundation is provided for subsequent cepstrum denoising and matrix construction. Specifically:
the stft function that calls MATLAB using these parameters is as follows:
[S, f, t] = stft(data, Fs, ...
'Window', win, ...
'OverlapLength', overlap, ...
'FFTLength', nfft);
The construction mode of using the Gaussian band-pass filter in the sixth step is as follows:
First, defining a frequency coordinate f as:
Where L is the spectral length and Fs is the signal sampling frequency.
Wherein, σ kou is the standard deviation of the Gaussian low-pass filter, σ kou is the standard deviation of the Gaussian high-pass filter, and f low and f high are the upper and lower cut-off frequencies of the Gaussian band-pass filter obtained in the step one.
According to the standard deviation obtained in the above steps, a gaussian low pass filter H low (f) is constructed as follows:
The gaussian high pass filter H high (f) is constructed as:
where θ is a coefficient for adjusting the attenuation degree of the high-pass filter.
The final gaussian bandpass filter H (f) is the product of the gaussian low-pass filter H low (f) and the high-pass filter H high (f):
The cepstrum matrix construction process comprises the steps of calculating magnitude spectrum I S j I of each column S j of a matrix S, carrying out frequency domain filtering on the magnitude spectrum by adopting a Gaussian band-pass filter H (f) with upper and lower cut-off frequencies set in the first step, taking logarithm of the magnitude spectrum after filtering, carrying out inverse Fourier transform and taking a real part of the magnitude spectrum to obtain cepstrum vectors of each window, and splicing the cepstrum vectors of all windows according to columns to form a two-dimensional cepstrum matrix C (i, j), wherein the calculation formula is as follows:
Where row index i represents the time delay on the inverse frequency axis (quefrency), column index j represents the corresponding time window position of the original signal on the time axis, S j represents a certain column of the matrix S, H (f) represents a gaussian bandpass filter, real { is represented by taking the real part, i.e. retaining only the real part of the inverse fourier transform result, IFFT (% is represented by the inverse fourier transform operation, ln (% is) is represented by the natural logarithm operation.
In the seventh step, the positive peak of the cepstrum corresponding to the larger bridge plug and the negative peak of the cepstrum corresponding to the smaller liquid inlet point are based on the wave impedance theory, and according to the wave impedance theory, the medium impedance Z c in the pipeline can be simplified into:
Where ρ is the fluid density, c is the wave velocity, and A is the tube cross-sectional area. The reflection coefficient R at the interface is defined as:
wherein R refers to a reflection coefficient, which represents the change direction and amplitude of energy when the pressure wave is reflected at the medium interface, AndRepresenting the impedance of the medium on both sides of the interface.
When the medium impedance is defined byTo increase to(E.g. a sudden decrease in cross-sectional area at the bridge plug), R > 0, appears as a positive peak in the cepstrum matrix, and when the impedance is changed fromIs reduced to(E.g., an increase in cross-sectional area at the fracture or perforation), R < 0 corresponds to a negative peak in the cepstrum. FIG. 3 shows a cepstrum cloud of the bridge plug signal, since the dielectric impedance at the bridge plug location is determined byTo increase toAccording to the theory of wave impedance reflection, the reflection coefficient R > 0 here, therefore, appears as a significant positive peak in the cepstrum matrix. The positive peak clearly seen in the graph corresponds to the distribution position of the bridge plug signal on the quefrency axis, and the effectiveness of accurately identifying the bridge plug position by the cepstrum analysis method is verified. Fig. 4 shows a bridge plug and a cepstrum cloud of two liquid inlet point signals, wherein the bridge plug position shows positive peaks in the cepstrum due to positive reflection caused by increased impedance (sudden decrease of sectional area), and the two liquid inlet point positions show negative peaks in the cepstrum due to reduced impedance (increase of sectional area) caused by cracks or perforations, and the reflection coefficient R is less than 0. In the figure, positive and negative peaks show clear space separation on quefrency axes, so that the bridge plug and the liquid inlet point can be effectively distinguished in a cepstrum domain, and the actual physical distance can be further calculated through the time interval between peaks, so that the accurate positioning of the liquid inlet point is realized.
Specifically:
first, the cepstral matrices are summed row by row, and the total energy E (i) of the i th row is noted as:
wherein, C (i, j) is the value of the ith row and the jth column in the cepstrum matrix, and N is the column number.
By maximizing all E (i), the row index r b with the greatest energy is determined, corresponding to the position of the bridge reflection:
Wherein argmax (#) is the maximum operation.
The algorithm formula for calculating the relative distance between the bridge plug and the liquid inlet point is as follows:
Where z r is the actual depth corresponding to the feed point, q plug is the known bridge plug depth, Δτ is the sampling time interval, v is the propagation velocity of the pressure wave in the medium, r b is the row index of the bridge plug, and r is the row index of the feed point. The wave velocity v can be derived from the bridge plug position and the signal period T, and the calculation formula is as follows:
where v is the wave velocity, q plug is the known bridge plug depth, and T is the water hammer signal period.
And taking the time vector t as a horizontal axis and the depth vector z r as a vertical axis, and taking the corresponding matrix C (i, j) as a contour map to obtain the accurate mapping between each row and the actual depth on the cepstrum cloud picture. The cepstrum cloud pattern drawing mode used in the seventh step is realized by a contourf method of MATLAB.
For further explanation of the above technical solutions, a specific embodiment of the present invention is provided as a specific example 1 for explanation. The actual data used was from a horizontal well. The well target stratum is positioned in a stratum group, and the well has a well completion depth of 3523.53m (vertical depth)/5640.00 m (inclined depth). The original design fracturing total section length of the well is 1963m, the main section length is 90-98 m, the average section length is 93.5m, and the fracturing section length is 21 sections. In this example, the 15 th hydraulic fracturing water hammer signal is selected for analysis, and the initial bridge plug setting depth is 4512 m, and finally is adjusted to 4227 m. The following is a 15 th stage perforation and bridge plug position data table:
TABLE 1
Table 1 lists the detailed parameters of the 15 th stage fracturing stage for analysis in the examples of the present invention, including bridge plug position, top and bottom depths of each perforation layer, and the like. The final setting depth of the bridge plug is 4227 m, the length of the bridge plug is 87 m, and the bridge plug comprises 12 layers of perforation, and the number of the perforation ranges from 37 layers to 48 layers. Each layer of perforation length is 0.45 m, and the total perforation number is 84.
The original waveform diagram of the hydraulic fracturing water hammer signal collected in the 15 th section is shown in fig. 5. Firstly, performing fast Fourier transform on an original signal of a 15 th-section hydraulic fracturing water hammer signal acquired on site, analyzing and positioning a main peak in a frequency spectrum, and adaptively determining upper and lower cut-off frequencies f low and f high of a Gaussian band-pass filter according to the main peak. As shown in the spectrum distribution diagram of the water hammer signal in fig. 6, the main peak frequency of the signal spectrum is about 0.08 Hz, however, considering that the bridge signal is accompanied by high frequency signal, the upper and lower cut-off frequency f low of the gaussian band-pass filter actually selected in this example is 0.1 hz, and f high is 35 Hz. To reduce the operational burden of subsequent processing and reduce the influence of random noise in the pipeline, the water hammer signal segment with the strongest energy after the pump is stopped is further intercepted for analysis, as shown in fig. 7.
Then, the original noise signal is subjected to alpha-trimmed filtering in the time domain, isolated peak noise is removed, and meanwhile the effective components of the water hammer signal are reserved. The window length of the alpha-trimmed filter is set to be 7, 1 point of each extreme value at two ends is removed in each window, namely the removal proportion alpha is 2/7, and middle effective samples are reserved for averaging, so that the robustness to local noise is improved, and meanwhile, the real waveform change of signals is reserved to the greatest extent. Fig. 8 shows waveform change of the original signal and the time-domain filtered signal in a comparison manner, wherein the upper graph is an unfiltered original signal, the lower graph is a time-domain filtered signal, peak noise is obviously weakened, and the overall waveform of the signal is clearer.
After time domain filtering is completed, the signal is subjected to periodic feature extraction by adopting an autocorrelation algorithm so as to obtain the period T of the water hammer signal. A normalized autocorrelation function is first calculated for the time-domain filtered signal x [ n ] and the non-negative delay component is truncated to form a vector R + [ k ]. Next, with 0.5max (R +) as a dynamic threshold, with the number of samples corresponding to the signal sampling frequency Fs (200 Hz in this example) as a minimum peak distance, a peak search function (findpeaks function of MATLAB is selected) is called to locate the first two significant peaks k 1 and k 2 on R + [ k ], and a sample point difference Δk=k 2 − k1 corresponding to the two is calculated, and the signal period T is calculated by combining the known signal sampling frequency Fs:
In this example, the period T of the water hammer signal extracted based on the autocorrelation algorithm is 11.71 seconds. As shown in fig. 9, the upper graph shows a normalized autocorrelation function curve (the horizontal axis shows the number of delay samples, and the vertical axis shows the autocorrelation value) corresponding to the time-domain filtered signal. The circles mark the detected signal main peak positions, the dashed lines are the periodic demarcations defined by the main peak spacing, the lower graph maps the extracted periodic structure back to the original pressure signal (the horizontal axis is the number of data points and the vertical axis is the signal amplitude), the curves are the processed signal waveforms, the dashed lines mark the "theoretical" start/end positions of each period (from the autocorrelation averaging period) which are used to divide the analysis window without having to precisely align each peak or trough. In this way, the original signal can be divided into a plurality of complete periodic units, and stable periodic references are provided for subsequent short-time Fourier transform parameter setting and cepstrum analysis.
After the period extraction is completed, a de-averaging operation is performed on the time-filtered signal, i.e. its average value is subtracted from the whole signal sequence. Then, based on the aforementioned extracted period t≡11.71 seconds, first, the window length peak_distance is determined according to the extracted water hammer signal period T, and three parameters of the window function win type, the window overlap length overlap, and the fourier transform point number nfft are adaptively set. In this example, the number of sample points corresponding to the period of 11.71 seconds is about 2300 points, and the window length peak_distance is 7200 points so that each analysis window covers a plurality of water hammer periods. In this example, the window function win adopts kaiser windows, so that the best balance can be achieved in terms of eliminating frequency spectrum leakage and achieving good frequency resolution, and the adjustable shape parameter beta is set to 7, so that sidelobe leakage is obviously reduced while the frequency resolution is highest. While in order to preserve more accurate signal period components, the window overlap length is set close to the window length, with the corresponding sliding step size remaining between 64 and 512 points. The fourier transform points nfft are chosen to be greater than the next integer power of 2 of the window length.
After the parameters are set, performing short-time Fourier transform on the time domain signals by utilizing a stft function of MATLAB, thereby obtaining a frequency domain two-dimensional matrix S.
Subsequently, the magnitude spectrum |s j | is calculated for each column of the matrix S, and frequency domain denoising processing is performed using the foregoing adaptively set gaussian band-pass filter (upper and lower cut-off frequencies in this example are f low = 0.1 Hz,fhigh =35 Hz, respectively). Fig. 10 shows the response curve of the gaussian band-pass filter used. And taking the logarithm of each column of the magnitude spectrum after filtering, performing inverse Fourier transform, taking the real part as a cepstrum vector, and finally splicing all window cepstrum vectors according to the columns to construct a two-dimensional cepstrum matrix C (i, j).
After the cepstrum matrix is built, summing the cepstrum matrix according to the rows to obtain total energy E (i) corresponding to each row, and selecting an index of the row with the maximum energy as a bridge plug reflection position. In this example, the pressure wave velocity is calculated as
The difference between the index and the other negative peak position on quefrency axis is multiplied by the sampling time interval delta tau (0.005 s in this example) and the wave velocity v for time-depth conversion. Finally, the actual depth difference between the bridge plug and each liquid inlet point can be obtained, and the accurate positioning of the liquid inlet points is finished according to the actual depth difference. And finally, drawing a cepstrum cloud picture by adopting contourf drawing functions of MATLAB, and realizing visual display of bridge plugs and liquid inlet points.
As shown in fig. 11, the cepstrum cloud chart intuitively shows the distribution condition of the bridge plug and a plurality of liquid inlet points in the depth direction, the cepstrum positive peak position corresponding to the bridge plug corresponds to 4227 m after deep conversion with time, is consistent with the actual bridge plug setting depth in construction, and verifies the accuracy of the method. The 5 fluid feed points are located at depths 4208m, 4201m, 4175m, 4168m, 4157m, respectively, and correspond approximately to perforation locations. The bridge plug and the liquid inlet point signals in the cepstrum cloud picture are clear and distinguishable, the positive and negative peaks are obviously separated, no obvious artifact interference exists in the image, the cepstrum cloud picture has good interpretability and positioning accuracy, and the adaptability and engineering application value of the cepstrum cloud picture under the actual complex working condition of the well site are verified.
The technical means disclosed by the scheme of the invention is not limited to the technical means disclosed by the embodiment, and also comprises the technical scheme formed by any combination of the technical features. It should be noted that modifications and adaptations to the invention may occur to one skilled in the art without departing from the principles of the present invention and are intended to be within the scope of the present invention.
Claims (4)
1. A method for determining a liquid inlet point position by utilizing a hydraulic fracturing water hammer signal is characterized by comprising the following steps:
1. Carrying out fast Fourier transform on the noisy time domain water hammer signal acquired by the wellhead of the hydraulic fracturing site, analyzing the main peak of the frequency spectrum of the signal, and adaptively determining the upper and lower cut-off frequencies of the Gaussian band-pass filter according to the main peak;
2. processing the original water hammer signal with noise in a time domain by adopting an alpha-trimmed filter to remove abnormal peak noise, and reserving an effective water hammer signal;
3. extracting the period T of the water hammer signal from the time domain filtered signal by adopting an autocorrelation algorithm;
4. subtracting the mean value of the signals filtered in the second step to eliminate the direct current component;
5. determining relevant parameters of short-time Fourier transform according to the water hammer signal period T obtained in the step three, and performing short-time Fourier transform operation on the signal windows to obtain a frequency domain two-dimensional matrix S;
6. Calculating magnitude spectrums for each column of the matrix S, filtering the magnitude spectrums by using Gaussian band-pass filters constructed by the upper and lower cut-off frequencies determined in the step one, taking logarithm of the filtered result, performing inverse Fourier transform, taking the real part to obtain cepstrum of each window, and splicing the cepstrum of all windows according to columns to form a cepstrum matrix;
7. the method comprises the steps of summing up the cepstrum matrix in the step six according to rows, selecting the row with the largest sum as the bridge plug position, calculating the time interval between the bridge plug and the liquid inlet points, multiplying the time interval by the wave velocity to finish time depth conversion, determining the actual position of each liquid inlet point, drawing a cepstrum cloud chart, distinguishing the bridge plug from the liquid inlet points based on different cepstrum value signs, wherein the bridge plug corresponds to a larger cepstrum positive peak and the liquid inlet point corresponds to a smaller cepstrum negative peak, determining the row index with the largest sum as the position of the bridge plug in the cepstrum by summing up the cepstrum matrix, multiplying the sampling time interval and the wave velocity according to the difference value between the row index of the bridge plug and the liquid inlet point, and calculating the actual depth of each liquid inlet point, wherein the algorithm formula for calculating the relative distance between the bridge plug and the liquid inlet point is as follows:
Where z r is the actual depth corresponding to the feed point, q plug is the known bridge plug depth, Δτ is the sampling time interval, v is the propagation velocity of the pressure wave in the medium, r b is the row index of the bridge plug, and r is the row index of the feed point.
2. The method of determining a position of a fluid intake point by using a hydraulic fracturing water hammer signal according to claim 1, wherein the period T in the third step is extracted by an autocorrelation algorithm, the filtered signal is subjected to autocorrelation calculation, a sample point difference corresponding to adjacent significant peaks is obtained by peak detection, and the difference is divided by a sampling frequency to obtain the period T.
3. The method for determining the position of a fluid inlet point by utilizing a hydraulic fracturing water hammer signal according to claim 1, wherein the short-time Fourier transform parameter in the fifth step is selected in a self-adaptive manner, and the method comprises the steps of determining a window length according to a water hammer signal period T extracted in the third step, and setting a window function type, a window overlapping length and a Fourier transform point number in a self-adaptive manner.
4. The method for determining the position of a liquid inlet point by utilizing a hydraulic fracturing water hammer signal according to claim 1, wherein the cepstrum matrix in the step six is constructed by calculating a magnitude spectrum for each column of a matrix S, filtering the magnitude spectrum by adopting a frequency domain Gaussian band-pass filter, wherein the frequency domain Gaussian band-pass filter is the product of a frequency domain Gaussian low-pass filter and a high-pass filter, performing Fourier inverse transformation on the filtered magnitude spectrum after taking logarithms, taking a real part, obtaining each window cepstrum, and sequentially splicing all cepstrum vectors according to columns to form a two-dimensional cepstrum matrix.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202510642276.8A CN120159397B (en) | 2025-05-19 | 2025-05-19 | Method for determining position of liquid inlet point by utilizing hydraulic fracturing water hammer signal |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202510642276.8A CN120159397B (en) | 2025-05-19 | 2025-05-19 | Method for determining position of liquid inlet point by utilizing hydraulic fracturing water hammer signal |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN120159397A CN120159397A (en) | 2025-06-17 |
| CN120159397B true CN120159397B (en) | 2025-07-22 |
Family
ID=96007093
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202510642276.8A Active CN120159397B (en) | 2025-05-19 | 2025-05-19 | Method for determining position of liquid inlet point by utilizing hydraulic fracturing water hammer signal |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN120159397B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN120871096A (en) * | 2025-09-23 | 2025-10-31 | 中国船舶集团有限公司第七一五研究所 | Interference equipment identification method based on forwarding harmonic characteristics |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114341462A (en) * | 2019-08-28 | 2022-04-12 | 斯伦贝谢技术有限公司 | Method for determining the position of a lowerable object in a wellbore |
| CN119025840A (en) * | 2024-08-21 | 2024-11-26 | 成都理工大学 | A method for calculating the main liquid inlet depth of a fracturing fracture based on wellhead water hammer signal |
Family Cites Families (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6847737B1 (en) * | 1998-03-13 | 2005-01-25 | University Of Houston System | Methods for performing DAF data filtering and padding |
| EP3480628A4 (en) * | 2016-07-01 | 2020-04-15 | Services Petroliers Schlumberger | METHOD AND SYSTEM FOR DETECTING IN THE WELLBORE OBJECTS REFLECTING A HYDRAULIC SIGNAL |
| WO2021016412A1 (en) * | 2019-07-23 | 2021-01-28 | Seismos, Inc. | Detecting operational anomalies for continuous hydraulic fracturing monitoring |
| CN111550230B (en) * | 2020-04-02 | 2021-03-02 | 中国石油大学(北京) | System for performing fracturing diagnosis based on water hammer pressure wave signal and fracturing diagnosis method |
| CN114239656B (en) * | 2021-12-17 | 2023-04-07 | 中国石油大学(北京) | Underground event positioning method and device based on pump stopping pressure signal |
| CN118423620A (en) * | 2024-04-25 | 2024-08-02 | 江南大学 | A method and system for detecting leakage of underground water supply pipe network based on audio signal |
| CN118349815A (en) * | 2024-04-29 | 2024-07-16 | 中国石油大学(北京) | Water hammer pressure wave signal filtering effect evaluation method and device |
-
2025
- 2025-05-19 CN CN202510642276.8A patent/CN120159397B/en active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114341462A (en) * | 2019-08-28 | 2022-04-12 | 斯伦贝谢技术有限公司 | Method for determining the position of a lowerable object in a wellbore |
| CN119025840A (en) * | 2024-08-21 | 2024-11-26 | 成都理工大学 | A method for calculating the main liquid inlet depth of a fracturing fracture based on wellhead water hammer signal |
Also Published As
| Publication number | Publication date |
|---|---|
| CN120159397A (en) | 2025-06-17 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN120159397B (en) | Method for determining position of liquid inlet point by utilizing hydraulic fracturing water hammer signal | |
| CN109285561B (en) | Ship propeller cavitation noise modulation spectrum feature fidelity enhancement method based on self-adaptive window length | |
| KR101294681B1 (en) | Apparatus and method for processing weather signal | |
| CN108548957A (en) | The double-spectrum analysis method being combined based on circular modulating frequency spectrum and segmentation cross-correlation | |
| CN110299141B (en) | An acoustic feature extraction method for recording playback attack detection in voiceprint recognition | |
| CN107329115B (en) | Parameter Estimation Method of LFM Signal Based on GCRBF Network | |
| CN116110417B (en) | A data enhancement method and device for ultrasonic voiceprint anti-counterfeiting | |
| CN108670291A (en) | The heart sound kind identification method of improved MFCC is combined based on EMD | |
| CN114563824A (en) | Identification method for second-order multiple synchronous extrusion polynomial chirp transform thin reservoir | |
| CN112364767A (en) | High-pressure circulating pump mechanical signal feature extraction method based on fast spectrum correlation | |
| CN113389541B (en) | High-precision extraction method for oil well working fluid level signal | |
| CN108682433A (en) | The heart sound kind identification method of first-order difference coefficient based on MFCC | |
| CN110838302B (en) | Audio frequency segmentation method based on signal energy peak identification | |
| CN111736222A (en) | Single-shot data signal-to-noise ratio determining method and device | |
| CN113740902A (en) | Method for identifying pinch-out point of geologic body based on generalized S transformation | |
| CN110275150B (en) | Variable acceleration moving target coherent accumulation method based on empirical mode decomposition and iterative endpoint fitting | |
| CN106019377A (en) | Two-dimensional seismic exploration noise removing method based on time-space-domain frequency reduction model | |
| CN112086105B (en) | Target identification method based on Gamma atom sub-band continuous spectrum characteristics | |
| CN116203636A (en) | A Seismic Identification Method for Deep Complex Fault System in Foreland Thrust Belt | |
| CN110806567B (en) | Ground clutter false alarm point eliminating method and system based on FIR Doppler filtering information | |
| CN111239814A (en) | Shallow profile data mechanical interference suppression method based on same-phase axis frequency division tracking smoothing | |
| EP2491412B1 (en) | Method for detecting a wavefront corresponding to an event in a signal received by a detector | |
| CN118244359A (en) | High-precision STC calculation method suitable for low signal-to-noise ratio acoustic logging data | |
| CN120877745B (en) | A method for optimal selection of DEMON demodulation frequency band under false target interference | |
| CN119511353B (en) | Method and device for estimating focus wavelet |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |