Subsea tunnel overlying strata fracture three-dimensional identification method

Abstract

The invention relates to the technical field of ocean engineering and structural health monitoring, in particular to a three-dimensional recognition method for overlying strata fractures of a subsea tunnel, and the method comprises the steps: when a propagation path between a sound wave signal transmitter and a sound wave signal receiver does not have a fracture, receiving signals are stable in performance; and when cracks exist in the path, the received signal is obviously decayed or suddenly dropped. According to the method, online, rapid and accurate detection of overlying strata fractures of the subsea tunnel is realized through a multi-source space superposition recognition algorithm based on sound wave decay characteristics, and a single-path measurement error can be eliminated and the positioning accuracy can be improved by utilizing a multi-path superposition principle; the automation degree of the system is high, signals can be rapidly collected and processed, and the on-line monitoring requirement in the construction process is met; equipment installation is flexible, the emitter only needs to be arranged outside the overlying strata, and transformation of the internal structure of the tunnel is reduced.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of ocean engineering and structural health monitoring, in particular to a three-dimensional identification method for overburden fissures of a submarine tunnel. BACKGROUND

[0002] With the increasing scale and depth of submarine tunnel construction, the geological structure in the overburden layer is complex and variable, and the fissure development degree has an important influence on the overall stability and safety of the tunnel. Specifically, the overburden fissures may weaken the self-stability of the surrounding rock, leading to uneven settlement of the foundation or deformation of the surrounding rock, and further causing uneven pressure on the tunnel lining, which may trigger the risk of cracks, benching and even collapse. In addition, fissures are the main channel for groundwater to seep into the tunnel, which may cause water accumulation behind the lining, concrete carbonation or steel corrosion, further expanding the cracks and reducing the structural durability. Once the lining is severely damaged due to the instability of the overburden layer, it may cause the tunnel to be interrupted, resulting in significant economic losses and social impact.

[0003] Currently, the detection methods for submarine tunnel fissures mainly rely on visual detection, geological drilling and empirical judgment based on mechanical reaction curve. However, these methods have the disadvantages of large workload, low efficiency and low precision, which cannot meet the detection requirements of three-dimensional full-parameter and high-precision online identification of fissures. SUMMARY

[0004] The present application aims to overcome the shortcomings of the prior art and provides a three-dimensional identification method for overburden fissures of a submarine tunnel. By using multi-source sound wave signal detection and spatial superposition analysis, the size and position of the fissures can be obtained in real time and accurately, thereby solving the problem of large detection workload, low efficiency and low precision in the prior art, which cannot meet the detection requirements of three-dimensional full-parameter and high-precision online identification of fissures.

[0005] To achieve the above purpose, the present application adopts the following technical solution: A three-dimensional identification method for overburden fissures of a submarine tunnel, comprising the following steps: S1. Intervals of sound wave signal transmitters are arranged in the overburden layer of the submarine tunnel; S2. Intervals of sound wave signal receivers are arranged on the inner wall of the submarine tunnel; S3. The sound wave signal transmitters sequentially emit sound wave signals to the inner wall of the submarine tunnel; S4. The sound wave signal receivers respectively collect the sound wave signals emitted by different sound wave signal transmitters, and convert the collected sound wave signals into digital signals; S5. The digital signals are denoised to obtain the envelope lines in the denoised digital signals, and decay detection is performed to determine the decay points and extract the coordinates of the decay boundary points; S6, moving all acoustic wave signal receivers along the axial direction of the seabed tunnel, continuously extracting the coordinates of new decay boundary points during the movement; S7, fitting the decay boundary points based on the same acoustic wave signal transmitter into a closed surface as the bottom surface of a cone, and taking the acoustic wave signal transmitter as the vertex of the cone, to construct a special-shaped bottom cone; S8, superimposing special-shaped bottom cones from different acoustic wave signal transmitters to form an intersection body; S9, calculating and outputting the geometric characteristics and spatial position of the intersection body; In S7, the process of fitting the decay boundary points into a closed surface includes: According to the decay amplitude, each decay boundary point is assigned a weight The calculation formula of the weight is as follows:

[0006] Wherein, the maximum decay amplitude , is the decay amplitude of the i-th decay boundary point; is the undetermined index, ; ; Select a spatial surface model , use a quadratic polynomial to calculate, and the calculation formula is as follows:

[0007] Wherein, is the undetermined parameter of the bottom surface of the cone, , is the transpose; is the Cartesian coordinate; Use the least square method to determine , and the calculation formula is as follows:

[0008] Wherein, is the coordinate of the i-th decay boundary point; The array of optimal parameters

[0009] can be obtained, and the closed surface equation is ; ; The process of constructing a special-shaped bottom cone includes: Taking the position of the corresponding acoustic wave signal transmitter as the vertex of the cone, connecting with each point on the closed surface equation , a special-shaped bottom cone with as the vertex and the closed surface as the bottom boundary can be constructed.

[0010] Further, in the S5, the de-noising processing adopts a wavelet de-noising method, including wavelet decomposition, threshold processing, up-sampling and wavelet inverse reconstruction.

[0011] Further, in the S4, the sound wave signal is converted into an original discrete signal wherein, is a sample index; in the S5, the process of wavelet decomposition includes: to the original discrete signal , let the approximation coefficient of the 0th layer be , let the detail coefficient of the 0th layer be , let the approximation coefficient of the 0th layer be , then the decomposition formula from the 0th layer to the 0th layer is as follows:

[0012] wherein, is the approximation coefficient of the 0th layer; is the detail coefficient of the 0th layer; is an analysis low-pass filter; is an analysis high-pass filter; is the tap number of the filter. Further, in the S5, the calculation formula of the threshold processing is as follows:

[0013]

[0014] wherein, is a sign function; is the detail coefficient of a certain layer obtained through wavelet decomposition; is a detail coefficient shrinkage threshold.

[0015] Further, in the S5, the calculation formula of the up-sampling is as follows:

[0016] wherein, is the detail coefficient after the up-sampling processing; is the approximation coefficient after the up-sampling processing.

[0017] Further, in the S5, the calculation formula of the wavelet inverse reconstruction is as follows:

[0018] wherein, is the 0th​​​​​​ Layer reconstruction signal; To synthesize a low-pass filter; To synthesize a high-pass filter.

[0019] Furthermore, in step S5, the process of obtaining the envelope of the denoised digital signal includes: The denoised digital signal is subjected to Hilbert transform to obtain the analytic signal, and the magnitude of the analytic signal is calculated. The magnitude of the analytic signal is the envelope of the denoised digital signal.

[0020] Furthermore, in S8, the calculation formula for the overlapping body is as follows:

[0021] in, It is an overlapping body. For the first A set of voxels for an irregularly shaped cone.

[0022] Furthermore, in S9, the process of calculating the geometric features and spatial position of the overlapping bodies includes: Overlapping bodies Coordinates of each voxel A set whose spatial location passes through the centroid. The calculation is as follows:

[0023] in, The coordinates are those of the centroid. The number of voxels in the overlapping body; The coordinates of each voxel of the overlapping volume; Overlapping volumes are represented using axis-aligned bounding boxes. The spatial range is calculated using the following formula:

[0024] in, The coordinates of the smallest voxel; The coordinates of the largest voxel; Perform principal axis analysis on the point set and construct the covariance matrix. as follows:

[0025] For covariance matrix Perform eigenvalue decomposition to obtain 3 eigenvalues. and the corresponding feature vector ; The eigenvalues ​​are sorted from largest to smallest as follows:

[0026] at this time, The direction of maximum variance; The direction of the second largest variance; The direction with the minimum variance; The point set is projected onto the principal axis coordinates as shown below:

[0027] in, The points are respectively at Components in the coordinate system; The three-dimensional scale of the crack is defined as follows:

[0028] in, The length along the principal axis of the crack. The width along the direction of the fracture propagation. The thickness is along the normal direction of the fracture surface; Finally, a set of parameters for the geometric features and spatial location of the overlapping body is obtained. .

