A rock mass mesoscopic fissure testing method and system

Abstract

The present application provides a kind of rock mass mesoscopic fissure test method and system, including the surface image in the process of loading rock mass is collected, and the acoustic emission data set containing three-dimensional coordinates, time, energy, the surface discrete strain field is solved using digital image correlation method, and the strain concentration area is extracted by the second invariant of partial curvature and threshold value;With the surface principal strain of characteristic time as reference, the anisotropic space-time distance of acoustic emission event is calculated and clustered, and independent event cluster is divided, the strain concentration area is spatially paired with acoustic emission cluster, respectively calculate spatial form compliance, time evolution synchronization degree, with spatial form compliance is modulated by synchronization degree, and combined with the maximum principal strain value is conditionally weighted correction, obtain matching confidence, when confidence exceeds the discrimination threshold, determine that the paired object is real rock mass mesoscopic fissure, can realize fissure reliable identification by synchronously coupling surface deformation and internal damage.

Description

Technical Field

[0001] This application belongs to the field of testing, and in particular relates to a method and system for testing microscopic fractures in rock masses. Background Technology

[0002] Rock masses, serving as a crucial foundation for numerous major engineering projects, are widely present in underground space development, mineral resource extraction, water conservancy and hydropower projects, and transportation tunnels. Under continuous external loads, rock masses typically undergo an evolutionary process from the accumulation of microscopic damage to macroscopic instability and failure. Primary and newly formed microscopic fractures within the rock mass continuously emerge, expand, and interconnect, potentially forming macroscopic fracture surfaces. The deterioration of rock mass structure caused by the evolution of microscopic fractures is a significant factor inducing sudden engineering geological disasters such as rockbursts, collapses, and landslides. Therefore, high-precision testing and monitoring of microscopic fractures during rock mass loading is not only a crucial foundation for revealing the mechanisms of rock mechanical damage and fracture but also a key prerequisite for identifying precursors to engineering disasters and providing early warnings of safety risks.

[0003] Digital image correlation (DIR) and acoustic emission (AE) techniques have been widely applied to the comprehensive monitoring of rock mass fracture processes. The former can acquire the full-field strain distribution of the rock mass through surface image tracking, characterizing surface fracture evolution; the latter can identify transient elastic waves released from internal fractures, achieving three-dimensional spatial localization of internal damage events. Existing methods for testing micro-fractures in rock masses typically analyze surface strain data and internal AE data as independent monitoring information in isolation, or simply at the level of data overlay. Due to the complex nonlinear and asynchronous relationships between the apparent deformation characteristics of rock masses and the internal micro-fracture behavior in three-dimensional spatial distribution and temporal evolution, existing methods struggle to establish a stable and accurate correspondence between surface strain concentration areas and internal AE events in the spatiotemporal dimensions. The lack of deep fusion between internal and external monitoring data makes actual testing susceptible to interference from environmental noise and spurious discrete signals, hindering the reliable identification and extraction of true micro-fractures in the rock mass, thus limiting the accurate assessment of rock mass damage status. Therefore, a testing method that can fuse surface deformation information and internal fracture information is urgently needed to overcome the limitations of spatiotemporal isolation of multi-source monitoring data and achieve reliable identification of micro-fractures in rock masses. Summary of the Invention

[0004] To improve the accuracy and reliability of identifying the evolution of internal and surface fractures in rock masses, this invention provides a method and system for testing microscopic fractures in rock masses.

[0005] In a first aspect, the present invention provides a method for testing microscopic fractures in rock mass, comprising the following steps: A sequence of surface images and a dataset of internal acoustic emission events are obtained during the loading process of the rock mass. The dataset of internal acoustic emission events includes the three-dimensional coordinates, occurrence time and energy value data of each internal acoustic emission event. The surface image sequence is processed based on digital image correlation to obtain the discrete strain field of the rock mass surface. The second invariant of the skew curvature corresponding to the discrete strain field is calculated, and the region with high strain concentration is identified by setting a threshold. The principal strain direction of the rock mass surface at the characteristic loading moment is extracted as the reference direction. The anisotropic spatiotemporal distance of each internal acoustic emission event in three-dimensional space and time interval is calculated. Based on the anisotropic spatiotemporal distance, data clustering operation is performed on all internal acoustic emission events to form multiple independent acoustic emission event clusters. The region with highly concentrated strain is spatially paired with the multiple independent acoustic emission event clusters. The spatial morphological conformity of each pair to be matched is calculated, and the temporal evolution synchronization of each pair to be matched is calculated. The temporal evolution synchronization is used as a modulation coefficient to perform data augmentation calculation on the spatial morphological conformity to obtain the matching confidence. When the matching confidence is greater than a preset discrimination threshold, it is determined that the pair to be matched corresponds to a real microscopic fracture in the rock mass.

[0006] Further, the step of processing the surface image sequence based on digital image correlation to obtain the discrete strain field of the rock mass surface, and calculating the second invariant of the skewness curvature corresponding to the discrete strain field, includes: The first image of the rock mass surface before loading is set as the reference image, and the subsequent loaded images are matched by digital image correlation to obtain the discrete displacement field. The discrete displacement field is smoothed and differentiated using a local surface fitting algorithm to obtain a continuous strain tensor field. The second-order partial derivatives of each independent strain component of the strain tensor field with respect to the surface spatial coordinates are calculated to obtain the second-order curvature matrix corresponding to each independent strain component. The inner product of the partial curvature parts of each second-order curvature matrix is ​​summed to obtain the second invariant of partial curvature, which represents the degree of non-uniformity of local deformation.

[0007] Furthermore, the step of identifying regions with highly concentrated strain by setting a threshold includes: The distribution characteristics of the second invariant of skewness curvature are statistically analyzed across the entire field, and the mean and standard deviation of the invariant are calculated. The sum of the mean and the standard deviation multiplied by a preset factor is set as the judgment threshold; Extract the connected components where the second invariant of the bias curvature is greater than the determination threshold, and mark the connected components as the regions where the strain is highly concentrated.

[0008] Furthermore, using the principal strain direction of the rock mass surface at the time of feature loading as the reference direction, the anisotropic spatiotemporal distance of each of the internal acoustic emission events in three-dimensional space and time interval is calculated, including: A three-dimensional anisotropic local coordinate system is constructed by taking the principal strain direction of the rock mass surface as the first coordinate axis, the direction perpendicular to the principal strain direction of the rock mass surface as the second coordinate axis, and the rock mass depth direction as the third coordinate axis. Calculate the coordinate difference vector between any two internal acoustic emission events in the three-dimensional space under the three-dimensional anisotropic local coordinate system, and assign different distance penalty weights to the corresponding spatial components according to the preset principal deformation direction weight, surface lateral weight and depth direction weight. Calculate the time difference between the occurrence of the two internal acoustic emission events; The magnitude of the weighted three-dimensional spatial coordinate difference vector is normalized by spatial scale, and the time difference of occurrence is normalized by temporal scale. Then, the anisotropic spatiotemporal distance is obtained by linearly combining the preset spatial contribution coefficient and temporal contribution coefficient.

