Seismic phase identification method of non-natural earthquake based on image processing

Abstract

The invention discloses a non-natural earthquake phase identification method based on image processing, and relates to the technical field of image processing, and the method comprises the following steps: capturing a continuous signal with the amplitude changing with time in a seismic wave propagation process in real time based on a sensor distributed by an earthquake monitoring station network, forming an original amplitude sequence, and carrying out the calculation of the original amplitude sequence; collected continuous signals are converted into three-channel images, and the image brightness is automatically adjusted according to the seismic wave data dynamic range. According to the method, a seismic wave multi-dimensional analysis system is constructed through a three-channel construction strategy, the first channel clearly shows the time domain difference between natural seismic progressive fluctuation and non-natural seismic pulse type sudden change, and the second channel accurately captures microsecond-level sudden change which is easy to neglect in traditional time domain analysis; the third channel clearly distinguishes natural earthquake low-frequency energy concentration and non-natural earthquake intermediate-frequency characteristic spectral lines, and the adaptive brightness adjustment technology aims at the intensity difference of multi-source earthquake signals, improves the distinguishing degree of non-natural and natural earthquake edges, and builds a recognition foundation firmly.

Description

Technical Field

[0001] The present invention relates to the technical field of image processing, and in particular to a method for identifying seismic phases of non-natural earthquakes based on image processing. Background Art

[0002] Non-natural earthquakes refer to earthquakes caused by human activities or geological engineering, including artificial blasting, underground nuclear tests, induced earthquakes and mine tremors. Natural earthquakes are mostly caused by plate movement and fault slip. The seismic phase characteristics are manifested as S-wave development, a wide spectrum and a deep focal depth. The phase identification of non-natural earthquakes plays a key role in earthquake monitoring, industrial safety, engineering assessment and scientific research. Its core value lies in the analysis of seismic wave characteristics.

[0003] At present, in the existing technology, seismic phase identification generally relies on a single time domain analysis method, which can only capture the amplitude information of seismic waves changing with time. It is difficult to simultaneously obtain the time domain waveform contour, frequency domain energy distribution and microsecond-level mutation characteristics, resulting in the subtle differences between natural and non-natural earthquakes being difficult to effectively distinguish. At the same time, natural earthquake signals are usually weak, and details are easily lost under fixed threshold processing; and non-natural earthquake signals are strong, which can easily cause pixel saturation problems, making it impossible to accurately present waveform characteristics, thereby reducing the accuracy of seismic phase identification.

[0004] Therefore, a seismic phase identification method for non-natural earthquakes based on image processing is proposed to solve the above problems. Summary of the Invention

[0005] The main purpose of the present invention is to provide a method for identifying seismic phases of non-natural earthquakes based on image processing to solve the problems raised in the above background.

[0006] To achieve the above-mentioned purpose, the technical solution adopted by the present invention is: a method for identifying seismic phases of non-natural earthquakes based on image processing, comprising the following steps: Step 1: Based on the sensors distributed across the seismic monitoring network, the system captures the continuous signal of the amplitude changing over time during the propagation of seismic waves in real time, forming an original amplitude sequence. The collected continuous signal is then converted into a three-channel image, and the image brightness is automatically adjusted based on the dynamic range of the seismic wave data. Step 2: Use the Canny operator to perform edge detection on the three-channel image. The Canny operator determines the edge binary map by calculating the gradient and double threshold screening. For the edge points in the edge binary map, the gradient direction is divided into multiple intervals based on the gradient direction. The proportion of edge points in each interval is counted. The difference between the maximum interval proportion and the average proportion is used to measure the direction concentration. The difference is compared with the set threshold and a preliminary mark is output. Step 3: Receive preliminary labels, perform STFT transform on the original amplitude sequence, generate a time-frequency diagram, extract frequency domain features to divide multiple frequency bands, extract frequency band concentration and energy peak features, and output the non-natural earthquake identification results through multi-feature weighted fusion judgment.

