A method and system for analyzing seismic activity based on Matlab software

Abstract

The present application provides a method and system for analyzing earthquake activity based on Matlab software. Among them, the method deploys a broadband acoustic wave sensor array to collect original data of focal depth; based on Matlab, the data is divided into equidistant grids, and a depth and frequency histogram is constructed, with the horizontal axis being the equidistant depth grid and the vertical axis recording the original frequency of earthquake events in each grid and retaining nonlinear characteristics; the acoustic wave propagation time difference of the same earthquake event is calculated, and if the time difference inversion depth is within the preset tolerance range, it is marked as valid data; based on Matlab, the valid data and the original frequency of the histogram are subjected to unsupervised classification processing, spatial earthquake clusters are iteratively divided, and classification boundaries are generated; the classification boundaries are combined with geographic coordinates to output a classification boundary map, and the dominant areas of shallow and deep earthquakes are distinguished by color blocks. The present application improves the accuracy of automatic division of shallow crust focal depth distribution patterns and the efficiency of geographic visualization expression.

Description

Technical Field

[0001] The present application relates to the technical field of seismic activity analysis, and in particular to a seismic activity analysis method and system based on Matlab software. Background Art

[0002] In shallow crustal earthquake monitoring, the automatic classification of focal depth distribution patterns is crucial for earthquake early warning and geological tectonic research. Because shallow and deep earthquakes exhibit significant differences in activity, this requires real-time processing of massive amounts of seismic data and automatic identification of the spatial distribution of earthquake clusters at different depths. Furthermore, the classification results must be combined with geographic coordinates to provide intuitive decision-making for disaster assessment.

[0003] The current mainstream approach uses a machine learning-based earthquake depth classification model. This model extracts the time-frequency characteristics of earthquake waveforms, combines them with regional crustal velocity models to invert focal depth, and then uses a clustering algorithm to group depth data. This approach trains on historical earthquake data to build a classification model, outputs depth distribution clustering results, and overlays the classification boundaries in a geographic information system.

[0004] Existing solutions have the following flaws: they rely on large amounts of labeled data to train models and are unable to handle scenarios where there is no historical earthquake data in the monitoring area; the velocity model ignores the heterogeneity of shallow crustal media when inverting depth, resulting in accumulated inversion errors; the superposition of classification results and geographic coordinates relies on manual adjustment of layer parameters, making it difficult to achieve automated boundary mapping, and the nonlinear characteristics of the original depth distribution are not retained, affecting the credibility of the classification. Summary of the Invention

[0005] The present application provides a seismic activity analysis method and system based on Matlab software to solve the problem of insufficient accuracy in automatic division of earthquake source depth distribution patterns in the prior art.

[0006] In a first aspect, the present application provides a seismicity analysis method based on Matlab software, comprising:

[0007] The original data of earthquake focal depth are collected synchronously by using detection devices deployed in the shallow crust monitoring area;

[0008] Based on Matlab software, the original focal depth data is converted into an equidistant grid at a preset depth resolution, and a depth-frequency histogram is directly constructed using the equidistant grid. The horizontal axis of the depth-frequency histogram is the equidistant depth grid, and the vertical axis records the original frequency of earthquake events within each depth grid, while retaining the nonlinear characteristics of the depth distribution.

[0009] Calculating the acoustic wave propagation time difference of the same earthquake event between adjacent sensor nodes based on the positional relationship of adjacent sensor nodes in the detection device, and marking the depth data of the earthquake event as valid data if the inversion depth of the acoustic wave propagation time difference is within a preset tolerance range;

[0010] Performing unsupervised classification processing on the valid data and the original frequencies recorded in the depth and frequency histogram based on the Matlab software, wherein the processing process automatically divides the spatially continuous earthquake clusters by iterative calculation, and generates classification boundaries according to the sudden change gradient of the central depth of the earthquake cluster;

[0011] After combining the classification boundaries with the geographic coordinates, an earthquake classification boundary map is output through the Matlab software.

[0012] Optionally, the Matlab software is used to perform unsupervised classification processing on the valid data and the original frequencies recorded in the depth and frequency histogram, wherein the processing process iteratively calculates and automatically divides spatially continuous earthquake clusters, and generates classification boundaries according to the sudden change gradient of the central depth of the earthquake cluster, including:

[0013] Merging the depth data in the valid data with the original frequencies recorded in the depth and frequency histogram according to corresponding equidistant grids to generate a joint input data set including the depth data and the number of event occurrences;

[0014] performing unlabeled grouping processing on the joint input data set based on the Matlab software, wherein the grouping processing process minimizes the distribution difference of depth data and event occurrence times within the same group by adjusting the group boundaries, and forms spatially continuous earthquake clusters after the grouping processing;

[0015] Performing geometric center calculation on all depth data points in the earthquake cluster to obtain a central depth value, wherein the geometric center calculation process is completed only based on the numerical distribution of the depth data points;

[0016] The central depth values of adjacent earthquake clusters are compared in sequence along the depth direction. If the absolute difference between the adjacent central depth values exceeds the mutation threshold, the boundary between the adjacent earthquake clusters is used as the classification boundary.

[0017] Optionally, performing unlabeled grouping processing on the joint input data set based on the Matlab software, wherein the grouping processing minimizes the distribution difference of depth data and event occurrence times within the same group by adjusting group boundaries, and forming spatially continuous earthquake clusters after grouping processing, includes:

[0018] Dividing the depth data intervals and event frequency intervals in the joint input data set into discrete candidate groups according to initial group boundaries based on the Matlab software;

[0019] Traversing the boundary positions of the candidate groups, calculating the weighted difference between the depth data distribution dispersion and the event frequency distribution dispersion within the same group;

[0020] Iteratively adjusting the positions of the group boundaries so that adjacent groups with the smallest weighted difference values are automatically merged to generate optimized spatially continuous groups;

[0021] The groups with continuous depth data and uniform event frequency distribution in the spatially continuous groups are marked as spatially continuous earthquake clusters.

[0022] Optionally, the step of calculating the acoustic wave propagation time difference of the same earthquake event between adjacent sensor nodes based on the positional relationship of adjacent sensor nodes in the detection device, and marking the depth data of the earthquake event as valid data if the inversion depth of the acoustic wave propagation time difference is within a preset tolerance range, includes:

[0023] determining a relative distance between adjacent sensor nodes based on the layout position of the detection device, and generating a sensor node group consisting of at least three adjacent sensor nodes according to the relative distance;

[0024] For the same earthquake event, the arrival time of the acoustic wave signal recorded by each sensor node in the sensor node group is extracted, and the acoustic wave propagation time difference between adjacent sensor nodes is calculated based on the arrival time;

[0025] According to the correspondence between the acoustic wave propagation time difference and the relative distance, converting and calculating the correspondence to obtain an inversion depth of the earthquake event in the sensor node group; if the difference between the inversion depth and a preset tolerance range is a non-positive number, determining that the inversion depth of the sensor node group is valid;

[0026] The valid number of the inversion depths in all sensor node groups is counted, and if the ratio of the valid number to the total number of sensor node groups exceeds a set threshold, the depth data of the earthquake event is marked as valid data.

[0027] Optionally, the combining of the classification boundaries with the geographic coordinates and outputting an earthquake classification boundary map through the Matlab software comprises:

[0028] Obtaining geographic coordinate information of the shallow crust monitoring area, wherein the geographic coordinate information includes latitude and longitude data of the boundary points of the location range of the area corresponding to each earthquake cluster;

[0029] Binding each boundary point in the classification boundary to the latitude and longitude data of the corresponding position in the geographic coordinate information to generate a classification boundary coordinate point set containing geographic coordinates;

[0030] Assigning a unique color identifier to the earthquake cluster based on the difference in the center depth value of the earthquake cluster, and defining the color identifier assignment rule as follows: if the center depth value of the earthquake cluster is less than a preset shallow source threshold, a first color identifier is assigned; if the center depth value of the earthquake cluster is greater than a preset deep source threshold, a second color identifier is assigned;

[0031] Based on the Matlab software, the classification boundary coordinate point set is regionally filled according to the geographic spatial position. During the filling process, the first color identifier and the second color identifier are respectively covered to the geographic coordinate range of the corresponding earthquake cluster to generate an earthquake classification boundary map containing color blocks.

