Real-time drilling depth measuring system of geological exploration drilling machine

Abstract

The invention relates to the technical field of geological prospecting, in particular to a geological prospecting drilling rig drilling real-time depth measuring system, which comprises a vibration acquisition module, a data processing module and a data processing module, the original vibration signals are divided into 256 time sequence segments, sampling interval values are set, the amplitude value and the phase value of each time sequence segment are calculated, and drilling dynamic characteristics are generated. According to the method, vibration sensors are arranged on the surface of a drilling tool to collect original vibration signals, the signals are segmented to calculate amplitude values and phase values, drilling dynamic characteristics are obtained, difference values of adjacent time sequence segments are calculated and compressed, Fourier transform is carried out to generate frequency spectrum parameters, and the frequency spectrum parameters are projected to Grassmann space to calculate characteristic vectors and characteristic values. And constructing a tangent plane projection function to calculate a geodesic distance to obtain a correction coordinate, and iteratively calculating a depth value based on a Riccati equation. And meanwhile, time sequence segmented compression is adopted, so that data redundancy is reduced, and the utilization efficiency of computing resources is optimized.

Description

Technical Field

[0001] The present invention relates to the technical field of geological exploration, and particularly to a real-time depth measurement system for drilling of geological exploration drills. Background Art

[0002] Geological exploration is a comprehensive field of earth science and technology, mainly studying the geological structure, distribution of mineral resources, and characteristics of geological environment on the earth's surface and its deep part. This field involves multiple sub-disciplines such as geophysical exploration, geochemical exploration, and drilling engineering. Through various exploration means and technical methods, systematic exploration and evaluation of underground resources are carried out to provide important geological data support for mineral development, engineering construction, environmental protection, etc. However, in the process of traditional geological exploration drilling, there are problems such as low sampling frequency and poor accuracy in the data acquisition link, which are difficult to accurately reflect the real-time state of the drill string; at the same time, the signal processing ability is limited, the quality of the collected data is unstable, and there are errors. Therefore, improvement is needed. Summary of the Invention

[0003] The purpose of the present invention is to solve the disadvantages existing in the prior art, and to propose a real-time depth measurement system for drilling of geological exploration drills.

[0004] In order to achieve the above purpose, the present invention adopts the following technical scheme: The real-time depth measurement system for drilling of geological exploration drills includes:

[0005] A vibration acquisition module, arranging vibration sensors on the surface of the drill string to obtain the original vibration signal; dividing the original vibration signal into 256 time series segments and setting the sampling interval value, calculating the amplitude and phase values for each time series segment, and generating the dynamic characteristics of drilling.

[0006] A signal transformation module, calculating the amplitude difference and phase difference between adjacent time series segments in the dynamic characteristics of drilling, compressing the time series segments according to the amplitude difference and phase difference to generate a compressed vibration sequence, and performing Fourier transform on the compressed vibration sequence to generate characteristic spectrum parameters.

[0007] A state correction module, projecting the characteristic spectrum parameters into the Grassmann space, calculating the projection eigenvector and eigenvalue to generate projection matrix parameters; constructing a tangent plane projection function for the projection matrix parameters and calculating the geodesic distance to generate the corrected space coordinates.

[0008] A depth determination module, based on the corrected space coordinates, calculating the position offset and angle offset of the drill string in the three-dimensional space through the tangent plane projection function, substituting the position offset and angle offset into the Riccati equation for iterative calculation and outputting the depth value to generate real-time depth data.

[0009] Preferably, the acquisition step of the original vibration signal is:

[0010] Vibration sensors are arranged on the surface of the drill string to collect the spatial position change amount and frequency change amount during the drilling process, and the original vibration signal is obtained;

[0011] Based on the original vibration signal, signal amplification, filtering and digitization are performed to obtain the preprocessed original vibration signal;

[0012] Based on the preprocessed original vibration signal, the continuous vibration signal is converted into discrete data points, time series analysis is performed on the discrete data points, and the dynamic characteristics of the drilling process are extracted to obtain the original vibration signal.

[0013] Preferably, the steps for obtaining the drilling dynamic characteristics are as follows:

[0014] The original vibration signal is divided into 256 time series segments, and a unified sampling interval value is set for each segment to obtain the segmented time series segments;

[0015] Based on the segmented time series segments, the amplitude and phase are calculated for each time series segment to obtain the amplitude-phase characteristics of each time series segment;

[0016] Based on the amplitude-phase characteristics of each time series segment, by analyzing the characteristics of all time series segments, statistical methods and pattern recognition are used to identify the dynamic change patterns during the drilling process, and the drilling dynamic characteristics are generated.

[0017] Preferably, the steps for obtaining the compressed vibration sequence are as follows:

[0018] Adjacent time series segments are extracted from the drilling dynamic characteristics, the amplitude difference and phase difference are calculated for the sequence segments, and by comparing the characteristics of each sequence segment, the vibration changes during the drilling process are captured to obtain the amplitude difference and phase difference data;

[0019] Based on the amplitude difference and phase difference data, data compression is performed on the time series segments, and adjacent and similar-featured sequence segments are merged to obtain the compressed time series segments;

[0020] Based on the compressed time series segments, re-integration is performed to construct a continuous time series to obtain the compressed vibration sequence.

[0021] Preferably, the steps for obtaining the characteristic spectrum parameters are as follows:

[0022] Data is extracted from the compressed vibration sequence to obtain the data ready for Fourier transform;

[0023] Based on the data ready for Fourier transform, Fourier transform is performed to obtain the spectrum data, and the formula is:

[0024]

[0025] Wherein, F(k) is the k-th frequency component in the frequency domain, x(j) is the j-th data point in the time series, and N is the total number of data points;

[0026] Based on the spectral data, the sparsity of the spectrum is optimized by applying the L1 norm constraint to obtain the characteristic spectrum parameters.

[0027] Preferably, the steps for obtaining the projection matrix parameters are as follows:

[0028] Based on the characteristic spectrum parameters, a frequency component matrix and frequency characteristic data are extracted, the frequency component matrix is decentralized to generate a decentralized frequency component matrix;

[0029] Based on the decentralized frequency component matrix, a projection matrix is calculated, and the expression is:

[0030]

[0031] Wherein, M ij is the element in the i-th row and j-th column of the projection matrix, f i and f j are respectively the amplitudes of the i-th and j-th frequency components in the frequency component matrix, u i and u j are respectively the eigenvectors of the i-th and j-th frequency components, λ i and λ j are respectively the eigenvalues of the i-th and j-th frequency components, θ i and θ j are respectively the phase angles of the i-th and j-th frequency components, φ i and φ j are respectively the amplitude angles of the i-th and j-th frequency components;

[0032] Based on the projection matrix, it is mapped to the Grassmann space, and geometric operations are performed through eigenvalue decomposition to obtain eigenvectors and eigenvalues, forming the projection matrix parameters.

[0033] Preferably, the steps for obtaining the corrected space coordinates are as follows:

[0034] Based on the projection matrix parameters, a tangent plane projection function is constructed, and the formula is:

[0035]

[0036] Wherein, x and y respectively represent the coordinate values of any two points in the projection matrix, and F(x, y) is the tangent plane projection function;

[0037] The distance of all points in the projection matrix is calculated using the tangent plane projection function, and the distance of each pair of points is mapped to a new distance value;

[0038] Calculate the geodesic distance between data points through the new distance value, and use the geodesic distance to correct the spatial coordinates to obtain the corrected spatial coordinates.

[0039] Preferably, the step of obtaining the real-time depth data is as follows:

[0040] Based on the corrected spatial coordinates, use the tangent plane projection function to calculate the position offset and angle offset of the drilling tool in three-dimensional space;

[0041] Based on the position offset and angle offset, calculate the depth adjustment value, and the calculation formula is:

[0042]

[0043] where, Δx and Δy represent the position offsets of the drilling tool in the horizontal and vertical directions, α is the total angle offset of the drilling tool, and β is the angle with the vertical axis;

[0044] Based on the depth adjustment value, substitute the depth adjustment value into the Riccati equation, perform iterative calculations to adjust the depth value of the drilling tool, and generate real-time depth data.

[0045] Compared with the prior art, the advantages and positive effects of the present invention are as follows:

[0046] In the present invention, vibration sensors are arranged on the surface of the drilling tool to collect the original vibration signals, and the signals are segmented to calculate the amplitude and phase values to obtain the drilling dynamic characteristics. The differences are calculated and compressed for adjacent time series segments, and Fourier transform is performed to generate spectral parameters, which are projected onto the Grassmann space to calculate the eigenvectors and eigenvalues. A tangent plane projection function is constructed to calculate the geodesic distance to obtain the corrected coordinates, and the depth value is iteratively calculated based on the Riccati equation. At the same time, time series segmentation compression is adopted to reduce data redundancy and optimize the utilization efficiency of computing resources; the introduction of Grassmann space projection and tangent plane projection function to establish a high-dimensional feature mapping relationship improves the depth calculation accuracy; the iterative operation based on the Riccati equation ensures the accuracy of spatial positioning and provides a reliable basis for real-time correction of the drilling trajectory. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] Figure 1 is the system flow chart of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

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