[0009] Furthermore, the step of performing data clustering operations on all the internal acoustic emission events based on the anisotropic spatiotemporal distance to form multiple independent acoustic emission event clusters includes: An acoustic emission event association graph is constructed using each of the aforementioned internal acoustic emission events as nodes and the anisotropic spatiotemporal distance as the edge weight. A density-based spatial clustering algorithm searches for core nodes in the acoustic emission event association graph and expands the clusters of adjacent nodes whose anisotropic spatiotemporal distance is less than a preset clustering radius, outputting the multiple independent acoustic emission event clusters.

[0010] Furthermore, the step of spatially pairing the region with the highly concentrated strain with the multiple independent acoustic emission event clusters, and calculating the spatial morphological conformity of each pair to be matched, includes: The geometric center and the first principal axis direction of the region with highly concentrated strain are extracted respectively, as well as the projection centroid of the event location point cloud on the rock surface and the projection of the first principal axis direction obtained based on principal component analysis in the multiple independent acoustic emission event clusters. Calculate the Euclidean distance between the geometric center and the projected centroid, and the spatial angle between the first principal axis direction of the region with high strain concentration and the projection of the first principal axis direction; After dividing the Euclidean distance and the spatial angle by preset reference distance and reference angle respectively to perform dimensionless processing, the results are input into a preset negative exponential decay evaluation function for calculation, to obtain the spatial morphology conformity with a value between 0 and 1.

[0011] Further, the calculation of the temporal evolution synchronization degree of each pair to be matched includes: Extract the effective strain rate of change over time sequence within the region of high strain concentration, and the cumulative acoustic emission energy rate of change over time sequence within the multiple independent acoustic emission event clusters; A sliding time window is set, and the two rate of change sequences are resampled and aligned within the sliding time window. Then, a fast Fourier transform is performed to convert them to the frequency domain. Calculate the cross power spectrum of the two rate-of-change sequences in the frequency domain, and obtain the time-delay cross-correlation function through inverse Fourier transform; The maximum positive correlation peak of the time-delay cross-correlation function is extracted within the preset allowable time delay range. When both rate of change sequences meet the preset effective fluctuation conditions, the maximum positive correlation peak is normalized and subjected to amplitude limiting processing as the time evolution synchronization degree. When any rate of change sequence does not meet the preset effective fluctuation condition, the time evolution synchronization degree of the corresponding window is set to 0 or the corresponding window is skipped.

[0012] Further, the step of using the temporal evolution synchronization degree as a modulation coefficient to perform data augmentation calculation on the spatial morphological conformity to obtain the matching confidence degree includes: The time evolution synchronization degree is used as a modulation coefficient with a value between 0 and 1; Multiplying the spatial morphology conformity by the modulation coefficient yields the preliminary confidence level; When both the spatial morphology conformity and the temporal evolution synchronization meet the corresponding minimum threshold conditions, the initial confidence level is weighted and corrected using the maximum principal strain value of the strain-concentrated region after normalization, and the correction result is limited so that the output matching confidence level is kept between 0 and 1. When the spatial morphology conformity or the temporal evolution synchronization does not meet the corresponding minimum threshold condition, the preliminary confidence level is used as the matching confidence level.

[0013] Secondly, the present invention also provides a rock mass micro-fracture testing system, comprising the following modules: The calculation module is used to acquire surface image sequences and internal acoustic emission event datasets during the loading process of the rock mass. The internal acoustic emission event dataset contains the three-dimensional coordinates, occurrence time and energy value data of each internal acoustic emission event. The surface image sequence is processed based on digital image correlation to obtain the discrete strain field of the rock mass surface. The second invariant of the eccentric curvature corresponding to the discrete strain field is calculated, and the region with high strain concentration is identified by setting a threshold. The extraction module is used to extract the principal strain direction of the rock mass surface at the characteristic loading time as the reference direction, calculate the anisotropic spatiotemporal distance of each internal acoustic emission event in three-dimensional space and time interval, and perform data clustering operation on all internal acoustic emission events according to the anisotropic spatiotemporal distance to form multiple independent acoustic emission event clusters. The determination module is used to spatially pair the region with the highly concentrated strain with the multiple independent acoustic emission event clusters, calculate the spatial morphological conformity of each pair to be matched, calculate the temporal evolution synchronization of each pair to be matched, use the temporal evolution synchronization as a modulation coefficient to perform data augmentation calculation on the spatial morphological conformity to obtain the matching confidence. When the matching confidence is greater than a preset discrimination threshold, it is determined that the pair to be matched corresponds to a real microscopic fracture in the rock mass.

[0014] Preferably, the step of processing the surface image sequence based on digital image correlation to obtain the discrete strain field of the rock mass surface, and calculating the second invariant of the skew curvature corresponding to the discrete strain field, includes: The first image of the rock mass surface before loading is set as the reference image, and the subsequent loaded images are matched by digital image correlation to obtain the discrete displacement field. The discrete displacement field is smoothed and differentiated using a local surface fitting algorithm to obtain a continuous strain tensor field. The second-order partial derivatives of each independent strain component of the strain tensor field with respect to the surface spatial coordinates are calculated to obtain the second-order curvature matrix corresponding to each independent strain component. The inner product of the partial curvature parts of each second-order curvature matrix is ​​summed to obtain the second invariant of partial curvature, which represents the degree of non-uniformity of local deformation.

[0015] Preferably, identifying regions of high strain concentration by setting a threshold includes: The distribution characteristics of the second invariant of skewness curvature are statistically analyzed across the entire field, and the mean and standard deviation of the invariant are calculated. The sum of the mean and the standard deviation multiplied by a preset factor is set as the judgment threshold; Extract the connected components where the second invariant of the bias curvature is greater than the determination threshold, and mark the connected components as the regions where the strain is highly concentrated.

[0016] Preferably, the extraction of the principal strain direction of the rock mass surface at the time of characteristic loading is used as the reference direction to calculate the anisotropic spatiotemporal distance of each of the internal acoustic emission events in three-dimensional space and time interval, including: A three-dimensional anisotropic local coordinate system is constructed by taking the principal strain direction of the rock mass surface as the first coordinate axis, the direction perpendicular to the principal strain direction of the rock mass surface as the second coordinate axis, and the rock mass depth direction as the third coordinate axis. Calculate the coordinate difference vector between any two internal acoustic emission events in the three-dimensional space under the three-dimensional anisotropic local coordinate system, and assign different distance penalty weights to the corresponding spatial components according to the preset principal deformation direction weight, surface lateral weight and depth direction weight. Calculate the time difference between the occurrence of the two internal acoustic emission events; The magnitude of the weighted three-dimensional spatial coordinate difference vector is normalized by spatial scale, and the time difference of occurrence is normalized by temporal scale. Then, the anisotropic spatiotemporal distance is obtained by linearly combining the preset spatial contribution coefficient and temporal contribution coefficient.

