A cement velocity inversion method based on array ultrasonic lamb waves

Abstract

The application discloses a cement velocity inversion method based on array ultrasonic Lamb waves, relates to the technical field of oil and gas well exploration and development, and solves the technical problems of low calculation efficiency and limited applicable velocity range of the cement velocity inversion algorithm. The method comprises the following steps: performing space-time domain preprocessing on an array signal to obtain a target guided wave mode, and extracting a complex dispersion curve of the target guided wave mode; a three-layer planar medium model is constructed, and a dispersion matrix is derived and established; acoustic parameters of in-well fluid and a casing are set, initial value sequences of multiple groups of to-be-inverted parameters are generated, and the value of the dispersion matrix determinant is calculated; a loss function is constructed, the to-be-inverted parameters are iteratively updated and cyclically executed until the iteration number reaches a preset iteration threshold; and the initial value sequences of all groups are traversed, and the iteration parameters corresponding to the round in which the loss function reaches a minimum value are extracted. The application widens the applicable working conditions and universality, is convenient for automatic and batch processing, and improves the calculation efficiency.

Description

Technical Field

[0001] This invention relates to the field of oil and gas well exploration and development technology, and in particular to a cement velocity inversion method based on array ultrasonic Lamb waves. Background Technology

[0002] In the exploration and development of oil and gas wells, cementing quality assessment is a crucial step in ensuring wellbore integrity and safe production. In recent years, ultrasonic Lamb wave-based logging technology has been widely used for evaluating the properties of cement in cased wells due to its high radial and axial resolution. In oil and gas cementing quality assessment, obtaining the P-wave and S-wave velocities of the cement sheath outside the casing is a core indicator for evaluating the early strength, sealing integrity, and mechanical properties of the cement.

[0003] Currently, the inversion of cement P-wave velocities mainly relies on two methods. The first is the ray theory method, which determines the velocity by extracting the propagation path and travel time of the sound wave. However, this requires the precise thickness of the cement to be known, which is often unknown under complex well conditions. The second is the mode splitting (frequency trap) method, which utilizes the mode splitting phenomenon that occurs in the A0 mode (Lamb wave fundamental antisymmetric mode) within a specific cement P-wave velocity range. By finding the frequency trap point in the spectrum of the first wave of A0 and locating the corresponding velocity on the dispersion curve, this method is heavily dependent on specific frequency trap physical phenomena and is only applicable to an extremely narrow cement velocity range, which cannot meet the inversion requirements of a wide velocity range.

[0004] Meanwhile, existing residual optimization inversion methods for dispersion curves are extremely time-consuming and cannot be automated. To achieve parameter inversion over a wide velocity range, the conventional approach is to calculate the theoretical dispersion curve using forward modeling, then calculate the residual between it and the measured dispersion curve, and iteratively optimize the parameters to be inverted until the residual is minimized. However, the forward modeling search process for complex dispersion curves is extremely complex. The algorithm needs to perform two-dimensional root-finding for each frequency point in the complex plane, which contains drastically changing exponential terms. The enormous computational overhead severely restricts inversion efficiency and cannot meet the rapid processing needs of massive amounts of data from well logging sites. Secondly, the poles of the multi-level dispersion equations in the complex plane are highly dependent on the initial search range and prior guesses, resulting in poor numerical stability and making the entire inversion evaluation process fragile. There is an urgent need for an efficient and automated inversion method that is free from the limitations of cement thickness and specific velocity ranges, and eliminates the need for a time-consuming two-dimensional forward modeling root-finding process in the complex plane, in order to achieve a stable solution for cement sound velocity.

[0005] In the process of realizing this invention, the inventors discovered at least the following problems in the prior art: The cement velocity inversion algorithm based on array ultrasonic Lamb waves has low computational efficiency and a limited applicable velocity range, making it difficult to provide automated and high-precision evaluation of cementing quality. Summary of the Invention

[0006] The purpose of this invention is to provide a cement velocity inversion method based on arrayed ultrasonic Lamb waves, to solve the technical problems of low computational efficiency, limited applicable velocity range, and difficulty in providing automated and high-precision evaluation of cementing quality in existing arrayed ultrasonic Lamb wave cement velocity inversion algorithms. The various technical effects of the preferred solutions among the many technical solutions provided by this invention are detailed below.

[0007] To achieve the above objectives, the present invention provides the following technical solution: This invention provides a cement velocity inversion method based on arrayed ultrasonic Lamb waves, comprising the following steps: S100: acquiring array signals collected by an ultrasonic Lamb wave logging instrument with array receiving capability, performing spatiotemporal preprocessing on the array signals to obtain the target guided wave mode, and extracting the complex dispersion curve of the target guided wave mode based on a dispersion curve extraction algorithm; S200: constructing a three-layer planar medium model of well fluid, casing, and cement based on displacement potential function and interlayer boundary conditions, and deriving and establishing the dispersion matrix; S300: setting the acoustic parameters of well fluid and casing, generating multiple sets of initial value sequences of parameters to be inverted within a preset physical constraint range of cement, and for each set of initial value sequences, extracting the complex dispersion curve... The data points and current medium parameters are substituted into the dispersion matrix. The dispersion matrix is ​​normalized row by row, and the determinant of the dispersion matrix is ​​calculated for all data points of the complex dispersion curves. S400: Construct a loss function based on the determinant of the dispersion matrix, Poisson's ratio, and acoustic impedance. Optimize the gradient of the objective function with respect to the longitudinal wave velocity, transverse wave velocity, and density of cement through automatic differentiation. Iterate and update the parameters to be inverted using an optimization algorithm and execute the loop until the number of iterations for a single initial value sequence reaches a preset iteration threshold. S500: Traverse all initial value sequences. In the historical iteration records, extract the iteration parameters corresponding to the round in which the loss function reaches its minimum value and output them as cement velocity inversion parameters.

