Method for determining position of liquid inlet point by using hydraulic fracturing water hammer signal

Abstract

The invention discloses a method for determining the position of a liquid inlet point by using a hydraulic fracturing water hammer signal. The method comprises the following steps: analyzing a main peak of a frequency spectrum through fast Fourier transform, and adaptively determining upper and lower cut-off frequencies of a Gaussian band-pass filter; an alpha-trimmed filter is used for removing peak noise, and an alpha-trimmed filter is used for removing peak noise; extracting a period T of the water hammer signal by adopting an autocorrelation algorithm; carrying out mean value removal processing on the signal after time domain filtering; short-time Fourier transform related parameters are adaptively determined according to the water attack signal period T, and short-time Fourier transform operation is executed to obtain a frequency domain two-dimensional matrix; performing Gaussian band-pass filtering and cepstrum transformation on each column of the two-dimensional matrix to construct a cepstrum matrix; and distinguishing the bridge plug and the liquid inlet points based on the symbolic 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 picture. Compared with the prior art, the method is improved in the aspects of adaptive filtering, period extraction accuracy and cepstrum cloud picture artifact removal, and has higher positioning precision 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 is stopped during the hydraulic fracturing process. Background Art

[0002] In the process of oil and gas reservoir development, 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 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 is stopped 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 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 typical methods: one is an edge computing data processing method based on high-frequency pressure fracture monitoring. This method installs a high-frequency pressure gauge at the wellhead four-way valve, and transmits and analyzes the pressure data in real time through edge computing equipment and a cloud platform 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 injection point. The other is a 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 is stopped, 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 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 not dynamically selecting 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 injection point position using hydraulic fracturing water hammer signals.

[0005] Specifically, the present invention provides a method for determining the liquid injection point position using hydraulic fracturing water hammer signals, and the method includes the following steps: 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; 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; Third, use the autocorrelation algorithm to extract the period T of the water hammer signal from the time-domain filtered signal; Fourth, subtract the mean value of the signal filtered in step two from the signal to eliminate the DC component; 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 the frequency-domain two-dimensional matrix S; 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; Seventh, sum the rows of the cepstrum matrix in step six, select the row with the largest sum as the bridge plug position, then calculate the time interval between the bridge plug and the liquid injection point, and multiply it by the wave velocity to complete the time-depth conversion, so as to determine the actual positions of each liquid injection point and draw a cepstrum cloud map.

[0006] Beneficial effects: The method of the present invention combines frequency-domain adaptive Gaussian band-pass filtering and non-linear time-domain filtering techniques. 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, combine the positive and negative cepstrum feature differences to distinguish the bridge plug from the liquid injection point, and calculate the liquid injection point position by using the bridge plug position, effectively improving the accuracy and real-time performance of the liquid injection point positioning, and having good engineering applicability and promotion value. Description of the Drawings

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

[0008] Figure 1It is a schematic flow chart of a method for determining the liquid inlet point position using hydraulic fracturing water hammer signals according to an embodiment of the present invention; Figure 2 It is a schematic diagram of the autocorrelation function graph of the autocorrelation algorithm and the signal period division of a method for determining the liquid inlet point position using hydraulic fracturing water hammer signals according to an embodiment of the present invention; Figure 3 It is the cepstrum cloud map of the bridge plug signal of a method for determining the liquid inlet point position using hydraulic fracturing water hammer signals according to an embodiment of the present invention; 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 using hydraulic fracturing water hammer signals according to an embodiment of the present invention; Figure 5 It is the original waveform diagram of the hydraulic fracturing water hammer signal collected in Example 1 of the present invention; Figure 6 It is the frequency spectrum distribution diagram of the water hammer signal in Example 1 of the present invention; Figure 7 It is the waveform diagram of the strong energy segment of the water hammer signal intercepted in Example 1 of the present invention; 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; Figure 9 It is the schematic diagram for extracting the signal period T by the autocorrelation algorithm in Example 1 of the present invention; Figure 10 It is the response curve graph of the frequency domain Gaussian band - pass filter used in Example 1 of the present invention; 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. Specific implementation mode

[0009] 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 conjunction with the implementation modes and the accompanying drawings. Here, the illustrative implementation modes of the present invention and their descriptions are used to explain the present invention, but do not limit the present invention.