[0017] Preferably, the step of performing data clustering operation on all the internal acoustic emission events based on the anisotropic spatiotemporal distance to form multiple independent acoustic emission event clusters includes: An acoustic emission event association graph is constructed using each of the aforementioned internal acoustic emission events as nodes and the anisotropic spatiotemporal distance as the edge weight. A density-based spatial clustering algorithm searches for core nodes in the acoustic emission event association graph and expands the clusters of adjacent nodes whose anisotropic spatiotemporal distance is less than a preset clustering radius, outputting the multiple independent acoustic emission event clusters.

[0018] Preferably, the step of spatially pairing the region with the highly concentrated strain with the plurality of independent acoustic emission event clusters, and calculating the spatial morphological conformity of each pair to be matched, includes: The geometric center and the first principal axis direction of the region with highly concentrated strain are extracted respectively, as well as the projection centroid of the event location point cloud on the rock surface and the projection of the first principal axis direction obtained based on principal component analysis in the multiple independent acoustic emission event clusters. Calculate the Euclidean distance between the geometric center and the projected centroid, and the spatial angle between the first principal axis direction of the region with high strain concentration and the projection of the first principal axis direction; After dividing the Euclidean distance and the spatial angle by preset reference distance and reference angle respectively to perform dimensionless processing, the results are input into a preset negative exponential decay evaluation function for calculation, to obtain the spatial morphology conformity with a value between 0 and 1.

[0019] Preferably, the calculation of the temporal evolution synchronization degree of each pair to be matched includes: Extract the effective strain rate of change over time sequence within the region of high strain concentration, and the cumulative acoustic emission energy rate of change over time sequence within the multiple independent acoustic emission event clusters; A sliding time window is set, and the two rate of change sequences are resampled and aligned within the sliding time window. Then, a fast Fourier transform is performed to convert them to the frequency domain. Calculate the cross power spectrum of the two rate-of-change sequences in the frequency domain, and obtain the time-delay cross-correlation function through inverse Fourier transform; The maximum positive correlation peak of the time-delay cross-correlation function is extracted within the preset allowable time delay range. When both rate of change sequences meet the preset effective fluctuation conditions, the maximum positive correlation peak is normalized and subjected to amplitude limiting processing as the time evolution synchronization degree. When any rate of change sequence does not meet the preset effective fluctuation condition, the time evolution synchronization degree of the corresponding window is set to 0 or the corresponding window is skipped.

[0020] Preferably, the step of using the temporal evolution synchronization degree as a modulation coefficient to perform data augmentation calculation on the spatial morphological conformity degree to obtain the matching confidence degree includes: The time evolution synchronization degree is used as a modulation coefficient with a value between 0 and 1; Multiplying the spatial morphology conformity by the modulation coefficient yields the preliminary confidence level; When both the spatial morphology conformity and the temporal evolution synchronization meet the corresponding minimum threshold conditions, the initial confidence level is weighted and corrected using the maximum principal strain value of the strain-concentrated region after normalization, and the correction result is limited so that the output matching confidence level is kept between 0 and 1. When the spatial morphology conformity or the temporal evolution synchronization does not meet the corresponding minimum threshold condition, the preliminary confidence level is used as the matching confidence level.

[0021] This invention effectively weakens irrelevant deformation interference and locates potential fracture initiation sites by fusing surface image sequences with internal acoustic emission event data and utilizing the second invariant of skew curvature to identify areas of high strain concentration. The anisotropic spatiotemporal distance of acoustic emission events is calculated and clustered based on the principal strain direction of the surface. This mechanism conforms to the directional and temporal characteristics of rock mass fracture under load, improving the accuracy of acoustic emission event classification and enabling independent acoustic emission clusters to more reliably characterize the development and propagation of individual fractures. In the matching stage, spatial morphological conformity and temporal evolution synchronization are integrated. The temporal evolution synchronization modulates the spatial morphological conformity, and conditional weighting correction is performed using the maximum principal strain value to obtain the matching confidence level. This processing achieves dual-dimensional verification of spatial morphology and temporal evolution, reducing the risk of misjudgment due to accidental spatial overlap, thereby achieving high-precision quantitative identification and extraction of true microscopic fractures in rock mass. Attached Figure Description

[0022] Figure 1A flowchart of a method for testing microscopic fractures in rock mass; Figure 2 This is a schematic diagram of the distribution of the second invariant; Figure 3 This is a schematic diagram of the projected centroid of the acoustic emission event cluster and the direction of the first principal axis. Detailed Implementation

[0023] To make the purpose, features, and advantages of this application more apparent and understandable, specific embodiments and accompanying drawings will be used to clearly and completely describe the technical solution protected by this application. Obviously, the embodiments described below are only some embodiments of this application, and not all embodiments. Based on the embodiments in this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0024] Reference Figure 1 This invention provides a method for testing microscopic fractures in rock masses, comprising the following steps: S1. Obtain the surface image sequence and internal acoustic emission event dataset of the rock mass during the loading process. The internal acoustic emission event dataset includes the three-dimensional coordinates, occurrence time and energy value data of each internal acoustic emission event. Process the surface image sequence based on digital image correlation to obtain the discrete strain field of the rock mass surface, calculate the second invariant of the skew curvature corresponding to the discrete strain field, and identify the region with high strain concentration by setting a threshold.

[0025] A high-resolution charge-coupled device (CCD) camera was used to acquire surface image sequences of the rock mass during uniaxial compression at a preset sampling frequency. Simultaneously, a network of acoustic emission sensors distributed on the rock mass surface, along with a multi-channel acoustic emission acquisition instrument, was used to record elastic wave signals in real time, calculating the three-dimensional coordinates, occurrence time, and energy values ​​of internal acoustic emission events. The surface image sequence was read, and a zero-mean normalized cross-correlation algorithm was used to calculate the matching degree of image sub-regions before and after deformation. The pixel displacement was then solved using the inverse combined Gauss-Newton iteration method. The displacement field was smoothed and differentiated to obtain the strain tensor field of the rock mass surface. The second-order partial derivatives of each independent strain component with respect to the surface spatial coordinates x and y were calculated, and a corresponding second-order curvature matrix was constructed for each independent strain component. The spherical portion formed by the trace average term was subtracted from each second-order curvature matrix to obtain the corresponding partial curvature matrix. The inner product of each partial curvature matrix with itself was accumulated as the second invariant of the partial curvature. Therefore, the object of calculation for the second invariant of the bias curvature is unified as the set of second-order curvature matrices of each independent component of the strain tensor field, rather than the Hessian matrix or principal curvature of a single scalar field.