[0049] Please refer to Figure 1, the present invention provides a technical solution: a real-time depth measurement system for geological exploration drilling rigs includes:

[0050] A vibration acquisition module arranges vibration sensors on the surface of the drill string to obtain the original vibration signal; divides the original vibration signal into 256 time series segments and sets the sampling interval value, calculates the amplitude and phase values for each time series segment, and generates the dynamic characteristics of drilling.

[0051] A signal transformation module calculates the amplitude difference and phase difference between adjacent time series segments in the dynamic characteristics of drilling, compresses the time series segments according to the amplitude difference and phase difference to generate a compressed vibration sequence, and performs Fourier transform on the compressed vibration sequence to generate characteristic spectrum parameters.

[0052] A state correction module projects the characteristic spectrum parameters onto the Grassmann space, calculates the projection eigenvector and eigenvalue to generate the projection matrix parameter; constructs a tangent plane projection function for the projection matrix parameter and calculates the geodesic distance to generate the corrected space coordinates.

[0053] A depth determination module, based on the corrected space coordinates, calculates the position offset and angle offset of the drill string in the three-dimensional space through the tangent plane projection function, substitutes the position offset and angle offset into the Riccati equation for iterative calculation and outputs the depth value to generate real-time depth data.

[0054] The steps for obtaining the original vibration signal are as follows:

[0055] Arrange vibration sensors on the surface of the drill string, collect the spatial position change amount and frequency change amount during the drilling process to obtain the original vibration signal;

[0056] Based on the original vibration signal, perform signal amplification, filtering, and digitization to obtain the preprocessed original vibration signal;

[0057] Based on the preprocessed original vibration signal, convert the continuous vibration signal into discrete data points, perform time series analysis on the discrete data points, extract the dynamic characteristics of the drilling process to obtain the original vibration signal.

[0058] Specifically, based on the vibration sensors arranged on the surface of the drilling tool, the spatial position change and frequency change during the drilling process are recorded. First, a sampling frequency of 100 times per second is set, and a corresponding time identifier is established for each sampling result. The collected amplitude and frequency information are uniformly saved for subsequent comparison. If the recorded amplitude exceeds the maximum allowable value of 3g obtained from prior test experience or the frequency data is higher than the 200Hz threshold determined by comprehensive indoor and outdoor tests, this data is summarized separately as abnormal information. At the same time, when comparing the spatial position change, the set range is 0 to 50 millimeters to determine whether there is an overlimit situation in the small-scale movement of the drilling tool in the horizontal and vertical directions. The comparison standard is determined by the data distribution obtained from long-term field observations and instrument measurements. When the change exceeds 50 millimeters, it is marked as out-of-range data in the summary. All data that meets the above collection range and threshold requirements is directly stored as normal records. After the stage-by-stage collection is completed, all records are checked for duplicates or omissions. The missing parts are resampled to make up for them, and duplicate records are marked to prevent double counting in subsequent statistics. Finally, a complete sequence of spatial position changes and frequency changes is summarized to obtain the original vibration signal.

[0059] Based on the original vibration signal, the recorded amplitude and frequency information are amplified, filtered, and digitized. First, the band-pass filter range is set to 10Hz to 300Hz with reference to the existing on-site noise level, and the frequency components outside this range are attenuated with an attenuation coefficient of 0.1. The filtered signal is digitized with a 16-bit resolution and each data point is matched with the corresponding timestamp. At this stage, it is compared whether the amplified value is still within the pre-set reference range. For example, the amplitude reference interval is set to 0.1g to 2.5g. When it exceeds 2.5g or is lower than 0.1g, it is marked as an abnormal amplification value or an abnormal attenuation value. These abnormal values are compared with the aforementioned original measurement values. If they are still not within the reasonable range after multiple comparisons, they are recorded as information to be investigated. For the normal data points after filtering, they are continued to be retained and arranged in chronological order. After confirming that the data is complete and error-free, they are concentrated to form a continuous preprocessed data sequence, obtaining the preprocessed original vibration signal.

[0060] Based on the preprocessed original vibration signal, to convert the continuous vibration signal into discrete data points and perform time series analysis, it is necessary to first segment the signal samples at a fixed time step. For example, each segment is intercepted every 0.01 seconds, and the average amplitude, peak amplitude, and corresponding main frequency within this period are extracted in sequence. Then, the above data are compared with the pre-set reference intervals respectively. The amplitude reference interval can be set from 0.1g to 2g, and the main frequency reference interval can be set from 10Hz to 200Hz. When it is detected that there is a situation beyond the range or below the lower limit, a suspected anomaly is marked additionally in the corresponding time period. Subsequently, the discrete data points within the normal range are statistically analyzed to obtain the mean and variance of each time period, and these means and variances are used to estimate the change trend of the drilling process. If the trend shows an obvious jump (such as three consecutive samplings all exceeding the mean plus twice the standard deviation calculated before), it can be noted in the result for subsequent detailed study. After finally summarizing all normal segments and abnormal segments, a time series result reflecting the dynamic characteristics of drilling is obtained, and the original vibration signal is obtained.

[0061] The steps to obtain the dynamic characteristics of drilling are as follows:

[0062] Divide the original vibration signal into 256 time series segments, and set a unified sampling interval value for each segment to obtain the segmented time series segments;

[0063] Based on the segmented time series segments, calculate the amplitude and phase of each time series segment to obtain the amplitude-phase characteristics of each time series segment;

[0064] Based on the amplitude-phase characteristics of each time series segment, by analyzing the characteristics of all time series segments, use statistical methods and pattern recognition to identify the dynamic change patterns during the drilling process and generate the dynamic characteristics of drilling.

[0065] Specifically, based on the previously acquired original vibration signals, according to the on-site observation results of the vibration duration and sampling stability during multiple drilling operations, it is determined that the total duration of the complete signal recording is divided into 256 time series segments within the current measurement range. First, calculate the total duration of this recording and divide it by 256 to obtain the time span required for a single sequence segment. Then, set this span as the unified sampling interval value and conduct an inspection. If it is found during the inspection that there are sequence segments with significantly shorter durations or containing burst pulse signals, refer to the critical threshold of 0.02 seconds obtained in laboratory tests. When the duration of a certain segment is lower than this threshold, it is merged with its adjacent segment. If the duration of a certain segment exceeds 0.1 seconds, it is further split into multiple sub-segments. For the situation where the split sub-segments are not sufficient to fill the set interval, zero values are used for padding. After completion of merging and padding, each sequence segment is re-numbered in turn and it is confirmed that the total number remains 256. At the same time, the start and end time tags of each sequence segment are compared to check for any time sequence overlap or discontinuity. When all sequence segments meet the preset length range, they can be regarded as the segmented time series segments.

[0066] Based on the segmented time series segments, several discrete sampling points are selected within each sequence segment to calculate the amplitude and phase. The root mean square amplitude of this sequence segment can be calculated first to represent the overall amplitude level. The specific operation is to sum the squares of the signal intensities of each sampling point within this sequence segment and divide by the number of sampling points, and then take the square root of the result to obtain the root mean square amplitude. If the root mean square amplitude is higher than the upper limit of the safe range summarized in multiple tests, such as 3g, then this segment is marked separately as a suspected overload segment. In addition, the phase information of this sequence segment can be estimated through the mapping of the time domain and the complex plane. For example, the phase angle is obtained according to the ratio of the sine and cosine components of the instantaneous signal. If the phase angle distribution shows extremely discontinuous phenomena, refer to the phase jump threshold set in the laboratory, such as 90 degrees. When the phase jump of a certain segment is greater than 90 degrees, it is recorded as an abnormal jump segment. Finally, the root mean square amplitudes and phase angles of all sequence segments are integrated and summarized and stored uniformly to obtain the amplitude-phase characteristics of each time series segment.