[0032] Optionally, generating an equidistant grid from the original focal depth data according to a preset depth resolution based on Matlab software, and directly constructing a depth and frequency histogram through the equidistant grid includes:

[0033] Determining the span of an equidistant grid according to a preset depth resolution, and generating an equidistant grid sequence covering all the original data of the focal depth according to the span, wherein the span is calculated by the ratio between the maximum detection depth of the monitoring area and the preset depth resolution;

[0034] Traversing the depth data points in the original data of the focal depth based on Matlab software, mapping the depth data points to the equidistant grid sequence, wherein the mapping process is completed by determining the grid span interval to which the depth data points belong;

[0035] Counting the number of depth data points mapped in the equidistant grid sequence, using the statistical results as the number of event occurrences, and generating an initial histogram structure with the equidistant grid as the horizontal axis and the number of event occurrences as the vertical axis;

[0036] The initial histogram structure is subjected to nonlinear feature retention processing based on Matlab software. The interpolation or smoothing operation on the number of event occurrences of adjacent equidistant grids is prohibited during the processing process, and the original statistical results are used as the vertical axis values of the processed histogram to form a depth and frequency histogram.

[0037] Optionally, converting and calculating the corresponding relationship between the acoustic wave propagation time difference and the relative distance to obtain the inversion depth of the seismic event in the sensor node group includes:

[0038] Obtaining the acoustic wave propagation time difference and the corresponding relative distance of adjacent nodes in the sensor node group, and generating a corresponding relationship set of time difference and distance;

[0039] Establishing a linear proportional relationship model between the acoustic wave propagation time difference and the relative distance according to each data point in the corresponding relationship set;

[0040] Based on the linear proportional relationship model, the acoustic wave propagation time difference is proportionally converted into a target depth, and during the conversion process, the product of the relative distance and the acoustic wave propagation time difference is kept as a fixed conversion coefficient;

[0041] An arithmetic mean is taken for the target depths converted from the seismic event in the sensor node group, and the inversion depth of the seismic event is output.

[0042] In a second aspect, the present application provides a seismic activity analysis system based on Matlab software, comprising:

[0043] An acquisition module is used to synchronously acquire raw data of earthquake focal depth through detection devices deployed in the shallow crust monitoring area;

[0044] A construction module is used to generate an equidistant grid according to a preset depth resolution from the original focal depth data based on Matlab software, and directly construct a depth and frequency histogram using the equidistant grid. The horizontal axis of the depth and frequency histogram is the equidistant depth grid, and the vertical axis records the original frequency of earthquake events in each depth grid, while retaining the nonlinear characteristics of the depth distribution;

[0045] a calculation module, configured to calculate the acoustic wave propagation time difference of the same earthquake event between adjacent sensor nodes according to the positional relationship of adjacent sensor nodes in the detection device, and mark the depth data of the earthquake event as valid data if the inversion depth of the acoustic wave propagation time difference is within a preset tolerance range;

[0046] a processing module for performing unsupervised classification processing on the valid data and the original frequencies recorded in the depth and frequency histogram based on the Matlab software, wherein the processing process automatically divides the spatially continuous earthquake clusters by iterative calculation and generates classification boundaries according to the sudden change gradient of the central depth of the earthquake cluster;

[0047] The output module is used to combine the classification boundary with the geographic coordinates and output the earthquake classification boundary map through the Matlab software.

[0048] In a third aspect, an embodiment of the present application provides a computing device comprising a processing component and a storage component; the storage component stores one or more computer instructions; the one or more computer instructions are used to be called and executed by the processing component to implement a seismic activity analysis method based on Matlab software as described in the first aspect above.

[0049] In a fourth aspect, an embodiment of the present application provides a computer storage medium storing a computer program. When the computer program is executed by a computer, it implements a seismic activity analysis method based on Matlab software as described in the first aspect.

[0050] In an embodiment of the present application, raw data of earthquake focal depth are synchronously collected by detection devices deployed in a shallow crust monitoring area; the raw data of earthquake focal depth are generated into an equidistant grid according to a preset depth resolution based on Matlab software, and a depth and frequency histogram is directly constructed through the equidistant grid. The horizontal axis of the depth and frequency histogram is the equidistant depth grid, and the vertical axis records the original frequency of earthquake events in each depth grid, and retains the nonlinear characteristics of the depth distribution; based on the positional relationship between adjacent sensor nodes in the detection device, the acoustic wave propagation time difference of the same earthquake event between adjacent sensor nodes is calculated, and if the inversion depth of the acoustic wave propagation time difference is within a preset tolerance range, the depth data of the earthquake event is marked as valid data; unsupervised classification processing is performed on the valid data and the raw frequency recorded in the depth and frequency histogram based on the Matlab software, and the processing process iteratively calculates and automatically divides spatially continuous earthquake clusters, and generates classification boundaries based on the sudden change gradient of the central depth of the earthquake cluster; after combining the classification boundaries with geographic coordinates, the earthquake classification boundary map is output through the Matlab software.

[0051] The technical solution of this application has the following beneficial effects:

[0052] Through the synchronous monitoring of the broadband acoustic wave sensor array, the temporal and spatial consistency of the original data of the earthquake source depth is ensured, and the coverage density and data reliability of the earthquake event detection in the shallow crust area are improved; based on the preset depth resolution, a histogram is directly constructed to retain the nonlinear characteristics of the original frequency distribution, avoid the distribution distortion caused by traditional data smoothing, and provide a real statistical benchmark for subsequent classification; through the inversion depth cross-validation of the multi-node acoustic wave propagation time difference, reliable depth data that conforms to the physical propagation laws is screened out, and the interference of abnormal data on the classification results is eliminated; the effective depth data and the histogram statistical features are integrated, and the spatial earthquake clusters are automatically divided through iterative optimization. The classification boundary is generated using the central depth mutation gradient to realize the depth distribution pattern recognition without human intervention; the classification boundary is dynamically bound to the geographic coordinates, and the dominant areas of shallow and deep earthquakes are intuitively displayed through color blocks, thereby improving the spatial interpretability of the earthquake source depth distribution pattern.

[0053] Furthermore, the effective depth data and the original frequency histogram records were merged into a joint dataset using an equidistant grid. Unlabeled grouping was then performed on the dataset using Matlab. Dynamic adjustment of group boundaries minimized differences in data distribution within groups, forming spatial earthquake clusters. The geometric center values of the depth data points in each cluster were calculated, and the absolute differences between adjacent center values were compared along the depth direction. If the absolute differences exceeded a mutation threshold, a classification boundary was generated. By combining the depth data with the statistical frequency and generating mutation gradient boundaries, the authors achieved automated classification of shallow crustal earthquake source depth distribution patterns, resolving the boundary ambiguity caused by traditional methods that rely on manually defined classification rules. This improved the physical consistency of depth clustering and the accuracy of geographic mapping.

[0054] These and other aspects of the present application will become more readily apparent from the description of the following embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0055] In order to more clearly illustrate the embodiments of the present application or the technical solutions in the prior art, a brief introduction will be given below to the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are some embodiments of the present application. For ordinary technicians in this field, other drawings can be obtained based on these drawings without any creative work.

[0056] Figure 1 A flowchart of a seismicity analysis method based on Matlab software provided by the present application is shown;

[0057] Figure 2 A scene diagram showing a seismic activity analysis method based on Matlab software provided by the present application is shown;

[0058] Figure 3 The present invention provides a schematic diagram of the structure of a seismic activity analysis system based on Matlab software;

[0059] Figure 4 A schematic structural diagram of a computing device provided by the present application is shown. DETAILED DESCRIPTION

[0060] In order to enable those skilled in the art to better understand the solution of the present application, the technical solution in the embodiments of the present application will be clearly and completely described below in conjunction with the drawings in the embodiments of the present application.