[0026] During threshold recognition, the second invariant matrix of skew curvature is traversed, and the mean and standard deviation of the entire field are calculated. The sum of the mean and standard deviation multiplied by a preset factor is set as the judgment threshold. Pixel connected regions exceeding the judgment threshold are marked as regions with highly concentrated strain. The distribution of the second invariant is as follows: Figure 2 As shown.

[0027] In one possible embodiment, the step of processing the surface image sequence based on digital image correlation to obtain the discrete strain field of the rock mass surface, and calculating the second invariant of the skew curvature corresponding to the discrete strain field, includes: The first image of the rock mass surface before loading is set as the reference image, and the subsequent loaded images are matched by digital image correlation to obtain the discrete displacement field. The discrete displacement field is smoothed and differentiated using a local surface fitting algorithm to obtain a continuous strain tensor field. The second-order partial derivatives of each independent strain component of the strain tensor field with respect to the surface spatial coordinates are calculated to obtain the second-order curvature matrix corresponding to each independent strain component. The inner product of the partial curvature parts of each second-order curvature matrix is ​​summed to obtain the second invariant of partial curvature, which represents the degree of non-uniformity of local deformation.

[0028] In the digital image correlation matching process, the image sub-region size is set to 31×31 pixels, and the calculation step size is 7 pixels to meet the displacement field resolution requirements. A 7×7 window quadratic or cubic polynomial local surface fitting algorithm is used to smooth and filter the discrete displacement field and differentiate it to generate a continuous strain tensor field. The second-order partial derivatives of each independent component of the strain tensor with respect to spatial coordinates x and y are calculated to obtain the second-order curvature matrix corresponding to each independent strain component. For each second-order curvature matrix, the spherical and partial components are separated, and the inner product of each partial curvature matrix with itself is accumulated to obtain the second invariant of the partial curvature. Specifically, for the q-th independent strain component, its second-order curvature matrix is ​​denoted as... The ball's mass is recorded as The biased part is denoted as Then the second invariant of the deviator curvature can be expressed as each The sum of its product with itself, or a proportional form of that sum; where q represents the index of the independent strain component, n is the dimension of the second-order curvature matrix, and I is the identity matrix. For example, when the surface strain distribution is in... When the magnitude is large, the magnitude of the second invariant of the partial curvature, obtained by second derivative, varies with the surface coordinate scale, sampling interval, smoothing window parameters, and second-order difference step size. Therefore, the judgment threshold can be adaptively determined by the statistical mean and standard deviation of the entire field. This invariant can amplify the local fluctuation characteristics of the strain field, reduce the interference of overall rigid body displacement and uniform deformation, and locate the small non-uniform deformation area caused by the initiation of microcracks in the rock mass.

[0029] In one possible embodiment, identifying regions of high strain concentration by setting a threshold includes: The distribution characteristics of the second invariant of skewness curvature are statistically analyzed across the entire field, and the mean and standard deviation of the invariant are calculated. The sum of the mean and the standard deviation multiplied by a preset factor is set as the judgment threshold; Extract the connected components where the second invariant of the bias curvature is greater than the determination threshold, and mark the connected components as the regions where the strain is highly concentrated.

[0030] The statistical analysis of the overall distribution characteristics requires traversing all calculation nodes in the discrete strain field to calculate the global statistical mean and standard deviation of the second invariant of eccentric curvature. The judgment threshold is set as the sum of the mean and standard deviation multiplied by a preset multiple, where the preferred range for the preset multiple is 1.5 to 3.0. For example, when the rock mass enters the nonlinear yielding stage, the calculated mean of the second invariant of eccentric curvature is... The standard deviation is If the preset multiplier is 2.0, then the threshold value is set to... When extracting connected components, an 8-neighborhood connected component labeling algorithm is applied to scan the full-field binarized image. To remove random noise caused by matching errors in digital image correlation methods, an area filtering threshold is set to automatically filter out tiny regions with a connected area of ​​less than 50 pixels. The retained continuous large-area high-strain gradient connected components are labeled and output as highly concentrated strain regions where microcracks may initiate and expand, serving as feature sources for subsequent spatial matching.

[0031] S2, extract the principal strain direction of the rock mass surface at the characteristic loading moment as the reference direction, calculate the anisotropic spatiotemporal distance of each internal acoustic emission event in three-dimensional space and time interval, and perform data clustering operation on all internal acoustic emission events according to the anisotropic spatiotemporal distance to form multiple independent acoustic emission event clusters.

[0032] Before calculating the anisotropic spatiotemporal distance, based on the pre-loading calibration results, the acoustic emission three-dimensional positioning coordinates are uniformly transformed to the specimen coordinate system consistent with the physical coordinates of the surface using the digital image correlation method. The two in-surface directions of the specimen coordinate system correspond to the physical coordinate directions of the image, and the depth direction is perpendicular to the rock mass observation surface. The moment when the rock mass stress-strain curve reaches its peak strength is selected as the characteristic loading moment. If multiple regions with highly concentrated strain exist at this moment, the strain tensors of the calculation nodes within each highly concentrated strain region are extracted, and a joint weighted average is performed according to the area of ​​each region and the peak value of the second invariant of the skew curvature, resulting in a full-field weighted average strain tensor characterizing the overall deformation trend at the current loading stage. The full-field weighted average strain tensor is then subjected to eigenvalue decomposition, and the eigenvector corresponding to the maximum principal strain is determined as the principal strain direction on the rock mass surface. The joint weight of each highly concentrated strain region is obtained by normalizing the product of the region's area and the peak value of the second invariant of the skew curvature. A three-dimensional anisotropic local coordinate system is constructed by taking the principal strain direction of the rock mass surface as the first coordinate axis, the direction perpendicular to the principal strain direction of the rock mass surface as the second coordinate axis, and the rock mass depth direction as the third coordinate axis.

[0033] For any two internal acoustic emission events, their three-dimensional coordinate difference vectors are first transformed to the three-dimensional anisotropic local coordinate system to obtain the coordinate difference components along the first, second, and third coordinate axes. Based on preset weights for the main deformation direction, surface transverse direction, and depth direction, the distance penalty weights for each spatial component are determined, and the weighted spatial distance is calculated. Since acoustic emission events corresponding to the same fracture are usually continuously distributed along the dominant deformation-related direction or its adjacent direction, the distance penalty weight in the first coordinate axis direction is lower than the distance penalty weights in the surface transverse and depth directions to avoid incorrectly splitting acoustic emission events extending along the main fracture direction into different clusters. Simultaneously, the absolute difference between the occurrence times of the two events is calculated as the time interval. The distance penalty weights in each direction are used to calculate the weighted spatial distance, and the weighted spatial distance is divided by a preset reference spatial scale for spatial normalization; the time interval is divided by a preset reference time scale for temporal normalization; and then a linear weighted sum is performed according to the set spatial contribution coefficient and temporal contribution coefficient to obtain the anisotropic spatiotemporal distance. The DBSCAN density clustering algorithm is invoked, taking all internal acoustic emission events as the input dataset, using the anisotropic spatiotemporal distance as the distance metric function, setting the neighborhood radius and the minimum number of core points, performing data clustering operations, and removing isolated noise points to form multiple independent acoustic emission event clusters.

