A method for separating multi-mode signals of casing well ultrasonic bending wave logging

Abstract

The application discloses a casing well ultrasonic bending wave logging multi-mode signal separation method, which comprises the following steps: inputting original mixed casing well ultrasonic bending wave logging data and a related parameter file; simulating A0 mode waves in the instrument centering condition; adopting a variational mode decomposition algorithm to decompose the original mixed signal into a set of multiple scale intrinsic mode functions; taking the correlation coefficient of the reconstructed A0 mode wave and the simulated standard A0 mode wave as an inversion objective function; adopting a simulated annealing method to optimize the optimal decomposition layer number and the reconstruction vector in the inversion objective function; and calculating a reconstructed A0 mode wave A' based on the optimal decomposition layer number k and the reconstruction vector v obtained by optimizing the inversion objective function. o The application can realize adaptive decomposition of A0 and S0 mode waves of a single detector in a non-axisymmetric casing well measuring environment, and can provide accurate A0 mode waves for a method for calculating wave impedance of a post-casing medium based on A0 first wave attenuation information.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of cementing quality evaluation, in particular to a method for separating multi-modal signals of casing well ultrasonic bending wave logging. BACKGROUND

[0002] Cementing quality evaluation is an essential part in the fields of oil and gas and geothermal resource exploitation, gas storage construction and carbon dioxide geological storage. Imaging of annulus medium outside the casing can quantitatively evaluate the cementing quality of casing-cement interface (I interface) and cement-formation interface (II interface), estimate the thickness of cement annulus medium, judge the integrity of wellbore, realize the sealing of oil, gas and water layers, prevent the channeling of interlayer fluid and combustible gas, and has important significance for accurate determination of perforated layer, improvement of oil and gas well exploitation efficiency and service life, and protection of ecological environment and exploitation safety.

[0003] Ultrasonic bending wave logging uses a specific angle of incidence to excite zero-order antisymmetric mode Lamb wave (A0 mode wave) propagating along the well axis direction in the casing. A0 primary wave attenuation and ultrasonic pulse echo calculation of post-casing medium acoustic impedance intersection analysis can effectively identify solid-liquid-gas three-phase medium. A0 mode wave leaks part of the energy to the annulus medium during propagation and is reflected back to the well by the formation interface, forming the third acoustic interface reflection echo (TIE) carrying the annulus medium-formation interface information. Analysis of the first arrival, dispersion and amplitude characteristics of TIE reflected wave can realize quantitative evaluation of II interface cementing quality.

[0004] However, in the measurement environment of high-deviation wells or horizontal wells, affected by factors such as instrument inclination, instrument and casing eccentricity, and change of incidence angle, the acoustic field excited by ultrasonic bending wave logging in the non-axisymmetric casing well environment contains A0 mode wave and S0 mode wave (zero-order symmetric mode Lamb wave) with different amplitudes, different propagation speeds but similar frequency components. The two kinds of mode waves are mixed together in the time domain and are difficult to distinguish, making it difficult for the existing method to estimate accurate A0 mode wave attenuation and introducing large errors in the quantitative evaluation of casing-cement interface cementing quality. Therefore, it is urgent to establish a casing well ultrasonic bending wave logging wave field separation method to suppress S0 mode wave and extract accurate A0 mode wave.

[0005] At present, the research on the casing well ultrasonic bending wave logging A0 mode wave and S0 mode wave separation method for eccentric sound field is less. In addition, the ultrasonic bending wave logging instrument only contains two detectors, far (source distance 0.35m) and near (source distance 0.25m), so that the conventional wellbore sound field separation algorithm based on array signal such as median filter (Wang Bing et al., 2011. Journal of China University of Petroleum, 35(2):57-63.), F-K filter (Liu Hao et al., 2021. Petroleum and Chemical Industry Application, 40(7).), Radon transform (Li Chao et al., 2014. Progress in Geophysics, 29(4):1678-1688.) and other methods cannot be applied to the ultrasonic bending wave logging instrument data with limited spatial sampling information, and it is difficult to realize the separation of A0 and S0 mode waves of a single detector. SUMMARY

[0006] In view of the above problems existing in the prior art, the purpose of the present application is to provide a casing well ultrasonic bending wave logging multi-modal signal separation method.

[0007] In order to achieve the above-mentioned target, the present application adopts the following technical solutions:

[0008] A casing well ultrasonic bending wave logging multi-modal signal separation method, comprising the following steps:

[0009] Step1: input the original mixed casing well ultrasonic bending wave logging waveform data, casing well structure and the elastic parameters of casing, fluid and cement;