[0061] In some of the processes described in the specification and claims of this application and the above-mentioned figures, multiple operations that appear in a specific order are included, but it should be clearly understood that these operations may not be executed in the order in which they appear in this document or may be executed in parallel. The serial numbers of the operations, such as 101, 102, etc., are only used to distinguish between different operations, and the serial numbers themselves do not represent any order of execution. In addition, these processes may include more or fewer operations, and these operations may be executed in sequence or in parallel. It should be noted that the descriptions of "first", "second", etc. in this document are used to distinguish different messages, devices, modules, etc., and do not represent a sequential order, nor do they limit "first" and "second" to being different types.

[0062] Among existing earthquake source depth classification technologies, solutions based on machine learning models face three key flaws: first, model training is highly dependent on well-labeled historical earthquake data and cannot adapt to new monitoring areas without historical data accumulation; second, when using the regional average wave velocity model for depth inversion, the non-uniform distribution characteristics of the density and velocity of the shallow crust medium are ignored, resulting in the accumulation of inversion errors as the complexity of the propagation path increases; third, the fusion of classification results with geographic information relies on manual adjustment of layer parameters, with a low degree of automation. In addition, the nonlinear characteristics of the original depth distribution are smoothed during the classification process, which weakens the geological significance of the depth cluster boundaries.

[0063] In response to the above problems, this application proposes a seismic activity classification method based on the collaborative verification of acoustic wave sensor arrays and the joint analysis of deep statistical characteristics. By deploying a broadband acoustic wave sensor array to synchronously collect source depth data, a depth-frequency histogram that retains nonlinear characteristics is constructed, and the reliability of the data is cross-validated using the multi-node acoustic wave propagation time difference; based on the unsupervised classification processing of the verification data and the histogram statistical characteristics, the spatial earthquake clusters are automatically divided and the sudden gradient classification boundaries are generated, and finally the automatic mapping output of the classification boundaries and geographic coordinates is realized. This method breaks through the dependence on historical data through the unsupervised collaborative calculation of the verification of the physical propagation characteristics of acoustic waves and the statistical distribution characteristics; through the multi-node cross-inversion depth verification mechanism, the error accumulation caused by inhomogeneous media is significantly reduced; through the retention of the nonlinear characteristics of the histogram and the automated geographic mapping, the credibility and visualization efficiency of the classification boundaries are improved, providing high-precision and highly adaptable technical support for the analysis of shallow crustal seismic activity.

[0064] The technical solution of the present application can be applied to the automatic classification scenario of earthquake source depth distribution pattern.

[0065] The following will be combined with the drawings in the embodiments of this application to clearly and completely describe the technical solutions in the embodiments of this application. Obviously, the embodiments described are only part of the embodiments of this application, not all of the embodiments. Based on the embodiments in this application, all other embodiments obtained by those skilled in the art without making creative efforts are within the scope of protection of this application.

[0066] Figure 1 A flowchart of a seismic activity analysis method based on Matlab software is provided for the embodiment of the present application. Figure 1 As shown, the method includes:

[0067] 101. Synchronously collect raw data on earthquake focal depth using detection devices deployed in the shallow crust monitoring area;

[0068] In this scheme, the detection device refers to a broadband acoustic sensor array deployed in the shallow crust monitoring area. It comprises multiple sensor nodes, each capable of collecting broadband acoustic signals. The raw focal depth data is the time series information of acoustic wave propagation synchronously recorded by the sensor array, including the propagation time from the earthquake source to each sensor node, waveform amplitude, and timestamp.

[0069] In an embodiment of the present application, first, a broadband acoustic wave sensor array is arranged in a triangular grid topology in the shallow crust monitoring area. Each sensor node has a built-in high-precision clock module, and a timestamp synchronization algorithm is used to calibrate the clock of each node to ensure global time consistency.

[0070] Secondly, when an earthquake event occurs, each sensor node captures the acoustic wave signal in real time through an event trigger mechanism. The trigger condition is that the acoustic wave amplitude exceeds the preset threshold. Then, bandpass filtering technology is used to filter out environmental noise outside the frequency band of 20Hz to 200Hz, and the acoustic wave arrival time, waveform amplitude and timestamp are synchronously recorded to form the original data set.

[0071] Again, the original data set is aligned by timestamp. The alignment process matches the signal segments of the same earthquake event at different nodes based on the dynamic time warping algorithm, and eliminates abnormal data with time differences exceeding the preset tolerance, generating a time series set centered on the earthquake event and storing it as structured data.

[0072] Finally, for the same earthquake event, the time difference of sound wave arrival in the structured data is extracted and combined with the relative distance between nodes and the sound wave propagation velocity model V=√(E / ρ), where √ is the square root, E is the elastic modulus of the medium, and ρ is the density; the focal depth H=ΔT×V×cosθ is calculated through geometric relationships, where ΔT is the time difference of sound wave arrival and θ is the angle between the sound wave propagation path and the vertical direction. The original data of the focal depth containing time series information is output to provide input for subsequent verification and classification.

[0073] 102. Using Matlab software, generate an equidistant grid from the original focal depth data at a preset depth resolution, and directly construct a depth-frequency histogram using the equidistant grid. The horizontal axis of the depth-frequency histogram represents the equidistant depth grids, and the vertical axis records the original frequency of earthquake events within each depth grid, while retaining the nonlinear characteristics of the depth distribution.

[0074] Optionally, step 102 may specifically include the following steps:

[0075] 1021. Determine a span of an equidistant grid according to a preset depth resolution, and generate an equidistant grid sequence covering all original data of focal depths according to the span, wherein the span is calculated by a ratio between a maximum detection depth of the monitoring area and the preset depth resolution;

[0076] 1022. Traversing the depth data points in the original focal depth data based on Matlab software, and mapping the depth data points to the equidistant grid sequence, wherein the mapping process is completed by determining the grid span interval to which the depth data point belongs;

[0077] 1023. Counting the number of depth data points mapped in the equidistant grid sequence, using the statistical result as the number of event occurrences, and generating an initial histogram structure with the equidistant grid as the horizontal axis and the number of event occurrences as the vertical axis;

[0078] 1024. Perform nonlinear feature retention processing on the initial histogram structure based on Matlab software, prohibit interpolation or smoothing operations on the number of event occurrences of adjacent equidistant grids through the processing process, and use the original statistical results as the vertical axis value of the processed histogram to form a depth and frequency histogram.

[0079] In the above scheme, the depth and frequency histogram is a statistical chart with the horizontal axis being the equidistant depth grid and the vertical axis being the original number of earthquake events within the corresponding grid. Non-linear feature retention prohibits interpolation or smoothing of the histogram data to maintain the irregular fluctuation characteristics of the original frequency distribution. The equidistant grid sequence is a set of depth intervals evenly divided by span, covering the entire focal depth data range. The depth data point is the focal depth value corresponding to a single earthquake event. The initial histogram structure is a two-dimensional chart with the grid as the horizontal axis and the number of events as the vertical axis. Non-linear feature retention processing prohibits interpolation or smoothing of the number of events in adjacent grids, and directly uses the original statistical value to fill the vertical axis.

[0080] In the embodiment of the present application, first, step 1021 is used to calculate the span of a single grid according to the preset depth resolution defined by the user. The formula is span = maximum detection depth / preset depth resolution. If the calculated result is a decimal, it is rounded up to an integer. Based on the span value, an equidistant grid sequence is generated. The sequence starts at 0 and increases in span until the maximum detection depth is covered. For example, when the maximum detection depth is 10 kilometers and the resolution is 100 meters, span = 10000 / 100 = 100, generating 100 grids, respectively, 0-100 meters, 100-200 meters...9900-10000 meters.

[0081] Next, step 1022 iterates through all depth data points in the original focal depth data. For each data point, the following operations are performed: the depth value is divided by the preset depth resolution and the result is rounded down to obtain the grid number, i.e., grid number = floor(depth value / resolution). The grid number is used to determine the equidistant grid interval to which it belongs: starting depth = grid number × resolution, ending depth = (grid number + 1) × resolution. For example, a depth value of 5230 meters and a resolution of 100 meters would result in a grid number = floor(5230 / 100) = 52, corresponding to grid number 52 (5200-5300 meters). This process uses a Matlab loop structure to iterate through all data points and complete the mapping labeling.

