A method for identifying gas-liquid interface using velocity spectrum of distributed acoustic wave sensing data

Abstract

The present invention proposes a method for identifying gas-liquid interface by using velocity spectrum of distributed acoustic wave sensor data, which belongs to the technical field of oil and gas well exploration; according to the mirror boundary condition, boundary artifacts on the input data are eliminated; first-order gradient filtering is applied to two-dimensional data, followed by frequency and wavenumber filtering; automatic gain control is applied to enhance weak signals in the data; speed-divided energy is extracted; according to the characteristics of F-K domain data, a two-dimensional frequency matrix and a two-dimensional wavenumber matrix are created, and a velocity spectrum is created by the frequency-wavenumber ratio of each corresponding unit in the two matrices; a velocity range is given to create a velocity shielding matrix; an inverse two-dimensional fast Fourier transform is applied to present the data in the time-space domain; and the autocorrelation coefficient is calculated based on the speed-divided energy data. The present invention can track the changes of the gas-liquid interface in real time and provide dynamic wellbore flow information.

Description

Technical Field

[0001] The invention provides a method for identifying a gas-liquid interface by utilizing a distributed acoustic wave sensing data velocity spectrum, and belongs to the technical field of oil and gas well exploration. Background Art

[0002] For the identification of gas-liquid interface in oil and gas wells, downhole multiphase flowmeters or pressure gauges can only provide information at discrete locations in the wellbore and cannot cover the entire wellbore interval at the same time. The use of logging technology for interface tracking and detection is costly, and production must be interrupted during logging, resulting in production losses. In recent years, it has become increasingly common to deploy optical fibers in oil and gas wells to analyze fluid flow characteristics using distributed acoustic sensing (DAS) technology. Compared with production logging based on traditional point sensors, a major advantage of distributed fiber optic sensors is that they can provide continuous information in space and time along the entire length of the wellbore or optical fiber without the need to install any additional electronic equipment. DAS can detect vibrations caused by fluid flow in the wellbore and changes in the speed of sound (SOS) in multiphase flow conditions. Summary of the invention

[0003] To overcome the above limitations and challenges, the present invention proposes a method for identifying the gas-liquid interface using the velocity spectrum of distributed acoustic sensor data. The method uses the proposed velocity-dependent energy extraction process to process the vibration data collected by DAS, and the velocity spectrum extracted from the DAS data is used to identify the gas-liquid interface.

[0004] The present invention considers the sound velocity of a certain phase fluid in a multiphase flow to extract the acoustic energy or signal unique to the phase fluid. This helps to estimate the flow characteristics of the phase fluid. Since the amount of DAS data may be very large (wellbore monitoring can reach several TB per day), the gas-liquid interface identification method based on velocity-dependent energy also has the practical advantage of compressing data without losing valuable information. Since DAS can provide spatially continuous information along the entire length of the optical fiber, this method is not limited by the flow state and is suitable for upward or downward vertical flow or inclined flow, as well as small and large diameter pipes. The method takes into account the sound velocity of a certain phase fluid to filter out the acoustic energy or signal contributed by the phase in the multiphase medium. The gas-liquid interface can be identified by processing DAS data through the workflow of the gas-liquid interface identification method based on velocity-dependent energy.

[0005] The specific technical solution is:

[0006] A method for identifying a gas-liquid interface using a velocity spectrum of distributed acoustic wave sensing data comprises the following steps:

[0007] (1) Eliminate boundary artifacts on the input data based on mirror boundary conditions.

[0008] The use of mirror boundary conditions is mainly used to deal with edge effects of finite-length data sets. The original data set is extended by creating virtual mirror data at the data boundary. For acquired 2D DAS data, mirror boundary conditions are applied to the depth dimension.

[0009] Specifically, mirrored data is added to the top and bottom boundaries of the data. At the top of the data: the first 25% of the original data is flipped upside down and added to the top of the data. At the bottom of the data: the last 25% of the original data is flipped upside down and added to the bottom of the data.