[0034] In one possible embodiment, the extraction of the principal strain direction of the rock mass surface at the characteristic loading moment is used as the reference direction, and the calculation of the anisotropic spatiotemporal distance of each of the internal acoustic emission events in three-dimensional space and time interval includes: A three-dimensional anisotropic local coordinate system is constructed by taking the principal strain direction of the rock mass surface as the first coordinate axis, the direction perpendicular to the principal strain direction of the rock mass surface as the second coordinate axis, and the rock mass depth direction as the third coordinate axis. Calculate the coordinate difference vector between any two internal acoustic emission events in the three-dimensional space under the three-dimensional anisotropic local coordinate system, and assign different distance penalty weights to the corresponding spatial components according to the preset principal deformation direction weight, surface lateral weight and depth direction weight. Calculate the time difference between the occurrence of the two internal acoustic emission events; The magnitude of the weighted three-dimensional spatial coordinate difference vector is normalized by spatial scale, and the time difference of occurrence is normalized by temporal scale. Then, the anisotropic spatiotemporal distance is obtained by linearly combining the preset spatial contribution coefficient and temporal contribution coefficient.

[0035] When constructing the local coordinate system, the direction of maximum principal strain on the surface of the loaded rock mass is set as the X' axis, the direction perpendicular to the X' axis within the rock mass surface is set as the Y' axis, and the depth direction perpendicular to the surface is set as the Z' axis. For the coordinate difference vector, distance penalty weights are set according to the deformation characteristics of the rock mass in each axis. The distance penalty weight for the direction of maximum principal strain is preferably 0.2, and the distance penalty weights for the lateral and depth directions of the surface are both preferably 0.4, so that acoustic emission events distributed along the dominant deformation-related direction or its adjacent direction maintain a small distance penalty, while giving higher penalties to lateral and depth deviations. The distance penalty weights are used to calculate the weighted spatial distance of the three-dimensional spatial coordinate difference vector; the spatial contribution coefficient is used to adjust the relative contribution between spatial distance and time interval in the anisotropic spatiotemporal distance calculation, and the two are weight parameters at different levels. The preferred range of the time contribution coefficient is 0.3 to 0.5, for example, set to 0.4. During the calculation, spatial normalization is achieved by dividing the square root of the weighted sum of squares in each direction by the maximum characteristic size of the rock mass, such as 100 mm; temporal normalization is achieved by dividing the absolute value of the time difference of occurrence by the total loading duration, such as 3600 s. For example, the normalized weighted spatial distance between two acoustic emission events is 0.05, and the normalized time difference is 0.02. Substituting the spatial contribution coefficient of 0.6 and the temporal contribution coefficient of 0.4, the anisotropic spatiotemporal distance obtained by linear combination is 0.038. This distance represents the spatiotemporal correlation between the events.

[0036] In one possible embodiment, the step of performing data clustering on all the internal acoustic emission events based on the anisotropic spatiotemporal distance to form multiple independent acoustic emission event clusters includes: An acoustic emission event association graph is constructed using each of the aforementioned internal acoustic emission events as nodes and the anisotropic spatiotemporal distance as the edge weight. A density-based spatial clustering algorithm searches for core nodes in the acoustic emission event association graph and expands the clusters of adjacent nodes whose anisotropic spatiotemporal distance is less than a preset clustering radius, outputting the multiple independent acoustic emission event clusters.

[0037] The number of nodes in the acoustic emission event correlation graph corresponds to the total number of acoustic emission events within a selected time period, for example, including 5000 microfracture events. The edge weight of the connection between any two nodes is assigned to the anisotropic spatiotemporal distance calculated above. The preset clustering radius is determined by calculating the inflection point of the K-distance distribution curve, preferably ranging from 0.02 to 0.08, for example, 0.035; the minimum number of nodes is preferably set to a range of 5 to 20, for example, 10. During the clustering expansion process, the correlation graph is traversed. When the number of neighboring nodes of an event node within a spatiotemporal distance radius of 0.035 is greater than or equal to 10, the node is promoted to a core node. All neighboring nodes of the node that meet the radius condition are included in the current cluster, and the neighborhood of the neighboring nodes is recursively searched to achieve continuous expansion of the cluster boundary until no new core nodes are added. Isolated nodes without neighbors or with insufficient numbers are set as noise points and removed, outputting multiple independent acoustic emission event clusters representing different fracture development concentration areas, for example, 3 to 5.

[0038] S3, spatially pair the region with the highly concentrated strain with the multiple independent acoustic emission event clusters, calculate the spatial morphological conformity of each pair to be matched, and calculate the temporal evolution synchronization of each pair to be matched. Use the temporal evolution synchronization as a modulation coefficient to perform data augmentation calculation on the spatial morphological conformity to obtain the matching confidence. When the matching confidence is greater than the preset discrimination threshold, determine that the pair to be matched corresponds to a real microscopic fracture in the rock mass.

[0039] An orthogonal projection algorithm is used to project the three-dimensional coordinates of each independent acoustic emission event cluster onto the two-dimensional plane containing the surface image sequence, forming a two-dimensional acoustic emission projection point cloud. Based on the coordinate transformation relationship between the acoustic emission coordinate system and the image physical coordinate system pre-calibrated before loading, the two-dimensional acoustic emission projection point cloud is transformed to the same surface physical coordinate system as the strain-concentrated region. In the spatial morphology matching calculation, fitting and registration that would change the relative spatial position relationship between the point cloud to be matched and the strain-concentrated region is not performed. The centroid of the acoustic emission event cluster projection and the direction of the first principal axis are as follows. Figure 3As shown. The projection centroid and the first principal axis projection of the two-dimensional acoustic emission projection point cloud are extracted, and the geometric center and the first principal axis direction of the strain-concentrated region are extracted; the Euclidean distance between the geometric center and the projection centroid, and the spatial angle between the first principal axis direction of the strain-concentrated region and the projection of the first principal axis direction are calculated; after dividing the Euclidean distance and the spatial angle by a preset reference distance and a preset reference angle respectively for dimensionless processing, the results are input into a negative exponential decay evaluation function to obtain the spatial morphology conformity.