[0010] As Figures 1 to 11 shown, an embodiment of the present invention discloses a method for determining the liquid inlet point position using hydraulic fracturing water hammer signals, including the following steps: Step 1: First, perform a fast Fourier transform on the original time - domain water hammer signal with noise collected at the wellhead of the hydraulic fracturing site, analyze its spectral main peak frequency f 0 , and taking the bandwidth adaptively selected on both sides of the main peak as a reference, adaptively set the upper and lower cut - off frequencies f low and f high, to match the signal characteristics under different well conditions and improve the robustness of subsequent analysis; 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; 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; 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; 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; 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 ); 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.

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

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

[0013] 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: First, calculate the normalized autocorrelation function of the time-domain filtered signal 𝑥[𝑛]: 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.

[0014] Then, the non-negative delayed part of the autocorrelation function is intercepted to form the vector 𝑅 + [𝑘], with 0.5max(𝑅 +) Using the dynamic threshold and taking the number of samples corresponding to the signal sampling frequency 𝐹𝑠 as the minimum peak distance, call the peak search function (such as the findpeaks function in MATLAB) to locate the first two significant peaks 𝑘1 and 𝑘2 on 𝑅 + [𝑘]. The difference in sample points corresponding to the two peaks is Δ𝑘 = 𝑘2 - 𝑘1, and finally dividing by 𝐹𝑠 gives the period T: where T is the signal period, Δ𝑘 is the difference in sample points corresponding to the two peaks, and 𝐹𝑠 is the signal sampling frequency.

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

[0016] Step four is to further perform a de-mean operation on the signal after time-domain filtering in step two. 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 false peaks in the subsequent short-time Fourier transform analysis, thereby improving the resolution and accuracy of the frequency-domain analysis.

[0017] Based on the water hammer signal period T extracted in step three, step five adaptively sets 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 covers at least 2 periods or more. The window function win type uses the kaiser window, which allows for a flexible trade-off between the main lobe width and sidelobe suppression, and its adjustable shape parameter βIt enables the highest frequency resolution while significantly reducing sidelobe leakage, and 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 FFT points nfft represents the number of transformation points used for discrete Fourier transform of each window, that is, zero-padding operation is performed on 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 the next power of 2 greater than the window length is taken. Such parameter selection not only ensures sufficient frequency resolution and anti-aliasing ability, but also takes into account time domain continuity and computational efficiency, laying a solid foundation for subsequent cepstrum matrix construction. In specific implementation, the parameters of the short-time Fourier transform mainly include: input signal data, signal sampling frequency Fs, window function win, window overlap length overlap, number of FFT points nfft. Using these parameters to call the stft function of MATLAB to perform short-time Fourier transform to obtain three parameters: two-dimensional frequency domain matrix S, frequency axis f, and time axis t, providing a basis for subsequent cepstrum denoising and matrix construction. Specifically: Call the stft function of MATLAB using these parameters as follows: [S, f, t] = stft(data, Fs,... 'Window', win,... 'OverlapLength', overlap,... 'FFTLength', nfft); The construction method of the Gaussian band-pass filter described in Step 6 is as follows: First, define the frequency coordinate f as: where 𝐿 is the spectrum length and 𝐹𝑠 is the signal sampling frequency.

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

[0019] Construct a Gaussian low-pass filter according to the standard deviation obtained in the above steps H low ( f ) is: Construct a Gaussian high-pass filter H high ( f ) is: where θ is the coefficient for adjusting the attenuation degree of the high-pass filter.

[0020] The final 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 ): The cepstrum matrix construction process is as follows: for each column of matrix S S j Calculate the magnitude spectrum | S j |, and use the Gaussian band-pass filter with the upper and lower cut-off frequencies set in Step 1 H ( f ) to perform frequency-domain filtering on it; take the logarithm of the filtered magnitude spectrum, then perform the inverse Fourier transform and take its real part to obtain the cepstrum vector 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: where the row index i represents the time delay on the quefrency axis, and the column index j represents the time window position corresponding to the original signal on the time axis, S j represents a certain column of 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 inverse Fourier transform result, and IFFT(⋅) represents the inverse Fourier transform operation, and ln(⋅) represents the operation of taking the natural logarithm.