[0008] Preferably, in step S100, the operation of performing spatiotemporal preprocessing on the array signal is as follows: along the trajectory of the target waveguide mode, the array signal is truncated using a Tukey window to suppress interference wave components other than the target waveguide mode.

[0009] Preferably, in step S100, the process of extracting the complex dispersion curve of the target guided wave mode is as follows: the array signal is converted to the frequency space domain, a prediction matrix characterizing the spatial propagation characteristics of the wave field is constructed, and the spatial poles are solved using a subspace algorithm. Finally, the phase velocity and spatial attenuation coefficient of the target guided wave mode at this frequency are extracted as the complex dispersion curve.

[0010] Preferably, in step S200, the process of deriving and establishing the dispersion matrix is ​​as follows: Without considering the sound source, the fluid in the well only contains reflected longitudinal waves. The longitudinal wave displacement potential function of the fluid in the well is... Represented as: in, For liquid along Directional wave number, For along Directional wave number, , For the speed of sound in liquid, Indicates time, Represents the reflection coefficient. Indicates along The coordinates of the axis, Indicates along The coordinates of the axis, Angular frequency, The imaginary unit; The casing is a solid in which both longitudinal and transverse waves exist. Both longitudinal and transverse waves include incident and reflected waves. The displacement potential of the casing includes the longitudinal wave displacement potential. transverse wave displacement potential Longitudinal wave displacement potential transverse wave displacement potential They are represented as follows: , These represent the longitudinal and transverse waves along the casing, respectively. The wave number in the direction, and satisfying , , , These represent the longitudinal wave velocity and the transverse wave velocity of the casing, respectively. , , , These represent the transmission coefficient and reflection coefficient of the longitudinal wave, and the transmission coefficient and reflection coefficient of the transverse wave, respectively. Cement is a solid that simultaneously contains longitudinal and transverse waves. The longitudinal wave displacement potential function in cement is... and transverse wave displacement potential function They are represented as follows: , , , These represent the longitudinal and transverse waves along the cement, respectively. The wave number in the axial direction, and satisfying , , , These represent the longitudinal wave velocity and transverse wave velocity of cement, respectively. , These represent the longitudinal wave transmission coefficient and the transverse wave transmission coefficient of cement, respectively. Based on the reflection coefficient of the fluid in the well The transmission coefficient of longitudinal waves in the casing and reflection coefficient Transmission coefficient of shear waves in the casing and reflection coefficient and the longitudinal wave transmission coefficient of cement and transverse wave transmission coefficient Seven unknown variables are used to construct a 7th order dispersion matrix.

[0011] Preferably, in step S300, a set of acoustic parameters of water are used as acoustic parameters of the fluid in the well, and acoustic parameters of steel are used as acoustic parameters of the casing. The acoustic parameters of the steel are obtained through the well construction record, and the thickness of the steel is extracted by ultrasonic pulse wave signal.

[0012] Preferably, in step S300, the preset physical constraint range of the cement includes the cement longitudinal wave velocity range. Cement shear wave velocity range Cement density range and the Poisson's ratio range of cement And based on the multi-starting point strategy, generate multiple sets of initial value sequences for the parameters to be inverted; Furthermore, a parameter mapping mechanism based on the Sigmoid function is constructed to map bounded parameters within the physical constraint interval to unconstrained latent variables. Arbitrary parameters to be inverted and their corresponding upper and lower limits With unconstrained potential In variables The mapping relationship is as follows: Initial values ​​of multiple latent variables are generated uniformly. The system automatically maps and obtains multiple sets of initial value sequences of cement longitudinal wave velocity, cement transverse wave velocity, and cement density with uniform spatial distribution.

[0013] Preferably, in step S400, the loss function The expression is: , in, This represents the mean of the determinant of the dispersion matrix. This represents the Poisson's ratio penalty term for cement. This represents the acoustic impedance penalty term; The mean expression for the determinant of the dispersion matrix is: , in , , These represent the angular frequency, real wavenumber, and imaginary wavenumber of the data points in the complex dispersion curve, respectively. Indicates the number of valid scatter points. , Indicates the first The dispersion matrix at each dispersion point Indicates the calculation modulus. This indicates the calculation of the determinant of a matrix; The current Poisson's ratio of the cement is calculated based on the input longitudinal wave velocity and transverse wave velocity. When the cement Poisson's ratio exceeds the set range... Calculate the square of the excess amount and assign a penalty weight as the cement Poisson's ratio penalty term. The expression for the cement Poisson's ratio penalty term is: , Represents Poisson's ratio. Indicates the penalty weighting factor; The expression for the acoustic impedance penalty term is: , , in, The weights for the impedance loss term, for Activation function This indicates the density of cement. This indicates the longitudinal wave velocity of the cement. This represents the acoustic impedance calculated through algorithmic inversion. This represents the reference acoustic impedance obtained by ultrasonic pulse-echo measurement. This indicates that the error between the acoustic impedance obtained by the inversion algorithm and the reference acoustic impedance is allowed. No penalty is imposed within this error range, and a penalty is imposed if the error exceeds this range.

[0014] Preferably, in step S400, after calculating the gradients of the longitudinal wave velocity, transverse wave velocity, and density of cement, an adaptive moment estimation optimization algorithm is used to dynamically adjust the update step size of each latent variable based on the calculated gradients. The updated latent variables are then remapped to the longitudinal wave velocity, transverse wave velocity, and density with physical constraints and incorporated into the next round of iteration calculation until the preset maximum number of iterations is reached.

[0015] Preferably, in step S400, the step of optimizing the gradient of the objective function with respect to the longitudinal wave velocity, transverse wave velocity, and density of cement by automatic differentiation specifically includes: dynamically constructing a computational graph from the parameters to be inverted to the scalar value of the loss function in the program memory, and calculating the gradient of each latent variable of cement in the loss function according to the chain rule using the backpropagation algorithm.

