A method for extracting features of reference signals of a drill string while drilling

Abstract

The application relates to the technical field of geophysical exploration, and particularly discloses a method for extracting a reference signal feature of a drill string in seismic drilling, which comprises the following steps: obtaining a drill string reference signal; performing autocorrelation processing on the drill string reference signal to obtain a correlation domain drill string reference signal; obtaining a spectrum of the correlation domain drill string reference signal; decomposing the correlation domain drill string reference signal into two modes of a middle-low frequency mode and a high frequency mode according to the spectrum of the correlation domain drill string reference signal, to obtain a middle-low frequency signal and a high frequency signal; performing sparse decomposition on the middle-low frequency signal, and further introducing a regularization parameter optimization strategy based on a spectrum energy physical constraint to extract a low frequency signal; and subtracting the low frequency signal from the middle-low frequency signal to obtain a middle frequency signal. Random noise is suppressed by autocorrelation, and then frequency band separation and sparse reconstruction are performed, so that the target frequency band feature of the drill string reference signal is effectively extracted.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of geophysical exploration, in particular to a method for extracting features of a drill string reference signal in seismic while drilling. BACKGROUND

[0002] In oil and gas field exploration and development, accurately drilling to the target layer is a complex project with high cost and high risk. At present, drilling design mainly relies on surface seismic data and adjacent well information. The uncertainty of the stratum leads to frequent accidents during drilling operations. Seismic while drilling technology uses the vibration of the drill bit at the bottom of the well as a seismic source. The core is to collect the drill bit vibration signal (i.e. "drill string reference signal") through the top of the drill string or downhole sensors, and to extract the stratum impulse response through cross-correlation processing with the data recorded by the surface geophone. The essence is the relative travel time information of the drill bit vibration signal and the surface record. This method can continuously record and control drilling risks in real time without interrupting drilling operations and without the need for additional surface seismic sources.

[0003] However, the drill string reference signal collected by the pilot sensor installed at the top of the drill string is subject to multiple interferences in complex drilling environments: energy attenuation of the drill string transmission, virtual reflection notch effect, multiple wave interference, and fixed frequency noise generated by mud pumps and other equipment. Therefore, the drill string reference signal in seismic while drilling has low signal-to-noise ratio and serious time-frequency aliasing, and the effective components are difficult to identify. It is still a core problem that restricts the accuracy of seismic while drilling to extract the drill bit vibration signal from the mixed noise and transmission effect signal with high fidelity.

[0004] Correlation processing is a standardized signal extraction process. The typical processing steps for the correlation domain signal include band-pass filtering, median filtering in the FX domain, noise attenuation in the FK domain, and deconvolution. However, the traditional band-pass filter is prone to boundary effects and waveform distortion, the processing in the FX / FK domain is highly sensitive to parameter selection, and the adaptability to low signal-to-noise ratio data is limited. Therefore, it is urgent to develop alternative methods based on stronger robustness and higher fidelity to improve the feature extraction effect of the drill string reference signal in seismic while drilling. SUMMARY

[0005] In view of the problems existing in the prior art, the present application provides a method for extracting features of a drill string reference signal in seismic while drilling to solve the problems of difficult feature extraction of the drill string reference signal in seismic while drilling and difficulty in obtaining high-fidelity drill bit vibration signals in the prior art.

[0006] The present application provides a method for extracting features of a drill string reference signal in seismic while drilling, comprising:

[0007] Step S101: collecting well section data using a pilot sensor installed at the top of the drill string to obtain a drill string reference signal;

[0008] Step S102: autocorrelation processing is performed on the drill string reference signal to obtain a correlation domain drill string reference signal;

[0009] Step S103: frequency scanning analysis is performed on the correlation domain drill string reference signal to obtain a correlation domain drill string reference signal spectrum, and according to the correlation domain drill string reference signal spectrum, it is determined that the correlation domain drill string reference signal has different characteristics in different frequency bands;

[0010] Step S104: according to the correlation domain drill string reference signal spectrum, the correlation domain drill string reference signal is decomposed into two modes of a low-frequency mode and a high-frequency mode to obtain a low-frequency signal and a high-frequency signal;

[0011] Step S105: a Rake wavelet over-complete atom dictionary is constructed, the low-frequency signal is extracted by sparse decomposition of the low-frequency signal, and a regularization parameter optimization strategy based on spectral energy physical constraint is further introduced;

[0012] Step S106: subtract the low-frequency signal from the low-frequency signal to obtain a medium-frequency signal.

[0013] Further, the step S102 comprises:

[0014] Assuming that a drill bit vibration signal generated by drill bit vibration is , a drill string transmission effect is , a drill string reference signal is , a drill bit vibration signal , a Z transform of the drill bit vibration signal is , a Z transform of the drill string transmission effect is , and a Z transform of the drill string reference signal is , , , t , where t represents a time variable, and Z represents a Z variable, the Z transform of the drill string reference signal autocorrelation function is represented as:

[0015]

[0016] , where represents the Z transform of the drill string reference signal, represents the Z transform of the drill bit vibration signal, represents the Z transform of the drill string transmission effect, represents the Z transform of the time-reversed signal of the drill string reference signal, represents the Z transform of the time-reversed signal of the drill bit vibration signal, represents the Z transform of the time-reversed signal of the drill string transmission effect.

[0017] Assuming that a drill bit vibration signal If it's white noise, then... In the frequency domain, it is a constant, let it be... Then the Z-transform of the autocorrelation function of the drill string reference signal Represented as:

[0018]

[0019] Z-transform of the autocorrelation function of the drill string reference signal Inverse transform to the time domain to obtain the correlation domain drill string reference signal after autocorrelation processing. .

[0020] Further, step S104 includes:

[0021] Reference signal of the relevant domain drill string Decomposed into the first modal component and the second modal component The first modal component The second modal component is a mid-to-low frequency mode. The high-frequency mode is a mid-to-low frequency signal; the mid-to-low frequency mode is a high-frequency signal; the variational optimization objective is to minimize the sum of the frequency bandwidths of each modal component, with the constraint that the superposition of the modal components equals the original signal, expressed as:

