Early disaster tracing method and system based on multi-scale feature fusion

Abstract

The application discloses a burst disaster early-stage tracing method and system based on multi-scale feature fusion, and relates to the field of disaster measurement. The method comprises the following steps: acquiring high-frequency signals and low-frequency signals; obtaining a high-frequency energy wave packet feature flow based on the high-frequency signals, obtaining a low-frequency band energy feature flow based on the low-frequency signals, and calculating the high-frequency propagation time and the low-frequency propagation time to each sensor; performing reverse translation compensation on the high-frequency energy wave packet feature flow and the low-frequency band energy feature flow on a time axis to obtain a high-frequency feature flow and a low-frequency feature flow; extracting burst peak features in the high-frequency feature flow, constructing a sliding time window, and projecting the low-frequency feature flow into the sliding time window; determining a residual time delay, fine-tuning the low-frequency feature flow according to the residual time delay, and obtaining multi-scale features; fusing the multi-scale features to obtain a fused feature map, and determining a disaster early-stage source position according to the fused feature map. The application can effectively improve the accuracy of the tracing direction judgment.

Description

Technical Field

[0001] This application relates to the field of disaster assessment, and in particular to a method and system for early source tracing of sudden disasters based on multi-scale feature fusion. Background Technology

[0002] With the acceleration of urbanization and the continuous expansion of infrastructure, early warning and source tracing of sudden disasters (such as landslides, underground pipeline leaks, and tunnel collapses) have become major challenges in the field of public safety. Disaster monitoring systems typically employ multiple types of sensors to monitor the area in real time, such as acoustic emission sensors and microseismic sensors. Through comprehensive analysis of these sensor signals, anomalies can be identified and the source location traced in the earliest stages of a disaster, buying valuable time for emergency response.

[0003] In existing technologies, early disaster source tracing methods mainly rely on single-scale signal analysis or simple multi-sensor data overlay. A typical approach is to first collect signal data from various sensors, and then estimate the source location using time-of-arrival (TOA) algorithms. These methods usually assume that signals at different scales remain synchronized during propagation or have only a fixed time delay. Therefore, after feature extraction from the sensor signals, they are directly fused, and the source location is inferred from the fusion result.

[0004] However, in the earliest stages of disaster development, weak anomalous signals at the source propagate outward through non-uniform media (such as rock strata, soil, and pipes). Due to the non-uniformity and dispersion characteristics of the media, sensor signals attenuate extremely rapidly and the propagation path is complex. This leads to asynchronous response times of the same source signal at different scales. Simultaneously, as the propagation distance increases, the causal signal at the source is interfered with and diluted by noise along the way, further degrading the signal quality. This results in misjudgments of the source tracing direction, severely impacting the accuracy and reliability of early disaster warnings. Summary of the Invention

[0005] This application provides a method and system for early source tracing of sudden disasters based on multi-scale feature fusion, which can effectively improve the accuracy of source direction judgment and ensure the accuracy and reliability of early disaster warning.

[0006] To achieve the above objectives, the embodiments of this application adopt the following technical solutions: Firstly, a method for early source tracing of sudden disasters based on multi-scale feature fusion is provided, which is applied to a disaster monitoring system. The disaster monitoring system includes an acoustic emission sensor array, a microseismic sensor array, and a data processing host. The acoustic emission sensor array and the microseismic sensor array are respectively communicatively connected to the data processing host. The method includes: Acquire the high-frequency signal output by the acoustic emission sensor array and the low-frequency signal output by the micro-vibration sensor array in the monitoring area; Short-time Fourier transform is performed on the high-frequency signal to obtain the high-frequency energy wave packet characteristic flow, and wavelet packet decomposition is performed on the low-frequency signal to obtain the low-frequency band energy characteristic flow. Based on a pre-built three-dimensional spatial database, spatial grid search points are set for high-frequency energy wave packet characteristic flow and low-frequency band energy characteristic flow, and the reverse acoustic wave propagation path from the spatial grid search points to each sensor, as well as the high-frequency propagation time and the low-frequency propagation time, are calculated. Based on the high-frequency propagation time and the low-frequency propagation time, the high-frequency energy wave packet characteristic flow and the low-frequency band energy characteristic flow are reverse-shifted and compensated on the time axis to obtain the spatiotemporally aligned high-frequency characteristic flow and the spatiotemporally aligned low-frequency characteristic flow. Extract burst spike features from the high-frequency feature stream, construct a sliding time window centered on the timestamp of the burst spike features, and project the low-frequency feature stream into the sliding time window. Calculate the normalized cross-correlation function of the high-frequency feature flow and the low-frequency feature flow, determine the residual time delay based on the peak position of the normalized cross-correlation function, and fine-tune the low-frequency feature flow based on the residual time delay to obtain multi-scale features; Multi-scale features are fused to obtain a fused feature map, and the location of the early source of the disaster is determined based on the fused feature map.

[0007] In one possible implementation of the first aspect, the calculation of the reverse acoustic wave propagation path from the spatial grid search point to each sensor, the high-frequency propagation time, and the low-frequency propagation time includes: Obtain the equivalent refractive index at each position on the path between the spatial grid search point and the sensor from the three-dimensional spatial database; Calculate the propagation path of sound waves in a non-homogeneous medium based on the equivalent refractive index; The propagation path is divided into multiple micro-segments, and the high-frequency equivalent propagation velocity and low-frequency equivalent propagation velocity corresponding to each micro-segment are obtained. Calculate the high-frequency propagation time and low-frequency propagation time of each micro-element, where the high-frequency propagation time is equal to the length of the micro-element divided by the high-frequency equivalent propagation speed, and the low-frequency propagation time is equal to the length of the micro-element divided by the low-frequency equivalent propagation speed. The high-frequency propagation time is obtained by summing the high-frequency propagation time of all micro-segments, and the low-frequency propagation time is obtained by summing the low-frequency propagation time of all micro-segments.

[0008] In another possible implementation of the first aspect, burst spike features are extracted from the high-frequency feature stream, including: Calculate the first-order time derivative of the high-frequency characteristic flow; Calculate the second time derivative of the first time derivative; Determine whether the second-order time derivative exceeds a preset mutation threshold; When the second time derivative exceeds a preset mutation threshold, the feature point at the corresponding time moment is determined to be a sudden spike feature.

[0009] In another possible implementation of the first aspect, the normalized cross-correlation function of the high-frequency feature stream and the low-frequency feature stream is calculated, and the residual time delay is determined based on the peak position of the normalized cross-correlation function, including: Obtain spatiotemporally aligned high-frequency feature stream segments and spatiotemporally aligned low-frequency feature stream segments within a sliding time window; Normalize the high-frequency characteristic flow segments and the low-frequency characteristic flow segments respectively; Calculate the cross-correlation function between the normalized high-frequency characteristic flow segment and the normalized low-frequency characteristic flow segment, and determine the peak position of the cross-correlation function; The residual time delay of the spatiotemporally aligned low-frequency feature stream segment relative to the spatiotemporally aligned high-frequency feature stream segment is determined based on the peak position.

[0010] In another possible implementation of the first aspect, multi-scale features are fused to obtain a fused feature map, including: Obtain the spatial coordinates of each sensor in the acoustic emission sensor array and the micro-vibration sensor array; A graph structure is constructed based on spatial coordinates, where multi-scale features serve as the initial vectors for the graph nodes in the graph structure. The attention coefficients between adjacent graph nodes are calculated based on the graph attention network, and the gradient of the attention coefficient of each graph node in space is calculated. The spatial receptive field radius of the graph node is adjusted according to the gradient. Multi-scale features are fused based on the adjusted spatial receptive field radius to output a fused feature map.

[0011] In another possible implementation of the first aspect, adjusting the spatial receptive field radius of the graph nodes according to the gradient includes: Calculate the magnitude of the gradient of the attention coefficient; Determine whether the magnitude of the gradient of the attention coefficient exceeds a preset gradient threshold; If the magnitude of the gradient of the attention coefficient exceeds the preset gradient threshold, the graph node is determined to be in the region of increased feature energy density, and the spatial receptive field radius of the graph node is shrunk to the preset minimum radius. If the magnitude of the gradient of the attention coefficient does not exceed the preset gradient threshold, the spatial receptive field radius of the graph node is kept at the preset standard radius.

[0012] In another possible implementation of the first aspect, multi-scale features are fused based on the adjusted spatial receptive field radius to output a fused feature map, including: Using the spatial coordinates of each graph node as the center and the corresponding adjusted spatial receptive field radius as the boundary, the spatial action sphere of each graph node is divided in three-dimensional space. Extract the high-frequency acoustic emission features and low-frequency microseismic features corresponding to each graph node, and map them into a virtual acoustic pressure scalar field and a virtual work vector field in the spatial action sphere, respectively. Within the overlapping region of the spatial action sphere of adjacent graph nodes, gradient divergence cross-coupling operation is performed on the virtual sound pressure scalar field and the virtual work vector field to generate a local spatiotemporal coupling tensor. Calculate the characteristic norm of the local spatiotemporal coupling tensor in three-dimensional space to quantitatively characterize the physical interaction strength of features at different scales in the overlapping space; The feature norm is interpolated radially along the spatial geometric path of the overlapping region to generate a continuous feature energy distribution field. Bilateral filtering is applied to the characteristic energy distribution field to output a fused feature map.

[0013] In another possible implementation of the first aspect, determining the early source location of a disaster based on the fused feature map includes: The fused feature map is mapped onto a spatial grid of the monitoring area, and the feature energy value of each spatial grid point is calculated. Calculate the gradient vector of the feature energy value in three-dimensional space and determine the temporal evolution direction of the multi-scale features; The causal flow direction of each spatial grid point is determined based on the temporal evolution direction; The gradient vector is combined with the causal flow direction to obtain the comprehensive gradient. Perform a reverse search in the opposite direction of the comprehensive gradient to determine the geometric center, and output the physical coordinates of the geometric center as the location of the early source of the disaster.

[0014] Secondly, this application provides a data processing host, comprising: The memory is configured to store instructions; and The processor is configured to retrieve the instructions from the memory and, when executing the instructions, to implement the aforementioned method for early source tracing of sudden disasters based on multi-scale feature fusion.

[0015] Thirdly, this application provides a disaster monitoring system, comprising: Acoustic emission sensor array; Micro-vibration sensor array; The data processing host is communicatively connected to both the acoustic emission sensor array and the micro-vibration sensor array.

[0016] By integrating acoustic emission sensor arrays and microseismic sensor arrays into the disaster monitoring system and employing a multi-scale feature fusion method, the system achieves accurate source tracing of early-stage sudden disasters, effectively overcoming the limitations of traditional single-scale signal analysis methods. Through short-time Fourier transform and wavelet packet decomposition of high-frequency and low-frequency signals respectively, energy feature flows at different scales are extracted, fully utilizing the complementary information of multi-scale signals and improving the comprehensiveness of feature representation. Based on a pre-constructed three-dimensional spatial database, the propagation paths and propagation times of different frequency signals in non-uniform media are accurately calculated, and inverse translation compensation is performed, effectively solving the time asynchrony problem caused by medium non-uniformity and dispersion characteristics, and achieving preliminary spatiotemporal alignment of multi-scale features. By extracting sudden spike features and calculating normalized cross-correlation functions, residual time delays are further eliminated, achieving precise time synchronization between high-frequency and low-frequency features, significantly improving the accuracy of feature alignment. Furthermore, a graph attention network is used for spatial fusion of multi-scale features, and the spatial receptive field radius is adaptively adjusted according to the attention coefficient gradient, effectively enhancing the spatial resolution of the source region and improving the targeting of feature fusion. Finally, by comprehensively considering spatial energy gradients and temporal causal flow, a reverse search was performed along the opposite direction of the comprehensive gradient, achieving precise location of the early source of the disaster. This technical solution is not only applicable to rockburst monitoring in underground mines, but can also be extended to early warning and source tracing of various sudden disasters such as landslides, underground pipeline leaks, and tunnel collapses, providing strong technical support for disaster prevention and control in the field of public safety and effectively improving the accuracy and reliability of early disaster warnings.

[0017] Other features and advantages of the embodiments of this application will be described in detail in the following detailed description section. Attached Figure Description

[0018] Figure 1 A flowchart illustrating an early source tracing method for sudden disasters based on multi-scale feature fusion, provided as an embodiment of this application; Figure 2 This is a schematic diagram of the structure of a disaster monitoring system provided in an embodiment of this application; Figure 3 This is a schematic diagram illustrating a process for dividing the spatial action sphere of each graph node in three-dimensional space, as provided in an embodiment of this application. Figure 4 This is a schematic diagram of a process for determining the early source location of a disaster based on a fused feature map, as provided in an embodiment of this application. Detailed Implementation

[0019] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. It should be understood that the specific embodiments described herein are only for illustration and explanation of the embodiments of this application and are not intended to limit the embodiments of this application. All other embodiments obtained by those skilled in the art based on the embodiments of this application without creative effort are within the scope of protection of this application.

[0020] It should be noted that if the embodiments of this application involve directional indicators (such as up, down, left, right, front, back, etc.), the directional indicators are only used to explain the relative positional relationship and movement of each component in a certain specific posture (as shown in the figure). If the specific posture changes, the directional indicators will also change accordingly.

[0021] Furthermore, if the embodiments of this application involve descriptions such as "first" or "second," these descriptions are for descriptive purposes only and should not be construed as indicating or implying their relative importance or implicitly specifying the number of technical features indicated. Therefore, features defined with "first" or "second" may explicitly or implicitly include at least one of those features. Additionally, the technical solutions of various embodiments can be combined with each other, but this must be based on the ability of those skilled in the art to implement them. If the combination of technical solutions is contradictory or impossible to implement, it should be considered that such a combination of technical solutions does not exist and is not within the scope of protection claimed in this application.

[0022] Figure 1 The illustration shows a flowchart of an early source tracing method for sudden disasters based on multi-scale feature fusion according to an embodiment of this application. Figure 1 As shown in the figure, this application provides an early source tracing method for sudden disasters based on multi-scale feature fusion, which is applied to a disaster monitoring system. The disaster monitoring system includes an acoustic emission sensor array, a microseismic sensor array, and a data processing host. The acoustic emission sensor array and the microseismic sensor array are respectively communicatively connected to the data processing host. The method may include the following steps.