[0016] Preferably, in step S500, after traversing the initial value sequence of all groups, the parameter combination with the minimum loss function is found from all historical iteration loss function values, the latent variable corresponding to the sequence is extracted, the latent variable is solved in reverse, and finally the cement velocity under the real physical dimensions is obtained.

[0017] Implementing one of the above-described technical solutions of the present invention has the following advantages or beneficial effects: This application utilizes complex dispersion data within a certain bandwidth range for direct velocity optimization. Unlike ray theory methods, it does not require pre-assuming cement thickness, nor is it limited by the specific narrow velocity range required by mode splitting (frequency trapping) methods, thus broadening the applicability and universality of the inversion algorithm. By directly substituting the measured complex dispersion data points into the dispersion matrix and constructing the mismatch objective function using the mean of the logarithmic determinant of this matrix, this strategy eliminates the extremely time-consuming process of two-dimensional root-finding for dispersion equations containing drastically changing exponential terms in the complex plane. The root-finding of complex multilayer medium equations is transformed into dispersion matrix determinant calculation and automatic differential gradient optimization, providing an algorithmic foundation for the automated and batch processing of ultrasonic logging data. By introducing multiple sets of preset physical constraint initial values ​​and parallel traversal, the algorithm can automatically escape local minima, exhibiting strong robustness and stably outputting P-wave and S-wave velocity inversion results that conform to physical common sense. Attached Figure Description

[0018] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort. In the drawings: Figure 1 This is a flowchart of a cement velocity inversion method based on array ultrasonic Lamb waves according to an embodiment of the present invention. Figure 1 ; Figure 2 This is a flowchart of a cement velocity inversion method based on array ultrasonic Lamb waves according to an embodiment of the present invention. Figure 2 ; Figure 3This is a waveform diagram of the ultrasonic Lamb wave array received in a forward simulation under a four-layer planar medium model in a cement velocity inversion method based on array ultrasonic Lamb waves according to an embodiment of the present invention. Figure 4 This is a schematic diagram of Tukey window processing of the array signal in a cement velocity inversion method based on array ultrasonic Lamb waves according to an embodiment of the present invention. Figure 5 This is a schematic diagram of the target guided wave signal after windowing the array signal in a cement velocity inversion method based on array ultrasonic Lamb waves according to an embodiment of the present invention. Figure 6 This is a frequency-real part curve of the complex dispersion curve of the Lamb wave extracted using the matrix pencil algorithm in a cement velocity inversion method based on array ultrasonic Lamb waves according to an embodiment of the present invention. Figure 7 This is a frequency-imaginary part curve of the complex dispersion curve of the Lamb wave extracted using the matrix pencil algorithm in a cement velocity inversion method based on array ultrasonic Lamb waves according to an embodiment of the present invention. Figure 8 This is a schematic diagram of the physical model of the water-steel-cement three-layer planar dielectric waveguide propagation and its wave field reflection and refraction path constructed in a cement velocity inversion method based on array ultrasonic Lamb waves according to an embodiment of the present invention. Figure 9 This is a comparison of the longitudinal wave velocity inversion results and theoretical true values ​​of a cement velocity inversion method based on array ultrasonic Lamb waves under 40 different cement models according to an embodiment of the present invention. Figure 10 This is a comparison of the transverse wave velocity inversion results and theoretical true values ​​of a cement velocity inversion method based on array ultrasonic Lamb waves under 40 different cement models, according to an embodiment of the present invention. Detailed Implementation

[0019] To make the objectives, technical solutions, and advantages of the present invention clearer, various exemplary embodiments described below will be referenced to the accompanying drawings, which form part of the exemplary embodiments, illustrating various exemplary embodiments that may be used to implement the present invention. Unless otherwise indicated, the same numbers in different drawings represent the same or similar elements. The embodiments described in the following exemplary embodiments do not represent all embodiments consistent with this disclosure. It should be understood that they are merely examples of processes, methods, and apparatuses consistent with some aspects of the present invention disclosed as detailed in the appended claims, and other embodiments may be used, or structural and functional modifications may be made to the embodiments listed herein without departing from the scope and spirit of the present invention.

[0020] In the description of this invention, it should be understood that the terms "center," "longitudinal," "lateral," etc., indicate the orientation or positional relationship based on the accompanying drawings, and are only for the convenience of describing the invention and simplifying the description, and do not indicate or imply that the referred element must have a specific orientation, or be constructed and operated in a specific orientation. The terms "first," "second," etc., are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of indicated technical features. The term "multiple" means two or more. The terms "connected" and "linked" should be interpreted broadly, for example, they can be fixed connections, detachable connections, integral connections, mechanical connections, electrical connections, communication connections, direct connections, indirect connections through an intermediate medium, and can be the internal connection of two elements or the interaction relationship between two elements. The term "and / or" includes any and all combinations of one or more of the related listed items. Those skilled in the art can understand the specific meaning of the above terms in this invention according to the specific circumstances.

[0021] To illustrate the technical solution described in this invention, specific embodiments are described below, showing only the parts related to the embodiments of this invention.