[0082] Next, in step 1023, the number of depth data points within each equidistant grid is cumulatively counted. The statistical method is to traverse the grid sequence and count the data points corresponding to each grid number. The statistical results are used as the vertical axis value to generate an initial histogram with the grid number as the horizontal axis and the number of event occurrences as the vertical axis. For example, if grid number 52 contains 85 data points, the vertical axis value is 85. This process is implemented using the Matlab histcounts function, which takes as input the depth data point set and the equidistant grid boundary sequence and directly outputs the count results for each grid.

[0083] Finally, step 1024 prohibits interpolation or smoothing of the initial histogram's vertical axis values. Specifically, Matlab's smoothdata function and any filtering algorithms are disabled, and the original statistical values are used directly to fill the histogram's vertical axis. For example, the vertical axis value of grid cell 52 remains at 85, even if the adjacent grid cell 51 has a value of 10 and the adjacent grid cell 53 has a value of 200. No sliding average or curve fitting is performed. The resulting depth and frequency histogram retains the nonlinear fluctuation characteristics of the original distribution, such as local frequency mutations and multi-peak distribution.

[0084] In practical applications, assume that in a certain monitoring area, the preset depth resolution is set to 100 meters in step 1021, and the maximum detection depth is 10 kilometers. The grid span is calculated according to the formula: span = 10,000 meters / 100 meters = 100. This generates 100 equally spaced grids in the order 0-100 meters, 100-200 meters, and so on, to 9,900-10,000 meters. If the actual maximum depth during a given detection is 10,230 meters, the grid is expanded upward to 10,300 meters to ensure full data coverage.

[0085] In step 1022, each collected focal depth data point is mapped to its corresponding grid. For example, at 5230 meters, the grid number is calculated as floor(5230 / 100) = 52, corresponding to grid number 52 (5200-5300 meters). At a depth of 5350 meters, the grid number is calculated as floor(5350 / 100) = 53, corresponding to grid number 53 (5300-5400 meters). Matlab loops through all 50,000 data points, labeling all data points as grids.

[0086] In step 1023, the number of marked data points within each grid is counted: Grid 52 contains 85 data points, such as those at 5230 meters and 5255 meters, with the vertical axis value set to 85; Grid 9 (900-1000 meters) contains 300 data points, such as those at 1980 meters and 995 meters, with the vertical axis value set to 300. Matlab's histcounts function inputs the data point set and the grid boundary sequence, outputs the count results for each grid, and generates an initial histogram with the grid number on the horizontal axis and the number of events on the vertical axis.

[0087] In step 1024, Matlab's smoothdata function and any filtering algorithms are disabled, and the histogram is populated directly with the original statistical values: Grid 52 retains 85 occurrences, while the adjacent grid 51 has 10 occurrences and the adjacent grid 53 has 200 occurrences, forming a steep "trough, peak, trough" distribution. Grid 9 has 300 occurrences, which is significantly different from the adjacent grid 8 (50 occurrences) and the adjacent grid 10 (20 occurrences), creating a significant frequency drop. The final histogram retains all mutation features, clearly reflecting the uneven distribution pattern of high frequency in the shallow 0-2000 meters region, low frequency in the mid-5000-6000 meters region, and medium frequency in the deep 8000-9000 meters region.

[0088] The overall solution of step 102 above dynamically generates an equidistant grid by presetting the depth resolution and strictly retaining the nonlinear characteristics of the original frequency distribution, thereby solving the distribution distortion problem caused by data smoothing in traditional methods; through grid mapping and statistics, it provides high-fidelity input for source depth classification, supports the unsupervised classification algorithm to accurately identify the distribution boundaries of shallow-source and deep-source earthquake clusters, and avoids the intervention of manually setting classification rules, significantly improving the efficiency of automated processing and the credibility of geological interpretation.

[0089] 103. Calculate the acoustic wave propagation time difference of the same earthquake event between adjacent sensor nodes based on the positional relationship of adjacent sensor nodes in the detection device. If the inversion depth of the acoustic wave propagation time difference is within a preset tolerance range, mark the depth data of the earthquake event as valid data.

[0090] Optionally, step 103 may specifically include the following steps:

[0091] 1031. Determine a relative distance between adjacent sensor nodes based on the layout position of the detection device, and generate a sensor node group consisting of at least three adjacent sensor nodes according to the relative distance;

[0092] 1032. For the same earthquake event, extract the arrival time of the acoustic wave signal recorded by each sensor node in the sensor node group, and calculate the acoustic wave propagation time difference between adjacent sensor nodes based on the arrival time;

[0093] 1033. According to the correspondence between the acoustic wave propagation time difference and the relative distance, convert and calculate the correspondence to obtain an inversion depth of the earthquake event in the sensor node group. If the difference between the inversion depth and a preset tolerance range is a non-positive number, determine that the inversion depth of the sensor node group is valid.

[0094] Among them, step 1033 may specifically include the following processes: obtaining the sound wave propagation time difference and the corresponding relative distance of adjacent nodes in the sensor node group, and generating a corresponding relationship set between the time difference and the distance; establishing a linear proportional relationship model between the sound wave propagation time difference and the relative distance according to each data point in the corresponding relationship set; based on the linear proportional relationship model, converting the sound wave propagation time difference into a target depth in proportion, and keeping the product of the relative distance and the sound wave propagation time difference as a fixed conversion coefficient during the conversion process; taking the arithmetic mean of the target depths obtained by converting the seismic event in the sensor node group, and outputting the inversion depth of the seismic event.

[0095] 1034. Count the valid numbers of the inversion depths in all sensor node groups. If the ratio of the valid numbers to the total number of sensor node groups exceeds a set threshold, mark the depth data of the earthquake event as valid data.

[0096] In the above scheme, adjacent sensor nodes refer to sensors that are physically adjacent and capable of signal synchronization. The acoustic propagation time difference is the time difference between the acoustic signals of the same earthquake event reaching adjacent sensors. The inversion depth is an estimate of the focal depth calculated based on the time difference and the sensor spacing. The relative distance is the straight-line distance between adjacent sensor nodes. A sensor node group is a verification unit consisting of at least three adjacent sensors. The arrival time is the precise timestamp of the acoustic signal arriving at the sensor. The time difference calculation is the arrival time difference of the same earthquake event at adjacent nodes. The conversion calculation is the inference of the focal depth based on the physical relationship between time difference and distance. A non-positive difference means that the absolute value of the error between the inversion depth and the true depth does not exceed the tolerance range.

[0097] In this embodiment, a sensor node group is first constructed in step 1031. Based on the sensor array's layout coordinates (three-dimensional coordinates x1, y1, z1; x2, y2, z2, etc.), the straight-line distance between adjacent nodes is calculated using the formula: relative distance = ((x2 - x1)^2 + (y2 - y1)^2 + (z2 - z1)^2)^(1 / 2). At least three adjacent nodes are selected to form a verification group, and the distance between nodes within the group must be less than a preset maximum verification distance. For example, node A (0, 0, 0), node B (1200, 0, 0), and node C (600, 1000, 0) are calculated to have distances of 1200 meters, 1118 meters, and 1000 meters, respectively, forming a triangular verification group.

[0098] Next, in step 1032, the acoustic propagation time difference is calculated. For the same earthquake event, the acoustic wave arrival timestamps recorded by each node within the sensor node group are extracted, and the acoustic wave propagation time difference between adjacent pairs of nodes is calculated: the time difference between adjacent nodes A and B is ΔTAB = TA - TB. This process relies on high-precision clock synchronization (for example, GPS synchronization error ≤ 1 microsecond) to ensure the accuracy of the time difference calculation. For example, if node A has a timestamp of T1 = 10:00:00.000, node B has a timestamp of T2 = 10:00:00.350, and node C has a timestamp of T3 = 10:00:00.520, then ΔT_AC = 0.52 seconds and ΔT_BC = 0.17 seconds.

[0099] Next, the validity of the inversion depth is determined in step 1033. Based on the acoustic wave propagation velocity model: velocity V = (elastic modulus E / density ρ)^(1 / 2), combined with the time difference ΔT, the focal depth H is calculated as ΔT × V × cos (angle θ between the propagation path and the vertical direction). The absolute value of the difference between H and the true depth is calculated. If the absolute value of the difference is less than or equal to the preset tolerance, the inversion depth is considered valid. Assuming V = 5000 m / s, the distance between nodes AB D = 1200 meters, ΔT = 0.35 seconds, and θ = 30 degrees, then HAB = 0.35 × 5000 × cos(30°) = 1515 meters. If the true depth is 1600 meters and the preset tolerance is ±200 meters, |1515 - 1600| = 85 ≤ 200, indicating validity.