[0023] S110. Acquire the high-frequency signal output by the acoustic emission sensor array and the low-frequency signal output by the micro-vibration sensor array in the monitoring area; S120. Perform short-time Fourier transform on the high-frequency signal to obtain the high-frequency energy wave packet characteristic flow, and perform wavelet packet decomposition on the low-frequency signal to obtain the low-frequency band energy characteristic flow. S130. Based on a pre-built three-dimensional spatial database, spatial grid search points are set for high-frequency energy wave packet characteristic flow and low-frequency band energy characteristic flow, and the reverse acoustic wave propagation path from the spatial grid search points to each sensor, as well as the high-frequency propagation time and the low-frequency propagation time are calculated. S140. Based on the high-frequency propagation time and the low-frequency propagation time, perform reverse translation compensation on the high-frequency energy wave packet characteristic flow and the low-frequency band energy characteristic flow on the time axis to obtain the spatiotemporally aligned high-frequency characteristic flow and the spatiotemporally aligned low-frequency characteristic flow. S150. Extract sudden spike features from the high-frequency feature stream, construct a sliding time window centered on the timestamp of the sudden spike features, and project the low-frequency feature stream into the sliding time window. S160. Calculate the normalized cross-correlation function of the high-frequency feature flow and the low-frequency feature flow, and determine the residual time delay based on the peak position of the normalized cross-correlation function. Then, fine-tune and shift the low-frequency feature flow based on the residual time delay to obtain multi-scale features. S170. The multi-scale features are fused to obtain a fused feature map, and the location of the early source of the disaster is determined based on the fused feature map.

[0024] In this embodiment, the data processing host can be a device with a processor, such as a tablet computer, desktop computer, laptop computer, handheld computer, wearable device, laptop computer, ultra-mobile personal computer (UMPC), or netbook. Of course, the data processing host can also be a server. This application embodiment does not impose any special limitations on the specific form of the data processing host.

[0025] In this embodiment, the method for constructing a three-dimensional spatial database may include the following steps: setting multiple active seismic source emission points at known locations within the monitoring area; controlling the active seismic source emission points to sequentially emit test signals; obtaining the arrival time of the test signals received by the sensor array; calculating the sound wave propagation speed of the medium in the monitoring area at different frequency bands based on the location of the active seismic source emission points, the sensor location, and the arrival time; mapping the sound wave propagation speed to an equivalent refractive index; and storing the equivalent refractive index according to three-dimensional spatial coordinates to form a three-dimensional spatial database.

[0026] In one embodiment of this example, determining the geometric center by performing a reverse search along the opposite direction of the comprehensive gradient may include: taking the spatial grid point with the highest feature energy value in the fused feature map as the search starting point; starting from the search starting point, performing an iterative search along the opposite direction of the comprehensive gradient; in each iteration, moving to the spatial grid point with the largest reverse component of the comprehensive gradient among adjacent spatial grid points; determining whether the magnitude of the comprehensive gradient of the current spatial grid point is less than a preset convergence threshold; if the magnitude of the comprehensive gradient is less than the preset convergence threshold, determining the current spatial grid point as the geometric center; if the magnitude of the comprehensive gradient is greater than or equal to the preset convergence threshold, continuing the iterative search until the magnitude of the comprehensive gradient is less than the preset convergence threshold.

[0027] Specifically, determining the geometric center by performing a reverse search along the opposite direction of the comprehensive gradient can include: taking the spatial grid point with the highest feature energy value in the fused feature map as the search starting point; starting from the search starting point, performing an iterative search along the opposite direction of the comprehensive gradient; in each iteration, moving to the spatial grid point with the largest component of the comprehensive gradient in the opposite direction among adjacent spatial grid points; determining whether the magnitude of the comprehensive gradient of the current spatial grid point is less than a preset convergence threshold; if the magnitude of the comprehensive gradient is less than the preset convergence threshold, determining the current spatial grid point as the geometric center; if the magnitude of the comprehensive gradient is greater than or equal to the preset convergence threshold, continuing the iterative search until the magnitude of the comprehensive gradient is less than the preset convergence threshold.

[0028] In each iteration, the process moves to the spatial grid point with the largest reverse component of the combined gradient among adjacent spatial grid points. This includes: calculating the second-order partial derivative of the eigenenergy value at the current spatial grid point in the current iteration step to construct a target matrix describing the local curvature; extracting the spatial local curvature eigenvalues ​​of the current grid region by performing eigenvalue decomposition on the target matrix; determining whether the current grid region is at a saddle point or a flat region caused by local extrema based on the spatial local curvature eigenvalues; if it is determined to be at a saddle point or a flat region, introducing the search direction and historical displacement of the previous iteration step to calculate the momentum correction term used to maintain search inertia; combining the combined gradient of the current step with the momentum correction term to obtain a corrected search direction vector used to eliminate local extrema interference; adjusting the current iteration search step size inversely proportionally based on the largest eigenvalue of the target matrix, so that a micro-step search is used in high curvature regions to prevent out-of-bounds errors, and a large-step search is used in low curvature regions to improve the convergence rate; and using the conjugate gradient direction instead of the original reverse component of the combined gradient for spatial coordinate iterative updates to eliminate sawtooth oscillating paths and determine the next moving grid point.

[0029] In this embodiment, the monitoring area can be an underground mining area, the acoustic emission sensor array can be deployed on the roadway perimeter, the microseismic sensor array can be deployed around the goaf, and the early source location of the disaster is the early lesion location of rockburst. Figure 2 This application provides a schematic diagram of the structure of a disaster monitoring system according to an embodiment of the present application. (Refer to...) Figure 2 The disaster monitoring system includes an acoustic emission sensor array, a microseismic sensor array, and a data processing host. The acoustic emission sensor array and the microseismic sensor array are communicatively connected to the data processing host. The acoustic emission sensor array is used to capture high-frequency signals, the microseismic sensor array is used to capture low-frequency signals, and the data processing host is used to execute the data processing flow from steps S110 to S170.

[0030] In practice, the data processing host first acquires monitoring data in real time from the acoustic emission sensor array and the microseismic sensor array. The sampling frequency of the acoustic emission sensor array is typically set in the range of 100kHz to 1MHz to capture high-frequency elastic wave signals generated by the propagation of microcracks in the rock mass. The frequency range of these signals is usually between 20kHz and 500kHz, which can reflect the earliest microscopic fracturing process in the formation of a disaster. The sampling frequency of the microseismic sensor array is typically set in the range of 1kHz to 10kHz to capture low-frequency vibration signals generated by stress adjustment in the rock mass. The frequency range of these signals is usually between 1Hz and 1000Hz, which can reflect a larger range of stress field changes.

[0031] After acquiring the raw monitoring data, time-frequency analysis needs to be performed on the high-frequency and low-frequency signals separately to extract feature flows at different scales. For the high-frequency signals, a short-time Fourier transform is used for processing. In practice, the time window length is set to 512 sampling points, the window overlap rate is 75%, and the Hanning window function is used to reduce spectral leakage.

[0032] The short-time Fourier transform converts the time-domain signal into a time-frequency domain representation, calculating the spectral energy distribution within each time window to obtain a time-varying spectral energy sequence, i.e., the high-frequency energy wave packet feature flow. This feature flow clearly reveals the energy abrupt change process of the high-frequency signal, reflecting the instantaneous propagation behavior of microcracks. For low-frequency signals, wavelet packet decomposition is used, selecting the db4 wavelet basis function and setting the decomposition level to 5, decomposing the signal into 32 frequency band sub-bands. The energy value of each frequency band sub-band in each time period is calculated, obtaining a matrix of multi-band energy changes over time, i.e., the low-frequency band energy feature flow. This feature flow can finely characterize the energy distribution features of low-frequency signals in different frequency bands, reflecting the multi-scale dynamic process of stress field adjustment. Through this step, the original time-domain signal is converted into a time-frequency domain feature representation, effectively extracting the energy evolution features of signals at different scales, providing a feature basis for subsequent spatiotemporal alignment processing.

[0033] After extracting the feature flow, the complexity of signal propagation in a non-uniform medium needs to be considered, and the propagation path and propagation time from the hypothetical source point to each sensor need to be calculated. First, spatial grid search points are set in the three-dimensional space of the monitoring area, with a grid spacing typically set to 5 to 10 meters, forming a set of candidate source points covering the entire monitoring area. For each spatial grid search point, the medium parameters along the path between that point and each sensor are queried from a pre-built three-dimensional spatial database.

[0034] A three-dimensional spatial database stores the equivalent refractive index of each spatial location within the monitoring area. This refractive index, obtained through active source calibration experiments, reflects the non-homogeneity of the medium. Based on the equivalent refractive index, a ray tracing algorithm is used to calculate the actual propagation path of sound waves in the non-homogeneous medium. This path is typically not a straight line, but rather the shortest time path determined according to Fermat's principle. The propagation path is divided into multiple micro-segments of 0.1 meters in length. For each micro-segment, the high-frequency equivalent propagation velocity at that location is obtained from the database. and low-frequency equivalent propagation speed Due to the dispersion characteristics of the medium, high-frequency and low-frequency signals propagate at different speeds in the same medium, with high-frequency signals typically propagating faster. Calculate the high-frequency propagation time for each micro-segment. and low-frequency propagation time ,in Let be the length of the infinitesimal segment. Sum the propagation times of all infinitesimal segments to obtain the high-frequency propagation time from the spatial grid search point to the sensor. and low-frequency propagation time For example, for a spatial grid search point 50 meters away from a sensor, ray tracing calculations show that the actual propagation path length is 52.3 meters, the high-frequency propagation time is 8.7 milliseconds, and the low-frequency propagation time is 10.2 milliseconds, a difference of 1.5 milliseconds. This step accurately calculates the propagation time difference of different frequency signals in a non-uniform medium, providing a theoretical basis for subsequent spatiotemporal alignment compensation.

[0035] After obtaining the propagation time, the feature stream needs to be spatiotemporally aligned to eliminate the time offset caused by the propagation delay. For each spatial grid search point, the calculated high-frequency propagation time is used as the basis for this process. and low-frequency propagation time Inverse translation compensation is performed on the high-frequency energy wave packet characteristic flow and the low-frequency band energy characteristic flow on the time axis.

[0036] In practical implementation, let's assume the current time is... The high-frequency feature stream received by the sensor is recorded at time... Then the actual time when the signal was emitted from the source should be Therefore, the high-frequency characteristic flow is shifted forward on the time axis. This aligns it with the source point's emission time.

[0037] Similarly, shift the low-frequency characteristic flow forward on the time axis. .because and Typically unequal, after reverse translation, the alignment reference for the high-frequency and low-frequency feature streams on the time axis becomes the assumed source emission time, rather than the sensor reception time. This reverse compensation method effectively eliminates the time delay caused by differences in propagation paths, enabling signals of different scales from the same source to be initially aligned on the time axis. For example, for a certain spatial grid search point, the high-frequency feature stream is shifted forward by 8.7 milliseconds, and the low-frequency feature stream is shifted forward by 10.2 milliseconds. After the shift, the start times of the energy mutations in the two feature streams tend to be consistent.

[0038] It is important to note that due to the complexity of medium parameter estimation errors and dispersion characteristics, inverse translation compensation can only achieve coarse alignment, and residual time delays still exist, requiring subsequent fine correction. This step achieves preliminary spatiotemporal alignment of multi-scale feature flows, creating conditions for subsequent fine alignment and feature fusion.

[0039] After completing the initial spatiotemporal alignment, it is necessary to identify key time markers of disaster events, namely sudden spike features, in the high-frequency feature stream. Sudden spike features correspond to the instantaneous propagation of microcracks in the rock mass or the sudden release of stress, and are key events in the disaster incubation process. The method for identifying sudden spike features is to calculate the time derivative of the high-frequency feature stream.

[0040] Specifically, the first time derivative of the high-frequency characteristic flow is first calculated, reflecting the rate of energy change. Then, the second time derivative of the first derivative is calculated, reflecting the rate of change of the energy change rate, i.e., the acceleration characteristic. When the second derivative exceeds a preset abrupt change threshold, it indicates a sharp energy change, and the corresponding feature point is marked as a sudden spike characteristic. The preset abrupt change threshold is usually determined based on historical data statistics and is set to 5 to 10 times the normal fluctuation level. For example, in a certain monitoring, the second derivative value of the high-frequency characteristic flow at time 15.3 seconds reached 8.2 times the normal level and was identified as a sudden spike characteristic.

[0041] After identifying the sudden spike characteristics, a sliding time window is constructed centered on the timestamp of this characteristic. The window length is typically set to 0.5 to 2 seconds, covering the time range before and after the sudden event. Low-frequency feature streams are projected into this sliding time window, extracting data segments of the low-frequency feature streams within that time period. This step identifies the key time nodes of the disaster event and establishes a time correlation window between high-frequency and low-frequency features, providing a local analytical scope for subsequent fine-grained time alignment.

[0042] After determining the sliding time window, fine-grained time alignment of the high-frequency and low-frequency feature flows within the window is required to eliminate residual delays. First, spatiotemporally aligned high-frequency and low-frequency feature flow segments are extracted within the sliding time window. These two segments are then normalized to eliminate the influence of amplitude differences, ensuring a mean of 0 and a standard deviation of 1. Finally, the cross-correlation function between the normalized high-frequency and low-frequency feature flow segments is calculated.

[0043] The cross-correlation function is defined as follows: ,in For the normalized high-frequency feature stream segment, The normalized low-frequency characteristic flow segment in time offset The value is determined by iterating through different time offsets. Calculate the cross-correlation function at each time offset and determine the peak position of the cross-correlation function. The peak position corresponds to the time offset where the high-frequency feature stream segment and the low-frequency feature stream segment have the highest similarity, i.e., the residual delay. For example, in a certain calculation, the cross-correlation function at... The peak value of 0.87 seconds indicates that the low-frequency characteristic flow segment lags behind the high-frequency characteristic flow segment by 0.15 seconds.

[0044] Based on residual delay Fine-tuning the low-frequency characteristic flow involves shifting it forward or backward along the time axis. This process ensures precise temporal alignment between the high-frequency and low-frequency feature streams. After fine-tuning and shifting, the high-frequency and low-frequency feature streams achieve precise alignment at critical moments of sudden events, with their energy mutation times highly consistent, forming multi-scale features. This step eliminates residual time delay errors after initial alignment, achieving precise time synchronization between high-frequency and low-frequency features and providing high-quality input data for subsequent feature fusion.

[0045] After obtaining precisely aligned multi-scale features, these features need to be spatially fused to generate a comprehensive fused feature map. First, the spatial coordinates of each sensor in the acoustic emission sensor array and the microseismic sensor array are obtained; these coordinates were measured and recorded during sensor deployment. A graph structure is constructed based on the sensor's spatial coordinates, where each node represents a sensor, and the edges between nodes represent the spatial proximity between sensors. The multi-scale features are used as the initial vectors for the graph nodes, meaning each node contains both high-frequency and low-frequency features corresponding to that sensor. A graph attention network is used to calculate the attention coefficients between adjacent graph nodes; these attention coefficients reflect the correlation strength between features from adjacent sensors.

[0046] The calculation formula for graph attention networks is as follows: ,in and They are nodes and nodes eigenvectors, This represents vector concatenation. This is a learnable weight matrix.