[0022] Example: like Figure 1As shown, this invention provides a cement velocity inversion method based on array ultrasonic Lamb waves, including the following steps: S100: Acquire array signals collected by an ultrasonic Lamb wave logging instrument with array receiving capability. The ultrasonic Lamb wave logging instrument is a casing well logging device, mainly used to accurately evaluate cementing quality. The array signal is a set of signals collected by arranging multiple transducers (array elements) according to specific rules. The array signal is subjected to spatiotemporal preprocessing (preprocessing is used to suppress interference wave components other than the target guided wave mode) to obtain the target guided wave mode. The complex dispersion curve of the target guided wave mode is extracted based on the dispersion curve extraction algorithm. The dispersion curve extraction algorithm accurately separates the wave velocity (phase velocity or group velocity) and frequency relationship of different modes (such as A0, S0, etc.) from the complex measured signal. The complex dispersion curve describes the wave propagation speed (dispersion characteristics) and energy attenuation (attenuation characteristics) at the same time. It includes phase velocity (which determines the speed of wave propagation and corresponds to the real part of the complex wave number) and spatial attenuation information (which determines the energy loss of the wave and corresponds to the imaginary part of the complex wave number). S200: Based on the displacement potential function and interlayer boundary conditions, a three-layer planar medium model of well fluid (liquid water in this embodiment), casing (solid steel in this embodiment) and cement is constructed, and the dispersion matrix is ​​derived. Since the actual single-layer casing well can be regarded as a four-layer medium model of water-steel-cement-formation, but the first wave of the Lamb wave is independent of the formation, the cement parameters can be inverted using only the three-layer medium model of water-steel-cement, thereby reducing the number of unknown parameters to be inverted and the dimension of the dispersion matrix, and improving the inversion efficiency of parameters. S300: Set the acoustic parameters of the well fluid and casing, and generate multiple sets of initial value sequences for the parameters to be inverted within a preset physical constraint range. For each set of initial value sequences, substitute the data points of the complex dispersion curves and the current medium parameters (i.e., the acoustic parameters of the well fluid and casing) into the dispersion matrix. Normalize the dispersion matrix row by row (e.g., ensure that the data in each row of the dispersion matrix has a sum of 1 or a modulus of 1, making different input parameters comparable and eliminating differences in total values), and calculate the determinant of the dispersion matrix for all data points of the complex dispersion curves. S400: Construct a loss function based on the determinant of the dispersion matrix, Poisson's ratio, and acoustic impedance. Optimize the gradient of the objective function with respect to the P-wave velocity, S-wave velocity, and density of cement through automatic differentiation calculation. Iterate and update the parameters to be inverted using an optimization algorithm and execute the process repeatedly until the number of iterations for a single set of initial value sequences reaches a preset iteration threshold. S500: Traverse the initial value sequence of all groups, extract the iteration parameters corresponding to the round in which the loss function reaches its minimum value from the historical iteration record, and output them as cement speed inversion parameters.The cement velocity inversion method in this embodiment directly optimizes velocity using complex dispersion data within a certain bandwidth range. Unlike ray theory methods, it does not require pre-assuming cement thickness, nor is it limited by the specific narrow velocity range required by mode splitting (frequency trapping) methods, thus broadening the applicability and universality of the inversion algorithm. The measured complex dispersion data points are directly substituted into the dispersion matrix, and the mismatch objective function is constructed using the mean of the logarithmic values ​​of the determinant of this matrix. This strategy eliminates the extremely time-consuming process of two-dimensional root finding for dispersion equations containing drastically changing exponential terms in the complex plane. The root finding of complex multilayer medium equations is transformed into dispersion matrix determinant calculation and automatic differential gradient optimization, providing an algorithmic foundation for the automated and batch processing of ultrasonic logging data. The introduction of multiple sets of preset physical constraint initial values ​​and parallel traversal enable the algorithm to automatically escape local minima, exhibiting strong robustness and stably outputting P-wave and S-wave velocity inversion results that conform to physical common sense.

[0023] As an optional implementation, in step S100, the spatiotemporal preprocessing of the array signal involves truncating the array signal along the trajectory of the target guided wave mode using a Tukey window to suppress interference wave components other than the target guided wave mode. Due to the complexity of ultrasonic waveforms, to suppress interference from volume waves and other modes, this embodiment constructs a two-dimensional spatiotemporal picking window within the spatiotemporal two-dimensional data matrix. Specifically, the array signal is truncated using a Tukey window along the trajectory of the target guided wave mode (such as the A0 mode, the lowest-order antisymmetric mode, a bending wave primarily vibrating in the thickness direction, which is extremely sensitive to changes in the surface and near-surface). Since the Tukey window is a time window that is flat in the middle and smoothly descends at both ends, it ensures that the target guided wave has no amplitude attenuation, thus reducing its impact on subsequent dispersion curve extraction. In this embodiment, the model settings are: sleeve thickness of 12mm, longitudinal wave velocity of 5860m / s, transverse wave velocity of 3130m / s, and density of 7850kg / m^3. The sound velocity of the fluid in the well is 1500 m / s, and the density is 1000 kg / m³. Forty different cement parameters were set: P-wave velocity ranged from 1960 to 4000 m / s, with 40 samples taken at uniform intervals; S-wave velocity ranged from 1120 to 2100 m / s, with 40 samples taken at uniform intervals; density ranged from 1200 to 2400 kg / m³, with 40 samples taken at uniform intervals; and cement thickness was 3 cm. The formation density was 2320 m / s, the P-wave velocity was 4500 m / s, and the S-wave velocity was 2455 m / s. The sound source used a Ricker wavelet with a dominant frequency of 220 kHz and a source width of 4 cm. Ten waveforms within a source distance of 15 cm to 24 cm were simulated using analytical methods, with a receiver spacing of 1 cm. The calculated array waveforms are shown below. Figure 3As shown. Due to the presence of the formation, the received signal contains many reflected signals from the cement-formation interface, which will affect the extraction of subsequent array signal dispersion. Therefore, it is necessary to suppress and eliminate these subsequent signals. For example... Figure 4 As shown, a Tukey window with a reasonable length and position can meet this requirement. The basic steps are: (1) Find the time of the first wave envelope peak of the Lamb wave. ,by The center of the time window is the frequency of the sound source, and the length of the time window can be determined based on the dominant frequency of the sound source. Therefore, the time window length used in this case is determined to be... The set time window is shown in the figure. After preprocessing, the following can be obtained: Figure 5 The array signal is shown.

[0024] As an optional implementation, in step S100, the process of extracting the complex dispersion curve of the target guided wave mode is as follows: the array signal is converted to the frequency space domain, which is convenient for obtaining the frequency components of the wave and their corresponding spatial rate of change (wave number), a prediction matrix characterizing the spatial transmission characteristics of the wave field is constructed, and the spatial poles are solved using a subspace algorithm (such as the matrix pencil algorithm), and finally the phase velocity and spatial attenuation coefficient of the target guided wave mode at this frequency are extracted as the complex dispersion curve.