[0022]

[0023]

[0024] in, Indicates the first One modal component, , For the first Each modal component corresponds to a center frequency; For time differential operators, For the Dirac function, This represents the convolution operation. The imaginary unit, The square of the L2 norm, t It is a time variable;

[0025] Introducing a penalty factor VMD-constrained Lagrange multipliers The constrained optimization problem is transformed into an unconstrained augmented Lagrangian function:

[0026]

[0027] in, Indicates the first One modal component, , is the center frequency of the kth modal component; is the center frequency of the kth modal component; is the time differential operator, is the Dirac function, denotes the convolution operation, is the imaginary unit, is the L2 norm square, t is the time variable, denotes the inner product operation, denotes the correlation domain drill string reference signal;

[0028] is the Fourier transform of , obtaining the kth modal component in the frequency domain is the Fourier transform of the correlation domain drill string reference signal , obtaining the frequency domain representation of the correlation domain drill string reference signal , then the kth modal component in the frequency domain at the nth iteration is:

[0029]

[0030] wherein, is the VMD constrained Lagrange multiplier in the frequency domain at the nth iteration, is the analytical solution after quadratic regularization of the bandwidth measure in the frequency domain; is the center frequency of the kth modal component at the nth iteration; is the frequency variable. Further, the step S104 further comprises: is obtained by using

[0031] to obtain discrete frequency points, the kth modal component in the frequency domain is updated at each discrete frequency point according to the following formula:

[0032]

[0033]

[0034] wherein, , , is the total number of discrete frequency points, is the sampling frequency, is the kth modal component in the frequency domain at the nth iteration, is the kth modal component in the frequency domain at the nth iteration, ​​​​​​​​is the frequency domain representation of the reference signal of the drill string in the correlation domain at the kth frequency point in the frequency domain;

[0035] the kth modal component at the nth iteration the center frequency of the kth modal component at the nth iteration The update formula at each discrete frequency point is:

[0036]

[0037] wherein, , , is the total number of discrete frequency points, is the sampling frequency, , denotes a set of discrete frequency points used for summation; is the power spectral density of the kth modal component at the nth iteration in the frequency domain Let

[0038] be the fidelity coefficient, the VMD constraint Lagrange multiplier of the kth iteration in the frequency domain The update formula of the VMD constraint Lagrange multiplier of the kth iteration in the frequency domain is:

[0039]

[0040] wherein, , , is the total number of discrete frequency points, is the sampling frequency, is the VMD constraint Lagrange multiplier of the kth iteration in the frequency domain, is the frequency domain representation of the reference signal of the drill string in the correlation domain at the kth frequency point in the frequency domain, is the kth modal component at the nth iteration in the frequency domain

[0041] When the following formula is satisfied, stop iteration, and obtain the kth modal component in the frequency domain ,

[0042]

[0043] wherein, , , is the total number of discrete frequency points, is the sampling frequency, , denotes a set of discrete frequency points used for summation,​​​​​​​ is the kth modal component in the frequency domain of the kth frequency point in the frequency domain, is the kth modal component in the frequency domain of the kth frequency point in the frequency domain, is the kth modal component in the frequency domain of the kth frequency point in the frequency domain, is the kth modal component in the frequency domain of the kth frequency point in the frequency domain, is the kth modal component in the frequency domain of the kth frequency point in the frequency domain, is the kth modal component in the frequency domain of the kth frequency point in the frequency domain, is a convergence threshold, is a frequency domain representation of a reference signal of the drill string in the relevant domain of the kth frequency point in the frequency domain;

[0044] performing Fourier inverse transform on the decomposed kth modal component in the frequency domain, to obtain the kth modal component in the time domain,

[0045]

[0046] wherein, is an imaginary unit, is a frequency variable, is a time variable, .

[0047] Further, the reference range of the penalty factor is , the fidelity coefficient is , and the convergence threshold is .

[0048] Further, the step S105 includes:

[0049] Step S1051, constructing a Ricker wavelet overcomplete atom dictionary H;

[0050] Step S1052, solving a constrained sparse reconstruction problem on the constructed Ricker wavelet overcomplete atom dictionary H:

[0051]

[0052] wherein, is an allowed reconstruction error, is a sparse coefficient, t is a time variable, is a medium-low frequency modal;

[0053] The above formula is changed into a second-order optimization problem without constraint:

[0054]

[0055] wherein, is a SPGL1 constraint Lagrange multiplier, H is a Ricker wavelet overcomplete atom dictionary, is a sparse coefficient,​​​ is a low frequency signal;

[0056] Step S1053, introducing a regularization parameter optimization strategy based on spectral energy physical constraints, to obtain the optimal regularization parameter ;

[0057] Step S1054, using the obtained optimal regularization parameter to solve the unconstrained second-order optimization problem:

[0058]

[0059] obtaining the optimal sparse coefficient , according to the optimal sparse coefficient , to obtain the reconstructed low frequency signal is:

[0060]

[0061] wherein, is a SPGL1 constraint Lagrange multiplier, H is a Ricker wavelet overcomplete atom dictionary, p is a sparse coefficient, is a low frequency signal.

[0062] Further, the step S1051 comprises:

[0063] The expression of the Ricker wavelet is:

[0064]

[0065] wherein, is the amplitude of the Ricker wavelet; is the center frequency; is the time delay of the Ricker wavelet, t is the time variable;

[0066] The Ricker wavelet overcomplete atom dictionary H is formed by adjusting and .

[0067] Further, the step S1053 comprises:

[0068] Let the target frequency band be , and define:

[0069] The low frequency signal in the frequency domain The energy in the target frequency band is:

[0070]

[0071] The reconstructed low frequency signal Energy in the same frequency band is:

[0072]

[0073] Target frequency band energy preservation rate is:

[0074]

[0075] Let the energy of the reconstructed low-frequency signal in the non-target frequency band be , and define the non-target frequency band energy leakage rate is:

[0076]

