A method for calculating acoustic slowness of array in laboratory

Abstract

The present invention discloses a method for calculating array acoustic wave slowness in a laboratory, comprising the following steps: Step 1, subjecting the collected monopole array waveform to a two-dimensional Fourier transform to obtain time-domain stratum longitudinal wave information; Step 2, determining the starting position of the first arrival of each column of waveform using a conditional fuzzy C-means clustering algorithm; Step 3, determining the exact position of the first arrival using the Bayesian information criterion; Step 4, performing anomaly detection on the first arrival point using a principal component analysis (PCA) method, and fitting the anomaly point to the true first arrival point using a norm; Step 5, calculating high-resolution longitudinal wave slowness using a slowness solution formula. The method disclosed by the present invention has little human intervention, simple parameter setting, and strong generalization ability, and is applicable to longitudinal wave slowness extraction of conventional digital acoustic wave instruments and array acoustic wave instruments; and the unsupervised machine learning method using conditional fuzzy clustering and the BIC information criterion has high computational efficiency.

Description

Technical Field

[0001] The invention belongs to the field of slowness extraction of longitudinal wave measurement signals of miniaturized array acoustic wave instruments in laboratories, and particularly relates to a method for calculating array acoustic wave slowness in laboratories in this field. Background Art

[0002] Acoustic slowness is the reciprocal of the propagation speed of an acoustic signal in a formation. As an essential component of the nine conventional logging curves, the acoustic slowness curve plays a crucial role in well logging measurements and evaluation. It can be used to identify lithology, analyze wellbore stability, calculate formation porosity, estimate formation permeability, and evaluate formation anisotropy.

[0003] Currently, acoustic slowness in well logging is primarily determined by correlating the arrival time of the detection wave with the slowness. Calculation methods can be categorized into two types: frequency-domain processing and time-domain processing. Frequency-domain processing, based on different principles, can be further categorized into methods such as Prony prediction and weighted spectral correlation. However, these methods suffer from poor noise immunity and are generally used as theoretical analysis tools without practical calculations. Time-domain processing methods primarily include threshold methods, long- and short-time window energy ratio methods, and waveform prediction methods. Among waveform prediction methods, the slowness-time correlation (STC) method and the Nth-root method have been studied. Due to the poor noise immunity of the threshold and short- and long-time window energy ratio methods, they suffer from severe distortion when processing low signal-to-noise ratio data, requiring manual intervention to achieve high accuracy. The slowness-time correlation (STC) method and the Nth-root method are commonly used slowness calculation methods in array acoustic data processing. These methods utilize correlation analysis techniques and offer high accuracy. However, for array acoustic data with severely undersampled spatial domain (wellbore axis), these methods suffer from reduced accuracy due to insufficient correlation. In addition, existing acoustic slowness calculation methods have difficulty identifying areas with low signal resolution, such as thin interlayers and interlayers. Summary of the Invention

[0004] The present invention provides a method for calculating array acoustic wave slowness in a laboratory, which is used to solve the problem of decreased slowness accuracy of miniaturized array acoustic wave instruments due to poor wave train signal-to-noise ratio and serious lack of spatial domain (well axial) sampling.

[0005] The present invention adopts the following technical solutions:

[0006] A method for calculating array acoustic wave slowness in a laboratory, the improvement of which is that it comprises the following steps:

[0007] Step 1: The collected monopole array waveform is transformed into the frequency-waveform domain X(k, ω) through a two-dimensional Fourier transform, and then filtered. Then, the time-domain formation compressional wave information X(z, t) is obtained through a two-dimensional inverse Fourier transform:

[0008] X(z,t)=∫∫X(k,ω)·Q(k,ω)e i(kz-ωt) dkdω

[0009] In the above formula, t represents the time variable, z represents the space variable, k represents the wave number, i is an imaginary number, ω represents the angle vector, and Q(k,ω) represents the filter factor;

[0010] Step 2: Use the conditional fuzzy C-means clustering algorithm to determine the starting position p0 of each column waveform:

[0011] The objective function of the conditional fuzzy C-means clustering algorithm is:

[0012]

[0013] In the above formula, v g is the g-th cluster center, x h (z,t) is the hth sample point, u gh is the membership value of the hth sample point of the gth cluster center, m is the fuzzy coefficient, C is the number of clusters, and N is the number of samples;

[0014] Then the membership value u is updated iteratively gh and cluster center v g To solve the objective function:

[0015]

[0016]

[0017] In the above formula, f h represents the conditional value of the fuzzy C-means clustering algorithm, which is defined as:

[0018]

[0019]