[0040] In the calculation of time evolution synchronization, the effective strain rate of change over time within the strain-concentrated region and the cumulative acoustic emission energy rate of change over time within the corresponding independent acoustic emission event clusters are extracted. The two rate of change sequences are resampled and aligned within a sliding time window and transformed to the frequency domain. The cross-power spectrum is calculated, and the time-delay cross-correlation function is obtained through inverse Fourier transform. The normalized maximum positive correlation peak is extracted and amplitude-limited to obtain the window synchronization of the current sliding time window. When multiple sliding time windows exist, the maximum value of the synchronization of each effective window is used as the time evolution synchronization of the pair to be matched. Subsequently, the time evolution synchronization is multiplied by the spatial morphology conformity as a modulation coefficient to obtain the preliminary confidence score. When both the spatial morphology conformity and the time evolution synchronization score meet the corresponding minimum threshold conditions, the preliminary confidence score is weighted and corrected using the normalized maximum principal strain value, and the correction result is amplitude-limited. When either index does not meet the corresponding minimum threshold condition, the preliminary confidence score is directly used as the matching confidence score, and the output is a matching confidence score between 0 and 1. The matching confidence level is compared with the discrimination threshold set to 0.70. When the value is greater than the discrimination threshold, it is determined that the current strain-concentrated area and the corresponding acoustic emission event cluster are of the same origin in terms of physical mechanism, and the matching pair is confirmed to correspond to a real micro-fracture in the rock mass.

[0041] In one possible embodiment, the step of spatially pairing the region with the highly concentrated strain with the plurality of independent acoustic emission event clusters, and calculating the spatial morphological conformity of each pair to be matched, includes: The geometric center and the first principal axis direction of the region with highly concentrated strain are extracted respectively, as well as the projection centroid of the event location point cloud on the rock surface and the projection of the first principal axis direction obtained based on principal component analysis in the multiple independent acoustic emission event clusters. Calculate the Euclidean distance between the geometric center and the projected centroid, and the spatial angle between the first principal axis direction of the region with high strain concentration and the projection of the first principal axis direction; After dividing the Euclidean distance and the spatial angle by preset reference distance and reference angle respectively to perform dimensionless processing, the results are input into a preset negative exponential decay evaluation function for calculation, to obtain the spatial morphology conformity with a value between 0 and 1.

[0042] During spatial feature extraction, a second-order image moment algorithm is used to calculate the geometric center coordinates and orientation angle of the two-dimensional strain-concentrated region. Simultaneously, the three-dimensional acoustic emission event cluster point cloud is projected perpendicularly along the Z' axis onto the surface, and the two-dimensional centroid coordinates of the projected point cloud are obtained. The maximum eigenvector of the covariance matrix is ​​then calculated to extract the projection angle along the first principal axis. The preset reference distance is typically set to 5% to 10% of the specimen size, for example, 10 mm; the preset reference angle is set to 10° to 30°, for example, 15°, to reduce the dimensional differences between size and angle. In the calculation of the negative exponential decay evaluation function, the spatial morphology conformity can be expressed as… , where d is the Euclidean distance between the geometric center and the projected centroid. The reference distance is set as a preset value, and θ is the spatial angle between the two primary axis directions. To preset the reference angle, and These are the weighting coefficients, and + =1. Since the first principal axis direction is equivalent in both directions, the spatial angle is determined by the inverse cosine of the absolute value of the dot product of the two unit direction vectors, i.e. Where u and v are the unit direction vectors of the two first principal axes, such that θ falls within the range of 0° to 90°. For example, if the calculated Euclidean distance between the geometric center and the projected centroid of a certain pair to be matched is 2mm, the spatial angle between the two principal axes is 3°, the reference distance is 10mm, and the reference angle is 15°, then... and Taking a value of 0.5 for all values, the spatial morphology conformity calculated using the negative exponent formula is approximately 0.819. This numerical result indicates the degree of agreement between the internal micro-fracture zone and the surface strain concentration zone in terms of macroscopic spatial location and extension direction.

[0043] In one possible embodiment, calculating the temporal evolution synchronization degree of each pair to be matched includes: Extract the effective strain rate of change over time sequence within the region of high strain concentration, and the cumulative acoustic emission energy rate of change over time sequence within the multiple independent acoustic emission event clusters; A sliding time window is set, and the two rate of change sequences are resampled and aligned within the sliding time window. Then, a fast Fourier transform is performed to convert them to the frequency domain. Calculate the cross power spectrum of the two rate-of-change sequences in the frequency domain, and obtain the time-delay cross-correlation function through inverse Fourier transform; The maximum positive correlation peak of the time-delay cross-correlation function is extracted within the preset allowable time delay range. When both rate of change sequences meet the preset effective fluctuation conditions, the maximum positive correlation peak is normalized and subjected to amplitude limiting processing as the time evolution synchronization degree. When any rate of change sequence does not meet the preset effective fluctuation condition, the time evolution synchronization degree of the corresponding window is set to 0 or the corresponding window is skipped.

[0044] The temporal feature extraction process divides the entire loading process into fixed time steps, such as 0.1 seconds. The area-weighted average of the absolute values ​​of the maximum principal strains at each calculation node within the strain-concentrated region is defined as the effective strain, and the effective strain rate of change sequence is calculated based on the difference between adjacent time steps. Simultaneously, the event energy within the acoustic emission event cluster is accumulated over time, and the accumulated acoustic emission energy rate of change sequence is calculated based on the difference between adjacent time steps. To address the inconsistency between the surface optical acquisition and internal acoustic acquisition frequencies, a sliding time window of 20 seconds is set. Within this window, a cubic spline interpolation algorithm is used to resample and align the two sequences to a standard sampling frequency of 100Hz. A Hamming window is applied to the aligned sequences to reduce spectral leakage, and a fast Fourier transform is performed to enter the frequency domain. In the frequency domain, the cross-power spectrum is obtained by multiplying the strain frequency domain signal and the energy frequency domain signal by their conjugates. An inverse Fourier transform is then performed to output a time-delay cross-correlation function containing time lead and lag information. Considering the reasonable delay of the physical response, a preset allowable time delay range is set to -2 seconds to +2 seconds. Search for the extreme value of the function within this range and extract the peak value of the maximum positive correlation. If the extracted peak value of the maximum positive correlation is 0.75× The maximum positive correlation peak is normalized by dividing it by the geometric mean of the peak values ​​at zero lag of the autocorrelation function of each of the two sequences. If the normalization result is less than 0, the time evolution synchronization is set to 0; if the normalization result is greater than 1, it is limited to 1. The peak value at zero lag of the autocorrelation function refers to the peak value obtained by correlation calculation between the rate of change sequence and itself in a time-off state, used to characterize the energy of the sequence itself and as a benchmark for cross-correlation normalization. When the variance of any rate of change sequence is less than a preset minimum threshold ε, or the peak value at zero lag of the autocorrelation function corresponding to any rate of change sequence is less than the preset minimum threshold ε, it is determined that there is no effective synchronization comparison condition within the current sliding time window, and the time evolution synchronization of that window is set to 0 or the window is skipped. When there are multiple sliding time windows, the maximum value of the synchronization of each effective window is taken as the time evolution synchronization of the pair to be matched. The dimensionless value of the pair to be matched, for example, 0.82, is the time evolution synchronization, which quantifies the transient synchronization characteristics of surface strain abrupt change and internal energy release.

