Geological stratification parameter calibration method based on acoustic profile and core drilling test

Abstract

The invention relates to the technical field of geological exploration and geotechnical engineering testing, in particular to a geological stratification parameter calibration method based on an acoustic profile and a core drilling test, comprising the following steps: S1, acquiring acoustic profile data of a target area, and performing time-frequency domain preprocessing on the acoustic profile; s2, dividing a stratum interface, and generating a sound wave layered profile with a boundary mark; s3, core drilling sampling is carried out, and three-dimensional CT scanning and mechanical parameter testing are carried out on a core drilling sample; s4, establishing a sound wave speed-mechanical parameter mapping relation model, and correcting an interlayer parameter transition value; and S5, constructing a three-dimensional stratigraphic parameter matrix based on the corrected parameter transition values, and outputting a geological stratigraphic parameter calibration report with credibility indexes. According to the method, continuous calibration and credibility quantification of the stratum mechanical parameters are achieved by fusing the acoustic profile and the core drilling test data, and the precision and reliability of geological stratified modeling are effectively improved.

Description

Technical Field

[0001] The invention relates to the technical field of geological exploration and geotechnical engineering testing, and in particular to a geological stratification parameter calibration method based on acoustic wave profile and core drilling test. Background Art

[0002] In the field of geological exploration and engineering investigation, acoustic profiling technology is widely used in underground structure detection and stratum identification due to its high resolution, non-destructiveness and continuous acquisition capabilities; however, acoustic profiling data essentially reflects the physical response characteristics of the stratum medium and is difficult to directly reflect the specific lithological mechanical properties; in order to accurately grasp the key stratum mechanical parameters in engineering design, it is usually necessary to conduct sample measurements in combination with core tests; however, due to the limited number of core samples and their distribution being restricted by cost and construction conditions, it is impossible to fully cover all stratified areas, resulting in difficulty in continuous calibration of stratum mechanical parameters, which seriously affects the integrity and accuracy of geological modeling.

[0003] In the existing technology, there is a lack of a systematic method to align and fuse the acoustic profile data with the drill core parameters with high precision, especially in the parameter space extension and multi-source data consistency verification. At the same time, the traditional parameter interpolation method lacks physical constraints, easily introduces non-geologically reasonable transition zones, and fails to establish a clear mapping relationship from acoustic response to rock mechanical characteristics. Therefore, a geological stratification parameter calibration method based on acoustic profile and drill core test is urgently needed to solve the problem of discontinuity, unreliability and unquantifiable evaluation of geological stratification parameters. Summary of the invention

[0004] Based on the above purpose, the present invention provides a geological stratification parameter calibration method based on acoustic wave profile and core drilling test.

[0005] A method for calibrating geological layering parameters based on acoustic wave profile and core drilling test comprises the following steps: S1: Acquire the acoustic wave profile data of the target area, and perform time-frequency domain preprocessing on the acoustic wave profile to generate a preprocessed acoustic wave profile containing time-frequency features; S2: Based on the time-frequency characteristic differences of the preprocessed acoustic wave profile, the sliding window energy ratio algorithm is used to divide the stratum interface and generate an acoustic wave layered profile with boundary marks; S3: Drill core sampling is carried out at the corresponding positions of the key layers of the acoustic layered section, and three-dimensional CT scanning and mechanical parameter tests are performed on the drill core samples to generate a drill core parameter data set containing spatial coordinates; S4: Establish a model of acoustic wave velocity-mechanical parameter mapping relationship, match the time-frequency characteristics of the drill core parameter data set with the acoustic wave layered profile, and correct the inter-layer parameter transition value by using a robust estimation algorithm; S5: Construct a three-dimensional formation parameter matrix based on the corrected parameter transition values, verify the parameter continuity using the wave impedance inversion method, and output a geological stratification parameter calibration report with credibility indicators.

[0006] Optionally, the S1 specifically includes: S11: Use a towed multi-channel acoustic wave detection device to carry out detection operations along the straight detection profile set in the target area. The detection device is equipped with 12 receiving channels and 1 high-energy sound source. The center frequency of the sound source is set to 60kHz, and the detection depth coverage range is 0 to 50 meters. In an environment where the water depth changes no more than ±2 meters, it advances uniformly along the profile direction at a stable towing speed of 0.5m / s, emits a pulse signal every 0.2 seconds, and synchronously collects the reflection signals of each channel to continuously obtain the original acoustic wave profile data covering the target area; S12: De-noising the acquired original acoustic wave profile data, using a bandpass filter to remove signal interference with a frequency less than 10 Hz or greater than 500 Hz, and outputting the de-noised pure acoustic wave profile data; S13: Perform Fourier transform on the pure acoustic wave profile data to obtain frequency domain distribution characteristics; then apply short-time Fourier transform to the frequency domain characteristics, set the short-time window length to 50ms and the overlap rate to 50%, extract the spectrum characteristics of the profile signal changing with time, and form a time-frequency feature matrix; S14: The envelope of the time-frequency feature matrix is ​​extracted using Hilbert transform, and the extracted envelope signal is normalized to generate a preprocessed acoustic wave profile containing complete time-frequency features.

[0007] Optionally, the S14 specifically includes: S141: for each time domain signal sequence in the time-frequency feature matrix obtained in S13, perform Hilbert transform respectively to obtain an analytical expression of the corresponding signal; S142: Calculate the envelope amplitude of each time domain signal sequence, and combine the Hilbert transform result obtained in S141 with the original signal to form an analytical signal; S143: normalizing the envelope amplitude obtained in S142 by using an extreme value normalization method to obtain a standardized time-frequency feature data matrix in a value range of [0,1]; S144: All normalized envelope signals are combined according to the original channel order to form a preprocessed sound wave profile containing complete time-frequency features.

