Rock fracture source acoustic emission positioning method and system based on multi-link collaborative optimization

Abstract

The application discloses a rock fracture source acoustic emission positioning method and system based on multi-link cooperative optimization, and the method comprises the following steps: collecting multi-channel acoustic emission signals generated by rock fracture, and screening effective channel signals; constructing an adaptive wavelet threshold function, taking a comprehensive index of signal-to-noise ratio and root mean square error, and taking signal overall variation as a constraint term, optimizing the wavelet threshold function through a particle swarm algorithm, and carrying out wavelet threshold denoising; performing complete set empirical mode decomposition on the denoised signal and reconstructing component signals containing high-frequency characteristics, locking a narrow time window containing the first arrival time on the component signals, picking up the first arrival time of each channel in the narrow time window based on a statistical information criterion; determining at least one reference sensor, and outputting the spatial coordinates of the fracture source based on the first arrival time difference of each effective channel and the reference sensor. Therefore, the positioning accuracy and robustness of the rock fracture source in a low signal-to-noise ratio and strong interference environment are improved.

Description

Technical Field

[0001] This invention relates to the field of rock acoustic emission monitoring and non-destructive testing technology, and in particular to a method and system for locating rock fracture sources based on multi-stage collaborative optimization. Background Technology

[0002] Major near-fault engineering projects and deep-earth resource development urgently require the understanding of the fracturing evolution mechanism of rock masses in fault zones. Earthquake fault zone dynamic disaster simulation experiments are a core means of reproducing fault displacement, rock fracturing, and energy release. Acoustic emission technology can capture micro-fracture signals in rocks in real time, making it an important technical approach for locating fracture sources and monitoring damage evolution.

[0003] In the field of rock acoustic emission signal processing, common existing processing methods include: threshold denoising methods based on wavelet transform, first-arrival time picking methods based on energy ratio methods or Akaike information criteria, and location inversion methods based on time difference methods. These methods are used to varying degrees in applications such as non-destructive testing in civil engineering and rock mechanics testing and monitoring. Furthermore, for the effectiveness evaluation of multi-channel sensor data, existing technologies have also attempted to use statistical analysis methods to make preliminary judgments on signal quality.

[0004] However, the acoustic emission signal from rock fracture has complex characteristics such as high frequency, non-stationarity, and low signal-to-noise ratio. Furthermore, the signal is accompanied by attenuation, dispersion, and waveform distortion during propagation in the rock mass. These factors place high demands on the accuracy of signal processing and positioning results. Summary of the Invention

[0005] This application provides a method for acoustic emission localization of rock fracture sources based on multi-stage collaborative optimization. Through multi-stage collaborative optimization of sensor selection, adaptive denoising, joint acquisition and localization inversion, the method significantly improves the localization accuracy and robustness of rock fracture sources in low signal-to-noise ratio and strong interference environments, and achieves highly reliable acoustic emission localization with full-process automation.

[0006] This application provides a method for locating acoustic emission sources of rock fractures based on multi-stage collaborative optimization, including: S101, acquiring multi-channel acoustic emission signals generated by rock fractures, extracting characteristic parameters of each channel signal, calculating the Pearson correlation coefficient between channels based on the characteristic parameters, and selecting effective channel signals according to the Pearson correlation coefficient;

[0007] S102, construct an adaptive wavelet threshold function that is continuously differentiable at the threshold and introduces an adjustable curvature parameter. Use the comprehensive index of signal-to-noise ratio and root mean square error as the fitness and the overall variation of the signal as the constraint term. Optimize the curvature parameter, threshold parameter and wavelet decomposition level of the wavelet threshold function through particle swarm optimization algorithm. Use the optimized parameters to perform wavelet threshold denoising on the effective channel signal to obtain the denoised signal. S103, perform complete set empirical mode decomposition on the denoised signal and reconstruct the component signal containing high-frequency features. Calculate the energy difference within the time window before and after the component signal to lock a narrow time window containing the arrival time of the first wave. Pick the arrival time of the first wave of each channel within the narrow time window based on statistical information criteria. S104, determine at least one reference sensor, and based on the arrival time difference of the first arrival wave between each effective channel and the reference sensor, construct a set of positioning equations using the time difference method and perform inversion to output the spatial coordinates of the rupture source.

[0008] Preferably, in step S104, determining at least one reference sensor is configured as follows: among all effective channels, the sensor with the shortest first arrival time is selected as the reference sensor.

[0009] Preferably, the characteristic parameters include at least one or more of the following: count, duration, rise time, peak frequency, frequency center, and the ratio of rise time to amplitude.

[0010] Preferably, the step of selecting effective channel signals based on the Pearson correlation coefficient specifically includes: For each channel i, calculate its average correlation coefficient with all other channels, select the channels with an average correlation coefficient greater than a preset empirical threshold and the top Q channels in descending order as valid channels, with the value of Q determined according to the positioning scenario; the data of the remaining channels are discarded.

[0011] Preferably, the adaptive wavelet threshold function in S102 is expressed as:

[0012] in, These are the original wavelet coefficients. These are the wavelet coefficients after thresholding. For threshold parameters, This is an adjustable curvature parameter.

[0013] Preferably, the complete set empirical mode decomposition in S103 includes: White noise of different amplitudes is added to the original signal multiple times. Empirical mode decomposition is performed on the signal after each addition of noise, and the decomposition results are averaged to obtain multiple intrinsic mode function components.

[0014] Preferably, the reconstructed component signal containing high-frequency features specifically comprises: The first M high-frequency intrinsic mode function components are selected and superimposed, where M is an integer greater than or equal to 2 and less than or equal to 4; Locking the narrow time window containing the arrival time of the first wave in S103 also includes: Simultaneously calculate the corrected energy ratio curve and the long-short time average ratio curve on the component signal. Determine multiple candidate narrow time windows based on the energy difference, the corrected energy ratio, and the long-short time average ratio, respectively. Use the fine search window obtained by fusing the multiple candidate narrow time windows as the final narrow time window.

