CN115407405A - Method for calculating slowness of array sound waves in laboratory - Google Patents

Abstract

The invention discloses a method for calculating slowness of array sound waves in a laboratory, which comprises the following steps: step 1, carrying out two-dimensional Fourier transform on the acquired monopole array waveform to obtain time-domain stratum longitudinal wave information: step 2, determining the position of the first arrival and the first jump of each row of waveforms by using a conditional fuzzy C-means clustering algorithm: step 3, determining the accurate position of the first arrival by applying a Bayesian information criterion: step 4, carrying out anomaly detection on the first-arrival points by applying a Principal Component Analysis (PCA) method, and fitting the anomaly points to real first-arrival points through norms; and 5, calculating the high-resolution longitudinal wave slowness by applying a slowness solving formula. The method disclosed by the invention has the advantages of less artificial interference, simple parameter setting and strong generalization capability, and is simultaneously suitable for extracting the longitudinal wave slowness of the conventional digital acoustic wave instrument and array acoustic wave instrument; the unsupervised machine learning method using the condition fuzzy clustering and the BIC information criterion has high calculation efficiency.

Description

Method for calculating slowness of array sound waves in laboratory

Technical Field

The invention belongs to the field of slowness extraction of longitudinal wave measurement signals of a miniaturized array acoustic wave instrument in a laboratory, and particularly relates to a slowness calculation method for an array acoustic wave in the laboratory in the field.

Background

The acoustic slowness is the inverse of the propagation speed of an acoustic signal in the formation. In well logging measurement and evaluation, the acoustic slowness curve is used as an indispensable well logging curve in nine conventional well logging lines, has very important application, and can be used for identifying lithology, analyzing borehole stability, calculating formation porosity, estimating formation permeability, evaluating formation anisotropy and the like.

Currently, the slowness of the acoustic wave in the well logging is mainly obtained by detecting the arrival time of the wave and the related value of the slowness, and the calculation method can be divided into 2 types of frequency domain processing method and time domain processing method. The frequency domain processing can be divided into a Prony prediction method, a weighted spectrum correlation method and the like based on different principles, but the noise resistance of the methods is low, and the methods are mostly used as theoretical analysis tools and do not carry out actual calculation. The time domain processing method mainly comprises a threshold value method, a long-time window energy ratio method, a short-time window energy ratio method and a waveform prediction method, wherein the waveform prediction method mainly researches a slowness-time correlation method (STC method) and N times of Fang Genfa. Because the threshold value method and the long and short time window energy ratio method have poor noise immunity, the problem of serious distortion exists in the low signal-to-noise ratio data processing, and the higher precision can be obtained only by manual intervention. Slowness-time correlation (STC) and square root N methods are commonly used slowness calculation methods in the array acoustic data processing process, the slowness calculation methods are processed by using correlation analysis technology and have high precision, but for array acoustic data with seriously insufficient spatial domain (well axial) sampling, the slowness-time correlation method and the square root N method cause precision reduction due to insufficient correlation. In addition, the existing acoustic wave slowness calculation method is difficult to identify the areas with low signal resolution such as thin interbed and interlayer.

Disclosure of Invention

The invention provides a method for calculating slowness of an array sound wave in a laboratory, which is used for solving the problem of slowness precision reduction caused by poor wave train signal-to-noise ratio and serious insufficient spatial domain (well axial) sampling of a miniaturized array sound wave instrument.

The invention adopts the following technical scheme:

the improvement of a method for calculating slowness of array sound waves in a laboratory is that the method comprises the following steps:

step 1, carrying out two-dimensional Fourier transform on the acquired unipolar subarray waveform to obtain a time-space domain waveform X ₀ (z, t) is converted into a frequency-waveform domain X (k, omega) for filtering, and then time-domain formation longitudinal wave information X (z, t) is obtained through two-dimensional inverse Fourier transform:

X(z,t)＝∫∫X(k,ω)·Q(k,ω)e ^i(kz-ωt) dkdω

in the above equation, t represents a time variable, z represents a space variable, k represents a wave number, i is an imaginary number, ω represents an angular vector, and Q (k, ω) represents a filter factor;

step 2, determining the position p of the first arrival and the first jump of each row of waveforms by using a conditional fuzzy C-means clustering algorithm ₀ ：

The objective function of the conditional fuzzy C-means clustering algorithm is as follows:

in the above formula, v _g Is the g-th cluster center, x _h (z, t) is the h sample point, u _gh Is the membership value of the h sample point of the g clustering center, m is a fuzzy coefficient, C is the clustering number, and N is the sample number;

the membership value u is then updated by iteration _gh And a clustering center v _g To solve the objective function:

in the above formula, f _h A condition value representing a fuzzy C-means clustering algorithm, which is defined as:

in the above formula, σ _g Represents the variance, σ, of the cluster g _max Represents σ in all clusters _g M (h), E (h) and R (h) are the absolute average, peak power spectral density and short-term to long-term average ratio, respectively, of a set of sequences d (h) in the X (z, t) waveform:

E(h)＝max(|D(h,ω)| ² )

in the above equation, the constant w represents half the length of the window around h, D (h, ω) is the modulus of the two-dimensional Fourier transform of D (h), and SW and LW are the lengths of the short-term and long-term windows, respectively;

obtaining a membership value capable of representing the wave train data characteristics through a target function, dividing a sample point corresponding to the membership value into a waveform signal class when the membership value is larger than a preset threshold value, and selecting a first component of the waveform signal class as a first arrival p of the wave train data ₀ ；

Step 3, determining the accurate position p of the first arrival by applying Bayesian information criterion ₁ ：

The Bayesian information criterion BIC function is defined as:

BIC(p ₀ )＝p ₀ ln(var{X(1,p ₀ )})+(N-p ₀ -1)ln(var{X(p ₀ +1,N)})-p ₀ ln(N)

in the above formula, X (1,p) ₀ ) Representing array waveform X front p ₀ Vector of data points, X (p) ₀ +1,N) represents a vector of array waveform X remaining data points;

applying BIC function to p ₀ Nearby region [ p ] ₀ -e，p ₀ +e]E represents a constant, and the point with the minimum BIC value is the accurate position p of the first arrival ₁ ；

Step 4, carrying out anomaly detection on the first-arrival points by applying a Principal Component Analysis (PCA) method, and fitting the anomaly points to real first-arrival points through an L1 norm;

in the above formula, Z ₁ Corresponding to the direction of least variance of the raw data, Z ₂ Corresponding to the direction of maximum variance, p, of the raw data _x And p _y Respectively representing the abscissa and ordinate of the first-to-point, p _xo And p _yo Respectively represents Z ₁ And Z ₂ Abscissa and ordinate of the intersection point, θ represents Z ₁ The included angle between the direction and the x direction; when the projection value of the abnormal data projected in the residual subspace is larger than a preset threshold value, judging the abnormal data as abnormal data q;

after the abnormal points are eliminated, the L1 norm is used for obtaining first arrival data at the positions of the abnormal points through straight line fitting, and the L1 norm used for fitting the straight lines is described as follows:

in the above formula, y _l And x _l The relationship is linear:

y(x)＝a+bx

in the above formula, a and b represent coefficients of variable x;

and 5, calculating the high-resolution longitudinal wave slowness by applying a slowness solving formula:

through the first to point p _g Obtaining corresponding head wave arrival time t _g And calculating the high-resolution longitudinal wave slowness by using a slowness formula:

in the above equation, s represents the compressional slowness and RR represents the receiver spacing.

The invention has the beneficial effects that:

the method disclosed by the invention has the advantages of less artificial interference, simple parameter setting and strong generalization capability, and is simultaneously suitable for extracting the longitudinal wave slowness of the conventional digital acoustic wave instrument and array acoustic wave instrument; the unsupervised machine learning method using the condition fuzzy clustering and the BIC information criterion has high calculation efficiency, and can effectively extract the longitudinal wave first arrival of the waveform data with lower signal-to-noise ratio; the principal component analysis method is used for detecting the abnormal points, and the L1 norm is adopted to fit the abnormal points to the real initial point, so that the calculation accuracy of the longitudinal wave slowness can be improved.

Drawings

FIG. 1 is a schematic flow diagram of the process of the present invention;

FIG. 2 is a frequency-waveform domain based filtering;

FIG. 3 is a diagram of first arrival extraction based on a fuzzy clustering algorithm;

FIG. 4 is a determination of a first arrival accurate position based on BIC criteria;

FIG. 5 is a principle of anomaly detection based on the PCA method;

FIG. 6 is anomaly detection based on the PCA method;

fig. 7 is a first arrival correction based on L1 norm fitting.

Detailed Description

In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention is further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.

The miniaturized array acoustic wave instrument in the laboratory has a single-pole transmitting transducer and 4 receiving transducers. The invention aims to provide a method for calculating the slowness of a longitudinal wave measurement signal of an array acoustic wave instrument in a laboratory, and solves the problems that a small array acoustic wave instrument has poor signal-to-noise ratio of a wave train due to the fact that emission interference and first wave first arrival mixing overlap caused by short source distance, and time difference precision is reduced due to low correlation caused by serious insufficient sampling of a spatial domain (well axial direction). In addition, the method can also be used for calculating the longitudinal wave slowness of a conventional array acoustic wave instrument and a digital acoustic wave instrument, and has a good effect on calculating the longitudinal wave slowness in acoustic logging of complex strata (such as a sand-shale thin interbed).