[0077] Low-frequency mode in frequency domain The average of the power spectral density in the target frequency band is:

[0078]

[0079] wherein, Low-frequency mode in frequency domain The power spectral density in the target frequency band;

[0080] Reconstructed low-frequency signal The average of the power spectral density in the target frequency band is:

[0081]

[0082] wherein, Reconstructed low-frequency signal The power spectral density in the target frequency band;

[0083] Target frequency band spectral correlation is expressed as:

[0084]

[0085] A representative segment of the low-frequency mode in the frequency domain is intercepted, and each SPGL1 constraint Lagrange multiplier is calculated by scanning the SPGL1 constraint Lagrange multiplier on a wider scale The corresponding target frequency band energy preservation rate , non-target frequency band energy leakage rate And target frequency band spectral correlation ;

[0086] Construct the objective function :

[0087]

[0088] wherein, , and are weight parameters, is the SPGL1 constraint Lagrange multiplier corresponding to the target frequency band energy preservation rate, is the SPGL1 constraint Lagrange multiplier corresponding to the non-target frequency band energy leakage rate, is the SPGL1 constraint Lagrange multiplier corresponding to the target frequency band spectral correlation, , M is the scanning number of the SPGL1 constraint Lagrange multiplier ; the SPGL1 constraint Lagrange multiplier that makes the minimum value of is taken as the optimal regularization parameter.

[0089] Further, the weight parameters satisfy the following conditions: .

[0090] Further, a relatively wide scale for scanning the SPGL1 constraint Lagrange multiplier is: , , M is the scanning number of the SPGL1 constraint Lagrange multiplier .

[0091] Compared with the prior art, the present application has the following remarkable technical effects.

[0092] (1) The traditional band-pass filter mainly realizes target frequency band signal extraction through passive attenuation of the spectrum, which will produce spectral boundary effect, phase distortion and waveform distortion. The present application first performs autocorrelation enhancement on the original record to produce periodic coherent components of the tri-cone bit, and then adaptively separates the modes (minimizes the mode bandwidth and estimates the mode center frequency) in the variational framework, so that the wide frequency energy produced by the tri-cone bit is concentrated into two independent modes, and the low frequency mode and the high frequency mode of the signal are separated; subsequently, based on the sparse representation of the Ricker wavelet dictionary matched with the seismic wavelet, the low frequency mode is processed. This processing method can better preserve the time-frequency characteristics of the original seismic source and reduce the waveform distortion and spectral boundary effect caused by the traditional band-pass filter.

[0093] (2) In the sparse reconstruction of the low and medium frequency mode, the traditional parameter selection method depends on experience or noise level estimation, and the reconstruction target frequency band drill string reference signal is not suitable for the case of denoising processing. The application proposes a physical constraint parameter optimization method based on spectral energy, which calculates the energy retention rate of the target frequency band, the energy leakage rate of the non-target frequency band, and the spectral correlation of the target frequency band, and selects the optimal regularization parameter by using scale scanning and weighted objective function. This method can improve the energy retention and spectral consistency of the reconstructed low frequency signal.

[0094] (3) The application combines autocorrelation preprocessing, mode separation and sparse reconstruction in order: autocorrelation enhances the drill bit coherent pulse, mode separation weakens the time-frequency aliasing problem to energy separation between high and low frequency modes, extracts high frequency components of strong interference, and sparse reconstruction accurately reconstructs strong energy low frequency signal on the over-complete wavelet dictionary. This method can reliably identify and extract the low, medium and high frequency components of the relevant domain drill string reference signal, reduce energy leakage and distortion between frequency bands, and facilitate subsequent joint analysis of drilling parameters, formation information and drill string reference signal. BRIEF DESCRIPTION OF DRAWINGS

[0095] In order to more clearly illustrate the technical solutions of the embodiments of the application, the following will briefly introduce the drawings needed to be used in the embodiments. Obviously, the drawings described in the following are only some embodiments of the application, and other drawings can also be obtained according to these drawings without creative labor for those skilled in the art.

[0096] Figure 1 A flowchart of a drill string reference signal feature extraction method provided by an embodiment of the application is shown in the figure.

[0097] Figure 2 A relevant domain drill string reference signal provided by an embodiment of the application is shown in the figure.

[0098] Figure 3 A relevant domain drill string reference signal spectrum provided by an embodiment of the application is shown in the figure.

[0099] Figure 4 A low frequency component of the relevant domain drill string reference signal extracted by an embodiment of the application is shown in the figure.

[0100] Figure 5 A low frequency component spectrum of the relevant domain drill string reference signal extracted by an embodiment of the application is shown in the figure.

[0101] Figure 6 A low frequency component of the relevant domain drill string reference signal extracted by a 4th order Butterworth bandpass filter provided by an embodiment of the application is shown in the figure.

[0102] Figure 7Low frequency component spectrum of the relevant domain drill string reference signal extracted by the fourth-order Butterworth band-pass filter provided in the embodiment of the present application;

[0103] Figure 8 Medium frequency component of the relevant domain drill string reference signal extracted in the embodiment of the present application;

[0104] Figure 9 Medium frequency component spectrum of the relevant domain drill string reference signal extracted in the embodiment of the present application;

[0105] Figure 10 High frequency component of the relevant domain drill string reference signal extracted in the embodiment of the present application;

[0106] Figure 11 High frequency component spectrum of the relevant domain drill string reference signal extracted in the embodiment of the present application. DETAILED DESCRIPTION

[0107] In order to better understand the technical solutions of the present application, the embodiments of the present application are described in detail below in combination with the drawings. It should be clear that the described embodiments are only some of the embodiments of the present application, but not all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without making creative efforts fall within the scope of protection of the present application.

[0108] The terms used in the embodiments of the present application are only for the purpose of describing specific embodiments, and are not intended to limit the present application. The singular forms "a", "an" and "the" used in the embodiments of the present application and the appended claims are also intended to include the plural forms, unless the context clearly indicates otherwise.