[0008] Optionally, the S2 specifically includes: S21: Based on the pre-processed sound wave profile, a sliding window with a length of 20 ms is set along the time axis direction of the sound wave propagation, and the sliding window is sequentially slid with a step length of 2 ms to intercept the signal data in each window; S22: Calculate the short-time energy value of the signal data in each sliding window, sum the square of the amplitude of each sampling point in the window, and obtain the energy value of the corresponding window; S23: sequentially calculating energy ratios between adjacent sliding windows; S24: Compare the energy ratio obtained in S23 with the preset formation interface determination threshold value 3.0. When the energy ratio exceeds the threshold value, the corresponding window boundary position is marked as the formation interface. The entire acoustic wave profile is scanned in sequence to obtain all formation interface positions that meet the conditions. S25: superimposing all the marked stratigraphic interface positions onto the pre-processed acoustic wave profile to generate an acoustic wave layered profile with boundary marks.

[0009] Optionally, the S3 specifically includes: S31: According to the acoustic layered profile obtained in S25, the layer with the highest energy ratio is selected as the key layer, and the spatial three-dimensional coordinate information of the key layer is recorded; S32: Arrange the drilling position at the spatial three-dimensional coordinate position recorded at the key layer, use the core drilling equipment to perform core sampling in the vertical direction, set the drilling speed to 0.2m / min, the core diameter to 50mm, and the target formation position to be drilled for 2m to obtain a columnar core sample; S33: performing an industrial-grade three-dimensional CT scan on the obtained columnar drill core sample, with the scanning parameters set to: scanning voltage 120 kV, scanning current 180 mA, and spatial scanning resolution 0.1 mm, to obtain three-dimensional structural image data inside the drill core sample; S34: The cylindrical drill core samples were cut into standard specimens with a diameter of 50 mm and a height of 100 mm. Uniaxial compressive strength tests were performed on the standard specimens at a constant loading rate of 0.5 MPa / s until the specimens were damaged. The compressive strength, elastic modulus and Poisson's ratio parameters of each specimen were measured. S35: The obtained spatial three-dimensional coordinates of the drill core sample, the structural image data obtained by CT scanning, and the mechanical parameter test data are integrated to generate a drill core parameter data set with spatial coordinate information.

[0010] Optionally, the S4 specifically includes: S41: extracting the sound wave propagation velocity data between the interfaces of various layers, obtaining the average sound wave propagation velocity in each layer segment, and forming a velocity vector set; S42: extracting the drill core sample parameters corresponding to the spatial position of the layered segment in S35, obtaining the compressive strength, elastic modulus and Poisson's ratio corresponding to each layer segment, and forming a mechanical parameter vector set; S43: Based on the correspondence between the velocity vector set and the mechanical parameter vector set, a mapping relationship model between the acoustic wave velocity and the mechanical parameters is constructed using a multivariate linear regression method, and a predicted mechanical parameter value corresponding to each acoustic wave velocity value is output; S44: applying the established mapping relationship model to all layered segments in the acoustic layered section, combining the average acoustic wave velocity thereof, calculating the initial mechanical parameter transition data, and performing spatial registration with the corresponding positions of the drill core data set; S45: For the deviations in the initial parameter transition data caused by local anomalies or measurement errors, a robust estimation algorithm is used to perform parameter correction to form a corrected mechanical parameter transition value.

[0011] Optionally, the S43 specifically includes: S431: Using the multivariate linear regression method, three linear mapping models were established with compressive strength, elastic modulus, and Poisson's ratio as dependent variables and acoustic wave velocity as independent variable. The regression equations are expressed as follows: ; ; ; In the formula, represents the sound wave propagation speed in the i-th layer; , , They represent the compressive strength, elastic modulus and Poisson’s ratio corresponding to the i-th layer segment respectively; , , is the slope coefficient of each regression model; , , is the intercept constant term of each regression model; S432: Use the least squares method to fit the regression coefficients, minimize the square error between the predicted value and the actual mechanical parameter value, and output the predicted mechanical parameter triple corresponding to each acoustic wave velocity value .

[0012] Optionally, the S44 specifically includes: S441: applying the established acoustic wave velocity and mechanical parameter mapping relationship model to the average acoustic wave velocity value of each layered segment in the acoustic wave layered section, substituting each acoustic wave velocity value into the prediction model to obtain the corresponding initial mechanical parameters, including compressive strength, elastic modulus and Poisson's ratio; S442: Associating each set of calculated initial mechanical parameters with the spatial coordinates in the acoustic wave layered section to form a preliminary parameter transition data set; S443: extracting the spatial position of each sample in the drill core parameter data set, and judging whether the registration condition is met according to the spatial distance between the center position of the acoustic layer segment and the center position of the drill core sample; when the spatial distance is within the preset error range, it is considered that the acoustic wave prediction parameter and the drill core measured parameter constitute a valid matching pair; S444: Output all registered predicted parameter and measured parameter pairs.