[0025] In this embodiment, the target guided wave time-domain signal extracted after Tukey window truncation is subjected to a Fast Fourier Transform along the time axis to obtain the array spatial response sequence at a specific angular frequency. ,in For the first The axial coordinates of each receiver, , Let be the total number of array receivers. At a single frequency, the spatial sequence of the wavefield propagating along the equally spaced array can be decomposed into: The superposition of individual modes and noise: in, For the first The complex amplitude of the target guided wave mode, To characterize the complex poles of the wavefield spatial propagation operator, it is defined as follows: ,in , Wave number The real part, Wave number The imaginary part, For receiver spacing, It is noise.

[0026] In this embodiment, the specific process of the matrix pencil algorithm is as follows.

[0027] At a single frequency Below, using spatial sequences Construct translation-invariant Hankel matrix ,in The matrix is: in These are parameters for the Hankel matrix, used to control the number of columns in the Hankel matrix. They are typically set to a certain value. .right Perform singular value decomposition: Set a threshold for singular values, which determines the number of modes in the signal. It equals the number of singular values ​​greater than the threshold. (Using a matrix...) Construct two subarrays and ,in It is a matrix Take before Column and remove the last row of the result. It is a matrix Take before The result of the column union with the first row removed, satisfying the following conditions: Solving for the poles using the above formula , The basic form is , Let the value be an undetermined positive integer, representing the distance between receivers when the requirements are met. , for The range, for The phase of the guided wave at that frequency. and attenuation Both can be obtained from the extreme point get: However, when the spacing between the receivers is too large, The solution requires careful attention to the phase selection to ensure that the phase velocity of the guided wave is within the theoretical range.

[0028] The calculated dispersion curve is as follows Figures 6-7As shown, the extracted dispersion curves reveal a significant difference between the S0 mode and the theoretical dispersion curve. This is primarily due to the larger group velocity distribution range of the S0 mode, resulting in a longer duration. Time window truncation distorts the dispersion curve. However, the dispersion curve of the A0 mode extracted by the algorithm closely matches the theoretical dispersion curve. Therefore, selecting the A0 mode dispersion curve for retrieving cement parameters is more appropriate.

[0029] As an optional implementation, in step S200, the process of deriving and establishing the dispersion matrix is ​​as follows: Without considering the sound source, the fluid in the well only contains reflected longitudinal waves. The longitudinal wave displacement potential function of the fluid in the well is... Represented as: in, For liquid along Directional wave number, For along Directional wave number, , For the speed of sound in liquid, Indicates time, Represents the reflection coefficient. Indicates along The coordinates of the axis, Indicates along The coordinates of the axis, Angular frequency, The imaginary unit; The casing is a solid in which both longitudinal and transverse waves exist. Both longitudinal and transverse waves include incident and reflected waves. The displacement potential of the casing includes the longitudinal wave displacement potential. transverse wave displacement potential Longitudinal wave displacement potential transverse wave displacement potential They are represented as follows: , These represent the longitudinal and transverse waves along the casing, respectively. The wave number in the direction, and satisfying , , , These represent the longitudinal wave velocity and the transverse wave velocity of the casing, respectively. , , , These represent the transmission coefficient and reflection coefficient of the longitudinal wave, and the transmission coefficient and reflection coefficient of the transverse wave, respectively. Cement is a solid that simultaneously contains longitudinal and transverse waves. The longitudinal wave displacement potential function in cement is... and transverse wave displacement potential function They are represented as follows: , , , These represent the longitudinal and transverse waves along the cement, respectively. The wave number in the axial direction, and satisfying , , , These represent the longitudinal wave velocity and transverse wave velocity of cement, respectively. , These represent the longitudinal wave transmission coefficient and the transverse wave transmission coefficient of cement, respectively. Based on the reflection coefficient of the fluid in the well The transmission coefficient of longitudinal waves in the casing and reflection coefficient Transmission coefficient of shear waves in the casing and reflection coefficient and the longitudinal wave transmission coefficient of cement and transverse wave transmission coefficient Seven unknown variables, namely For unknown coefficients, based on 7 boundary conditions, at each angular frequency and A simultaneous equation can be formed to obtain a statement about the variable. The homogeneous linear equation system, i.e. ,in , It is a 7×7 dispersion matrix.

[0030] In this embodiment, the displacement and stress are calculated as follows: The two-dimensional (xz plane) problems involved include: in , They represent along shaft and Displacement in the axial direction, Γ represents the longitudinal wave displacement potential, and Γ represents the transverse wave displacement potential. The method for calculating solid stress is as follows: in , Indicates normal stress and tangential stress. , The Latin American coefficient is calculated as follows: , ,in For the density of the medium, For the longitudinal wave velocity of the medium, The transverse wave velocity of the medium; For fluid media, the stress calculation method is as follows: In this embodiment, the displacement potential functions of each medium layer are related through boundary conditions. Interface I and interface II are as follows: Figure 8 As shown, the first layer is water, the second layer is steel, and the third layer is cement. ) indicates the first Longitudinal wave velocity of the layered medium ) indicates the first Transverse wave velocity of the layered medium ) indicates the first The density of the medium layer, and thus the specific boundary conditions are as follows.

[0031] At the water-steel boundary, the displacement potential function satisfies the fluid-solid boundary condition: Normal displacement continuity: in This represents the normal displacement of the first layer of medium (water) at interface I. This represents the normal displacement of the second layer of medium (steel) at interface I; Normal stress continuity: in This represents the normal stress at interface I of the first layer of medium (water). This represents the normal stress at interface I of the second layer of medium (steel); The tangential stress is 0. in This represents the tangential stress at interface I of the second layer of medium (steel).