[0007] Preferably, the continuous signal collected in step 1 is converted into a three-channel image, and the specific implementation steps are as follows: S1: The first channel is constructed, the original amplitude sequence is mapped, and the original amplitude sequence is normalized to the interval [0, 255] as the pixel value of the first channel image; S2: Second channel construction, performing first-order difference operation on the original amplitude sequence to generate a difference sequence, performing the same normalization on the difference sequence, and mapping it into the second channel image; S3: The third channel is constructed. The seismic wave signal is subjected to short-time Fourier transform in different time periods. The Hanning window and window length are set to obtain the frequency spectrum sequence in different time periods. The frequency spectrum is normalized to the interval [0, 255] and mapped into the third channel image.

[0008] Preferably, in step 1, the image brightness is automatically adjusted according to the dynamic range of the seismic wave data, and the specific implementation steps are as follows: S4: For the original amplitude sequence, difference sequence and frequency spectrum sequence of the three-channel image, traverse all elements in the sequence respectively, and select the internal maximum and minimum values ​​as the global maximum and global minimum values ​​of the entire sequence; S5: A dynamic range is set based on a global maximum value and a global minimum value, and a linear scaling algorithm is used to remap the original pixel values ​​according to the dynamic range.

[0009] Preferably, in step 2, edge detection is performed on the three-channel image using the Canny operator, and the Canny operator determines the edge binary image by calculating the gradient and performing double threshold screening. The specific implementation steps are as follows: S6: Based on the three-channel image in step 1, perform Canny edge detection on each channel image independently; S7: Convolution of each channel image with a Gaussian kernel; S8: Calculate the gradient magnitude map in the X and Y directions using the Sobel operator, and calculate the gradient magnitude and gradient direction; S9: traverse the gradient magnitude map, and for each pixel, compare the magnitude of the two adjacent pixels in the gradient direction. If it is greater than the magnitude of the adjacent pixel, it is retained as a candidate edge, otherwise it is suppressed; S10: By setting a high threshold T1 and a low threshold T2, and comparing them with the gradient amplitude, if the gradient amplitude is greater than or equal to the high threshold T1 and is a candidate edge, it is marked as a strong edge and directly retained. If the gradient amplitude is less than or equal to the high threshold T1, greater than or equal to the low threshold T2 and is a candidate edge, and there is a strong edge in the field, it is retained to obtain the edge binary map of each channel image.

[0010] Preferably, in step 2, based on the gradient direction, the gradient direction is divided into multiple intervals, the edge point ratio of each interval is counted, and the direction concentration is measured by the difference between the maximum interval ratio and the average ratio. The specific steps are as follows: S11: Based on the edge binary image of each single channel, the edge point coordinates are extracted, and based on the gradient direction in S8, the corresponding gradient direction value is extracted; S12: Divide the gradient direction into N intervals, count the number of edge points in each interval, and calculate the interval ratio; S13: Based on the interval proportion, calculate the average proportion and traverse all interval proportions at the same time to find the maximum interval proportion, and use the difference between the maximum interval proportion and the average proportion as the maximum proportion difference.

[0011] Preferably, in step 2, the data is compared with a set threshold value and a preliminary mark is output. The specific steps are as follows: S14: setting a directional concentration threshold T3 based on experiments; S15: If the maximum proportion difference is greater than the directional concentration threshold T3, the edge of the channel is preliminarily marked as a non-natural earthquake feature edge. If the maximum proportion difference is less than or equal to the directional concentration threshold T3, it is marked as a natural earthquake or noise edge.

[0012] Preferably, in step three, the energy peak characteristics include calculating peak energy, average energy and peak significance.

[0013] Preferably, in step three, the preliminary mark is received, the original amplitude sequence is subjected to STFT transformation, and the window length is set to generate a time-frequency diagram. The specific steps are as follows: S16: receiving preliminary marks, inputting the corresponding original amplitude sequence, and adjusting the window length, Hanning window, and overlap rate according to the main frequency of the seismic wave; S17: Generate a time-frequency diagram using the short-time Fourier transform formula.