[0100] Finally, in step 1034, the effective ratio is calculated and the data is marked. The number of valid inversion depths in all sensor node groups is counted, and the effective ratio is calculated as: effective number / total number of node groups × 100%. If the effective ratio is less than the set threshold, the depth data for the earthquake event is marked as valid. For example, if the effective ratio is 70%, the data is retained if the threshold is set to 65%, otherwise it is discarded.

[0101] In practical applications, assuming that in a sensor array deployed in a certain volcanic area, step 1031 selects node A (0, 0, 0), node B (1500, 0, 0), and node C (750, 1299, 0), and calculates the spacing D_AB = 1500 meters, D_AC = 1500 meters, and D_BC = 1500 meters to form an equilateral triangle verification group.

[0102] In step 1032, node A records the arrival time as T1 = 08:30:15.000, node B as T2 = 08:30:15.300, and node C as T3 = 08:30:15.519, and calculates ΔT_AB = 0.3 seconds and ΔT_AC = 0.519 seconds.

[0103] In step 1033, assuming the medium velocity V = 5200 m / s, the inversion depth at node AB is H_AB = 0.3 × 5200 × cosine(25 degrees) = 1410 meters, and the inversion depth at node AC is H_AC = 0.519 × 5200 × cosine(40 degrees) = 2065 meters. If the true depth is 1500 meters and the tolerance is ±250 meters, then H_AB is valid (with an error of 90 meters) and H_AC is invalid (with an error of 565 meters).

[0104] In step 1034, the earthquake event verified 15 node groups, of which 11 had an inversion depth error of ≤ 250 meters, for a validity rate of 73.3%. A threshold of 70% was set, marking the data as valid and inputting it into the classification algorithm.

[0105] The overall solution of the above-mentioned step 103, through the inversion verification of the acoustic wave propagation time difference of multiple node groups and the dynamic judgment of the tolerance threshold, screens out reliable data that conforms to physical laws, and solves the error accumulation problem caused by medium heterogeneity in traditional methods; automatically marks data through the effective ratio threshold to avoid manual intervention, improves the data quality of the input classification algorithm, provides high-confidence input for the automatic classification of the source depth distribution pattern, and enhances the physical consistency of the classification boundary and the reliability of the geological interpretation.

[0106] 104. Performing unsupervised classification processing on the valid data and the original frequencies recorded in the depth and frequency histogram based on the Matlab software, wherein the processing process automatically divides the earthquake clusters into spatially continuous earthquake clusters through iterative calculations, and generates classification boundaries based on the sudden change gradient of the central depth of the earthquake clusters;

[0107] Optionally, step 104 may specifically include the following steps:

[0108] 1041. Merge the depth data in the valid data and the original frequencies recorded in the depth and frequency histogram according to corresponding equidistant grids to generate a joint input data set including depth data and event occurrence counts;

[0109] 1042. Performing unlabeled grouping processing on the joint input data set based on the Matlab software, wherein the grouping processing process minimizes the distribution difference of depth data and event occurrence counts within the same group by adjusting group boundaries, and forming spatially continuous earthquake clusters after the grouping processing;

[0110] Among them, step 1042 may specifically include the following processes: based on the Matlab software, the depth data intervals and event frequency intervals in the joint input data set are divided into discrete candidate groups according to the initial group boundaries; the boundary positions of the candidate groups are traversed, and the weighted difference values of the discreteness of the depth data distribution and the discreteness of the event frequency distribution within the same group are calculated; the adjacent groups with the smallest weighted difference values are automatically merged by iteratively adjusting the positions of the group boundaries to generate optimized spatially continuous groups; the groups with continuous depth data and uniform event frequency distribution in the spatially continuous groups are marked as spatially continuous earthquake clusters.

[0111] 1043. Performing geometric center calculation on all depth data points in the earthquake cluster to obtain a central depth value, wherein the geometric center calculation process is completed only based on the numerical distribution of the depth data points;

[0112] 1044. Compare the center depth values of adjacent earthquake clusters in sequence along the depth direction. If the absolute difference between the adjacent center depth values exceeds a mutation threshold, use the boundary between the adjacent earthquake clusters as a classification boundary.

[0113] In the above scheme, the combined input dataset is a fused dataset formed by aligning the effective depth data with the original frequency histogram records on a grid. Label-free grouping is a data grouping method that automatically partitions the data set through an algorithm without requiring manual labeling. Geometric center calculation calculates the center value based on the mean of the numerical distribution of the data points, without relying on weights or density. The mutation threshold is the critical value for determining the depth difference between adjacent earthquake clusters.

[0114] In this embodiment of the present application, each equidistant grid is first traversed through step 1041, and all depth data points belonging to that grid in the valid data are merged with the original frequencies recorded in the depth and frequency histogram according to the grid number, generating a structured joint input data set with the data structure: {grid number, depth data point set, original frequency}. For example, grid No. 52 contains 85 data points such as 5230 meters and 5255 meters. The frequency of grid No. 52 is 85 times, so the corresponding entry for grid 52 is {52, [5230, 5255, ...], 85}. This process uses Matlab to loop through the grid numbers and use array indexing to ensure strict data alignment.

[0115] Next, in step 1042, based on the combined input dataset, the group boundaries are dynamically adjusted using an iterative optimization algorithm. The objective function is to minimize the difference between the dispersion and frequency distribution of depth data points within the same group. Specifically, the following operations are performed: The standard deviation of the depth data for the candidate group is calculated as the square root of (Σ(depth value - mean)^2 / number of data points); the frequency coefficient of variation is calculated as standard deviation / mean frequency; the optimization goal is to minimize the sum of the depth standard deviation and the frequency coefficient of variation, and the group boundary positions are iteratively adjusted using a gradient descent method. For example, the 0-2000 meter grid is divided into shallow source clusters with a depth standard deviation of 150 meters and a frequency coefficient of variation of 0.2; the 5000-7000 meter grid is divided into deep source clusters with a standard deviation of 800 meters and a frequency coefficient of variation of 0.8.

[0116] Next, in step 1043, the arithmetic mean of all depth data points within each earthquake cluster is calculated using the formula: Center depth = Σ(depth data point) / number of data points. For example, if an earthquake cluster contains depth data points [1200, 1250, 1300], the calculated center value is (1200 + 1250 + 1300) / 3 = 1250 meters. This process is based solely on depth values and does not include frequency weights or density parameters.

[0117] Finally, in step 1044, the center depths of adjacent earthquake clusters are compared sequentially along the depth direction, and the absolute difference is calculated as |center depth value 1 - center depth value 2|. For example, if the center of the shallow source cluster is 1500 meters and the center of the deep source cluster is 2200 meters, the absolute difference is |1500 - 2200| = 700 meters. If the difference exceeds the preset mutation threshold, a classification boundary is generated at the intersection of the two clusters. For example, if the center of the shallow source cluster is 1500 meters and the center of the medium source cluster is 2200 meters, and the difference is 700 meters > 500 meters, a boundary is generated between 1500 and 2200 meters. If the difference in the center values of adjacent clusters is 300 meters < 500 meters, they are merged into the same cluster.

[0118] In practical applications, assuming that in a certain earthquake zone monitoring, in step 1041, the valid data includes 120 depth data points of 0-2000 meter grid, with a frequency of 120 times, and 30 data points of 5000-7000 meter grid, with a frequency of 30 times, to generate a joint data set entry.

[0119] Through step 1042, Matlab iteratively optimizes the grouping boundaries, dividing 0-2000 meters into shallow source clusters with a depth standard deviation of 180 meters and a frequency variation coefficient of 0.15; 5000-7000 meters into deep source clusters with a standard deviation of 950 meters and a frequency variation coefficient of 0.75; the middle 2000-5000 meters are not clustered separately because the frequency is less than 10 times and the distribution is discrete.