[0013] Optionally, the S45 specifically includes: S451: for each set of predicted parameters and measured parameters that have been aligned in S444, calculating the difference between the two; S452: judging each deviation value according to the set residual threshold; when the deviation amplitude is within the allowable range, retaining the deviation; when it exceeds the range, adjusting the deviation to the threshold value and retaining its positive and negative signs; S453: subtracting the adjusted deviation from the original prediction parameter to obtain a corrected mechanical parameter transition value; S454: Summarize all corrected mechanical parameter transition values ​​in stratigraphic order to form a complete transition parameter sequence.

[0014] Optionally, the S5 specifically includes: S51: Using the corrected mechanical parameter transition value as input, the three-dimensional space covered by the acoustic layered section is discretized into an equilateral voxel grid, and the compressive strength, elastic modulus and Poisson's ratio are filled in each voxel according to the spatial coordinates of the voxel center point using the layered Kriging interpolation method to construct a three-dimensional formation parameter matrix with a resolution of 0.25m×0.25m×0.25m; S52: For each voxel in the three-dimensional stratigraphic parameter matrix, the wave impedance value is calculated according to the acoustic wave propagation velocity and rock density constant corresponding to the voxel; then, the synthetic acoustic wave record is generated by forward convolution modeling, and compared with the original acoustic wave profile one by one, and the residual energy distribution between the two is calculated; S53: Apply the Gauss-Newton iterative inversion method to minimize the residual energy, and iteratively update the voxel wave impedance value and the gradient between adjacent voxels until the residual energy converges to below the preset threshold of 2%; S54: Calculate the credibility index of each layer segment according to the final residual energy and the inter-voxel gradient smoothness, and classify the credibility into three levels: A, B, and C; S55: Generate a geological stratification parameter calibration report, which includes three-dimensional stratigraphic parameter matrix visualization slices and the credibility level of each layer segment, and is output and archived in an electronically signed manner.

[0015] Beneficial effects of the present invention: The present invention achieves high-precision prediction and spatial extension of formation parameters by establishing a multivariate linear mapping relationship model between acoustic wave propagation velocity and rock mechanical parameters, and combines the measured data of core samples for alignment and correction, thereby effectively solving the parameter missing problem caused by sparse core data. At the same time, the formation interface is accurately identified through a sliding window energy ratio algorithm, and the layered identification capability of the acoustic wave profile is improved by combining Hilbert envelope extraction and normalization processing.

[0016] The present invention introduces robust estimation and wave impedance inversion algorithms to eliminate anomalies and optimize continuity of initial prediction parameters, and quantitatively evaluates the calibration results through a credibility grading mechanism, thereby enhancing the engineering usability and interpretability of the parameter matrix. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] In order to more clearly illustrate the technical solutions in the present invention or the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings in the following description are only for the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0018] Figure 1 A schematic diagram of a geological layering parameter calibration method according to an embodiment of the present invention; Figure 2 It is a schematic diagram of a mapping relationship model construction and transition value correction method according to an embodiment of the present invention. DETAILED DESCRIPTION

[0019] The present invention is described in detail below in conjunction with the accompanying drawings and specific embodiments. At the same time, it is explained here that in order to make the embodiments more detailed, the following embodiments are the best and preferred embodiments, and those skilled in the art may also adopt other alternatives to implement some known technologies; and the accompanying drawings are only for more specific description of the embodiments, and are not intended to specifically limit the present invention.

[0020] like Figure 1-Figure 2 As shown, a geological layer parameter calibration method based on acoustic profile and core test includes the following steps: S1: Acquire the acoustic wave profile data of the target area, and perform time-frequency domain preprocessing on the acoustic wave profile to generate a preprocessed acoustic wave profile containing time-frequency features; S2: Based on the time-frequency characteristic differences of the preprocessed acoustic wave profile, the sliding window energy ratio algorithm is used to divide the stratum interface and generate an acoustic wave layered profile with boundary marks; S3: Drill core sampling is carried out at the corresponding positions of the key layers of the acoustic layered section, and three-dimensional CT scanning and mechanical parameter tests are performed on the drill core samples to generate a drill core parameter data set containing spatial coordinates; S4: Establish a model of acoustic wave velocity-mechanical parameter mapping relationship, match the time-frequency characteristics of the drill core parameter data set with the acoustic wave layered profile, and correct the inter-layer parameter transition value by using a robust estimation algorithm; S5: Construct a three-dimensional formation parameter matrix based on the corrected parameter transition values, verify the parameter continuity using the wave impedance inversion method, and output a geological stratification parameter calibration report with credibility indicators.

[0021] S1 specifically includes: S11: Use a towed multi-channel acoustic wave detection device to carry out detection operations along the straight detection profile set in the target area. The detection device is equipped with 12 receiving channels and 1 high-energy sound source. The center frequency of the sound source is set to 60kHz, and the detection depth coverage range is 0 to 50 meters. In an environment where the water depth changes no more than ±2 meters, it advances uniformly along the profile direction at a stable towing speed of 0.5m / s, emits a pulse signal every 0.2 seconds, and synchronously collects the reflection signals of each channel to continuously obtain the original acoustic wave profile data covering the target area; S12: De-noising the acquired original acoustic wave profile data, using a bandpass filter to remove signal interference with a frequency less than 10 Hz or greater than 500 Hz, and outputting the de-noised pure acoustic wave profile data; S13: Perform Fourier transform on the pure acoustic wave profile data to obtain frequency domain distribution characteristics; then apply short-time Fourier transform to the frequency domain characteristics, set the short-time window length to 50ms and the overlap rate to 50%, extract the spectrum characteristics of the profile signal changing with time, and form a time-frequency feature matrix; S14: Use Hilbert transform to extract the envelope of the time-frequency feature matrix, and normalize the extracted envelope signal to generate a preprocessed acoustic wave profile containing complete time-frequency features; through the specific implementation of the above steps S11 to S14, the interference factors are effectively removed, the time-frequency domain characteristics of the acoustic wave data are clarified, and the accuracy and reliability of the subsequent stratigraphic division are ensured.