[0014] Preferably, in step three, frequency domain features are extracted to divide multiple frequency bands, and frequency band concentration and energy peak features are extracted. The specific steps are as follows: S18: According to the physical characteristics of seismic waves, it is divided into low frequency band, medium frequency band and high frequency band; S19: For each frequency band, calculate the proportion of energy in the total energy, and calculate the difference between the proportion of low-frequency band energy in the total energy and the proportion of high-frequency band energy in the total energy as the frequency band concentration; S20: Calculate the peak energy and the average energy by finding the frequency of the energy peak in the time-frequency graph, and calculate the peak significance based on the peak energy and the average energy.

[0015] Preferably, in step 3, the weighted fusion judgment is performed by frequency band concentration, peak significance and direction concentration, and the specific steps are as follows: S21: By setting weights for frequency band concentration, peak significance, and direction concentration, the combined scores are obtained. S22: Based on the experimental calibration of the comprehensive threshold T4, if the comprehensive score is greater than the comprehensive threshold T4, it is judged to be a non-natural earthquake. If the comprehensive score is less than the comprehensive threshold T4, it is judged to be a natural earthquake or noise.

[0016] The present invention has the following beneficial effects: 1. In the present invention, a three-channel construction strategy is used to construct a multi-dimensional analysis system for seismic waves. The first channel completely retains the original amplitude waveform, clearly presenting the time domain difference between the gradual fluctuation of natural earthquakes and the pulse-like mutation of non-natural earthquakes. The second channel converts the steep peaks and troughs produced by non-natural events such as blasting into high-contrast pixel areas through first-order difference operations, accurately capturing microsecond-level mutations that are easily overlooked by traditional time domain analysis. The third channel uses STFT time-frequency transformation to convert frequency components into spatial distribution, clearly distinguishing the low-frequency energy concentration of natural earthquakes from the medium-frequency characteristic spectrum lines of non-natural earthquakes, laying the foundation for frequency band concentration analysis. The adaptive brightness adjustment technology targets the intensity differences of multi-source seismic signals, and dynamically adapts through global maximum statistics and linear scaling to improve the edge distinction between non-natural and natural earthquakes, laying a solid foundation for identification.

[0017] 2. In the present invention, the Canny operator is used to independently detect the edge of three-channel images, deeply mining the multidimensional features of seismic waves. In the preprocessing stage, Gaussian kernel convolution effectively filters sensor noise and environmental interference. The Sobel operator calculates the gradient amplitude map in the X and Y directions, converts the grayscale change into gradient intensity information, accurately locates the edge, divides the gradient direction interval and calculates the proportion, and uses the maximum proportion difference as a quantitative indicator to convert the visual edge into a mathematical feature. When the channel edge direction concentration difference is compared with the set threshold T3, the corresponding operation is output. Multiple links work together to improve the recognition accuracy, provide high-quality candidate areas for subsequent analysis, reduce computational redundancy, and enhance system robustness.

[0018] 3. In the present invention, a time-frequency diagram is generated by STFT transformation, breaking through the time domain limitations and presenting a panoramic view of the three-dimensional distribution of seismic waves. The low-frequency band, medium-frequency band and high-frequency band accurately capture the energy distribution law of the earthquake source. The frequency band concentration quantifies the aggregation tendency through the difference in the proportion of high and low frequency energy. The energy peak feature calculates the peak energy, average energy and significance, highlights the mutation details, and weightedly integrates the direction concentration, frequency band concentration and peak significance to construct a multi-dimensional discrimination system. Quantitative discrimination is based on the experimental calibration threshold T4, and multi-domain features are integrated throughout the process to improve recognition accuracy, avoid the risk of misjudgment of a single feature, build a solid barrier for non-natural earthquake identification, and promote the accuracy of seismic phase identification and differentiation. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] Figure 1 This is a flow chart of a method for identifying seismic phases of non-natural earthquakes based on image processing according to the present invention; Figure 2 This is a flowchart of step three of a method for identifying seismic phases of non-natural earthquakes based on image processing according to the present invention. DETAILED DESCRIPTION