[0067] Based on the amplitude-phase characteristics of each time series segment, it is necessary to statistically summarize and perform pattern recognition on the amplitude and phase data of all segments. First, arrange these amplitude and phase values in sequence by segment and generate a two-dimensional feature matrix. Then, set several reference ranges with reference to the amplitude-phase distribution rules that have appeared in historical drilling data. For example, the amplitude reference range is defined as 0.2g to 2.5g, and the phase reference range is set as 0 degrees to 180 degrees. If the amplitude and phase of some segments fall into the specified reference range at the same time, they are marked as normal segments. If the amplitude exceeds 2.5g or the phase jump exceeds 180 degrees, it is regarded as an abnormal segment and is given key attention in subsequent calculations. After the preliminary classification of normal segments and abnormal segments is completed, a pattern recognition model trained with existing sample data can be used for further recognition. This model adopts a supervised learning method. In the training stage, a sufficient number of labeled samples are collected first. The features of each sample, including amplitude values and phase values, are input. By iteratively updating the model weights, the recognition accuracy reaches the empirically set standard on the validation set. In the inference stage, the amplitude and phase characteristics of each segment are input into the trained model, and the model outputs the specific dynamic change categories. If it is recognized that there is a continuous high-amplitude area or a long-time phase jump, it is listed separately in the final classification result. After summarizing all recognition results, a dynamic change pattern covering the entire drilling process is obtained, and the drilling dynamic characteristics are generated.

[0068] The steps for obtaining the compressed vibration sequence are as follows:

[0069] Extract adjacent time series segments from the drilling dynamic characteristics, calculate the amplitude difference and phase difference for the segments, and capture the vibration changes during the drilling process by comparing the characteristics of each segment to obtain the amplitude difference and phase difference data;

[0070] Based on the amplitude difference and phase difference data, compress the time series segments, merge adjacent and similar-featured segments, and obtain the compressed time series segments;

[0071] Based on the compressed time series segments, re-integrate and construct a continuous time series to obtain the compressed vibration sequence.

[0072] Specifically, extract adjacent time series segments from the drilling dynamic characteristics. First, screen out the time series segments that can form a continuous comparison relationship according to the amplitude and phase distribution information obtained previously, number each segment and read its amplitude and phase values, and then calculate the amplitude difference and phase difference between adjacent segments. The amplitude difference can be expressed as ΔA = A i+1 -A i , and the phase difference can be expressed as Δφ = φ i+1 -φ i, when ΔA or Δφ exceeds the reference range set according to the actual drilling fluctuation degree. For example, the amplitude difference threshold is set to 1.0 g and the phase difference threshold is set to 45 degrees. These reference ranges are obtained by comparing field tests with historical data. If ΔA is greater than 1.0 g, it is marked as a large change. If Δφ exceeds 45 degrees, it is marked as a severe phase jump. After calculating the differences for each sequence segment one by one, these differences are recorded in a unified dataset. During this process, if it is found that there are repeated time periods or missing acquisitions, backtracking and filling are performed according to the original records. For segments with strong noise, they are removed according to the noise threshold obtained from multiple measurements to prevent extreme differences. After completing the difference statistics for all adjacent segments, the amplitude difference and phase difference data are summarized.

[0073] Based on the amplitude difference and phase difference data, data compression is performed on the time series segments according to the pre-established merging criteria. This criterion first divides the sequence segments into two categories: high similarity or low similarity. The judgment of similarity depends on the situation where both the amplitude difference ΔA and the phase difference Δφ fall within the set range. For example, ΔA less than 0.5 g and Δφ lower than 20 degrees are regarded as similar features. If adjacent segments meet the above judgment conditions, they are merged into a new time series segment. When merging, the amplitude and phase values of adjacent segments are weighted and averaged, and the average coefficient is determined by the effective sampling points of each segment. At the same time, the start and end time identifiers of the merged segment need to be retained for subsequent tracing. If the cumulative length of some segments far exceeds the maximum recommended single-segment duration obtained by experience, such as 2 seconds, then it is judged again whether to split them into shorter sub-segments before merging. After merging, they are re-numbered according to the new sequence order, and it is checked whether the amplitude and phase of each merged segment still conform to the reference range. After all checks are correct, the compressed time series segments are obtained.

[0074] Based on the compressed time series segments, these time series are re-integrated and a continuous time series is constructed. The specific method is to first arrange the merged and numbered sequence segments in chronological order, compare the start time of each segment with the end time of the previous segment. If a time stamp jump is found, interpolation is attempted by referring to the original recorded amplitude and phase curves. For example, for the missing amplitude part, the average value of the upper and lower segments is taken, and for the phase, the smoothed transition value calculated from adjacent segments is taken. If there is still a situation that does not conform to the previously obtained amplitude reference interval after interpolation, it is marked as a segment that needs to be reconfirmed and corrected later. When all sequence segments are smoothly connected, they are sequentially spliced to form a complete time series. The amplitude and phase information recorded previously is retained in this time series and made to correspond one by one with the time coordinates. Finally, it is checked whether the sequence length is consistent with the duration of the entire drilling process. If correct, the constructed continuous time series is output to obtain the compressed vibration sequence.

[0075] The steps for obtaining the characteristic spectrum parameters are as follows:

[0076] Extract data from the compressed vibration sequence to obtain data ready for Fourier transform;

[0077] Based on the data ready for Fourier transform, perform Fourier transform to obtain spectral data. The formula is:

[0078]

[0079] where F(k) is the k-th frequency component in the frequency domain, x(j) is the j-th data point in the time series, and N is the total number of data points;

[0080] Based on the spectral data, apply L1 norm constraint to optimize the sparsity of the spectrum to obtain characteristic spectral parameters.

[0081] Specifically, to extract data from the compressed vibration sequence, first read the previously obtained compressed vibration sequence and confirm whether its time coverage is consistent with the total duration recorded during the previous drilling process. If there is missing data or time period gaps, interpolation processing is performed according to the previously recorded sequence segment information. For example, the amplitude is filled with the average value of adjacent sampling segments, and the phase is set with a smooth connection value according to the phase distribution of adjacent segments. After confirming the integrity of the time series, synchronize it with the vibration intensity distribution collected by the sensor according to the time tag corresponding to each record, and screen the records with reference to the effective range determined in advance based on multiple analyses of the drilling rig vibration data. For example, compare the amplitude in the range of 0.1g to 3g, and compare the phase in the range of 0 degrees to 180 degrees. If the amplitude or phase of a record is not within this range, further trace the data of the previous drilling link to locate the possible abnormal segment. After correcting or removing the abnormal segment, form a relatively continuous vibration data sequence without obvious abnormal jitter. Finally, select the corresponding number of data according to the requirements of the number of sampling points for Fourier transform. When the actual available data volume in the compressed vibration sequence is more than the target requirement, relatively stable and representative continuous segments are selected according to the average energy density. When the available data just meets or is slightly less than the target number, it is filled in the way of zero-value extension. After summarizing all the data that meet the requirements, the data ready for Fourier transform can be obtained.

[0082] The benefit of the formula is that by quantifying each frequency component in the frequency domain, the vibration signal originally in the time domain is transformed into the energy distribution in different frequency ranges. This method can, together with the subsequent sparsity optimization step, show the concentration degree of vibration characteristics on various frequency components, making it easier to identify key vibration modes and detect potential abnormalities in the subsequent processes.

[0083] The steps to obtain the x(j) parameter are:

[0084] Here, \(x(j)\) represents the vibration amplitude at the \(j\)-th time sampling point, which is usually measured by a vibration sensor arranged on the surface of the drill string. The sensor obtains the original vibration at a sampling frequency higher than 200 Hz in the previous process, and after band-pass filtering and digital processing, a discretized amplitude value is obtained. In practice, the sensor will be calibrated to obtain the output voltage range corresponding to each 1 g, and then combined with the hardware acquisition frequency and range conversion coefficient, the voltage value measured each time is converted into a vibration amplitude expressed in the physical unit g. Then, combined with the time series record that has been compressed and interpolated and supplemented before, they are uniformly numbered 0, 1, 2, …, N - 1. Finally, in this formula, \(x(0), x(1), …, x(N - 1)\) are read in sequence to complete the subsequent frequency domain analysis. For example, in a certain measurement, \(x(0)=0.25g\), \(x(1)=0.27g\), …, \(x(1023)=1.32g\) can be obtained. These data are obtained from the long-term operation of the sensor and are formed by the accumulation of multiple vibration samples per second of the drill string. The parts with abnormal peaks have been removed or corrected in the previous process.

[0085] The steps to obtain the parameter \(j\) are as follows:

[0086] Here, \(j\) is the index of the sampling point in the time series, ranging from 0 to N - 1, and its value is determined by the number of sampling points in the compressed vibration sequence. This index does not directly represent the time length, but is used to correspond each sampling point to the corresponding amplitude value one by one. For example, when the sensor samples at 512 Hz and records data for 2 seconds, about 1024 sampling points will be obtained. At this time, N can be taken as 1024, and then \(j\) will range from 0 to 1023, and each \(j\) corresponds to the vibration amplitude collected instantaneously within this time period. In the specific implementation, the sampling time will be marked in milliseconds and then remapped into an integer index value \(j\) to facilitate directly using an integer to quickly locate the sampling point during the formula operation.