[0029] This invention achieves online, rapid, and accurate detection of overburden fractures in submarine tunnels through a multi-source spatial superposition identification algorithm based on acoustic wave decay characteristics. Compared with existing technologies, the advantages of this invention are: (1) High detection accuracy: By utilizing the principle of multi-path superposition, the measurement error of a single path can be eliminated, thereby improving the positioning accuracy; (2) Online monitoring: The overall automation of the identification system is high, and it can quickly collect and process signals to meet the online monitoring needs during the construction process; (3) Easy to deploy: The equipment is flexible to install. The transmitter only needs to be placed outside the overburden layer, reducing the need to modify the internal structure of the tunnel. Attached Figure Description

[0030] Figure 1 This is a schematic diagram of the structure of the submarine overlying rock strata in the yOz coordinate system. Figure 2 This is a schematic diagram of the structure of the seafloor overlying rock layer in the yOx coordinate system; Figure 3 A schematic diagram showing the superposition of two irregularly shaped bottom cones; Figure 4 This is a schematic diagram of the signal receiving unit.

[0031] The attached figures are labeled as follows: 1. Acoustic signal transmitter; 2. Acoustic signal receiver; 3. Seabed overburden; 4. Submarine tunnel; 5. Test stand. Detailed Implementation

[0032] To enable those skilled in the art to better understand the present application, the technical solutions in the embodiments of the present application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present application, and not all embodiments. Based on the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present application.

[0033] This embodiment provides a method for three-dimensional identification of overburden fractures in submarine tunnels. Based on the principle of acoustic tomography, this method analyzes multipath signals and combines spatial reconstruction and signal superposition algorithms to achieve high-precision three-dimensional reconstruction of overburden fractures. The technical principle of this method is as follows: when there are no fractures in the propagation path between the acoustic signal transmitter and receiver, the received signal is stable; however, when fractures exist in the path, the received signal will exhibit significant attenuation or a sharp drop. The degree of signal attenuation is closely related to the size and spatial location of the fracture and the acoustic propagation path.

[0034] A method for three-dimensional identification of overburden fissures in submarine tunnels includes the following steps: Step S1: Place multiple (≥2) acoustic signal transmitters 1 into pre-drilled holes in the seabed overburden 3 along the preset distribution locations.

[0035] The distribution locations of the acoustic signal transmitters 1 are calculated based on the geometric parameters of the submarine tunnel 4 and the expected fracture distribution. Specifically, the spacing and depth of the acoustic signal transmitters 1 are determined using the following formula:

[0036] in, The spacing is set for the acoustic signal transmitter 1; The width of the undersea tunnel is 4. The depth of the sound wave signal transmitter 1; Empirical coefficient (usually) The above formula ensures that the acoustic signal transmitter 1 can cover the entire cross-section of the undersea tunnel 4, while avoiding strong scattering and interference near the boundary of the undersea tunnel 4.

[0037] For example, when the width of the undersea tunnel At that time, the spacing of the acoustic signal transmitters Depth of acoustic signal transmitter .

[0038] Furthermore, the pre-drilled holes in the seabed overburden strata 3 are filled with an acoustic coupling agent (silicone-based material) to enhance the transmission of acoustic signals. In actual operation, the distribution positions of the acoustic signal transmitters 1 are first calculated, then pre-drilled holes are made at the distribution positions, and then multiple acoustic signal transmitters 1 are installed in the pre-drilled holes respectively. The installation angle of the acoustic signal transmitters 1 is perpendicular to the axis of the seabed tunnel 4 to ensure signal coverage.

[0039] Step S2: Install signal receiving units on the inner wall of the undersea tunnel 4.