Firstly, preprocessing the acquired monopole array waveform on a frequency-waveform domain, and separating to obtain time domain longitudinal wave information of a stratum; secondly, determining the rough range of the take-off of the primary point of each row of waveforms by using a clustering algorithm; determining the accurate position of the initial point by applying Bayesian information criterion again; then, applying a principal component analysis method to detect abnormal points of the first arrival, and fitting the abnormal points to real first arrival points through an L1 norm; and finally, calculating the high-resolution longitudinal wave slowness by applying a slowness solving formula.

Embodiment 1, this embodiment discloses a method for calculating slowness of array sound waves in a laboratory, as shown in fig. 1, including the following steps:

step 1, the acquired monopole array waveform is processed by two-dimensional Fourier transform, as shown in figure 2, and the time-space domain waveform X is processed ₀ (z, t) is converted into a frequency-waveform domain X (k, omega) for filtering, and then time-domain formation longitudinal wave information X (z, t) is obtained through two-dimensional inverse Fourier transform:

X(z,t)＝∫∫X(k,ω)·Q(k,ω)e ^i(kz-ωt) dkdω

in the above equation, t represents a time variable, z represents a space variable, k represents a wave number, i is an imaginary number, ω represents an angular vector, and Q (k, ω) represents a filter factor;

step 2, as shown in fig. 3, determining the approximate position p of the first arrival and the first jump of each row of waveforms by using a conditional fuzzy C-means clustering algorithm ₀ ：

The objective function of the conditional fuzzy C-means clustering algorithm is as follows:

in the above formula, v _g Is the g-th cluster center, x _h (z, t) is the h sample point, u _gh Is the membership value of the h sample point of the g clustering center, m is a fuzzy coefficient, C is the clustering number, and N is the sample number;

the membership value u is then updated by iteration _gh And a clustering center v _g To solve the objective function:

in the above formula, f _h A condition value representing a fuzzy C-means clustering algorithm, which is defined as:

in the above formula, σ _g Represents the variance, σ, of the cluster g _max Represents σ in all clusters _g M (h), E (h) and R (h) are the absolute average, peak power spectral density and short-term to long-term average ratio, respectively, of a set of sequences d (h) in the X (z, t) waveform:

E(h)＝max(|D(h,ω)| ² )

in the above equation, the constant w represents half the length of the window around h, D (h, ω) is the modulus of the two-dimensional Fourier transform of D (h), and SW and LW are the lengths of the short-term and long-term windows, respectively;

the method comprises the following steps that wave train information of the array sound wave can be divided into 2 types of wave train signals and noise signals, a membership value capable of representing wave train data characteristics is obtained through continuously optimizing a target function according to characteristic difference between the wave train information and the noise signals, when the membership value is larger than a preset threshold value, a sample point corresponding to the membership value is divided into the wave train signal type, and the first component of the wave train signal type is selected as the first arrival p of the wave train data ₀ ；

Step 3, as shown in FIG. 4, determining the accurate position p of the first arrival by applying the Bayesian information criterion ₁ ：

The Bayesian information criterion BIC function is defined as:

BIC(p ₀ )＝p ₀ ln(var{X(1,p ₀ )})+(N-p ₀ -1)ln(var{X(p ₀ +1,N)})-p ₀ ln(N)

in the above formula, X (1,p) ₀ ) Representing array waveform X front p ₀ Vector of data points, X (p) ₀ +1,N) represents a vector of remaining data points of array waveform X;

applying BIC function to p ₀ Nearby region [ p ] ₀ -e，p ₀ +e]E represents a constant, and the point with the minimum BIC value is the accurate position p of the first arrival ₁ ；

Step 4, as shown in fig. 5-6, performing anomaly detection on the first-arrival point by applying a Principal Component Analysis (PCA) method, and fitting the anomaly point to the real first-arrival point through an L1 norm; the principal component analysis technology can effectively find the characteristic direction Z of the most main element in the data ₁ ；

In the above formula, Z ₁ Corresponding to the direction of least variance of the raw data, Z ₂ Corresponding to the direction of maximum variance, p, of the raw data _x And p _y Respectively representing the abscissa and ordinate of the first to point, p _xo And p _yo Respectively represents Z ₁ And Z ₂ Abscissa and ordinate of the intersection point, θ represents Z ₁ The included angle between the direction and the x direction; when the projection value of the abnormal data projected in the residual subspace is larger than a preset threshold value, judging the abnormal data as abnormal data q;

as shown in fig. 7, after the outliers are eliminated, the first arrival data at the location of the outliers is obtained by fitting a straight line using the L1 norm, and the L1 norm for fitting the straight line is described as follows:

in the above formula, y _l And x _l The relationship is linear:

y(x)＝a+bx

in the above formula, a and b represent coefficients of variable x;

and 5, calculating the high-resolution longitudinal wave slowness by applying a slowness solving formula:

through the first to point p _g Obtaining corresponding head wave arrival time t _g And calculating the high-resolution longitudinal wave slowness by using a slowness formula:

in the above equation, s represents the compressional slowness and RR represents the receiver spacing.

Claims (1)

1. A method for calculating slowness of an array sound wave in a laboratory is characterized by comprising the following steps of:

step 1, the acquired monopole array waveform is subjected to two-dimensional Fourier transform, and time-space is processedDomain waveform X ₀ (z, t) is converted into a frequency-waveform domain X (k, omega) for filtering, and then time-domain formation longitudinal wave information X (z, t) is obtained through two-dimensional inverse Fourier transform:

X(z,t)＝∫∫X(k,ω)·Q(k,ω)e ^i(kz-ωt) dkdω

in the above equation, t represents a time variable, z represents a space variable, k represents a wave number, i is an imaginary number, ω represents an angular vector, and Q (k, ω) represents a filter factor;

step 2, determining the position p of the first arrival and the first jump of each row of waveforms by using a conditional fuzzy C-means clustering algorithm ₀ ：

The objective function of the conditional fuzzy C-means clustering algorithm is as follows:

in the above formula, v _g Is the g-th cluster center, x _h (z, t) is the h sample point, u _gh Is the membership value of the h sample point of the g clustering center, m is a fuzzy coefficient, C is the clustering number, and N is the sample number;

the membership value u is then updated by iteration _gh And a clustering center v _g To solve the objective function:

in the above formula, f _h A condition value representing a fuzzy C-means clustering algorithm, which is defined as:

in the above formula, σ _g Represents the variance, σ, of the cluster g _max Represents σ in all clusters _g M (h), E (h) and R (h) are the absolute average, peak power spectral density and short-term to long-term average ratio, respectively, of a set of sequences d (h) in the X (z, t) waveform:

E(h)＝max(|D(h,ω)| ² )

in the above equation, the constant w represents half the length of the window around h, D (h, ω) is the modulus of the two-dimensional Fourier transform of D (h), and SW and LW are the lengths of the short-term and long-term windows, respectively;

obtaining a membership value capable of representing the wave train data characteristics through a target function, dividing a sample point corresponding to the membership value into a waveform signal class when the membership value is larger than a preset threshold value, and selecting a first component of the waveform signal class as a first arrival p of the wave train data ₀ ；

Step 3, determining the accurate position p of the first arrival by applying Bayesian information criterion ₁ ：

The Bayesian information criterion BIC function is defined as:

BIC(p ₀ )＝p ₀ ln(var{X(1,p ₀ )})+(N-p ₀ -1)ln(var{X(p ₀ +1,N)})-p ₀ ln(N)

in the above formula, X (1,p) ₀ ) Representing the array waveform Xfront p ₀ Vector of data points, X (p) ₀ +1,N) represents a vector of remaining data points of array waveform X;

applying BIC function to p ₀ Nearby region [ p ] ₀ -e，p ₀ +e]E represents a constant, and the point with the minimum BIC value is the accurate position p of the first arrival ₁ ；

Step 4, carrying out anomaly detection on the first-arrival points by applying a Principal Component Analysis (PCA) method, and fitting the anomaly points to real first-arrival points through an L1 norm;

in the above formula, Z ₁ Corresponding to the direction of least variance of the raw data, Z ₂ Corresponding to the direction of maximum variance, p, of the raw data _x And p _y Respectively representing the abscissa and ordinate of the first to point, p _xo And p _yo Respectively represent Z ₁ And Z ₂ Abscissa and ordinate of the intersection point, θ represents Z ₁ The included angle between the direction and the x direction; when the projection value of the abnormal data projected in the residual subspace is larger than a preset threshold value, judging the abnormal data to be abnormal data q;

after the abnormal points are eliminated, the L1 norm is used for obtaining first arrival data at the positions of the abnormal points through straight line fitting, and the L1 norm used for fitting the straight lines is described as follows:

in the above formula, y _l And x _l The relationship is linear:

y(x)＝a+bx

in the above formula, a and b represent coefficients of variable x;

and 5, calculating the high-resolution longitudinal wave slowness by applying a slowness solving formula:

through the first to point p _g Obtaining corresponding head wave arrival time t _g And calculating the high-resolution longitudinal wave slowness by using a slowness formula:

in the above equation, s represents the compressional slowness and RR represents the receiver spacing.

Images