[0020] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0021] See also Figure 1 and Figure 2 The present invention provides a technical solution: a method for identifying seismic phases of non-natural earthquakes based on image processing, comprising the following steps: Step 1: Based on the sensors distributed across the seismic monitoring network, the system captures the continuous signal of the amplitude changing over time during the propagation of seismic waves in real time, forming an original amplitude sequence. The collected continuous signal is then converted into a three-channel image, and the image brightness is automatically adjusted based on the dynamic range of the seismic wave data. Step 2: Use the Canny operator to perform edge detection on the three-channel image. The Canny operator determines the edge binary map by calculating the gradient and double threshold screening. For the edge points in the edge binary map, the gradient direction is divided into multiple intervals based on the gradient direction. The proportion of edge points in each interval is counted. The difference between the maximum interval proportion and the average proportion is used to measure the direction concentration. The difference is compared with the set threshold and a preliminary mark is output. Step 3: Receive preliminary labels, perform STFT transform on the original amplitude sequence, generate a time-frequency diagram, extract frequency domain features to divide multiple frequency bands, extract frequency band concentration and energy peak features, and output the non-natural earthquake identification results through multi-feature weighted fusion judgment.

[0022] The continuous signal collected in step 1 is converted into a three-channel image. The specific implementation steps are as follows: S1: The first channel is constructed, the original amplitude sequence is mapped, and the original amplitude sequence is normalized to the interval [0, 255] as the pixel value of the first channel image; S2: Second channel construction, performing first-order difference operation on the original amplitude sequence to generate a difference sequence, performing the same normalization on the difference sequence, and mapping it into the second channel image; Specifically, the first-order difference calculation formula is: ; in, represents the amplitude difference between adjacent moments, Represents the first The amplitude value at the moment; Generate difference sequence , and perform normalization.

[0023] S3: The third channel is constructed. The seismic wave signal is subjected to short-time Fourier transform in different time periods. The Hanning window and window length are set to obtain the frequency spectrum sequence in different time periods. The frequency spectrum is normalized to the interval [0, 255] and mapped into the third channel image.

[0024] Specifically, the calculation formula is: ; in, Represents the window index, represents the frequency point index, represents the Hanning window, represents the complex exponential term of the Fourier transform, Represents the spectrum amplitude; by calculating the spectrum amplitude, a frequency spectrum sequence is generated and normalized.

[0025] In step 1, the image brightness is automatically adjusted based on the dynamic range of the seismic wave data. The specific implementation steps are as follows: S4: For the original amplitude sequence, difference sequence and frequency spectrum sequence of the three-channel image, traverse all elements in the sequence respectively, and select the internal maximum and minimum values ​​as the global maximum and global minimum values ​​of the entire sequence; S5: A dynamic range is set based on a global maximum value and a global minimum value, and a linear scaling algorithm is used to remap the original pixel values ​​according to the dynamic range.

[0026] Specifically, the linear scaling algorithm formula is: ; in, represents the original pixel value, Represents the pixel value after scaling, represents the difference between the global maximum and the global minimum, represents the global minimum.

[0027] In step 2, the Canny operator is used to perform edge detection on the three-channel image, and the Canny operator determines the edge binary image by calculating the gradient and double threshold screening. The specific implementation steps are as follows: S6: Based on the three-channel image in step 1, perform Canny edge detection on each channel image independently; S7: Convolution of each channel image with a Gaussian kernel; S8: Calculate the gradient magnitude map in the X and Y directions using the Sobel operator, and calculate the gradient magnitude and gradient direction; Specifically, the gradient amplitude calculation formula is: ; in, and Represents the Sobel operator to calculate the gradient magnitude map in the X and Y directions, represents the gradient amplitude; The gradient direction calculation formula is: ; in, Represents the gradient direction, where the four main directions are 0 degrees, 45 degrees, 90 degrees, and 135 degrees.

