Lithology identification-oriented while-drilling signal optimization reconstruction method

Abstract

The application discloses a lithology identification-oriented while-drilling signal optimization reconstruction method, which is used for processing vibration signals during drilling rock samples by a drill bit and belongs to the technical field of geological drilling. The method comprises the following steps: obtaining vibration signals; obtaining a frequency spectrum of the vibration signals through Fourier transform; adaptively generating a frequency band boundary on the frequency spectrum; constructing an initial matrix A and a Hankel matrix; and sequentially reconstructing the Hankel matrix, the initial matrix and the signal after singular value decomposition and singular value screening of the Hankel matrix. The method improves the empirical wavelet transform method, avoids the problem of inaccurate frequency band division in the empirical wavelet transform, accurately extracts an intrinsic mode function, ensures that the function can better capture the local characteristics of the signal at different frequencies, makes the reconstructed signal have stronger feature reservation capability and noise resistance, and comprehensively captures time domain, frequency domain and entropy domain information of the signal, so that efficient representation of the vibration signal is realized.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of geological drilling, and particularly relates to a lithology identification-oriented while-drilling signal optimization reconstruction method. BACKGROUND

[0002] Lithology identification has an important influence on drilling efficiency and cost-effectiveness of resource development. Although traditional lithology identification methods such as rock chip logging and drilling coring are widely used in actual operation, these methods rely on a large amount of manpower and material resources, and have a large lag in data acquisition, which is difficult to meet the needs of modern fast-paced drilling operations.

[0003] While-drilling lithology identification technology is different from traditional methods. It collects drilling rig drilling parameter data in real time and applies deep learning algorithms for interpretation to output the strength and structure information of the rock formation. According to the different types of signals, while-drilling lithology identification research mainly divides into three methods: mechanical parameter-based, acoustic parameter-based and vibration parameter-based. The mechanical parameter-based method monitors the running state of the drilling rig, but is easily affected by the noise of the equipment itself, resulting in a decrease in identification accuracy. The acoustic parameter-based method uses acoustic signals generated during drilling. However, in a high-noise environment, acoustic signals are easily disturbed, affecting the identification effect. The vibration parameter-based method identifies lithology by analyzing the vibration characteristics of the drill pipe. Compared with mechanical signals, vibration signals are not easily affected by equipment noise; compared with acoustic signals, vibration signals have more significant differences in distinguishing rock properties and have better potential. Although vibration signals perform well in lithology identification, their processing usually needs to be converted into time-frequency images. However, vibration signals have time-varying and non-stationary characteristics, making feature extraction and classification complex. At the same time, in the case of small samples, it will lead to insufficient training or overfitting of the deep learning model, affecting the classification accuracy. Therefore, researching efficient feature extraction methods and deep learning models suitable for small samples is of great significance for dealing with the time-varying and non-stationary problems of while-drilling vibration signals and improving the accuracy of lithology identification.

[0004] In recent years, while-drilling lithology identification technology has been increasingly widely applied in geological engineering. Researchers have proposed various innovative methods in signal processing and feature extraction to improve the feature extraction accuracy and classification accuracy of vibration signals. Current research on while-drilling lithology identification mainly focuses on improving signal processing methods, applying deep learning technology for classification, multi-domain feature fusion and dealing with the time-varying and non-stationary characteristics of vibration signals. These research results not only improve the accuracy and real-time performance of lithology identification, but also lay a foundation for the further application of while-drilling vibration signal identification technology in geological engineering. However, there are still deficiencies in dealing with the time-varying and non-stationary problems of while-drilling vibration signals.

[0005] In view of the time-varying non-stationary problem of the while-drilling vibration signal, scholars have proposed a variety of solutions. Although these methods have made certain progress in dealing with time-varying non-stationary signals, there are still limitations such as insufficient frequency domain segmentation and inaccurate signal component extraction, which affect the effectiveness of the subsequent feature reconstruction process.