[0120] The central depth values are calculated in step 1043 : the central depth of the shallow source cluster = Σ (data points from 0 to 2000 meters) / 120 = 1050 meters; the central depth of the deep source cluster = Σ (data points from 5000 to 7000 meters) / 30 = 6200 meters.

[0121] The classification boundary is generated through step 1044. The center difference of 5150 meters between the shallow source cluster of 1050 meters and the deep source cluster of 6200 meters is greater than the mutation threshold of 500 meters, and the classification boundary is generated between 2000-5000 meters.

[0122] The overall solution of step 104 above solves the boundary ambiguity problem caused by the traditional method's reliance on artificial rules by jointly grouping depth data and statistical frequency; improves the geological interpretability of the classification results by generating sudden gradient boundaries; strictly retains the nonlinear characteristics of the original distribution to ensure that the classification boundaries are consistent with the actual geological activity pattern, significantly improving the automation level and accuracy of shallow crust earthquake source depth classification.

[0123] 105. After combining the classification boundary with the geographic coordinates, output an earthquake classification boundary map through the Matlab software.

[0124] Optionally, step 105 may specifically include the following steps:

[0125] 1051. Obtain geographic coordinate information of the shallow crust monitoring area, wherein the geographic coordinate information includes latitude and longitude data of the boundary points of the location range of the area corresponding to each earthquake cluster;

[0126] 1052. Bind each boundary point in the classification boundary to the latitude and longitude data of the corresponding position in the geographic coordinate information to generate a classification boundary coordinate point set including geographic coordinates;

[0127] 1053. Assign a unique color identifier to the earthquake cluster based on the difference in the center depth values of the earthquake cluster, and define a rule for assigning the color identifier as follows: if the center depth value of the earthquake cluster is less than a preset shallow source threshold, assign a first color identifier; if the center depth value of the earthquake cluster is greater than a preset deep source threshold, assign a second color identifier;

[0128] 1054. Based on the Matlab software, the classification boundary coordinate point set is regionally filled according to the geographic spatial position. During the filling process, the first color identifier and the second color identifier are respectively covered within the geographic coordinate range of the corresponding earthquake cluster to generate an earthquake classification boundary map containing color blocks.

[0129] In the above scheme, geographic coordinate information refers to the longitude and latitude data set of the geographic boundary points of the region corresponding to the earthquake cluster. The classification boundary coordinate point set is a dataset that binds the depth values of the classification boundary points to their actual geographic coordinates. The color code assignment rule defines the mapping rules for different color blocks based on the comparison between the depth value of the earthquake cluster center and the preset threshold. Area filling is the process of covering the color code with the geographic area corresponding to the earthquake cluster.

[0130] In this embodiment of the present application, step 1051 first extracts the longitude and latitude data of the earthquake cluster boundary points in the monitoring area from a geographic information database or map service. This is achieved by parsing a geographic information file, such as a Shapefile or GeoJSON. Matlab's shaperead function is used to read the file and extract the coordinate set of the boundary points. For example, the boundary points of a shallow source cluster are four coordinate points (longitude 116.3 to 116.5, latitude 39.6 to 39.8), which are stored as a structure array geo_data.Boundary = [116.3, 39.8; 116.5, 39.8; 116.5, 39.6; 116.3, 39.6].

[0131] Next, in step 1052, the classification boundary points are spatially mapped to geographic coordinates using a linear interpolation algorithm. Based on the depth boundary values and the crustal vertical layering model, the corresponding longitude and latitude coordinates are calculated. The formula is: longitude = base longitude + (depth - base depth) × longitude scaling factor; latitude = base latitude + (depth - base depth) × latitude scaling factor. For example, a depth of 2000 meters corresponds to a base longitude of 116.4 and a latitude of 39.7, resulting in the entry (2000 meters, 116.4, 39.7).

[0132] Next, define the color assignment rules in step 1053. Shallow source threshold: If the depth of the earthquake cluster center is less than or equal to the shallow source threshold, assign a red marker (RGB value 255, 0, 0). Deep source threshold: If the depth of the earthquake cluster center is greater than or equal to the deep source threshold, assign a blue marker (RGB value 0, 0, 255). The intermediate area: Assign a yellow marker (RGB value 255, 255, 0). For example, the center of the shallow source cluster at 1500 meters is assigned red, and the center of the deep source cluster at 6000 meters is assigned blue.

[0133] Finally, in step 1054, based on Matlab's m_map toolbox, polygonal areas are drawn according to geographic coordinates and filled with colors: the m_polygon function is used to draw the shallow source cluster area and fill it with red; the deep source cluster area is drawn and filled with blue; the intermediate area is drawn and filled with yellow; and finally, a classification boundary map with superimposed geographic grids and color blocks is output.

[0134] In practical applications, it is assumed that the coordinates of the boundary points of the shallow source cluster in a potential earthquake area are (115.8, 40.2), (116.2, 40.2), (116.2, 39.8), and (115.8, 39.8), which are obtained by reading the GeoJSON file using Matlab.

[0135] In step 1052 , the bounding coordinates are generated and the 2000-meter depth boundary is mapped to the geographic coordinates (116.0, 40.0), and the 5000-meter depth boundary is mapped to (115.9, 39.9), generating entries (2000 meters, 116.0, 40.0) and (5000 meters, 115.9, 39.9).

[0136] The color identification is assigned in step 1053 , where the center of the shallow source cluster at 1800 meters is assigned red, the center of the deep source cluster at 5200 meters is assigned blue, and the center of the intermediate cluster at 3500 meters is assigned yellow.

[0137] Generate a map through step 1054, draw the map using Matlab, fill in red (115.8-116.2, 39.8-40.2), fill in blue (115.6-115.8, 39.6-39.8), fill in the transition area with yellow, and output the map.

[0138] The overall solution of step 105 above seamlessly connects the classification results to the actual geographic space through automated geographic coordinate mapping and color identification assignment, solving the inefficiency problem caused by the traditional method's reliance on manual parameter adjustment. It intuitively distinguishes the dominant areas of shallow-source and deep-source earthquake activity through color blocks, improves the efficiency of visual expression of the source depth distribution pattern, and provides highly recognizable technical support for earthquake risk assessment and emergency decision-making.

[0139] The following is a complete embodiment of steps 101 to 105:

[0140] like Figure 2 As shown in the figure, in a hypothetical earthquake monitoring area, a broadband acoustic sensor array is deployed in step 101. Fifty sensor nodes are deployed in a triangular grid topology, with a spacing of 1.5 kilometers between nodes and microsecond-level time synchronization achieved via GPS. During an earthquake event, the array synchronously collects raw data on the focal depth, including acoustic wave propagation time series information: Node A records the timestamp 10:00:00.000, and Node B records 10:00:00.350.

[0141] Step 102 divides the maximum exploration depth of 10 kilometers into 100 equally spaced grids, based on a preset depth resolution of 100 meters: 0-100 meters, 100-200 meters, and so on, to 9900-10000 meters. The 50,000 collected depth data points are mapped to the corresponding grids. For example, 5230 meters is mapped to grid 52, between 5200 and 5300 meters. The frequency of each grid is counted to generate a histogram, disabling interpolation to preserve nonlinear characteristics. Grid 50 has a frequency of 85, and the adjacent grid 51 has a frequency of 10.

[0142] Step 103 selects three adjacent sensor nodes, spaced 1.2 km, 1.5 km, and 1.8 km apart, to form a validation group. The acoustic propagation time differences for the same earthquake event, ΔTAB = 0.35 seconds and ΔTAC = 0.52 seconds, are calculated. Based on a wave velocity V = 5000 m / s, the inversion depths H_AB = 0.35 × 5000 × cos(30°) = 1515 meters, and H_AC = 0.52 × 5000 × cos(45°) = 1838 meters. With a preset tolerance of ±200 meters and a true depth of 1600 meters, an error of 85 meters deems HAB valid, while an error of 238 meters deems HAC invalid. Seven of the ten node groups are statistically valid, representing 70% of the total, and are therefore marked as valid data.