[0047] The gradient of the attention coefficient for each graph node in space is calculated. The gradient reflects the rate of change of the attention coefficient in space. The spatial receptive field radius of the graph node is adjusted based on the gradient. Specifically, the magnitude of the attention coefficient gradient is calculated, and it is determined whether this magnitude exceeds a preset gradient threshold. If it exceeds the threshold, it indicates that the node is in a region of increased feature energy density, possibly close to a disaster source. In this case, the spatial receptive field radius of the node is reduced to a preset minimum radius, such as 3 meters, to improve spatial resolution. If it does not exceed the threshold, it indicates that the node is in a region of flat features, and the spatial receptive field radius is maintained at a preset standard radius, such as 10 meters.

[0048] After adjusting the spatial receptive field radius, a spatial action sphere is delineated in three-dimensional space, centered on the spatial coordinates of each graph node and bounded by the corresponding adjusted radius. High-frequency acoustic emission features and low-frequency microseismic features corresponding to each graph node are extracted and mapped to a virtual sound pressure scalar field and a virtual work vector field within the spatial action sphere, respectively. Gradient divergence cross-coupling operations are performed on the virtual sound pressure scalar field and the virtual work vector field within the overlapping region of the spatial action spheres of adjacent graph nodes to generate a local spatiotemporal coupling tensor. The feature norm of the local spatiotemporal coupling tensor in three-dimensional space is calculated; this norm quantitatively characterizes the physical interaction intensity of features at different scales within the overlapping space. Radial basis function interpolation is performed on the feature norm along the spatial geometric path of the overlapping region to generate a continuous feature energy distribution field. Bilateral filtering is applied to the feature energy distribution field to smooth noise while preserving edge features, outputting a fused feature map. This step achieves deep fusion of multi-scale features in the spatial dimension, generating a fused feature map containing rich spatiotemporal information, providing a high-quality feature representation for final source localization.

[0049] After obtaining the fused feature map, the early source location of the disaster can be determined based on this map. First, the fused feature map is mapped onto a spatial grid of the monitoring area. Each spatial grid point corresponds to a feature energy value, which reflects the probability that the location is a disaster source. The gradient vector of the feature energy value in three-dimensional space is calculated; the gradient vector points in the direction of the fastest growth of the feature energy. Simultaneously, based on the temporal evolution characteristics of the multi-scale features, the temporal evolution direction of the features is determined, reflecting the causal relationship of energy propagation. The causal flow direction of each spatial grid point is determined based on the temporal evolution direction, indicating the direction of energy propagation outward from the source.

[0050] The gradient vector is synthesized with the causal flow direction to obtain the comprehensive gradient, which considers both spatial energy distribution and temporal causality. A reverse search is then performed along the opposite direction of the comprehensive gradient, i.e., searching from regions of high feature energy to regions of low energy, gradually approaching the origin of energy. The reverse search process employs an iterative algorithm. First, the spatial grid point with the highest feature energy value in the fused feature map is used as the search starting point. Then, starting from the search starting point, in each iteration, the search moves to the spatial grid point with the largest reverse component of the comprehensive gradient among adjacent spatial grid points.

[0051] To avoid getting trapped in local optima, the second-order partial derivatives of the eigenenergy values ​​at the current spatial grid point are calculated in the current iteration step to construct a target matrix describing the local curvature. Eigenvalue decomposition is performed on the target matrix to extract spatial local curvature eigenvalues, determining whether the current grid region is at a saddle point or a flat region caused by a local extremum. If it is at a saddle point or flat region, the search direction and historical displacement from previous iterations are introduced to calculate a momentum correction term. The current step's combined gradient and the momentum correction term are vector-synthesized to obtain a corrected search direction vector. Based on the largest eigenvalue of the target matrix, the current iteration search step size is adjusted inversely proportionally, using a micro-step search for high-curvature regions and a large-step search for low-curvature regions. The conjugate gradient direction is used instead of the original combined gradient in the opposite direction for spatial coordinate iterative updates to determine the next moving grid point.

[0052] The algorithm determines whether the comprehensive gradient magnitude of the current spatial grid point is less than a preset convergence threshold. If it is less than the threshold, the current spatial grid point is determined as the geometric center, i.e., the location of the early source of the disaster. If it is greater than or equal to the threshold, the iterative search continues. The final output is the physical coordinates of the geometric center as the location of the early source of the disaster. For example, in tracing the source of a rockburst event, the reverse search converged after 23 iterations, determining the source location as coordinates (125.3 meters, 87.6 meters, -450.2 meters). This location differs from the rockburst lesion location confirmed by subsequent on-site investigation by only 2.8 meters. This step achieves precise location of the early source of the disaster, providing accurate location information for emergency response and disaster prevention.

[0053] This embodiment integrates acoustic emission sensor arrays and microseismic sensor arrays into a disaster monitoring system and employs a multi-scale feature fusion method to achieve accurate source tracing of sudden disasters in their early stages, effectively overcoming the limitations of traditional single-scale signal analysis methods. By performing short-time Fourier transform and wavelet packet decomposition on high-frequency and low-frequency signals respectively, energy feature flows at different scales are extracted, fully utilizing the complementary information of multi-scale signals and improving the comprehensiveness of feature representation. Based on a pre-constructed three-dimensional spatial database, the propagation paths and propagation times of signals of different frequencies in non-uniform media are accurately calculated, and inverse translation compensation is performed, effectively solving the time asynchrony problem caused by medium non-uniformity and dispersion characteristics, and achieving preliminary spatiotemporal alignment of multi-scale features. By extracting sudden spike features and calculating normalized cross-correlation functions, residual time delays are further eliminated, achieving precise time synchronization between high-frequency and low-frequency features, significantly improving the accuracy of feature alignment. In addition, a graph attention network is used to spatially fuse multi-scale features, and the spatial receptive field radius is adaptively adjusted according to the attention coefficient gradient, effectively enhancing the spatial resolution of the source region and improving the targeting of feature fusion. By comprehensively considering spatial energy gradients and temporal causal flow, a reverse search is performed along the opposite direction of the comprehensive gradient. Momentum correction and conjugate gradient optimization are introduced to effectively avoid local extremum interference and achieve precise location of the early source of disasters. This technical solution is not only applicable to rockburst monitoring in underground mines, but can also be extended to early warning and source tracing of various sudden disasters such as landslides, underground pipeline leaks, and tunnel collapses. It provides strong technical support for disaster prevention and control in the field of public safety, effectively improving the accuracy and reliability of early disaster warnings.

[0054] In one embodiment of this invention, calculating the reverse acoustic wave propagation path from the spatial grid search point to each sensor, the high-frequency propagation time, and the low-frequency propagation time includes the following steps: S210. Obtain the equivalent refractive index at each position on the path between the spatial grid search point and the sensor from the three-dimensional spatial database; S220. Calculate the propagation path of sound waves in a non-uniform medium based on the equivalent refractive index; S230. Divide the propagation path into multiple micro-segments and obtain the high-frequency equivalent propagation velocity and low-frequency equivalent propagation velocity corresponding to each micro-segment. S240. Calculate the high-frequency propagation time and low-frequency propagation time of each micro-element segment, where the high-frequency propagation time is equal to the length of the micro-element segment divided by the high-frequency equivalent propagation speed, and the low-frequency propagation time is equal to the length of the micro-element segment divided by the low-frequency equivalent propagation speed. S250. The high-frequency propagation time is obtained by summing the high-frequency propagation time of all micro-segments, and the low-frequency propagation time is obtained by summing the low-frequency propagation time of all micro-segments.

[0055] In practical applications, achieving multi-scale feature spatiotemporal alignment requires accurate calculation of the propagation path and propagation time of sound waves in non-homogeneous media. Traditional methods typically assume a homogeneous medium, employing a straight-line propagation model and a fixed propagation velocity. This simplification leads to significant time delay estimation errors in non-homogeneous geological environments, resulting in inaccurate feature alignment and misjudgment of the source direction. This implementation effectively addresses the impact of non-homogeneous media and dispersion effects on time delay calculations by introducing a three-dimensional spatial database and a frequency-dependent propagation velocity model, providing a reliable theoretical foundation for subsequent accurate spatiotemporal alignment.

[0056] Specifically, the equivalent refractive index at each location along the path between the spatial grid search point and the sensor is first obtained from a three-dimensional spatial database. The three-dimensional spatial database is established through a pre-conducted active seismic source calibration experiment. Multiple active seismic source emission points at known locations are set up within the monitoring area, and these emission points are controlled to sequentially emit test signals, recording the arrival time of the received test signals by the sensor array. Based on the location of the active seismic source emission points, the sensor location, and the arrival time, a tomographic inversion algorithm is used to calculate the distribution of sound wave propagation velocity in the monitored area at different frequency bands. The sound wave propagation velocity... Mapped to equivalent refractive index ,in For reference velocity, the speed of sound in air is usually taken. The equivalent refractive index is then expressed according to three-dimensional spatial coordinates. The data is stored to form a discrete three-dimensional spatial database.

[0057] In actual queries, for the line connecting the spatial grid search point and the sensor, samples are taken at 0.5-meter intervals along the direction of this line. At each sampling location, the equivalent refractive index at that location is obtained from the database through trilinear interpolation. This step yields spatial distribution parameters describing the non-uniformity of the medium. These parameters directly reflect the degree to which the sound wave propagation path deviates from a straight line, providing necessary input data for subsequent ray tracing calculations and effectively resolving the systematic errors caused by the traditional homogeneous medium assumption.

[0058] After obtaining the equivalent refractive index distribution, the actual propagation path of the sound wave in the non-uniform medium is calculated based on the equivalent refractive index. The propagation of sound waves in a non-uniform medium follows Fermat's principle, meaning the sound wave propagates along the path with the shortest propagation time. A ray tracing algorithm is used to solve for this shortest time path; in practice, an improved iterative algorithm is employed. Starting from a search point on the spatial grid, the initial propagation direction is set towards the sensor position, and the propagation path is discretized into a series of small segments. At the endpoint of each segment, the propagation path is calculated based on the equivalent refractive index gradient at that location. Calculate the change in the refraction angle of the sound wave. The refraction angle satisfies the following relationship: ,in and The equivalent refractive index of adjacent segments. and This refers to the corresponding communication angle.

[0059] Through iterative calculations, the propagation direction is gradually adjusted so that the path converges to the sensor position while satisfying the law of refraction. During the iteration process, a variable step-size strategy is employed: a smaller step-size is used in regions with a large equivalent refractive index gradient to improve accuracy, and a larger step-size is used in regions with a small gradient to improve efficiency. After multiple iterations, the actual propagation path from the spatial grid search point to the sensor is obtained. For example, for a spatial grid search point and the sensor 80 meters apart, the propagation path is a straight line with a length of 80 meters under the assumption of a homogeneous medium. However, considering non-homogeneity, the actual propagation path exhibits a curved shape, increasing the path length to 83.5 meters, and the maximum deviation from a straight line reaches 4.2 meters. This step accurately describes the actual propagation trajectory of sound waves in complex geological environments. This trajectory considers the refraction effect caused by medium non-homogeneity, significantly improving the realism of the propagation path modeling and laying the geometric foundation for subsequent frequency-dependent propagation time calculations. It effectively overcomes the inapplicability of the straight-line propagation assumption in non-homogeneous media.

[0060] After determining the propagation path, the path is divided into multiple micro-segments, and the high-frequency equivalent propagation velocity and low-frequency equivalent propagation velocity corresponding to each micro-segment are obtained. The segment length is usually set to 0.1 meters to 0.2 meters to ensure that the change in medium parameters within each micro-segment is small enough to approximate a homogeneous medium.

[0061] For each micro-segment, the spatial coordinates of its center are first determined, and then the medium parameters at that location are queried from a three-dimensional spatial database. Because geological media such as rocks have significant dispersion characteristics, sound waves of different frequencies propagate at different speeds within the same medium. The high-frequency equivalent propagation velocity at that location is then obtained from the database. and low-frequency equivalent propagation speed The high-frequency equivalent propagation velocity corresponds to the dominant frequency range of the acoustic emission signal, typically between 100kHz and 500kHz, while the low-frequency equivalent propagation velocity corresponds to the dominant frequency range of the micro-vibration signal, typically between 10Hz and 1000Hz.

[0062] In rock media, high-frequency sound waves typically propagate faster than low-frequency sound waves. This is because high-frequency sound waves mainly propagate along the interior of grains, while low-frequency sound waves are affected by structures such as grain boundaries and fissures. For example, in a certain granite region, the equivalent propagation speed of high frequency is 5800 meters per second, while the equivalent propagation speed of low frequency is 5200 meters per second, a difference of approximately 11.5%.

[0063] For each micro-segment along the propagation path, its corresponding high-frequency and low-frequency propagation velocities are obtained to form a velocity sequence. and ,in The total number of micro-segments. Through this step, a frequency-dependent propagation velocity model was established. This model finely characterizes the differentiated impact of medium dispersion characteristics on the propagation of signals at different scales, providing a parameter basis for accurately calculating the propagation time difference between high-frequency and low-frequency signals, and effectively solving the problem of time delay estimation bias caused by neglecting dispersion effects in traditional methods.

[0064] After obtaining the propagation speed of each micro-segment, the high-frequency propagation time and low-frequency propagation time of each micro-segment are calculated. For the th... Each infinitesimal segment has a length of _____. The high-frequency equivalent propagation speed is The low-frequency equivalent propagation speed is Based on the fundamental relationship that time equals distance divided by velocity, the high-frequency propagation time of this infinitesimal segment is calculated. and low-frequency propagation time For example, for a micro-segment with a length of 0.15 meters, the high-frequency equivalent propagation speed is 5800 meters per second, and the low-frequency equivalent propagation speed is 5200 meters per second. Therefore, the high-frequency propagation time is... microseconds, low-frequency propagation time is The high-frequency and low-frequency propagation time difference of this micro-segment is 2.99 microseconds. This calculation process is repeated for all micro-segments along the propagation path to obtain the high-frequency propagation time series for each micro-segment. and low-frequency propagation time series This step decomposes the overall propagation path time calculation into fine calculations at the micro-segment level, fully considering the spatial variations of medium parameters and frequency-related propagation speed differences along the propagation path. This significantly improves the accuracy of propagation time calculation and provides high-quality basic data for subsequent summation.