[0006] Feature extraction and feature reconstruction are key steps in vibration signal analysis. The main purpose of feature reconstruction is to further strengthen the useful features of the signal and reduce noise and redundant information, thereby improving the classification accuracy and robustness of the model. Current feature reconstruction research mainly focuses on multi-modal reconstruction based on classification tasks, which combines different feature domain information to improve model performance. Existing methods for identifying lithology changes mainly analyze the time domain and frequency domain features of the vibration signal, but the time domain features are easily affected by external interference, and the frequency domain features cannot directly reflect the vibration intensity. For the combination of time domain and frequency domain, the amplitude variation degree and frequency distribution cannot be quantified. Information entropy is a measure of signal complexity and disorder, which can quantify the difference of a single component or reconstructed signal, and can provide an effective supplement to time-frequency analysis methods. However, existing multi-modal reconstruction methods mainly focus on time domain and frequency domain features, ignoring the entropy domain features in the vibration signal.

[0007] Empirical Wavelet Transform (EWT) is a signal processing method that combines wavelet analysis and Empirical Mode Decomposition (EMD). Wavelet transform achieves multi-resolution analysis of signals by translating and scaling wavelet basis functions at different scales and positions. Empirical Mode Decomposition can adaptively decompose complex signals into a finite number of intrinsic mode functions, also known as Intrinsic Mode Function (IMF) or IMF components. Unlike traditional wavelet transform, which uses predefined basis functions, EWT adaptively constructs wavelet filters based on signal spectral characteristics, providing more accurate time-frequency analysis. EWT has shown broad application prospects in fault diagnosis, biomedical signal processing, image analysis, and other fields. The method flow is as follows Figure 1: (1) The original signal is transformed into frequency domain by FFT to obtain a frequency spectrum; (2) Band division: the band boundaries are automatically divided by detecting the local maximum value of the Fourier spectrum of the signal, avoiding manual selection of the wavelet basis and improving adaptability; (3) Constructing an empirical wavelet filter bank: constructing a set of orthogonal wavelet filters (empirical wavelet function, EWF) with a tight support on each frequency band to form a set of band-pass filters, which are close to ideal filters; (4) Obtaining component signals: decomposing the signal into a plurality of empirical mode functions (EMF) by the filter bank, each EMF representing the time-frequency characteristics of a frequency band. The signal components obtained in this step can also be regarded as obtaining a plurality of intrinsic mode functions.

[0008] Singular Value Decomposition (SVD) is an important matrix decomposition in linear algebra, which can decompose any real or complex matrix (whether it is a square matrix or not, whether it is full rank or not) into the product of three simple matrices, thereby revealing the internal structure, key features and energy distribution of the original matrix. For any matrix A ∈ R m×n , it can be decomposed into A = UΣV ⊤ , where U is the left singular matrix, Σ is the diagonal matrix, and V is the right singular matrix. The functions of SVD include (1) feature dimension reduction: reducing the computational load and visualizing high-dimensional data; (2) data reconstruction: such as reconstructing signals and images (lossy compression can be achieved, and the smaller k is, the higher the compression rate is, but the greater the image quality loss is); (3) noise reduction, etc. SUMMARY

[0009] The purpose of the present application is to provide a lithology identification-oriented while-drilling signal optimization reconstruction method that solves the problems of inaccurate downstream tasks caused by inaccurate band division and incomplete signal feature extraction in empirical wavelet transform.

[0010] To achieve the above purpose, the technical scheme adopted by the present application is as follows: a lithology identification-oriented while-drilling signal optimization reconstruction method, wherein the while-drilling signal is a vibration signal of a rock sample in a drilling area drilled by a drill bit, and the method comprises the following steps:

[0011] S1, obtaining a vibration signal x, and obtaining the frequency spectrum f(x) of x by Fourier transform, wherein the frequency spectrum length is p1;

[0012] S2, adaptively generating a group of band boundaries on f(x);

[0013] S3, processing the band boundaries based on the empirical wavelet transform method to obtain N1 intrinsic mode functions imf1~imf N1 ;

[0014] S4, Construct the initial matrix A, and construct the Hankel matrix H based on A. n×k1 ;

[0015] A=(imf1,imf2,…imf N-1 ,imf N ) T T represents the transpose operation;

[0016] ;

[0017] In the formula, n and k1 are respectively H n×k1 The number of rows and columns, 1≤n≤k1, k1=N1-n+1;

[0018] S5, for H n×k1 Perform singular value decomposition and select the first k2 singular values ​​to reconstruct the Hankel matrix H';