[0045] In one possible embodiment, the step of using the temporal evolution synchronization degree as a modulation coefficient to perform data augmentation calculation on the spatial morphological conformity to obtain a matching confidence degree includes: The time evolution synchronization degree is used as a modulation coefficient with a value between 0 and 1; Multiplying the spatial morphology conformity by the modulation coefficient yields the preliminary confidence level; When both the spatial morphology conformity and the temporal evolution synchronization meet the corresponding minimum threshold conditions, the initial confidence level is weighted and corrected using the maximum principal strain value of the strain-concentrated region after normalization, and the correction result is limited so that the output matching confidence level is kept between 0 and 1. When the spatial morphology conformity or the temporal evolution synchronization does not meet the corresponding minimum threshold condition, the preliminary confidence level is used as the matching confidence level.

[0046] The calculated temporal evolution synchronization degree, such as 0.82, is used as the modulation coefficient and multiplied by the spatial morphological conformity degree, such as 0.818, to obtain a preliminary confidence level of 0.67. This multiplication mechanism constitutes a double check; under the condition of highly overlapping spatial geometry and synchronized fracture temporal evolution, the pair to be matched can obtain a high baseline confidence score, filtering out false positive results due to coincidental overlap. In the weighted correction stage, the absolute value of the maximum principal strain at the nodes in the highly concentrated strain region is extracted, and the absolute value is input into a preset Sigmoid nonlinear mapping function, such as setting the strain reference value as 0.818. The smoothing coefficient is set to 1000 and converted to a local intensity weighting factor between 0 and 1. When the spatial morphology conformity is not lower than a preset minimum spatial threshold and the temporal evolution synchronization is not lower than a preset minimum synchronization threshold, the formula is used... Perform augmented computation, in which The initial confidence level is defined as W, the local intensity weighting factor is defined as λ, and the correction coefficient is defined as min{1,·}, which represents the upper limit truncation of the enhancement results to ensure that the final matching confidence level does not exceed 1. When the spatial morphological conformity or temporal evolution synchronization is lower than the corresponding minimum threshold, the maximum principal strain enhancement correction is not performed, and the initial confidence level is used as the matching confidence level. The correction coefficient is defined as follows: The preferred value is 0.15. If the weighting factor W is 0.9, the output matching confidence score is corrected and enhanced to 0.76. When this confidence score is greater than the preset discrimination threshold of 0.70, it can be determined that the surface macroscopic strain and the internal acoustic emission clusters originate from the development of the same real rock mass microscopic fractures.

[0047] This invention also proposes a rock mass micro-fracture testing system, comprising the following modules: The calculation module is used to acquire surface image sequences and internal acoustic emission event datasets during the loading process of the rock mass. The internal acoustic emission event dataset contains the three-dimensional coordinates, occurrence time and energy value data of each internal acoustic emission event. The surface image sequence is processed based on digital image correlation to obtain the discrete strain field of the rock mass surface. The second invariant of the eccentric curvature corresponding to the discrete strain field is calculated, and the region with high strain concentration is identified by setting a threshold. The extraction module is used to extract the principal strain direction of the rock mass surface at the characteristic loading time as the reference direction, calculate the anisotropic spatiotemporal distance of each internal acoustic emission event in three-dimensional space and time interval, and perform data clustering operation on all internal acoustic emission events according to the anisotropic spatiotemporal distance to form multiple independent acoustic emission event clusters. The determination module is used to spatially pair the region with the highly concentrated strain with the multiple independent acoustic emission event clusters, calculate the spatial morphological conformity of each pair to be matched, calculate the temporal evolution synchronization of each pair to be matched, use the temporal evolution synchronization as a modulation coefficient to perform data augmentation calculation on the spatial morphological conformity to obtain the matching confidence. When the matching confidence is greater than a preset discrimination threshold, it is determined that the pair to be matched corresponds to a real microscopic fracture in the rock mass.

[0048] The above description of the disclosed embodiments enables those skilled in the art to make or use this application. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of this application. Therefore, this application is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for testing microscopic fractures in rock mass, characterized in that, Includes the following steps: A sequence of surface images and a dataset of internal acoustic emission events are obtained during the loading process of the rock mass. The dataset of internal acoustic emission events includes the three-dimensional coordinates, occurrence time and energy value data of each internal acoustic emission event. The surface image sequence is processed based on digital image correlation to obtain the discrete strain field of the rock mass surface. The second invariant of the skew curvature corresponding to the discrete strain field is calculated, and the region with high strain concentration is identified by setting a threshold. The principal strain direction of the rock mass surface at the characteristic loading moment is extracted as the reference direction. The anisotropic spatiotemporal distance of each internal acoustic emission event in three-dimensional space and time interval is calculated. Based on the anisotropic spatiotemporal distance, data clustering operation is performed on all internal acoustic emission events to form multiple independent acoustic emission event clusters. The region with highly concentrated strain is spatially paired with the multiple independent acoustic emission event clusters. The spatial morphological conformity of each pair to be matched is calculated, and the temporal evolution synchronization of each pair to be matched is calculated. The temporal evolution synchronization is used as a modulation coefficient to perform data augmentation calculation on the spatial morphological conformity to obtain the matching confidence. When the matching confidence is greater than the preset discrimination threshold, it is determined that the pair to be matched corresponds to a real microscopic fracture in the rock mass. The process of using the principal strain direction of the rock mass surface at the time of feature loading as the reference direction to calculate the anisotropic spatiotemporal distance of each internal acoustic emission event in three-dimensional space and time interval includes: constructing a three-dimensional anisotropic local coordinate system with the principal strain direction of the rock mass surface as the first coordinate axis, the direction perpendicular to the principal strain direction of the rock mass surface as the second coordinate axis, and the rock mass depth direction as the third coordinate axis; calculating the coordinate difference vector of any two internal acoustic emission events in three-dimensional space under the three-dimensional anisotropic local coordinate system, and assigning different distance penalty weights to the corresponding spatial components according to the preset principal deformation direction weight, surface lateral weight, and depth direction weight; calculating the time difference of occurrence of the two internal acoustic emission events; normalizing the magnitude of the weighted three-dimensional spatial coordinate difference vector in terms of spatial scale, normalizing the time difference of occurrence in terms of time scale, and then linearly combining them according to the preset spatial contribution coefficient and time contribution coefficient to obtain the anisotropic spatiotemporal distance; The calculation of the time evolution synchronization degree of each pair to be matched includes: extracting the effective strain rate of change sequence over time in the region where the strain is highly concentrated, and the cumulative acoustic emission energy rate of change sequence over time in the multiple independent acoustic emission event clusters; setting a sliding time window, resampling and aligning the two rate of change sequences within the sliding time window, and then performing a fast Fourier transform to convert to the frequency domain; calculating the cross power spectrum of the two rate of change sequences in the frequency domain, and obtaining the time-delay cross-correlation function through inverse Fourier transform; extracting the maximum positive correlation peak of the time-delay cross-correlation function within a preset allowable time delay range; when both rate of change sequences meet the preset effective fluctuation condition, normalizing and limiting the maximum positive correlation peak as the time evolution synchronization degree; when either rate of change sequence does not meet the preset effective fluctuation condition, setting the time evolution synchronization degree of the corresponding window to 0 or skipping the corresponding window; The process of spatially pairing the strain-concentrated region with the multiple independent acoustic emission event clusters and calculating the spatial morphological conformity of each pair includes: extracting the geometric center and first principal axis direction of the strain-concentrated region, and the projection centroid of the event location point cloud on the rock surface and the projection of the first principal axis direction obtained based on principal component analysis; calculating the Euclidean distance between the geometric center and the projection centroid, and the spatial angle between the first principal axis direction of the strain-concentrated region and the projection of the first principal axis direction; dividing the Euclidean distance and the spatial angle by a preset reference distance and reference angle respectively, performing dimensionless processing, and inputting them into a preset negative exponential decay evaluation function for calculation to obtain the spatial morphological conformity with values ​​between 0 and 1; in the calculation of the negative exponential decay evaluation function, the spatial morphological conformity can be expressed as... , where d is the Euclidean distance between the geometric center and the projected centroid. The reference distance is set as a preset value, and θ is the spatial angle between the two primary axis directions. To preset the reference angle, and These are the weighting coefficients, and + =1.