[0022] S14 specifically includes: S141: For each time domain signal sequence in the time-frequency feature matrix obtained in S13, perform Hilbert transform respectively to obtain an analytical expression of the corresponding signal, where the Hilbert transform calculation formula is: , where represents the imaginary part of the analytical signal obtained by Hilbert transform; Indicates that the time domain signal to be processed is in the integral variable The signal amplitude at ; t represents the target time point for solving the Hilbert transform result; is the integral variable, which indicates the integration time of the signal; is the integration variable The differential of ; PV is the Cauchy principal value integral operator; S142: Calculate the envelope amplitude of each time domain signal sequence, combine the Hilbert transform result obtained in S141 with the original signal to form an analytical signal. The specific calculation formula is as follows: , where Represents the instantaneous envelope amplitude of the time domain signal; is the original amplitude of the time domain signal at time t; is the imaginary part of the analytical signal obtained by Hilbert transform; S143: The envelope amplitude obtained in S142 is normalized by using the extreme value normalization method to obtain a standardized time-frequency feature data matrix in the range of [0,1]. The specific calculation formula is as follows: , where is the normalized envelope amplitude; is the maximum value of the envelope amplitude during the entire analysis period; is the minimum value of the envelope amplitude in the entire analysis period; S144: All normalized envelope signals are combined in the original channel order to form a preprocessed acoustic wave profile containing complete time-frequency characteristics; through the implementation of the above steps S141 to S144, the envelope characteristics extracted by the Hilbert transform make the stratigraphic structure characteristics hidden in the acoustic wave signal more prominent, and at the same time, the normalization processing eliminates the interference of the signal amplitude difference between different measurement channels, thereby improving the accuracy and reliability of subsequent geological boundary identification.

[0023] S2 specifically includes: S21: Based on the pre-processed sound wave profile, a sliding window with a length of 20 ms is set along the time axis direction of the sound wave propagation, and the sliding window is sequentially slid with a step length of 2 ms to intercept the signal data in each window; S22: Calculate the short-time energy value of the signal data in each sliding window, sum the square of the amplitude of each sampling point in the window, and obtain the energy value of the corresponding window. The calculation formula is as follows: , where Indicates the energy value of the signal in the current sliding window; represents the signal amplitude of the i-th sampling point in the current sliding window; N represents the number of sampling points in the current sliding window; S23: Calculate the energy ratios between adjacent sliding windows in sequence, and use the energy value of the window with larger energy divided by the energy value of the window with smaller energy as the energy ratio. The calculation formula is as follows: , where: Represents the energy ratio between the jth window and the j+1th window; Represents the energy value in the jth sliding window; Represents the energy value in the j+1th sliding window; S24: Compare the energy ratio obtained in S23 with the preset formation interface determination threshold value 3.0. When the energy ratio exceeds the threshold value, the corresponding window boundary position is marked as the formation interface. The entire acoustic wave profile is scanned in sequence to obtain all formation interface positions that meet the conditions. S25: superimpose all marked stratigraphic interface positions onto the preprocessed acoustic profile to generate an acoustic layered profile with boundary marks; through the implementation of the above steps S21 to S25, the sliding window energy ratio algorithm is used to accurately identify the positions with significant energy differences between adjacent strata, effectively improving the recognition accuracy of stratigraphic interfaces, thereby providing a clear spatial reference basis for subsequent core sampling and formation parameter calibration.

[0024] S3 specifically includes: S31: According to the acoustic layered profile obtained in S25, a layer with the highest energy ratio and the most obvious continuity is selected as a key layer, and the spatial three-dimensional coordinate information of the key layer is recorded; S32: Arrange the drilling position at the spatial three-dimensional coordinate position recorded at the key layer, use the core drilling equipment to perform core sampling in the vertical direction, set the drilling speed to 0.2m / min, the core diameter to 50mm, and the target formation position to be drilled for 2m to obtain a columnar core sample; S33: performing an industrial-grade three-dimensional CT scan on the obtained columnar drill core sample, with the scanning parameters set to: scanning voltage 120 kV, scanning current 180 mA, and spatial scanning resolution 0.1 mm, to obtain three-dimensional structural image data inside the drill core sample; S34: The cylindrical drill core samples were cut into standard specimens with a diameter of 50 mm and a height of 100 mm. Uniaxial compressive strength tests were performed on the standard specimens at a constant loading rate of 0.5 MPa / s until the specimens were damaged. The compressive strength, elastic modulus and Poisson's ratio parameters of each specimen were measured. S35: The spatial three-dimensional coordinates of the obtained core samples, the structural image data obtained by CT scanning, and the mechanical parameter test data are integrated to generate a core parameter data set with spatial coordinate information; through the implementation of the above steps S31 to S35, the core sampling positions corresponding to the key layers are clarified, the internal structure of the formation and the detailed mechanical performance parameters are accurately obtained, and the precise spatial alignment of the core parameters and the acoustic layering profile is achieved, providing high-precision basic data guarantee for the subsequent geological parameter calibration.