[0109] It should be understood that the term "and / or" used herein is only to describe the association relationship of the associated objects, which means that there can be three relationships, for example, A and / or B, which can represent the three cases of A alone, A and B together, and B alone. In addition, the character " / " in this paper generally represents that the front and rear associated objects are a "or" relationship.

[0110] Referring to Figure 1 A flowchart of a drill string reference signal feature extraction method for seismic while drilling provided in the embodiment of the present application is shown in FIG. 1. As shown in FIG. 1, it mainly includes the following steps. Figure 1

[0111] Step S101: Collecting well section data by using a pilot sensor installed at the top end of the drill string to obtain a drill string reference signal.

[0112] ​The data used in this embodiment comes from well G3 in an oilfield in eastern China. A pilot sensor installed at the top of the drill string is used to collect data on the well section with a drill bit depth of 2850m to 3290m as a reference signal for the drill string.

[0113] Step S102: Perform autocorrelation processing on the drill string reference signal to obtain the correlation domain drill string reference signal.

[0114] Assume the drill bit vibration signal generated by the drill bit vibration is The drill string transmission effect is The drill string reference signal is Drill bit vibration signal The 𝑍 transform is represented as Drill string transmission effect The 𝑍 transform is represented as Drill string reference signal The 𝑍 transform is represented as , t Let Z represent the time variable and Z represent the Z variable.

[0115] Drill string reference signal autocorrelation function Defined as:

[0116]

[0117] in, t Represents a time variable. This indicates the time delay of the signal's autocorrelation.

[0118] The above autocorrelation function Written in convolution form:

[0119]

[0120] in, Indicates drill string reference signal Time reversal signal, t Represents a time variable. This represents temporal convolution.

[0121] According to the convolution theorem of the Z-transform: the Z-transform of the convolution of two signals is equal to the product of their Z-transforms. Therefore, the Z-transform of the autocorrelation function mentioned above... for:

[0122]

[0123] in, Represents the Z-transform. Indicates drill string reference signal Time reversal signal, tdenotes a time variable, denotes a drill string reference signal of the Z-transform.

[0124] If the Z-transform of the signal is:

[0125]

[0126] wherein, denotes a Z-transform, denotes a drill string reference signal, t denotes a time variable, Z denotes Z a variable;

[0127] the Z-transform of the time reversed signal of the drill string reference signal is:

[0128]

[0129] wherein, denotes a Z-transform, t denotes a time variable, Z denotes a Z variable;

[0130] Let then i.e. the Z-transform of the time reversed signal of the drill string reference signal is:

[0131]

[0132] wherein, denotes a Z-transform, t denotes a time variable, Z denotes Z a variable;

[0133] The drill string reference signal is obtained from the bit vibration signal by a modification of the drill string transmission effect and is denoted as:

[0134]

[0135] The Z-transform of the drill string reference signal is denoted as:

[0136]

[0137] wherein, denotes a Z-transform, tRepresents a time variable. Z express Z variable; Indicates drill bit vibration signal The 𝑍 transformation, Indicates drill string transmission effect Z-transform.

[0138] Therefore, the Z-transform of the autocorrelation function of the drill string reference signal Represented as:

[0139]

[0140] in, This represents the α transform of the drill string reference signal. This represents the α transform of the drill bit vibration signal. The α transform represents the transmission effect of the drill string. The Z-transform of the time-reversed drill string reference signal. The Z-transform of the time-reversed drill bit vibration signal. The Z-transform of the time-reversed signal representing the transmission effect of the drill string.

[0141] Assuming drill bit vibration signal If it's white noise, then... In the frequency domain, it is a constant, let it be... Then the Z-transform of the autocorrelation function of the drill string reference signal Represented as:

[0142]

[0143] Z-transform of the autocorrelation function of the drill string reference signal Inverse transform to the time domain to obtain the correlation domain drill string reference signal after autocorrelation processing. .

[0144] After autocorrelation processing, both the drill bit vibration signal and the drill string transmission effect component are enhanced. When the drill bit vibration signal... When the noise level is white, autocorrelation preserves only the drill string transmission effect, allowing for the calculation of the period of drill string multiples and further estimation of the propagation speed of the drill bit vibration signal. Autocorrelation also highlights the overlapping coherent components in the drill string reference signal, suppressing random noise.

[0145] For the well section data collected from the drill bit depth of G3 well from 2850m to 3290m, one record was extracted every 10m interval. After autocorrelation processing, the first 3000ms record was taken as the reference signal for the correlation domain drill string in this embodiment. Figure 2 As shown.

[0146] Step S103: performing frequency scanning analysis on the correlation domain drill string reference signal to obtain a correlation domain drill string reference signal spectrum, and determining that the correlation domain drill string reference signal has different characteristics in different frequency bands according to the correlation domain drill string reference signal spectrum.

[0147] The frequency scanning analysis on the correlation domain drill string reference signal obtains a correlation domain drill string reference signal spectrum as shown in FIG. 3. Figure 3 As can be seen from FIG. 3, the correlation domain drill string reference signal has strong low-frequency and high-frequency component energy, and the signal has different characteristics in three frequency bands of a low frequency band of 1-12 Hz, a medium frequency band of 13-30 Hz and a high frequency band of 31-120 Hz.

[0148] Step S104: decomposing the correlation domain drill string reference signal into two modes of a low-medium frequency mode and a high frequency mode according to the correlation domain drill string reference signal spectrum to obtain a low-medium frequency signal and a high frequency signal.

[0149] The correlation domain drill string reference signal is decomposed into a first modal component and a second modal component , the first modal component is the low-medium frequency mode, and the second modal component is the high frequency mode, the low-medium frequency mode is the low-medium frequency signal, and the high frequency mode is the high frequency signal. The variational optimization objective is to minimize the sum of the frequency bandwidths of the modal components, and the constraint condition is that the superposition of the modal components is equal to the original signal, and the expression is as follows:

[0150]

[0151]

[0152] wherein, represents the first modal component, , , represents the second modal component, , is the center frequency corresponding to the first modal component; is a time differential operator, is a Dirac function, represents a convolution operation, is an imaginary unit, t is an L2 norm square, is a time variable.