[0143] In step 104, the valid data and histogram frequency data are combined into a joint dataset. The clustering boundaries are iteratively optimized to create shallow source clusters ranging from 0 to 2000 meters, with a frequency of 120 and a standard deviation of 180 meters; and deep source clusters ranging from 5000 to 7000 meters, with a frequency of 30 and a standard deviation of 950 meters. Cluster center depths are calculated: 1050 meters for shallow source clusters and 6200 meters for deep source clusters. The depth difference between adjacent cluster centers is compared along the depth direction, and the difference is found to be greater than the mutation threshold of 500 meters. This generates a classification boundary between 2000 and 5000 meters.

[0144] Step 105 extracts the monitoring area boundary coordinates from the geographic database: longitudes 116.3-116.5 and latitudes 39.6-39.8 represent shallow source clusters; latitudes 116.1-116.3 and 39.2-39.4 represent deep source clusters. The classification boundary points (2000 meters) at 116.4, 39.7 and 5000 meters at 116.2, 39.3) are bound to the geographic coordinates. Colors are assigned according to the rules: shallow source clusters are red, deep source clusters are blue, and intermediate areas are yellow. A map is drawn using Matlab's m_map toolbox, outputting a classification boundary map containing color blocks: red areas cover major urban areas, and blue areas are located in uninhabited mountainous areas.

[0145] This application uses a broadband acoustic sensor array to synchronously collect raw data on the depth of the earthquake source, and combines it with a multi-node acoustic wave propagation time difference cross-validation mechanism to screen out valid data with inversion depth errors within a preset tolerance range, significantly reducing the error accumulation problem caused by the heterogeneity of the shallow crust medium; constructs a depth-frequency histogram based on a preset depth resolution and strictly retains the nonlinear characteristics of the original distribution, providing a high-fidelity statistical benchmark for classification processing; through unsupervised classification processing that integrates valid data and histogram frequency, dynamically divides spatial earthquake clusters and automatically generates classification boundaries based on the central depth mutation gradient, breaking through the boundary ambiguity and efficiency bottleneck caused by traditional methods relying on artificial rules; finally, the classification boundaries are dynamically bound to geographic coordinates, and the dominant areas of shallow and deep earthquakes are intuitively distinguished through color blocks, realizing the geographical visualization expression of the earthquake source depth distribution pattern, improving the decision-making support efficiency and geological interpretation reliability of earthquake activity analysis, and providing high-precision technical support for disaster risk assessment and emergency response.

[0146] Figure 3 A schematic diagram of a seismic activity analysis system based on Matlab software is provided for the embodiment of the present application. Figure 3 As shown, the system includes:

[0147] The acquisition module 31 is used to synchronously acquire raw data of earthquake focal depth through detection devices deployed in the shallow crust monitoring area;

[0148] A construction module 32 is configured to generate an equidistant grid according to a preset depth resolution from the original focal depth data using Matlab software, and directly construct a depth-frequency histogram using the equidistant grid. The horizontal axis of the depth-frequency histogram is the equidistant depth grid, and the vertical axis records the original frequency of earthquake events within each depth grid, while retaining the nonlinear characteristics of the depth distribution.

[0149] a calculation module 33 for calculating the acoustic wave propagation time difference of the same earthquake event between adjacent sensor nodes according to the positional relationship of adjacent sensor nodes in the detection device, and marking the depth data of the earthquake event as valid data if the inversion depth of the acoustic wave propagation time difference is within a preset tolerance range;

[0150] a processing module 34 for performing unsupervised classification processing on the valid data and the original frequencies recorded in the depth and frequency histogram based on the Matlab software, wherein the processing process iteratively calculates and automatically divides the earthquake clusters into spatially continuous clusters and generates classification boundaries based on the sudden change gradient of the central depth of the earthquake clusters;

[0151] The output module 35 is used to combine the classification boundary with the geographic coordinates and output the earthquake classification boundary map through the Matlab software.

[0152] Figure 3 The seismic activity analysis system based on MATLAB software can be executed Figure 1 The implementation principle and technical effects of the MATLAB software-based seismicity analysis method described in the illustrated embodiment will not be elaborated on here. The specific manner in which the various modules and units of the Matlab software-based seismicity analysis system in the above embodiment perform operations has been described in detail in the embodiments of the method and will not be elaborated on here.

[0153] In one possible design, Figure 3 The seismicity analysis system based on Matlab software of the embodiment shown can be implemented as a computing device, such as Figure 4 As shown, the computing device may include a storage component 41 and a processing component 42;

[0154] The storage component 41 stores one or more computer instructions, wherein the one or more computer instructions are called and executed by the processing component 42 .

[0155] The processing component 42 is used for the above Figure 1 The embodiment provides a seismic activity analysis method based on Matlab software.

[0156] The processing component 42 may include one or more processors to execute computer instructions to complete all or part of the steps in the above method. Of course, the processing component may also be implemented as one or more application-specific integrated circuits (ASICs), digital signal processors (DSPs), digital signal processing devices (DSPDs), programmable logic devices (PLDs), field programmable gate arrays (FPGAs), controllers, microcontrollers, microprocessors, or other electronic components to perform the above method.

[0157] The storage component 41 is configured to store various types of data to support operations at the terminal. The storage component can be implemented by any type of volatile or non-volatile memory device, or a combination thereof, such as static random access memory (SRAM), electrically erasable programmable read-only memory (EEPROM), erasable programmable read-only memory (EPROM), programmable read-only memory (PROM), read-only memory (ROM), magnetic memory, flash memory, magnetic disk, or optical disk.

[0158] Of course, a computing device may also include other components, such as input / output interfaces, display components, communication components, etc.

[0159] The input / output interface provides an interface between the processing component and the peripheral interface module, which can be an output device, an input device, etc.

[0160] The communication component is configured to facilitate, among other things, wired or wireless communications between the computing device and other devices.

[0161] Among them, the computing device can be a physical device or an elastic computing host provided by a cloud computing platform, etc. In this case, the computing device can refer to a cloud server, and the above-mentioned processing components, storage components, etc. can be basic server resources rented or purchased from the cloud computing platform.

[0162] The present application also provides a computer storage medium storing a computer program, wherein the computer program can achieve the above-mentioned Figure 1 The embodiment shown is a method for analyzing seismicity based on MATLAB software.

[0163] Those skilled in the art will clearly understand that, for the convenience and brevity of description, the specific working processes of the systems, devices and units described above can refer to the corresponding processes in the aforementioned method embodiments and will not be repeated here.

[0164] The device embodiments described above are merely illustrative. The units described as separate components may or may not be physically separate, and the components shown as units may or may not be physical units, i.e., they may be located in one location or distributed across multiple network units. Some or all of the modules may be selected based on actual needs to achieve the objectives of the present embodiment. Persons of ordinary skill in the art will be able to understand and implement the present invention without inventive effort.

[0165] Through the above description of the embodiments, those skilled in the art will clearly understand that each embodiment can be implemented using software plus a necessary general-purpose hardware platform, or of course, hardware. Based on this understanding, the essence of the above technical solution, or the portion that contributes to the prior art, can be embodied in the form of a software product. This computer software product can be stored in a computer-readable storage medium, such as ROM / RAM, a magnetic disk, or an optical disk, and includes a number of instructions for causing a computer device (such as a personal computer, server, or network device) to execute the methods described in each embodiment or certain portions of the embodiments.

[0166] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present application, rather than to limit them. Although the present application has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present application.

Claims