[0025] S4 specifically includes: S41: extracting the sound wave propagation velocity data between the interfaces of various layers, obtaining the average sound wave propagation velocity in each layer segment, and forming a velocity vector set; S42: extracting the drill core sample parameters corresponding to the spatial position of the layered segment in S35, obtaining the compressive strength, elastic modulus and Poisson's ratio corresponding to each layer segment, and forming a mechanical parameter vector set; S43: Based on the correspondence between the velocity vector set and the mechanical parameter vector set, a mapping relationship model between the acoustic wave velocity and the mechanical parameters is constructed using a multivariate linear regression method, and a predicted mechanical parameter value corresponding to each acoustic wave velocity value is output; S44: applying the established mapping relationship model to all layered segments in the acoustic layered section, combining the average acoustic wave velocity thereof, calculating the initial mechanical parameter transition data, and performing spatial registration with the corresponding positions of the drill core data set; S45: In view of the deviations in the initial parameter transition data caused by local anomalies or measurement errors, a robust estimation algorithm is used to correct the parameters to form a corrected mechanical parameter transition value; through the implementation of the above steps S41 to S45, a stable parameter mapping relationship is constructed by combining the velocity information with the mechanical characteristics of the drill core, and the robust estimation is introduced to improve the robustness of the model to abnormal data, laying an accurate data foundation for the construction of the three-dimensional formation parameter matrix.

[0026] S43 specifically includes: S431: Using the multivariate linear regression method, three linear mapping models were established with compressive strength, elastic modulus, and Poisson's ratio as dependent variables and acoustic wave velocity as independent variable. The regression equations are expressed as follows: ; ; ; In the formula, represents the sound wave propagation speed in the i-th layer; , , They represent the compressive strength, elastic modulus and Poisson’s ratio corresponding to the i-th layer segment respectively; , , is the slope coefficient of each regression model; , , is the intercept constant term of each regression model; S432: Use the least squares method to fit the regression coefficients, minimize the square error between the predicted value and the actual mechanical parameter value, and output the predicted mechanical parameter triple corresponding to each acoustic wave velocity value , which is used for the subsequent parameter transition construction; through the implementation of the above steps S431 to S432, a linear response relationship between the acoustic wave velocity and the key mechanical parameters is established, and a stable and interpretable geomechanical parameter prediction model is constructed, which significantly improves the parameter compensation capability of the acoustic profile in the coreless section area, and provides a complete prediction basis for anti-error correction and three-dimensional parameter matrix construction.

[0027] S44 specifically includes: S441: The established acoustic wave velocity and mechanical parameter mapping model is applied to the average acoustic wave velocity value of each layer segment in the acoustic wave layered section. Substitute into the prediction model to obtain the corresponding initial mechanical parameters, including compressive strength , elastic modulus and Poisson's ratio , and its calculation formula is as follows: ; ; ; In the formula, is the average sound wave propagation velocity in the i-th layer; , , are the predicted compressive strength, elastic modulus and Poisson’s ratio, respectively; , , is the slope coefficient in the mapping model; , , is the intercept constant term; S442: Associating each set of calculated initial mechanical parameters with the spatial coordinates in the acoustic wave layered section to form a preliminary parameter transition data set; S443: extracting the spatial position of each sample in the drill core parameter data set, and judging whether the registration condition is met according to the spatial distance between the center position of the acoustic layer segment and the center position of the drill core sample; when the spatial distance is within the preset error range, it is considered that the acoustic wave prediction parameter and the drill core measured parameter constitute a valid matching pair; Specifically, the Euclidean distance matching method is used, and the coordinates of the center of gravity of the acoustic wave segment are set as , the center coordinates of the drill core sample are , when the following distance constraints are met: , then the acoustic wave prediction parameters are spatially aligned with the drill core sample; is the spatial distance between the acoustic layer and the core sample, is the registration distance threshold (set to 0.5m); S444: Output all the registered predicted parameter and measured parameter pairs, which are used as the input basis for the subsequent robust estimation process to ensure that the predicted data is highly consistent with the actual drill core position; through the implementation of the above steps S441 to S444, combined with the average acoustic wave velocity and the established mapping relationship, the prediction generation of mechanical parameters is realized, and it is aligned with the measured drill core data with high precision through spatial registration technology, providing a unified and reliable input basis for parameter anomaly identification and model correction.

[0028] S45 specifically includes: S451: for each set of registered predicted parameters and measured parameters in S444, the difference between the two is calculated to measure the model prediction deviation; S452: judging each deviation value according to the set residual threshold; when the deviation amplitude is within the allowable range, retaining the deviation; when it exceeds the range, adjusting the deviation to the threshold value and retaining its positive and negative signs to limit the influence of abnormal values; S453: deducting the adjusted deviation from the original predicted parameter to obtain a corrected transition value of the mechanical parameter, thereby ensuring that the parameter of each layer section is closer to the measured value; S454: All corrected transition values ​​of mechanical parameters are summarized in stratigraphic order to form a complete transition parameter sequence, which is directly used for the subsequent construction of the three-dimensional parameter matrix. Through the above steps, the abnormal fluctuations in the initial prediction parameters are effectively suppressed, and high-precision alignment of the model output and the measured data is achieved, and the stability and reliability of the formation parameter calibration are significantly improved.