[0010] Step2: simulate A0 mode wave according to the casing well structure and the elastic parameters of the current casing well;

[0011] Step3: use the variational mode decomposition algorithm to decompose the original mixed signal into a set of multiple scale intrinsic mode functions;

[0012] Step4: based on the correlation coefficient of the reconstructed A0 mode wave and the simulated standard A0 mode wave, construct an inversion objective function;

[0013] Step5: use the simulated annealing method to optimize the optimal decomposition layer number and the reconstruction vector in the inversion objective function;

[0014] Step6: based on the optimal decomposition layer number k o and the reconstruction vector v obtained by optimizing the inversion objective function, calculate the reconstructed A0 mode wave A'.

[0015] Further, in Step1, the maximum decomposition layer number of variational mode decomposition is set.

[0016] Generally, the maximum decomposition layers are set to 5-8 layers, and too small decomposition layers cannot obtain ideal reconstruction results, and too large decomposition layers will introduce unnecessary calculation amount for subsequent inversion of the objective function.

[0017] Further, in Step 2, according to the well diameter, casing thickness and elastic information of casing, well fluid and cement, a high-order staggered grid finite difference method is used to simulate the A0 mode wave under the condition of instrument strict centering as a standard A0 wave form.

[0018] Further, in Step 3, the original signal is decomposed into a set of multiple scale intrinsic functions with compact support by the variational mode decomposition:

[0019]

[0020]

[0021] wherein u k and w k are the intrinsic mode function of the kth scale and its corresponding center frequency, δ(t) is the Dirac function, j is a complex number, e is the natural logarithm, π is the circular constant, M is the original aliasing signal, t is the time sampling point, and n is the maximum decomposition layer.

[0022] Further, in Step 4, based on the correlation coefficient of the reconstructed A0 mode wave and the simulated standard A0 mode wave, an inversion objective function is constructed:

[0023]

[0024] wherein k o is the optimal decomposition layer, v is the multi-scale reconstruction vector, A is the A0 signal of instrument strict centering, A' is the A0 signal reconstructed by the variational mode decomposition, Cov is the covariance calculation, and Var is the variance calculation.

[0025] Further, in Step 5, the initial temperature T0 and the maximum temperature drop times n max are needed to be preset for the algorithm, in order to ensure the stability of the optimization process, a linear cooling method is adopted:

[0026] T n+1 = λT n

[0027] wherein T n and T n+1 are the temperature of the nth time and the n+1th time, respectively, and λ is the temperature drop rate. For a certain temperature T n , an isothermal equilibrium method is adopted, and based on the Metropolis criterion, the new solution acceptance probability of the inversion objective function of N cycles at the current temperature is calculated.

[0028]

[0029] where L is the inversion objective function, and v n is the optimal decomposition level and reconstruction vector at the n-th iteration, and v n+1 is the optimal decomposition level and reconstruction vector at the n+1-th iteration, and the new solution is generated randomly in the neighborhood of the previous solution.

[0030] In a preferred embodiment of the present application, the initial temperature is set to 200, and the temperature reduction rate is set to 0.95. The algorithm has two layers of loops, the first layer of loops is linear temperature reduction, and the second layer of loops is to generate a new solution by random disturbance for N times at a certain temperature. v n+1 and calculate the change of the objective function value to determine whether to be accepted. Since the initial temperature of the algorithm is relatively high, the new solution with an increased inversion objective function L may also be accepted at the beginning, so the algorithm has the ability to jump out of the local minimum value. Then by slowly reducing the temperature, the algorithm finally converges to the global optimal solution. Further, in Step 6, the optimal decomposition level k o and the reconstruction vector v calculated based on the inversion objective function optimization are used to calculate the reconstructed A0 signal A':

[0031]

[0032] where n v is the total number of elements of the reconstruction vector v, i is the element index of the reconstruction vector v, k is the scale index, and u k is the intrinsic mode function of the k-th scale.

[0033] Compared with the prior art, the present application has the following beneficial effects

[0034] (1) The method of the present application applies the variational mode decomposition algorithm suitable for multi-scale decomposition of transient signals to non-axisymmetric casing well ultrasonic bending wave logging, can realize multi-scale decomposition of the aliasing signals of A0 and S0 mode waves, can give the optimal decomposition level and reconstruction vector according to the established inversion objective function and optimization algorithm, and can adaptively separate A0 mode waves and S0 mode waves from the aliasing signals of a single geophone to extract accurate A0 mode waves, thereby laying a foundation for the method of calculating the wave impedance of the medium behind the casing based on A0 first wave attenuation information.