[0028] S9: traverse the gradient magnitude map, and for each pixel, compare the magnitude of the two adjacent pixels in the gradient direction. If it is greater than the magnitude of the adjacent pixel, it is retained as a candidate edge, otherwise it is suppressed; S10: By setting a high threshold T1 and a low threshold T2, and comparing them with the gradient amplitude, if the gradient amplitude is greater than or equal to the high threshold T1 and is a candidate edge, it is marked as a strong edge and directly retained. If the gradient amplitude is less than or equal to the high threshold T1, greater than or equal to the low threshold T2 and is a candidate edge, and there is a strong edge in the field, it is retained to obtain the edge binary map of each channel image.

[0029] In step 2, based on the gradient direction, the gradient direction is divided into multiple intervals. The proportion of edge points in each interval is counted, and the difference between the maximum interval proportion and the average proportion is used to measure the direction concentration. The specific steps are as follows: S11: Based on the edge binary image of each single channel, the edge point coordinates are extracted, and based on the gradient direction in S8, the corresponding gradient direction value is extracted; S12: Divide the gradient direction into N intervals, count the number of edge points in each interval, and calculate the interval ratio; Specifically, the interval ratio calculation formula is: ; in, Represents the interval ratio, Representative The number of edge points contained in the interval, Represents the total number of edge points; S13: Based on the interval proportion, calculate the average proportion and traverse all interval proportions at the same time to find the maximum interval proportion, and use the difference between the maximum interval proportion and the average proportion as the maximum proportion difference.

[0030] In step 2, the data is compared with the set threshold and a preliminary mark is output. The specific steps are as follows: S14: setting a directional concentration threshold T3 based on experiments; S15: If the maximum proportion difference is greater than the directional concentration threshold T3, the edge of the channel is preliminarily marked as a non-natural earthquake feature edge. If the maximum proportion difference is less than or equal to the directional concentration threshold T3, it is marked as a natural earthquake or noise edge.

[0031] In step three, energy peak features include calculating peak energy, average energy and peak significance.

[0032] In step 3, the preliminary markup is received, the original amplitude sequence is subjected to STFT transformation, and the window length is set to generate a time-frequency diagram. The specific steps are as follows: S16: receiving preliminary marks, inputting the corresponding original amplitude sequence, and adjusting the window length, Hanning window, and overlap rate according to the main frequency of the seismic wave; S17: Generate time-frequency diagram through short-time Fourier transform formula; Specifically, the calculation formula is: ; in, Represents the result of short-time Fourier transform, reflecting the spectrum characteristics of the signal in different time windows and different frequency points. represents the Hanning window, represents the original amplitude sequence, represents the window length, represents the time window index, Represents the frequency point index; right Take the square of the modulus, with the time window index as the horizontal axis and the frequency point index as the vertical axis. Map to color and generate time-frequency graph.

[0033] In step 3, frequency domain features are extracted to divide the frequency bands into multiple layers, and the frequency band concentration and energy peak features are extracted. The specific steps are as follows: S18: According to the physical characteristics of seismic waves, it is divided into low frequency band, medium frequency band and high frequency band; Specifically, the low frequency band is 0-10 Hz, the middle frequency band is 10-30 Hz, and the high frequency band is above 30 Hz.

[0034] S19: For each frequency band, calculate the proportion of energy in the total energy, and calculate the difference between the proportion of low-frequency band energy in the total energy and the proportion of high-frequency band energy in the total energy as the frequency band concentration; Specifically, the formula for calculating the proportion of energy to total energy is: ; in, Represents the ratio of energy to total energy, Representative The set of frequency points corresponding to the frequency band, represents the result of short-time Fourier transform, Represents the frequency point, Represents a time window.

[0035] S20: Calculate the peak energy and the average energy by finding the frequency of the energy peak in the time-frequency graph, and calculate the peak significance based on the peak energy and the average energy.