[0087] The steps to obtain the parameter \(N\) are as follows:

[0088] Here, \(N\) represents the total number of data points, which is used to control the length of the discrete Fourier transform. When the sampling frequency and sampling duration are given, \(N\) can be determined. For example, during the repeated measurement at the drill rig site, sampling 512 times per second and lasting for 2 seconds can obtain 1024 points, so \(N = 1024\). At this time, if in order to adapt to the fast Fourier transform algorithm for subsequent processing, the data can also be interpolated or reduced to ensure that \(N\) is an integer power of 2. At the same time, \(N\) cannot be changed arbitrarily in the whole data segment and needs to correspond one by one to the previously obtained vibration amplitude sequence. If only 1000 effective sampling points are obtained through calculation in the previous compression or interpolation stage, interpolation operations will also be used to expand it to 1024 points. In this way, each newly inserted point will follow the average value of adjacent measured values or a smoothly transitional value to keep the data coherent.

[0089] The steps to obtain the parameter \(k\) are as follows:

[0090] Here, k increases from 0 to N - 1. Specifically, it represents the frequency index. The corresponding actual frequency can be calculated through When the sampling frequency f s is known, the physical frequency position corresponding to each k can be obtained. For example, when the sampling frequency f s = 512 Hz, k = 0 corresponds to 0 Hz, k = 1 corresponds to 0.5 Hz, and up to k = 1023 corresponds to approximately 511.5 Hz. Here, k is not simply a virtual serial number, but an important index for subsequent judgment of which frequency ranges the main vibration energy is in. When monitoring the vibration characteristics of the drill string, if it is found that the amplitude corresponding to a certain k is significantly higher than other frequencies, then in-depth statistics or troubleshooting can be focused on this frequency segment. The accuracy of k is jointly determined by N and f s .

[0091] Calculation process:

[0092] First step, select N = 1024, and read the sampling amplitudes numbered from 0 to 1023 into x(j). For example, select 5 data points for illustration:

[0093] x(0) = 0.25g, x(1) = 0.27g, x(2) = 0.31g, x(3) = 0.36g, x(4) = 0.42g

[0094] Second step, for the specified k, calculate the exponential factor e -i2πkj / N . For example, when k = 4, there is

[0095] e -i2π·4·j / 1024

[0096] Third step, multiply the corresponding x(j) by this exponential factor term by term and accumulate. If only the first five terms are illustrated, we get

[0097]

[0098] Fourth step, continue to perform the same operation for all points where j ranges from 0 to 1023. The accumulated result is F(4).

[0099] Fifth step, if the entire frequency spectrum needs to be calculated, let k traverse from 0 to 1023 to obtain F(0), F(1), …, F(1023) respectively. These results constitute all the components in the frequency domain.

[0100] The results show that for a specific k, the intensity and phase information of the vibration signal can be obtained in the frequency domain through such continuous accumulation and complex exponential multiplication operations. After all the actual operations are completed, an amplitude-frequency distribution curve and a phase-frequency distribution curve will be formed. If the amplitude corresponding to a certain k far exceeds the other components, it means that the energy distribution in this frequency range is relatively concentrated. Subsequently, the contribution of this frequency component to the vibration and whether there is vibration abnormality can be judged in combination with the actual working state of the drill string.

[0101] Based on the spectral data, first arrange and statistically analyze all k values according to the frequency domain components and phase information obtained previously, and at the same time read the corresponding amplitude data to establish a frequency-amplitude mapping sequence. In order to implement the L1 norm constraint, a certain proportion of data points will be selected as the key points for sparsification processing. For example, compare the amplitude values in the interval from 0 dB to 60 dB and mark the points with significant increase. If a certain point exceeds the threshold of 40 dB formed by multiple tests, it will be added to the key sparse list. At this time, combined with the noise energy determination distribution obtained by frequency domain analysis in the previous text, weaken the weight of the relatively low-energy and discrete frequency points in the way of L1 constraint. When implementing this constraint, the corresponding parameters will be gradually adjusted according to the compared noise distribution, and record the current amplitude correction result for each frequency point. After completing the sparsification adjustment of all frequency points, a new spectral distribution data will be generated uniformly and interpolated or smoothly transitioned. Finally, all the frequency components that meet the set threshold conditions after adjustment will be summarized to generate characteristic spectral parameters.

[0102] The steps to obtain the projection matrix parameters are as follows:

[0103] Based on the characteristic spectral parameters, extract the frequency component matrix and frequency characteristic data, and perform a centering process on the frequency component matrix to generate a centered frequency component matrix;

[0104] Based on the centered frequency component matrix, calculate the projection matrix, and the expression is:

[0105]

[0106] where, M ij is the element in the i-th row and j-th column of the projection matrix, f i and f j are the amplitudes of the i-th and j-th frequency components in the frequency component matrix respectively, u i and u j are the eigenvectors of the i-th and j-th frequency components respectively, λ i and λ j are the eigenvalues of the i-th and j-th frequency components respectively, θ i and θ j are the phase angles of the i-th and j-th frequency components respectively, φ iand φ j are the amplitude angles of the i-th and j-th frequency components respectively;

[0107] Based on the projection matrix, it is mapped to the Grassmann space, and geometric operations are performed through eigenvalue decomposition to obtain eigenvectors and eigenvalues, forming the projection matrix parameters.

[0108] Specifically, based on the characteristic spectrum parameters, first extract all frequency indices and corresponding amplitude values from the previously obtained characteristic frequency components and phase data, organize these amplitude values and corresponding frequency indices into a preliminary two-dimensional matrix form, and then check whether there is any missing or inconsistent amplitude and phase information recorded previously. If complete data cannot be matched at a certain frequency, interpolation is performed with reference to the previously obtained amplitude or phase sequence. For example, at the frequency point where the amplitude is missing, it is filled with the average amplitude of the two adjacent frequency components, or at adjacent frequency points, observe their amplitude change trends and calculate the smooth transition values. When all frequencies are filled, a frequency component matrix can be formed. Subsequently, according to the reference standard established in multiple on-site measurements, perform a centering process on this frequency component matrix, subtract the average value of each row or column in the matrix. If the average value exceeds the amplitude threshold of 0.5 obtained by comprehensive analysis of multiple data analyses, it is recorded as a high-value offset, and if it is lower than 0.1, it is recorded as a low-value offset. During this process, these offset values will be recorded for subsequent comparison. Finally, after completing the centering, a new centered frequency component matrix is obtained, and each element of this matrix represents the numerical result after removing the average offset at the amplitude and phase levels of the corresponding frequency.

[0109] The advantage of the formula is that through the simultaneous operations of multiple factors such as amplitude, eigenvector, phase angle, and amplitude angle, a more comprehensive correlation measure can be established between frequency components. And with the addition of eigenvalues, the weights of different frequency components in the overall vibration mode can be more intuitively reflected, which has more comprehensive reference significance for the subsequent evaluation of the vibration characteristics of the drill in three-dimensional space or higher dimensions.

[0110] |f i The steps for obtaining the | parameter are as follows:

[0111] This parameter represents the amplitude of the i-th frequency component, and its value is extracted from the characteristic frequency components obtained through Fourier transform and sparsity optimization in the previous stage. In the specific acquisition process, first, the vibration of the drill is cumulatively measured based on multiple sampling points, and the time-domain signal is processed through the previously mentioned discrete Fourier transform to form a spectrum. Then, the L1 norm constraint is applied to the amplitude of each frequency component in the spectrum to reduce the influence of noise and insignificant frequencies. Subsequently, while retaining the main energy components, a relatively clear set of amplitude values is obtained. Finally, the amplitude of the i-th frequency component is directly read from this set as |f i|, in the vibration scenario of an industrial drill rig, after retrieval, the common amplitude can fluctuate within the range of 0.1g to 5g. When used for further calculations, the unit needs to be kept consistent with other parameters or necessary dimensional conversions are carried out. For example, convert 0.1g to 5g to a dimensionless value for formula operations.

[0112] |f j |The steps for obtaining the parameter are as follows:

[0113] Same as |f i |, |f j | represents the amplitude of the j-th frequency component. The acquisition method used is the same as that of the i-th component. First, perform a frequency-domain transformation, then use the amplitude information recorded previously, and then combine the upper and lower amplitude limits during on-site measurement for secondary verification. If the amplitude value of a certain frequency component in continuous sampling is higher than 5g, the on-site vibration record will be traced back for inspection to confirm whether there is an instantaneous peak caused by an impact load. Only after determining that it is an effective measurement will it be included in the final data. After determining this value, |f j | can be used as an amplitude quantity of the same status as |f i | to enter the subsequent formula operations.

[0114] u i The steps for obtaining the parameter are as follows:

[0115] This parameter represents the eigenvector of the i-th frequency component. Its source is often obtained by vectorizing the eigenmodes corresponding to different frequency components when performing principal component or singular value decomposition on the spectral matrix. For example, in on-site measurement, each frequency component corresponds to dozens of sampling results. By performing correlation analysis and covariance calculation on these results, the corresponding eigenvector can be decomposed to reveal the common distribution of this frequency in multiple measurement scenarios. This vector is usually recorded in real number form, and the number ranges from 1 to several hundred components, depending on the sampling dimension. In actual data processing, the eigenvectors of the same frequency after multiple measurements are superimposed and averaged or weighted averaged, and the weight values come from the signal-to-noise ratios of each measurement round. Finally, u is determined. i And normalized to, for example, between 0.0 and 1.0 for operations in the formula.

[0116] u j The steps for obtaining the parameter are as follows:

[0117] As the eigenvector of the j-th frequency component, the acquisition method is the same as that of u i |. It is necessary to extract the response vectors of the j-th frequency component under multiple measurements and perform the same normalization process. The numerical range is usually the same as that of u iKeep them unified, for example, all within the range of 0.0 to 1.0. If it is found in actual operation that a certain eigenvector shows an extremely high value relative to other frequency components, further investigation will be carried out in combination with the noise distribution and instrument calibration data recorded previously to determine whether the eigenvector is real and effective. After confirmation, it will be incorporated into the subsequent calculation to finally form a complete u j numerical sequence.

[0118] λ i The steps to obtain the parameter are as follows:

[0119] This parameter is the eigenvalue of the i-th frequency component, which is mainly obtained by performing eigenvalue decomposition or singular value decomposition on the frequency component matrix. When the matrix dimension is large, multiple groups of measurement data will be integrated to obtain the eigenvalue, and the stability of the eigenvalue under different measurement dates or working conditions will be compared. In the case of industrial drilling sites, common λ i will be in the range of 0.1 to 10.0. The larger this value is, the more significant the proportion of this frequency component in the corresponding eigenmode usually means. The acquisition of the eigenvalue is inseparable from the covariance or correlation measurement results corresponding to each frequency component. By combining the eigenvector and the eigenvalue, the energy distribution of the spectrum and its main influencing factors can be quantified.

[0120] λ j The steps to obtain the parameter are as follows:

[0121] Similar to λ i λ j is the eigenvalue of the j-th frequency component, which is obtained by a symmetric covariance matrix or other mathematical means in the same decomposition process. When the value of λ j is significantly higher than other frequency components, it is necessary to check whether there is noise concentration or mechanical resonance, etc. After determining its effectiveness, λ j will also be normalized or kept in actual dimension for calculation. In actual industrial scenarios, if multiple rounds of measurement results are used to superimpose to obtain a more stable eigenvalue, the multiple rounds of data will be spliced, and a variance-based weighting process will be performed on each frequency entry to form the final λ j sequence and used in subsequent projection matrix operations.

[0122] θ i The steps to obtain the parameter are as follows:

[0123] This parameter represents the phase angle of the \(i\)-th frequency component. The value is obtained by performing the arctangent operation on the complex form of the frequency-domain component during the discrete Fourier transform mentioned above. At the same time, combined with the noise clipping strategy, it is ensured that the phase distribution is within the range of 0 degrees to 180 degrees or 0 degrees to 360 degrees. In industrial drill rig vibration monitoring, it is more convenient to analyze when the common phase angle is concentrated between 0 degrees and 180 degrees. Therefore, the range of 0 to 180 degrees will be selected as the recording interval. If the phase angle is greater than 180 degrees during a certain measurement period, appropriate mapping or period reduction processing will be performed to keep the data of all measurement periods comparable, and finally \(\theta\) is formed. i and saved at the position corresponding to \(f\) i

[0124] \(\theta\) j The steps to obtain the parameter are as follows:

[0125] Similar to \(\theta\), \(\theta\) i represents the phase angle of the \(j\)-th frequency component. The acquisition method is similar. First, the phase of the frequency-domain component is obtained, and then the phase is normalized. If the measurement result exceeds 180 degrees, it is mapped back to the range of 0 to 180 degrees. In the final component list, \(\theta\) j is usually used together with \(\theta\) j In the calculation of the projection matrix, the difference in phase between these two components can be reflected by \(\cos(\theta\) i -\(\theta\) i ). When the difference is too large, it indicates that the peak points of the two frequency components are inconsistent during the vibration process or there is an obvious phase misalignment. j

[0126] \(\varphi\) i The steps to obtain the parameter are as follows:

[0127] This parameter represents the amplitude angle of the \(i\)-th frequency component. Different from the amplitude value \(|f|\) i , the amplitude angle is more inclined to quantify the directionality of the amplitude distribution. In drill string vibration analysis, sometimes the amplitude vector is projected in polar coordinates or spherical coordinates to obtain an angle with geometric meaning. For example, when recording vibration components in three-dimensional space, the polar angle and azimuth angle are calculated based on the amplitude ratio of each axis, and the combination of the two forms the amplitude angle. The value is usually in the range of 0 degrees to 90 degrees. When the amplitude angle is too high, it indicates that the vibration is mainly concentrated in a specific direction. When obtaining \(\varphi\) i , it is necessary to first sum up the values of the sensors on each axis, then calculate through trigonometric functions, and finally summarize at the position corresponding to the \(i\)-th frequency component.

[0128] \(\varphi\) j The steps to obtain the parameter are as follows:

[0129] Similar to \(\varphi\), \(\varphi\) i j ​​​Represents the amplitude angle of the j-th frequency component. In the calculation, it is also necessary to synthesize the multi-axis vibration measurement data, and then use the geometric projection method to obtain the angle relative to a certain main axis or reference plane. There will be certain fluctuations in the values that appear in multiple measurements, and a more stable φ can be determined by statistical median or weighted average. j For feature analysis. In the feature parameter list, this amplitude angle and θ j are used together so that both the phase difference and amplitude distribution information of the projection matrix can be incorporated into the calculation.

[0130] Calculation process:

[0131] First step, given the specific values of each parameter:

[0132] |f i | = 1.20, |f j | = 0.95, u i = 0.65, u j = 0.72,

[0133] λ i = 2.10, λ j = 3.20, θ i = 50°, θ j = 70°,

[0134] φ i = 20°, φ j = 30°

[0135] Second step, calculate the denominator

[0136]

[0137] Third step, calculate the numerator |f i |·|f j |:

[0138] 1.20 × 0.95 = 1.14

[0139] Fourth step, calculate cos(θ i - θ j ):

[0140] θ i - θ j = 50° - 70° = -20°, cos(-20°) = cos(20°) ≈ 0.94

[0141] Fifth step, calculate φ i ·φ j and divide it by λ i + λ j |:

[0142] φ i ·φ j = 20° × 30°

[0143] Since the amplitude angle is often converted to radians or a conversion ratio is set in actual calculations, it can be regarded as a numerical ratio of 20 × 30 = 600 here;

[0144] λ i +λ j = 2.10 + 3.20 = 5.30

[0145]

[0146] Step 6, substitute the above results into the original formula:

[0147]

[0148] This result indicates that between the selected frequency components i and j, the projection matrix element M calculated by the combined action of each parameter ij is approximately 107.596. The larger this value, the more significant the amplitude product and phase coupling of these two frequency components are. When performing eigenvalue decomposition on the entire projection matrix subsequently, the frequency pair with a larger M ij may dominate. Subsequently, integrating this value into the geometric operations mapped to the Grassmann space can ultimately better extract or distinguish the characteristic patterns of different frequency components.

[0149] Based on the projection matrix, first read the value ranges of each element in the previously calculated projection matrix, summarize them by rows and columns and ensure that the number of rows and columns is consistent with the number of frequency components constructed previously. Then, for a projection matrix with a large number of rows and columns, eigenvalue decomposition will be used to analyze the main eigenvalues and corresponding eigenvectors one by one. At this time, it is necessary to ensure that the input matrix is symmetric or can be symmetrized. If it is detected that some row and column components exceed the upper threshold of 2000 summarized through multiple measurements, then trace back to check the corresponding frequency components, and include them in the calculation after confirming that the data is correct. Arrange the identified main eigenvalues and secondary eigenvalues item by item, and record the corresponding eigenvectors for subsequent geometric mapping in the Grassmann space. When mapping, the eigenvectors will be represented as the coordinate vectors of the basis vectors in this space, and each row and column will be mapped into corresponding subspace points. Then, perform geometric calculations on these points through geodesic or other distance algorithms. During the processing stage, for high-dimensional matrices, block operations may be required. The matrix is separated into manageable sub-blocks by rows or columns and then eigenvalue decomposition is performed. After the block operation is completed, the overall eigenvalues and eigenvector list are unified and merged. Finally, these data are summarized to obtain the new projection matrix parameters, which usually include several numerical tables in row and column form and the corresponding main eigenvalue entries, record their order and numerical size, so as to complete the geometric operation and analysis of the mapping results of the entire frequency component.