[0020] In the above formula, σ g represents the variance of cluster g, σ max Indicates the σ in all clusters g The maximum value of M(h), E(h) and R(h) are the absolute mean, peak power spectral density and short-term to long-term average ratio of a set of sequences d(h) in the X(z,t) waveform respectively:

[0021]

[0022] E(h)=max(|D(h,ω)| 2 )

[0023]

[0024] In the above formula, 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), SW and LW are the lengths of the short-term and long-term windows respectively;

[0025] The membership value that can characterize the characteristics of the wave train data is obtained through the objective function. When the membership value is greater than the preset threshold, the sample point corresponding to the value is divided into the waveform signal class, and the first component of the waveform signal class is selected as the initial arrival p0 of the wave train data;

[0026] Step 3: Apply the Bayesian Information Criterion to determine the exact position p1 of the first arrival:

[0027] The Bayesian Information Criterion BIC function is defined as:

[0028] BIC(p0)=p0ln(var{X(1,p0)})+(N-p0-1)ln(var{X(p0+1,N)})-p0ln(N)

[0029] In the above formula, X(1,p0) represents the vector consisting of the first p0 data points of the array waveform X, and X(p0+1,N) represents the vector consisting of the remaining data points of the array waveform X;

[0030] Apply the BIC function to the interval [p0-e, p0+e] near p0, where e represents a constant. The point with the minimum BIC value is the exact position of the first arrival, p1.

[0031] Step 4: Apply the principal component analysis (PCA) method to detect anomalies of the first arrival points and fit the anomaly points to the true first arrival points using the L1 norm.

[0032]

[0033] In the above formula, Z1 corresponds to the direction with the smallest variance of the original data, Z2 corresponds to the direction with the largest variance of the original data, and p x and p y They represent the horizontal and vertical coordinates of the initial arrival point, respectively, p xo and p yo represent the horizontal and vertical coordinates of the intersection of Z1 and Z2, respectively, and θ represents the angle between the Z1 direction and the x direction; when the projection value of the abnormal data in the residual subspace is greater than the preset threshold, it is judged as abnormal data q;

[0034] After removing the outliers, the L1 norm is used to obtain the first arrival data at the outlier position through straight line fitting. The L1 norm used for fitting the straight line is described as follows:

[0035]

[0036] In the above formula, y l with x l There is a linear relationship between them:

[0037] y(x)=a+bx

[0038] In the above formula, a and b represent the coefficients of variable x;

[0039] Step 5: Apply the slowness solution formula to calculate the high-resolution P-wave slowness:

[0040] Through the initial point p g Get the corresponding first wave arrival time t g , use the slowness formula to calculate the high-resolution longitudinal wave slowness:

[0041]

[0042] In the above formula, s represents the longitudinal wave slowness and RR represents the receiver spacing.

[0043] The beneficial effects of the present invention are:

[0044] The method disclosed in the present invention has little human intervention, simple parameter setting, and strong generalization ability, and is applicable to the extraction of longitudinal wave slowness from conventional digital acoustic wave instruments and array acoustic wave instruments. The unsupervised machine learning method using conditional fuzzy clustering and BIC information criterion has high computational efficiency and can effectively extract the first arrival of longitudinal waves from waveform data with a low signal-to-noise ratio. The principal component analysis method is used to detect outliers, and the L1 norm is used to fit the outliers to the true first arrival points, which can improve the calculation accuracy of the longitudinal wave slowness. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] Figure 1 It is a schematic flow diagram of the method of the present invention;

[0046] Figure 2 It is filtering based on frequency-waveform domain;

[0047] Figure 3 It is based on fuzzy clustering algorithm to extract the first arrival;

[0048] Figure 4 The accurate position of the first arrival is determined based on the BIC criterion;

[0049] Figure 5 It is an anomaly detection principle based on the PCA method;

[0050] Figure 6 It is an outlier detection based on PCA method;

[0051] Figure 7 It is a first arrival correction based on L1 norm fitting. DETAILED DESCRIPTION

[0052] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention.

[0053] The miniaturized array acoustic instrument used in the laboratory has a monopole transmitting transducer and four receiving transducers. The purpose of the present invention is to provide a method for calculating the slowness of the P-wave measurement signal of the laboratory array acoustic instrument. This method addresses the problems of the miniaturized array acoustic instrument, such as the short source distance causing the transmission interference to overlap with the first arrival of the first wave, resulting in a poor signal-to-noise ratio of the wave train, and the serious lack of spatial domain (well axial) sampling, resulting in low correlation and reduced time difference accuracy. Furthermore, the present invention can also be used for calculating P-wave slowness in conventional array acoustic instruments and digital acoustic instruments, and is also effective in calculating P-wave slowness in acoustic logging of complex formations (such as thin interbedded sandstone and mudstone).

