A method for determining the liquid inlet point position by using hydraulic fracturing water hammer signals

Abstract

The present invention discloses a method for determining the liquid inlet point position by using hydraulic fracturing water hammer signals. The method comprises the following steps: adaptively determining the upper and lower cut-off frequencies of a Gaussian band-pass filter by analyzing the main peak of the spectrum through fast Fourier transform; removing spike noise by using an alpha-trimmed filter; extracting the period T of the water hammer signal by using an autocorrelation algorithm; performing a de-mean processing on the signal filtered in the time domain; adaptively determining the relevant parameters of the short-time Fourier transform according to the water hammer signal period T, and performing a short-time Fourier transform operation to obtain a two-dimensional frequency domain matrix; performing Gaussian band-pass filtering and cepstrum transformation on each column of the two-dimensional matrix to construct a cepstrum matrix; distinguishing the bridge plug and the liquid inlet point based on the sign feature 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 has improvements in adaptive filtering, the accuracy of period extraction, and artifact removal in the cepstrum cloud map, and has higher positioning accuracy and better real-time performance.

Description

Technical Field

[0001] The present invention relates to the technical field of oil and gas reservoir development, and particularly relates to a method for determining the liquid injection point position by using the water hammer signal generated when the wellhead pump stops during the hydraulic fracturing process. Background Art

[0002] During the development of oil and gas reservoirs, the hydraulic fracturing technology, as an important means to improve the single-well productivity of low-permeability and tight oil and gas reservoirs, has been widely used. In order to improve the fracturing effect, accurately grasping the actual position where the fracturing fluid enters the formation, that is, the liquid injection point position, is a key issue in fracturing monitoring and effect evaluation. At present, the commonly used fracturing monitoring methods in the industry include microseismic monitoring, perforation acoustic emission monitoring, distributed optical fiber monitoring, and surface pressure curve analysis, etc. Microseismic monitoring can reflect the fracture propagation trajectory, but the equipment cost is high, the construction is complex, and the positioning accuracy is limited; distributed optical fiber monitoring requires pre-burying optical fibers, the construction difficulty is large, and the data interpretation is complex; although the surface pressure curve analysis is simple, it can only qualitatively infer the fracture behavior and it is difficult to accurately locate the liquid injection point.

[0003] However, the existing methods for determining the liquid injection point position by using the water hammer signal generated when the wellhead pump stops during the hydraulic fracturing process mostly only rely on one-time global frequency domain analysis and fixed bandwidth filtering, lacking effective suppression of time-domain spike noise and dynamic bandwidth selection. At the same time, it is difficult to balance the time-frequency resolution and the calculation efficiency, and it is difficult to achieve high-precision and real-time liquid injection point positioning in a complex noise environment. At present, the existing technologies for liquid injection point positioning during the hydraulic fracturing process of oil and gas reservoirs mainly include two types of typical methods: one is the edge computing data processing method based on high-frequency pressure fracture monitoring. This method installs high-frequency pressure gauges at the wellhead four-way valve, and transmits and analyzes the pressure data in real time through edge computing devices and cloud platforms to check the perforation quality, packer leakage, and fracture initiation, etc. However, this method requires a large number of sensors and network devices to be arranged, the system deployment is complex, the construction and debugging cost is high, and it is difficult to achieve precise positioning and rapid response of the liquid injection point. The other is the method for calculating the depth of the main liquid injection point of the fracturing fracture based on the wellhead water hammer signal. This method uses the water hammer pressure wave generated when the pump stops, first removes the direct current from the original time-domain signal and performs Gaussian or masking filtering, then uses window functions such as Hamming window and Hanning window to intercept the water hammer time-domain waveform, and finally calculates the liquid injection point depth through global Fourier transform and frequency difference method. Although this technology avoids the high cost of large-scale microseismic monitoring, it only relies on one-time removal of the direct current component and unified Gaussian / masking filtering, lacking effective suppression of time-domain spike noise, and does not dynamically select the filter bandwidth parameter according to the spectral characteristics of the on-site signal; in addition, the pure frequency domain global analysis method has limitations in terms of time-frequency resolution and calculation efficiency, and it is difficult to obtain high-precision positioning while ensuring real-time performance. Summary of the Invention

[0004] To solve the problems of poor real-time performance and low positioning accuracy in complex noise environments in the prior art, the present invention provides a method for determining the liquid inlet point position using hydraulic fracturing water hammer signals.