[0029] The steps for parameter correction using the robust estimation algorithm are as follows: First, for each set of registered predicted mechanical parameters and measured parameters, the residual value is calculated and the predicted parameters are expressed as , the measured parameters are expressed as , then the i-th residual It is calculated as the difference between the predicted value and the measured value, that is, the formula is: ; Then, for each residual , use the Huber influence function to calculate the correction component , where the threshold is set , when the absolute value of the residual does not exceed the threshold, take ; When the absolute value of the residual is greater than the threshold, take ; The specific expression is: ; Where: is the correction component of the i-th residual; is the preset residual threshold; for The symbolic function of Next, according to the Huber correction component Update the predicted parameters to form the i-th corrected mechanical parameter transition value , which is calculated as the original predicted value minus the correction component, that is, ; Finally, all The layers are summarized in order to form a complete sequence of corrected mechanical parameter transition values.

[0030] S5 specifically includes: S51: Using the corrected mechanical parameter transition value as input, the three-dimensional space covered by the acoustic layered section is discretized into an equilateral voxel grid, and the compressive strength, elastic modulus and Poisson's ratio are filled in each voxel according to the spatial coordinates of the voxel center point using the layered Kriging interpolation method to construct a three-dimensional formation parameter matrix with a resolution of 0.25m×0.25m×0.25m; S52: For each voxel in the three-dimensional stratigraphic parameter matrix, the wave impedance value is calculated according to the acoustic wave propagation velocity and rock density constant corresponding to the voxel; then, the synthetic acoustic wave record is generated by forward convolution modeling, and compared with the original acoustic wave profile one by one, and the residual energy distribution between the two is calculated; The specific steps of the above S52 are as follows: S521: For each voxel in the three-dimensional formation parameter matrix, extract the sound wave propagation velocity value at its center position and the set rock density constant, and calculate the corresponding wave impedance value. The formula is as follows: In the formula, represents the wave impedance value of the kth voxel; represents the rock density constant of the kth voxel (in units of ); represents the longitudinal wave propagation velocity of the k-th voxel (in m / s); S522: According to the wave impedance sequence Arrange along the detection path direction and calculate the reflection coefficient sequence between adjacent voxels , the reflection coefficient is calculated as follows: , where represents the reflection coefficient on the kth interface; , are the wave impedance values ​​on the upper and lower sides of the interface respectively; S523: Perform one-dimensional convolution on the reflection coefficient sequence and the standard wavelet signal to generate a synthetic sound wave record on the path. The calculation expression is as follows: ; In the formula, represents a synthetic sound wave record; is the discrete form of the reflection coefficient sequence; is the standard wavelet used in the acquisition system; Represents the convolution operation; S524: For each detection path node, the generated synthetic sound wave record is compared with the measured original sound wave profile data at the same position to calculate the residual signal , and calculate the residual energy E in the form of time-domain square error, the expression is as follows: , where E is the residual energy of the path node; To synthesize the sound wave signal; is the original measured sound wave data; t is the sampling time point.

[0031] S53: Apply the Gauss-Newton iterative inversion method to minimize the residual energy, and iteratively update the voxel wave impedance value and the gradient between adjacent voxels until the residual energy converges to below the preset threshold of 2%; The specific steps of S53 are as follows: S531: Define the objective function. For each path node, set the objective to minimize the residual energy between the synthetic sound wave record and the measured sound wave record. Let the synthetic sound wave be , the measured sound wave is , the wave impedance parameter vector is , then the residual vector is: ; The objective function is defined as the sum of the squares of the residual vector: ; S532: Derivative the wave impedance vector Z and construct the Jacobian matrix of the residual function with respect to the parameters ,Right now: ,in, represents the partial derivative of the i-th residual with respect to the j-th wave impedance parameter; S533: Use Gauss-Newton iteration method to minimize optimization. Let the current iteration step be k, then the parameter update amount Solve the following linear equations: ; The updated wave impedance parameters are: ; S534: In the iterative process, in order to avoid sudden changes in the wave impedance values ​​of adjacent voxels, a gradient constraint term is added to penalize the square of the impedance difference between adjacent voxels, so that the solution has spatial smoothness, forming an improved objective function: ,in, is the smoothing coefficient (value ranges from 0.1 to 1.0). This regularization term is included in the Jacobian matrix solution in each iteration; S535: Repeat the iteration step, if the residual energy change rate between two consecutive iterations satisfies: , that is, the residual energy change rate is less than 2%, it is judged to have converged, and the current wave impedance distribution is output as the final inversion result.

[0032] S54: The credibility index of each layer segment is calculated based on the final residual energy and the gradient smoothness between voxels, and the credibility is graded into three levels: A, B, and C. Level A corresponds to residual energy less than 1% and gradient smoothness greater than 95%; Level B corresponds to residual energy between 1% and 2% and gradient smoothness between 90% and 95%; Level C corresponds to residual energy greater than 2% or gradient smoothness less than 90%; S55: Generate a geological stratification parameter calibration report, the report content includes a three-dimensional stratigraphic parameter matrix visualization slice and each layer segment credibility level, and output it for archiving in the form of an electronic signature; through the implementation of steps S51 to S55, high-precision interpolation of correction parameters in three-dimensional space is achieved, the continuity of parameters is strictly checked using wave impedance inversion, and the results are output in the form of quantitative credibility, thereby ensuring the integrity and reliability of geological stratification parameter calibration.

[0033] The present invention covers any substitution, modification, equivalent method and scheme made on the essence and scope of the present invention. In order to make the public have a thorough understanding of the present invention, specific details are described in detail in the following preferred embodiments of the present invention, but those skilled in the art can fully understand the present invention without the description of these details. In addition, in order to avoid unnecessary confusion about the essence of the present invention, well-known methods, processes, procedures, components and circuits are not described in detail.