[0036] Specifically, peak energy The calculation formula is as follows: ; Average energy The calculation formula is as follows: ; Peak significance The calculation formula is as follows: ; In step 3, the weighted fusion judgment is performed through the frequency band concentration, peak significance and direction concentration. The specific steps are as follows: S21: By setting weights for frequency band concentration, peak significance, and direction concentration, the combined scores are obtained. S22: Based on the experimental calibration of the comprehensive threshold T4, if the comprehensive score is greater than the comprehensive threshold T4, it is judged to be a non-natural earthquake. If the comprehensive score is less than the comprehensive threshold T4, it is judged to be a natural earthquake or noise.

[0037] In the present invention, a method for identifying seismic phases of non-natural earthquakes based on image processing is provided. The three-channel construction strategy achieves a multi-dimensional representation of seismic wave characteristics. The first channel preserves the original amplitude waveform, fully recording the time-domain profile of the seismic wave. This channel allows for a direct contrast between the gradual fluctuations of natural earthquakes and the pulse-like sudden changes of unnatural earthquakes. The second channel enhances the amplitude change rate through first-order difference operations, converting the sharp peaks and troughs produced by unnatural events such as explosions into high-contrast pixel regions, effectively highlighting microsecond-level sudden changes that are easily overlooked in traditional time-domain analysis. The third channel utilizes the STFT time-frequency transform to convert frequency components into spatial distributions, forming a characteristic distinction between the low-frequency energy concentration of natural earthquakes and the mid-frequency characteristic spectral lines of unnatural earthquakes, providing a basis for subsequent frequency band concentration analysis. Adaptive brightness adjustment technology, through global maximum statistics and linear scaling, effectively addresses the identification challenges posed by intensity differences between multiple sources of seismic signals. In natural earthquake scenarios, edge details are preserved after brightness stretching of weak signals, avoiding the information loss caused by traditional fixed threshold processing. In strong earthquake signals, dynamic range compression prevents pixel saturation and ensures that sudden changes are not submerged. This adaptive mapping based on the dynamic range of the data itself enables the subsequent dual-threshold screening of Canny edge detection to be performed in a more reasonable gradient space. The improvement in the quality of the gradient amplitude map directly translates into improved edge positioning accuracy, significantly enhancing the distinction between the regular edges unique to non-natural earthquakes and the random edges of natural earthquakes.

[0038] The Canny operator's independent edge detection mechanism for three-channel images realizes the multi-dimensional and refined capture of seismic wave characteristics. In the preprocessing stage, Gaussian kernel convolution effectively suppresses sensor noise and environmental interference. The Sobel operator converts the grayscale changes of the seismic wave image into gradient intensity information by calculating the gradient amplitude map in the X and Y directions. The gradient direction interval division and proportion statistics technology further deepen the quantitative expression of edge features. By dividing the gradient direction space into N intervals (such as 16 intervals, each interval is 22.5°), the system can capture the distribution law of edge direction and use the maximum proportion difference as a quantitative indicator to make the edge Edge features are transformed from vague visual impressions into precise mathematical features. When the directional concentration difference of a channel edge is greater than T3, its edge direction distribution shows a significant directional preference, which is consistent with the characteristic pattern of non-natural earthquakes. Otherwise, it is more likely to be a natural earthquake or noise. This data-driven threshold setting avoids the uncertainty of subjective judgment and elevates seismic phase identification from traditional empirical judgment to scientific decision-making based on quantitative features. The synergistic effect of multiple factors improves the accuracy of non-natural earthquake identification, provides high-quality candidate regions for time-frequency feature extraction, reduces computational redundancy, and enhances the robustness and efficiency of the entire identification system. From the perspective of feature mining, the STFT transform generates time-frequency diagrams, breaking the limitations of time-domain analysis and presenting the three-dimensional distribution of seismic waves. The layered frequency band division is consistent with the physical characteristics of seismic waves. The low-frequency, medium-frequency, and high-frequency bands accurately capture the energy distribution patterns of different earthquake sources. The frequency band concentration quantifies the energy accumulation tendency by the difference in the proportion of high and low frequency energy. The energy peak feature calculates the peak energy, average energy, and peak significance, highlighting energy mutations. Multiple features are extracted in parallel to fully cover the unique characteristics of non-natural earthquakes in the frequency domain. By weightedly integrating directional concentration, frequency band concentration, and peak significance, a multi-dimensional feature discrimination system is constructed. Based on the experimental calibration of the comprehensive threshold T4, quantitative discrimination is achieved, integrating multi-domain features, improving recognition accuracy, reducing the risk of misjudgment of a single feature, building a reliable technical barrier for non-natural earthquake identification, and facilitating the efficient distinction of seismic phase types.