[0150] The steps to obtain the corrected space coordinates are as follows:

[0151] Based on the projection matrix parameters, construct a tangent plane projection function, and the formula is:

[0152]

[0153] where x and y respectively represent the coordinate values of any two points in the projection matrix, and F(x, y) is the tangent plane projection function;

[0154] Use the tangent plane projection function to calculate the distances of all points in the projection matrix, and map the distance of each pair of points to a new distance value;

[0155] Through the new distance value, calculate the geodesic distance between data points, and use the geodesic distance to correct the space coordinates to obtain the corrected space coordinates.

[0156] Specifically, the benefit of the formula is that by combining the product of x and y in the numerator part, and including The two factors of exp(-|x - y|) can incorporate the coupling relationship between the magnitude and the relative difference within the same function framework. This can take into account the influence of both the magnitude and the difference in various industrial scenarios such as spatial mapping or data measurement. It can also be further combined with geodesic distance or other geometric calculations to achieve the purpose of comprehensively considering the deviation degree and the amplitude weight, so as to accurately reflect the gap degree between data points in subsequent projection and correction.

[0157] The steps to obtain the x parameter are as follows:

[0158] This parameter represents the coordinate value of a certain point in the projection matrix. The value comes from the definition of the coordinate axes in the previously obtained projection matrix parameters. In the vibration analysis of industrial drills, the mapped data is usually distributed in a coordinate system of one or more dimensions, and x is the coordinate on the horizontal or any arbitrarily selected axis. To obtain x, it is necessary to read the previously recorded projection matrix row by row or column by column, then correspond to the specific point position according to its row and column indices, and then extract the matrix value at that position. This value may be derived from the projection coordinates obtained by combining operations such as eigenvalues and eigenvectors. In actual monitoring, the common coordinate range usually fluctuates between -10 and 10. This interval can be set with reference to the large vibration dimension at the site. Then, it is necessary to check whether there is overflow or missing. If it is found that some coordinates fall within the interval of -1.0 to 1.0 after normalization in the previous steps, then x at this time also needs to be processed with the same normalization standard. For example, if the actual dimension is between -8.5 and +9.2, it will be mapped to between -0.85 and 0.92, and finally an x value adjusted by a unified scale is obtained.

[0159] The steps to obtain the y parameter are as follows:

[0160] Similar to x, the y parameter represents the position of the same point in the projection matrix in another coordinate dimension, also from the previously obtained projection matrix parameters. After projection or eigen-decomposition, each point will have multi-dimensional coordinates, where y can be regarded as the coordinate in the vertical or orthogonal direction. The acquisition method is usually to retrieve the same matrix element at another row and column index. In the application scenario of industrial drills, if the vibration is mainly concentrated in certain frequency bands or some eigenvalues are particularly prominent, it will also cause the value of y to be relatively large. Through multiple on-site records, it is found that the value fluctuates within the range of -10 to 10 relatively commonly. It should also be noted that after normalization, its interval may become -1.0 to 1.0. To ensure consistent operations, the same coordinate transformation standard as x needs to be adopted. If it is found that the y value at a certain time has exceeded the reasonable range formed by the previous cumulative analysis, such as being higher than 10 or lower than -10, then it is necessary to trace back to the previous steps and verify its validity in combination with the sensor range and the eigenvector decomposition results. Finally, after retaining or revising the y value, subsequent operations can be carried out.

[0161] The steps to obtain the F(x, y) parameter are as follows:

[0162] This parameter is the output value of the tangent plane projection function. It is not directly read from the projection matrix but is obtained by substituting both x and y into the formula. F(x, y) is often regarded as a type of modified distance or similarity metric in spatial mapping and measure. After the previous confirmation of x and y, the item-by-item calculation of the numerator and denominator needs to be performed. In industrial applications, to prevent the denominator from being too small or too large, observation intervals are set for and exp(-|x - y|) respectively. For example, if is less than 0.01 or greater than 20, the previously obtained sensor calibration data will be called for troubleshooting and adjustment. Similarly, exp(-|x - y|) may decay to a minimum value when |x - y| is too large and approach 1 when |x - y| is extremely small, and it needs to be adapted to the actual coordinate range. The finally obtained F(x, y) can directly enter the subsequent distance or geodesic calculation.

[0163] Calculation process:

[0164] First step, select specific values. For example, let x = 2.5 and y = 3.2. Then, according to the known environmental measurement, the values of the corresponding points of these two coordinates in the projection matrix are confirmed to be 2.5 and 3.2, and it is confirmed that these two values are between -10 and 10, meeting the previously set coordinate range.

[0165] Second step, calculate the numerator x·y:

[0166] x·y = 2.5 × 3.2 = 8.0

[0167] Third step, calculate the denominator part

[0168] x 2 + y 2 = (2.5) 2 +(3.2) 2 = 6.25 + 10.24 = 16.49

[0169]

[0170] Fourth step, calculate exp(-|x - y|):

[0171] |x - y| = |2.5 - 3.2| = 0.7, exp(-0.7) ≈ 0.4966

[0172] Fifth step, add the two terms to get the overall denominator:

[0173] 4.06 + 0.4966 ≈ 4.5566

[0174] Sixth step, finally obtain F(x, y):

[0175]

[0176] This result shows that in the case of the projection matrix coordinates x = 2.5 and y = 3.2, using this tangent plane projection function, a numerical value of approximately 1.756 can be obtained. The practical significance of this numerical value is that, compared with the situation of simply using the Euclidean distance or only using the absolute difference, this function comprehensively considers factors such as the amplitude multiplication and difference attenuation between points. When the obtained F(x, y) is greater than, for example, 2.0, it means that the relative relationship between these two points tends to be highly coupled under this projection function. When it is less than 1.0, it indicates that the two points are not prominent in terms of amplitude product or difference degree. By combining this result with the subsequent geodesic distance, the spatial coordinates can be further corrected and a more refined characterization of the inter-point correlation can be obtained.

[0177] Use the tangent plane projection function to calculate the distances for all points in the projection matrix. During the expansion process, it is necessary to first read the list of projection matrix parameters formed previously and confirm that the number of rows and columns is consistent with the number of frequency components or coordinate components. For each point, mark its coordinates x and y, and then pair them with the coordinates of other points in the matrix in sequence according to the same coordinate index order to form all possible point pairs. For each point pair, it is necessary to determine the difference in the horizontal value and the difference in the vertical value between the two coordinates. Here, all the absolute differences of x and y are compared with the range established in multiple rounds of on-site observations. For example, if |x a -x b | is greater than 10 and the two points are in extremely distant positions in the projection space, they are marked as a height difference pair in the record. If |x a -x b | is less than 0.5 and |y a -y b|If it is also less than 0.5, it is marked as an approximate pair. Then, each pair of coordinates is substituted into the tangent plane projection function for term-by-term operations of the numerator and denominator, and combined with the exponential mapping of the difference part for attenuation or amplification. If some differences significantly exceed the previous empirical threshold, such as 15, this point pair needs to be checked again for anomalies or missing data in the original projection matrix. After confirmation, the projection function output of the current point pair can be obtained. Then, the results of all point pairs are stored to form a preliminary distance mapping table. During this process, the result values appearing in the distance mapping table are also compared with the pre-established valid range. For example, when the distance mapping result is between 0 and 10, it is regarded as the normal range. If it exceeds 10, it is recorded as an ultra-high mapping value and needs to be reviewed separately during subsequent classification. Finally, after traversing all point pairs, a new set of distance values can be obtained, which covers the projection function relationship degrees between all points in the projection matrix. Subsequently, these new distance values are normalized, compressing the highest value to, for example, 10, and raising the lowest value to, for example, 0.1, and combined with the original coordinates of each point in the projection matrix to form a multi-dimensional table. This table can be used to continue performing geometric calculations or aggregating geodesic distances, thereby presenting the distance situation between points in a more refined manner in the projection space, and obtaining a new distance mapping distribution.