[0032] There are two types of boundary conditions. There are two interfaces in the third-layer planar medium of water-casing-cement: water-casing and casing-cement. The former is a fluid-solid interface (interface I), and the latter is a solid-solid interface (interface II).

[0033] At the steel-cement boundary, the displacement potential function satisfies the solid-solid boundary condition: Normal displacement continuity: in This indicates the normal displacement of the second layer of medium (steel) at interface II. This indicates the normal displacement of the third layer of medium (cement) at interface II; Tangential displacement is continuous: in This indicates the tangential displacement of the second layer of medium (steel) at interface II. This indicates the tangential displacement of the third layer of medium (cement) at interface II; Normal stress continuity: in This represents the normal stress at interface II of the second layer of medium (steel). This represents the normal stress at interface II of the third layer of medium (steel). Tangential stress continuity: in This represents the tangential stress at interface II of the second layer of medium (steel). This represents the tangential stress at interface II of the third layer of medium (steel).

[0034] Based on the boundary conditions, displacement and stress calculation methods, and combined with the displacement potential function, at each angular frequency... and A simultaneous equation can be formed to obtain a statement about the variable. The homogeneous linear equation system: in The dispersion matrix is ​​7×7, and therefore each term in the matrix is ​​expressed as: , , , , , , ; , , , , ， ， ； ， ， ， ， ， ， ； ， ， ， ， ， ， ； ， ， ， ， ， ， ； ， ， ， ， ， ， ； ， ， ， ， ， ， ； in , , , For the first Lami coefficient of the layered medium These correspond to the fluid inside the well, the casing, and the cement, respectively. The thickness of the sleeve is denoted by the Lame coefficient, which is one of the two important material parameters used to describe the elastic properties of homogeneous isotropic materials (mediums).

[0035] As an optional implementation, in step S300, a set of acoustic parameters of water are used as the acoustic parameters of the fluid in the well. In actual casing wells, the properties of the fluid are usually close to those of water, since water is usually used as the fluid in the well. In actual engineering, if the properties of the fluid in the well are known, the known fluid parameters are used instead. The acoustic parameters of steel are used as the acoustic parameters of the casing. The acoustic parameters of steel are obtained through the well construction records, and the thickness of the steel is extracted through ultrasonic pulse wave signals.

[0036] As an optional implementation, in step S300, the preset physical constraint range of the cement includes the cement longitudinal wave velocity range. Cement shear wave velocity range Cement density range and the Poisson's ratio range of cement The algorithm generates multiple sets of initial value sequences for the parameters to be inverted based on a multi-starting-point strategy. This strategy increases the probability of finding the global optimum by exploring from multiple different initial points and avoids the risk of gradient-based optimization algorithms getting trapped in local minima during the inversion process. In this embodiment, a three-layer planar medium of water-steel-cement (structure as shown) is derived using the displacement potential function and boundary conditions. Figure 8 The dispersion matrix is ​​shown in the figure. The parameters of water and casing are set to be consistent with the theoretical values, and the iteration range of cement parameters is set. The iteration range of longitudinal wave velocity of cement is 1000-4500m / s, the iteration range of transverse wave velocity is 1000-2500m / s, and the iteration range of density is 1000-2500kg / m^3. In order to ensure that the calculated parameters are consistent with the characteristics of cement, the iteration range of Poisson's ratio of cement is additionally set to 0.1-0.4. If the Poisson's ratio exceeds the set iteration range, an additional penalty loss will be added.

[0037] To avoid the inversion results getting trapped in local optima, this embodiment first uses 1000 randomly distributed sequences as the preset physical constraint range for cement. Simultaneously, additional physical constraints are applied to the values ​​of the three parameters: the Poisson's ratio of cement is limited to 0.1-0.4, and the cement density must not exceed the longitudinal wave velocity. Then, calculations are performed on these 1000 sets of initial parameters. ,Pick The 100 sets of parameters with the smallest values ​​are used as the final initial value points and then incorporated into the subsequent iteration steps.

[0038] To prevent parameters from exceeding the preset physical constraint range during the optimization process, a parameter mapping mechanism based on the Sigmoid function is constructed to map bounded parameters within the physical constraint range to unconstrained latent variables. Arbitrary parameters to be inverted and their corresponding upper and lower limits With unconstrained latent variables The mapping relationship is as follows: Initial values ​​of multiple latent variables are generated uniformly. The system automatically maps and obtains multiple sets of initial value sequences for cement longitudinal wave velocity, cement transverse wave velocity, and cement density, all with uniform spatial distribution. Subsequent optimization parameters are no longer the physical parameters themselves, but rather the parameters they map to. The Sigmoid mapping is globally smooth and differentiable. As physical parameters approach set upper and lower limits, the derivative of the Sigmoid function naturally decreases, strictly limiting the parameters within a given range while ensuring gradient continuity. This is achieved by using the mapping parameters... The goal is to map the cement longitudinal wave velocity, cement transverse wave velocity, and cement density so that the optimizer no longer directly optimizes these three dimensionless latent variables that fluctuate around 0 and have the same magnitude. This allows the optimizer to use a single learning rate to efficiently drive the three physical quantities to converge simultaneously.