[0153] A penalty factor and a VMD constraint Lagrange multiplier are introduced, and the constrained optimization problem is converted into an unconstrained augmented Lagrange function:

[0154]

[0155] in, Indicates the first One modal component, , For the first Each modal component corresponds to a center frequency; For time differential operators, For the Dirac function, This represents the convolution operation. The imaginary unit, The square of the L2 norm, t For time variables, This represents the inner product operation. Indicates the relevant domain drill string reference signal;

[0156] right Perform a Fourier transform to obtain the frequency domain... Modal components For the relevant domain drill string reference signal Perform a Fourier transform to obtain the frequency domain representation of the drill string reference signal in the correlation domain. Then the frequency domain The iteration of the ... Modal components for:

[0157]

[0158] in, For the frequency domain VMD-constrained Lagrange multipliers in the next iteration This is the analytical solution for the bandwidth metric after quadratic regularization in the frequency domain. For the first The modal component in the ... The center frequency of the next iteration; For frequency variables.

[0159] use To obtain discrete frequency points, the first frequency point in the frequency domain Modal components The update formula at each discrete frequency point is:

[0160]

[0161] in, , , This represents the total number of discrete frequency points. Sampling frequency, For the k-th frequency point in the frequency domain, the first... The iteration of the ... a modal component, a frequency domain representation of the correlation domain drill string reference signal at the kth frequency point in the frequency domain;

[0162] a modal component, a modal component at the kth frequency point in the frequency domain at the nth iteration a central frequency of the modal component at the nth iteration The update formula at each discrete frequency point is:

[0163]

[0164] wherein, , , the total number of discrete frequency points, a sampling frequency, , a set of discrete frequency points used for summation; a power spectral density of the modal component at the kth frequency point in the frequency domain at the nth iteration Let

[0165] be a fidelity coefficient, a VMD constraint Lagrange multiplier at the kth frequency point in the frequency domain at the nth iteration The update formula of the VMD constraint Lagrange multiplier at the kth frequency point in the frequency domain at the nth iteration is:

[0166]

[0167] wherein, , , the total number of discrete frequency points, a sampling frequency, a VMD constraint Lagrange multiplier at the kth frequency point in the frequency domain at the nth iteration, a frequency domain representation of the correlation domain drill string reference signal at the kth frequency point in the frequency domain, a modal component at the kth frequency point in the frequency domain at the nth iteration

[0168] When the following formula is satisfied, stop iteration, and obtain the kth modal component in the frequency domain ,

[0169]

[0170] wherein, , , the total number of discrete frequency points, a sampling frequency, ,​​​​​​ a set of discrete frequency points for summation, is the kth modal component of the 1th iteration of the frequency domain kth frequency point, is the kth modal component of the 1th iteration of the frequency domain kth frequency point, is the kth modal component of the 1th iteration of the frequency domain kth frequency point, is the kth modal component of the 1th iteration of the frequency domain kth frequency point, is the kth modal component of the 1th iteration of the frequency domain kth frequency point, is the kth modal component of the 1th iteration of the frequency domain kth frequency point, is a convergence threshold, is a frequency domain representation of the correlation domain drill string reference signal of the kth frequency point in the frequency domain;

[0171] performing inverse Fourier transform on the decomposed kth modal component of the frequency domain, to obtain the kth modal component of the solution,

[0172]

[0173] wherein, is an imaginary unit, is a frequency variable, is a time variable, .

[0174] Parameter selection in the above algorithm: the number of modes is selected as 2 according to the spectrum of the correlation domain drill string reference signal; the reference range of the penalty factor is , since the spectrum of the correlation domain drill string reference signal is wide, a lower penalty factor can relax the modal bandwidth constraint, thereby ensuring that the decomposition result covers a wide frequency band and avoiding energy leakage caused by excessive shrinkage, and therefore is selected in an embodiment; a smaller fidelity coefficient can ensure that the decomposition result accurately matches the real components of the signal, and therefore is selected in an embodiment; as shown in Figure 2 , the correlation domain drill string reference signal does not contain a direct current component, and therefore the DC component is set to 0 to remove the influence of direct current drift; the convergence threshold adopts a strict convergence condition to ensure the stability and accuracy of modal decomposition and avoid the problem of incomplete frequency bands caused by premature stopping.

[0175] Step S105: constructing a Ricker wavelet overcomplete atom dictionary, performing sparse decomposition on the medium-low frequency signal, and further introducing a regularization parameter optimization strategy based on spectral energy physical constraints to extract a low frequency signal, comprising:

[0176] Step S1051, constructing a Ricker wavelet overcomplete atom dictionary H, and the expression of the Ricker wavelet is:

[0177] ​​​

[0178] wherein, is the amplitude of the Ricker wavelet; is the center frequency; is the time delay of the Ricker wavelet, t is the time variable;

[0179] by regulating and to form a Ricker wavelet overcomplete atom dictionary H.

[0180] Since the Ricker wavelet frequency band range is uncertain, the full width at half maximum of the Ricker wavelet amplitude spectrum is used to screen the atoms on the target frequency band at an interval of 1 Hz, and the Ricker wavelet with energy mainly concentrated in the target frequency band is selected for constructing the atom dictionary. The specific method is as follows: first, mirror symmetry is performed on the first 500 ms of the decomposed low-frequency mode to prevent boundary effects, and then a Ricker wavelet dictionary is constructed using Ricker wavelets with different center frequencies and time delays. The finally constructed overcomplete atom dictionary has a center frequency range of [1, 6] with a step size of 0.5 Hz and a sampling rate of 1000 Hz, covering the entire signal in the time scale with a range of [0, 3499] and a step size of 2 ms. Therefore, the size of the constructed atom dictionary is 3500 x 19250.

[0181] Step S1052, solve the sparse reconstruction problem with constraints on the constructed Ricker wavelet overcomplete atom dictionary H:

[0182]

[0183] wherein, is the allowed reconstruction error, is the sparse coefficient, t is the time variable, is the low-frequency mode;