[0034] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principle of the present invention. These improvements and modifications should also be regarded as the scope of protection of the present invention.

Claims

1. A method for calibrating geological layering parameters based on acoustic profile and core drilling test, characterized in that: The following steps are involved: S1: Acquire the acoustic wave profile data of the target area, and perform time-frequency domain preprocessing on the acoustic wave profile to generate a preprocessed acoustic wave profile containing time-frequency features; S2: Based on the time-frequency characteristic differences of the preprocessed acoustic wave profile, the sliding window energy ratio algorithm is used to divide the stratum interface and generate an acoustic wave layered profile with boundary marks; S3: Drill core sampling is carried out at the corresponding positions of the key layers of the acoustic layered section, and three-dimensional CT scanning and mechanical parameter tests are performed on the drill core samples to generate a drill core parameter data set containing spatial coordinates; S4: Establish a model of acoustic wave velocity-mechanical parameter mapping relationship, match the time-frequency characteristics of the drill core parameter data set with the acoustic wave layered profile, and correct the inter-layer parameter transition value by using a robust estimation algorithm; S5: Construct a three-dimensional formation parameter matrix based on the corrected parameter transition values, verify the parameter continuity using the wave impedance inversion method, and output a geological stratification parameter calibration report with credibility indicators.

2. The method for calibrating geological layering parameters based on acoustic wave profile and core drilling test according to claim 1 is characterized in that: The S1 specifically includes: S11: Use a towed multi-channel acoustic wave detection device to carry out detection operations along the straight detection profile set in the target area. The detection device is equipped with 12 receiving channels and 1 high-energy sound source. The center frequency of the sound source is set to 60kHz, and the detection depth coverage range is 0 to 50 meters. In an environment where the water depth changes no more than ±2 meters, it advances uniformly along the profile direction at a stable towing speed of 0.5m / s, emits a pulse signal every 0.2 seconds, and synchronously collects the reflection signals of each channel to continuously obtain the original acoustic wave profile data covering the target area; S12: De-noising the acquired original acoustic wave profile data, using a bandpass filter to remove signal interference with a frequency less than 10 Hz or greater than 500 Hz, and outputting the de-noised pure acoustic wave profile data; S13: Perform Fourier transform on the pure acoustic wave profile data to obtain frequency domain distribution characteristics; then apply short-time Fourier transform to the frequency domain characteristics, set the short-time window length to 50ms and the overlap rate to 50%, extract the spectrum characteristics of the profile signal changing with time, and form a time-frequency feature matrix; S14: The envelope of the time-frequency feature matrix is ​​extracted using Hilbert transform, and the extracted envelope signal is normalized to generate a preprocessed acoustic wave profile containing complete time-frequency features.

3. The method for calibrating geological stratification parameters based on acoustic profile and core drilling test according to claim 2 is characterized in that: The S14 specifically includes: S141: for each time domain signal sequence in the time-frequency feature matrix obtained in S13, perform Hilbert transform respectively to obtain an analytical expression of the corresponding signal; S142: Calculate the envelope amplitude of each time domain signal sequence, and combine the Hilbert transform result obtained in S141 with the original signal to form an analytical signal; S143: normalizing the envelope amplitude obtained in S142 by using an extreme value normalization method to obtain a standardized time-frequency feature data matrix in a value range of [0,1]; S144: All normalized envelope signals are combined according to the original channel order to form a preprocessed sound wave profile containing complete time-frequency features.

4. The method for calibrating geological layering parameters based on acoustic profile and core drilling test according to claim 1 is characterized in that: The S2 specifically includes: S21: Based on the pre-processed sound wave profile, a sliding window with a length of 20 ms is set along the time axis direction of the sound wave propagation, and the sliding window is sequentially slid with a step length of 2 ms to intercept the signal data in each window; S22: Calculate the short-time energy value of the signal data in each sliding window, sum the square of the amplitude of each sampling point in the window, and obtain the energy value of the corresponding window; S23: sequentially calculating energy ratios between adjacent sliding windows; S24: Compare the energy ratio obtained in S23 with the preset formation interface determination threshold value 3.

0. When the energy ratio exceeds the threshold value, the corresponding window boundary position is marked as the formation interface. The entire acoustic wave profile is scanned in sequence to obtain all formation interface positions that meet the conditions. S25: superimposing all the marked stratigraphic interface positions onto the pre-processed acoustic wave profile to generate an acoustic wave layered profile with boundary marks.

5. The method for calibrating geological layering parameters based on acoustic wave profile and core drilling test according to claim 4 is characterized in that: The S3 specifically includes: S31: According to the acoustic layered profile obtained in S25, the layer with the highest energy ratio is selected as the key layer, and the spatial three-dimensional coordinate information of the key layer is recorded; S32: Arrange the drilling position at the spatial three-dimensional coordinate position recorded at the key layer, use the core drilling equipment to perform core sampling in the vertical direction, set the drilling speed to 0.2m / min, the core diameter to 50mm, and the target formation position to be drilled for 2m to obtain a columnar core sample; S33: performing an industrial-grade three-dimensional CT scan on the obtained columnar drill core sample, with the scanning parameters set to: scanning voltage 120 kV, scanning current 180 mA, and spatial scanning resolution 0.1 mm, to obtain three-dimensional structural image data inside the drill core sample; S34: The cylindrical drill core samples were cut into standard specimens with a diameter of 50 mm and a height of 100 mm. Uniaxial compressive strength tests were performed on the standard specimens at a constant loading rate of 0.5 MPa / s until the specimens were damaged. The compressive strength, elastic modulus and Poisson's ratio parameters of each specimen were measured. S35: The obtained spatial three-dimensional coordinates of the drill core sample, the structural image data obtained by CT scanning, and the mechanical parameter test data are integrated to generate a drill core parameter data set with spatial coordinate information.