1. A seismicity analysis method based on Matlab software, characterized in that: include: The original data of earthquake focal depth are collected synchronously by using detection devices deployed in the shallow crust monitoring area; Based on Matlab software, the original focal depth data is converted into an equidistant grid at a preset depth resolution, and a depth-frequency histogram is directly constructed using the equidistant grid. The horizontal axis of the depth-frequency histogram is the equidistant depth grid, and the vertical axis records the original frequency of earthquake events within each depth grid, while retaining the nonlinear characteristics of the depth distribution. Calculating the acoustic wave propagation time difference of the same earthquake event between adjacent sensor nodes based on the positional relationship of adjacent sensor nodes in the detection device, and marking the depth data of the earthquake event as valid data if the inversion depth of the acoustic wave propagation time difference is within a preset tolerance range; Performing unsupervised classification processing on the valid data and the original frequencies recorded in the depth and frequency histogram based on the Matlab software, wherein the processing process automatically divides the spatially continuous earthquake clusters by iterative calculation, and generates classification boundaries according to the sudden change gradient of the central depth of the earthquake cluster; After combining the classification boundary with the geographic coordinates, the earthquake classification boundary map is output through the Matlab software; The Matlab software is used to perform unsupervised classification processing on the valid data and the original frequencies recorded in the depth and frequency histogram. The processing process automatically divides spatially continuous earthquake clusters by iterative calculation and generates classification boundaries based on the sudden change gradient of the center depth of the earthquake cluster, including: Merging the depth data in the valid data with the original frequencies recorded in the depth and frequency histogram according to corresponding equidistant grids to generate a joint input data set including the depth data and the number of event occurrences; performing unlabeled grouping processing on the joint input data set based on the Matlab software, wherein the grouping processing process minimizes the distribution difference of depth data and event occurrence times within the same group by adjusting the group boundaries, and forms spatially continuous earthquake clusters after the grouping processing; Performing geometric center calculation on all depth data points in the earthquake cluster to obtain a central depth value, wherein the geometric center calculation process is completed only based on the numerical distribution of the depth data points; The central depth values of adjacent earthquake clusters are compared in sequence along the depth direction. If the absolute difference between the adjacent central depth values exceeds the mutation threshold, the boundary between the adjacent earthquake clusters is used as the classification boundary.

2. The method according to claim 1, characterized in that The unlabeled grouping process is performed on the combined input data set based on the Matlab software. The grouping process minimizes the distribution difference of depth data and event occurrence times within the same group by adjusting the group boundaries. After the grouping process, spatially continuous earthquake clusters are formed, including: Dividing the depth data intervals and event frequency intervals in the joint input data set into discrete candidate groups according to initial group boundaries based on the Matlab software; Traversing the boundary positions of the candidate groups, calculating the weighted difference between the depth data distribution dispersion and the event frequency distribution dispersion within the same group; Iteratively adjusting the positions of the group boundaries so that adjacent groups with the smallest weighted difference values are automatically merged to generate optimized spatially continuous groups; The groups with continuous depth data and uniform event frequency distribution in the spatially continuous groups are marked as spatially continuous earthquake clusters.

3. The method according to claim 1, characterized in that The method further comprises calculating the acoustic wave propagation time difference between adjacent sensor nodes of the same earthquake event based on the positional relationship between adjacent sensor nodes in the detection device, and marking the depth data of the earthquake event as valid data if the inversion depth of the acoustic wave propagation time difference is within a preset tolerance range. determining a relative distance between adjacent sensor nodes based on the layout position of the detection device, and generating a sensor node group consisting of at least three adjacent sensor nodes according to the relative distance; For the same earthquake event, the arrival time of the acoustic wave signal recorded by each sensor node in the sensor node group is extracted, and the acoustic wave propagation time difference between adjacent sensor nodes is calculated based on the arrival time; According to the correspondence between the acoustic wave propagation time difference and the relative distance, converting and calculating the correspondence to obtain an inversion depth of the earthquake event in the sensor node group; if the difference between the inversion depth and a preset tolerance range is a non-positive number, determining that the inversion depth of the sensor node group is valid; The valid number of the inversion depths in all sensor node groups is counted, and if the ratio of the valid number to the total number of sensor node groups exceeds a set threshold, the depth data of the earthquake event is marked as valid data.

4. The method according to claim 1, wherein The method of combining the classification boundary with the geographic coordinates and outputting an earthquake classification boundary map through the Matlab software comprises: Obtaining geographic coordinate information of the shallow crust monitoring area, wherein the geographic coordinate information includes latitude and longitude data of the boundary points of the location range of the area corresponding to each earthquake cluster; Binding each boundary point in the classification boundary to the latitude and longitude data of the corresponding position in the geographic coordinate information to generate a classification boundary coordinate point set containing geographic coordinates; Assigning a unique color identifier to the earthquake cluster based on the difference in the center depth value of the earthquake cluster, and defining the color identifier assignment rule as follows: if the center depth value of the earthquake cluster is less than a preset shallow source threshold, a first color identifier is assigned; if the center depth value of the earthquake cluster is greater than a preset deep source threshold, a second color identifier is assigned; Based on the Matlab software, the classification boundary coordinate point set is regionally filled according to the geographic spatial position. During the filling process, the first color identifier and the second color identifier are respectively covered to the geographic coordinate range of the corresponding earthquake cluster to generate an earthquake classification boundary map containing color blocks.

5. The method according to claim 1, wherein The method comprises generating an equidistant grid based on the original focal depth data according to a preset depth resolution based on Matlab software, and directly constructing a depth and frequency histogram through the equidistant grid, including: Determining the span of an equidistant grid according to a preset depth resolution, and generating an equidistant grid sequence covering all the original data of the focal depth according to the span, wherein the span is calculated by the ratio between the maximum detection depth of the monitoring area and the preset depth resolution; Traversing the depth data points in the original data of the focal depth based on Matlab software, mapping the depth data points to the equidistant grid sequence, wherein the mapping process is completed by determining the grid span interval to which the depth data points belong; Counting the number of depth data points mapped in the equidistant grid sequence, using the statistical results as the number of event occurrences, and generating an initial histogram structure with the equidistant grid as the horizontal axis and the number of event occurrences as the vertical axis; The initial histogram structure is subjected to nonlinear feature retention processing based on Matlab software. The interpolation or smoothing operation on the number of event occurrences of adjacent equidistant grids is prohibited during the processing process, and the original statistical results are used as the vertical axis values of the processed histogram to form a depth and frequency histogram.

6. The method according to claim 3, characterized in that The converting and calculating the corresponding relationship between the acoustic wave propagation time difference and the relative distance to obtain the inversion depth of the seismic event in the sensor node group includes: Obtaining the acoustic wave propagation time difference and the corresponding relative distance of adjacent nodes in the sensor node group, and generating a corresponding relationship set of time difference and distance; Establishing a linear proportional relationship model between the acoustic wave propagation time difference and the relative distance according to each data point in the corresponding relationship set; Based on the linear proportional relationship model, the acoustic wave propagation time difference is proportionally converted into a target depth, and during the conversion process, the product of the relative distance and the acoustic wave propagation time difference is kept as a fixed conversion coefficient; An arithmetic mean is taken for the target depths converted from the seismic event in the sensor node group, and the inversion depth of the seismic event is output.

7. A Matlab-based seismicity analysis system, used to execute the Matlab-based seismicity analysis method according to any one of claims 1 to 6, characterized in that: include: An acquisition module is used to synchronously acquire raw data of earthquake focal depth through detection devices deployed in the shallow crust monitoring area; A construction module is used to generate an equidistant grid according to a preset depth resolution from the original focal depth data based on Matlab software, and directly construct a depth and frequency histogram using the equidistant grid. The horizontal axis of the depth and frequency histogram is the equidistant depth grid, and the vertical axis records the original frequency of earthquake events in each depth grid, while retaining the nonlinear characteristics of the depth distribution; a calculation module, configured to calculate the acoustic wave propagation time difference of the same earthquake event between adjacent sensor nodes according to the positional relationship of adjacent sensor nodes in the detection device, and mark the depth data of the earthquake event as valid data if the inversion depth of the acoustic wave propagation time difference is within a preset tolerance range; a processing module for performing unsupervised classification processing on the valid data and the original frequencies recorded in the depth and frequency histogram based on the Matlab software, wherein the processing process automatically divides the spatially continuous earthquake clusters by iterative calculation and generates classification boundaries according to the sudden change gradient of the central depth of the earthquake cluster; The output module is used to combine the classification boundary with the geographic coordinates and output the earthquake classification boundary map through the Matlab software.

8. A computing device, characterized in that It comprises a processing component and a storage component; the storage component stores one or more computer instructions; the one or more computer instructions are used to be called and executed by the processing component to implement a seismic activity analysis method based on Matlab software as described in any one of claims 1 to 6.

9. A computer storage medium, characterized in that A computer program is stored, and when the computer program is executed by a computer, the seismic activity analysis method based on Matlab software as described in any one of claims 1 to 6 is implemented.

Images