[0015] Preferably, in step S104, determining at least one reference sensor may further include: The positioning equations are solved sequentially using each effective channel as a reference sensor to obtain multiple candidate coordinates. The final rupture source coordinates are determined based on the posterior residuals of each candidate coordinate.

[0016] Preferably, the determination of the final fracture source coordinates based on the posterior residuals of each candidate coordinate specifically involves: Based on the posterior residuals, several candidate coordinates are selected and weighted averaged to obtain the final rupture source coordinates:

[0017] Where S is the selected set of candidate coordinate indices, and the weights of each index are inversely proportional to the square of the posterior residuals: , The coordinates of the final rupture source, These are candidate coordinates.

[0018] This application also proposes a rock fracture source acoustic emission localization system based on multi-stage collaborative optimization, including: an extraction module, a noise reduction module, a decomposition and reconstruction module, and a localization module; The extraction module is used to collect multi-channel acoustic emission signals generated by rock fracturing, extract the characteristic parameters of each channel signal, calculate the Pearson correlation coefficient between channels based on the characteristic parameters, and filter out the effective channel signals according to the Pearson correlation coefficient. The denoising module is used to construct an adaptive wavelet threshold function that is continuously differentiable at the threshold and introduces an adjustable curvature parameter. Using a comprehensive index of signal-to-noise ratio and root mean square error as the fitness, the curvature parameter, threshold parameter, and wavelet decomposition level of the wavelet threshold function are optimized by particle swarm optimization algorithm. The optimized parameters are then used to perform wavelet threshold denoising on the effective channel signal to obtain a denoised signal. The decomposition and reconstruction module is used to perform complete set empirical mode decomposition on the denoised signal and reconstruct the component signal containing high-frequency features, calculate the energy difference within the time window before and after the component signal to lock a narrow time window containing the arrival time of the first wave, and pick the arrival time of the first wave of each channel within the narrow time window based on statistical information criteria. The positioning module is used to determine at least one reference sensor, and based on the arrival time difference of the first wave between each effective channel and the reference sensor, construct a set of positioning equations using the time difference method and perform inversion to output the spatial coordinates of the rupture source.

[0019] One or more technical solutions provided in this application have at least the following technical effects or advantages: By performing a sensor selection step based on Pearson correlation coefficient before positioning inversion, characteristic-level correlation analysis is conducted on the acoustic emission signals of each channel to actively screen out high-quality effective channel signals. This ensures that the signals entering the subsequent processing link have reliable data sources, thereby avoiding positioning errors caused by sensor failure or waveform distortion from the source. This significantly enhances the adaptability of the entire positioning system to sensor deployment conditions and on-site noise environment.

[0020] By constructing a continuously differentiable wavelet threshold function and introducing a particle swarm optimization algorithm for global parameter optimization, high-fidelity processing of low signal-to-noise ratio signals is achieved in the denoising stage. This threshold function is continuously differentiable at the threshold, and does not have the discontinuity problem of traditional hard threshold functions, nor the constant deviation of soft threshold functions. After optimization by the particle swarm optimization algorithm, the denoised signal can effectively retain the high-frequency abrupt change components in the rock fracture signal, while suppressing pseudo-Gibbs oscillations and over-smoothing phenomena.

[0021] By organically integrating complete set empirical mode decomposition, improved energy difference method, and statistical information criterion, a joint strategy of "coarse positioning + fine positioning" is formed in the arrival time picking stage. Complete set empirical mode decomposition effectively suppresses mode aliasing, and the reconstructed high-frequency components enhance the abrupt change characteristics of the first arrival wave. Improved energy difference method quickly locks the narrow time window containing the accurate arrival time, greatly reducing the search range. Statistical information criterion performs fine picking within this narrow window, avoiding the problem of being easily disturbed by local extrema during global search. The synergistic effect of the three enables stable and low-error first arrival wave arrival times to be obtained even in environments with strong tailwave interference and low signal-to-noise ratio.

[0022] The four stages of sensor selection, adaptive noise reduction, accurate acquisition, and positioning inversion form a collaborative link that supports each other and transmits information step by step. Through the collaboration of the entire process, the overall optimal solution is achieved, the positioning error is significantly reduced, and strong robustness is maintained under complex working conditions such as signal attenuation, dispersion, and data loss. Attached Figure Description

[0023] Figure 1 This is a flowchart illustrating the acoustic emission localization method for rock fracture sources based on multi-stage collaborative optimization, according to an embodiment of the present invention. Figure 2 This is a structural block diagram of a rock fracture source acoustic emission localization system based on multi-stage collaborative optimization, according to an embodiment of the present invention. Detailed Implementation

[0024] To facilitate understanding of the present invention, a more complete description of this application will be given below with reference to the accompanying drawings, which illustrate preferred embodiments of the invention. However, the invention can be implemented in many different forms and is not limited to the embodiments described herein. Rather, these embodiments are provided to enable a more thorough and complete understanding of the disclosure of the present invention.

[0025] It should be noted that the terms "vertical," "horizontal," "up," "down," "left," "right," and similar expressions used in this article are for illustrative purposes only and do not represent the only possible implementation.

[0026] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains; the terminology used herein in the description of the invention is for the purpose of describing particular embodiments only and is not intended to limit the invention; the term "and / or" as used herein includes any and all combinations of one or more of the associated listed items.

[0027] Example 1: In conventional acoustic emission monitoring of rock fractures, the acquired signals are often mixed with environmental noise and wake interference, and the signal quality of different sensor channels varies. Traditional methods process each stage independently, which can easily lead to noise reduction distortion and misjudgment, thus significantly increasing the positioning error.