[0065] Finally, the high-frequency propagation time is obtained by summing the high-frequency propagation times of all infinitesimal segments, and the low-frequency propagation time is obtained by summing the low-frequency propagation times of all infinitesimal segments. The formula for calculating the high-frequency propagation time is as follows: The formula for calculating the low-frequency propagation time is: For example, for a propagation path containing 550 micro-segments, the high-frequency propagation time, when summed up, is 14.23 milliseconds, while the low-frequency propagation time, when summed up, is 15.87 milliseconds, a difference of 1.64 milliseconds. This time difference reflects the cumulative effect of dispersion along the entire propagation path. The high-frequency and low-frequency propagation times will serve as key parameters in the subsequent spatiotemporal alignment step, used for reverse translation compensation of the high-frequency energy packet characteristic flow and the low-frequency band energy characteristic flow. Through this step, the complete propagation time calculation from the spatial grid search point to the sensor is completed. This calculation comprehensively considers the path bending effect caused by medium inhomogeneity and the velocity difference effect caused by dispersion characteristics, significantly improving the accuracy of propagation delay estimation.

[0066] This implementation achieves accurate calculation of the propagation path and propagation time of sound waves in non-homogeneous media, effectively solving the systematic errors caused by the traditional homogeneous medium assumption and single propagation velocity model. By obtaining the equivalent refractive index distribution from a three-dimensional spatial database, the non-homogeneous characteristics of the monitoring area are accurately described, providing realistic medium parameter input for ray tracing. The actual propagation path of sound waves is calculated using a ray tracing algorithm, fully considering the refraction effect caused by medium non-homogeneity, making the propagation path model more consistent with physical reality. By dividing the propagation path into micro-segments and obtaining the frequency-dependent propagation velocity of each micro-segment, the differentiated influence of medium dispersion characteristics on signals at different scales is finely characterized, significantly improving the accuracy of propagation velocity modeling. In addition, through the calculation and summation of propagation time at the micro-segment level, accurate estimation of high-frequency and low-frequency propagation time is achieved, providing reliable time delay compensation parameters for subsequent multi-scale feature spatiotemporal alignment. This technical solution effectively overcomes the complex influence of non-homogeneous media and dispersion effects on time delay calculation, significantly improves the accuracy of spatiotemporal alignment, lays a solid theoretical and computational foundation for the accuracy of early disaster source tracing, and strongly supports the reliable operation of early warning systems for sudden disasters.

[0067] In one embodiment of this invention, extracting burst spike features from the high-frequency feature stream includes the following steps: S310. Calculate the first-order time derivative of the high-frequency characteristic flow; S320. Calculate the second time derivative of the first time derivative; S330. Determine whether the second-order time derivative exceeds the preset mutation threshold. S340. When the second time derivative exceeds the preset mutation threshold, the feature point at the corresponding time is determined to be a sudden spike feature.

[0068] In practical implementation, accurately identifying sudden spikes in high-frequency feature streams facilitates the subsequent establishment of multi-scale feature temporal correlations. Traditional methods typically employ simple amplitude threshold judgments or first-order derivative detection. These methods are easily affected by background noise fluctuations, leading to misidentification of normal fluctuations as sudden events or omission of real sudden events with small amplitudes but drastic changes. This implementation introduces second-order time derivative analysis to effectively capture the acceleration characteristics of energy changes, significantly improving the accuracy and anti-interference capability of sudden spike feature identification, and providing a reliable time reference point for subsequent sliding time window construction and fine time alignment.

[0069] Specifically, we first calculate the first-order time derivative of the high-frequency characteristic flow, which reflects the rate of change of high-frequency energy with time. Let the high-frequency energy wave packet characteristic flow be... ,in If the variable is time, then the formula for calculating the first-order time derivative is: In the discretization implementation, the central difference method is used for numerical calculation, that is... ,in The sampling time interval, For the first Each sampling time.

[0070] The physical meaning of the first-order time derivative is the rate of energy increase or decrease; a positive value indicates an energy increase, and a negative value indicates an energy decrease. The absolute value reflects the degree of drastic change. For example, in a certain monitoring, the energy value of the high-frequency characteristic stream suddenly jumped from 0.35 to 0.82 at time 15.2 seconds, and the corresponding first-order derivative jumped from 0.02 at the background level to 1.85, showing a significant energy surge. However, the magnitude of the first-order derivative alone is insufficient to distinguish between gradual energy growth and abrupt energy jumps, as both can produce large first-order derivative values. This step obtains information on the energy change rate of the high-frequency characteristic stream, which provides basic data for subsequent abrupt change detection, but is not enough to accurately identify sudden spikes; further analysis of the change rate itself is needed.

[0071] After obtaining the first-order time derivative, the second-order time derivative is calculated. This second-order derivative reflects the rate of change of energy, i.e., the acceleration characteristic of energy change. Let the first-order time derivative be... The formula for calculating the second-order time derivative is: In the discretization implementation, the central difference method is also used, i.e. .

[0072] The physical meaning of the second derivative is the degree of acceleration or deceleration of energy change. A positive value indicates that the rate of energy growth is accelerating, while a negative value indicates that the rate of energy growth is slowing down or the rate of energy decay is accelerating. The essential characteristic of a sudden spike is that energy changes abruptly within a very short time. This jump corresponds to a sudden increase in the first derivative, which in turn corresponds to the maximum value of the second derivative. In contrast, although the first derivative may be larger in a gradual energy growth pattern, the change in the first derivative is gradual, resulting in a smaller second derivative.

[0073] For example, in the energy surge event at 15.2 seconds, the first derivative jumps from 0.02 to 1.85, and the corresponding second derivative reaches a peak of 91.5 at 15.2 seconds, far exceeding the background level of 3.2. In contrast, for a gradual energy increase process, the first derivative slowly increases from 0.02 to 1.20, with a maximum second derivative of only 8.6. This step extracts the energy change acceleration feature of the high-frequency characteristic stream. This feature can effectively distinguish between abrupt and gradual changes, providing a key criterion for accurately identifying sudden spikes and significantly improving the sensitivity and specificity of sudden event detection.

[0074] After calculating the second time derivative, it is determined whether the second time derivative exceeds a preset mutation threshold in order to identify the true sudden spike characteristics. The determination of the preset mutation threshold needs to comprehensively consider the background noise level of the monitoring area and the statistical characteristics of historical sudden events.

[0075] In practice, background monitoring data is first collected over a period of time, and the mean of the second derivative is calculated for that period. and standard deviation These statistics reflect the acceleration characteristics of energy fluctuations under normal conditions. The preset mutation threshold is set to... ,in This is a multiplier, typically ranging from 5 to 10. Smaller... A higher value increases detection sensitivity but may increase false alarms; a larger value... The value reduces false alarms but may miss weak emergencies.

[0076] For example, in the monitoring of an underground mine, the mean of the background second derivative was 3.2, and the standard deviation was 1.8. [The following is a separate, unrelated sentence:] Setting... The preset mutation threshold is For each sampling time Determine the second derivative at that moment. Does it exceed the preset mutation threshold? .like If so, it is determined that there is a sudden spike characteristic at that moment; if If the signal is clear, then that moment is considered a normal fluctuation. This step fully considers the randomness of background noise and the salience of sudden events, effectively reducing the false positive rate and ensuring that the identified sudden spike features have real physical meaning.

[0077] After determining that the second-order time derivative exceeds a preset abrupt change threshold, the feature point at the corresponding time is identified as a sudden spike feature, and the timestamp and related parameters of this feature point are recorded. For each identified sudden spike feature, its occurrence time is recorded. The corresponding energy value First derivative value and second derivative value These parameters fully describe the temporal location and intensity characteristics of the sudden spike features.

[0078] In actual monitoring, multiple moments exceeding the threshold may be detected within a short period of time. These moments may correspond to different stages of the same sudden event. To avoid duplicate identification, a temporal clustering method is used to group multiple abrupt change moments with a time interval less than a preset clustering window (e.g., 0.1 seconds) into a single sudden peak feature, and the moment with the largest second derivative is selected as the representative timestamp of that sudden peak feature.

[0079] For example, if five times exceeding the threshold are detected between 15.18 and 15.24 seconds, with the largest second derivative value (91.5) at 15.21 seconds, then 15.21 seconds is determined as the timestamp of this sudden spike feature. The identified sudden spike feature will serve as the center time for constructing the subsequent sliding time window, used to establish the temporal correlation between high-frequency and low-frequency feature flows. This step achieves precise localization and parameter extraction of the sudden spike feature, providing an accurate time reference point for subsequent multi-scale feature fine alignment.

[0080] This implementation achieves accurate identification of sudden spike features in high-frequency feature streams, effectively solving the problems of traditional simple thresholding methods being susceptible to noise interference and difficulty in distinguishing between abrupt changes and gradual changes. Energy change rate information is obtained by calculating the first-order time derivative, providing basic data for subsequent acceleration feature extraction. Acceleration features of energy change are extracted by calculating the second-order time derivative; this feature can effectively distinguish between sudden jumps and gradual increases, significantly improving the discriminative ability of sudden event detection. An adaptive sudden event judgment criterion is established by setting a preset abrupt change threshold based on background noise statistical features, effectively reducing the false positive and false negative rates. Furthermore, temporal clustering avoids duplicate identification, ensuring that each sudden spike feature corresponds to a real physical event, improving the accuracy of feature identification. This technical solution provides a reliable time reference point for subsequent sliding time window construction and multi-scale feature fine alignment, strongly supporting the accuracy and reliability of early disaster source tracing.

[0081] In one embodiment of this invention, the normalized cross-correlation function of the high-frequency feature stream and the low-frequency feature stream is calculated, and the residual time delay is determined based on the peak position of the normalized cross-correlation function, including the following steps: S410. Obtain the spatiotemporally aligned high-frequency feature stream segment and the spatiotemporally aligned low-frequency feature stream segment within the sliding time window; S420. Normalize the high-frequency characteristic flow segments and the low-frequency characteristic flow segments respectively; S430. Calculate the cross-correlation function between the normalized high-frequency characteristic flow segment and the normalized low-frequency characteristic flow segment, and determine the peak position of the cross-correlation function. S440. Determine the residual time delay of the spatiotemporally aligned low-frequency characteristic stream segment relative to the spatiotemporally aligned high-frequency characteristic stream segment based on the peak position.

[0082] Before acquiring the spatiotemporally aligned high-frequency and low-frequency feature stream segments within the sliding time window, the process may further include: estimating the approximate propagation distance from the burst source corresponding to the high-frequency feature stream to the sensor array; calculating the distribution range of the physical time delay difference between the high-frequency and low-frequency signals due to dispersion effects at the current propagation distance based on the acoustic dispersion physical model; determining the upper limit of the physical time delay difference distribution range as the reference physical width of the sliding time window; extracting the energy accumulation rate of change near the burst peak features in the high-frequency feature stream in real time, and determining the time-series inflection point where the energy accumulation rate of change changes from increasing to flat; adaptively adjusting the forward and backward time-series boundaries of the sliding time window according to the time delay deviation corresponding to the time-series inflection point, so that the sliding time window encompasses the complete low-frequency energy packet; introducing window functions with smooth transition characteristics at both ends of the adjusted sliding time window to smoothly attenuate the feature amplitude at the window edge, thereby eliminating the spectral leakage effect caused by time-domain truncation; and outputting the adjusted and smoothed sliding time window for truncation of the spatiotemporally aligned high-frequency and low-frequency feature stream segments.

[0083] In practice, although preliminary spatiotemporal alignment has been achieved through reverse translation compensation based on a 3D spatial database, slight residual time delays still exist after the initial alignment due to errors in medium parameter estimation, simplification of the dispersion model, and discretization errors in propagation path calculation. Traditional methods typically ignore these residual time delays or use fixed compensation values, resulting in insufficient alignment accuracy of multi-scale features on the time axis, which in turn affects the quality of feature fusion and the accuracy of source tracing. This implementation introduces normalized cross-correlation function analysis to accurately measure the time offset between high-frequency and low-frequency feature flows within a local time window, effectively eliminating residual time delay errors and significantly improving the time synchronization accuracy of multi-scale features, providing high-quality aligned data for subsequent feature fusion.

[0084] Specifically, before calculating the cross-correlation function, a suitable sliding time window needs to be constructed to extract the feature stream segments. First, the approximate propagation distance from the burst source corresponding to the high-frequency feature stream to the sensor array is estimated. This distance can be roughly estimated by analyzing the time differences between the burst signals received by multiple sensors and using a simplified time-difference localization algorithm. For example, if the time differences between the burst signals received by the three sensors are 3.2 milliseconds and 5.8 milliseconds, respectively, and the sensor spacing is 20 meters, then the distance from the burst source to the sensor array can be estimated to be approximately 45 to 60 meters.

[0085] Based on a physical model of acoustic wave dispersion, the distribution range of the physical time delay difference between high-frequency and low-frequency signals due to dispersion effects at the current propagation distance is calculated. The dispersion physical model considers the difference in propagation speed of sound waves of different frequencies in rock media. For a typical granite medium, at a propagation distance of 50 meters, the time delay difference between the high-frequency signal (200kHz) and the low-frequency signal (100Hz) is approximately 0.8 milliseconds to 1.5 milliseconds. The upper limit of the physical time delay difference distribution range is determined as the baseline physical width of the sliding time window, for example, set to 2.0 milliseconds, to ensure that the window can cover the maximum possible time delay difference.

[0086] However, a fixed-width window may not be suitable for the energy release characteristics of different burst events, thus requiring adaptive adjustment. The rate of change of energy accumulation near burst spikes in the high-frequency feature stream is extracted in real time; this rate of change reflects the duration of energy release. The energy accumulation curve is then calculated. ,in The timestamp of the sudden spike characteristic is used, and then the cumulative rate of change of energy is calculated. Determine the time-series inflection point where the rate of change of energy accumulation transitions from increasing to leveling off; this inflection point corresponds to the end of the energy release from the sudden event.

[0087] For example, if the rate of change of energy accumulation is at time... The moment when the latency drops from the peak of 0.85 to 0.15 milliseconds is considered the timing inflection point. Based on the latency deviation corresponding to the inflection point, the forward and backward timing boundaries of the sliding time window are adaptively adjusted. The forward boundary is set to... The backward boundary is set to the time inflection point plus 0.5 milliseconds to ensure that the sliding time window encompasses the complete low-frequency energy packet.

[0088] Window functions with smooth transition characteristics, such as Hanning or Gaussian windows, are introduced at both ends of the adjusted sliding time window to smoothly attenuate the characteristic amplitudes at the window edges. The attenuation coefficient gradually decreases from 1.0 at the center of the window to 0.0 at the edges, thus eliminating the spectral leakage effect caused by time-domain truncation. The adjusted and smoothed sliding time window is output, with a time range of [missing information]. The total width is 2.0 milliseconds. This step constructs an adaptive sliding time window that fully considers the physical constraints of dispersion effects and the energy release characteristics of sudden events. It effectively solves the problem that fixed windows cannot adapt to different event characteristics, providing a reasonable time range for subsequent feature segment extraction and cross-correlation calculations.