[0054] The present invention first preprocesses the collected monopole array waveforms in the frequency-waveform domain to separate the time-domain P-wave information of the formation; secondly, a clustering algorithm is used to determine the approximate range of the first arrival point of each column of waveforms; thirdly, the Bayesian information criterion is applied to determine the precise location of the first arrival point; then, the principal component analysis method is used to detect outliers on the first arrivals, and the outliers are fitted to the true first arrival points through the L1 norm; finally, the slowness solution formula is used to calculate the high-resolution P-wave slowness.

[0055] Example 1: This example discloses a method for calculating the slowness of an array acoustic wave in a laboratory. Figure 1 As shown, the following steps are included:

[0056] Step 1: The collected monopole array waveform is transformed by two-dimensional Fourier transform, as shown in Figure 2 As shown in Figure 1, the time-space domain waveform X0(z,t) is converted to the frequency-waveform domain X(k,ω) for filtering, and then the time-domain formation compressional wave information X(z,t) is obtained through a two-dimensional inverse Fourier transform:

[0057] X(z,t)=∫∫X(k,ω)·Q(k,ω)e i(kz-ωt) dkdω

[0058] In the above formula, t represents the time variable, z represents the space variable, k represents the wave number, i is an imaginary number, ω represents the angle vector, and Q(k,ω) represents the filter factor;

[0059] Step 2, such as Figure 3 As shown in the figure, the conditional fuzzy C-means clustering algorithm is used to determine the approximate position p0 of the first arrival of each column of waveform:

[0060] The objective function of the conditional fuzzy C-means clustering algorithm is:

[0061]

[0062] In the above formula, v g is the gth cluster center, x h (z,t) is the hth sample point, u gh is the membership value of the hth sample point of the gth cluster center, m is the fuzzy coefficient, C is the number of clusters, and N is the number of samples;

[0063] Then the membership value u is updated iteratively gh and cluster center v g To solve the objective function:

[0064]

[0065]

[0066] In the above formula, f h represents the conditional value of the fuzzy C-means clustering algorithm, which is defined as:

[0067]

[0068]

[0069] In the above formula, σ g represents the variance of cluster g, σ max Indicates the σ in all clusters g The maximum value of M(h), E(h) and R(h) are the absolute mean, peak power spectral density and short-term to long-term average ratio of a set of sequences d(h) in the X(z,t) waveform respectively:

[0070]

[0071] E(h)=max(|D(h,ω)| 2 )

[0072]

[0073] In the above formula, 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), SW and LW are the lengths of the short-term and long-term windows respectively;

[0074] The wave train information of the array acoustic wave can be divided into two categories: waveform signals and noise signals. In this step, based on the characteristic differences between the two, the objective function is continuously optimized to obtain the membership value that can characterize the wave train data characteristics. When the membership value is greater than a preset threshold, the sample point corresponding to the value is classified as the waveform signal class, and the first component of the waveform signal class is selected as the first arrival p0 of the wave train data.

[0075] Step 3, such as Figure 4 As shown, the Bayesian Information Criterion is applied to determine the exact position p1 of the first arrival:

[0076] The Bayesian Information Criterion BIC function is defined as:

[0077] BIC(p0)=p0ln(var{X(1,p0)})+(N-p0-1)ln(var{X(p0+1,N)})-p0ln(N)

[0078] In the above formula, X(1,p0) represents the vector consisting of the first p0 data points of the array waveform X, and X(p0+1,N) represents the vector consisting of the remaining data points of the array waveform X;

[0079] Apply the BIC function to the interval [p0-e, p0+e] near p0, where e represents a constant. The point with the minimum BIC value is the exact position of the first arrival, p1.

[0080] Step 4, such as Figure 5 As shown in Figure 6, the principal component analysis (PCA) method is used to detect anomalies of the initial arrival points, and the anomaly points are fitted to the true initial arrival points through the L1 norm. The principal component analysis technology can effectively find the most important element characteristic direction Z1 in the data.

[0081]

[0082] In the above formula, Z1 corresponds to the direction with the smallest variance of the original data, Z2 corresponds to the direction with the largest variance of the original data, and p x and p y They represent the horizontal and vertical coordinates of the initial arrival point, respectively, p xo and p yo represent the horizontal and vertical coordinates of the intersection of Z1 and Z2, respectively, and θ represents the angle between the Z1 direction and the x direction; when the projection value of the abnormal data in the residual subspace is greater than the preset threshold, it is judged as abnormal data q;