[0010] (2) Applying a first-order gradient filter to the two-dimensional data to remove the horizontal noise observed in the data, followed by frequency and wavenumber filtering. In the present invention, the FK filter calculation is to perform frequency filtering and wavenumber filtering in sequence. For frequency filtering, a high-pass Butterworth filter is applied in the time direction to remove DC noise (i.e., at zero frequency) and submersible pump noise in the signal.

[0011] (3) After FK filtering, automatic gain control is applied to enhance weak signals in the data; the automatic gain control calculation involves multiplying the DAS data by a scaling factor h(x), which is estimated in equation (1) as the inverse of the root mean square of all DAS data at that depth in each considered time window, expressed as:

[0012]

[0013] Where K i is the DAS data for a given row in the two-dimensional matrix, and M is the total number of samples at each depth within the recording time window.

[0014] (4) Perform velocity-dependent energy extraction, including phase velocity filtering. The phase velocity filtering step uses the AGC filtered data as input. Then, a two-dimensional fast Fourier transform is applied, as shown in formula (2), to convert the AGC processed data into the frequency-wavenumber (FK) domain.

[0015] G(ω,λ)=∫∫f(τ,ξ)e -(λξ-ωτ) dξdτ (2)

[0016] Where ω represents the angular frequency, λ represents the angular wave number, τ represents the time, and ξ represents the spatial coordinate.

[0017] (5) According to the characteristics of the FK domain data, a two-dimensional frequency matrix and a two-dimensional wavenumber matrix are created, and a velocity spectrum is created by the frequency-wavenumber ratio of each corresponding unit in the two matrices. Then, a velocity range is given to create a velocity shielding matrix. The selected velocity range depends on the target fluid. If the focus is on a mixture of nitrogen and water, in order to detect the gas-water two-phase region, a velocity range of 270 to 360 m / s is specified based on the sound speed in the nitrogen-water mixture. In order to detect the single-phase water region, a velocity range of 1500 to 1600 m / s is specified. The FK domain data is multiplied by the velocity shielding data to obtain velocity filtered data (in the FK space), that is, the velocity-divided energy data.

[0018] (6) Apply inverse two-dimensional fast Fourier transform to present the data in the time-space domain. Then calculate the autocorrelation coefficient based on the speed-energy data. The autocorrelation function is used to obtain the similarity between the time series and itself, which helps to obtain the periodic information embedded in the data. The autocorrelation function of the analysis function s(t) is as shown in formula (3):

[0019]

[0020] Where s(t) represents the signal in the time domain and δ represents the time delay.

[0021] Calculating the autocorrelation of a signal s(t) is equivalent to calculating the square of the amplitude spectrum in the frequency domain. In other words, calculating the power spectrum of x(t) in the frequency domain is another way to calculate the autocorrelation function, as shown in formula (4)

[0022]

[0023] Where A(δ) is the autocorrelation function in formula 3, S(f) is the Fourier transform of the signal s(t), and f is the angular frequency.

[0024] The present invention screens out the acoustic energy or signal contributed by a specific phase according to the sound velocity range of the phase, which is called "differential velocity energy". The differential velocity energy is then used for characteristic analysis of multiphase flow. The present invention demonstrates the application of the proposed differential velocity energy workflow in real-time tracking of the gas-liquid interface along the entire length of the wellbore. The differential velocity energy technology uses a series of signal processing and adjustment steps to reduce noise and enhance the target signal. The results obtained from this new method will further help to verify DAS-based flow monitoring algorithms, leak detection and quantification, and reservoir characteristic analysis.

[0025] The technical effects of the technical solution of the present invention are as follows:

[0026] The present invention can track the changes in the gas-liquid interface in real time and provide dynamic wellbore flow information. It helps to timely discover and respond to downhole flow problems. By accurately identifying the gas-liquid interface, production parameters such as pump speed and gas injection volume can be optimized.

[0027] The entire wellbore can be continuously monitored, not just at discrete points. A more comprehensive wellbore flow image is provided. Abnormal gas influx can be detected in time to prevent safety accidents such as blowouts. BRIEF DESCRIPTION OF THE DRAWINGS

[0028] Figure 1 A schematic diagram of the optical fiber layout for collecting data in an embodiment;

[0029] Figure 2 is a flow chart of the present invention;