[0021] In Step 7, the bridge plug corresponds to a larger cepstrum positive peak, and the liquid inlet point corresponds to a smaller cepstrum negative 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: where 𝜌 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: where R denotes the reflection coefficient, representing 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.

[0022] When the medium impedance changes from to (such as a sudden decrease in the cross-sectional area at the bridge plug), 𝑅 > 0, which appears as a positive peak in the cepstrum matrix; while when the impedance changes 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 map 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 maps of the bridge plug and the signals of two liquid inlet points. At the bridge plug position, due to the increase in impedance (sudden decrease in cross-sectional area), positive reflection occurs, and it appears as a positive peak in the cepstrum; while at 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.

[0023] Specifically: First, sum the rows of the cepstrum matrix, and denote the total energy i of the E ( i ) row as: Among them, C ( i,j ) is the value at the i -th row and the j -th column of the cepstrum matrix, N is the number of columns.

[0024] By finding the maximum value of all E ( i ), the row index r b with the maximum energy is determined, corresponding to the position of the packer reflection: Among them, argmax(⋅) is the operation of finding the maximum value.

[0025] The algorithm formula for calculating the relative distance between the packer and the liquid inlet point is as follows: Among them, z r is the actual depth corresponding to the liquid inlet point, q plug is the known depth of the packer, Δτ 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 packer is located, r is the row index where the liquid inlet point is located. The wave speed v can be derived from the packer position and the signal period T, and its calculation formula is: Among them, v is the wave speed, q plug is the known depth of the packer, and T is the water hammer signal period.

[0026] 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 ), a contour map is made, and the precise mapping between each row of the cepstrum cloud map and the actual depth can be obtained. The method of making the cepstrum cloud map used in Step 7 is implemented through the contourf method of MATLAB.

[0027] To further explain the above technical solutions, a specific example 1 is provided in the specific embodiments 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 length of the fracturing sections 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 hydraulic fracturing water hammer signal of the 15th 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: Table 1 Table 1 lists the detailed parameters of the 15th fracturing section used for analysis in the embodiments of the present invention, including data such as the bridge plug position, the top depth and bottom depth of each perforated 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 a total of 12 layers of perforations, numbered from layer 37 to layer 48. The perforation length of each layer is 0.45 meters, and the total number of perforation holes is 84.

[0028] The original waveform diagram of the hydraulic fracturing water hammer signal collected in the 15th section is as shown in Figure 5 Figure. First, perform a fast Fourier transform on the original signal of the hydraulic fracturing water hammer signal collected on site in the 15th section, 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 shown in the spectrum distribution diagram of the water hammer signal, the main peak frequency of the signal is approximately around 0.08 Hz. However, considering that high-frequency signals are often accompanied in the bridge plug signal, the actually selected upper and lower cut-off frequencies Figure 6 of the Gaussian band-pass filter 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, a water hammer signal segment with the strongest energy after pump shutdown is further intercepted for analysis, as shown in Figure 7 Figure.

[0029] Subsequently, apply an alpha-trimmed filter 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 values is removed in each window, that is, the removal ratio α is 2 / 7, and the effective samples in the middle are retained for averaging to improve the robustness to local noise while maximizing the retention of the true waveform changes of the signal. Figure 8The comparison 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 noise is significantly reduced, and the overall waveform of the signal is clearer.

[0030] After completing the time-domain filtering, the autocorrelation algorithm is used 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. Call the peak search function (select the findpeaks function in MATLAB) to locate the first two significant peaks 𝑘1 and 𝑘2 on 𝑅 + [𝑘]. The difference in the corresponding sample points between the two, Δ𝑘 = 𝑘2− 𝑘1, combined with the known signal sampling frequency 𝐹𝑠, is used to calculate the signal period T: 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 detected main signal peaks, and the dashed line is the period demarcation line defined based on the main peak spacing; the lower figure 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 curve is the waveform of the processed signal, and the dashed lines indicate the "theoretically" start / end positions of each period (obtained by averaging the autocorrelation period), which are used to divide the analysis window and do not need to be precisely aligned with each peak or trough. In this way, the original signal can be divided into several complete periodic units, providing a stable periodic reference for the subsequent short-time Fourier transform parameter setting and cepstrum analysis.