[0039] As an optional implementation, in step S400, the loss function The expression is: , in, This represents the mean of the determinant of the dispersion matrix. This represents the Poisson's ratio penalty term for cement. This represents the acoustic impedance penalty term; The mean expression for the determinant of the dispersion matrix is: , in , , These represent the angular frequency, real wavenumber, and imaginary wavenumber of the data points in the complex dispersion curve, respectively. Indicates the number of valid scatter points. , Indicates the first The dispersion matrix at each dispersion point Indicates the calculation modulus. This represents the calculation of the matrix determinant; extracting the effective complex dispersion data points (angular frequencies) Real wavenumber Imaginary wavenumber is the attenuation coefficient. The media parameters of the current iteration are substituted into the dispersion matrix. Then, the logarithm of the determinant of the dispersion matrix is ​​calculated to base 10, and the average of the logarithms of the determinants for all dispersion data points is taken as the optimization target for the parameters. The current Poisson's ratio of the cement is calculated based on the input longitudinal wave velocity and transverse wave velocity. When the cement Poisson's ratio exceeds the set range... Calculate the square of the excess amount and assign a penalty weight as the cement Poisson's ratio penalty term. The expression for the cement Poisson's ratio penalty term is: , Represents Poisson's ratio. This represents the penalty weighting factor, the value of which can be adjusted according to the actual situation; The expression for the acoustic impedance (which reflects the resistance encountered by sound waves when propagating in a medium) penalty term is: , , in, The weights for the impedance loss term, for Activation function This represents the acoustic impedance calculated through algorithmic inversion. This represents the reference acoustic impedance obtained by ultrasonic pulse-echo measurement. This indicates that the error between the acoustic impedance obtained by the inversion algorithm and the reference acoustic impedance is allowed; no penalty is imposed within this error range, and a penalty is imposed if the error exceeds this range. In this embodiment, the preferred method is... , It is the unit of acoustic impedance. The value is 10^6. The cement velocity obtained by this algorithm from the array waveform inversion of 40 models, and the inversion error are as follows: Figures 9-10 As shown, inversion was achieved by adjusting the shear wave velocities of 40 different types of cement (e.g., Figure 9 The average error shown is only 3.83%, and effective tracking of the P-wave velocity is also achieved (e.g. Figure 10 The average error shown is approximately 3.68%, indicating extremely high potential for industrial automation applications and engineering value.

[0040] As an optional implementation, in step S400, after calculating the gradients of the P-wave velocity, S-wave velocity, and density of cement, an adaptive moment estimation optimization algorithm (such as the Adam optimization algorithm; in this case, the optimizer learning rate is set to 0.002 and the iteration step size is 500 steps. To improve the inversion efficiency of the algorithm, GPU-based parallel computing can be considered in the actual inversion process) is used to dynamically adjust the update step size of each latent variable based on the calculated gradients. The updated latent variables are then remapped to P-wave velocity, S-wave velocity, and density with physical constraints and substituted into the next round of iteration calculation until the preset maximum number of iterations is reached.

[0041] As an optional implementation, in step S400, the step of automatically differentiating and calculating the gradient of the objective function with respect to the P-wave velocity, S-wave velocity, and density of cement specifically includes: dynamically constructing a loss function from the parameters to be inverted (unconstrained latent variable sequence) in the program memory. The computation graph of the scalar value of the final optimization objective function is used to calculate the loss function using the backpropagation algorithm following the chain rule (for easy automatic and accurate calculation). The gradients of various latent variables related to cement.

[0042] As an optional implementation, in step S500, after traversing the initial value sequence of all groups, the parameter combination with the smallest loss function is found from all historical iteration loss function values, the latent variable corresponding to the sequence is extracted, and the latent variable is substituted into the preset Sigmoid mapping function for reverse calculation, and finally the cement velocity under the real physical dimensions is obtained.

[0043] The embodiment is merely a specific example and does not indicate that this is the only way to implement the present invention.

[0044] The above description is merely a preferred embodiment of the present invention. Those skilled in the art will understand that various changes or equivalent substitutions can be made to these features and embodiments without departing from the spirit and scope of the present invention. Furthermore, under the teachings of the present invention, these features and embodiments can be modified to adapt to specific situations and materials without departing from the spirit and scope of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed herein, and all embodiments falling within the scope of the claims of this application are within the protection scope of the present invention.

Claims

1. A cement velocity inversion method based on array ultrasonic Lamb waves, characterized in that, Includes the following steps: S100: Acquire array signals collected by an ultrasonic Lamb wave logging instrument with array receiving capability, perform spatiotemporal preprocessing on the array signals to obtain the target guided wave mode, and extract the complex dispersion curve of the target guided wave mode based on the dispersion curve extraction algorithm. S200: Based on the displacement potential function and interlayer boundary conditions, a three-layer planar medium model of well fluid, casing and cement is constructed, and the dispersion matrix is ​​derived and established. S300: Set the acoustic parameters of the fluid and casing in the well, generate multiple sets of initial value sequences of parameters to be inverted within the preset physical constraint range of cement, and for each set of initial value sequences, substitute the data points of the complex dispersion curve and the current medium parameters into the dispersion matrix, normalize the dispersion matrix by row, and calculate the value of the determinant of the dispersion matrix under the data points of all complex dispersion curves. S400: Construct a loss function based on the determinant of the dispersion matrix, Poisson's ratio, and acoustic impedance value. Optimize the gradient of the objective function with respect to the longitudinal wave velocity, transverse wave velocity, and density of cement through automatic differentiation calculation. Iteratively update the parameters to be inverted using an optimization algorithm and execute the process repeatedly until the number of iterations for a single set of initial value sequences reaches the preset iteration threshold. S500: Traverse the initial value sequence of all groups, extract the iteration parameters corresponding to the round in which the loss function reaches the minimum value in the historical iteration record, and output them as cement speed inversion parameters. In step S100, the process of extracting the complex dispersion curve of the target waveguide mode is as follows: the array signal is converted to the frequency space domain, a prediction matrix characterizing the spatial propagation characteristics of the wave field is constructed, and the spatial poles are solved using the subspace algorithm. Finally, the phase velocity and spatial attenuation coefficient of the target waveguide mode at this frequency are extracted as the complex dispersion curve. In step S200, the process of deriving and establishing the dispersion matrix is ​​as follows: Without considering the sound source, the fluid in the well only contains reflected longitudinal waves. The longitudinal wave displacement potential function of the fluid in the well is... Represented as: in, For liquid along Directional wave number, For along Directional wave number, , For the speed of sound in liquid, Indicates time, Represents the reflection coefficient. Indicates along The coordinates of the axes, Indicates along The coordinates of the axes, Angular frequency, The imaginary unit; The casing is a solid in which both longitudinal and transverse waves exist. Both longitudinal and transverse waves include incident and reflected waves. The displacement potential of the casing includes the longitudinal wave displacement potential. transverse wave displacement potential Longitudinal wave displacement potential transverse wave displacement potential They are represented as follows: , These represent the longitudinal and transverse waves along the casing, respectively. The wave number in the direction, and satisfying , , , These represent the longitudinal wave velocity and the transverse wave velocity of the casing, respectively. , , , These represent the transmission coefficient and reflection coefficient of the longitudinal wave, and the transmission coefficient and reflection coefficient of the transverse wave, respectively. Cement is a solid that simultaneously contains longitudinal and transverse waves. The longitudinal wave displacement potential function in cement is... and transverse wave displacement potential function They are represented as follows: ， ， , These represent the longitudinal and transverse waves along the cement, respectively. The wave number in the axial direction, and satisfying , , , These represent the longitudinal wave velocity and transverse wave velocity of cement, respectively. , These represent the longitudinal wave transmission coefficient and the transverse wave transmission coefficient of cement, respectively. Based on the reflection coefficient of the fluid in the well The transmission coefficient of longitudinal waves in the casing and reflection coefficient Transmission coefficient of shear waves in the casing and reflection coefficient and the longitudinal wave transmission coefficient of cement and transverse wave transmission coefficient Seven unknown variables are used to construct a 7th order dispersion matrix.

2. The cement velocity inversion method based on array ultrasonic Lamb waves according to claim 1, characterized in that, In step S100, the operation of performing spatiotemporal preprocessing on the array signal is as follows: along the trajectory of the target waveguide mode, the array signal is truncated using a Tukey window to suppress interference wave components other than the target waveguide mode.

3. The cement velocity inversion method based on array ultrasonic Lamb waves according to claim 1, characterized in that, In step S300, a set of acoustic parameters of water are used as the acoustic parameters of the fluid in the well, and the acoustic parameters of steel are used as the acoustic parameters of the casing. The acoustic parameters of steel are obtained through the well construction records, and the thickness of steel is extracted through ultrasonic pulse wave signals.

4. The cement velocity inversion method based on array ultrasonic Lamb waves according to claim 1, characterized in that, In step S300, the preset physical constraint range for cement includes the cement longitudinal wave velocity range. Cement shear wave velocity range Cement density range and the Poisson's ratio range of cement And based on the multi-starting point strategy, generate multiple sets of initial value sequences for the parameters to be inverted; Furthermore, a parameter mapping mechanism based on the Sigmoid function is constructed to map bounded parameters within the physical constraint interval to unconstrained latent variables. , Arbitrary parameters to be inverted and their corresponding upper and lower limits With unconstrained latent variables The mapping relationship is as follows: Initial values ​​of multiple latent variables are generated uniformly. The system automatically maps and obtains multiple sets of initial value sequences of cement longitudinal wave velocity, cement transverse wave velocity, and cement density with uniform spatial distribution.

5. The cement velocity inversion method based on array ultrasonic Lamb waves according to claim 1, characterized in that, In step S400, the loss function The expression is: ， in, This represents the mean of the determinant of the dispersion matrix. This represents the Poisson's ratio penalty term for cement. This represents the acoustic impedance penalty term; The mean expression for the determinant of the dispersion matrix is: ， in , , These represent the angular frequency, real wavenumber, and imaginary wavenumber of the data points in the complex dispersion curve, respectively. Indicates the number of valid scatter points. , Indicates the first The dispersion matrix at each dispersion point Indicates the calculation modulus. This indicates the calculation of the determinant of a matrix; The current Poisson's ratio of the cement is calculated based on the input longitudinal wave velocity and transverse wave velocity. When the cement Poisson's ratio exceeds the set range... Calculate the square of the excess amount and assign a penalty weight as the cement Poisson's ratio penalty term. The expression for the cement Poisson's ratio penalty term is: ， Represents Poisson's ratio. This represents the penalty weighting factor; The expression for the acoustic impedance penalty term is: ， ， in, The weights for the impedance loss term, for Activation function This indicates the density of cement. This indicates the longitudinal wave velocity of the cement. This represents the acoustic impedance calculated through algorithmic inversion. This represents the reference acoustic impedance obtained by ultrasonic pulse-echo measurement. This indicates that the error between the acoustic impedance obtained by the inversion algorithm and the reference acoustic impedance is allowed. No penalty is imposed within this error range, and a penalty is imposed if the error exceeds this range.

6. The cement velocity inversion method based on array ultrasonic Lamb waves according to claim 1, characterized in that, In step S400, after calculating the gradients of the longitudinal wave velocity, transverse wave velocity, and density of cement, an adaptive moment estimation optimization algorithm is used to dynamically adjust the update step size of each latent variable based on the calculated gradients. The updated latent variables are then remapped to the longitudinal wave velocity, transverse wave velocity, and density with physical constraints and substituted into the next round of iteration calculation until the preset maximum number of iterations is reached.

7. The cement velocity inversion method based on array ultrasonic Lamb waves according to claim 1, characterized in that, In step S400, the step of optimizing the gradient of the objective function with respect to the longitudinal wave velocity, transverse wave velocity, and density of cement through automatic differentiation calculation specifically includes: dynamically constructing a computational graph from the parameters to be inverted to the scalar value of the loss function in the program memory, and calculating the gradient of each latent variable of cement in the loss function according to the chain rule through the backpropagation algorithm.

8. The cement velocity inversion method based on array ultrasonic Lamb waves according to claim 1, characterized in that, In step S500, after traversing the initial value sequence of all groups, the parameter combination with the minimum loss function is found from all historical iteration loss function values. The latent variable corresponding to this sequence is extracted, and the latent variable is solved in reverse to finally obtain the cement velocity under the real physical dimensions.

Images