[0184] The above formula is changed to solve the unconstrained second-order optimization problem:

[0185]

[0186] wherein, is the SPGL1 constraint Lagrange multiplier, H is the Ricker wavelet overcomplete atom dictionary, is the sparse coefficient, is the low-frequency mode.

[0187] Step S1053, introduce a regularization parameter optimization strategy based on the physical constraint of spectral energy, and set the target frequency band as , and define:

[0188] low-frequency mode in frequency domain Energy in the target band is:

[0189]

[0190] Reconstructed low frequency signal Energy in the same band is:

[0191]

[0192] Target band energy preservation rate is:

[0193]

[0194] Let the energy of the reconstructed low frequency signal in the non-target band be , define the non-target band energy leakage rate is:

[0195]

[0196] Low frequency mode in the frequency domain Mean value of the power spectral density in the target band is:

[0197]

[0198] Wherein, Low frequency mode in the frequency domain Power spectral density in the target band

[0199] Reconstructed low frequency signal Mean value of the power spectral density in the target band is:

[0200]

[0201] Wherein, Reconstructed low frequency signal Power spectral density in the target band

[0202] Target band spectral correlation is expressed as:

[0203]

[0204] A representative segment of the low frequency mode in the frequency domain is intercepted, and each SPGL1 constraint Lagrange multiplier is calculated by scanning the SPGL1 constraint Lagrange multiplier on a wider scale Corresponding target band energy preservation rate , non-target band energy leakage rate and target band spectrum correlation . The one wider scale is: , M is the SPGL1 constraint Lagrange multiplier the number of scans.

[0205] In one embodiment, the value range of is: , .

[0206] Constructing the objective function :

[0207]

[0208] wherein, , and are weight parameters, M is the SPGL1 constraint Lagrange multiplier corresponding to the target band energy preservation rate, M is the SPGL1 constraint Lagrange multiplier corresponding to the non-target band energy leakage rate, M is the SPGL1 constraint Lagrange multiplier corresponding to the target band spectrum correlation, M is the SPGL1 constraint Lagrange multiplier the number of scans. The SPGL1 constraint Lagrange multiplier that makes the minimum value of is taken as the optimal regularization parameter.

[0209] The target of the present application is to extract the low frequency signal of the target band, the target band energy preservation rate, the non-target band energy leakage rate and the contribution of the target band spectrum correlation to the target band, which is target band spectrum correlation > target band energy preservation rate > energy target band energy leakage rate. Therefore, the selection of the weight parameter should satisfy In one embodiment, the weight parameter is set as .

[0210] Step S1054, the obtained optimal regularization parameter is used to solve the unconstrained second-order optimization problem:

[0211]

[0212] The optimal sparse coefficient is obtained, according to the optimal sparse coefficient , the reconstructed low frequency signal is obtained as:

[0213]

[0214] wherein, is the SPGL1 constraint Lagrange multiplier, H is the overcomplete dictionary of Ricker wavelets, p is the sparse coefficient, is the low frequency component of the relevant domain drill string reference signal extracted by the embodiment of the present application,

[0215] Step S106: subtracting the low frequency signal from the medium-low frequency signal to obtain a medium frequency signal.

[0216] Figure 4 is the low frequency component of the relevant domain drill string reference signal extracted by the embodiment of the present application, Figure 5 is the low frequency component spectrum of the relevant domain drill string reference signal extracted by the embodiment of the present application; wherein Figure 4 and Figure 5 As can be seen, the low frequency component spectrum extracted by the embodiment of the present application is concentrated in 1~12Hz, the time domain waveform is relatively smooth, and is in good agreement with the target frequency band of the original relevant domain drill string reference signal (as shown in Figure 2 ); Figure 6 is the low frequency component of the relevant domain drill string reference signal extracted by the 4th order Butterworth band-pass filter provided by the embodiment of the present application, Figure 7 is the low frequency component spectrum of the relevant domain drill string reference signal extracted by the 4th order Butterworth band-pass filter provided by the embodiment of the present application, and by comparison Figure 4-7 , it can be seen that the method of the present application effectively suppresses the spectrum boundary effect and waveform distortion, and better preserves the time-frequency characteristics of the drill bit source. After removing the low frequency component, the medium frequency component of 12~30Hz of the relevant domain drill string reference signal is successfully extracted, Figure 8 is the medium frequency component of the relevant domain drill string reference signal extracted by the embodiment of the present application, Figure 9 is the medium frequency component spectrum of the relevant domain drill string reference signal extracted by the embodiment of the present application; wherein Figure 8 As can be seen, each trace record periodically appears a peak value, which shows that the signal has a periodicity of 200~210ms. This is related to the mechanical rotating speed of the three-cone drill bit, and the rotating speed of the drill bit is basically maintained at 98~101rpm during the acquisition process. The multiple wave development in the medium frequency record is obvious, including short period drill tool assembly multiple wave and long period drill string multiple wave, which can be used for monitoring the drill tool assembly state; based on the characteristics of the time delay of the first order and second order multiple wave, the propagation speed of the energy of the drill bit vibration along the drill string calculated by the embodiment of the present application is 5084.7m / s, which is consistent with the range of the longitudinal wave speed of the steel drill string, and the extracted relevant domain medium frequency drill string reference signal can reflect the propagation characteristics of the drill string and estimate the propagation time delay of the energy of the drill bit vibration along the drill string. The high frequency component of 31~120Hz is directly separated by step S104, Figure 10High frequency component of the relevant domain drill string reference signal extracted for the embodiment of the present application, Figure 11 High frequency component spectrum of the relevant domain drill string reference signal extracted for the embodiment of the present application, Figure 10 And Figure 11 It can be seen that the time domain signal extreme value presents a periodic distribution related to the drill bit mechanical rotating speed, and appears at the position of 630 ms and its integer multiples, which can provide a basis for monitoring the drill bit mechanical rotating speed and bearing mechanical state.

[0217] Compared with the prior art, the present application has the following remarkable technical effects.