[0030] Figure 3 is a normalized data graph of the embodiment;

[0031] In the figure, the gas-water mixing curve is mainly on the left side, with the peak on the right side, and the single-phase water curve is mainly on the right side, with the peak on the left side. DETAILED DESCRIPTION

[0032] In this embodiment, the data acquisition optical fiber is coupled outside the oil pipe, and the optical fiber is connected to the DAS demodulator, and the collected data enters the data processing flow method proposed in the present invention. The schematic diagram of the optical fiber layout is shown in FIG. Figure 1 shown.

[0033] The technical roadmap of the present invention is as follows Figure 2 As shown, from receiving the data collected by the optical fiber to the final calculation of the position of the gas-liquid interface.

[0034] The specific steps include:

[0035] (1) The DAS data collected by optical fiber is two-dimensional data, which exists in two dimensions: time and depth. First, the boundary artifacts on the input data are eliminated according to the mirror boundary condition. In image or two-dimensional data analysis, if the boundary value problem is not handled properly, it will usually cause ring or circular artifacts in the processed data. Here, the use of mirror boundary conditions is mainly used to deal with the edge effects of finite-length data sets. It extends the original data set by creating virtual mirror data at the data boundary. For the collected two-dimensional DAS data (time-depth two-dimensional data), the mirror boundary condition is applied to the depth dimension. Specifically, mirror data is added to the top and bottom boundaries of the data. At the top of the data: the first 25% of the original data is flipped upside down and added to the top of the data. At the bottom of the data: the last 25% of the original data is flipped upside down and added to the bottom of the data. The purpose of using the mirror boundary method is to reduce edge effects. By adding mirror data, the distortion that may occur when performing frequency-wavenumber (FK) transformation at the data boundary can be reduced.

[0036] (2) Next, a first order gradient filter is applied to the two-dimensional data to remove the horizontal noise observed in the data used in the present invention, followed by frequency and wavenumber filtering. In the present invention, the FK filter is calculated as frequency filtering and wavenumber filtering in sequence. For frequency filtering, a high pass Butterworth filter is applied in the time direction to remove DC noise (i.e., zero frequency) and submersible pump noise from the signal. The purpose of wavenumber filtering is to improve the signal-to-noise ratio of a specific signal frequency, and a Butterworth filter is applied in the depth direction to ensure that long wavelengths are removed from the data to improve the signal-to-noise ratio of the expansion bubble signal in the mixed fluid.

[0037] (3) After FK filtering, automatic gain control is applied to enhance weak signals in the data to improve the signal-to-noise ratio. Weak signals are introduced by bubbles, which are more attenuating than flowing water. The automatic gain control calculation involves multiplying the DAS data by a scaling factor h(x), which is estimated in equation (1) as the inverse of the root mean square of all DAS data at that depth in each considered time window, expressed as:

[0038]

[0039] Where K i is the DAS data for a given row in the two-dimensional matrix, and M is the total number of samples at each depth within the recording time window.

[0040] (4) The next step is to extract the energy by velocity. The key step is phase velocity filtering. The phase velocity filtering step takes the AGC filtered data as input. Then a two-dimensional fast Fourier transform is applied, as shown in formula (2), to convert the AGC processed data into the frequency-wavenumber (FK) domain.

[0041] G(ω,λ)=∫∫f(τ,ξ)e -(λξ-ωτ) dξdτ (2)

[0042] Where ω represents the angular frequency, λ represents the angular wave number, τ represents the time, and ξ represents the spatial coordinate.

[0043] (5) According to the characteristics of the FK domain data, a two-dimensional frequency matrix and a two-dimensional wavenumber matrix are created, and a velocity spectrum is created by the frequency-wavenumber ratio of each corresponding unit in the two matrices. Then, a velocity range is given to create a velocity shielding matrix. The selected velocity range depends on the target fluid. If the focus is on a mixture of nitrogen and water, in order to detect the gas-water two-phase region, based on the sound speed in the nitrogen-water mixture, a velocity range of 270 to 360 m / s can be specified. Similarly, in order to detect the single-phase water region, a velocity range of 1500 to 1600 m / s can be specified. The FK domain data is multiplied by the velocity shielding data to obtain velocity filtered data (in the FK space), that is, the velocity-divided energy data.

