A method for identifying three-dimensional cracks in overburden of a submarine tunnel

Abstract

This invention relates to the field of marine engineering and structural health monitoring technology, specifically a three-dimensional identification method for overburden fissures in submarine tunnels. When there are no fissures in the propagation path between the acoustic signal transmitter and receiver, the received signal is stable; however, when fissures exist in the path, the received signal will show significant attenuation or a sharp drop. This invention achieves online, rapid, and accurate detection of overburden fissures in submarine tunnels through a multi-source spatial superposition identification algorithm based on acoustic attenuation characteristics. Utilizing the principle of multi-path superposition, it can eliminate single-path measurement errors and improve positioning accuracy. The system has a high degree of automation, can quickly acquire and process signals, and meets the online monitoring needs during construction. The equipment is flexible in installation; the transmitter only needs to be placed outside the overburden layer, reducing modifications to the internal structure of the tunnel.

Description

Technical Field

[0001] This invention relates to the field of marine engineering and structural health monitoring technology, specifically a three-dimensional identification method for overburden fractures in submarine tunnels. Background Technology

[0002] As the scale and depth of undersea tunnel construction continue to increase, the geological structure of the overburden strata becomes increasingly complex and variable. The degree of fissure development has a significant impact on the overall stability and safety of the tunnel. Specifically, fissures in the overburden strata may weaken the self-stabilizing ability of the surrounding rock, leading to uneven settlement of the foundation or deformation of the surrounding rock. This, in turn, causes uneven pressure on the tunnel lining, resulting in cracks, misalignment, or even the risk of collapse. Furthermore, fissures are the main channels for groundwater to seep into the tunnel, which may lead to water accumulation behind the lining, concrete carbonization, or steel corrosion, further enlarging the cracks and reducing structural durability. Moreover, if the undersea tunnel lining is severely damaged due to overburden instability, it may cause the tunnel to be closed to traffic, resulting in significant economic losses and social impact.

[0003] Currently, methods for detecting cracks in submarine tunnels mainly rely on visual inspection, geological drilling, and empirical judgment based on mechanical reaction curves. These methods suffer from drawbacks such as high workload, low efficiency, and low accuracy, making it difficult to meet the requirement for high-precision online identification of cracks with full three-dimensional parameters. Summary of the Invention

[0004] The purpose of this invention is to address the shortcomings of existing technologies by proposing a three-dimensional identification method for overburden fractures in submarine tunnels. This method uses multi-source acoustic signal detection and spatial superposition analysis to obtain the size and location of fractures in real time and with high accuracy. This solves the problems of large workload, low efficiency, and low accuracy in existing technologies, which make it difficult to meet the detection requirements for three-dimensional full-parameter and high-precision online identification of fractures.

[0005] To achieve the above objectives, the present invention adopts the following technical solution:

[0006] A method for three-dimensional identification of overburden fissures in submarine tunnels includes the following steps:

[0007] S1. Multiple acoustic signal transmitters are deployed at intervals within the overlying rock layer of the submarine tunnel.

[0008] S2. Multiple acoustic signal receivers are installed at intervals on the inner wall of the submarine tunnel.

[0009] S3. Transmit sound wave signals sequentially into the inner wall of the undersea tunnel using a sound wave signal transmitter;

[0010] S4. Acquire sound wave signals emitted by different sound wave signal transmitters through a sound wave signal receiver, and convert the acquired sound wave signals into digital signals.

[0011] S5. Perform denoising on the digital signal, obtain the envelope of the denoised digital signal, perform attenuation detection, determine the attenuation point, and extract the coordinates of the attenuation boundary point.

[0012] S6. Move all acoustic signal receivers axially along the seabed tunnel and continuously extract the coordinates of new decay boundary points during the movement.

[0013] S7. Fit the decay boundary points based on the same acoustic signal transmitter into a closed surface as the bottom surface of the cone, and use the acoustic signal transmitter as the vertex of the cone to construct an irregular bottom cone.

[0014] S8. Stack the irregularly shaped bottom cones from different acoustic signal transmitters to form an overlapping body;

[0015] S9. Calculate and output the geometric features and spatial position of the overlapping body;

[0016] In step S7, the process of fitting the decay boundary points into a closed surface includes:

[0017] Weights are assigned to each decay boundary point based on the magnitude of the decay amplitude. Weight The calculation formula is as follows:

[0018]

[0019] Among them, the maximum decay amplitude , For the first The decay amplitude at each decay boundary point; The index is yet to be determined. ;

[0020] Selecting a spatial surface model The calculation is performed using a quadratic polynomial, and the formula is as follows:

[0021]

[0022] in, For the undetermined parameters of the base of the cone, , For transpose; Cartesian coordinates;

[0023] Determined using the least squares method The calculation formula is as follows:

[0024]

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

[0026] This will give you the array of optimal parameters. At this point, the equation of the closed surface is ;

[0027] The process of constructing an irregularly shaped cone includes:

[0028] 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.

[0029] Furthermore, in S5, the denoising process employs a wavelet denoising method, including wavelet decomposition, thresholding, upsampling, and inverse wavelet reconstruction.

[0030] Furthermore, in step S4, the acoustic wave signal is converted into the original discrete signal. ,in, For sample indexing; in S5, the wavelet decomposition process includes:

[0031] 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:

[0032]

[0033] in, For the first Approximation coefficients of the layer; For the first The detail factor of the layer; To analyze low-pass filters; To analyze high-pass filters; This is the tap number of the filter.

[0034] Furthermore, in S5, the calculation formula for the threshold processing is as follows:

[0035]

[0036] in, It is a symbolic function; These are the detail coefficients of a certain layer obtained through wavelet decomposition; The threshold for shrinking the detail coefficients.

[0037] Furthermore, in S5, the calculation formula for upsampling is as follows:

[0038]

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

[0040] Furthermore, in S5, the calculation formula for wavelet inverse reconstruction is as follows:

[0041]

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

[0043] Furthermore, in step S5, the process of obtaining the envelope of the denoised digital signal includes:

[0044] 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.

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

[0046]

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

[0048] Furthermore, in S9, the process of calculating the geometric features and spatial position of the overlapping bodies includes:

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

[0050]

[0051] 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;

[0052] Overlapping volumes are represented using axis-aligned bounding boxes. The spatial range is calculated using the following formula:

[0053]

[0054] in, The coordinates of the smallest voxel; The coordinates of the largest voxel;

[0055] Perform principal axis analysis on the point set and construct the covariance matrix. as follows:

[0056]

[0057] For covariance matrix Perform eigenvalue decomposition to obtain 3 eigenvalues. and the corresponding feature vector ;

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

[0059]

[0060] at this time, The direction of maximum variance; The direction of the second largest variance; The direction with the minimum variance;

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

[0062]

[0063] in, The points are respectively at Components in the coordinate system;

[0064] The three-dimensional scale of the crack is defined as follows:

[0065]

[0066] 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;

[0067] Finally, a set of parameters for the geometric features and spatial location of the overlapping body is obtained. .

[0068] 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:

[0069] (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;

[0070] (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;

[0071] (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

[0072] Figure 1 This is a schematic diagram of the structure of the submarine overlying rock strata in the yOz coordinate system.

[0073] Figure 2 This is a schematic diagram of the structure of the seafloor overlying rock layer in the yOx coordinate system;

[0074] Figure 3 A schematic diagram showing the superposition of two irregularly shaped bottom cones;

[0075] Figure 4 This is a schematic diagram of the signal receiving unit.

[0076] The attached figures are labeled as follows:

[0077] 1. Acoustic signal transmitter; 2. Acoustic signal receiver; 3. Seabed overburden; 4. Submarine tunnel; 5. Test stand. Detailed Implementation

[0078] 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.

[0079] 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.

[0080] A method for three-dimensional identification of overburden fissures in submarine tunnels includes the following steps:

[0081] Step S1: Place multiple (≥2) acoustic signal transmitters 1 into pre-drilled holes in the seabed overburden 3 along the preset distribution locations.

[0082] 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:

[0083]

[0084] 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.

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

[0086] 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.

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

[0088] 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.

[0089] 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.

[0090] 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.

[0091] 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.

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

[0093] 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:

[0094]

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

[0096] 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.

[0097] 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.

[0098] 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.

[0099] 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.

[0100] 1. The specific process of wavelet decomposition is as follows:

[0101] 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:

[0102]

[0103]

[0104] 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 .

[0105] 2. The specific process of threshold processing is as follows:

[0106] Thresholding is applied to the detail coefficients of each layer, and the calculation formula is as follows:

[0107]

[0108] 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:

[0109]

[0110] 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:

[0111]

[0112] 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).

[0113] 3. The specific process of upsampling is as follows:

[0114] Use the detail factor after thresholding. and the approximate coefficients without thresholding. Upsampling is performed, and the calculation formula is as follows:

[0115]

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

[0117] 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.

[0118] 4. The specific process of wavelet inverse reconstruction is as follows:

[0119] Use the detail coefficients after upsampling. and the approximation coefficients after upsampling. The inverse wavelet reconstruction is calculated using the following formula:

[0120]

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

[0122] 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. .

[0123] 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. .

[0124] 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:

[0125]

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

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

[0128]

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

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

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

[0132] 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:

[0133] ① 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 ;

[0134] ② Set a threshold The calculation formula is as follows:

[0135]

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

[0137] 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.

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

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

[0140] 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).

[0141] 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.

[0142] 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.

[0143] 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.

[0144] 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.

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

[0146] 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:

[0147]

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

[0149] 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.

[0150] 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.

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

[0152] 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:

[0153]

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

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

[0156]

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

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

[0159] 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. .

[0160] 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.

[0161] 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.

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

[0163] 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

[0164] What form?

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

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

[0167]

[0168] 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, .

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

[0170]

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

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

[0173]

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

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

[0176]

[0177] 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.

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

[0179]

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

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

[0182]

[0183] 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).