[0218] (1) The traditional band-pass filter mainly realizes target frequency band signal extraction through passive attenuation of the frequency spectrum, which will produce spectrum boundary effect, phase distortion and waveform distortion. The present application first performs autocorrelation enhancement on the original record to produce periodic coherent components of the three-cone drill bit, and then adaptively separates the modes (minimizes the mode bandwidth and estimates the mode center frequency) in the variational framework, so that the wide frequency energy produced by the three-cone drill bit is concentrated into two independent modes, and the low frequency mode and the high frequency mode of the signal are separated; then, on the low frequency mode, sparse representation is performed based on the Ricker wavelet dictionary matched with the seismic wavelet. This processing method can better preserve the time-frequency characteristics of the original seismic source and reduce the waveform distortion and spectrum boundary effect caused by the traditional band-pass filter.

[0219] (2) In the sparse reconstruction of the low frequency mode, the traditional parameter selection method mainly depends on experience or noise level estimation, and the method is no longer applicable to the case of reconstructing the target frequency band drill string reference signal instead of denoising. The present application proposes a physical constraint parameter optimization method based on spectrum energy, which quantitatively calculates the energy retention rate of the target frequency band, the energy leakage rate of the non-target frequency band, and the spectrum correlation of the target frequency band, and selects the optimal regularization parameter by using scale scanning and weighted objective function. This method can improve the energy retention and spectral consistency of the reconstructed low frequency signal.

[0220] (3) The present application combines autocorrelation preprocessing, mode separation and sparse reconstruction in an orderly manner: autocorrelation enhancement of drill bit coherent pulse, mode separation weakens the time-frequency aliasing problem to energy separation between high and low frequency modes, extracts the high frequency component of strong interference, and sparse reconstruction accurately reconstructs the strong energy low frequency signal in a sparse manner on the Ricker wavelet over-complete atomic dictionary. This method can reliably identify and extract the low frequency, medium frequency and high frequency components of the relevant domain drill string reference signal, reduce the energy leakage and distortion between frequency bands, and is beneficial to the subsequent joint analysis of drilling parameters, formation information and drill string reference signal.

[0221] In summary, the cooperative processing procedure adopted by the embodiments of the present application performs excellently in reducing inter-band energy leakage, reducing boundary effect and maintaining waveform physical consistency, and the processing result is ideal, which can provide reliable data basis for subsequent drilling parameter analysis, drilling tool state diagnosis and drilling formation information prediction.

[0222] In the embodiments of the present application, "at least one" means one or more, and "multiple" means two or more. The "and / or" describes the association relationship of the associated objects, which means that there can be three kinds of relationships, for example, A and / or B can represent the cases of A alone, A and B together, and B alone. Wherein A and B can be singular or plural. The character " / " generally represents that the front and rear associated objects are in an "or" relationship. "At least one of the following" and similar expressions mean any combination of these items, including any combination of single or multiple items. For example, at least one of a, b and c can represent: a, b, c, a-b, a-c, b-c or a-b-c, wherein a, b and c can be single or multiple.

[0223] The above description is merely specific embodiments of the present application, and any skilled in the art within the technical scope disclosed by the present application can easily think of changes or replacements, which shall be covered within the protection scope of the present application.

Claims

1. A method for extracting features of a reference signal of a drillstring while drilling, characterized in that, The method comprises the steps of: Step S101: collecting well section data by using a pilot sensor installed at the top of a drill string to obtain a drill string reference signal; Step S102: performing autocorrelation processing on the drill string reference signal to obtain a correlation domain drill string reference signal; The step S102 comprises: Assume the drill bit vibration signal generated by the drill bit vibration is The drill string transmission effect is The drill string reference signal is Drill bit vibration signal The 𝑍 transform is represented as Drill string transmission effect The 𝑍 transform is represented as Drill string reference signal The 𝑍 transform is represented as , t Let Z represent the time variable and Z represent the Z-variable. Then, the Z-transform of the autocorrelation function of the drill string reference signal is... Represented as: wherein, Z-transform of a drill string reference signal, Z-transform of a drill bit vibration signal, Z-transform of a drill string transmission effect, Z-transform of a time-reversed signal of a drill string reference signal, Z-transform of a time-reversed signal of a drill bit vibration signal, Z-transform of a time-reversed signal of a drill string transmission effect; Assume that the drill bit vibration signal is white noise, then is constant in the frequency domain, set as , then the Z transform of the drill string reference signal autocorrelation function is expressed as: Z-transform of the drill string reference signal autocorrelation function inverse transform to the time domain to obtain the autocorrelation processed correlation domain drill string reference signal ; Step S103: performing frequency scanning analysis on the correlation domain drill string reference signal to obtain a correlation domain drill string reference signal spectrum, and determining that the correlation domain drill string reference signal has different characteristics in different frequency bands according to the correlation domain drill string reference signal spectrum; Step S104: decomposing the correlation domain drill string reference signal into two modes of a medium-low frequency mode and a high frequency mode according to the correlation domain drill string reference signal spectrum to obtain a medium-low frequency signal and a high frequency signal; The step S104 comprises: Correlating domain drill string reference signals Decompose into 1st modal component And 2nd modal component 1st modal component Is mid-low frequency modal, 2nd modal component Is high frequency modal, mid-low frequency modal is mid-low frequency signal; high frequency modal is high frequency signal; variational optimization objective is to minimize sum of frequency bandwidth of each modal component, constraint condition is that modal component superposition equals to original signal, expression is: wherein, represents the th modal component, , is the th modal component corresponding center frequency; is the time differential operator, is the Dirac function, represents the convolution operation, is the imaginary unit, is the L2 norm square, t is the time variable; Introducing a penalty factor and the VMD constraint Lagrange multiplier transforming the constrained optimization problem into an unconstrained augmented Lagrangian function: wherein, represents the th modal component, , is the th modal component corresponding center frequency; is the time differential operator, is the Dirac function, represents the convolution operation, is the imaginary unit, is the L2 norm square, t is the time variable, represents the inner product operation, represents the correlation domain drill string reference signal; ​​​​​​​​​ wherein, VMD is the VMD constraint Lagrange multiplier for the frequency domain, is the analytical solution of the bandwidth metric after quadratic regularization in the frequency domain; is the analytical solution of the bandwidth metric after quadratic regularization in the frequency domain; is the center frequency of the th modal component at the th iteration; is the frequency variable; Step S105: constructing a Ricker wavelet over-complete atom dictionary, performing sparse decomposition on the medium-low frequency signal, and further introducing a regularization parameter optimization strategy based on a spectrum energy physical constraint to extract a low frequency signal; Step S106: subtracting the low frequency signal from the medium-low frequency signal to obtain a medium frequency signal.