6. The method for calibrating geological layering parameters based on acoustic profile and core drilling test according to claim 1 is characterized in that: The S4 specifically includes: S41: extracting the sound wave propagation velocity data between the interfaces of various layers, obtaining the average sound wave propagation velocity in each layer segment, and forming a velocity vector set; S42: extracting the drill core sample parameters corresponding to the spatial position of the layered segment in S35, obtaining the compressive strength, elastic modulus and Poisson's ratio corresponding to each layer segment, and forming a mechanical parameter vector set; S43: Based on the correspondence between the velocity vector set and the mechanical parameter vector set, a mapping relationship model between the acoustic wave velocity and the mechanical parameters is constructed using a multivariate linear regression method, and a predicted mechanical parameter value corresponding to each acoustic wave velocity value is output; S44: applying the established mapping relationship model to all layered segments in the acoustic layered section, combining the average acoustic wave velocity thereof, calculating the initial mechanical parameter transition data, and performing spatial registration with the corresponding positions of the drill core data set; S45: For the deviations in the initial parameter transition data caused by local anomalies or measurement errors, a robust estimation algorithm is used to perform parameter correction to form a corrected mechanical parameter transition value.

7. A method for calibrating geological layering parameters based on acoustic profile and core drilling test according to claim 6, characterized in that: The S43 specifically includes: S431: Using the multivariate linear regression method, three linear mapping models were established with compressive strength, elastic modulus, and Poisson's ratio as dependent variables and acoustic wave velocity as independent variable. The regression equations are expressed as follows: ; ; ; In the formula, represents the sound wave propagation speed in the i-th layer; , , They represent the compressive strength, elastic modulus and Poisson’s ratio corresponding to the i-th layer segment respectively; , , is the slope coefficient of each regression model; , , is the intercept constant term of each regression model; S432: Use the least squares method to fit the regression coefficients, minimize the square error between the predicted value and the actual mechanical parameter value, and output the predicted mechanical parameter triple corresponding to each acoustic wave velocity value .

8. The method for calibrating geological layering parameters based on acoustic profile and core drilling test according to claim 7 is characterized in that: The S44 specifically includes: S441: applying the established acoustic wave velocity and mechanical parameter mapping relationship model to the average acoustic wave velocity value of each layered segment in the acoustic wave layered section, substituting each acoustic wave velocity value into the prediction model to obtain the corresponding initial mechanical parameters, including compressive strength, elastic modulus and Poisson's ratio; S442: Associating each set of calculated initial mechanical parameters with the spatial coordinates in the acoustic wave layered section to form a preliminary parameter transition data set; S443: extracting the spatial position of each sample in the drill core parameter data set, and judging whether the registration condition is met according to the spatial distance between the center position of the acoustic layer segment and the center position of the drill core sample; when the spatial distance is within the preset error range, it is considered that the acoustic wave prediction parameter and the drill core measured parameter constitute a valid matching pair; S444: Output all registered predicted parameter and measured parameter pairs.

9. A method for calibrating geological layering parameters based on acoustic profile and core drilling test according to claim 8, characterized in that: The S45 specifically includes: S451: for each set of predicted parameters and measured parameters that have been aligned in S444, calculating the difference between the two; S452: judging each deviation value according to the set residual threshold; when the deviation amplitude is within the allowable range, retaining the deviation; when it exceeds the range, adjusting the deviation to the threshold value and retaining its positive and negative signs; S453: subtracting the adjusted deviation from the original prediction parameter to obtain a corrected mechanical parameter transition value; S454: Summarize all corrected mechanical parameter transition values ​​in stratigraphic order to form a complete transition parameter sequence.

10. The method for calibrating geological layering parameters based on acoustic wave profile and core drilling test according to claim 1, characterized in that: The S5 specifically includes: S51: Using the corrected mechanical parameter transition value as input, the three-dimensional space covered by the acoustic layered section is discretized into an equilateral voxel grid, and the compressive strength, elastic modulus and Poisson's ratio are filled in each voxel according to the spatial coordinates of the voxel center point using the layered Kriging interpolation method to construct a three-dimensional formation parameter matrix with a resolution of 0.25m×0.25m×0.25m; S52: For each voxel in the three-dimensional stratigraphic parameter matrix, the wave impedance value is calculated according to the acoustic wave propagation velocity and rock density constant corresponding to the voxel; then, the synthetic acoustic wave record is generated by forward convolution modeling, and compared with the original acoustic wave profile one by one, and the residual energy distribution between the two is calculated; S53: Apply the Gauss-Newton iterative inversion method to minimize the residual energy, and iteratively update the voxel wave impedance value and the gradient between adjacent voxels until the residual energy converges to below the preset threshold of 2%; S54: Calculate the credibility index of each layer segment according to the final residual energy and the inter-voxel gradient smoothness, and classify the credibility into three levels: A, B, and C; S55: Generate a geological stratification parameter calibration report, which includes three-dimensional stratigraphic parameter matrix visualization slices and the credibility level of each layer segment, and is output and archived in an electronically signed manner.

Images