[0035] (2) The variational mode decomposition method proposed by the present application does not require array waveform data, reduces the requirements on the acquisition conditions, is a data-driven signal separation method, has fast calculation speed, and meets the needs of rapid evaluation of casing well ultrasonic bending wave cementing quality. BRIEF DESCRIPTION OF DRAWINGS

[0036] Other features, objects, and advantages of the application will become more apparent from the following detailed description of non-limiting embodiments thereof, when read in connection with the following drawings:

[0037] Figure 1 is a flow chart of the method provided by the present application;

[0038] Figure 2 is the simulated standard A0 mode wave (top), S0 mode wave (middle) and the aliasing waveform M0 (bottom);

[0039] Figure 3 is the intrinsic function at each scale obtained by the variational mode decomposition of the aliasing waveform;

[0040] Figure 4 is the similarity coefficient of the intrinsic function at each scale and the standard A0 mode wave;

[0041] Figure 5 is the comparison of the separated A0 mode wave and the simulated standard A0 mode wave and the comparison of the aliasing signal M0 and the simulated standard A0 mode wave;

[0042] Figure 6 is the measured casing well ultrasonic bending wave aliasing waveform under the condition of instrument eccentricity;

[0043] Figure 7 is the intrinsic function at each scale obtained by the variational mode decomposition of the measured aliasing waveform;

[0044] Figure 8 is the A0 mode wave separated from the measured aliasing waveform. DETAILED DESCRIPTION

[0045] The present application will be described in detail below with specific embodiments. The following examples will help those skilled in the art to further understand the present application, but do not limit the present application in any form. It should be noted that for those skilled in the art, without departing from the concept of the present application, a number of modifications and improvements can be made. These all belong to the protection scope of the present application.

[0046] The application applies a variational mode decomposition algorithm with signal adaptivity, shift invariance and compact support to casing well ultrasonic bending wave logging, aims at the problem of A0 mode wave and S0 mode wave aliasing caused by non-axisymmetric acoustic field environment, decomposes the aliasing signal into multiple time-frequency subintervals with limited bandwidth according to the oscillation characteristics thereof, constructs an optimization objective function based on the similarity of the effective scale reconstructed signal and the standard A0 mode wave, obtains the optimal decomposition layer number and the reconstructed vector through optimization of the objective function, and adaptively realizes A0 signal reconstruction of a single detector, suppresses the interference of the S0 mode wave, and provides accurate A0 mode wave for quantitative evaluation of the cementing quality of the casing-cement interface (I interface) and the cement-formation interface (I interface) in the oil exploitation process.

[0047] To solve the problems of the prior art, the purpose of the present application is to provide a casing well ultrasonic bending wave logging multi-modal signal separation method, which uses a multi-scale variational mode decomposition algorithm to adaptively separate A0 mode wave and S0 mode wave signals of a single detector in casing well ultrasonic bending wave logging.

[0048] Embodiment 1: Simulation of ultrasonic bending wave logging data implementation case

[0049] As shown in Figure 1 , a casing well ultrasonic bending wave logging multi-modal signal separation method comprises the following steps:

[0050] Step 1: read in the original aliasing simulated ultrasonic bending wave logging waveform data (SimulatedWave.dat), input the casing well bore diameter, casing size and elastic information of the casing, well fluid and cement, and set the maximum decomposition layer number of the variational mode decomposition to 8 layers.

[0051] Step 2: use a high-order staggered grid finite difference method to simulate the A0 mode wave when the casing well instrument is strictly centered. Figure 2 The simulated standard A0 mode wave, S0 mode wave and the aliasing signal M0 of the two.

[0052] Step 3: use the variational mode decomposition to decompose the original aliasing signal M0 into each scale domain. Taking 8 layers of decomposition as an example, the intrinsic function results of each scale decomposition are as shown in Figure 3 , in which the solid line is the intrinsic function of each scale, and the dashed line is the simulated standard A0 mode wave. The similarity of the intrinsic function of each scale and the simulated standard A0 mode wave is as shown in Figure 4 . As can be seen from the figure, the intrinsic functions of the 7th scale and the 8th scale have higher similarity with the simulated A0 mode wave.

[0053] In Step 3, the variational mode decomposition decomposes the original signal into a set of multiple scale intrinsic functions with compact support:

[0054]

[0055]

[0056] Among them, u k and w k These are the intrinsic mode function and its corresponding center frequency at the k-th scale, respectively. δ(t) is the Dirac function, j is a complex number, e is the natural logarithm, π is pi, M is the original aliased signal, t is the time sampling point, and n is the maximum decomposition level.