[0089] After constructing the sliding time window, spatiotemporally aligned high-frequency feature stream segments and spatiotemporally aligned low-frequency feature stream segments are obtained within the sliding time window. For the high-frequency feature stream, the extraction time range is [missing information]. The energy value sequence within a millisecond is denoted as ,in This represents the number of sampling points within the time window. For example, if the sampling rate of a high-frequency signal is 100kHz, then a 2.0-millisecond window contains 200 sampling points.

[0090] For low-frequency characteristic flows, the energy value sequence within the same time range is also extracted, denoted as... ,in This represents the number of sampling points within the time window. Since low-frequency signals typically have a low sampling rate, such as 10kHz, a 2.0ms window contains 20 sampling points. To facilitate subsequent cross-correlation calculations, the low-frequency feature stream segment needs to be interpolated and upsampled to match the number of sampling points of the high-frequency feature stream segment. Using cubic spline interpolation, the low-frequency feature stream segment is interpolated from 20 sampling points to 200 sampling points, resulting in... This step extracts high-frequency and low-frequency feature stream segments within the sliding time window. These segments contain the complete energy release process of the sudden event, providing raw data for subsequent normalization and cross-correlation calculations.

[0091] After acquiring the feature flow segments, high-frequency and low-frequency feature flow segments are normalized separately to eliminate the influence of amplitude differences on cross-correlation calculations. The purpose of normalization is to make feature flows of different scales comparable and to prevent feature flows with larger amplitudes from dominating the cross-correlation results.

[0092] For high-frequency characteristic flow segments First, calculate its mean. and standard deviation Then, a normalization transformation is performed on each sampling point to obtain the normalized high-frequency feature stream segment. The normalized characteristic current has a mean of 0 and a standard deviation of 1, eliminating the influence of the original amplitude.

[0093] For low-frequency characteristic flow segments The mean was calculated using the same normalization method. and standard deviation The normalized low-frequency characteristic flow segment is obtained. For example, in a certain calculation, the original amplitude range of the high-frequency feature flow segment is 0.15 to 0.92, and the normalized amplitude range is -1.57 to 2.10. The original amplitude range of the low-frequency feature flow segment is 0.08 to 0.35, and the normalized amplitude range is -1.38 to 2.00. After normalization, the two feature flow segments have the same statistical properties, allowing the cross-correlation function to accurately reflect the waveform similarity between them, unaffected by the original amplitude differences. This step achieves amplitude standardization of multi-scale feature flows, providing preprocessing assurance for accurate calculation of the cross-correlation function.

[0094] After normalization, the cross-correlation function between the normalized high-frequency characteristic flow segment and the normalized low-frequency characteristic flow segment is calculated, and the peak position of the cross-correlation function is determined. The cross-correlation function is defined as follows: ,in The time offset represents the time shift of a low-frequency feature stream segment relative to a high-frequency feature stream segment.

[0095] By iterating through different time offsets Calculate the cross-correlation function value at each offset. The search range for the time offset is typically set to... Milliseconds, the search step size is set to the sampling interval, for example, 0.01 milliseconds. For each time offset Calculate the cross-correlation function value This value reflects the time offset. The waveform similarity between high-frequency characteristic flow segments and low-frequency characteristic flow segments is then analyzed.

[0096] Among all time offsets, determine the maximum value of the cross-correlation function and its corresponding time offset, i.e., the peak position. For example, in a certain calculation, the cross-correlation function at a time offset... The peak value of 0.89 is reached at milliseconds, indicating that the low-frequency feature flow segment has a lead of 0.12 milliseconds over the high-frequency feature flow segment.

[0097] To improve the estimation accuracy of the peak position, a parabolic interpolation method is used to fit the cross-correlation function values ​​near the peak to obtain a sub-sampling accuracy peak position. Let three points near the peak position be... , and The precise peak position is obtained through parabolic fitting. ,in The sampling interval is defined as follows. This step accurately calculates the time offset between the high-frequency and low-frequency feature flows. This offset reflects the residual time delay after the initial spatiotemporal alignment, providing precise compensation parameters for subsequent fine-tuning and translation.

[0098] After determining the peak position of the cross-correlation function, the residual time delay of the spatiotemporally aligned low-frequency feature stream segment relative to the spatiotemporally aligned high-frequency feature stream segment is determined based on the peak position. Peak position This is the estimated residual time delay; the sign and magnitude of this value reflect the direction and magnitude of the time adjustment required for the low-frequency characteristic flow. If This indicates that the low-frequency feature flow segment lags behind the high-frequency feature flow segment, and the low-frequency feature flow needs to be shifted forward on the time axis. ;like This indicates that the low-frequency feature flow segment leads the high-frequency feature flow segment, and the low-frequency feature flow needs to be shifted backward on the time axis. .

[0099] Based on residual delay Fine-tuning the translation of the entire low-frequency feature stream, that is, for each moment in the low-frequency feature stream. Shift its corresponding energy value to time 1. In the discretization implementation, time-domain interpolation is used to shift non-integer sampling points. For example, if the residual delay is -0.12 milliseconds and the sampling interval is 0.01 milliseconds, the low-frequency feature stream needs to be shifted backward by 12 sampling points. For shifting non-integer multiple sampling points, Sinc interpolation or fractional delay filters are used to achieve high-precision time-domain shifting.

[0100] After fine-tuning and shifting, the low-frequency and high-frequency feature streams are precisely aligned on the time axis, achieving maximum similarity at critical moments of sudden events. For example, in one processing iteration, after initial spatiotemporal alignment, the energy peak of the high-frequency feature stream occurs at 15.21 seconds, while the energy peak of the low-frequency feature stream occurs at 15.33 seconds, a difference of 0.12 milliseconds. After residual delay compensation, the energy peak of the low-frequency feature stream is adjusted to 15.21 seconds, achieving precise synchronization with the high-frequency feature stream. The fine-tuned and shifted high-frequency and low-frequency feature streams together constitute a multi-scale feature, which is highly aligned in the time dimension, providing high-quality input data for subsequent spatial fusion. This step achieves accurate estimation and compensation of residual delay, effectively eliminating residual errors from the initial spatiotemporal alignment and realizing fine-grained time synchronization of multi-scale features.

[0101] This implementation achieves accurate estimation and compensation of residual time delay between high-frequency and low-frequency feature flows, effectively solving the residual error problem in initial spatiotemporal alignment. By constructing an adaptive window based on a dispersion physics model and energy release characteristics, the integrity and rationality of feature flow segment extraction are ensured, providing a high-quality data range for cross-correlation calculation. Normalization eliminates amplitude differences between feature flows of different scales, enabling the cross-correlation function to accurately reflect waveform similarity and significantly improving the accuracy of time delay estimation. By calculating the normalized cross-correlation function and using parabolic interpolation to determine the peak position, sub-sampling precision residual time delay estimation is achieved, effectively improving the accuracy of time alignment. Furthermore, by fine-tuning the translation of the low-frequency feature flow based on the peak position, fine-grained time synchronization of multi-scale features is achieved, providing high-quality aligned data for subsequent feature fusion. This technical solution effectively overcomes the residual time delay problems caused by medium parameter estimation errors and dispersion model simplification, significantly improving the accuracy of multi-scale feature spatiotemporal alignment and providing a solid data foundation for accurate early disaster source tracing.

[0102] In one embodiment of this invention, multi-scale features are fused to obtain a fused feature map, including the following steps: S510: Obtain the spatial coordinates of each sensor in the acoustic emission sensor array and the micro-vibration sensor array; S520. Construct a graph structure based on spatial coordinates, where multi-scale features serve as the initial vectors for the graph nodes in the graph structure. S530. Calculate the attention coefficients between adjacent graph nodes based on the graph attention network, and calculate the spatial gradient of the attention coefficient of each graph node. Adjust the spatial receptive field radius of the graph node according to the gradient. S540. Based on the adjusted spatial receptive field radius, multi-scale features are fused to output a fused feature map.

[0103] Traditional methods typically fuse multi-sensor data using simple weighted averaging or direct superposition. These methods ignore the spatial correlation between sensors and the spatial distribution characteristics of feature energy, resulting in fusion results that fail to accurately reflect the spatial location characteristics of the disaster source. Especially in the early stages of a disaster, the signal energy density is high at the source and low in areas far from the source. Simple fusion methods cannot adaptively adjust the fusion strategy for different spatial regions, easily generating false peaks or obscuring the true source signal. This implementation effectively solves the spatial fusion problem of multi-scale features by introducing a graph attention network and a spatial receptive field adaptive adjustment mechanism, significantly improving the ability of the fused feature map to represent the location of the disaster source.

[0104] Specifically, the spatial coordinates of each sensor in the acoustic emission sensor array and the microseismic sensor array are first obtained. During the sensor deployment phase, a high-precision total station or laser rangefinder is used to measure the three-dimensional spatial position of each sensor, achieving centimeter-level accuracy. The acoustic emission sensor array contains... There are 1 sensor node, and the spatial coordinates of each node are as follows: ,in The microseismic sensor array contains There are 1 sensor node, and the spatial coordinates of each node are denoted as . ,in .

[0105] For example, in a monitoring system for an underground mine, the acoustic emission sensor array contains 32 nodes, and the microseismic sensor array contains 16 nodes. The spatial coordinates of all sensors were accurately recorded and stored in the system database during deployment. The spatial coordinates of the sensors not only describe their physical locations but also implicitly contain their spatial proximity and relative geometric configurations. This spatial information is crucial for subsequent graph structure construction and spatial feature fusion because the propagation of disaster source signals in space follows physical laws; signals received by adjacent sensors have strong correlations, while signals received by sensors farther away have weaker correlations. Through this step, the basic spatial topology data of the sensor network is obtained, providing geometric parameters for subsequent graph structure construction.

[0106] After acquiring the spatial coordinates of the sensors, a graph structure is constructed based on these coordinates, using multi-scale features as the initial vectors for the graph nodes. A graph structure is a data structure that effectively represents spatial relationships, consisting of nodes and edges. Each sensor is mapped to a node in the graph structure, with a total number of nodes. For each node Its initial eigenvector It consists of the multi-scale features corresponding to this sensor.

[0107] Specifically, for acoustic emission sensor nodes, their feature vectors contain high-frequency energy packet features after spatiotemporal alignment and residual delay compensation; for microseismic sensor nodes, their feature vectors contain low-frequency band energy features after spatiotemporal alignment and residual delay compensation. The dimension of the feature vectors is typically set to 64 to 128 dimensions, encompassing the temporal statistical characteristics, frequency domain distribution characteristics, and energy evolution characteristics of the features.

[0108] In the graph structure, edges represent the spatial proximity between nodes, and the K-nearest neighbor method is used to construct edge connections. For each node, the Euclidean distance between it and all other nodes is calculated. Choose the one closest to you. Each node is designated as a neighbor node, and an edge connection is established between the node and its neighbor nodes. The value is usually set to 5 to 10 to ensure that the graph structure captures local spatial relationships without being too dense.