[0031] After completing the cycle extraction, a de-mean operation is performed on the time-domain filtered signal, that is, the mean value is subtracted from the entire signal sequence. Subsequently, based on the previously extracted period T ≈ 11.71 s, 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. In this example, the number of sample points corresponding to the period of 11.71 s 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 that while the frequency resolution is the highest, the sidelobe leakage is also significantly reduced. In order to retain more accurate 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.

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

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

[0034] 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 as Combined with the difference between this index and the positions of other negative peaks on the quefrency axis, multiply it 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 based on this, the precise positioning of the liquid inlet point can be completed. Finally, use the contourf plotting function of MATLAB to draw the cepstrum cloud diagram to realize the visual display of the positions of the bridge plug and the liquid inlet points.

[0035] As Figure 11 shown, the cepstrum cloud diagram intuitively shows the distribution of the bridge plug and multiple liquid inlet points in the depth direction. The position of the positive peak of the cepstrum 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 during 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 diagram are clearly distinguishable, the positive and negative peaks are clearly separated, and there is no obvious artifact interference 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.

[0036] 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, 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 location of a fluid inlet point using a hydraulic fracturing water hammer signal, characterized in that: The method comprises the following steps: First, perform 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 cutoff frequencies of the Gaussian bandpass filter based on it; Second, the original noisy water hammer signal is processed in the time domain using an alpha-trimmed filter to remove abnormal spike noise and retain the effective water hammer signal; 3. The period T of the water hammer signal is extracted by using the autocorrelation algorithm on the signal after time domain filtering; 4. Subtract the mean value of the signal filtered in step 2 to eliminate the DC component; 5. Determine the relevant parameters of short-time Fourier transform according to the water hammer signal period T obtained in step 3, and perform short-time Fourier transform operation on the signal window to obtain a two-dimensional matrix S in the frequency domain; 6. Calculate the amplitude spectrum for each column of the matrix S, apply the Gaussian bandpass filter constructed by the upper and lower cutoff frequencies determined in step 1 to filter it, then take the logarithm of the filtered result and perform inverse Fourier transform, take the real part to obtain the cepstrum of each window, and splice the cepstrum of all windows by column to form a cepstrum matrix; 7. Sum the inverse spectrum matrix in step 6 row by row, select the row with the largest sum as the bridge plug position, and 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-to-depth conversion, so as to determine the actual position of each liquid inlet point and draw the inverse spectrum cloud diagram.

2. The method for determining the location of the liquid inlet point using the hydraulic fracturing water hammer signal according to claim 1, characterized in that: The period T described in step three is extracted through an 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 location of the liquid inlet point using the hydraulic fracturing water hammer signal according to claim 1, characterized in that: The short-time Fourier transform parameters described in step five are selected adaptively, including: determining the window length according to the water hammer signal period T extracted in step three, and adaptively setting the window function type, window overlap length and Fourier transform points.

4. The method for determining the location of the liquid inlet point using the hydraulic fracturing water hammer signal according to claim 1, characterized in that: The cepstrum matrix described in step six is ​​constructed in the following manner: the amplitude spectrum is calculated for each column of the matrix S, and it is filtered using a frequency domain Gaussian bandpass filter, where the frequency domain Gaussian bandpass filter is the product of a frequency domain Gaussian low-pass and high-pass filter; the logarithm of the filtered amplitude spectrum is taken, and then an inverse Fourier transform is performed and the real part is taken to obtain the cepstrum of each window; all cepstrum vectors are concatenated in sequence by column to form a two-dimensional cepstrum matrix.

5. The method for determining the location of the liquid inlet point using the hydraulic fracturing water hammer signal according to claim 1, characterized in that: The method for distinguishing the bridge plug from the liquid inlet point in step seven is based on the difference in the signs of the cepstrum values, where the bridge plug corresponds to a larger cepstrum positive peak and the liquid inlet point corresponds to a smaller cepstrum negative peak; by summing the cepstrum matrix row by row, the row index with the largest total energy is determined as the position of the bridge plug in the cepstrum, and the actual depth of each liquid inlet point is calculated based on the difference between it and other row indices multiplied by the sampling time interval and wave velocity.

Images