[0039] It should be noted that, in this document, relational terms such as first and second, etc., are used only to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any actual relationship or order between these entities or operations. Moreover, the terms "comprises," "comprising," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that includes a list of elements includes not only those elements but also other elements not explicitly listed, or elements inherent to such process, method, article, or apparatus.

[0040] Although embodiments of the present invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and variations may be made to these embodiments without departing from the principles and spirit of the present invention, and the scope of protection claimed in the present invention is defined by the appended claims and their equivalents.

Claims

1. A method for identifying seismic phases of non-natural earthquakes based on image processing, characterized in that: The following steps are involved: Step 1: Based on the sensors distributed across the seismic monitoring network, the system captures the continuous signal of the amplitude changing over time during the propagation of seismic waves in real time, forming an original amplitude sequence. The collected continuous signal is then converted into a three-channel image, and the image brightness is automatically adjusted based on the dynamic range of the seismic wave data. Step 2: Use the Canny operator to perform edge detection on the three-channel image. The Canny operator determines the edge binary map by calculating the gradient and double threshold screening. For the edge points in the edge binary map, the gradient direction is divided into multiple intervals based on the gradient direction. The proportion of edge points in each interval is counted. The difference between the maximum interval proportion and the average proportion is used to measure the direction concentration. The difference is compared with the set threshold and a preliminary mark is output. Step 3: Receive preliminary labels, perform STFT transform on the original amplitude sequence, generate a time-frequency diagram, extract frequency domain features to divide multiple frequency bands, extract frequency band concentration and energy peak features, and output the non-natural earthquake identification results through multi-feature weighted fusion judgment.

2. The method for identifying seismic phases of non-natural earthquakes based on image processing according to claim 1, characterized in that: The continuous signal collected in step 1 is converted into a three-channel image. The specific implementation steps are as follows: S1: The first channel is constructed, the original amplitude sequence is mapped, and the original amplitude sequence is normalized to the interval [0, 255] as the pixel value of the first channel image; S2: Second channel construction, performing first-order difference operation on the original amplitude sequence to generate a difference sequence, performing the same normalization on the difference sequence, and mapping it into the second channel image; S3: The third channel is constructed. The seismic wave signal is subjected to short-time Fourier transform in different time periods. The Hanning window and window length are set to obtain the frequency spectrum sequence in different time periods. The frequency spectrum is normalized to the interval [0, 255] and mapped into the third channel image.

3. The method for identifying seismic phases of non-natural earthquakes based on image processing according to claim 1, characterized in that: In step 1, the image brightness is automatically adjusted according to the dynamic range of the seismic wave data. The specific implementation steps are as follows: S4: For the original amplitude sequence, difference sequence and frequency spectrum sequence of the three-channel image, traverse all elements in the sequence respectively, and select the internal maximum and minimum values ​​as the global maximum and global minimum values ​​of the entire sequence; S5: A dynamic range is set based on a global maximum value and a global minimum value, and a linear scaling algorithm is used to remap the original pixel values ​​according to the dynamic range.