[0040] The signal receiving unit includes multiple acoustic signal receivers (2 acoustic sensors) spaced apart circumferentially along the inner wall of the seabed tunnel 4. Further, to facilitate the simultaneous axial movement of the multiple acoustic signal receivers 2 along the seabed tunnel 4, in this embodiment, the signal receiving unit also includes a platform 5 for mounting the multiple acoustic signal receivers. Specifically, the main body of the platform 5 is made of high-strength aluminum alloy and adopts a multi-section telescopic curved rod design. The curved rod is made of an elastic material (such as carbon fiber composite material) and has a built-in electric telescopic mechanism (such as a linear motor or hydraulic cylinder). The length of the curved rod is adjusted by a control system, ensuring that the multiple acoustic signal receivers 2 are spaced apart while remaining close to the inner wall of the seabed tunnel 4. Guide rails are installed along the axial direction of the curved rods, allowing the acoustic signal receivers 2 to move along the guide rails, i.e., the acoustic signal receivers 2 can move circumferentially along the inner wall of the seabed tunnel 4. The movement accuracy is controlled by an encoder (accuracy up to ±1mm). Furthermore, the bottom of the platform 5 is equipped with casters and a locking mechanism for easy movement and stopping along the inner wall of the seabed tunnel 4. The spacing of the acoustic signal receivers 2 is determined according to the wavelength of the acoustic waves. Determine the spacing of the acoustic signal receiver 2. To avoid spatial overlap, the control system adjusts the extension and retraction of the crank to ensure that the acoustic signal receiver 2 is in close contact with the rock wall (a pressure sensor ensures constant contact force). The initial position of the acoustic signal receiver 2 is close to the working face of the undersea tunnel 4.

[0041] Step S3: Activate all acoustic signal transmitters 1 to sequentially transmit acoustic signals (frequency range 1-20kHz, pulse width adjustable) into the inner wall of the underwater tunnel 4. The transmission sequence is synchronized by the control unit to avoid interference.

[0042] Step S4: Collect sound wave signals emitted by different sound wave signal transmitters 1 through the sound wave signal receiver 2, and convert the collected sound wave signals into digital signals.

[0043] Specifically, the sound wave signal is converted into a digital signal through the data acquisition module. , For the original discrete signal (discrete samples of the time-domain signal), where, This is the sample index, and its value range is... ; This represents the total number of signal samples.

[0044] Furthermore, the intensity and frequency domain characteristics of the signals received by each acoustic signal receiver 2 are collected through the data acquisition module.

[0045] For digital signals Frequency domain features can be obtained by applying the Fast Fourier Transform (FFT) algorithm. The calculation formula for the FFT is as follows:

[0046] in, For the first in the frequency domain Complex representation of each frequency component It is the imaginary unit.

[0047] Frequency domain characteristics include spectral amplitude Phase Among them, the spectral amplitude Used to calculate bandwidth attenuation and as a spatial reconstruction weight, phase Used for propagation delay estimation and as a constraint for geometric fitting. Data is transmitted to the signal processing unit in real time.

[0048] More specifically, the data acquisition module includes high-speed ADC (analog-to-digital converter) and signal conditioning circuitry, with a sampling rate no less than twice the sound wave frequency (e.g., when the sound wave frequency is 10kHz, the sampling rate is ≥20kHz). Time-domain characteristics include signal amplitude and envelope. Frequency-domain characteristics are obtained through FFT (Fast Fourier Transform), including spectral energy and center frequency shift.

[0049] Step S5: The digital signal is denoised by the signal processing unit to obtain the envelope of the denoised digital signal, and decay detection is performed to determine the decay point and extract the coordinates of the decay boundary point.

[0050] Step S51: Use wavelet denoising method to denoise the original discrete signal after conversion by the data acquisition module. The noise reduction process is performed and the envelope signal is extracted. Specifically, the noise reduction process includes wavelet decomposition, thresholding, upsampling, and inverse wavelet reconstruction.

[0051] 1. The specific process of wavelet decomposition is as follows: For the original discrete signal Let the first The approximation coefficient of the layer is , No. The detail factor of the layer is ,make That is, the approximation coefficients of the 0th layer (initial layer) approximate the original discrete-time series. The low-pass filter is analyzed as follows: Analyze the high-pass filter as The filter length is Then the first layer to the first The decomposition formula for the layer is as follows:

[0052]

[0053] in, For the first Approximation coefficients of the layer; For the first The detail factor of the layer, This refers to the tap number of the filter, which is also the index variable for the convolution summation. In this embodiment, The range of values ​​is .

[0054] 2. The specific process of threshold processing is as follows: Thresholding is applied to the detail coefficients of each layer, and the calculation formula is as follows:

[0055] in, Let be the detail coefficients of a certain layer obtained by wavelet decomposition, i.e., the th Layer detail factor Or perhaps the first The detail factor of the layer is ; This is the sign function, used to preserve the sign or phase direction of detail coefficients; This refers to the detail factor shrinkage threshold; specifically, the detail factor shrinkage threshold. The calculation formula is as follows:

[0056] in, The standard deviation is the coefficient of detail. The total number of signal samples. For sample indexing, more specifically, The possible values ​​for are as follows:

[0057] The thresholding calculation formula shifts the absolute values ​​of all coefficients towards zero. Units, for small noise figure ( The coefficients are directly set to zero to achieve noise reduction, preserving meaningful coefficients but slightly reducing their amplitude (introducing bias). This helps avoid abrupt artifacts near the threshold during noise reduction. Due to the approximate coefficients... Characterizing the low-frequency components of a signal (i.e., the overall trend), noise is typically distributed in the high-frequency components (i.e., the detail coefficients). If the low-frequency approximation coefficients remain unchanged, the high-frequency detail coefficients will be thresholded, thus affecting the approximation coefficients. Thresholding can destroy the low-frequency structure and trends of a signal, and affect the detail coefficients. Denoising can suppress high-frequency noise while preserving low-frequency information of the signal, and usually performs better in terms of mean square error (MSE).

[0058] 3. The specific process of upsampling is as follows: Use the detail factor after thresholding. and the approximate coefficients without thresholding. Upsampling is performed, and the calculation formula is as follows:

[0059] in, These are the detail coefficients after upsampling. These are the approximate coefficients after upsampling.

[0060] By upsampling, the sequence length is doubled (the number of points is doubled), zero values ​​are inserted at odd positions, and the coefficients at even positions are the original values, resulting in a higher resolution sequence.

[0061] 4. The specific process of wavelet inverse reconstruction is as follows: Use the detail coefficients after upsampling. and the approximation coefficients after upsampling. The inverse wavelet reconstruction is calculated using the following formula:

[0062] in, For the first Layer reconstruction signal; To synthesize a low-pass filter; To synthesize a high-pass filter.

[0063] The above formula uses a synthetic filter bank and Perform inverse wavelet transform, i.e., from the th Layer by layer reconstruction, the reconstruction signal of layer 0 is finally obtained. .

[0064] Let the denoised discrete time series The denoised discrete time series Considered as a continuous time signal At discrete time The sampled value at that point, i.e. .

[0065] right Perform Hilbert transform (here) It is a "hypothetical" continuous signal; the actual calculation is performed on a denoised discrete-time series. The Hilbert transform is performed on the above (the principle is given using a continuous formula), and the calculation formula is as follows:

[0066] in, To analyze the signal; for Hilbert transform; It is the imaginary unit.

[0067] The magnitude of the analytic signal is calculated using the following formula:

[0068] in, To analyze the magnitude of the signal, also known as the envelope signal (signal envelope).

[0069] Step S52: Determine the decay point by performing decay detection based on the threshold method.

[0070] The entire envelope signal With decay event detection threshold In comparison, any The point is determined to be the decay point.

[0071] Furthermore, the decay event detection threshold It is not a fixed absolute value, but rather it is adaptively determined based on the background noise level. The specific determination method is as follows: ① Calculate the statistical characteristics of the signal amplitude (envelope value), such as the mean signal strength, in the initial segment of the signal or in a segment known to have no decay events. and standard deviation ; ② Set a threshold The calculation formula is as follows:

[0072] in, It is a constant, chosen according to the required confidence level (usually 3 to 5).

[0073] For example, With a confidence level of 99.7% (assuming the noise follows a Gaussian distribution), this means that signals exceeding this threshold have a high probability of being genuine decay events rather than noise.

[0074] Step S53: Extract the coordinates of the decay boundary points based on the decay points.

[0075] The decay boundary points are extracted from the transmit-receive path. The spatial coordinates of the decay boundary points are: The envelope amplitude corresponding to the decay boundary point is As an index of decay intensity at the decay boundary point, where, To determine the decay boundary time on this path using threshold determination, specifically, .

[0076] Specifically, the coordinates of the decay boundary points are recorded by an acoustic signal receiver and their absolute coordinates are obtained in conjunction with a bench positioning system (such as a laser rangefinder or IMU).

[0077] More specifically, it can be understood as comparing the envelope with a threshold. Points on the envelope that are below the threshold are decay points. The curve (arc) formed by connecting multiple decay points is the decay line. The two ends of the decay line are the decay boundary points. A decay line has two decay boundary points. The specific coordinates of the decay boundary points are determined by the acoustic signal receiver that is closest to the decay boundary point.

[0078] Furthermore, the phase information from step S4 (i.e., phase) is utilized. Accurately calculating the propagation delay of sound waves allows for a more precise determination of the spatial location of the decay point. Since the propagation delay is distance-dependent, phase information can help calibrate the position of the decay point relative to the sound wave transmitter and receiver, reducing positioning errors.

[0079] Step S6: Move all acoustic signal receivers 2 along the seabed tunnel 4 axially, and continuously extract the coordinates of new decay boundary points during the movement.

[0080] The platform 5 is moved towards the tunnel face, causing all the acoustic signal receivers 2 to move synchronously along the underwater tunnel 4 axially. Specifically, the movement mechanism of the platform 5 is as follows: the platform 5 moves along the inner wall track of the underwater tunnel 4 via an electric pulley system, and the movement speed is... To ensure sufficient sampling.

[0081] Repeat steps S4 and S5, obtaining the coordinates of the new decay boundary point after each move.

[0082] By repeating the coordinate extraction process described above for all decay point locations corresponding to the same acoustic signal transmitter, a set of discrete decay boundary points can be obtained, the mathematical expression of which is:

[0083] in, It is the first Spatial coordinates of the decay boundary point For the first The decay amplitude at each decay boundary point.

[0084] The coordinate extraction process described above is performed on all acoustic signal transmitters 1 to obtain the discrete set of decay boundary points within the entire seabed overburden 3.

[0085] Step S7: Perform weighted curve fitting on the discrete points of the decay boundary obtained in step S6 to reconstruct the smooth closed boundary line at the bottom of the cone, and construct multiple irregular bottom cones with the corresponding acoustic signal transmitter 1 as the cone vertex.

[0086] Weights are assigned to each decay boundary point based on the magnitude of the decay amplitude. Weight , It is an index to be determined. .

[0087] To describe the decay boundary surrounding the acoustic signal transmitter, a class of closable spatial surface models is selected. The calculation is performed using a quadratic polynomial, and the formula is as follows:

[0088] in, For the parameters to be determined (to be solved) of the base of the cone, specifically... , This is a transpose.

[0089] Determined by the weighted least squares criterion That is, to solve:

[0090] in, For the first The coordinates of the decay boundary points;

[0091] This will give you the array of optimal parameters. ,Right now Make the function At its minimum, the surface equation is: .

[0092] Based on the location of the sound wave signal transmitter As the vertex of the cone, With curved surfaces By connecting the points above, a network can be constructed. An irregularly shaped cone with its vertex at the top and the fitted surface as its bottom boundary. .

[0093] Repeat the above process for all acoustic signal transmitters 1 to obtain a family of cone combinations. The final number of irregularly shaped bottom cones is the same as the number of acoustic signal transmitters 1.

[0094] Using the spectral amplitude information obtained in step S4 (i.e., spectral amplitude) Weights are assigned to each decay boundary point. Decay boundary points with larger amplitudes typically indicate stronger signal decay and may correspond to more significant crack features. When fitting a closed surface, a weighted algorithm is used to give points with larger amplitudes a greater impact on the surface shape, thereby improving reconstruction accuracy.

[0095] Step S8: Stack the irregularly shaped bottom cones from different acoustic signal transmitters 1 to form an overlapping body.

[0096] The superposition uses a voxel-based Boolean intersection operation: Let the first... The voxel set of the irregularly shaped base cone is Then overlapping bodies Overlapping bodies That is, the number of cracks What form?

[0097] Step S9: The crack identification unit calculates and outputs the length, width, and specific spatial location of the identified crack.

[0098] Overlapping bodies Coordinates of each voxel A set whose spatial location passes through the centroid. The calculation is as follows:

[0099] in, The coordinates are those of the centroid. The number of voxels in the overlapping body; For the coordinates of each voxel of the overlapping volume, specifically, .

[0100] Overlapping bodies The spatial extent can be represented using an axis-aligned bounding box as follows:

[0101] in, The coordinates of the smallest voxel; The coordinates of the largest voxel are given.

[0102] Perform principal component analysis (PMAC) on the point set to construct the covariance matrix. As shown below:

[0103] By analyzing the covariance matrix By performing eigenvalue decomposition, we can obtain three eigenvalues. and the corresponding feature vector .

[0104] The eigenvalues ​​are sorted from largest to smallest as follows:

[0105] at this time, The direction of maximum variance can be used to describe the principal axis of crack propagation. The direction of the second largest variance can be used to describe the secondary axis of fracture distribution. The direction with the minimum variance can be used as the normal to the fracture surface.

[0106] The point set is projected onto the principal axis coordinates as shown below:

[0107] in, The points are respectively at Components in the coordinate system.

[0108] The three-dimensional scale of the fracture is defined as follows:

[0109] in, The length along the principal axis of the crack. The width along the direction of the fracture propagation. The thickness is along the normal direction of the fracture surface (which can characterize the thickness of the fracture zone).

[0110] Finally, the parameter set of the output fracture is obtained. Including orientation, tilt angle, and other attitude information, and using 3D rendering software (such as OpenGL) to generate a 3D visualization report, the overlapping bodies... Visualize the parameters and output a parameter report.

[0111] Although the present invention has been described using the above preferred embodiments, it is not intended to limit the scope of protection of the present invention. Any changes and modifications made by those skilled in the art to the above embodiments without departing from the spirit and scope of the present invention shall still fall within the scope of protection of the present invention.

Claims

1. A method for identifying three-dimensional cracks in overburden rock of a submarine tunnel, characterized in that, The method comprises the following steps: S1, arranging multiple acoustic signal transmitters in the overburden layer of the submarine tunnel at intervals; S2, arranging multiple acoustic signal receivers on the inner wall of the submarine tunnel at intervals; S3, sequentially transmitting acoustic signals to the inner wall of the submarine tunnel by the acoustic signal transmitters; S4, collecting the acoustic signals emitted by different acoustic signal transmitters by the acoustic signal receivers respectively, and converting the collected acoustic signals into digital signals; S5, performing denoising processing on the digital signals, obtaining the envelope line in the denoised digital signals, and performing decay detection to determine the decay point and extract the coordinates of the decay boundary point; S6, moving all the acoustic signal receivers along the axial direction of the submarine tunnel, and continuously extracting the coordinates of new decay boundary points during the movement; S7, fitting the decay boundary points based on the same acoustic signal transmitter into a closed surface as the bottom surface of a cone, taking the acoustic signal transmitter as the apex of the cone, and constructing a special-shaped bottom cone; S8, superimposing the special-shaped bottom cones from different acoustic signal transmitters to form an overlapping body, and the overlapping body is the geometric morphology of the fracture; S9, identifying the geometric features and spatial positions of the overlapping body; In the S7, the process of fitting the decay boundary points into a closed surface comprises: The weight is assigned to each decay boundary point according to the decay amplitude size , the weight is calculated according to the following formula: ; wherein the maximum decay amplitude , is the decay amplitude of the first decay boundary point; is the undetermined exponent, ; Selecting a spatial surface model , using a quadratic polynomial, the calculation formula is as follows: ; wherein is a pending parameter for the base of the cone, , is the transpose; are Cartesian coordinates; The least square method is used to determine The calculation formula is as follows: ； wherein is the coordinate of the first decay boundary point; is the coordinate of the second decay boundary point; An array of optimal parameters is obtained The closed surface equation is ; The process of constructing a special-shaped bottom cone comprises: According to the position of the corresponding sound wave signal transmitter As the vertex of the cone, equations of closed surfaces By connecting the points above, a network can be constructed. An irregularly shaped cone with a vertex and a closed surface as its bottom boundary.

2. The method according to claim 1, wherein, In the S5, the denoising processing adopts a wavelet denoising method, which comprises wavelet decomposition, threshold processing, upsampling, and wavelet inverse reconstruction.

3. The method according to claim 2, wherein, In the S4, the acoustic wave signal is converted into the original discrete signal wherein, is the sample index; In the S5, the process of wavelet decomposition comprises: For the original discrete signal Let the first The approximation coefficient of the layer is , No. The detail factor of the layer is Let the approximation coefficient of the 0th layer be... Then the first layer to the first The decomposition formula for the layer is as follows: ； ； wherein is the approximation coefficient of the first layer; is the approximation coefficient of the first layer; is the detail coefficient of the first layer; is the detail coefficient of the first layer; is an analysis low-pass filter; is an analysis high-pass filter; is the tap number of the filter.

4. The method according to claim 3, wherein, In the S5, the calculation formula of threshold processing is as follows: ； wherein is a symbol function; is a detail coefficient of a certain layer obtained by wavelet decomposition; is a detail coefficient shrinkage threshold.

5. The method according to claim 4, wherein, In the S5, the calculation formula of upsampling is as follows: ; ; wherein, is a detail coefficient after up-sampling processing; is an approximation coefficient after up-sampling processing.

6. The method according to claim 5, wherein, In the S5, the calculation formula of wavelet inverse reconstruction is as follows: ； wherein is a first reconstructed signal of a layer; is a synthetic low-pass filter; is a synthetic high-pass filter.

7. The method according to claim 1, wherein, In the S5, the process of obtaining the envelope line in the denoised digital signal comprises: Performing Hilbert transform on the denoised digital signal to obtain an analytical signal, and calculating the modulus value of the analytical signal, which is the envelope line in the denoised digital signal.

8. The method according to claim 1, wherein, In the S8, the calculation formula of the overlapping body is as follows: ； wherein, is an overlap, is the set of voxels of the i-th irregular base cone.

9. The method according to claim 8, wherein, In the S9, the process of calculating the geometric features and spatial positions of the overlapping body comprises: Interleaver is a set of voxel coordinates whose spatial position is calculated by the centroid The calculation formula is as follows: ； wherein, is the centroid coordinate; is the number of voxels of the intersection; is the coordinate of each voxel of the intersection; Representing an intersection using axis-aligned bounding boxes The spatial extent of the intersection body is calculated as follows: ； ； wherein is the minimum voxel coordinate; is the maximum voxel coordinate; Performing principal component analysis on the point set, constructing a covariance matrix As follows: ; On the covariance matrix Do eigen decomposition, get 3 eigenvalues And the corresponding eigenvectors ; The characteristic values are sorted from large to small as follows: ; At this time, is the direction of the largest variance; is the direction of the second largest variance; is the direction of the smallest variance; The point set is projected onto the principal axis coordinates as follows: ； ； ; wherein are the components of the point in the coordinate system, respectively; coordinate system. The three-dimensional scale of the fracture is defined as follows: ; ; ； wherein L is the length along the main axis of the fissure, W is the width along the direction of the fissure, T is the thickness along the normal to the fissure plane; The final set of parameters of the geometry and spatial position of the overlap .

Images