[0109] For example, for a certain coordinate ( The acoustic emission sensor nodes are located at distances of 8.5 meters, 12.3 meters, 15.7 meters, 18.2 meters, and 21.6 meters, respectively, and edge connections are established with these five nodes. This step transforms the dispersed sensor network into a graph structure with well-defined spatial topological relationships. This structure preserves the multi-scale feature information of each sensor while establishing spatial relationships between sensors, providing structured input data for subsequent graph attention network processing.

[0110] After constructing the graph structure, the attention coefficients between adjacent graph nodes are calculated using the graph attention network, and the gradient of the attention coefficient of each graph node in space is calculated. The spatial receptive field radius of the graph nodes is then adjusted based on the gradient. The graph attention network is a neural network architecture that can adaptively learn the importance weights between nodes. Its core idea is to dynamically calculate the contribution weights of adjacent nodes to the current node through the attention mechanism.

[0111] For nodes and its neighboring nodes Calculate the attention coefficient The formula is ,in , The weight matrix is ​​a learnable matrix. For attention weight vectors, This represents a vector concatenation operation. For nodes The set of neighboring nodes. Attention coefficient. Reflects neighboring nodes Features of nodes The greater the coefficient, the stronger the correlation between the features of the neighboring node and the features of the current node.

[0112] Calculate the gradient of the attention coefficient in space for each graph node; this gradient reflects the rate of change of the attention coefficient in three-dimensional space. For each node... The spatial gradient of its attention coefficient is calculated as follows: The gradient magnitude is approximated using the finite difference method. This reflects the degree of spatial variation in the attention coefficient near the node.

[0113] A large gradient magnitude indicates that the node is located in a region of rapidly changing feature energy density, potentially near the source of a disaster; a small gradient magnitude indicates that the node is located in a region of relatively flat feature density, far from the source of a disaster. The spatial receptive field radius of the graph node is adjusted based on the gradient magnitude. Specifically, this is done by determining whether the gradient magnitude exceeds a preset gradient threshold. .like Determine the node In regions with increased characteristic energy density, the spatial receptive field radius is reduced to a preset minimum radius. For example, 3 meters, to improve the spatial resolution of the area; if Keep the node The spatial receptive field radius is a preset standard radius. For example, 10 meters, to cover a larger area.

[0114] For example, in one processing step, the gradient magnitudes of the attention coefficients of the five nodes closest to the disaster source were 0.082, 0.095, 0.118, 0.091, and 0.076, respectively, all exceeding the preset gradient threshold of 0.070, thus shrinking their spatial receptive field radius to 3 meters. Meanwhile, the gradient magnitudes of the other nodes farther from the source were all less than 0.070, maintaining a standard radius of 10 meters. This step achieves adaptive spatial receptive field adjustment based on the spatial gradient of the attention coefficients. This mechanism can dynamically adjust the fusion range according to the spatial distribution characteristics of the feature energy, employing a fine-grained fusion strategy in the source region and a coarse-grained fusion strategy in the far-field region.

[0115] After adjusting the spatial receptive field radius, multi-scale features are fused based on the adjusted spatial receptive field radius to output a fused feature map. For each graph node... With its spatial coordinates Centered on the adjusted spatial receptive field radius Using the boundary as the boundary, the spatial action sphere of this node is divided in three-dimensional space. .

[0116] Extract Nodes The corresponding high-frequency acoustic emission characteristics and low-frequency microseismic characteristics are denoted as follows: and High-frequency features Mapped to a virtual sound pressure scalar field within a spatial action sphere low-frequency features Mapped to virtual work vector field The mapping method uses radial basis function interpolation, i.e. ,in The distance from a spatial point to the center of a node. Parameters used to control the decay rate.

[0117] Within the overlapping region of the spatial action sphere of adjacent graph nodes, gradient divergence cross-coupling operations are performed on the virtual sound pressure scalar field and the virtual work vector field to generate a local spatiotemporal coupling tensor. For spatial points within the overlapping region... Calculate the scalar gradient of the sound pressure field at all relevant nodes at that point. and virtual work vector field divergence Then perform cross-coupling operations. ,in The set of nodes that cover this spatial point. This represents the tensor outer product.

[0118] Calculate the characteristic norm of the local spatiotemporal coupling tensor in three-dimensional space. This norm quantitatively characterizes the physical interaction strength of features at different scales within the overlapping space. By interpolating the feature norm along the spatial geometric path of the overlapping region using radial basis functions, a continuous feature energy distribution field is generated. .

[0119] A bilateral filter is applied to the characteristic energy distribution field, taking into account both spatial distance and energy value differences. The filtering formula is as follows: ,in For spatial distance, Due to differences in energy values, The normalization coefficient is... This serves as the filtering window. Bilateral filtering smooths noise while preserving energy abrupt changes, outputting a fused feature map. .

[0120] A fused feature map is a scalar field defined on a three-dimensional spatial grid, where the value at each grid point reflects the likelihood of that location being a hazard source. For example, in a certain processing, the fused feature map is located at coordinates ( The energy value at point (0.92) reached a peak, while the energy values ​​in the surrounding area were all below 0.35, clearly indicating the location of the disaster source. Through this step, deep spatial fusion of multi-scale features was achieved. The fusion process fully considered the spatial topological relationship of the sensors, the spatial distribution characteristics of feature energy, and the physical interaction mechanism of features at different scales, generating a high-quality fused feature map, which provides accurate spatial energy distribution information for subsequent source location determination.

[0121] This implementation achieves effective spatial fusion of multi-scale features, effectively solving the problems of traditional simple fusion methods being unable to adapt to spatial heterogeneity and ignoring sensor spatial correlations. By acquiring sensor spatial coordinates and constructing a graph structure, a spatial topological representation of the sensor network is established, providing a structured framework for subsequent spatial correlation modeling. Attention coefficients between adjacent nodes are calculated using a graph attention network, enabling adaptive evaluation of sensor feature importance and effectively capturing spatial correlations. Spatial adaptation of the fusion strategy is achieved by calculating the spatial gradient of the attention coefficients and adjusting the spatial receptive field radius based on the gradient. Fine-grained fusion is used to improve resolution in the source region, while coarse-grained fusion is used in the far-field region to improve coverage. Furthermore, deep physical fusion of multi-scale features is achieved through virtual field mapping based on the adjusted receptive field radius, gradient divergence cross-coupling, and bilateral filtering, generating a high-quality fused feature map. This technical solution significantly improves the fused feature map's ability to represent the location of disaster sources and its spatial resolution, providing key technical support for the accuracy of early disaster source tracing.

[0122] In one embodiment of this invention, adjusting the spatial receptive field radius of a graph node based on a gradient includes the following steps: S610, Calculate the magnitude of the gradient of the attention coefficient; S620. Determine whether the magnitude of the gradient of the attention coefficient exceeds the preset gradient threshold. S630. If the magnitude of the gradient of the attention coefficient exceeds the preset gradient threshold, determine that the graph node is in the region of increased feature energy density, and shrink the spatial receptive field radius of the graph node to the preset minimum radius. S640. If the magnitude of the gradient of the attention coefficient does not exceed the preset gradient threshold, keep the spatial receptive field radius of the graph node at the preset standard radius.

[0123] In the early stages of a disaster, the signal energy in the source region exhibits a sharp spatial attenuation characteristic, requiring a smaller fusion range for accurate source location. Conversely, the signal energy in areas far from the source changes gradually, necessitating a larger fusion range to fully utilize spatial information. This implementation effectively addresses the issues of insufficient spatial resolution and wasted computational resources caused by a fixed fusion range by introducing an adaptive adjustment mechanism for the spatial receptive field based on attention coefficient gradients. This significantly improves the sensitivity and positioning accuracy of the fused feature map to the disaster source location.

[0124] Specifically, the magnitude of the gradient of the attention coefficient is first calculated, which quantitatively characterizes the degree of drastic change of the attention coefficient in three-dimensional space. For graph nodes... Its attention coefficient It is calculated using a graph attention network and reflects the strength of feature correlation between the node and its neighboring nodes. The spatial gradient of the attention coefficient is defined as follows: This indicates that the attention coefficient is in , , Rates of change in the three spatial directions.

[0125] In the discretization implementation, the central difference method is used to calculate the partial derivatives in each direction. For Direction, calculated using the following formula: ,in The spatial step size is typically set to 1 to 2 meters. Indicates at node of Coordinates increase Then, the attention coefficient values ​​are obtained through interpolation or recalculation. For and The direction is calculated using the same method.

[0126] The formula for calculating the magnitude of the gradient is: This modulus is a scalar that reflects the intensity of spatial variation of the attention coefficient near the node, independent of specific directions. For example, in a certain calculation, the node... Located at coordinates ( Its attention coefficient is 0.78, in The partial derivative in the direction is 0.065. The partial derivative in the direction is 0.082. If the partial derivative in the direction is 0.041, then the gradient magnitude is... .

[0127] A large gradient magnitude indicates drastic changes in the attention coefficient near the node, implying a significant spatial gradient in the characteristic energy density of the region, potentially indicating proximity to the disaster source. Conversely, a small gradient magnitude indicates a gradual change in the attention coefficient near the node, suggesting a more uniform distribution of characteristic energy density in the region, away from the disaster source. This step yields a quantitative indicator of the spatial variation intensity of the attention coefficient for each graph node, providing a basis for subsequent adjustments to the spatial receptive field radius.

[0128] After calculating the gradient magnitude, it is determined whether the magnitude of the gradient of the attention coefficient exceeds a preset gradient threshold to identify graph nodes located in regions where feature energy density changes rapidly. Preset gradient threshold. The setting needs to take into account the spatial scale of the monitoring area, the density of sensor deployment, and the statistical characteristics of historical disaster events.

[0129] In practice, background monitoring data is first collected over a period of time. Under normal conditions without disasters, the gradient magnitude of the attention coefficients of all graph nodes is calculated, their distribution characteristics are statistically analyzed, and the mean of the background gradient magnitude is obtained. and standard deviation The preset gradient threshold is set to... ,in This is a multiplier, typically ranging from 3 to 5. Smaller... A higher value results in more nodes being identified as high-gradient regions, improving spatial resolution but increasing computational burden; a larger value... This reduces the number of nodes identified as high-gradient regions, decreasing computational burden but potentially reducing resolution.

[0130] For example, in the monitoring system of an underground mine, the mean background gradient magnitude is 0.025, and the standard deviation is 0.015. The system is set to... The preset gradient threshold is For each graph node Determine its gradient magnitude Does it exceed the preset gradient threshold? .like If the node is located in a region of rapidly changing feature energy density, the spatial receptive field radius needs to be reduced to improve fusion accuracy; if If the node is located in a region with a flat characteristic energy density, the standard spatial receptive field radius is maintained to cover a sufficient spatial range.

[0131] For example, in one monitoring session, 7 out of 48 graph nodes had gradient magnitudes exceeding the threshold of 0.085, specifically 0.109, 0.095, 0.118, 0.091, 0.102, 0.088, and 0.096. These 7 nodes were identified as high-gradient region nodes. The remaining 41 nodes had gradient magnitudes less than 0.085 and were identified as low-gradient region nodes. This step established a spatial region classification criterion based on statistical principles. This criterion can automatically identify key regions with rapidly changing feature energy density, providing a spatial division basis for subsequent differentiated fusion strategies.

[0132] After determining the region type of a node, for nodes whose gradient magnitude exceeds a preset gradient threshold, they are identified as being in a region of increased feature energy density, and their spatial receptive field radius is reduced to a preset minimum radius. The spatial receptive field radius defines the spatial range considered by the node during feature fusion. A smaller radius indicates a more local fusion range and higher spatial resolution; a larger radius indicates a wider fusion range and stronger spatial coverage.

[0133] For nodes located in regions of increased feature energy density, the spatial variation in feature energy in their vicinity is drastic. Using a large fusion range would cause features at different energy levels to be mixed and averaged, obscuring the peak features at the source location and reducing positioning accuracy. Therefore, it is necessary to shrink the spatial receptive field radius, considering only feature information within the node's nearest neighborhood to precisely characterize the spatial location of the energy peak. A preset minimum radius is required. The setting needs to take into account the sensor deployment density and the spatial resolution requirements of the monitoring area, and is usually set to 0.3 to 0.5 times the distance between adjacent sensors.

[0134] For example, if the average sensor spacing is 10 meters, the preset minimum radius is set to 3 to 5 meters. In one implementation, the preset minimum radius was set to 3 meters, and for seven nodes in high-gradient regions, their spatial receptive field radius was reduced from the standard value of 10 meters to 3 meters. After reduction, when these nodes are used for feature fusion, only spatial information within a 3-meter radius of their center is considered, effectively eliminating interference from low-energy features in the far field and highlighting high-energy features near the source.

[0135] For example, nodes Located at coordinates ( Before contraction, the spatial receptive field of the high-gradient region nodes covered a spherical area with a radius of 10 meters and a volume of approximately 4189 cubic meters, containing a large amount of low-energy information far from the source. After contraction, the spatial receptive field of the high-gradient region nodes covered only a spherical area with a radius of 3 meters and a volume of approximately 113 cubic meters, mainly containing high-energy information near the source. The fusion result is more focused on the true source location. Through this step, the spatial receptive field of the high-gradient region nodes was contracted. This mechanism effectively improved the spatial resolution of the source region and enhanced the sensitivity of the fused feature map to the source location.

[0136] For nodes whose gradient magnitude does not exceed a preset gradient threshold, their spatial receptive field radius is maintained at a preset standard radius. These nodes are located in regions with flat feature energy density, where the feature energy in their vicinity changes slowly in space, and there are no abrupt energy gradients. In these regions, using a larger fusion range can fully utilize information from surrounding sensors, improving the stability and noise resistance of feature estimation. (Preset standard radius) The settings need to take into account the spatial scale of the monitoring area and the coverage of the sensor network, and are usually set to 1 to 1.5 times the average spacing between sensors.

[0137] For example, if the average sensor spacing is 10 meters, the preset standard radius is set to 10 to 15 meters. In one implementation, the preset standard radius was set to 10 meters, and for 41 nodes in low-gradient regions, their spatial receptive field radius was maintained at 10 meters. By maintaining the standard radius, when these nodes perform feature fusion, they consider spatial information within a 10-meter radius of their center, which can integrate features from multiple neighboring sensors, smooth local noise fluctuations, and improve the robustness of the fusion results.

[0138] For example, nodes Located at coordinates ( The node, located approximately 70 meters from the disaster source, has a gradient magnitude of 0.032, which does not exceed the threshold of 0.085, maintaining a spatial receptive field radius of 10 meters. This node's spatial action sphere covers a spherical region with a radius of 10 meters, encompassing feature information from 5 to 8 surrounding sensors. The fusion result accurately reflects the average energy level of this region, providing stable background information for the overall feature map. This step preserves the spatial receptive field of nodes in low-gradient regions. This mechanism ensures sufficient spatial coverage in the far-field region, avoids information loss due to excessive contraction, and guarantees the integrity and continuity of the fused feature map.

[0139] This implementation achieves intelligent adjustment of the graph node fusion range, effectively solving the problem that a fixed fusion range cannot adapt to spatial heterogeneity. By calculating the gradient magnitude of the attention coefficient, the intensity of spatial variation of feature energy near each graph node is quantitatively evaluated, providing a quantitative basis for adjusting the fusion strategy. By judging whether the gradient magnitude exceeds a preset gradient threshold, automatic identification and classification of high-gradient and low-gradient regions are achieved, establishing a spatial region division criterion based on statistical principles. By shrinking the spatial receptive field radius of nodes in high-gradient regions to a preset minimum radius, the spatial resolution of the source region is significantly improved, enhancing the sensitivity of the fused feature map to the source location and the positioning accuracy. In addition, by maintaining the spatial receptive field radius of nodes in low-gradient regions at a preset standard radius, sufficient spatial coverage and feature estimation stability in the far-field region are ensured. This technical solution achieves spatial adaptation of the fusion strategy, using fine fusion in the source region to improve positioning accuracy and coarse-grained fusion in the far-field region to improve robustness, significantly improving the quality of multi-scale feature fusion and the accuracy of early disaster source tracing.

[0140] In one embodiment of this invention, multi-scale features are fused based on the adjusted spatial receptive field radius to output a fused feature map, including the following steps: S710. Using the spatial coordinates of each graph node as the center and the corresponding adjusted spatial receptive field radius as the boundary, divide the spatial action sphere of each graph node in three-dimensional space. S720. Extract the high-frequency acoustic emission features and low-frequency micro-vibration features corresponding to each graph node, and map them as virtual acoustic pressure scalar field and virtual work vector field in the spatial action sphere, respectively. S730. Within the overlapping region of the spatial action sphere of adjacent graph nodes, perform gradient divergence cross-coupling calculations of the virtual sound pressure scalar field and the virtual work vector field to generate a local spatiotemporal coupling tensor. S740. Calculate the characteristic norm of the local spatiotemporal coupling tensor in three-dimensional space to quantitatively characterize the physical interaction strength of features at different scales in the overlapping space. S750. The characteristic norm is interpolated radially along the spatial geometric path of the overlapping region to generate a continuous characteristic energy distribution field. S760 performs bilateral filtering on the characteristic energy distribution field and outputs a fused feature map.

[0141] In early disaster monitoring, high-frequency acoustic emission signals reflect the instantaneous propagation of microcracks, while low-frequency microseismic signals reflect the adjustment of macroscopic stress fields. These two phenomena are inherently causally related and energy-coupled. Simple fusion methods cannot capture this deep-seated physical interaction, resulting in fusion results lacking physical meaning and failing to accurately reflect the true location of the disaster source. This implementation method establishes a physical fusion framework for multi-scale features by introducing methods such as virtual physical field mapping, gradient divergence cross-coupling, and radial basis function interpolation. This effectively solves the problem of missing physical mechanisms in traditional methods and significantly improves the physical realism of the fused feature map and the accuracy of source location.

[0142] Reference Figure 3 Specifically, firstly, using the spatial coordinates of each graph node as the center and the corresponding adjusted spatial receptive field radius as the boundary, the spatial domain of action for each graph node is divided in three-dimensional space. For graph nodes... Its spatial coordinates are The adjusted spatial receptive field radius is The radius is adaptively determined based on the gradient magnitude of the node's attention coefficient; the radius of nodes in high-gradient regions is the preset minimum radius. The radius of nodes in the low gradient region is a preset standard radius. .

[0143] The spatial action sphere is defined as , indicating that by node Centered on, with radius A spherical region. All spatial points within this spherical region are affected by nodes. The characteristics of the nodes The characteristic information will be spatially diffused and mapped within this spherical domain.

[0144] For example, for a location at coordinate ( The high-gradient region node, after adjustment, has a radius of 3 meters and a spatial action sphere covering a volume of approximately 113 cubic meters; for nodes located at coordinates ( The low-gradient region nodes, after adjustment, have a radius of 10 meters, and their spatial interaction sphere covers a volume of approximately 4189 cubic meters. Within the monitoring area, the spatial interaction spheres of all graph nodes collectively cover the entire 3D space, with overlapping areas between the spheres of adjacent nodes. These overlapping areas represent the spatial extent of feature interaction among multiple nodes and are key regions for achieving multi-scale feature fusion.

[0145] For example, in one implementation, 48 graph nodes formed 127 overlapping regions in 3D space, with each overlapping region involving 2 to 5 nodes. This step established the spatial scope of the graph node features, which adaptively adjusts according to the feature energy gradient of the region where the node is located, providing a geometric basis for subsequent feature space mapping and fusion.

[0146] After dividing the spatial action sphere, the high-frequency acoustic emission features and low-frequency microseismic features corresponding to each graph node are extracted and mapped to virtual sound pressure scalar field and virtual work vector field within the spatial action sphere, respectively. For graph nodes... The multi-scale feature vector is extracted, which contains high-frequency energy features after spatiotemporal alignment and residual time delay compensation. and low-frequency energy characteristics The high-frequency energy characteristics reflect the energy intensity of the acoustic emission signal detected by the node, and are a scalar value; the low-frequency energy characteristics reflect the energy intensity of the microseismic signal detected by the node, and are also a scalar value.

[0147] To represent the spatial distribution of these features in three-dimensional space, high-frequency features are... Mapped to a sphere of spatial action Virtual sound pressure scalar field within low-frequency features Mapped to virtual work vector field The virtual sound pressure scalar field represents the spatial distribution of high-frequency energy, and is mapped using radial basis functions. The mapping formula is as follows: ,in For spatial points Distance to the node center The parameter used to control the decay rate is usually set to... This causes the sound pressure at the boundary of the sphere to decrease to about 5% of the central value.

[0148] The virtual work vector field represents the spatial propagation direction and intensity of low-frequency energy. Its direction points towards the node center, and its amplitude decreases with distance. The mapping formula is: ,in Let be the vector pointing from a point in space to the center of the node. The parameter used to control the decay rate is usually set to... .

[0149] For example, for high-frequency eigenvalues ,radius The virtual sound pressure at a node 1 meter away from the center is... At a distance of 2 meters from the center This step transforms discrete node features into a continuous spatial field distribution. This mapping method fully considers the spatial decay characteristics of feature energy, providing a spatialized feature representation for subsequent physical field interaction calculations.

[0150] After completing the virtual physics mapping, gradient divergence cross-coupling operations are performed on the virtual sound pressure scalar field and the virtual work vector field within the overlapping region of the spatial action spheres of adjacent graph nodes to generate a local spatiotemporal coupling tensor. The overlapping region is the intersection of the spatial action spheres of multiple nodes, and multiple virtual physics fields mapped by these nodes exist simultaneously within this region.

[0151] For any spatial point within the overlapping region First, identify all graph nodes that cover the point, and denote them as the node set. ,in The number of nodes covering this point is typically 2 to 5. For each node... Calculate the position of the node in space. Virtual sound pressure scalar field value at the location and its spatial gradient Simultaneously calculate the virtual work vector field and its divergence .

[0152] gradient It is a three-dimensional vector that reflects the direction and rate of spatial variation of high-frequency energy at that point; divergence It is a scalar that reflects the degree of convergence or divergence of low-frequency energy at that point. Gradient divergence cross-coupling operations are performed to generate a local spatiotemporal coupled tensor. ,in Representing the tensor cross product, it transforms a three-dimensional vector... With scalar Perform the outer product operation to obtain a The tensor is then summed over the tensors of all covered nodes.

[0153] The physical meaning of this tensor is the coupling strength between the spatial gradient of high-frequency energy and the spatial convergence / divergence of low-frequency energy, reflecting the degree of physical interaction between features of different scales at this spatial point. For example, at a spatial point in an overlapping region ( At point ), there are 3 nodes covering the area, and the gradient of node 1 is ( The divergence of node 1 is 0.35; the gradient of node 2 is ( The divergence of node 3 is 0.28; the gradient of node 3 is ( The divergence is 0.42. Through this step, a physical coupling relationship between high-frequency and low-frequency features is established. This coupling operation fully considers the physical propagation mechanism and energy interaction characteristics of signals at different scales, providing a physically meaningful tensor expression for subsequent feature intensity quantization.

[0154] After generating the local spatiotemporal coupling tensor, the eigennorm of the tensor in three-dimensional space is calculated to quantitatively characterize the physical interaction strength of features at different scales within the overlapping space. The eigennorm is a measure of a tensor, reflecting its overall size or energy. For the local spatiotemporal coupling tensor... The calculation is performed using a preset norm, and the formula is as follows: ,in For the tensor's first Each component.

[0155] The predefined norm is the square root of the sum of squares of all elements in the tensor. It is a non-negative scalar; a larger value indicates a stronger multi-scale feature interaction at that spatial point, suggesting it is closer to the disaster source; a smaller value indicates a weaker interaction and greater distance from the disaster source. For example, for the aforementioned spatial point ( The tensor at () has a preset norm of 0.251.

[0156] In the three-dimensional space of the entire monitoring area, the characteristic norm is calculated for spatial points within all overlapping regions, forming a discrete characteristic norm distribution. Regions with higher characteristic norms are typically concentrated near the disaster source because the high-frequency and low-frequency signal energy at the source is strong, and the spatial gradient and divergence are large, resulting in a larger norm for the coupling tensor. For example, in one calculation, the maximum characteristic norm of 0.856 appeared at coordinates (…). The location of the disaster source closely matches the location identified later. This step transforms the complex tensor representation into a concise scalar index, which quantitatively reflects the physical interaction strength of multi-scale features, providing core data for subsequent spatial interpolation and fusion feature map generation.

[0157] After calculating the characteristic norm, radial basis function interpolation is performed along the spatial geometric path of the overlapping region to generate a continuous characteristic energy distribution field. Since the characteristic norm is defined only at discrete spatial points within the overlapping region, while the fused feature map needs to be continuously defined on a three-dimensional spatial grid across the entire monitoring area, spatial interpolation is necessary. Radial basis function interpolation is a commonly used multidimensional spatial interpolation method that can generate smooth interpolation surfaces between irregularly distributed data points.

[0158] For any spatial point within the monitoring area Its characteristic energy value is calculated by radial basis function interpolation as follows: ,in The number of sampling points with known feature norms. For the first The coordinates of each sampling point This represents the weighting coefficient corresponding to that sampling point. Radial basis functions, commonly used forms include Gaussian functions. or multiple quadratic functions , The shape parameter controls the decay rate of the function.

[0159] Weighting coefficient The interpolation function is determined by solving a system of linear equations, ensuring that it is exactly equal to the known characteristic norm at all sampling points. For example, in a certain interpolation, a Gaussian radial basis function is used, and the shape parameter... The feature norms of the center points of 127 overlapping regions are interpolated to generate a continuous feature energy distribution field covering the entire monitoring area. The interpolated energy distribution field is exactly equal to the calculated feature norm at the sampling points, and transitions smoothly between sampling points, avoiding abrupt changes and discontinuities. Through this step, the discrete feature norm is extended into a continuous spatial field. This interpolation method ensures the spatial smoothness and physical rationality of the energy distribution, providing a continuous data foundation for subsequent filtering processing and feature map fusion output.

[0160] After generating a continuous feature energy distribution field, bilateral filtering is applied to the feature energy distribution field to output the final fused feature map. Bilateral filtering is an edge-preserving smoothing filtering method that can smooth noise while retaining the edge features of energy abrupt changes. It is particularly suitable for disaster source location scenarios because the source location usually corresponds to abrupt energy jumps.

[0161] Bilateral filtering considers both spatial distance and energy difference; the filtering formula is as follows: ,in For The filter window is centered on the center, and the window size is usually set to 5 to 10 meters; For spatial domain Gaussian functions, , The standard deviation of the spatial domain controls the degree of spatial smoothness and is usually set to 2 to 3 meters. It is a Gaussian function in the energy field. , The standard deviation of the energy range controls the tolerance for energy differences, and is typically set to 10% to 20% of the characteristic energy range. This is the normalization coefficient.

[0162] The core idea of ​​bilateral filtering is to assign greater weight to points that are spatially close and have similar energy values ​​for smoothing, while assigning less weight to points that are spatially close but have significantly different energy values ​​to avoid smoothing across energy boundaries. For example, near the source of a disaster, the energy rapidly decays from a peak of 0.856 to the surrounding area of ​​0.35. Bilateral filtering will preserve this energy boundary and will not perform cross-boundary smoothing, thus accurately preserving the peak characteristics of the source location. In the far-field flat region, the energy value changes slowly between 0.10 and 0.15. Bilateral filtering will sufficiently smooth this area and eliminate noise fluctuations.

[0163] After bilateral filtering, the fused feature map is output. This feature map is a scalar field defined on a three-dimensional spatial grid, where the value of each grid point reflects the likelihood of that location being a source of a hazard. For example, in a certain processing step, the fused feature map is located at coordinates ( The peak value of 0.842 was reached at a certain location, which was identified as the early source of the disaster. This step completed the final generation of the fused feature map. Bilateral filtering effectively suppressed noise interference while preserving the peak features of the source, significantly improving the quality of the fused feature map and the reliability of source localization.

[0164] This implementation achieves deep physical fusion of multi-scale features, effectively solving the problems of traditional simple fusion methods lacking physical mechanisms and failing to preserve edge features. By dividing the spatial action sphere, the spatial range of features is established, providing a geometric basis for the spatialized representation of features. By mapping high-frequency and low-frequency features to virtual sound pressure scalar fields and virtual work vector fields, the conversion from discrete features to a continuous spatial field is realized, providing a spatialized representation for physical field interaction calculations. A local spatiotemporal coupling tensor is generated through gradient divergence cross-coupling operations, establishing physical coupling relationships between features of different scales and fully considering the physical mechanisms of signal propagation. Furthermore, the physical interaction strength is quantitatively characterized by calculating the feature norm of the tensor, converting the complex tensor into a concise energy index. A continuous feature energy distribution field is generated through radial basis function interpolation, achieving a smooth conversion from discrete data to a continuous field. Finally, bilateral filtering smooths noise while preserving source peak features, significantly improving the quality of the fused feature map. This technical solution establishes a multi-scale feature fusion framework with clear physical meaning, significantly improving the fused feature map's ability to represent and locate disaster sources.

[0165] In one embodiment of this example, refer to Figure 4 Determining the early source location of a disaster based on the fused feature map includes the following steps: S810. Map the fused feature map onto the spatial grid of the monitoring area and calculate the feature energy value of each spatial grid point; S820. Calculate the gradient vector of the feature energy value in three-dimensional space and determine the temporal evolution direction of the multi-scale feature. S830. Determine the causal flow direction of each spatial grid point based on the temporal evolution direction; S840. Combine the gradient vector with the causal flow direction to obtain the comprehensive gradient. S850. Perform a reverse search in the opposite direction of the comprehensive gradient to determine the geometric center, and output the physical coordinates of the geometric center as the location of the early source of the disaster.

[0166] This implementation method establishes a source location mechanism based on spatiotemporal causal relationships by introducing methods such as energy gradient analysis, temporal evolution direction identification, and causal flow integration. This effectively solves the problem of insufficient anti-interference capability of traditional simple peak detection methods and significantly improves the accuracy and reliability of early disaster source location.

[0167] Specifically, the fused feature map is first mapped onto a spatial grid of the monitoring area, and the feature energy value of each spatial grid point is calculated. The fused feature map is a continuous three-dimensional scalar field generated through the aforementioned steps. The feature map is defined across the entire three-dimensional space of the monitoring area. To facilitate subsequent discretization calculations and source search, the continuous fused feature map needs to be discretized onto a regular spatial grid.

[0168] A three-dimensional spatial grid is established within the monitoring area, with the grid spacing typically set between 0.5 meters and 2 meters, determined based on the spatial scale of the monitoring area and the required positioning accuracy. For example, for a monitoring area of ​​100 meters × 80 meters × 60 meters, using a grid spacing of 1 meter, the total number of spatial grid points would be approximately 480,000. For each spatial grid point... Extract the feature energy value of the point from the fused feature map. If different spatial resolutions are used in the fused feature maps, the energy values ​​of the grid points are calculated using the trilinear interpolation method.

[0169] The characteristic energy value reflects the likelihood of a spatial location being a source of a disaster; the higher the energy value, the more likely the location is to be a source. For example, in a certain processing, spatial grid points ( The feature energy value of the fused feature map is 0.835, while the energy values ​​of the surrounding grid points are all below 0.40, showing obvious energy peak characteristics. Through this step, the continuous fused feature map is converted into a discrete spatial grid representation. This discretization facilitates subsequent gradient calculation and back-search, providing a structured data foundation for source localization.

[0170] After obtaining the feature energy values ​​on the spatial grid, the gradient vector of the feature energy values ​​in three-dimensional space is calculated, and the temporal evolution direction of the multi-scale features is determined. For spatial grid points... Its energy gradient vector is calculated as follows In the discretization implementation, the central difference method is used to calculate the partial derivatives in each direction, for example... ,in This represents the grid spacing.

[0171] The energy gradient vector points in the direction of the fastest energy growth; physically, it points from a low-energy region to a high-energy region, i.e., from the far field to the source. The magnitude of the gradient vector... This reflects the severity of energy change at that location; the gradient magnitude is typically larger near the source. For example, at spatial grid points ( At point ), the energy gradient vector is ( The modulus is 0.28, while at the far-field grid points ( At position ), the gradient vector is ( The modulus is only 0.04.

[0172] Simultaneously, the temporal evolution direction of multi-scale features is determined, reflecting the evolutionary trend of features along the time axis. The temporal direction of energy propagation is determined by analyzing the chronological order of high-frequency and low-frequency feature flows along the time axis. For each spatial grid point, the timestamp sequence of the corresponding multi-scale features is extracted, and the mean and standard deviation of the timestamps are calculated. Earlier timestamps correspond to the starting point of energy propagation, and later timestamps correspond to the ending point. The temporal evolution direction is defined as pointing from a later timestamp to an earlier timestamp, i.e., in reverse time, towards the origin of energy.

[0173] For example, at grid points near the source, the average timestamp of multi-scale features is 15.21 seconds, while at grid points in the far field, the average timestamp is 15.35 seconds, with the temporal evolution direction pointing from the far field towards the source. This step obtains spatial gradient information and temporal evolution information of the energy. These two types of information reflect the location characteristics of the source from different dimensions, providing a dual basis for subsequent determination of causal flow direction.

[0174] After obtaining the energy gradient vector and temporal evolution direction, the causal flow direction of each spatial grid point is determined based on the temporal evolution direction. The causal flow direction reflects the causal relationship of energy propagation, i.e., the direction of propagation from the source outwards. Physically, the disaster source is the starting point of energy release; energy radiates and propagates from the source into the surrounding space, forming a causal flow direction from the source to the periphery. The temporal evolution direction already reflects this causal relationship, i.e., from a later position to an earlier position, corresponding to the direction of tracing back from the receiving point to the transmitting point. Therefore, the causal flow direction can be directly set as the temporal evolution direction, denoted as […]. .

[0175] To improve the accuracy of causal flow, corrections are made using the energy gradient vector. If the energy gradient vector... With the direction of temporal evolution If the angle between the two is less than 90 degrees, it indicates that they point in the same direction, and the causal flow remains the direction of temporal evolution. If the angle is greater than 90 degrees, it indicates that there is a conflict between the two. A weighted average method is used to combine the two, and the corrected causal flow is: ,in and The weighting coefficient is usually set to... , This gives higher weight to time-series information because time-series causal relationships are more physically reliable.

[0176] For example, at a certain grid point, the temporal evolution direction is (0.6, 0.8, 0.0), and the energy gradient vector is ( The angle between the two is approximately 10 degrees, and the causal flow direction remains the temporal evolution direction; while at another grid point, the temporal evolution direction is (0.6, 0.8, 0.0), and the energy gradient vector is ( The angle between the two is approximately 160 degrees, and the corrected causal flow direction is... This step establishes the causal flow direction for each spatial grid point. This flow direction integrates temporal evolution information and energy gradient information, accurately reflecting the causal relationship of energy propagation from the source outward, and providing a causal basis for subsequent comprehensive gradient calculations.

[0177] After determining the causal flow direction, the gradient vector and the causal flow direction are vector-synthesized to obtain the comprehensive gradient. The comprehensive gradient is a combined expression of the energy spatial gradient and the causal flow direction, considering both the spatial distribution characteristics of energy and the temporal causal relationship of energy propagation. For spatial grid points... Its comprehensive gradient calculation is as follows ,in and These are weighting coefficients that control the relative contributions of the energy gradient and causal flow direction to the overall gradient. The weighting coefficients need to be adjusted based on the specific application scenario; they are typically set to [value missing]. , We assign a slightly higher weight to the energy gradient because the energy gradient directly reflects the spatial location characteristics of the source.

[0178] Comprehensive gradient It is a three-dimensional vector pointing in the direction of the most likely source under the combined effects of energy and causality. For example, at a certain grid point, the energy gradient vector is ( The causal flow is ( The combined gradient is The direction of the integrated gradient points towards the source, and its magnitude reflects the certainty of the direction from that location to the source; a larger magnitude indicates a more definitive direction. This step achieves deep fusion of energy spatial information and temporal causal information. The integrated gradient fully utilizes multi-dimensional information, significantly improving the accuracy of source direction determination and providing a reliable guiding vector for subsequent back-search.

[0179] After obtaining the comprehensive gradient, a reverse search is performed in the opposite direction of the comprehensive gradient to determine the geometric center, and the physical coordinates of the geometric center are output as the early source location of the disaster. The comprehensive gradient points in the direction of the source, so the search is performed in the opposite direction of the comprehensive gradient, that is, tracing back from the high-energy region to the low-energy region, and from the later time position to the earlier time position, gradually approaching the origin point of energy.

[0180] The reverse search employs an iterative algorithm. It first takes the spatial grid point with the highest feature energy value in the fused feature map as the starting point for the search, denoted as . Starting from the search origin, in each iteration, the comprehensive gradient of the current grid point is calculated. It moves to the next grid point in the opposite direction of the combined gradient. The direction of movement is... The movement step size is set to 0.5 to 1 times the grid spacing, for example, 0.5 meters. The coordinates of the next grid point are calculated as follows: ,in This is the step size parameter.

[0181] After each iteration, determine whether the comprehensive gradient magnitude of the current grid point is less than the preset convergence threshold. This threshold is typically set to 2 to 3 times the background gradient magnitude, for example, 0.05. If This indicates that the gradient minimum region has been reached, i.e., the origin of energy. The current grid point is then determined as the geometric center, and the iteration terminates. Continue iterative searching.

[0182] To avoid getting trapped in local optima or saddle points, a momentum term is introduced during the iteration process, i.e. ,in This is the momentum coefficient, usually set to 0.2 to 0.3, to maintain a certain inertia in the search path and avoid sawtooth-like oscillations.

[0183] For example, in a search, starting from the starting point ( Starting from ( ), after 18 iterations, it reaches the coordinates ( ). The combined gradient magnitude at this location is 0.032, which is less than the convergence threshold of 0.05, thus it is determined to be the geometric center. Output the physical coordinates of the geometric center ( This location serves as the early source of the disaster. The error between this location and the confirmed rockburst lesion location from subsequent on-site investigation is only 2.1 meters, verifying the accuracy of the location method. Through this step, the precise determination of the early source location of the disaster was achieved. The reverse search method fully utilizes the guiding effect of the comprehensive gradient, effectively avoiding local peak interference and realizing an accurate mapping from energy distribution to the source location.

[0184] This implementation achieves precise early disaster source location based on spatiotemporal causality, effectively solving the problems of insufficient anti-interference capability and neglect of temporal causal information in traditional simple peak detection methods. By mapping the fused feature map to a spatial grid and calculating energy values, a discretized data foundation for source location is established. By calculating the energy gradient vector and determining the temporal evolution direction, feature information of the source location is extracted from both spatial and temporal dimensions. By determining the causal flow direction based on the temporal evolution direction, a causal relationship expression for energy propagation is established, fully utilizing temporal information. Furthermore, by vector synthesis of the energy gradient and the causal flow direction to obtain a comprehensive gradient, deep fusion of spatial and temporal information is achieved, significantly improving the accuracy of source direction determination. Finally, by performing a reverse search along the opposite direction of the comprehensive gradient and introducing a momentum term to avoid local extrema, precise determination of the early disaster source location is achieved. This technical solution significantly improves the accuracy and reliability of the location results, providing precise location information for early disaster warning and emergency response.

[0185] This application embodiment also provides a data processing host, including: The memory is configured to store instructions; and The processor is configured to retrieve the instructions from the memory and, when executing the instructions, to implement the aforementioned method for early source tracing of sudden disasters based on multi-scale feature fusion.

[0186] Reference Figure 2 This application also provides a disaster monitoring system, including: Acoustic emission sensor array; Micro-vibration sensor array; The data processing host is communicatively connected to both the acoustic emission sensor array and the micro-vibration sensor array.

[0187] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this application can take the form of a computer program product embodied on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0188] This application is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this application. It will be understood that each block of the flowchart illustrations and / or block diagrams, as well as combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create a machine for implementing the flowchart illustrations. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0189] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0190] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0191] In a typical configuration, a computing device includes one or more processors (CPU), input / output interfaces, network interfaces, and memory.

[0192] Memory may include non-persistent memory in computer-readable media, such as random access memory (RAM) and / or non-volatile memory, such as read-only memory (ROM) or flash RAM. Memory is an example of computer-readable media.

[0193] Computer-readable media includes both permanent and non-permanent, removable and non-removable media that can store information using any method or technology. Information can be computer-readable instructions, data structures, modules of programs, or other data. Examples of computer storage media include, but are not limited to, phase-change memory (PRAM), static random access memory (SRAM), dynamic random access memory (DRAM), other types of random access memory (RAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other memory technologies, CD-ROM, digital versatile optical disc (DVD) or other optical storage, magnetic tape, magnetic disk storage or other magnetic storage devices, or any other non-transferable medium that can be used to store information accessible by a computing device. As defined herein, computer-readable media does not include transient computer-readable media, such as modulated data signals and carrier waves.

[0194] It should also be noted that the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such process, method, article, or apparatus. Unless otherwise specified, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes that element.

[0195] The above are merely embodiments of this application and are not intended to limit the scope of this application. Various modifications and variations can be made to this application by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this application should be included within the scope of the claims of this application.

Claims

1. A method for early source tracing of sudden disasters based on multi-scale feature fusion, characterized in that, Applied to a disaster monitoring system, the disaster monitoring system includes an acoustic emission sensor array, a microseismic sensor array, and a data processing host. The acoustic emission sensor array and the microseismic sensor array are respectively communicatively connected to the data processing host. The method includes: Acquire the high-frequency signal output by the acoustic emission sensor array and the low-frequency signal output by the micro-vibration sensor array in the monitoring area; Short-time Fourier transform is performed on the high-frequency signal to obtain the high-frequency energy wave packet characteristic flow, and wavelet packet decomposition is performed on the low-frequency signal to obtain the low-frequency band energy characteristic flow. Based on a pre-built three-dimensional spatial database, spatial grid search points are set for high-frequency energy wave packet characteristic flow and low-frequency band energy characteristic flow, and the equivalent refractive index at each position on the path between the spatial grid search points and the sensor is obtained from the three-dimensional spatial database. Calculate the propagation path of sound waves in a non-homogeneous medium based on the equivalent refractive index; The propagation path is divided into multiple micro-segments, and the high-frequency equivalent propagation velocity and low-frequency equivalent propagation velocity corresponding to each micro-segment are obtained. Calculate the high-frequency propagation time and low-frequency propagation time of each micro-element, where the high-frequency propagation time is equal to the length of the micro-element divided by the high-frequency equivalent propagation speed, and the low-frequency propagation time is equal to the length of the micro-element divided by the low-frequency equivalent propagation speed. The high-frequency propagation time is obtained by summing the high-frequency propagation time of all infinitesimal segments, and the low-frequency propagation time is obtained by summing the low-frequency propagation time of all infinitesimal segments. Based on the high-frequency propagation time and the low-frequency propagation time, the high-frequency energy wave packet characteristic flow and the low-frequency band energy characteristic flow are reverse-shifted and compensated on the time axis to obtain the spatiotemporally aligned high-frequency characteristic flow and the spatiotemporally aligned low-frequency characteristic flow. Extract burst spike features from the high-frequency feature stream, construct a sliding time window centered on the timestamp of the burst spike features, and project the low-frequency feature stream into the sliding time window. Calculate the normalized cross-correlation function of the high-frequency feature flow and the low-frequency feature flow, determine the residual time delay based on the peak position of the normalized cross-correlation function, and fine-tune the low-frequency feature flow based on the residual time delay to obtain multi-scale features; Obtain the spatial coordinates of each sensor in the acoustic emission sensor array and the micro-vibration sensor array; A graph structure is constructed based on spatial coordinates, where multi-scale features serve as the initial vectors for the graph nodes in the graph structure. The attention coefficients between adjacent graph nodes are calculated based on the graph attention network, and the gradient of the attention coefficient of each graph node in space is calculated. The magnitude of the gradient of the attention coefficient is also calculated. Determine whether the magnitude of the gradient of the attention coefficient exceeds a preset gradient threshold; If the magnitude of the gradient of the attention coefficient exceeds the preset gradient threshold, the graph node is determined to be in the region of increased feature energy density, and the spatial receptive field radius of the graph node is shrunk to the preset minimum radius. If the magnitude of the gradient of the attention coefficient does not exceed the preset gradient threshold, the spatial receptive field radius of the graph node is kept at the preset standard radius. Multi-scale features are fused based on the adjusted spatial receptive field radius to output a fused feature map, and the early source location of the disaster is determined based on the fused feature map.

2. The method according to claim 1, characterized in that, Extracting burst spike features from high-frequency feature streams includes: Calculate the first-order time derivative of the high-frequency characteristic flow; Calculate the second time derivative of the first time derivative; Determine whether the second-order time derivative exceeds a preset mutation threshold; When the second time derivative exceeds a preset mutation threshold, the feature point at the corresponding time moment is determined to be a sudden spike feature.

3. The method according to claim 1, characterized in that, Calculate the normalized cross-correlation function of the high-frequency and low-frequency characteristic flows, and determine the residual time delay based on the peak position of the normalized cross-correlation function, including: Obtain spatiotemporally aligned high-frequency feature stream segments and spatiotemporally aligned low-frequency feature stream segments within a sliding time window; Normalize the high-frequency characteristic flow segments and the low-frequency characteristic flow segments respectively; Calculate the cross-correlation function between the normalized high-frequency characteristic flow segment and the normalized low-frequency characteristic flow segment, and determine the peak position of the cross-correlation function; The residual time delay of the spatiotemporally aligned low-frequency feature stream segment relative to the spatiotemporally aligned high-frequency feature stream segment is determined based on the peak position.

4. The method according to claim 1, characterized in that, Multi-scale features are fused based on the adjusted spatial receptive field radius, and a fused feature map is output, including: Using the spatial coordinates of each graph node as the center and the corresponding adjusted spatial receptive field radius as the boundary, the spatial action sphere of each graph node is divided in three-dimensional space. Extract the high-frequency acoustic emission features and low-frequency microseismic features corresponding to each graph node, and map them into a virtual acoustic pressure scalar field and a virtual work vector field in the spatial action sphere, respectively. Within the overlapping region of the spatial action sphere of adjacent graph nodes, gradient divergence cross-coupling operation is performed on the virtual sound pressure scalar field and the virtual work vector field to generate a local spatiotemporal coupling tensor. Calculate the characteristic norm of the local spatiotemporal coupling tensor in three-dimensional space to quantitatively characterize the physical interaction strength of features at different scales in the overlapping space; The feature norm is interpolated radially along the spatial geometric path of the overlapping region to generate a continuous feature energy distribution field. Bilateral filtering is applied to the characteristic energy distribution field to output a fused feature map.

5. The method according to claim 4, characterized in that, The location of the early source of the disaster is determined based on the fused feature map, including: The fused feature map is mapped onto a spatial grid of the monitoring area, and the feature energy value of each spatial grid point is calculated. Calculate the gradient vector of the feature energy value in three-dimensional space and determine the temporal evolution direction of the multi-scale features; The causal flow direction of each spatial grid point is determined based on the temporal evolution direction; The gradient vector is combined with the causal flow direction to obtain the comprehensive gradient. Perform a reverse search in the opposite direction of the comprehensive gradient to determine the geometric center, and output the physical coordinates of the geometric center as the location of the early source of the disaster.

6. A data processing host, characterized in that, include: The memory is configured to store instructions; as well as The processor is configured to retrieve the instructions from the memory and, when executing the instructions, to implement the early source tracing method for sudden disasters based on multi-scale feature fusion according to any one of claims 1 to 5.

7. A disaster monitoring system, applied to the early source tracing method for sudden disasters based on multi-scale feature fusion as described in any one of claims 1-5, characterized in that, include: Acoustic emission sensor array; Micro-vibration sensor array; The data processing host is communicatively connected to both the acoustic emission sensor array and the micro-vibration sensor array.

Images