[0083] like Figure 7 As shown in the figure, after removing the outliers, the L1 norm is used to obtain the first arrival data at the outlier position through straight line fitting. The L1 norm used for fitting the straight line is described as follows:

[0084]

[0085] In the above formula, y l with x l There is a linear relationship between them:

[0086] y(x)=a+bx

[0087] In the above formula, a and b represent the coefficients of variable x;

[0088] Step 5: Apply the slowness solution formula to calculate the high-resolution P-wave slowness:

[0089] Through the initial point p g Get the corresponding first wave arrival time t g , use the slowness formula to calculate the high-resolution longitudinal wave slowness:

[0090]

[0091] In the above formula, s represents the longitudinal wave slowness and RR represents the receiver spacing.

Claims

1. A method for calculating array acoustic wave slowness in a laboratory, characterized in that: The steps include: Step 1: The collected monopole array waveform is transformed into the frequency-waveform domain X(k, ω) through a two-dimensional Fourier transform, and then filtered. Then, the time-domain formation compressional wave information X(z, t) is obtained through a two-dimensional inverse Fourier transform: X(z,t)=∫∫X(k,ω)·Q(k,ω)e i(kz-ωt) dkdω In the above formula, t represents the time variable, z represents the space variable, k represents the wave number, i is an imaginary number, ω represents the angle vector, and Q(k,ω) represents the filter factor; Step 2: Use the conditional fuzzy C-means clustering algorithm to determine the starting position p0 of each column waveform: The objective function of the conditional fuzzy C-means clustering algorithm is: In the above formula, v g is the g-th cluster center, x h (z,t) is the hth sample point, u gh is the membership value of the hth sample point of the gth cluster center, m is the fuzzy coefficient, C is the number of clusters, and N is the number of samples; Then the membership value u is updated iteratively gh and cluster center v g To solve the objective function: In the above formula, f h represents the conditional value of the fuzzy C-means clustering algorithm, which is defined as: In the above formula, σ g represents the variance of cluster g, σ max Indicates the σ in all clusters g The maximum value of M(h), E(h) and R(h) are the absolute mean, peak power spectral density and short-term to long-term average ratio of a set of sequences d(h) in the X(z,t) waveform respectively: E(h)=max(|D(h,ω)| 2 ) In the above formula, d h represents the central element of the hth sequence d(h) in the X(z,t) waveform, 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), SW and LW are the lengths of the short-term and long-term windows respectively, d j represents the central element d of the hth group sequence h The j-th element value in the constructed short-term or long-term window. The membership value that can characterize the characteristics of the wave train data is obtained through the objective function. When the membership value is greater than the preset threshold, the sample point corresponding to the value is divided into the waveform signal class, and the first component of the waveform signal class is selected as the initial arrival p0 of the wave train data; Step 3: Apply the Bayesian Information Criterion to determine the exact position p1 of the first arrival: The Bayesian Information Criterion BIC function is defined as: BIC(p0)=p0ln(var{X(1,p0)})+(N-p0-1)ln(var{X(p0+1,N)})-p0ln(N) In the above formula, X(1,p0) represents the vector consisting of the first p0 data points of the array waveform X, and X(p0+1,N) represents the vector consisting of the remaining data points of the array waveform X; Apply the BIC function to the interval [p0-e, p0+e] near p0, where e represents a constant. The point with the minimum BIC value is the exact position of the first arrival, p1. Step 4: Apply the principal component analysis (PCA) method to detect anomalies of the first arrival points and fit the anomaly points to the true first arrival points using the L1 norm. In the above formula, Z1 corresponds to the direction with the smallest variance of the original data, Z2 corresponds to the direction with the largest variance of the original data, and p x and p y They represent the horizontal and vertical coordinates of the initial arrival point, respectively, p xo and p yo represent the horizontal and vertical coordinates of the intersection of Z1 and Z2, respectively, and θ represents the angle between the Z1 direction and the x direction; when the projection value of the abnormal data in the residual subspace is greater than the preset threshold, it is judged as abnormal data q; After removing the outliers, the L1 norm is used to obtain the first arrival data at the outlier position through straight line fitting. The L1 norm used for fitting the straight line is described as follows: In the above formula, y l with x l There is a linear relationship between them: y(x)=a+bx In the above formula, a and b represent the coefficients of variable x; Step 5: Apply the slowness solution formula to calculate the high-resolution P-wave slowness: Through the initial point p g Get the corresponding first wave arrival time t g , use the slowness formula to calculate the high-resolution longitudinal wave slowness: In the above formula, s represents the longitudinal wave slowness and RR represents the receiver spacing.

Images