2. The method of claim 1, wherein, The step S104 further comprises: Adopting obtained discrete frequency points, the first modal component in the frequency domain The update formula at each discrete frequency point is: wherein, , , is the total number of discrete frequency points, is the sampling frequency, is the kth modal component of the kth frequency point in the frequency domain at the nth iteration, is the kth modal component of the kth frequency point in the frequency domain at the nth iteration, is the kth modal component of the kth frequency point in the frequency domain at the nth iteration, is the frequency domain representation of the associated domain tubular reference signal at the kth frequency point in the frequency domain; The first modal component is at the center frequency of the first iteration The update formula at each discrete frequency point is: in, , , This represents the total number of discrete frequency points. Sampling frequency, , This represents the set of discrete frequency points used for summation; For the k-th frequency point in the frequency domain, the first... Power spectral density of modal components in the next iteration; set up For fidelity coefficients, the k-th frequency point in the frequency domain is the first... VMD-constrained Lagrange multipliers in the next iteration The update formula is: in, , , This represents the total number of discrete frequency points. Sampling frequency, For the k-th frequency point in the frequency domain, the first... VMD-constrained Lagrange multipliers in the next iteration Let be the frequency domain representation of the correlation domain drill string reference signal at the k-th frequency point in the frequency domain. For the k-th frequency point in the frequency domain, the first... The iteration of the ... One modal component; When the following formula is satisfied, the iteration is stopped, and the frequency domain first modal component is obtained in, , , This represents the total number of discrete frequency points. Sampling frequency, , This represents the set of discrete frequency points used for summation. For the k-th frequency point in the frequency domain, the first... The iteration of the ... One modal component, For the k-th frequency point in the frequency domain, the first... The iteration of the ... One modal component, The convergence threshold, This is the frequency domain representation of the correlation domain drill string reference signal at the k-th frequency point in the frequency domain; The decomposed frequency domain first modal component The decomposed frequency domain first modal component The decomposed frequency domain first modal component The decomposed frequency domain first modal component The decomposed frequency domain first modal component wherein is the imaginary unit, is the frequency variable, is the time variable, .

3. The method of claim 2, wherein, a reference range for the penalty factor is a fidelity coefficient a convergence threshold .

4. The method of claim 2, wherein, The step S105 comprises: Step S1051: constructing a Ricker wavelet over-complete atom dictionary H; Step S1052: solving a sparse reconstruction problem with constraints on the constructed Ricker wavelet over-complete atom dictionary H: wherein, is the allowed reconstruction error, is the sparse coefficient, t is the time variable, is the mid-low frequency mode; Change the above formula into a second-order optimization problem without constraints: wherein, is the SPGL1 constraint Lagrange multiplier, H is the Rake sub-wave over-complete atom dictionary, is the sparse coefficient, is the middle-low frequency mode; Step S1053, a regularization parameter optimization strategy based on spectrum energy physical constraint is introduced to obtain the optimal regularization parameter ; Step S1054, using the obtained optimal regularization parameter solving the unconstrained second-order optimization problem: obtaining optimal sparse coefficients , reconstructing a low frequency signal , reconstructing a low frequency signal is wherein, is the SPGL1 constraint Lagrange multiplier, H is the Rake sub-wave over-complete atom dictionary, p is the sparse coefficient, is the mid-low frequency mode.

5. The method of claim 4, wherein, The step S1051 comprises: The expression of the Ricker wavelet is: wherein, is the amplitude of the Ricker wavelet; is the center frequency; is the Ricker wavelet time delay, t is the time variable; By regulating and Forming a rake wavelet overcomplete atom dictionary H.

6. The method of claim 4, wherein, The step S1053 comprises: Target frequency band is set as Definition: Low frequency mode in frequency domain Energy within target frequency band is: Reconstructing low frequency signals Energy within the same frequency band is: Then the target frequency band energy reservation rate is: Let the energy of the reconstructed low frequency signal in the non-target frequency band be , and define the non-target frequency band energy leakage rate as: Low frequency mode in frequency domain Mean value of the power spectral density in the target frequency band is: wherein is a low frequency mode in the frequency domain power spectral density in the target frequency band; Reconstructing low frequency signals Mean value of the power spectral density in the target frequency band is: wherein, to reconstruct the low frequency signal power spectral density in the target frequency band; The target band spectrum correlation is represented as: A representative segment of the low frequency mode in the frequency domain is intercepted, and the SPGL1 constraint Lagrange multipliers are scanned in a wide scale, and each SPGL1 constraint Lagrange multiplier is calculated Corresponding target frequency band energy retention rate , non-target frequency band energy leakage rate And target frequency band spectral correlation ; Constructing the objective function : wherein, , and are weight parameters, is the SPGL1 constraint Lagrange multiplier corresponding to the target band energy preservation rate, is the SPGL1 constraint Lagrange multiplier corresponding to the non-target band energy leakage rate, is the SPGL1 constraint Lagrange multiplier corresponding to the target band spectral correlation, , M is the SPGL1 constraint Lagrange multiplier the number of scans; the SPGL1 constraint Lagrange multiplier that minimizes is taken as the optimal regularization parameter.

7. The method of claim 6, wherein, The weight parameters satisfy the following conditions: .

8. The method of claim 6, wherein, A wider scale of scanning the SPGL1 constraint Lagrange multiplier is: , , M is the number of scanning the SPGL1 constraint Lagrange multiplier .

Images