[0044] (6) Finally, the inverse two-dimensional fast Fourier transform is applied to present the data in the time-space domain. The autocorrelation coefficient is then calculated based on the speed-energy data. The autocorrelation function is used to obtain the similarity between the time series and itself, which helps to obtain the periodic information embedded in the data. The autocorrelation function of the analysis function s(t) is shown in formula (3):

[0045]

[0046] Where s(t) represents the signal in the time domain and δ represents the time delay.

[0047] Calculating the autocorrelation of a signal s(t) is equivalent to calculating the square of the amplitude spectrum in the frequency domain. In other words, calculating the power spectrum of x(t) in the frequency domain is another way to calculate the autocorrelation function, as shown in formula (4)

[0048]

[0049] Where A(δ) is the autocorrelation function in formula 3, S(f) is the Fourier transform of the signal s(t), and f is the angular frequency.

[0050] In the present invention, the autocorrelation is calculated for each depth in the frequency domain, and then the inverse Fourier transform is applied to convert the result back to the time-space domain, and then superimposed, and the root mean square amplitude is calculated. The root mean square amplitude value is then normalized using the minimum-maximum method to enhance phase detection. According to the calculated normalized root mean square amplitude data curve, a high value of the root mean square amplitude of water indicates the presence of single-phase water, while a high value of the root mean square amplitude of gas-water indicates the presence of a gas-water (two-phase) region. Based on the root mean square amplitude corresponding to water and gas-water mixtures, the "half-maximum full width" can be used to infer the interface between the single-phase water and gas-water mixture regions and the height of the two-phase gas-water region. The half-maximum full width value is mainly used to describe the peak value or width. Simply put, it refers to the width measured at half the height of the maximum value on a peak or distribution curve. This width is the height of the gas-water two-phase, and the two height surfaces are at the gas-liquid interface.

[0051] In this embodiment, downhole pressure sensors and temperature sensors are set in a specific test well. The test well is a cased well with a total length of 1573m and a depth, where the injection gas (radon gas) test is conducted. The test well is equipped with four downhole pressure sensors and temperature sensors, located at depths of 147m, 616m, 1067m and 1531m respectively. An armored optical fiber is fixed to the outside of the oil pipe by a steel clamp.

[0052] In the test well described, the wellbore (including tubing and annulus) was initially filled with water. A fixed volume of nitrogen was then injected through the tubing at a constant rate, with the surface pressure set at 1.4 MPa. During the gas injection process, the amount of water displaced by the gas was recorded once a minute. When the cumulative return volume reached a predetermined volume, the gas injection was stopped. Subsequently, injection water was injected through the tubing at a rate of 1 cubic meter per minute to displace the gas to the bottom of the well. Once the gas plug reached the bottom of the tubing, the circulation was stopped and the gas was allowed to rise to the surface through the annulus. During the water injection process, pressure, temperature, and DAS data were recorded in real time. It was injected through the tubing and then bubbled up to the surface through the annulus as a two-phase (nitrogen-water) mixture.

[0053] The collected data are processed by the method of the present invention, and two normalized RMS data are calculated, such as Figure 3 shown.

[0054] To verify the accuracy of data processing, the pressure difference between 616m and 1067m was analyzed. The moment when the pressure difference began to decrease (t1) indicated that the gas-water interface first crossed the sensor position. Finally, the pressure difference stabilized (at time t2), indicating that the gas-water mixing area had completely passed the sensor at 1067m.

[0055] To estimate the interface in the gas-water mixing zone, we take the time difference (t2 minus t1) and multiply it by the average gas rise velocity estimated experimentally and verified by numerical simulations, which in this test is 0.15 m / s. This gives an average column height of 235 m, which agrees well with the estimate using the workflow in this paper.

Claims