[0019] S6, Construct H based on the initial matrix A n×k1 The method is to perform an inverse transformation to obtain the initial reconstruction matrix A' from H', and then use each row element of A' as an intrinsic mode function to reconstruct the signal, thus obtaining the reconstructed signal x' of x.

[0020] As a preferred option, S2 specifically includes Sa1~Sa3;

[0021] Sa1, with a preset step size SL, where SL is the Nyquist frequency, and f(x) is divided into p1 segments L1~L according to the following formula. p1 The j-th segment is L j The maximum amplitude is A j , ;

[0022] Sa2, preset amplitude threshold ∂, select A j The segments of >∂ are arranged in descending order of their maximum amplitude values ​​to form sequence L1;

[0023] Sa3, for L1, generates a frequency band boundary by grouping every three adjacent segments, where one frequency band boundary [ω k ,ω h The generation method is as follows: obtain the frequencies corresponding to the maximum amplitude values ​​of a set of segments, and label them sequentially as ω. a ω b ω c ω k =(ω a +ω b ) / 2, ω h =(ω b +ω c ) / 2.

[0024] As preferred, S5 specifically comprises S51-S53:

[0025] S51, according to the following formula, H n×k1 is singular value decomposed;

[0026] ;

[0027] In the formula, U n×n is an n x n left singular matrix, S is a diagonal matrix, V k1×k1 is a k1 x k1 right singular matrix, and n singular values are contained in S;

[0028] S52, by accumulating the contribution rate, the first k2 singular values greater than 0 and their corresponding left singular matrix and right singular matrix are retained to obtain a reconstructed Hankel matrix H'.

[0029] As preferred, the following steps are further included:

[0030] S7, Fourier transform of x' is performed to obtain its frequency spectrum f(x'), and the frequency spectrum length is p2;