[0184] 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.

[0185] 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 three-dimensional identification of overburden fissures in submarine tunnels, characterized in that, Includes the following steps: S1. Multiple acoustic signal transmitters are deployed at intervals within the overlying rock layer of the submarine tunnel. S2. Multiple acoustic signal receivers are installed at intervals on the inner wall of the submarine tunnel. S3. Transmit sound wave signals sequentially into the inner wall of the undersea tunnel using a sound wave signal transmitter; S4. Acquire sound wave signals emitted by different sound wave signal transmitters through a sound wave signal receiver, and convert the acquired sound wave signals into digital signals. S5. Perform denoising on the digital signal, obtain the envelope of the denoised digital signal, perform attenuation detection, determine the attenuation point, and extract the coordinates of the attenuation boundary point. S6. Move all acoustic signal receivers axially along the seabed tunnel and continuously extract the coordinates of new decay boundary points during the movement. S7. Fit the decay boundary points based on the same acoustic signal transmitter into a closed surface as the bottom surface of the cone, and use the acoustic signal transmitter as the vertex of the cone to construct an irregular bottom cone. S8. The irregularly shaped cones from different acoustic signal transmitters are superimposed to form an overlapping body, which is the geometry of the crack. S9. Identify the geometric features and spatial location of the overlapping bodies; In step S7, the process of fitting the decay boundary points into a closed surface includes: Weights are assigned to each decay boundary point based on the magnitude of the decay amplitude. Weight The calculation formula is as follows: ; Among them, the maximum decay amplitude , For the first The decay amplitude at each decay boundary point; The index is yet to be determined. ; Selecting a spatial surface model The calculation is performed using a quadratic polynomial, and the formula is as follows: ; in, For the undetermined parameters of the base of the cone, , For transpose; Cartesian coordinates; Determined using the least squares method The calculation formula is as follows: ； in, For the first The coordinates of the decay boundary points; This will give you the array of optimal parameters. At this point, the equation of the closed surface is ; The process of constructing an irregularly shaped cone includes: 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 for three-dimensional identification of overburden fissures in submarine tunnels according to claim 1, characterized in that, In S5, the denoising process employs a wavelet denoising method, including wavelet decomposition, thresholding, upsampling, and inverse wavelet reconstruction.

3. The method for three-dimensional identification of overburden fissures in submarine tunnels according to claim 2, characterized in that, In step S4, the acoustic signal is converted into the original discrete signal. ,in, For sample index; In S5, the wavelet decomposition process includes: 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: ； ； in, For the first Approximation coefficients of the layer; For the first The detail factor of the layer; To analyze low-pass filters; To analyze high-pass filters; This is the tap number of the filter.

4. The method for three-dimensional identification of overburden fissures in submarine tunnels according to claim 3, characterized in that, In step S5, the calculation formula for the threshold processing is as follows: ； in, It is a symbolic function; These are the detail coefficients of a certain layer obtained through wavelet decomposition; The threshold for shrinking the detail coefficients.

5. The method for three-dimensional identification of overburden fissures in submarine tunnels according to claim 4, characterized in that, In S5, the calculation formula for upsampling is as follows: ; ; in, These are the detail coefficients after upsampling. These are the approximate coefficients after upsampling.

6. The method for three-dimensional identification of overburden fissures in submarine tunnels according to claim 5, characterized in that, In S5, the calculation formula for wavelet inverse reconstruction is as follows: ； in, For the first Layer reconstruction signal; To synthesize a low-pass filter; To synthesize a high-pass filter.

7. The method for three-dimensional identification of overburden fissures in submarine tunnels according to claim 1, characterized in that, 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.

8. The method for three-dimensional identification of overburden fissures in submarine tunnels according to claim 1, characterized in that, In S8, the calculation formula for the overlapping body is as follows: ； in, It is an overlapping body. For the first A set of voxels for an irregularly shaped base cone.

9. The method for three-dimensional identification of overburden fissures in submarine tunnels according to claim 8, characterized in that, In step 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: ； 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: ； ； 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: ; 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: ; 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: ； ； ; in, The points are respectively at Components in the coordinate system; The three-dimensional scale of the crack is defined as follows: ; ; ； 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. .

Images