2. The method according to claim 1, characterized in that, The step of processing the surface image sequence based on digital image correlation to obtain the discrete strain field of the rock mass surface, and calculating the second invariant of the skew curvature corresponding to the discrete strain field, includes: The first image of the rock mass surface before loading is set as the reference image, and the subsequent loaded images are matched by digital image correlation to obtain the discrete displacement field. The discrete displacement field is smoothed and differentiated using a local surface fitting algorithm to obtain a continuous strain tensor field. The second-order partial derivatives of each independent strain component of the strain tensor field with respect to the surface spatial coordinates are calculated to obtain the second-order curvature matrix corresponding to each independent strain component. The inner product of the partial curvature parts of each second-order curvature matrix is ​​summed to obtain the second invariant of partial curvature, which represents the degree of non-uniformity of local deformation.

3. The method according to claim 2, characterized in that, The method of identifying regions with highly concentrated strain by setting a threshold includes: The distribution characteristics of the second invariant of skewness curvature are statistically analyzed across the entire field, and the mean and standard deviation of the invariant are calculated. The sum of the mean and the standard deviation multiplied by a preset factor is set as the judgment threshold; Extract the connected components where the second invariant of the bias curvature is greater than the determination threshold, and mark the connected components as the regions where the strain is highly concentrated.

4. The method according to claim 1, characterized in that, The step of performing data clustering operations on all the internal acoustic emission events based on the anisotropic spatiotemporal distance to form multiple independent acoustic emission event clusters includes: An acoustic emission event association graph is constructed using each of the aforementioned internal acoustic emission events as nodes and the anisotropic spatiotemporal distance as the edge weight. A density-based spatial clustering algorithm searches for core nodes in the acoustic emission event association graph and expands the clusters of adjacent nodes whose anisotropic spatiotemporal distance is less than a preset clustering radius, outputting the multiple independent acoustic emission event clusters.

5. The method according to claim 1, characterized in that, The step of using the temporal evolution synchronization degree as a modulation coefficient to perform data augmentation calculation on the spatial morphological conformity to obtain the matching confidence degree includes: The time evolution synchronization degree is used as a modulation coefficient with a value between 0 and 1; Multiplying the spatial morphology conformity by the modulation coefficient yields the preliminary confidence level; When both the spatial morphology conformity and the temporal evolution synchronization meet the corresponding minimum threshold conditions, the initial confidence level is weighted and corrected using the maximum principal strain value of the strain-concentrated region after normalization, and the correction result is limited so that the output matching confidence level is kept between 0 and 1. When the spatial morphology conformity or the temporal evolution synchronization does not meet the corresponding minimum threshold condition, the preliminary confidence level is used as the matching confidence level.

6. A rock mass micro-fracture testing system according to any one of claims 1-5, characterized in that, Includes the following modules: The calculation module is used to acquire surface image sequences and internal acoustic emission event datasets during the loading process of the rock mass. The internal acoustic emission event dataset contains the three-dimensional coordinates, occurrence time and energy value data of each internal acoustic emission event. The surface image sequence is processed based on digital image correlation to obtain the discrete strain field of the rock mass surface. The second invariant of the eccentric curvature corresponding to the discrete strain field is calculated, and the region with high strain concentration is identified by setting a threshold. The extraction module is used to extract the principal strain direction of the rock mass surface at the characteristic loading time as the reference direction, calculate the anisotropic spatiotemporal distance of each internal acoustic emission event in three-dimensional space and time interval, and perform data clustering operation on all internal acoustic emission events according to the anisotropic spatiotemporal distance to form multiple independent acoustic emission event clusters. The determination module is used to spatially pair the region with the highly concentrated strain with the multiple independent acoustic emission event clusters, calculate the spatial morphological conformity of each pair to be matched, calculate the temporal evolution synchronization of each pair to be matched, use the temporal evolution synchronization as a modulation coefficient to perform data augmentation calculation on the spatial morphological conformity to obtain the matching confidence. When the matching confidence is greater than a preset discrimination threshold, it is determined that the pair to be matched corresponds to a real microscopic fracture in the rock mass.

7. The system according to claim 6, characterized in that, The step of processing the surface image sequence based on digital image correlation to obtain the discrete strain field of the rock mass surface, and calculating the second invariant of the skew curvature corresponding to the discrete strain field, includes: The first image of the rock mass surface before loading is set as the reference image, and the subsequent loaded images are matched by digital image correlation to obtain the discrete displacement field. The discrete displacement field is smoothed and differentiated using a local surface fitting algorithm to obtain a continuous strain tensor field. The second-order partial derivatives of each independent strain component of the strain tensor field with respect to the surface spatial coordinates are calculated to obtain the second-order curvature matrix corresponding to each independent strain component. The inner product of the partial curvature parts of each second-order curvature matrix is ​​summed to obtain the second invariant of partial curvature, which represents the degree of non-uniformity of local deformation.

Images