[0031] S8, f(x') is processed according to S2-S5 to obtain the intrinsic mode function corresponding to x' and the singular value, wherein the intrinsic mode function is N2, and in sequence ~ , and the singular value is k3, and in sequence σ1-σ k3 ;

[0032] S9, based on x', the intrinsic mode function and the singular value of x', a plurality of time domain features, frequency domain features and entropy domain features are extracted, and all the features constitute a multi-domain feature fusion vector F.

[0033] As preferred, the entropy domain features include power spectrum entropy H eE , singular spectrum entropy H svdpE , and energy entropy H psdE , which are obtained according to the following formulas respectively.

[0034] ,

[0035] ,

[0036] ,

[0037] ,

[0038] In the formula, y m is the value of f(x') when the frequency is m, 1≤m≤P2, σ z is the zth singular value of x', k3 is the total number of singular values corresponding to x', 1≤z≤k3, is the qth intrinsic mode function of x' E is the energy value of x' at time t q is the energy value of x' at time t is the energy value of x' at time t

[0039] Preferably, in S9, the time domain features include the maximum value, minimum value, mean value, variance, standard deviation, kurtosis, skewness, root mean square value, waveform factor, peak factor, impact factor and pulse gap factor of the reconstructed signal x';

[0040] The frequency domain features include the center frequency, mean square frequency, root mean square frequency, frequency variance and frequency standard deviation of the spectrum f(x');

[0041] The entropy domain features include power spectrum entropy, singular spectrum entropy and energy entropy.

[0042] Compared with the prior art, the present application has the following advantages:

[0043] (1) The empirical wavelet transform method is improved to more accurately divide the frequency band: the present application proposes a new method for adaptively generating frequency band boundaries, which replaces the original frequency band boundary generation method of the empirical wavelet transform method with the improved method in steps Sa1-Sa3. According to steps Sa1-Sa3, the present application generates and filters segments, and generates a frequency band boundary for every three segments as an adaptively generated window. The window can dynamically adjust the size and shape of the window according to the spectral characteristics of the signal, so that the empirical wavelet transform can more accurately divide the appropriate spectral boundaries, avoiding the problem of inaccurate frequency band division in the empirical wavelet transform, and ensuring that the intrinsic mode functions can better capture the local characteristics of the signal at different frequencies.

[0044] (2) Improve the feature extraction capability of the reconstructed signal: in order to reduce the influence of invalid components in the intrinsic mode function on the subsequent lithology classification task, the adaptive empirical wavelet transform is integrated with the singular value decomposition method, which aims to filter out components that have a significant influence on the main structure and dynamic characteristics of the vibration signal from the intrinsic mode function, thereby improving the feature extraction capability of the reconstructed signal.

[0045] (3) A multi-domain feature fusion vector F is constructed: the present application introduces entropy domain feature analysis to comprehensively capture the time domain, frequency domain and entropy domain information of the signal and construct a multi-domain feature fusion vector F, which realizes efficient representation of the time-frequency characteristics of the vibration signal, and in the feature extraction and classification experiments in the embodiments of the present application, it shows a higher classification accuracy, especially in distinguishing complex categories, with a lower error rate. BRIEF DESCRIPTION OF DRAWINGS

[0046] Figure 1 is the flowchart of the empirical wavelet transform in the prior art;

[0047] Figure 2 Flow chart for generating a reconstructed signal of the present application;

[0048] Figure 3 Schematic diagram for generating a band boundary with three adjacent segments in step Sa3;

[0049] Figure 4 Schematic diagram for the empirical mode decomposition of a vibration signal;

[0050] Figure 5 Schematic diagram for the empirical wavelet transform of a vibration signal;

[0051] Figure 6 Schematic diagram for the empirical wavelet transform autocorrelation analysis of a vibration signal;

[0052] Figure 7 Schematic diagram for the adaptive empirical wavelet transform of a vibration signal according to the present application. DETAILED DESCRIPTION

[0053] The present application will be further described below in conjunction with embodiments and drawings.

[0054] Embodiment 1: see Figures 1-3 A lithology-identifying while-drilling signal optimized reconstruction method, the while-drilling signal being a vibration signal of a drill bit when drilling a rock sample in a work area, comprising the following steps:

[0055] S1, obtaining a vibration signal x, and obtaining its frequency spectrum f(x) by Fourier transform, the frequency spectrum length being p1;

[0056] S2, adaptively generating a set of band boundaries on f(x);

[0057] S3, processing the band boundaries based on an empirical wavelet transform method to obtain N1 intrinsic mode functions imf1~imf N1 ;

[0058] S4, constructing an initial matrix A, and constructing a Hankel matrix H n×k1 based on A;

[0059] A=(imf1, imf2, …imf N-1 , imf N ) T , T being a transpose operation;

[0060] ;

[0061] In the formula, n and k1 are the row number and column number of H n×k1 , respectively, 1≤n≤k1, k1=N1-n+1;

[0062] S5, performing singular value decomposition on H n×k1Singular value decomposition is performed, and the first k2 singular values are selected to reconstruct the Hankel matrix H';

[0063] S6, constructing H from the initial matrix A n×k1 , and the reconstructed signal x' of x is obtained by taking each row of A' as an intrinsic mode function to reconstruct the signal.

[0064] S2 specifically includes Sa1-Sa3;

[0065] Sa1, a preset step size SL, the step size SL is the Nyquist frequency, and f(x) is divided into p1 segments L1-L p1 , according to the following formula j , where the jth segment is L j , the amplitude maximum value is A j , ;

[0066] Sa2, a preset amplitude threshold ∂, selecting segments with A >∂, and arranging them in descending order of amplitude maximum value to obtain a sequence L1;

[0067] Sa3, generating a frequency band boundary for every three adjacent segments in L1, where the generation method of a frequency band boundary [ω k ,ω h ] is as follows: obtaining the frequencies corresponding to the amplitude maximum values of a group of segments, and sequentially marking them as ω a ,ω b ,ω c ,ω k =(ω a +ω b ) / 2,ω h =(ω b +ω c ) / 2. Taking Figure 3 for example, it is assumed that a spectrum f(x') is as shown in Figure 3 , for ease of distinction, the overall background is represented as a shaded area, and the obtained segments are displayed in white areas, and it can be seen that the three white intervals in Figure 3 correspond to ω a ,ω b ,ω c ,ω k ,ω h is located between ω a ,ω b ,ω c .

[0068] S5 specifically includes S51-S53:

[0069] S51, singular value decomposition is performed on H n×k1 according to the following formula;

[0070] ;

[0071] wherein U n×n is an n x n left singular matrix, S is a diagonal matrix, V k1×k1 is a k1 x k1 right singular matrix, S contains n singular values;

[0072] S52, by accumulating the contribution rate, the first k2 singular values greater than 0 and the corresponding left singular matrix and right singular matrix are retained to obtain a reconstructed Hankel matrix H'.

[0073] Embodiment 2: refer to Figures 1-3 On the basis of embodiment 1, in addition to obtaining the reconstructed signal x' of x, the reconstructed signal x' can be further feature extracted to generate fusion features for downstream tasks. For example, on the basis of embodiment, in addition to steps S1-S6, the following steps are included:

[0074] S7, Fourier transform of x' is performed to obtain its frequency spectrum f(x'), and the frequency spectrum length is p2;

[0075] S8, f(x') is processed according to S2-S5 to obtain the intrinsic mode function and singular value corresponding to x', wherein the intrinsic mode function is N2, and the intrinsic mode functions are ~ in turn, and the singular values are k3, and the singular values are σ1-σ k3 ;

[0076] S9, based on x', the intrinsic mode function and the singular value of x', a plurality of time domain features, frequency domain features and entropy domain features are extracted, and all the features constitute a multi-domain feature fusion vector F.

[0077] The entropy domain features include power spectrum entropy H eE , singular spectrum entropy H svdpE , and energy entropy H psdE , which are obtained according to the following formulas respectively.

[0078] ,

[0079] ,

[0080] ,

[0081] ,

[0082] wherein y m is the value of f(x') when the frequency is m, 1≤m≤P2, σ z is the zth singular value of x', k3 is the total number of singular values corresponding to x', 1≤z≤k3, Let x' be the q-th intrinsic mode function. At time t, E q for The energy value.

[0083] In S9, the time-domain features include the maximum value, minimum value, mean, variance, standard deviation, kurtosis, skewness, root mean square value, waveform factor, peak factor, impulse factor, and pulse gap factor of the reconstructed signal x'.

[0084] The frequency domain features include the center frequency, mean square frequency, root mean square frequency, frequency variance, and frequency standard deviation of the spectrum f(x');

[0085] The entropy domain features include power spectral entropy, singular spectral entropy, and energy entropy.

[0086] Example 3: See Figures 1-7 In this embodiment, the ASK-101 drilling rig's indoor drilling equipment was used to collect 700 vibration signal data points from shale, marl, limestone, coarse yellow sandstone, feldspathic sandstone, marble, and granite, with 100 vibration signal data points from each of the seven rock types. Following the method in Embodiment 1, a corresponding reconstructed signal was generated for each vibration signal according to steps S1-S6, and then a corresponding multi-domain feature fusion vector was generated according to steps S7-S9, resulting in a total of 700 multi-domain feature fusion vectors. In step S9, when extracting multiple time-domain, frequency-domain, and entropy-domain features, this embodiment extracted 12 time-domain features, labeled as f1~f... 12 The five frequency domain features are labeled as f 13 ~f 17 The three entropy domain features are labeled as f. 18 ~f 20 Then a multi-domain feature fusion vector F can be expressed as F=[f1,⋅⋅⋅,f 12 ,⋅⋅⋅f 20 ].

[0087] To demonstrate the advantages of the reconstructed signal of this invention compared to the original vibration signal, and the advantages of the multi-domain feature fusion vector compared to other feature vectors, the following two experiments are conducted:

[0088] I. Ablation Experiment:

[0089] Let F = [f1, ⋅⋅⋅, f 12 ,⋅⋅⋅f 20 F is the time-frequency entropy fusion feature. Three entropy domain features are removed from F to obtain F. 12 =[f1,⋅⋅⋅,f 12 ,⋅⋅⋅f 17] as a time-frequency fusion feature. In the ablation experiment, the time-frequency entropy fusion feature outperforms the time-frequency fusion feature in the lithology classification task. In terms of accuracy and recall rate, the model based on the time-frequency entropy fusion feature exhibits higher classification accuracy, especially when facing complex lithology samples. The time-frequency entropy fusion feature effectively captures key information in the signal, making the model more capable of distinguishing between different rock classes.

[0090] II. Comparative experiment

[0091] For the 700 vibration signals mentioned in the foregoing of the embodiment, the method S1-S6 of the application is used to reconstruct each vibration signal to generate 700 reconstructed signals as a self-built data set. The method of steps S7-S9 of the application is used as the feature extraction method of the application. Three existing feature extraction methods, EWTKC-SVD, IEWT-SVD, and IEWT-ISVD, are selected for feature extraction.

[0092] EWT-KC-SVD: Empirical Wavelet Transform-Hankel Matrix Construction-Singular Value Decomposition. Wu Y, Cheng T, Yi Q, et al. Research on Mine Microseismic Signal Recognition Based on Empirical Wavelet Transform[J]. Nonferrous Metal Science and Engineering, 2023, 14(3): 45-53.

[0093] IEWT-SVD: Improved Empirical Wavelet Transform-Singular Value Decomposition. Leng J, Diao K, Pang Z, Feng H. Automatic Identification Method of Offshore Platform Structural Modal Parameters Based on Improved Empirical Wavelet Transform[J]. Vibration and Shock, 2024, 43(7): 196-204.

[0094] IEWT-ISVD: Improved Empirical Wavelet Transform-Improved Singular Value Decomposition. See patent application No. 202310554003.9.

[0095] Then, the same classifier is used to classify rocks based on the extracted features, i.e., the 7 types of rocks mentioned in the foregoing. The performance indicators of the four methods, including accuracy, precision, recall, and F1 score, are analyzed to obtain the following Table 1:

[0096] Table 1 Comparison of different methods on self-built data set

[0097] Feature extraction method Accuracy rate Precision rate Recall rate F1 EWTKC-SVD 0.867 0.874 0.882 0.874 IEWT-SVD 0.889 0.892 0.894 0.872 IEWT-ISVD 0.896 0.897 0.904 0.902 The present invention 0.902 0.908 0.913 0.903

[0098] From Table 1, it can be seen that the application has obvious advantages in adaptability of spectral division and accuracy of feature extraction.

[0099] To further verify the generalization ability of the application, the public motor fault vibration data set is used for feature extraction and classification experiment, and the following Table 2 is obtained:

[0100] Table 2 Comparison of different methods on public datasets

[0101] Feature extraction method Accuracy rate Precision rate Recall rate F1 EWTKC-SVD 0.858 0.853 0.855 0.865 IEWT-SVD 0.861 0.859 0.860 0.874 IEWT-ISVD 0.872 0.861 0.863 0.882 The present invention 0.886 0.867 0.874 0.889

[0102] As shown in Table 2, through the processing of motor vibration signals, the present invention continues to demonstrate significant advantages in key indicators such as accuracy, precision, recall, and F1 score.

[0103] This invention selects an original vibration signal and uses empirical mode decomposition, empirical wavelet transform, and empirical wavelet transform autocorrelation analysis to decompose the signal and obtain... Figures 4-6 Then, the intrinsic mode functions in the reconstructed initial matrix A' are obtained by performing S1~S6 on the original vibration signal using the method of this invention, as shown below. Figure 7 As shown. The method for generating intrinsic mode functions according to steps S1-S6 of this invention is also called adaptive empirical wavelet transform. In this method, the vibration signal is decomposed by empirical mode decomposition to obtain 10 intrinsic mode functions. Figure 4 The functions, labeled IMF1 to IMF10, are decomposed by empirical wavelet transform to obtain six intrinsic mode functions. Figure 5 The signals, labeled IMF1~IMF6, were decomposed using empirical wavelet transform autocorrelation analysis to obtain six intrinsic mode functions. Figure 6 The intrinsic mode functions (IMFs) are labeled IMF1 to IMF6 and generated by this invention are 6 intrinsic mode functions. Figure 7 They are labeled as IMF1 to IMF6 respectively.

[0104] Combination Figures 4-7 Analysis of Tables 1 and 2 shows that the time-domain, frequency-domain, and entropy-domain feature values ​​extracted using the method of this invention have a certain degree of uniqueness, and some feature indicators fluctuate within a specific range, which can reflect the physical properties and structural characteristics of rocks to a certain extent. It also further demonstrates that there is a nonlinear relationship between the multi-domain statistical characteristics of vibration signals and the corresponding rock masses. The multi-domain feature vectors extracted from different rocks have different forms of expression, which can better reflect the characteristic differences of different rocks. This further verifies the effectiveness of the drilling vibration signal feature extraction and multi-domain fusion feature vector construction method of this invention.

[0105] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for optimizing and reconstructing drilling signals for lithology identification, wherein the drilling signal is the vibration signal of a rock sample drilled by the drill bit in the drilling area, and the rock sample is shale, marl, limestone, coarse yellow sandstone, feldspathic sandstone, marble, or granite, characterized in that, Includes the following steps: S1, obtain the vibration signal x, and obtain its spectrum f(x) by Fourier transform of x, with a spectrum length of p1; S2, adaptively generates a set of frequency band boundaries on f(x); S3, based on the empirical wavelet transform method to process the frequency band boundary, obtain N1 intrinsic mode functions imf1~imf N1 ; S4, Construct the initial matrix A, and construct the Hankel matrix H based on A. n×k1 ; A=(imf1,imf2,…imf N-1 ,imf N ) T T represents the transpose operation; ； In the formula, n and k1 are respectively H n×k1 The number of rows and columns, 1≤n≤k1, k1=N1-n+1; S5, for H n×k1 Perform singular value decomposition and select the first k2 singular values ​​to reconstruct the Hankel matrix H'; S6, Construct H based on the initial matrix A n×k1 The method is to perform an inverse transformation to obtain the initial reconstruction matrix A' from H', and then use each row element of A' as an intrinsic mode function to reconstruct the signal, thus obtaining the reconstructed signal x' of x; S7, obtain the spectrum f(x') of x' by Fourier transform, with a spectrum length of p2; S8, process f(x') according to S2~S5 to obtain the intrinsic mode functions and singular values ​​corresponding to x', where there are N2 intrinsic mode functions, in order as follows: ~ There are k3 singular values, numbered σ1 to σ2 respectively. k3 ; S9, based on the intrinsic mode functions and singular values ​​of x' and x', extracts multiple time-domain features, frequency-domain features and entropy-domain features, and constructs a multi-domain feature fusion vector F from all features; The entropy domain features include the power spectral entropy H. eE Singular spectral entropy H svdpE Energy entropy H psdE The results are obtained from the following formulas; ， ， ， ， In the formula, y m Let f(x') be the value of frequency m, 1≤m≤P2, σ z Let x' be the z-th singular value, and k3 be the total number of singular values ​​corresponding to x', where 1 ≤ z ≤ k3. Let x' be the q-th intrinsic mode function. At time t, E q for The energy value.

2. The drilling signal optimization and reconstruction method for lithology identification according to claim 1, characterized in that, S2 specifically includes Sa1~Sa3; Sa1, with a preset step size SL, where SL is the Nyquist frequency, and f(x) is divided into p1 segments L1~L according to the following formula. p1 The j-th segment is L j The maximum amplitude is A j , ; Sa2, preset amplitude threshold ∂, select A j The segments of >∂ are arranged in descending order of their maximum amplitude values ​​to form sequence L1; Sa3, for L1, generates a frequency band boundary by grouping every three adjacent segments, where one frequency band boundary [ω k ,ω h The generation method is as follows: obtain the frequencies corresponding to the maximum amplitude values ​​of a set of segments, and label them sequentially as ω. a ω b ω c ω k =(ω a +ω b ) / 2, ω h =(ω b +ω c ) / 2.

3. The drilling signal optimization and reconstruction method for lithology identification according to claim 1, characterized in that, S5 specifically includes S51~S53: S51, according to the following formula for H n×k1 Perform singular value decomposition; ； In the formula, U n×n Let S be an n×n left singular matrix, and V be a diagonal matrix. k1×k1 S is a k1×k1 right singular matrix containing n singular values; S52, by retaining the first k2 singular values ​​greater than 0 and their corresponding left and right singular matrices through the cumulative contribution rate, the reconstructed Hankel matrix H' is obtained.

4. The drilling signal optimization and reconstruction method for lithology identification according to claim 1, characterized in that, In S9, the time-domain features include the maximum value, minimum value, mean, variance, standard deviation, kurtosis, skewness, root mean square value, waveform factor, peak factor, impulse factor, and pulse gap factor of the reconstructed signal x'. The frequency domain features include the center frequency, mean square frequency, root mean square frequency, frequency variance, and frequency standard deviation of the spectrum f(x'); The entropy domain features include power spectral entropy, singular spectral entropy, and energy entropy.

Images