[0057] Step 4: Construct the inversion objective function based on the correlation coefficient between the reconstructed A0 mode wave and the simulated A0 mode wave.

[0058]

[0059] Where, k o denoted as the optimal decomposition level, v as the multi-scale reconstruction vector, A as the A0 signal strictly centered by the instrument, A' as the A0 signal reconstructed through variational mode decomposition, Cov as the covariance calculation, and Var as the variance calculation.

[0060] Step 5: Simulated annealing is used to optimize the inversion objective function to obtain the optimal number of decomposition layers and the reconstructed vector. The initial temperature is set to 200°C, and the temperature decrease rate is set to 0.95%.

[0061] In Step 5, the initial temperature T0 and the maximum number of cooling cycles n need to be preset for the algorithm. max To ensure the stability of the optimization process, a linear cooling method is adopted:

[0062] T n+1 =λT n

[0063] Among them, T n and T n+1 Let T be the temperatures at the nth and (n+1)th times, respectively, and λ be the rate of temperature decrease. For a given temperature T... n Using an isothermal equilibrium approach and based on the Metropolis criterion, the probability of receiving a new solution to the inversion objective function after N cycles at the current temperature is calculated.

[0064]

[0065] Where L is the inversion objective function, and v n Let be the optimal decomposition level and reconstruction vector for the nth iteration. and v n+1 Let be the optimal decomposition level and reconstruction vector for the (n+1)th iteration. The new solution is generated randomly within the neighborhood of the previous solution.

[0066] The algorithm consists of two nested loops. The first loop linearly reduces the temperature, while the second loop generates a new solution by performing N random perturbations on a given temperature. v n+1 The algorithm calculates the change in the objective function value and determines whether it is accepted. Since the initial temperature of the algorithm is relatively high, a new solution that increases L may also be accepted initially, thus escaping local minima. Then, by slowly decreasing the temperature, the algorithm eventually converges to the global optimum.

[0067] Step 6: The optimization result of the inversion objective function is: the optimal decomposition level is 8, and the reconstruction vectors are 7 and 8. Therefore, the sum of the intrinsic functions at the 7th and 8th scales is used as the reconstructed A0 mode wave.

[0068] The optimal decomposition level k is obtained based on the optimization of the inversion objective function. o Calculate the reconstructed A0 signal A' from the reconstructed vector v:

[0069]

[0070] Where, n v Let i be the total number of elements in the reconstructed vector v, i be the element index of the reconstructed vector v, k be the scale index, and u be the value of the vector. k It is the intrinsic mode function at the k-th scale.

[0071] Figure 5 Comparison of the reconstructed A0 mode waveform with the standard A0 mode waveform (dashed line) (top); comparison of the original aliased signal M0 with the standard A0 mode waveform (dashed line) (bottom). Figure 5 It can be seen that after multi-scale decomposition and reconstruction, the separated A0 mode wave is very similar to the standard A0 mode wave, which proves the effectiveness of the method of the present invention.

[0072] Example 2: Case Study of Measured Ultrasonic Curved Wave Logging Data

[0073] Step 1: Read in the original aliased simulated ultrasonic bending wave logging waveform data (3_middle sub_log.dat), input the casing well diameter, casing size, and elastic information of the casing, well fluid, and cement, and set the maximum number of variational mode decomposition layers to 8 layers. Figure 6 The image shows the original aliased ultrasonic bending wave logging waveform generated by the instrument's tilt in a high-angle well measurement environment. Figure 6 It can be seen that after the first wave A0, there is a strong overlap between the S0 mode wave and the A0 mode wave.

[0074] Step 2: Since the standard A0 mode wave cannot be obtained from the actual well data, the high-order staggered grid finite difference method is still used to simulate the A0 mode wave of the casing well when the instrument is strictly centered, which is used as the standard A0 wave form for establishing the inversion objective function.

[0075] Step 3: The original mixed signal M0 is decomposed into each scale domain by using the variational mode decomposition. Taking 8-layer decomposition as an example, the intrinsic function results of each scale decomposition are shown in FIG. 1, wherein the solid line is the intrinsic function of each scale, and the dashed line is the simulated A0 mode wave. As shown in the figure, the intrinsic functions of the first scale and the second scale have higher similarity with the simulated A0 mode wave. The third scale and the fourth scale mainly contain S0 component. The fifth scale to the eighth scale mainly contain borehole environmental noise. Figure 7