[0005] Specifically, the present invention provides a method for determining the liquid inlet point position using hydraulic fracturing water hammer signals, and the method includes the following steps:

[0006] First, perform a fast Fourier transform on the noisy time-domain water hammer signal collected at the wellhead of the hydraulic fracturing site, analyze the main peak of its spectrum, and adaptively determine the upper and lower cut-off frequencies of the Gaussian band-pass filter accordingly;

[0007] Second, process the original noisy water hammer signal in the time domain using an alpha-trimmed filter to remove abnormal spike noise and retain the effective water hammer signal;

[0008] Third, use the autocorrelation algorithm to extract the period T of the water hammer signal from the time-domain filtered signal;

[0009] Fourth, subtract the mean value from the signal filtered in step two to eliminate the DC component;

[0010] Fifth, determine the relevant parameters of the short-time Fourier transform according to the water hammer signal period T obtained in step three, and perform short-time Fourier transform operations on the signal in sub-windows to obtain a two-dimensional frequency-domain matrix S;

[0011] Sixth, calculate the amplitude spectrum for each column of the matrix S, filter it using the Gaussian band-pass filter constructed with the upper and lower cut-off frequencies determined in step one, then take the logarithm of the filtered result and perform the inverse Fourier transform, take the real part to obtain the cepstrum of each window, and splice the cepstra of all windows by column to form a cepstrum matrix;

[0012] Seventh, sum the cepstrum matrix in step six by row, select the row with the largest sum as the bridge plug position, then calculate the time interval between the bridge plug and the liquid inlet point, and multiply it by the wave velocity to complete the time-depth conversion, so as to determine the actual positions of each liquid inlet point and draw a cepstrum cloud map.

[0013] Beneficial effects: The method of the present invention combines frequency-domain adaptive Gaussian band-pass filtering and non-linear time-domain filtering technologies. It can not only dynamically adjust the cut-off frequency of the Gaussian band-pass filter according to different well conditions to adapt to noise changes, but also efficiently remove abnormal spikes and retain the key features of the water hammer signal; use the autocorrelation algorithm to determine the water hammer signal period T, distinguish the bridge plug from the liquid inlet point by combining the positive and negative characteristics of the cepstrum, and calculate the liquid inlet point position by using the bridge plug position, effectively improving the accuracy and real-time performance of the liquid inlet point positioning, and having good engineering applicability and promotion value. Description of the Drawings

[0014] To more clearly illustrate the implementation of the present invention or the existing technical solutions, the following will briefly introduce the drawings required for the description of the embodiments or the existing technologies.

[0015] Figure 1 It is a schematic flow chart of a method for determining the liquid inlet point position by using hydraulic fracturing water hammer signals in an embodiment of the present invention;

[0016] Figure 2 It is a schematic diagram of the autocorrelation function graph and signal period division of the autocorrelation algorithm of a method for determining the liquid inlet point position by using hydraulic fracturing water hammer signals in an embodiment of the present invention;

[0017] Figure 3 It is the cepstrum cloud map of the bridge plug signal of a method for determining the liquid inlet point position by using hydraulic fracturing water hammer signals in an embodiment of the present invention;

[0018] Figure 4 It is the cepstrum cloud map of the bridge plug and the signals of two liquid inlet points of a method for determining the liquid inlet point position by using hydraulic fracturing water hammer signals in an embodiment of the present invention;

[0019] Figure 5 It is the original waveform diagram of the hydraulic fracturing water hammer signal collected in Example 1 of the present invention;

[0020] Figure 6 It is the frequency spectrum distribution diagram of the water hammer signal in Example 1 of the present invention;

[0021] Figure 7 It is the waveform diagram of the strong energy section of the intercepted water hammer signal in Example 1 of the present invention;

[0022] Figure 8 It is the comparison diagram of the original signal and the signal after time-domain filtering after interception in Example 1 of the present invention;

[0023] Figure 9 It is the schematic diagram of extracting the signal period T by the autocorrelation algorithm in Example 1 of the present invention;

[0024] Figure 10 It is the response curve graph of the frequency-domain Gaussian band-pass filter used in Example 1 of the present invention;

[0025] Figure 11 It is the cepstrum cloud map of the water hammer signal obtained based on cepstrum analysis in Example 1 of the present invention. Detailed implementation manners

[0026] To make the objectives, technical solutions and advantages of the present invention clearer and more understandable, the present invention will be further described in detail below in combination with the implementation manners and the drawings. Herein, the illustrative implementation manners of the present invention and their descriptions are used to explain the present invention, but do not limit the present invention.