[0178] With the new distance values, first decompose the obtained distance mapping table in a row-column manner, store and summarize the distance data corresponding to each point index. If the distance value between any pair of points is higher than the upper limit of the segmented range formed previously, for example, 10, then check whether there is a significant coordinate difference between this pair of points, and when the true deviation is confirmed, classify this pair of points into the long-distance group. If the distance value between any pair of points is lower than 0.1, it is regarded as the adjacent group. When the preliminary division of the distances between all points is completed, it enters the geodesic distance calculation link. Here, the connection paths of each point will be statistically analyzed in the projection space in the form of a surface or a manifold. Accumulate the total distance of each point communicating with other points on multiple possible connection lines and use multiple iterations to find the overall shortest path on this basis. During the process, first mark the long-distance pairs higher than 10 or the adjacent pairs lower than 0.1 respectively to accelerate the investigation. Record the path of each pair of point connections and compare it with the reference value obtained by laser ranging or ultrasonic measurement on-site to ensure that there is no dimensional mismatch. Then extract the geodesic according to the principle of minimum energy consumption on the surface or manifold. When the corresponding geodesic length is found, it can be used as the true measurement value between this point and the target point. If the lengths of some geodesics deviate too much from the external reference value, return to modify the previous normalization step or check whether abnormal coordinates need to be excluded. Finally, form the geodesic distance matrix of each point in the projection space. After obtaining all the geodesic distances, reposition the correction values of each point in the three-dimensional or higher-dimensional space coordinates according to these distances. Compare the original projection coordinates with the correction values one by one and perform coordinate updates. For points with significant deviations, the vibration characteristics or sensor status can be further traced back in combination with the actual measurement log. For points whose corrected position offset reaches the set threshold, they can also be included in the next round of tracking records. After completing the overall coordinate correction, the final corrected space coordinates can be obtained.

[0179] The steps for obtaining real-time depth data are as follows:

[0180] Based on the corrected space coordinates, use the tangent plane projection function to calculate the position offset and angle offset of the drill string in the three-dimensional space;

[0181] Based on the position offset and angle offset, calculate the depth adjustment value, and the calculation formula is:

[0182]

[0183] where, Δx and Δy represent the position offsets of the drill string in the horizontal and vertical directions, α is the total angle offset of the drill string, and β is the angle with the vertical axis;

[0184] Based on the depth adjustment value, substitute the depth adjustment value into the Riccati equation, perform iterative calculations to adjust the depth value of the drill string, and generate real-time depth data.

[0185] Specifically, based on the calibrated spatial coordinates, first read the coordinate values obtained through projection and geodesic calculation in the previous text. Determine the specific position of the drill string in the three-dimensional space according to these coordinates. Then, compare with the reference intervals established in advance through multiple on-site measurements and industrial analysis to evaluate whether there is a large deviation between each coordinate point and the actual position of the drill string. If it is found that the deviation exceeds the critical threshold set by empirical statistics, such as 0.5 meters, an offset record will be noted in the result. Next, read the azimuth angle information recorded at this coordinate. These angle information can be obtained by the angle measurement units arranged in different azimuths of the drill string. In the statistical process, each azimuth data will be checked item by item to exclude abnormal acquisition values caused by factors such as sensor interference or cable faults. If the value is within the common range, such as 0 degrees to 90 degrees, it is considered normal. If an angle higher than 90 degrees appears, continue to retrieve the sensor distribution and drill string attitude mentioned above to confirm its validity. Once the angle value does not exceed the feasible range in the existing records, it will be incorporated into the subsequent summary. Subsequently, combine the collected position coordinates and angle information, and perform label matching on them to correspond to the timestamp and operation sequence of the same drilling. On this basis, use the tangent plane projection function to evaluate the comprehensive difference degree of these positions and angles. First, extract the relative displacement between the coordinates for each measurement period, and record whether it exceeds the reference standard summarized through multiple on-site measurements in the horizontal and vertical axes. For example, when the horizontal displacement is higher than 0.2 meters or the vertical displacement exceeds 0.3 meters, the corresponding period will be listed as a high-offset interval. Compare with the angle information to check the included angle with the reference direction to see if it exceeds the 15-degree safety upper and lower limits defined in the drill rig design. Finally, through these comparisons, the actual position offset and angle offset of the drill string can be refined, and complete reference data can be provided for the subsequent depth adjustment link in the system.

[0186] The benefit of the formula is that by simultaneously incorporating the horizontal and vertical displacement Δx and Δy, the total angle offset α, and the angle β with the vertical axis, multiple factors such as position and attitude can be comprehensively considered to evaluate the correction amplitude for depth, and trigonometric operations are incorporated into it, making the depth adjustment value Δz not only reflect the offset degree of the drill string in the plane but also consider the influence of the attitude angle on the actual depth, thus being more in line with the multi-dimensional motion state in the actual drilling environment.

[0187] The steps to obtain the Δx parameter are as follows:

[0188] This parameter represents the actual displacement of the drill string in the horizontal direction. The value is sourced from the three-dimensional space correction coordinates constructed previously and verified against the comparison records of on-site laser rangefinders or ground surveying instruments. To form an accurate Δx, the lateral coordinates of the drill string recorded on the same horizontal reference line during multiple rounds of operations are compared item by item. For example, at the initial stage, the origin is aligned when the equipment is in the calibrated position, and subsequent horizontal direction changes at each time period are measured based on this origin. If it is found that any time period exceeds the horizontal safety threshold of 0.5 meters obtained from the summary of construction experience and equipment specifications during the process, the instrument output curve is checked and the displacement data collected by other auxiliary sensors is referred to. When it is confirmed to be valid, it is incorporated into the Δx sequence. The values measured multiple times are also denoised or sensor instantaneous jumps are filtered to form a smoother average or median as the final Δx. When all horizontal coordinate data is differentiated from the reference position, a series of Δx values can be aggregated. Usually, at the operation site, Δx fluctuates within the range of 0 meters to 1 meter, specifically varying depending on geological conditions and drill string models.

[0189] The steps to obtain the Δy parameter are as follows:

[0190] This parameter is the displacement of the drill string in the vertical direction. The method of obtaining it is similar to that of Δx. First, the value of this point on the vertical axis is queried by referring to the three-dimensional correction space coordinates formed previously, and then this value is differentiated from the initial reference position to obtain Δy. In industrial drilling operations, vertical deviation is often affected by factors such as well inclination and formation structure. Therefore, an acceptable vertical deviation range is specified during equipment joint debugging. For example, when the logging depth is 1000 meters, a vertical deviation of 0 meters to 5 meters is generally regarded as the normal range. If the collected Δy is greater than 5 meters, it is compared with historical cases, and the reliability of this value is confirmed based on the long-term stability test results of the instrument or the feedback of other downhole positioning systems. If it meets the authenticity requirements, this value is written into the Δy data sequence. In subsequent data processing, the vibration characteristics obtained previously may also be combined to determine whether this deviation is related to formation changes, and finally a more robust set of vertical displacement values is formed.

[0191] The steps to obtain the α parameter are as follows:

[0192] This parameter is the total angular offset of the drill string. Its value is measured by a downhole triaxial gyroscope or a multi-azimuth inclination sensor, and is comprehensively obtained by combining the attitude information when mapping to the three-dimensional space before. In on-site implementation, angle data is often collected every 10 to 20 meters of footage. These data lists are corresponding to timestamps. After stagewise smoothing filtering and calibration, they are unified into the same reference coordinate system. Then, after dimension matching, when the output unit of the gyroscope is inconsistent with the projection coordinate system, the angle value needs to be converted. Thus, a unified α distribution list is formed. Then, compare with the reference direction calibrated in the projection and geodesic distance stage in the previous text, calculate the angle between this direction and the current gyroscope measurement direction and further verify it. If there is an obvious jump in a certain measurement, for example, greater than 10 degrees, then retrieve the load and vibration records in the relevant period to exclude abnormal interference. If it is confirmed to be correct, then summarize it into the final sequence of α. In engineering practice, α is mostly in the range of 0 degrees to 30 degrees, depending on the geology and the design of the drill string. Once it exceeds this range, it indicates that there is an excessive attitude deflection, and at the same time, it is entered into the subsequent depth calculation link.

[0193] The steps to obtain the β parameter are as follows:

[0194] This parameter represents the angle between the current attitude of the drill string and the vertical axis. The value mainly comes from the inclination measurement device installed near the drill bit or drill collar. Each measurement is distributed at multiple points around to improve accuracy. The inclination information of these points is summarized and the average angle relative to the ground vertical line is calculated, and then coordinate transformation is performed with the previously recorded coordinate system and stored as β. In on-site acquisition, β is often between 0 degrees and 90 degrees. Exceeding 90 degrees usually means that the direction of the drill string has been severely inclined or reversed, and it is necessary to immediately check the formation conditions and the surrounding support status. The β data accumulated by this method will be continuously refreshed during multiple rounds of drilling, and provide the basis for attitude correction when finally entering the depth adjustment calculation. If there are omissions in the previous stage of measurement or attitude estimation, additional acquisition is required to improve the β sequence to ensure that continuous and reliable angle information can be read in subsequent formula operations.

[0195] Calculation process:

[0196] First step, select specific values. For example, let Δx = 0.20 m, Δy = 0.15 m, and confirm that it is within the allowable range of 0 m to 1 m with reference to the previous instrument measurement data. Take α = 12°, β = 8°. These two angles are both obtained by summarizing the records of the aforementioned triaxial gyroscope and inclination sensor, and the actual geological conditions and equipment attitude have been verified.

[0197] Second step, calculate Δx 2 +Δy 2 :

[0198] 0.20 2 +0.152 = 0.04 + 0.0225 = 0.0625

[0199] Step 3, calculate tan 2 (α). When α = 12°,

[0200] tan(12°) ≈ 0.2126, tan 2 (12°) ≈ 0.0452

[0201] Step 4, calculate cos(β). When β = 8°,

[0202] cos(8°) ≈ 0.9903

[0203] Step 5, combine the product of these three parts:

[0204] (0.0625) × (0.0452) × (0.9903) ≈ 0.00280

[0205] Step 6, take the square root of the result:

[0206]

[0207] This result indicates that under the above measured conditions of Δx = 0.20 m, Δy = 0.15 m, α = 12°, β = 8°, etc., the depth adjustment value Δz is approximately 0.0529 m. This value represents the amplitude of the additional correction required for the depth of the drill string in the current attitude. If this value is greater than, for example, 0.1 m, it is considered that the drill string has a significant deviation and may require further correction in subsequent operations. If it is less than 0.05 m, it is classified as a displacement within a small range and is usually considered within the normal tolerance. Thus, it can provide an exact starting value for the depth correction for the subsequent iteration of the Riccati equation.

[0208] Based on the depth adjustment value, first select the corresponding Riccati equation form and perform multiple iterations in combination with the drilling process data recorded previously. First, for each drilling stage, read the current depth information and add it to the previously obtained depth adjustment value Δz. If the accumulated depth value at this time is greater than the depth threshold proposed by the construction party during the preliminary exploration, such as 1000 meters, it is noted in the record that the planned drilling limit is approaching. Then, in the equation operation, refine the iteration step size of each parameter one by one. For example, subtract the previous round of error correction term from Δz in each round of iteration to obtain a new estimated value. The error correction term can be weighted according to the drill bit force, formation resistance, and vibration characteristic distribution, etc. If excessive vibration or obvious change in well inclination is found during some periods, modify the weighting coefficient to make the equation iteration converge more quickly. After each iteration, compare whether there is an obvious gap between the depth output by the current equation and the actual bottom hole position measured by the sensor. If the gap is greater than the preset limit, such as 0.1 meters, continue with the next iteration. If the deviation is within a small range for several consecutive iterations, the depth can be regarded as a relatively reliable real-time depth. After all acquisition cycles are completed, the real-time depth data will be centralized and made into an overall distribution corresponding to the time stamp. In the subsequent drilling monitoring stage, the depth data can be retrieved at any time for synchronous display or control scheduling, so that the depth information during the entire drilling process is in a dynamically updated state, and finally form real-time depth data.

Claims

1. Real-time depth measurement system for geological exploration drilling rigs, characterized in that: The system comprises: The vibration acquisition module arranges vibration sensors on the surface of the drilling tool to obtain the original vibration signal; divides the original vibration signal into 256 time series segments and sets the sampling interval value, calculates the amplitude and phase value for each time series segment, and generates the drilling dynamic characteristics; A signal conversion module calculates the amplitude difference and phase difference of adjacent time series segments in the drilling dynamic characteristics, compresses the time series segments according to the amplitude difference and phase difference to generate a compressed vibration sequence, and performs Fourier transformation on the compressed vibration sequence to generate characteristic spectrum parameters; A state correction module projects the characteristic spectrum parameters into the Grassmann space, calculates the projection eigenvector and eigenvalue, and generates projection matrix parameters; constructs a tangent plane projection function for the projection matrix parameters and calculates the geodesic distance to generate correction space coordinates; The depth measurement module calculates the position offset and angle offset of the drilling tool in three-dimensional space through the tangent plane projection function based on the corrected space coordinates, substitutes the position offset and angle offset into the Riccati equation for iterative calculation and outputs the depth value to generate real-time depth data.

2. The real-time depth measurement system for geological exploration drilling rigs according to claim 1 is characterized in that: The steps of obtaining the original vibration signal are: Arrange vibration sensors on the surface of the drilling tool to collect spatial position changes and frequency changes during drilling to obtain original vibration signals; Based on the original vibration signal, signal amplification, filtering and digitization are performed to obtain a preprocessed original vibration signal; Based on the pre-processed original vibration signal, the continuous vibration signal is converted into discrete data points, and the discrete data points are subjected to time series analysis to extract the dynamic characteristics of the drilling process to obtain the original vibration signal.

3. The real-time depth measurement system for geological exploration drilling rigs according to claim 1 is characterized in that: The steps for acquiring the drilling dynamic characteristics are as follows: The original vibration signal is divided into 256 time series segments, and a uniform sampling interval value is set for each segment to obtain the divided time series segments; Based on the segmented time series segments, the amplitude and phase of each time series segment are calculated to obtain the amplitude and phase characteristics of each time series segment; Based on the amplitude and phase characteristics of each time series segment, by analyzing the characteristics of all time series segments, statistical methods and pattern recognition are used to identify the dynamic change pattern during the drilling process and generate drilling dynamic characteristics.

4. The real-time depth measurement system for geological exploration drilling rigs according to claim 1 is characterized in that: The steps of obtaining the compressed vibration sequence are: Extracting adjacent time series segments from the drilling dynamic characteristics, calculating amplitude difference and phase difference for the series segments, capturing vibration changes during the drilling process by comparing the characteristics of each series segment, and obtaining amplitude difference and phase difference data; Based on the amplitude difference and phase difference data, data compression is performed on the time series segment, and adjacent sequence segments with similar characteristics are merged to obtain compressed time series segments; Based on the compressed time series segments, a continuous time series is re-integrated to obtain a compressed vibration sequence.

5. The real-time depth measurement system for geological exploration drilling rigs according to claim 1 is characterized in that: The steps of obtaining the characteristic spectrum parameters are as follows: Extracting data from the compressed vibration sequence to obtain data ready for Fourier transformation; Based on the data to be Fourier transformed, Fourier transform is performed to obtain spectrum data, and the formula is: Where F(k) is the kth frequency component in the frequency domain, x(j) is the jth data point in the time series, and N is the total number of data points; Based on the spectrum data, L1 norm constraint is applied to optimize the sparsity of the spectrum to obtain characteristic spectrum parameters.

6. The real-time depth measurement system for geological exploration drilling rigs according to claim 1 is characterized in that: The steps for obtaining the projection matrix parameters are: Based on the characteristic spectrum parameters, extract the frequency component matrix and the frequency characteristic data, perform decentralized processing on the frequency component matrix, and generate a decentralized frequency component matrix; Based on the decentralized frequency component matrix, the projection matrix is ​​calculated, and the expression is: Among them, M ij is the element in the i-th row and j-th column of the projection matrix, f i and f j are the amplitudes of the i-th and j-th frequency components in the frequency component matrix, u i and u j are the eigenvectors of the i-th and j-th frequency components, respectively, and λ i and λ j are the eigenvalues ​​of the i-th and j-th frequency components, θ i and θ j are the phase angles of the i-th and j-th frequency components, φ i and φ j are the amplitude angles of the i-th and j-th frequency components respectively; Based on the projection matrix, it is mapped to the Grassmann space, and geometric operations are performed through eigendecomposition to obtain eigenvectors and eigenvalues ​​to form projection matrix parameters.

7. The real-time depth measurement system for geological exploration drilling rigs according to claim 1 is characterized in that: The steps for obtaining the correction space coordinates are: Based on the projection matrix parameters, a cutting plane projection function is constructed, and the formula is: Where x and y represent the coordinate values ​​of any two points in the projection matrix, and F(x,y) is the tangent plane projection function; Use the tangent plane projection function to calculate the distance of all points in the projection matrix, and map the distance of each pair of points to a new distance value; The geodesic distance between the data points is calculated using the new distance value, and the spatial coordinates are corrected using the geodesic distance to obtain the corrected spatial coordinates.

8. The real-time depth measurement system for geological exploration drilling rigs according to claim 1 is characterized in that: The steps for acquiring the real-time depth data are as follows: Based on the corrected spatial coordinates, the position offset and angle offset of the drilling tool in three-dimensional space are calculated using a tangent plane projection function; Based on the position offset and angle offset, the depth adjustment value is calculated using the following formula: Among them, Δx and Δy represent the position offset of the drilling tool in the horizontal and vertical directions, α is the total angular offset of the drilling tool, and β is the angle with the vertical axis; Based on the depth adjustment value, the depth adjustment value is substituted into the Riccati equation, and iterative calculation is performed to adjust the depth value of the drilling tool to generate real-time depth data.

Images