[0076] Step 4: Based on the correlation coefficient of the reconstructed A0 mode wave and the simulated A0 mode wave, the inversion objective function is constructed.

[0077] Step 5: The simulated annealing method is used to optimize the inversion objective function to obtain the optimal decomposition layer number and the reconstruction vector. The initial temperature is set to 200, and the temperature drop rate is set to 0.95.

[0078] Step 6: The optimization result of the inversion objective function is that the optimal decomposition layer number is 8, and the reconstruction vector is 1 and 2. Therefore, the sum of the intrinsic functions of the first scale and the second scale is used as the reconstructed A0 mode wave. Figure 8 The comparison between the reconstructed A0 mode wave and the original mixed signal M0 is shown in FIG. 2. As shown in the figure, after the multi-scale decomposition and reconstruction, the amplitude of the S0 mode wave is effectively suppressed, which proves the effectiveness of the method of the present application.

[0079] The specific embodiments of the present application are described above. It should be understood that the present application is not limited to the above specific embodiments, and various modifications or changes can be made by those skilled in the art within the scope of the claims, which does not affect the essential content of the present application.​

Claims

1. A method of separating multi-modal signals in ultrasonic bending wave logging of a cased well, characterized by, Comprising the following steps: Step1: input the original mixed casing well ultrasonic bending wave logging waveform data, casing well structure and elastic parameters of casing, fluid and cement; Step2: simulate A0 mode wave of instrument centering according to the current casing well structure and elastic parameters; Step3: use the variational mode decomposition algorithm to decompose the original mixed signal into a set of multiple scale intrinsic mode functions; Step4: based on the correlation coefficient of the reconstructed A0 mode wave and the simulated standard A0 mode wave, construct the inversion objective function L: ， wherein, is the optimal number of decomposition levels, v is the multiscale reconstruction vector, A is the instrument strictly centered A0 signal, A' is the A0 signal reconstructed by the variational mode decomposition, Cov is the covariance calculation, and Var is the variance calculation; Step5: use simulated annealing method to optimize the optimal decomposition layer number and reconstruction vector in the inversion objective function; Step 6: the optimal decomposition layer number obtained based on inversion of the objective function and the reconstructed A0 mode wave A' is calculated based on the reconstruction vector v: ， where n v is the total number of elements of the reconstruction vector v, i is the index of an element of the reconstruction vector v, k is the scale index, u k is the intrinsic mode function of the kth scale.

2. The method of separation of the casing ultrasonic bending wave well logging multimodal signal according to claim 1, characterized in that, In Step1, the maximum decomposition layer number of variational mode decomposition is set.

3. The method of separating the multi-modal signals of a cased hole ultrasonic bending wave log according to claim 2, wherein, The maximum decomposition layer number is set to 5-8 layers.

4. The method of separating the multi-modal signals of a cased hole ultrasonic bending wave log according to claim 1, wherein, In Step2, according to the casing well diameter, casing thickness and elastic information of casing, fluid and cement in the well, high-order staggered grid finite difference method is used to simulate A0 mode wave under the condition of strict instrument centering, as the standard A0 waveform.

5. The method of separating the multi-modal signals of a cased hole ultrasonic bending wave log according to claim 1, wherein, In Step3, the original signal is decomposed into a set of multiple scale intrinsic functions with compact support by variational mode decomposition: ， where u k and w k are the intrinsic mode function and its corresponding central frequency of the kth scale, respectively, is the Dirac function, j is the imaginary unit, e is the natural logarithm, π is the circle constant, M is the original mixed signal, t is the time sampling point, and n is the maximum decomposition level.

6. The method of separating the casing well ultrasonic bending wave logging multimodal signals according to claim 1, characterized in that, In Step 5, the initial temperature T0 and the maximum number of temperature drops n need to be preset for the algorithm max A linear cooling method is adopted: ， wherein, and Tnand Tn+1are the nthand (n+1)thtemperatures, is the temperature decrease rate; for a certain temperature The new solution acceptance probability of the inversion objective function for the current temperature N cycles is calculated in an isothermal equilibrium manner based on the Metropolis criterion. ， where L is the inversion objective function, and is the optimal decomposition level and reconstruction vector at the n-th iteration, and is the optimal decomposition level and reconstruction vector at the n+1-th iteration, and the new solution is generated randomly in the neighborhood of the previous solution.

7. The method of separating the multi-modal signal of a cased hole ultrasonic bending wave log according to claim 6, characterized in that, The initial temperature is set to 100-200℃, and the temperature drop rate is set to 0.95-0.98.

Images