[0028] Figure 1 This is a flowchart illustrating the acoustic emission localization method for rock fracture sources based on multi-stage collaborative optimization, according to an embodiment of the present invention.

[0029] like Figure 1 As shown, a method for locating acoustic emission sources of rock fractures based on multi-stage collaborative optimization includes the following steps: S101, collect multi-channel acoustic emission signals generated by rock fracturing, extract characteristic parameters of each channel signal, calculate the Pearson correlation coefficient between channels based on the characteristic parameters, and filter out effective channel signals according to the Pearson correlation coefficient.

[0030] Specifically: A1. Multiple acoustic emission sensors (e.g., Nano30) are non-coplanarly arranged on the rock surface to form a multi-channel acoustic emission sensor array. The spatial coordinates of each sensor are predetermined. The multi-channel acoustic emission full waveform signal generated by rock fracturing is continuously collected at a preset sampling rate (2MHz in this embodiment).

[0031] It should be noted that the trigger threshold can be set as follows: when there is no rupture event, continuously collect 0.5 seconds of silent signal as background noise samples and calculate its root mean square value in the time domain. Set trigger threshold When the amplitude of any channel signal exceeds Th, the complete waveform recording of the event is initiated to eliminate environmental and equipment-inherent noise.

[0032] A2. For the acoustic emission signal received by each sensor channel, extract the feature parameter vector, including but not limited to: count, duration, rise time, peak frequency, frequency center, and the ratio of rise time to amplitude (RA value).

[0033] Among them, count is defined as the number of pulses exceeding the preset threshold, duration is defined as the time interval from the first time exceeding the preset threshold to the time when it drops below the preset threshold, rise time is defined as the time from the first time exceeding the threshold to the maximum amplitude, peak frequency is defined as the frequency corresponding to the largest amplitude value in the amplitude spectrum, frequency center is defined as the centroid frequency of the spectrum, and RA value is defined as the ratio of rise time to amplitude.

[0034] A3. Normalize the feature parameter vector of each channel and calculate the Pearson correlation coefficient between any two channels i and j:

[0035] Where K is the total number of parameters in the feature parameter vector (K is 6 in this embodiment), and k is the feature parameter index. and These are the k-th eigenvalues ​​of the two channels, and These are the average values ​​of the characteristic parameters of the two channels, respectively.

[0036] A4. For each channel i, calculate its average correlation coefficient with all other channels. , where N is the total number of channels (in this embodiment, N is preferably 8), and the effective channels are determined based on the average correlation coefficient.

[0037] Optionally, the quartiles with an average correlation coefficient greater than a preset empirical threshold and sorted in descending order (rounded up if not integer) are selected as valid channels, and the data of the remaining channels are discarded. In other embodiments, the channels with an average correlation coefficient greater than a preset empirical threshold and sorted in descending order, the top Q channels, are also selected as valid channels. The value of Q is determined according to the positioning scenario, that is, the top 3 channels with an average correlation coefficient greater than the preset empirical threshold (3 for two-dimensional positioning scenarios, Q value is preferably 3 in this embodiment) or the top 4 channels (4 for three-dimensional positioning scenarios) are selected as valid channels, and the data of the remaining channels are discarded. If the number of channels that meet the preset empirical threshold is insufficient, all channels that meet the threshold are selected.

[0038] Among them, the preset empirical threshold Pre-calibration can be performed as follows: Specifically, conduct 10 repeated lead-breaking tests at the same location on the rock surface (using HB pencil leads with a diameter of 0.5 mm and an extension length of 2.5 mm, at a 30° angle to the rock surface). Collect signals from 8 channels, calculate the average correlation coefficient between each channel and the remaining channels in each test, and take the minimum average correlation coefficient of all effective channels (with no significant signal distortion) from the 10 tests. Multiply this by 0.9 as a safety margin. For example, the preset empirical threshold is set to 0.6. It should be noted that if the average correlation coefficient is lower than the preset empirical threshold, the data for that channel is severely distorted and should not be used.

[0039] S102, construct an adaptive wavelet threshold function that is continuously differentiable at the threshold and introduces an adjustable curvature parameter. Using the combined index of signal-to-noise ratio and root mean square error as the fitness, optimize the curvature parameter, threshold parameter and wavelet decomposition level of the wavelet threshold function through particle swarm optimization algorithm. Then, use the optimized parameters to perform wavelet threshold denoising on the effective channel signal to obtain a denoised signal.

[0040] The fitness function of the particle swarm optimization algorithm is a weighted combination of signal-to-noise ratio and root mean square error.

[0041] In some embodiments, S102 specifically includes: S201, Perform discrete wavelet transform on each selected valid channel signal to obtain the wavelet coefficients of each layer. , where j is the scale (number of layers) and k is the time index.

[0042] For example, the Daubechies wavelet basis (order db4) is selected, and the number of decomposition layers L is initially set to 3-8 layers to be optimized.

[0043] S202, Construct a continuously differentiable adaptive wavelet threshold function:

[0044] in, These are the original wavelet coefficients. To process the wavelet coefficients, The threshold parameter to be optimized (using the general threshold form) , Where is the noise standard deviation, and N is the signal length. (It can also be scaled up as needed). For adjustable curvature parameters ( This function is in It is continuous at the point of origin and its derivative is continuous, thus overcoming the discontinuity of the hard threshold and the constant deviation of the soft threshold.

[0045] S203 uses a preset particle swarm optimization algorithm to optimize the threshold parameters. Given the curvature parameter k and the number of decomposition layers L, the optimization objective is set as: maximizing the signal-to-noise ratio (SNR) and minimizing the root mean square error (RMSE). A fitness function is constructed (a weighted combination of SNR and RMSE).

[0046] in, , The weighting coefficient (in this embodiment, we take...) ), Defined as: , The signal after denoising. The original noisy signal, Defined as: ; Using a pre-defined particle swarm optimization algorithm in the parameter space Global optimization is performed in the middleware; as an example, in the particle swarm optimization algorithm, the particle swarm settings are as follows: number of particles 30, maximum number of iterations 25, and the parameter search range is determined based on the following prior knowledge: (lower limit) Retaining too much noise, upper limit (leading to oversmoothing) (When k=0.5, the function approaches the hard threshold; when k=5, it approaches the soft threshold.) (Lower limit of 3 allows for noise separation, upper limit of 8 avoids excessive decomposition); Inertia weight is calculated according to the formula. The value decreases linearly from 0.9 to 0.4, where t is the current iteration number. The acceleration factor is 25. Each iteration updates the particle's velocity and position, recording the individual optimal and global optimal. Based on experiments, the fitness converges after approximately 15 iterations, outputting the optimal parameters. ).