4. The method for identifying seismic phases of non-natural earthquakes based on image processing according to claim 1, characterized in that: In the second step, the Canny operator is used to perform edge detection on the three-channel image, and the Canny operator determines the edge binary image by calculating the gradient and double threshold screening. The specific implementation steps are as follows: S6: Based on the three-channel image in step 1, perform Canny edge detection on each channel image independently; S7: Convolution of each channel image with a Gaussian kernel; S8: Calculate the gradient magnitude map in the X and Y directions using the Sobel operator, and calculate the gradient magnitude and gradient direction; S9: traverse the gradient magnitude map, and for each pixel, compare the magnitude of the two adjacent pixels in the gradient direction. If it is greater than the magnitude of the adjacent pixel, it is retained as a candidate edge, otherwise it is suppressed; S10: By setting a high threshold T1 and a low threshold T2, and comparing them with the gradient amplitude, if the gradient amplitude is greater than or equal to the high threshold T1 and is a candidate edge, it is marked as a strong edge and directly retained. If the gradient amplitude is less than or equal to the high threshold T1, greater than or equal to the low threshold T2 and is a candidate edge, and there is a strong edge in the field, it is retained to obtain the edge binary map of each channel image.

5. The method for identifying seismic phases of non-natural earthquakes based on image processing according to claim 1, characterized in that: In step 2, based on the gradient direction, the gradient direction is divided into multiple intervals, the edge point ratio of each interval is counted, and the difference between the maximum interval ratio and the average ratio is used to measure the direction concentration. The specific steps are as follows: S11: Based on the edge binary image of each single channel, the edge point coordinates are extracted, and based on the gradient direction in S8, the corresponding gradient direction value is extracted; S12: Divide the gradient direction into N intervals, count the number of edge points in each interval, and calculate the interval ratio; S13: Based on the interval proportion, calculate the average proportion and traverse all interval proportions at the same time to find the maximum interval proportion, and use the difference between the maximum interval proportion and the average proportion as the maximum proportion difference.

6. The method for identifying seismic phases of non-natural earthquakes based on image processing according to claim 1, characterized in that: In step 2, the comparison is performed with the set threshold value and a preliminary mark is output. The specific steps are as follows: S14: setting a directional concentration threshold T3 based on experiments; S15: If the maximum proportion difference is greater than the directional concentration threshold T3, the edge of the channel is preliminarily marked as a non-natural earthquake feature edge. If the maximum proportion difference is less than or equal to the directional concentration threshold T3, it is marked as a natural earthquake or noise edge.

7. The method for identifying seismic phases of non-natural earthquakes based on image processing according to claim 1, characterized in that: In step three, the energy peak feature includes calculating the peak energy, average energy and peak significance.

8. The method for identifying seismic phases of non-natural earthquakes based on image processing according to claim 1, characterized in that: In step 3, the preliminary mark is received, the original amplitude sequence is subjected to STFT transformation, and the window length is set to generate a time-frequency diagram. The specific steps are as follows: S16: receiving preliminary marks, inputting the corresponding original amplitude sequence, and adjusting the window length, Hanning window, and overlap rate according to the main frequency of the seismic wave; S17: Generate a time-frequency diagram using the short-time Fourier transform formula.

9. The method for identifying seismic phases of non-natural earthquakes based on image processing according to claim 1, characterized in that: In step 3, frequency domain features are extracted to divide multiple frequency bands, and frequency band concentration and energy peak features are extracted. The specific steps are as follows: S18: According to the physical characteristics of seismic waves, it is divided into low frequency band, medium frequency band and high frequency band; S19: For each frequency band, calculate the proportion of energy in the total energy, and calculate the difference between the proportion of low-frequency band energy in the total energy and the proportion of high-frequency band energy in the total energy as the frequency band concentration; S20: Calculate the peak energy and the average energy by finding the frequency of the energy peak in the time-frequency graph, and calculate the peak significance based on the peak energy and the average energy.

10. The method for identifying seismic phases of non-natural earthquakes based on image processing according to claim 7, characterized in that: In step 3, the weighted fusion judgment of frequency band concentration, peak significance and direction concentration is performed. The specific steps are as follows: S21: By setting weights for frequency band concentration, peak significance, and direction concentration, the combined scores are obtained. S22: Based on the experimental calibration of the comprehensive threshold T4, if the comprehensive score is greater than the comprehensive threshold T4, it is judged to be a non-natural earthquake. If the comprehensive score is less than the comprehensive threshold T4, it is judged to be a natural earthquake or noise.

Images