1. A method for identifying a gas-liquid interface using a velocity spectrum of distributed acoustic wave sensing data, characterized in that: The following steps are involved: (1) Eliminate boundary artifacts on the input data according to the mirror boundary condition; the mirror boundary condition is mainly used to deal with the edge effect of a finite length data set; the original data set is extended by creating virtual mirror data at the data boundary; for the acquired two-dimensional DAS data, the mirror boundary condition is applied to the depth dimension; (2) applying a first-order gradient filter to the 2D data to remove the horizontal noise observed in the data, followed by frequency and wavenumber filtering; (3) After filtering, automatic gain control is applied to enhance weak signals in the data; (4) Perform velocity-dependent energy extraction, including phase velocity filtering; (5) According to the characteristics of FK domain data, a two-dimensional frequency matrix and a two-dimensional wave number matrix are created, and a velocity spectrum is created by calculating the frequency and wave number ratio of each corresponding unit in the two matrices; Then, a velocity range is given to create a velocity masking matrix; (6) Applying inverse two-dimensional fast Fourier transform to present the data in the time-space domain; The autocorrelation coefficient is then calculated based on the velocity-energy data. The velocity-energy data is used to analyze the characteristics of multiphase flow, track the changes in the gas-liquid interface in real time, and provide dynamic wellbore flow information.

2. The method for identifying the gas-liquid interface using the velocity spectrum of distributed acoustic wave sensing data according to claim 1, characterized in that: The specific method of step (1) is as follows: mirror data is added to the top and bottom boundaries of the data; at the top of the data: the first 25% of the original data is flipped upside down and added to the top of the data; at the bottom of the data: the last 25% of the original data is flipped upside down and added to the bottom of the data.

3. The method for identifying the gas-liquid interface using the velocity spectrum of distributed acoustic wave sensing data according to claim 1, characterized in that: The specific method of step (2) is: The FK filtering calculation is to perform frequency filtering and wave number filtering in sequence; for frequency filtering, a high-pass Butterworth filter is applied in the time direction to remove DC noise and submersible pump noise in the signal.

4. The method for identifying a gas-liquid interface using a velocity spectrum of distributed acoustic wave sensing data according to claim 1, characterized in that: The specific method of step (3) is as follows: the automatic gain control calculation multiplies the DAS data by a scaling factor h(x), which is estimated in formula (1) as the inverse of the root mean square of all DAS data at that depth in each considered time window, expressed as: Where K i is the DAS data for a given row in the two-dimensional matrix, and M is the total number of samples at each depth within the recording time window.

5. The method for identifying the gas-liquid interface using the velocity spectrum of distributed acoustic wave sensing data according to claim 1, characterized in that: The phase velocity filtering step in step (4) uses the AGC filtered data as input; then a two-dimensional fast Fourier transform is applied, as shown in formula (2), to convert the AGC processed data into the frequency-wavenumber (FK) domain; G(ω,λ)=∫∫f(τ,ξ)e -(λξ-ωτ) dξdτ (2) Where ω represents the angular frequency, λ represents the angular wave number, τ represents the time, and ξ represents the spatial coordinate.

6. The method for identifying a gas-liquid interface using a velocity spectrum of distributed acoustic wave sensing data according to claim 1, characterized in that: In step (5), the selected velocity range depends on the target fluid; if the focus is on a mixture of nitrogen and water, in order to detect the gas-water two-phase region, a velocity range of 270 to 360 m / s is specified based on the speed of sound in the nitrogen-water mixture; in order to detect the single-phase water region, a velocity range of 1500 to 1600 m / s is specified; the FK domain data is multiplied by the velocity shielding data to obtain velocity filtering data, i.e., velocity-divided energy data.

7. The method for identifying a gas-liquid interface using a velocity spectrum of distributed acoustic wave sensing data according to claim 1, characterized in that: In step (6), the autocorrelation function is used to obtain the similarity between the time series and itself. The autocorrelation function of the analysis function s(t) is as shown in formula (3): Where s(t) represents the signal in the time domain and δ represents the time delay; Calculating the autocorrelation of a signal s(t) is equivalent to calculating the square of the amplitude spectrum in the frequency domain. Calculating the power spectrum of x(t) in the frequency domain is another way to calculate the autocorrelation function, as shown in formula (4): Where A(δ) is the autocorrelation function in formula 3, S(f) is the Fourier transform of the signal s(t), and f is the angular frequency.

Images