[0047] S204 uses the optimal parameters to perform thresholding on the wavelet coefficients of each layer, and then performs inverse wavelet transform to reconstruct the signal, resulting in a high-quality denoised signal.

[0048] S103, perform complete ensemble empirical mode decomposition on the denoised signal and reconstruct component signals containing high-frequency features, calculate the energy difference within the time window before and after the component signals to lock a narrow time window containing the arrival time of the first wave, and pick the arrival time of the first wave of each channel within the narrow time window based on statistical information criteria.

[0049] The complete set of empirical mode decomposition includes: adding white noise of different amplitudes to the original signal multiple times, performing empirical mode decomposition on the signal after each addition of noise, and averaging the decomposition results to obtain multiple intrinsic mode function components; the reconstruction of the component signal containing high-frequency features specifically involves: selecting the first M high-frequency intrinsic mode function components and superimposing them, where M is an integer greater than or equal to 2 and less than or equal to 4; the energy difference is the difference in signal energy within the preceding and following windows.

[0050] In some embodiments, S103 specifically includes: S301, for noise-reducing signals Perform complete set empirical mode decomposition (CEEMDAN). The specific steps are as follows: To the original signal Add white noise of different amplitudes in group Q (in this embodiment, Q=30, which is the empirically optimal value of the CEEMDAN algorithm; when it is greater than 30, the computational cost increases but the improvement in decomposition accuracy is not significant; when it is less than 20, the mode mixing suppression is insufficient; the noise amplitude is selected as follows). The standard deviation is 0.2 times (this value can achieve a balance between decomposition accuracy and convergence speed). Empirical mode decomposition (EMD) is performed on each group of noisy signals to obtain a set of intrinsic mode functions (IMFs). The average of the Q sets of results is taken to obtain the final IMF components, which are multiple intrinsic mode function components. and residual r(n).

[0051] Therefore, a complete set of empirical modes can effectively suppress mode aliasing and make the physical meaning of each IMF clearer.

[0052] S302, Reconstructing high-frequency characteristic components to form component signals: Calculate the center frequency of each IMF component, and select the top M high-frequency intrinsic mode function components with the highest center frequencies (M is an integer greater than or equal to 2 and less than or equal to 4; in this embodiment, M=3, i.e., the first 3 IMFs, denoted as...). The components are superimposed to obtain the component signals: ,Should It highlights the high-frequency abrupt change characteristics of the first arrival wave and suppresses low-frequency noise and wake waves.

[0053] S303 uses CIMF-corrected energy difference (MED) locking to lock a narrow time window. The window length w depends on the sampling rate. and the expected duration of the first wave Sure: ,in The width of the main peak of the first arrival wave is approximately 25. This embodiment Since the frequency is 2MHz, w=50. The specific operation is as follows: Define the front and rear window lengths w (where w represents 50 sampling points) on CIMF, and calculate the energy difference for each sampling point n: Find the point on the MED curve where it first exceeds the mean. And the maximum point (extreme point) of the MED curve. Define a narrow time window containing the first arrival wave as ,in, ; The narrow window extension is set to 200 sampling points (corresponding to 100μs), which is greater than the typical arrival time difference range of acoustic emission signals in rock mass (50-80μs). This ensures that the first arrival wave falls within the window while avoiding the introduction of noise by making the window too large. The length of this narrow window is usually 200-300 sampling points, which significantly reduces the search range for subsequent accurate picking.

[0054] S304, in a narrow time window Within the system, the arrival time of the first wave is accurately determined using preset statistical information criteria. In this embodiment, the preset statistical information criterion is set as the Akaike Information Criterion (AIC), and calculations are performed point by point:

[0055] in, For window Maximum likelihood estimation variance of internal signal ( , (mean of the signals within this window). For window The maximum likelihood estimate variance within ( , (where C is the sample mean of the signal within the corresponding window), and C is a constant term, typically taken as... (p is the order of the autoregressive model, p=2 in this embodiment), but since C does not change with n and has no effect on the location of the minimum value, it can be omitted in the actual calculation; calculate the AIC value of each sampling point n, and take... The sampling point corresponding to the minimum value is taken as the arrival time of the first wave of that channel. .

[0056] S104, determine at least one reference sensor, and based on the arrival time difference of the first arrival wave between each effective channel and the reference sensor, construct a set of positioning equations using the time difference method and perform inversion to output the spatial coordinates of the rupture source.

[0057] The reference sensor is the first sensor in all valid channels to trigger the first arrival wave.

[0058] In some embodiments, S104 specifically includes: S401, Determine the reference sensor. Among all valid channels, select the sensor with the shortest arrival time of the first arrival wave as the reference sensor (i.e., the first trigger sensor). Let its arrival time be... Spatial coordinates are .

[0059] S402, using the arrival time of the first arrival wave from the reference sensor as a benchmark, calculate the arrival time difference for each of the other valid channels i. .

[0060] S403, let the coordinates of the fracture source be (x, y, z), and the wave velocity be v (pre-calibrated through a lead-breaking test, 3800 m / s in this embodiment), and construct the hyperboloid equation system:

[0061] Where i traverses all valid channels except the reference sensor, and (x, y, z) are the coordinates of the rupture source to be determined. , , Let be the coordinates of the i-th sensor, and v be the equivalent wave velocity of the rock mass. For reference sensor coordinates; it should be noted that the wave velocity v can be pre-calibrated in the following way: with known coordinates on the rock mass surface. A lead breakage test was conducted at the location (HB pencil lead diameter 0.5mm, protrusion length 2.5mm, angle with surface 30°), and the arrival time of sensor i was recorded. According to the formula Calculation, where The trigger time for lead breakage is N, where N is the number of sensors. The calibration result in this embodiment is 3800 m / s.

[0062] S404, the positioning result is optimized using a preset iterative optimization algorithm, the initial position of which is set as the geometric center of the sensor array. Output the final rupture source coordinates.

[0063] The iterative optimization algorithm is the Geiger iterative method, which minimizes the sum of squares of the arrival time residuals using the Levenberg-Marquardt algorithm. Specifically, in each iteration, the nonlinear equations are linearized (using a first-order Taylor expansion approximation), and the correction is solved using the least squares method. To prevent the step size from being too large, the Levenberg-Marquardt algorithm is used to add a damping factor. The iteration stops when the correction amount is less than 0.1 mm or after 20 iterations, and the final fracture source coordinates are output.

[0064] The technical solutions described in the embodiments of this application have at least the following technical effects or advantages: By extracting six-dimensional feature parameters such as count, duration, rise time, peak frequency, frequency center, and RA value from multi-channel acoustic emission signals, the Pearson correlation coefficient between channels is calculated, and the three (two-dimensional) or four (three-dimensional) channels with the highest average correlation are selected as valid inputs. This mechanism can actively identify and eliminate abnormal channel data caused by poor sensor coupling, signal attenuation distortion, or hardware failure, thus avoiding contamination of subsequent processing by inferior signals. To address the issues of function discontinuity (hard thresholding) or constant deviation (soft thresholding) in traditional wavelet soft and hard thresholding methods, an adaptive thresholding function with an adjustable curvature parameter is constructed. This function can continuously transition between hard and soft thresholds by adjusting the curvature parameter k. Simultaneously, the particle swarm optimization (PSO) algorithm is used with a weighted combination of signal-to-noise ratio and root mean square error as the fitness to automatically search for the globally optimal threshold parameter λ, curvature parameter k, and decomposition level L. This frees the denoising process from the limitations of manual parameter tuning and enables an adaptive balance between noise suppression and signal detail preservation based on the characteristics of the signal itself. A four-step joint strategy (CMAIC) is adopted, consisting of "CEEMDAN decomposition → high-frequency component reconstruction → energy difference correction narrow window locking → AIC fine pickup". CEEMDAN decomposition decomposes the signal into multiple intrinsic mode functions, effectively suppressing mode aliasing; high-frequency component reconstruction highlights the abrupt change characteristics of the first arrival wave; energy difference correction (MED) quickly locks the narrow time window containing the arrival time by calculating the energy difference between the front and rear windows, narrowing the search range to 200-300 sampling points; finally, the AIC criterion is applied within the narrow window for fine positioning. This strategy balances anti-interference capability (suppressing noise and wake waves through mode decomposition and high-frequency reconstruction) with pickup accuracy (precise AIC positioning within the narrow window). Using the first trigger sensor as the reference node, the arrival time difference between other sensors and the reference is calculated, a hyperboloid equation system is constructed, and the coordinates of the rupture source are solved. This positioning scheme has a simple structure and high computational efficiency, realizing single-reference node positioning inversion based on the time difference method.

[0065] Example 2: In simulation tests of deep rock masses or strongly disturbed fault zones, acoustic emission signals may drop below -10dB, and some sensor channels may completely fail due to attenuation, dispersion, or hardware malfunction. The fixed parameter selection and single-reference-node positioning strategies described in Example 1 may encounter problems under such extreme conditions, such as generally low correlation coefficients, premature convergence of the particle swarm optimization algorithm, narrow-window positioning failure, or overall positioning offset caused by the reference node's own picking errors.

[0066] Based on Example 1, this embodiment introduces an adaptive threshold adjustment mechanism, particle swarm optimization with a restart strategy, multi-window fusion picking, and multi-reference node weighted localization, aiming to further improve the localization robustness under extreme low signal-to-noise ratio and partial sensor failure conditions, making the method applicable to more demanding engineering scenarios.

[0067] In some embodiments, after S101, the method may further include: If the maximum value of the average Pearson correlation coefficient for all channel combinations If so, then adaptive threshold reduction will be initiated.

[0068] Specifically, the adaptive threshold reduction is set as follows: Adaptive threshold reduction uses step size The arithmetic progression strategy is as follows: Based on a preset empirical threshold, the threshold is decreased step by step (for example, the threshold is set to 0.5, 0.4, and 0.3 in sequence) until at least 3 effective channels are selected; if this is still not possible, the current event is determined to be an "extremely low quality event", triggering the backup sensor (if any) to replace the channel with the largest signal variance (the largest variance usually means the most severe noise or distortion), and recalculating the correlation coefficient of all channels. For the selected effective channels, the energy attenuation factor of each channel signal is further calculated. ,in Let i be the energy of the signal in channel i (i.e., the sum of squares of the signal); if a certain channel's < ( The preset attenuation threshold is set to 0.3. Its physical meaning is: when the channel energy is less than 30% of the strongest channel, the signal-to-noise ratio of the channel is usually lower than 0dB, and the positioning contribution is significantly reduced. This threshold is an empirical value and can be adjusted in the range of 0.2 to 0.4 according to the actual monitoring environment. Then it is marked as a "semi-failed channel" and given a degrading weight in the subsequent positioning inversion (the degrading weight coefficient is set to 0.5, and the weight of the normal channel is 1), instead of being directly eliminated.

[0069] In some embodiments, in S102, an enhancement strategy is introduced to address the premature convergence problem of PSO that may be caused by extremely low signal-to-noise ratio signals. The enhancement strategy is set as follows: Latin hypercube sampling is used during particle swarm initialization to ensure parameter space. Uniform coverage to avoid initial concentration; In each iteration, the population diversity index is calculated: ,in, Let be the position vector of the i-th particle. This is the current globally optimal position. The threshold is 30; when D is less than the diversity threshold and the fitness value does not improve after 5 consecutive iterations (5 consecutive iterations account for approximately 20% of the maximum number of iterations 25), a restart is triggered: the current globally optimal particle is retained, the positions of the remaining 70% of the particles (i.e., 21 particles) are randomly initialized, and the remaining 30% of the particles retain their original positions to maintain convergence; where, the diversity threshold is... The physical meaning is: when the average behavior of all particles in the particle swarm to the optimal position is less than 1% of the parameter space ( Spatial span (k spatial span 4.5, L spatial span 5), which suggests that the group is highly clustered.

[0070] The fitness function is modified to use the overall signal variation as a constraint term. The overall signal variation is defined as the sum of the absolute differences between adjacent sampling points of the denoised signal, and is used to suppress excessive smoothing and spurious oscillations during the denoising process.

[0071] in, , and These are the corresponding weighting coefficients, preferably 1, 0.1, and 0.05. The overall variation of the signal, i.e., the overall variation of the signal after denoising, is defined as the sum of the absolute differences between adjacent sampling points: Introducing the TV term can suppress spurious oscillations and make the denoised signal smoother.

[0072] In some embodiments, in step S103, to address the problem that a single narrow window may not be accurately positioned under extremely low signal-to-noise ratios, the method may further include: simultaneously calculating a corrected energy ratio curve and a long-short time average ratio curve on the component signal; determining multiple candidate narrow time windows based on the energy difference, the corrected energy ratio, and the long-short time average ratio, respectively; and using the finely searched window obtained by fusing the multiple candidate narrow time windows as the final narrow time window. Specifically: B1. After performing CEEMDAN and CIMF reconstruction in S103, three characteristic curves are calculated simultaneously on the CIMF signal: MED curve (same as S103), corrected energy ratio (MER) curve, and short-time average / long-time average ratio (STA / LTA) curve.

[0073] The corrected energy ratio (MER) curve is defined as follows: The logarithm is base 10, with units in dB, to facilitate the identification of energy abrupt change points; the short-term average / long-term average ratio (STA / LTA) curve is defined as follows: short window length 10, long window length 100, calculated... .

[0074] B2. Extract candidate narrow windows from the three curves respectively: first narrow window, second narrow window, and third narrow window: From the MED curve: Take the extreme point and the first time exceeding the mean point Get the first narrow window ; From the MER curve: take the rising edge position that exceeds the preset threshold (the sum of the mean and 0.5 times the standard deviation). and peak position Get the second narrow window ; From the STA / LTA curve: A pre-set STA / LTA threshold of 2.5 is used. This value is based on statistical results from typical granite acoustic emission lead-cutting tests (signal-to-noise ratio above 0dB, mean STA / LTA peak value of 3.2 and standard deviation of 0.4 for the first arrival wave in 100 tests). The mean value minus 1.75 times the standard deviation is used as the detection threshold to ensure detection rate while reducing false alarms. In practical applications, this can be adjusted appropriately according to lithology. The starting point of the interval continuously exceeding the STA / LTA threshold of 2.5 is taken. and the end point The third narrow window ].

[0075] B3. Use the intersection of the three candidate windows as the fine-grained search window. If there is no intersection, then take the union of the three and expand outward by 50 sampling points.

[0076] In this embodiment, the statistical information criteria used to pick the arrival times of the first arrival waves for each channel within the narrow time window include the Akaike information criterion and the Bayesian information criterion. Specifically, the picking process involves calculating the arrival times of the first arrival waves using both the Akaike and Bayesian information criteria, and then taking the median of the two results as the final picking result. B4. Within the refined search window, simultaneously calculate the time using both AIC and Bayesian Information Criterion (BIC):

[0077] and Each is a fine search window and The variance of the internal signal (maximum likelihood estimation, same as the variance in AIC); take the sampling points corresponding to the minimum values ​​of AIC and BIC respectively. Take the median of the two results: As the arrival time of the final first arrival wave, it reduces the outlier risk of a single criterion.

[0078] In some embodiments, to avoid overall positioning offset caused by the picking error of a single reference sensor (such as the first trigger sensor), this embodiment adopts a multi-reference node strategy instead of the original single-reference node scheme. In S104, determining at least one reference sensor may further include: solving the positioning equation system sequentially with each effective channel as a reference sensor to obtain multiple candidate coordinates, and determining the final rupture source coordinates based on the posterior residuals of each candidate coordinate; the determination of the final rupture source coordinates based on the posterior residuals of each candidate coordinate specifically involves: selecting the candidate coordinate with the smallest posterior residual as the final rupture source coordinate, or taking a weighted average of the first few candidate coordinates after ascending order as the final rupture source coordinate, with the weight inversely proportional to the square of the posterior residual. Specifically: C1. Using each valid channel as a reference sensor, sequentially solve the positioning equations in S104 (but do not use weights in this solution; first obtain candidate coordinates) to get multiple candidate coordinates. ,in (M is the number of effective channels); C2. Calculate the posterior residual for each candidate coordinate:

[0079] , That is, the sum of the absolute values ​​of the differences between the observation time difference and the calculated time difference; C3. Select the three candidate coordinates with the smallest posterior residuals (if M<3, then take all of them), and calculate the weighted average as the final rupture source coordinates:

[0080] Where S is the selected set of three indices, and the weights are... (The smaller the residual, the greater the weight). Simultaneously, output the positional ellipsoid: with the weighted average coordinates as the ellipsoid center, and the semi-axis length taken as the weighted standard deviation of the three candidate coordinates from the center; C4. For the “semi-failed channel” marked in Example 1, its weight is multiplied by 0.5 when calculating the residual to reduce its impact on the positioning result.

[0081] The technical solutions described in the embodiments of this application have at least the following technical effects or advantages: When the average Pearson correlation coefficient of all channels is lower than a preset empirical threshold, an arithmetic threshold reduction strategy is initiated to ensure that effective channels can still be selected even in extreme noise environments. If enough channels still cannot be selected after threshold reduction, a backup sensor is triggered to replace the channel with the largest signal variance. At the same time, an energy attenuation factor is introduced to evaluate the effectiveness of each channel: channels with less than 30% of the energy of the strongest channel are marked as "semi-failed channels" and given a degraded weight in subsequent positioning. This solves the failure problem of not being able to select effective channels under low signal-to-noise ratio conditions. To address the premature convergence issue of PSO under extremely low signal-to-noise ratios, Latin hypercube sampling initialization is introduced to ensure uniform coverage of the parameter space. During iteration, the population diversity index D is calculated, and a restart is triggered when D < 0.1 and the fitness has not improved for five consecutive iterations. Simultaneously, the fitness function is modified by adding a total variational (TV) penalty term to suppress excessive denoising. These measures effectively improve the convergence rate of PSO under extremely low signal-to-noise ratios. Not limited to a single MED narrow window, this method simultaneously calculates three characteristic curves: MED, Modified Energy Ratio (MER), and STA / LTA, extracts candidate narrow windows for each, and takes the intersection of the three (or the union expanded by 50 points if there is no intersection) as the fine search window. When calculating the results using both AIC and Bayesian Information Criterion (BIC) within the fine window, the median of the two is taken as the final result. This solves the problems of single narrow window localization failure and high outlier risk under extremely low signal-to-noise ratios. For the selection of reference sensors, the positioning equations are solved independently for each effective channel as a reference sensor to obtain multiple candidate coordinates. The posterior residual of each candidate is calculated, and the three candidate coordinates with the smallest residuals are selected for weighted average (the weight is inversely proportional to the square of the residual) as the final result. At the same time, the positioning confidence ellipsoid is output. For the half-failed channel, its weight is multiplied by 0.5 in the residual calculation. This solves the problem of overall positioning offset caused by the picking error of the first trigger sensor itself, and provides users with positioning uncertainty assessment through the confidence ellipsoid.

[0082] In summary, this application addresses common problems in rock fracture acoustic emission localization, such as low signal fidelity, poor first-arrival pickup accuracy, and significant interference from anomalous channels, by constructing a multi-stage collaborative architecture (sensor optimization, adaptive denoising, joint pickup, and localization inversion). Specifically, the channel optimization mechanism based on Pearson correlation coefficient actively eliminates distorted or faulty signals, ensuring data quality from the source; the continuous wavelet threshold denoising method optimized by particle swarm optimization achieves high signal fidelity at low signal-to-noise ratios and is free of pseudo-Gibbs oscillations; the CMAIC joint pickup strategy significantly improves anti-interference capability and pickup accuracy at first-arrival through deep fusion of mode decomposition, energy difference narrow window locking, and statistical information criteria; and by determining at least one reference sensor and performing localization inversion based on the time-difference method, a stable solution for the fracture source coordinates is achieved.

[0083] Furthermore, this embodiment of the invention also provides a rock fracture source acoustic emission localization system based on multi-stage collaborative optimization.

[0084] Figure 2 This is a structural block diagram of the acoustic emission localization system for rock fracture sources based on multi-stage collaborative optimization, according to an embodiment of the present invention.

[0085] like Figure 2 As shown, the acoustic emission localization system for rock fracture sources based on multi-stage collaborative optimization includes: an extraction module, a denoising module, a decomposition and reconstruction module, and a localization module.

[0086] The extraction module is used to acquire multi-channel acoustic emission signals generated by rock fracturing, extract feature parameters of each channel signal, calculate the Pearson correlation coefficient between channels based on the feature parameters, and filter out effective channel signals according to the Pearson correlation coefficient. The denoising module is used to construct an adaptive wavelet threshold function that is continuously differentiable at a threshold and introduces an adjustable curvature parameter. Using a combination of signal-to-noise ratio and root mean square error as the fitness index, the curvature parameter, threshold parameter, and wavelet decomposition level of the wavelet threshold function are optimized using a particle swarm optimization algorithm. The optimized parameters are then used to denoise the effective channel signals. The wave threshold denoising method is used to obtain a denoised signal. The decomposition and reconstruction module is used to perform complete ensemble empirical mode decomposition on the denoised signal and reconstruct the component signal containing high-frequency features. The energy difference within the time window before and after the component signal is calculated to lock a narrow time window containing the arrival time of the first wave. The arrival time of the first wave of each channel is picked up within the narrow time window based on statistical information criteria. The positioning module is used to determine at least one reference sensor. Based on the arrival time difference of the first wave of each effective channel and the reference sensor, the positioning equation set is constructed using the time difference method and inverted to output the spatial coordinates of the rupture source.

[0087] It should be noted that other specific implementation details of the embodiments of the present invention can refer to the above-described method for locating rock fracture sources based on multi-stage collaborative optimization.

[0088] The above description is merely a preferred embodiment of the present invention and is not intended to limit the invention. For those skilled in the art, the present invention can have various modifications and variations. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for locating acoustic emission sources of rock fractures based on multi-stage collaborative optimization, characterized in that, include: S101: Collect multi-channel acoustic emission signals generated by rock fracture, extract characteristic parameters of each channel signal, calculate the Pearson correlation coefficient between channels based on the characteristic parameters, and select effective channel signals according to the Pearson correlation coefficient. S102, construct an adaptive wavelet threshold function that is continuously differentiable at the threshold and introduces an adjustable curvature parameter. Use the comprehensive index of signal-to-noise ratio and root mean square error as the fitness and the overall variation of the signal as the constraint term. Optimize the curvature parameter, threshold parameter and wavelet decomposition level of the wavelet threshold function through particle swarm optimization algorithm. Use the optimized parameters to perform wavelet threshold denoising on the effective channel signal to obtain the denoised signal. S103, perform complete set empirical mode decomposition on the denoised signal and reconstruct the component signal containing high-frequency features. Calculate the energy difference within the time window before and after the component signal to lock a narrow time window containing the arrival time of the first wave. Pick the arrival time of the first wave of each channel within the narrow time window based on statistical information criteria. S104, determine at least one reference sensor, and based on the arrival time difference of the first arrival wave between each effective channel and the reference sensor, construct a set of positioning equations using the time difference method and perform inversion to output the spatial coordinates of the rupture source.

2. The method for locating rock fracture sources based on multi-stage collaborative optimization according to claim 1, characterized in that, In step S104, determining at least one reference sensor is configured as follows: among all effective channels, the sensor with the shortest first arrival time is selected as the reference sensor.

3. The method for locating rock fracture sources based on multi-stage collaborative optimization according to claim 2, characterized in that, The characteristic parameters include at least one or more of the following: count, duration, rise time, peak frequency, frequency center, and the ratio of rise time to amplitude.

4. The method for locating rock fracture sources based on multi-stage collaborative optimization according to claim 3, characterized in that, The process of selecting effective channel signals based on the Pearson correlation coefficient specifically includes: For each channel i, calculate its average correlation coefficient with all other channels, select the channels with an average correlation coefficient greater than a preset empirical threshold and the top Q channels in descending order as valid channels, with the value of Q determined according to the positioning scenario; the data of the remaining channels are discarded.

5. The method for locating rock fracture sources based on multi-stage collaborative optimization according to claim 4, characterized in that, The adaptive wavelet threshold function described in step S102 is expressed as follows: ， in, These are the original wavelet coefficients. These are the wavelet coefficients after thresholding. For threshold parameters, This is an adjustable curvature parameter.

6. The method for locating rock fracture sources based on multi-stage collaborative optimization according to claim 5, characterized in that, The complete set empirical mode decomposition mentioned in step S103 includes: White noise of different amplitudes is added to the original signal multiple times. Empirical mode decomposition is performed on the signal after each addition of noise, and the decomposition results are averaged to obtain multiple intrinsic mode function components.

7. The method for locating rock fracture sources based on multi-stage collaborative optimization according to claim 6, characterized in that, The reconstructed component signal containing high-frequency features specifically includes: The first M high-frequency intrinsic mode function components are selected and superimposed, where M is an integer greater than or equal to 2 and less than or equal to 4; Step S103, which involves locking a narrow time window that includes the arrival time of the first arrival wave, also includes: Simultaneously calculate the corrected energy ratio curve and the long-short time average ratio curve on the component signal. Determine multiple candidate narrow time windows based on the energy difference, the corrected energy ratio, and the long-short time average ratio, respectively. Use the fine search window obtained by fusing the multiple candidate narrow time windows as the final narrow time window.

8. The method for locating rock fracture sources based on multi-stage collaborative optimization according to claim 7, characterized in that, In step S104, determining at least one reference sensor may further include: The positioning equations are solved sequentially using each effective channel as a reference sensor to obtain multiple candidate coordinates. The final rupture source coordinates are determined based on the posterior residuals of each candidate coordinate.

9. The method for locating rock fracture sources based on multi-stage collaborative optimization according to claim 8, characterized in that, The determination of the final fracture source coordinates based on the posterior residuals of each candidate coordinate is specifically as follows: Based on the posterior residuals, several candidate coordinates are selected and weighted averaged to obtain the final rupture source coordinates: ， Where S is the selected set of candidate coordinate indices, and the weights of each index are inversely proportional to the square of the posterior residuals: , The coordinates of the final rupture source, These are candidate coordinates.

10. A rock fracture source acoustic emission localization system based on multi-stage collaborative optimization, applied to the rock fracture source acoustic emission localization method based on multi-stage collaborative optimization as described in any one of claims 1 to 9, characterized in that, The system includes: an extraction module, a noise reduction module, a decomposition and reconstruction module, and a localization module; The extraction module is used to collect multi-channel acoustic emission signals generated by rock fracturing, extract the characteristic parameters of each channel signal, calculate the Pearson correlation coefficient between channels based on the characteristic parameters, and filter out the effective channel signals according to the Pearson correlation coefficient. The denoising module is used to construct an adaptive wavelet threshold function that is continuously differentiable at the threshold and introduces an adjustable curvature parameter. The fitness is based on a comprehensive index of signal-to-noise ratio and root mean square error, and the overall variation of the signal is used as a constraint term. The curvature parameter, threshold parameter and wavelet decomposition level of the wavelet threshold function are optimized by particle swarm optimization algorithm. The optimized parameters are then used to perform wavelet threshold denoising on the effective channel signal to obtain a denoised signal. The decomposition and reconstruction module is used to perform complete set empirical mode decomposition on the denoised signal and reconstruct the component signal containing high-frequency features, calculate the energy difference within the time window before and after the component signal to lock a narrow time window containing the arrival time of the first wave, and pick the arrival time of the first wave of each channel within the narrow time window based on statistical information criteria. The positioning module is used to determine at least one reference sensor, and based on the arrival time difference of the first wave between each effective channel and the reference sensor, construct a set of positioning equations using the time difference method and perform inversion to output the spatial coordinates of the rupture source.

Images