[0027] Such as Figures 1 to 11As shown, the present invention discloses a method for determining the location of a fluid inlet point using a hydraulic fracturing water hammer signal, comprising the following steps:

[0028] Step 1: First, perform fast Fourier transform on the original noisy time-domain water hammer signal collected at the wellhead of the hydraulic fracturing site to analyze the main peak frequency of its spectrum. f 0 The upper and lower cutoff frequencies of the Gaussian bandpass filter are adaptively set based on the bandwidth adaptively selected on both sides of the main peak. f low and f high , to match the signal characteristics under different well conditions and improve the robustness of subsequent analysis;

[0029] Step 2: Use alpha-trimmed filter to filter the time domain signal, remove isolated spike noise, and retain the effective high-frequency and low-frequency components of the water hammer signal, so as to provide a pure input signal for the next period extraction;

[0030] Step 3: Use the autocorrelation algorithm to extract the period T of the water hammer signal from the time-domain filtered signal, and use the normalized autocorrelation function, combined with the dynamic threshold and the minimum peak distance, to extract the sample point difference Δ between adjacent significant peaks. k , and calculate the water hammer signal period T = Δ k / Fs , using period T to guide the selection of subsequent operation parameters;

[0031] Step 4: Subtract the mean value of the output signal of step 2 as a whole to remove the DC component and avoid low-frequency leakage affecting the accuracy of frequency domain and cepstrum analysis;

[0032] Step 5: Based on the period T calculated in step 3, the window length of the short-time Fourier transform is adaptively determined, and the window function type (preferably Kaiser window), overlap length, and number of Fourier transform points are adaptively determined, so that each analysis window covers at least 2 water hammer signal periods and retains the important periodic characteristics of the water hammer signal; the short-time Fourier transform operation is performed on the signal window to obtain the frequency domain matrix S;

[0033] Step 6: For each column of matrix S S j Compute Magnitude Spectrum | S j |, use the Gaussian bandpass filter with upper and lower cutoff frequencies set in step 1 H ( f ) Perform frequency domain filtering on it; take the logarithm of the filtered amplitude spectrum, then perform inverse Fourier transform and take its real part to obtain the cepstrum vector of each window; the cepstrum vectors of all windows are concatenated column by column to form a two-dimensional cepstrum matrix C ( i, j);

[0034] Step 7: The cepstrum matrix obtained in step 6 C ( i, j ) Sum by row, select energy E i The largest row is taken as the location of the bridge plug signal (corresponding to the larger positive peak of the cepstrum); then the difference between the corresponding row indexes between the bridge plug and the inlet point (corresponding to the smaller negative peak of the cepstrum) is calculated, and the difference is multiplied by the sampling time interval corresponding to the adjacent row indexes in the cepstrum. Δτ and wave speed v , complete the time-to-depth conversion, and thus determine the actual depth position of each liquid inlet point; finally, draw the inverse spectrum cloud map to realize the visual display of the bridge plug and the liquid inlet point position.

[0035] Step 1 is to perform fast Fourier transform on the original noisy time-domain water hammer signal, normalize it and extract the one-sided amplitude spectrum, and then locate the main peak frequency in the spectrum. f 0 ; Then take the adaptive bandwidth on both sides of the main peak to determine the upper and lower cutoff frequencies of the Gaussian bandpass filter f low and f high , complete the determination of the parameters of the adaptive Gaussian bandpass filter, and prepare for the subsequent steps to achieve adaptive frequency domain filtering of the target frequency band.

[0036] The alpha-trimmed filter described in step 2 is a nonlinear time domain filter. In the time domain, a nonlinear filter is used to remove isolated peaks from the original noisy signal without losing the effective frequency components of the water hammer signal. In each sliding window, the alpha-trimmed filter first removes a preset ratio ( α ), and then calculate the average of the remaining samples as the output. α value, which can flexibly balance the noise reduction depth and signal fidelity. α →0, it degenerates into a mean filter. α When the half window length is reached, it degenerates into a median filter, which can not only remove the peaks, but also better retain the effective high-frequency and low-frequency components of the signal, providing a clean time domain signal input for subsequent period extraction.

[0037] The autocorrelation algorithm described in step 3 is a time domain analysis method based on the similarity calculation between the signal and its different delayed versions. The principle is as follows:

[0038] First, calculate the normalized autocorrelation function of the time-domain filtered signal 𝑥[𝑛]:

[0039]

[0040] in,k is the number of delayed sample points, n is the time index, R [ k ] indicates delay k The normalized autocorrelation function value when the samples are x [ n + k ] for delay k The signal after samples, R [0] represents the autocorrelation value at zero lag.

[0041] Then, the non-negative delayed part of the autocorrelation function is intercepted to form the vector 𝑅 + [𝑘], with 0.5max(𝑅 + ) as the dynamic threshold, the number of samples corresponding to the signal sampling frequency 𝐹𝑠 as the minimum peak distance, and call the peak search function (such as MATLAB's findpeaks function) at 𝑅 + [𝑘] locates the first two significant peaks 𝑘1 and 𝑘2. The sample point difference between the two peaks Δ𝑘 = 𝑘2−𝑘1, and finally divided by 𝐹𝑠 to get the period T:

[0042]

[0043] Where T is the signal period, Δ𝑘 is the corresponding sample point difference between the two peaks, and 𝐹𝑠 is the signal sampling frequency.

[0044] Figure 2 The autocorrelation function diagram and signal period division diagram of the autocorrelation algorithm of the present invention are shown. The upper figure is the autocorrelation function diagram (the horizontal axis is the number of delayed samples, and the vertical axis is the autocorrelation value), and the lower figure is the signal period division diagram according to the autocorrelation algorithm (the horizontal axis is the number of data points, and the vertical axis is the signal amplitude). The autocorrelation extraction method has strong robustness to random noise and baseline drift, and can accurately identify repeated water hammer pulse periods in complex noise environments. The accurately extracted period T is an important basis for the adaptive setting of short-time Fourier transform parameters in the subsequent step five.

[0045] Step 4 is to further perform a de-averaging operation on the signal after time domain filtering in step 2. This operation eliminates the DC component or baseline drift by subtracting its average value from the signal sequence. Eliminating the DC component helps to avoid strong low-frequency leakage and pseudo peaks in the subsequent short-time Fourier transform analysis, thereby improving the resolution and accuracy of frequency domain analysis.

[0046] Step 5: Based on the water hammer signal period T extracted in Step 3, adaptively set the parameters of the short-time Fourier transform. Specifically, first determine the window length peak_distance according to the extracted water hammer signal period T, and adaptively set the three parameters of the window function win type, window overlap length overlap, and the number of Fourier transform points nfft. The selection principle of the window length peak_distance is that each window should cover at least 2 periods or more. The window function win type uses the kaiser window, which enables a flexible trade-off between the main lobe width and sidelobe suppression, and its adjustable shape parameter β enables the highest frequency resolution while significantly reducing sidelobe leakage, which is particularly effective in eliminating spectral aliasing and spectral leakage. The window overlap length overlap refers to the number of overlapping samples when two adjacent windows slide. Appropriate overlap can not only ensure smooth transition between adjacent windows and avoid losing mutation information due to window segmentation. The number of Fourier transform points nfft represents the number of transform points used for the discrete Fourier transform of each window, that is, zero-padding the window. It can be equal to the window length, or the next power of 2 greater than the window length, or the next power of 2 greater than the window length + 1. Zero-padding can achieve higher-resolution frequency interpolation. Increasing nfft can improve the frequency domain sampling density, making the spectrum smoother and finer, but the corresponding computational amount will also increase. In the present invention, according to the experimental results, generally take the next power of 2 greater than the window length. Such parameter selection not only ensures sufficient frequency resolution and anti-aliasing ability, but also takes into account the time domain continuity and computational efficiency, laying a solid foundation for the subsequent cepstrum matrix construction. In the specific implementation, the parameters of the short-time Fourier transform mainly include: the input signal data, the signal sampling frequency Fs, the window function win, the window overlap length overlap, and the number of Fourier transform points nfft. Use these parameters to call the stft function of MATLAB to perform the short-time Fourier transform to obtain three parameters: the two-dimensional frequency domain matrix S, the frequency axis f, and the time axis t, providing a basis for subsequent cepstrum denoising and matrix construction. Specifically:

[0047] Call the stft function of MATLAB with these parameters as follows:

[0048] [S, f, t] = stft(data, Fs,...

[0049] 'Window', win,...

[0050] 'OverlapLength', overlap,...

[0051] 'FFTLength', nfft);

[0052] The construction method of the Gaussian band-pass filter described in Step 6 is as follows:

[0053] First, define the frequency coordinate f as:

[0054]

[0055] where 𝐿 is the spectrum length and 𝐹𝑠 is the signal sampling frequency.

[0056]

[0057] where 𝜎 𝐿 is the standard deviation of the Gaussian low-pass filter, 𝜎 𝐻 is the standard deviation of the Gaussian high-pass filter, f low and f high are the upper and lower cut-off frequencies of the Gaussian band-pass filter obtained in Step 1.

[0058] Construct the Gaussian low-pass filter according to the standard deviation obtained from the above steps H low ( f ) is:

[0059]

[0060] Construct the Gaussian high-pass filter H high ( f ) is:

[0061]

[0062] where θ is the coefficient for adjusting the attenuation degree of the high-pass filter.

[0063] Finally, the Gaussian band-pass filter H ( f ) is the product of the Gaussian low-pass filter H low ( f ) and the high-pass filter H high ( f ):

[0064]

[0065] The process of constructing the cepstrum matrix is as follows: For each column of matrix S S j calculate the magnitude spectrum | S j|, a Gaussian band - pass filter with the upper and lower cut - off frequencies set in Step 1 H ( f ) performs frequency - domain filtering on it; take the logarithm of the filtered amplitude spectrum, then perform the inverse Fourier transform and take its real part to obtain the cepstrum vectors of each window; the cepstrum vectors of all windows are concatenated by columns to form a two - dimensional cepstrum matrix C ( i, j ), and its calculation formula is:

[0066]

[0067] Among them, the row index i represents the time delay on the quefrency axis, and the column index j represents the position of the time window corresponding to the original signal on the time axis, S j represents a certain column of the matrix S, H ( f ) represents the Gaussian band - pass filter, real{⋅} represents the operation of taking the real part, that is, only retaining the real part of the result of the inverse Fourier transform, IFFT(⋅) represents the inverse Fourier transform operation, and ln(⋅) represents the operation of taking the natural logarithm.

[0068] In Step 7, the bridge plug corresponds to a relatively large positive cepstrum peak, and the liquid inlet point corresponds to a relatively small negative cepstrum peak. Based on the wave impedance theory, according to the wave impedance theory, the impedance of the medium in the pipeline Z c can be simplified as:

[0069]

[0070] Among them, 𝜌 is the fluid density, 𝑐 is the wave speed, and 𝐴 is the cross - sectional area of the pipe. The reflection coefficient at the interface R is defined as:

[0071]

[0072] Among them, R refers to the reflection coefficient, indicating the direction and amplitude of the energy change when the pressure wave is reflected at the medium interface, and represent the impedances of the medium on both sides of the interface.

[0073] When the medium impedance increases from to (such as a sudden decrease in the cross - sectional area at the bridge plug), 𝑅 > 0, which is manifested as a positive peak in the cepstrum matrix; when the impedance decreases from to (such as an increase in the cross - sectional area at the crack or perforation), 𝑅 < 0, corresponding to the negative peak in the cepstrum.Figure 3 shows the cepstrum cloud diagram of the bridge plug signal. Since the medium impedance at the bridge plug position changes from to , according to the wave impedance reflection theory, the reflection coefficient 𝑅 > 0 here, so it appears as a significant positive peak in the cepstrum matrix. The clearly visible positive peak in the figure corresponds to the distribution position of the bridge plug signal on the quefrency axis, verifying the effectiveness of accurately identifying the bridge plug position by the cepstrum analysis method of the present invention. Figure 4 shows the cepstrum cloud diagram of the bridge plug and the signals of two liquid inlet points. At the bridge plug position, due to the increase in impedance (sudden reduction in cross-sectional area), positive reflection occurs, which appears as a positive peak in the cepstrum; while at the positions of the two liquid inlet points, due to cracks or perforations, the impedance decreases (increase in cross-sectional area), and the reflection coefficient 𝑅 < 0, which appears as a negative peak in the cepstrum. The positive and negative peaks in the figure show a clear spatial separation on the quefrency axis, not only proving that the bridge plug and the liquid inlet points can be effectively distinguished in the cepstrum domain, but also the actual physical distance can be further calculated through the time interval between the peaks to achieve accurate positioning of the liquid inlet points.

[0074] Specifically:

[0075] First, sum the cepstrum matrix by rows, and denote the total energy i of the E ( i ) row as:

[0076]

[0077] Among them, C ( i,j ) is the value of the i row and the j column in the cepstrum matrix, and N is the number of columns.

[0078] By finding the maximum value of all E ( i ), determine the row index r b of the maximum energy, corresponding to the position of the bridge plug reflection:

[0079]

[0080] Among them, argmax(⋅) is the operation of finding the maximum value.

[0081] The algorithm formula for calculating the relative distance between the bridge plug and the liquid inlet point is as follows:

[0082]

[0083] Among them, z ris the actual depth corresponding to the liquid inlet point, q plug is the known bridge plug depth, Δτ is the sampling time interval, v is the propagation speed of the pressure wave in the medium, r b is the row index where the bridge plug is located, r is the row index where the liquid inlet point is located. Wave speed v can be derived from the bridge plug position and the signal period T, and its calculation formula is:

[0084]

[0085] Among them, v is the wave speed, q plug is the known bridge plug depth, and T is the water hammer signal period.

[0086] Taking the time vector 𝑡 as the horizontal axis and the depth vector z r as the vertical axis, corresponding to the matrix C ( i,j ) to make a contour map, the precise mapping between each row on the cepstrum cloud map and the actual depth can be obtained. The method of making the cepstrum cloud map used in Step Seven is implemented through the contourf method of MATLAB.

[0087] In order to further explain the above technical solution, a specific example 1 is provided in the specific embodiment of the present invention for illustration. The actual data used comes from a horizontal well. The target formation of this well is located in a certain formation group. The total measured depth of this well is 3523.53 m (vertical depth) / 5640.00 m (slant depth). The original designed total fracturing section length of this well is 1963 m, the main section length is 90 - 98 m, the average section length is 93.5 m, and there are 21 fracturing sections. In this example, the water hammer signal of the 15th hydraulic fracturing section is selected for analysis. The initial setting depth of the bridge plug in this section is 4512 m, and it is finally adjusted to 4227 m. The following is the data table of the perforation and bridge plug positions of the 15th section:

[0088] Table 1

[0089]

[0090] Table 1 lists the detailed parameters of the 15th fracturing section used for analysis in the embodiment of the present invention, including data such as the bridge plug position, the top depth and bottom depth of each perforation layer, etc. The final setting depth of the bridge plug in this section is 4227 meters, the section length is 87 meters, and it contains 12 layers of perforations, numbered from layer 37 to layer 48. The length of each layer of perforation is 0.45 meters, and the total number of perforation holes is 84.

[0091] The original waveform diagram of the water hammer signal collected in the 15th section is asFigure 5 As shown in the figure. First, perform a fast Fourier transform on the original signal of the 15th hydraulic fracturing water hammer signal collected on site, analyze and locate the main peak in the spectrum, and adaptively determine the upper and lower cut-off frequencies of the Gaussian band-pass filter accordingly. f low and f high . As Figure 6 shown in the spectrum distribution diagram of the water hammer signal, the main peak frequency of the signal spectrum is approximately around 0.08 Hz. However, considering that high-frequency signals are also often accompanied in the bridge plug signal, the upper and lower cut-off frequencies of the Gaussian band-pass filter actually selected in this example f low are 0.1 Hz, f high and 35 Hz. To reduce the computational burden of subsequent processing and reduce the influence of random noise in the pipeline, further intercept the water hammer signal segment with the strongest energy after the pump is stopped for analysis, as Figure 7 shown in the figure.

[0092] Subsequently, apply alpha-trimmed filtering to the original noise signal in the time domain to remove isolated spike noises while retaining the effective components of the water hammer signal. Among them, the window length of the alpha-trimmed filter is set to 7, and 1 point at each end of the extreme value is removed in each window, that is, the removal ratio α is 2 / 7, and the remaining effective samples in the middle are averaged to improve the robustness to local noise while retaining the true waveform changes of the signal to the greatest extent. Figure 8 The figure compares and shows the waveform changes of the original signal and the signal after time-domain filtering: the upper figure is the unfiltered original signal, and the lower figure is the signal after time-domain filtering. It can be seen that the spike noises are significantly weakened, and the overall waveform of the signal is clearer.

[0093] After completing the time-domain filtering, use the autocorrelation algorithm to extract the periodic characteristics of the signal to obtain the period T of the water hammer signal. First, calculate the normalized autocorrelation function of the signal 𝑥[𝑛] after time-domain filtering, and intercept the non-negative delay part to form a vector 𝑅 + [𝑘]. Then, use 0.5max(𝑅 + ) as the dynamic threshold and the number of samples corresponding to the signal sampling frequency 𝐹𝑠 (200 Hz in this example) as the minimum peak distance, and call the peak search function (select the findpeaks function of MATLAB) to locate the first two significant peaks 𝑘1 and 𝑘2 on 𝑅 + [𝑘]. The difference Δ𝑘 = 𝑘2 - 𝑘1 between the corresponding sample points of the two, combined with the known signal sampling frequency 𝐹𝑠, calculate the signal period T:

[0094]

[0095] In this example, the period T of the water hammer signal extracted based on this autocorrelation algorithm is 11.71 seconds. As Figure 9 shown, the upper figure is the normalized autocorrelation function curve corresponding to the signal after time-domain filtering (the horizontal axis is the number of delay samples, and the vertical axis is the autocorrelation value). The circles mark the positions of the main peaks of the detected signals, and the dashed lines are the period dividing lines defined based on the main peak spacing; the lower figure maps the extracted period 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 curve is the waveform of the processed signal, and the dashed lines indicate the "theoretical" start / end positions of each period (obtained from the average period of the autocorrelation), which are used to divide the analysis window, rather than precisely aligning each peak or valley. In this way, the original signal can be divided into several complete period units, providing a stable period reference for subsequent short-time Fourier transform parameter setting and cepstrum analysis.

[0096] After completing the period extraction, perform a mean removal operation on the signal after time-domain filtering, that is, subtract its mean from the entire signal sequence. Subsequently, based on the previously extracted period T ≈ 11.71 seconds, first determine the window length peak_distance according to the extracted period T of the water hammer signal, and adaptively set the three parameters of the window function win type, window overlap length overlap, and the number of Fourier transform points nfft. In this example, the number of samples corresponding to the period of 11.71 seconds is approximately 2300 points. To make each analysis window cover multiple water hammer periods, the window length peak_distance is selected as 7200 points. In this example, the window function win type uses the kaiser window, which can achieve the best balance in eliminating spectral leakage and achieving good frequency resolution, and its adjustable shape parameter β is set to 7 ， so as to significantly reduce the sidelobe leakage while achieving the highest frequency resolution. And in order to retain more precise signal period components, the window overlap length is set close to the window length, and the corresponding sliding step is kept between 64 and 512 points. The number of Fourier transform points nfft is selected as the next power of 2 greater than the window length.

[0097] After setting the above parameters, use the built-in stft function in MATLAB to perform a short-time Fourier transform on the time-domain signal to obtain the two-dimensional frequency-domain matrix S.

[0098] Subsequently, calculate the magnitude spectrum | S j | for each column of the matrix S, and use the adaptively set Gaussian band-pass filter mentioned above (in this example, the upper and lower cut-off frequencies are respectively f low = 0.1 Hz, fhigh Perform frequency-domain denoising processing at a frequency of 35 Hz). Figure 10 The response curve of the Gaussian band-pass filter used is shown. After filtering, take the logarithm of the amplitude spectrum of each column, perform the inverse Fourier transform, and take its real part as the cepstrum vector. Finally, splice all the window cepstrum vectors by column to construct a two-dimensional cepstrum matrix C ( i,j ).

[0099] After completing the construction of the cepstrum matrix, sum it by row to obtain the total energy 𝐸(𝑖) corresponding to each row, and select the index of the row with the maximum energy as the bridge plug reflection position. In this example, the pressure wave velocity is calculated according to the formula

[0100]

[0101] Then, combine the difference between this index and the positions of other negative peaks on the quefrency axis, multiply by the sampling time interval Δτ (0.005 s in this example) and the wave velocity 𝑣 for time-depth conversion. Finally, the actual depth difference between the bridge plug and each liquid inlet point can be obtained, and the precise positioning of the liquid inlet point can be completed accordingly. Finally, use the contourf plotting function of MATLAB to draw the cepstrum cloud map to realize the visual display of the positions of the bridge plug and the liquid inlet points.

[0102] As Figure 11 shown, the cepstrum cloud map intuitively shows the distribution of the bridge plug and multiple liquid inlet points in the depth direction. The position of the positive cepstrum peak corresponding to the bridge plug corresponds to 4227 m after time-depth conversion, which is consistent with the actual setting depth of the bridge plug in the construction, verifying the accuracy of the method. The 5 liquid inlet points are located at depths of 4208 m, 4201 m, 4175 m, 4168 m, and 4157 m respectively, corresponding roughly to the perforation positions. The signals of the bridge plug and the liquid inlet points in the cepstrum cloud map are clearly distinguishable, the positive and negative peaks are clearly separated, and there are no obvious artifacts interfering in the image, with good interpretability and positioning accuracy, verifying the adaptability and engineering application value of the present invention under the complex working conditions of the actual well site.

[0103] The technical means disclosed in the solution of the present invention are not limited to the technical means disclosed in the above embodiments, but also include technical solutions composed of any combination of the above technical features. It should be noted that for those of ordinary skill in the art in this technical field, without departing from the principle of the present invention, several improvements and refinements can be made, and these improvements and refinements are also regarded as the protection scope of the present invention.

Claims

1. A method for determining the liquid inlet point position by using hydraulic fracturing water hammer signals, characterized in that: The method includes the following steps:

1. Perform a fast Fourier transform on the noisy time-domain water hammer signal collected at the wellhead of the hydraulic fracturing site, analyze the main peak of its spectrum, and adaptively determine the upper and lower cut-off frequencies of the Gaussian band-pass filter accordingly; 2. Process the original noisy water hammer signal in the time domain using an alpha-trimmed filter to remove abnormal spike noise and retain the effective water hammer signal; 3. Use the autocorrelation algorithm to extract the period T of the water hammer signal from the time-domain filtered signal; 4. Subtract the mean value from the signal filtered in step 2 to eliminate the DC component; 5. Determine the relevant parameters of the short-time Fourier transform according to the water hammer signal period T obtained in step 3, and perform short-time Fourier transform operations on the signal in windows to obtain the frequency-domain two-dimensional matrix S; 6. Calculate the amplitude spectrum for each column of the matrix S, filter it using the Gaussian band-pass filter constructed with the upper and lower cut-off frequencies determined in step 1, then take the logarithm of the filtered result, perform the inverse Fourier transform, take the real part to obtain the cepstrum of each window, and splice the cepstra of all windows by columns to form a cepstrum matrix; 7. Sum the rows of the cepstrum matrix in step 6, select the row with the largest sum as the position of the bridge plug, then calculate the time interval between the bridge plug and the liquid injection point, multiply it by the wave speed to complete the time-depth conversion, thereby determining the actual positions of each liquid injection point and drawing a cepstrum cloud map; The method for distinguishing the bridge plug and the liquid injection point is based on the different signs of the cepstrum values, where the bridge plug corresponds to a larger positive cepstrum peak and the liquid injection point corresponds to a smaller negative cepstrum peak; By summing the rows of the cepstrum matrix, determine the row index with the largest sum as the position of the bridge plug in the cepstrum, and calculate the actual depth of each liquid injection point according to the difference between its row index and the row index of the liquid injection point, multiplied by the sampling time interval and the wave speed; The algorithm formula for calculating the relative distance between the bridge plug and the liquid injection point is as follows: Wherein, z r is the actual depth corresponding to the liquid inlet point, q plug is the known bridge plug depth, Δτ is the sampling time interval, v is the propagation speed of the pressure wave in the medium, r b is the row index where the bridge plug is located, r is the row index where the liquid inlet point is located.

2. The method for determining the liquid inlet point position by using hydraulic fracturing water hammer signals according to claim 1, wherein: The period T described in step 3 is extracted by the autocorrelation algorithm, including performing autocorrelation calculation on the filtered signal, obtaining the sample point difference corresponding to adjacent significant peaks through peak detection, and then dividing the difference by the sampling frequency to obtain the period T.

3. The method for determining the liquid inlet point position by using hydraulic fracturing water hammer signals according to claim 1, wherein: The short-time Fourier transform parameters described in step 5 are adaptively selected, including: determining the window length according to the water hammer signal period T extracted in step 3, and adaptively setting the window function type, window overlap length, and the number of Fourier transform points.

4. The method for determining the liquid inlet point position by using hydraulic fracturing water hammer signals according to claim 1, wherein: The cepstrum matrix described in step 6 is constructed as follows: Calculate the amplitude spectrum for each column of the matrix S, filter it using a frequency-domain Gaussian band-pass filter, which is the product of a frequency-domain Gaussian low-pass filter and a high-pass filter; Take the logarithm of the filtered amplitude spectrum, perform the inverse Fourier transform, and take the real part to obtain the cepstrum of each window; Splice all the cepstrum vectors by columns in sequence to form a two-dimensional cepstrum matrix.

Images