Seismic observation data automatic calibration method based on median filtering

Abstract

The invention discloses a seismic observation data automatic calibration method based on median filtering, and the method comprises the steps: firstly analyzing an original signal, determining a time window of a seismic signal, extracting a signal segment in the time window of the seismic signal, carrying out the band-pass filtering of the signal segment according to the type of the original signal, and obtaining an initial band-pass filtering signal; processing the initial band-pass filtering signal by using median filtering windows with different lengths, calculating indexes such as a signal-to-noise ratio and an effective signal retention rate, and screening an optimal initial window through second-order difference and sliding standard deviation; and processing the initial band-pass filtering signal by using an optimal signal-to-noise ratio window to obtain an initial median filtering signal, calculating a point-by-point correlation coefficient of the initial band-pass filtering signal and a residual error, carrying out regional dynamic window adjustment according to the correlation coefficient, carrying out median filtering processing again, and finally carrying out band-pass filtering processing again to obtain a final calibration signal. According to the design, the selected length of the window is not fixed, signals of different sampling rates and observation frequency bands can be adapted, the universality is high, and the noise reduction effect and effective signal retention are considered at the same time.

Description

TECHNICAL FIELD

[0001] The present application relates to a calibration method of seismic observation data, belonging to the field of seismic observation data processing and calibration, and particularly relates to an automatic calibration method of seismic observation data based on median filtering. BACKGROUND

[0002] In seismic observation, due to the complexity of the observation environment, a large amount of interference is often mixed in the original signal recorded by the seismic observation equipment. In order to eliminate the interference in the original signal, the original signal needs to be calibrated, and the filtering method is the main choice for such original signal calibration. Common filtering methods include band-pass filtering, high-pass and low-pass filtering, median filtering, etc. Among them, the core principle of median filtering is to set a fixed length of filtering window near the original signal data, and replace the original data with the median value in the filtering window to eliminate accidental errors. This traditional median filtering method has certain limitations, because the length of the filtering window selected in this method is relatively fixed, and the universality is poor, which cannot adapt to the original signals produced by observation equipment with different sampling rates and different observation frequency ranges, and the length of the filtering window is not appropriate, which may lead to incomplete noise reduction of the original signal or loss of effective signal.

[0003] The Chinese patent application with the application number CN202110102152.2 and the application date of January 26, 2021 discloses a square wave suppression method for seismic records, and the core technical solution is as follows: setting a mean value constraint unit, determining a filtering threshold by slidingly calculating the local signal mean value of the seismic record to avoid excessive filtering of effective signals; configuring a median filtering unit to slide and filter the seismic record containing square wave interference in segments with a preset fixed window, replacing the center data with the median value of the data in the window to achieve square wave noise suppression; a symmetrical extension device is provided to symmetrically expand the data at the beginning and end of the seismic record to prevent edge distortion during filtering; at the same time, noise energy detection (quantifying noise suppression effect) and waveform fidelity detection (ensuring effective signal retention) modules are provided, and finally the filtered waveform and related effect data are output. Although this invention can effectively suppress square wave interference, reduce edge distortion, and improve the recognition degree of effective seismic signals, it still has the following defects:

[0004] The length of the filtering window is relatively fixed, and the universality is poor.

[0005] The information disclosed in this BACKGROUND section is only intended to increase an understanding of the general context of the present application, and is not to be taken as an acknowledgement or any form of suggestion that this information forms prior art that is publicly known. SUMMARY

[0006] The application aims to overcome the problem that the length of the existing filter window is fixed and the universality is poor, and provide a seismic observation data automatic calibration method based on median filtering with a variable filter window length and strong universality.

[0007] To achieve the above object, the technical solution of the application is:

[0008] A seismic observation data automatic calibration method based on median filtering, the method comprising the following steps:

[0009] S1, first, select a paragraph with an amplitude variance lower than a preset threshold in an original signal as a background noise paragraph to calculate a background noise level; then detect a time point at which signal energy first continuously exceeds the background noise level by a first set threshold, and record the time point as a trigger time point; then, based on a seismic time and a distance from a seismic center provided by a seismic network, calibrate the trigger time point to determine an effective signal starting point; then, from the effective signal starting point, track the signal energy backward, and when the signal energy falls within the background noise level range and continuously falls below a second set threshold for a specified time length, record an end time point that meets the continuous time length condition as an effective signal ending point; finally, expand the effective signal starting point forward by a first set time length, and expand the effective signal ending point backward by a second set time length, to finally form a seismic signal time window containing pre-seismic background, complete seismic phase and post-seismic tail wave;

[0010] S2, first, extract a signal segment in the seismic signal time window in S1, and then perform band-pass filtering processing on the signal segment according to the type of the original signal in S1 to obtain an initial band-pass filtered signal;

[0011] S3, first, perform median filtering processing on the initial band-pass filtered signal obtained in S2 using median filtering windows with different window lengths in sequence to obtain filtered signals corresponding to each median filtering window; then, calculate the signal-to-noise ratio, effective signal retention rate and noise attenuation rate corresponding to each median filtering window; the signal-to-noise ratio is obtained by calculating the root mean square amplitude ratio of the effective segment and the noise segment of the signal after median filtering; the effective signal retention rate is obtained by calculating the root mean square amplitude ratio of the energy of the effective signal segment before and after median filtering; and the noise attenuation rate is obtained by calculating the root mean square amplitude ratio of the energy of the noise segment before and after median filtering;

[0012] S4, first, calculate the difference sequence of the signal-to-noise ratio, effective signal retention rate and noise attenuation rate corresponding to each median filtering window in S3, respectively:

[0013] For the signal-to-noise ratio, calculate the first-order difference sequence and the second-order difference sequence to obtain the first-order difference sequence of the signal-to-noise ratio and the second-order difference sequence of the signal-to-noise ratio;

[0014] For the effective signal retention rate and the noise attenuation rate, second-order difference sequences are calculated to obtain the effective signal retention rate second-order difference sequence and the noise attenuation rate second-order difference sequence;

[0015] Then, the second-order difference sequences of the signal-to-noise ratio, the effective signal retention rate and the noise attenuation rate are calculated to obtain the signal-to-noise ratio sliding standard deviation sequence, the effective signal retention rate sliding standard deviation sequence and the noise attenuation rate sliding standard deviation sequence; and based on the sliding standard deviation sequences of the signal-to-noise ratio, the effective signal retention rate and the noise attenuation rate, screening is performed to preliminarily obtain the signal-to-noise ratio optimal initial window, the effective signal retention optimal initial window and the noise attenuation optimal initial window;

[0016] S5, the signal-to-noise ratio optimal initial window obtained in S4 is used to perform median filtering on the initial band-pass filtered signal obtained in S2 to obtain a preliminary median filtered signal, and then a difference between the initial band-pass filtered signal and the preliminary median filtered signal is calculated to obtain a residual signal;

[0017] S6, the signal-to-noise ratio optimal initial window obtained in S4 is used as a sliding unit, and then a point-by-point correlation coefficient calculation is performed on the initial band-pass filtered signal obtained in S2 and the residual signal obtained in S5 to obtain a correlation coefficient corresponding to each data point in the initial band-pass filtered signal;

[0018] S7, according to the numerical value of the correlation coefficient of each data point obtained in S6, the initial band-pass filtered signal obtained in S2 is again subjected to point-by-point median filtering:

[0019] If the correlation coefficient of a data point is greater than 0.99, the median filtering window adjusted based on the noise attenuation optimal initial window obtained in S4 is used to perform median filtering on the data point; if the correlation coefficient of a data point is less than or equal to 0.99, the median filtering window adjusted based on the effective signal retention optimal initial window obtained in S4 is used to perform median filtering on the data point, and after each data point in the initial band-pass filtered signal is processed, a secondary median filtered signal is obtained;

[0020] S8, the secondary median filtered signal obtained in S7 is subjected to band-pass filtering again to obtain a final calibration signal.

[0021] In S2, the specific manner of performing band-pass filtering on the signal segment according to the type of the original signal in S1 is that: for different original signal types with different propagation characteristics, corresponding band-pass filtering frequency bands are preset, and the original signal types at least include near-shock signals, far-shock signals and ordinary seismic signals.

[0022] The presetting manner of the corresponding band-pass filtering frequency bands is that:

[0023] When the original signal type is a near-shock signal, the frequency band of the band-pass filtering is set to be 0.05 Hz to 30 Hz;

[0024] When the original signal type is a far-shock signal, the frequency band of the band-pass filtering is set to be 0.025 Hz to 10 Hz;

[0025] When the original signal type is a common earthquake signal, the frequency band of the band-pass filtering is set to be 0.05 Hz to 10 Hz.

[0026] In the S3, the parameters of the median filtering window of different lengths are determined by the following method:

[0027] The minimum window length, the window step and the maximum window length are calculated based on the sampling rate of the seismic observation equipment and the high-frequency cutoff frequency of the band-pass filtering in the S2, and the minimum window length and the maximum window length are both odd numbers.

[0028] The calculation formula of the minimum window length is:

[0029]

[0030] In the formula, bw is the minimum window length, f is the sampling rate of the seismic observation equipment, h is the high-frequency cutoff frequency of the band-pass filtering in the S2, and round is the rounding function. If the calculation result of bw is even, the final value of bw is: the calculation result of bw plus 1. If the calculation result of bw is less than 3, the final value of bw is equal to 3. f

[0031] The calculation formula of the window step is:

[0032]

[0033] In the formula, sw is the window step, f is the sampling rate of the seismic observation equipment, h is the high-frequency cutoff frequency of the band-pass filtering in the S2, and round is the rounding function. If the calculation result of sw is odd, the final value of sw is: the calculation result of sw plus 1. If the calculation result of sw is less than 2, the final value of sw is equal to 2. f

[0034] The calculation formula of the maximum window length is:

[0035] ew=round(f×1),

[0036] ​​In the formula, ew is the maximum window length, f is the sampling rate of the seismic observation device; if the calculation result of ew is an even number, the final value of ew is: the calculation result of ew plus 1, and the final value of ew is greater than or equal to the final value of bw plus 20 times 2.

[0037] In the S3, the calculation formula of the signal-to-noise ratio corresponding to each median filtering window is:

[0038]

[0039] In the formula, SNR i represents the signal-to-noise ratio corresponding to the i th median filtering window, rms represents the root mean square calculation function, SG i represents the effective signal component in the filtered signal corresponding to the i th median filtering window, NS i represents the noise component in the filtered signal corresponding to the i th median filtering window.

[0040] The calculation formula of the effective signal retention rate corresponding to each median filtering window is:

[0041]

[0042] In the formula, SR i represents the effective signal retention rate corresponding to the i th median filtering window, rms represents the root mean square calculation function, SG i represents the effective signal component in the filtered signal corresponding to the i th median filtering window, SG0 represents the effective signal component in the initial band-pass filtered signal obtained in S2;

[0043] The calculation formula of the noise attenuation rate corresponding to each median filtering window is:

[0044]

[0045] In the formula, NR i represents the noise attenuation rate corresponding to the i th median filtering window, rms represents the root mean square calculation function, NS i represents the noise component in the filtered signal corresponding to the i th median filtering window, NS0 represents the noise component in the initial band-pass filtered signal obtained in S2.

[0046] In the S4, the second-order difference sequences of the signal-to-noise ratio, the effective signal retention rate and the noise attenuation rate are respectively calculated to obtain the signal-to-noise ratio sliding standard deviation sequence, the effective signal retention rate sliding standard deviation sequence and the noise attenuation rate sliding standard deviation sequence.

[0047] For the second-order difference sequence of signal-to-noise ratio, the second-order difference sequence of effective signal retention rate and the second-order difference sequence of noise attenuation rate, each is independently processed according to the following steps:

[0048] Step one, first define a sliding window, each sliding window is composed of 5 second-order difference results in the second-order difference sequence in position continuity, specifically: take the kth second-order difference result in the second-order difference sequence as the starting point, containing the kth, k+1th, k+2th, k+3th, k+4th second-order difference result;

[0049] Step two, starting from the first second-order difference result of the second-order difference sequence, when k=1, the window contains the first to the fifth second-order difference result, calculate the standard deviation of the 5 second-order difference results in the sliding window to get 1 sliding standard deviation data; then slide the sliding window by 1 position, when k=2, take the second to the sixth second-order difference result to form a new sliding window and calculate the standard deviation, when the starting point k of the window satisfies k+4 is greater than the total length of the second-order difference sequence, that is, the sliding window cannot completely contain 5 second-order difference results, stop calculating;

[0050] Step three, arrange all the calculated sliding standard deviation data in order of their corresponding sliding window starting positions to get the sliding standard deviation sequence corresponding to the second-order difference sequence.

[0051] In S4, the sliding standard deviation sequences based on the signal-to-noise ratio, the effective signal retention rate and the noise attenuation rate are respectively screened to obtain the specific ways of the signal-to-noise ratio optimal initial window, the effective signal retention optimal initial window and the noise attenuation optimal initial window.

[0052] For the sliding standard deviation sequence of signal-to-noise ratio, the sliding standard deviation sequence of effective signal retention rate and the sliding standard deviation sequence of noise attenuation rate, each is independently processed according to the following steps:

[0053] First, calculate the average value of the sliding standard deviation data ranked in the last 40% of the sliding standard deviation sequence, and then take twice the average value as the stationary threshold corresponding to the sliding standard deviation sequence;

[0054] Second, starting from the first sliding standard deviation data in the sliding standard deviation sequence, sequentially find out whether the three consecutive sliding standard deviation data are all less than the stationary threshold obtained in the first step, that is, whether the jth, j+1th and j+2th sliding standard deviation data are all less than the stationary threshold;

[0055] Third, among the three sliding standard deviation data that first satisfy the condition that the three consecutive sliding standard deviation data are all less than the stationary threshold found in the second step, the first, i.e. the jth sliding standard deviation data, is determined as the optimal initial window of the median filtering window corresponding to the index of the sliding standard deviation sequence.

[0056] If the correlation coefficient of the data point is greater than 0.99, the data point is subjected to median filtering processing by using the median filtering window adjusted based on the noise attenuation optimal initial window obtained in S4, and the adjustment manner is that the length of the median filtering window for the median filtering processing of the data point is adjusted to 1.1 times of the length of the noise attenuation optimal initial window.

[0057] If the correlation coefficient of the data point is less than or equal to 0.99, the data point is subjected to median filtering processing by using the median filtering window adjusted based on the effective signal reservation optimal initial window obtained in S4, and the adjustment manner is that the length of the median filtering window for the median filtering processing of the data point is adjusted to 0.9 times of the length of the effective signal reservation optimal initial window.

[0058] In S4, after the signal-to-noise ratio optimal initial window, the effective signal reservation optimal initial window and the noise attenuation optimal initial window are obtained, the signal-to-noise ratio first-order difference sequence is further judged.

[0059] If in the signal-to-noise ratio first-order difference sequence, there are 8 or more values less than 0 in 10 continuous first-order difference results, the previous median filtering window corresponding to the first first-order difference result less than 0 in the 10 continuous first-order difference results is determined as a critical window at which the signal-to-noise ratio starts to decrease with the increase of the window length; if the window length of the signal-to-noise ratio optimal initial window obtained initially is greater than the window length of the critical window, the final value of the window length of the signal-to-noise ratio optimal initial window is corrected to the window length of the critical window; if the window length of the effective signal reservation optimal initial window obtained initially is greater than the window length of the critical window, the final value of the window length of the effective signal reservation optimal initial window is corrected to the window length of the critical window.

[0060] Compared with the prior art, the present application has the following beneficial effects:

[0061] 1. The application is a kind of automatic calibration method of seismic observation data based on median filtering, first, the original signal is analyzed, the time window of effective seismic signal is determined, then the signal segment in the time window is extracted, then the signal segment is band pass filtered according to the type of original signal, and the initial band pass filtered signal is obtained;Then, for the initial band pass filtered signal, different length median filtering windows are processed in turn, then the signal-to-noise ratio, effective signal retention rate and noise attenuation rate of each median filtering window are calculated, then the second order difference sequence and the sliding standard deviation sequence are calculated through the three indicators, and the three best initial windows corresponding to the signal-to-noise ratio, effective signal retention and noise attenuation are selected, and the preliminary selection of core median filtering window is completed;Subsequent dynamic calibration stage, first, the signal-to-noise ratio of the best initial window is applied to the initial band pass filtered signal to do preliminary median filtering, and the initial band pass filtered signal is obtained and the residual signal is calculated, then the signal-to-noise ratio of the best initial window is taken as the sliding unit, the point-by-point correlation coefficient of the initial signal and the residual signal is calculated, then the correlation coefficient is divided into scenes according to the correlation coefficient, if the correlation coefficient is greater than 0.99 (mostly noise intensive area), the median filtering window adjusted based on the noise attenuation best initial window is used to process to strengthen noise reduction;If the correlation coefficient is less than or equal to 0.99 (mostly effective signal intensive area), the median filtering window adjusted based on the best initial window of effective signal retention is used to process to reduce signal loss, and the secondary median filtering signal is obtained after point-by-point processing;Finally, the secondary median filtering signal is band pass filtered again, and the final calibration signal is obtained.

[0062] First, the correlation coefficient is calculated point by point to distinguish the noise intensive area and the effective signal intensive area, and the filtering window is dynamically matched according to different areas: the length of the noise intensive area is adjusted based on the noise attenuation best initial window, and the length of the effective signal intensive area is adjusted based on the best initial window of effective signal retention, which avoids the disadvantages of "incomplete noise reduction" or "effective signal loss" caused by the fixed length of the filtering window in the traditional method.

[0063] Second, the length of the filtering window is not fixed in the whole process, but the window range is dynamically calculated combined with the sampling rate of the observation device and the signal observation frequency band, and then the best initial window is selected through multiple indicators such as signal-to-noise ratio and effective signal retention rate, and the band pass filtering is designed combined with the type of original signal, which can adapt to different sampling rates, different observation frequency bands and different types of seismic signals, and solves the poor universality problem caused by the fixed length of the filtering window in the existing one-dimensional median filtering improvement scheme.

[0064] Therefore, the application breaks through the limitation of the fixed length of the filtering window in the traditional median filtering method, the length of the filtering window is not fixed, the universality is strong, and the original signal produced by the observation device with different sampling rates and different observation frequency band ranges can be adapted, and the accurate noise reduction of seismic observation data and the reduction of effective signal loss are considered.

[0065] 2. In the automatic calibration method for seismic observation data based on median filtering, when performing band-pass filtering on the original signal, the type of the original signal needs to be judged in combination with the observation scene. If it is a near earthquake monitoring or artificial earthquake exploration scene (such as a regional seismic network capturing short-distance earthquake signals), after identifying the near earthquake signal, the signal segment within the effective earthquake signal time window is extracted, and a frequency band of 0.05Hz-30Hz is used for filtering to retain the rich high-frequency details of the near earthquake signal. If it is a teleseismic monitoring (such as a global seismic network capturing interplate earthquake signals), after identifying the teleseismic signal, a frequency band of 0.025Hz-10Hz is used for filtering to filter the naturally attenuated redundant high-frequency interference in the teleseismic signal. If it is a medium-short or medium-long distance earthquake monitoring, after identifying the ordinary earthquake signal, a frequency band of 0.05Hz-10Hz is used for filtering to adapt to its concentrated frequency characteristics. This kind of "signal type-quantitative frequency band" corresponding processing can not only avoid the loss of high-frequency details of near earthquake signals and the residual of invalid noise of teleseismic signals in general frequency bands, but also improve the purity of the initial band-pass filtered signal from the source, lay a high-quality foundation for the selection of the subsequent median filtering window, eliminate the subjectivity of artificial debugging through clear frequency band parameters, ensure the consistency of different observers and different equipment processing the same type of signal, and reduce the operation complexity. Therefore, by matching the clear band-pass filtering frequency band according to the type of the original signal, the accuracy of the original signal processing and the standardization and reproducibility of the method are effectively improved.

[0066] 3、The method for automatically calibrating seismic observation data based on median filtering, when processing the second-order difference sequence of signal-to-noise ratio, effective signal retention rate and noise attenuation rate, first calculates the sliding standard deviation: taking the second-order difference sequence of signal-to-noise ratio as an example, define five continuous second-order difference results as a sliding window, start calculating the standard deviation from k=1 (the window contains the first to fifth second-order difference results), slide bit by bit to k+4 greater than the total length of the sequence, and generate a complete sliding standard deviation sequence; then take twice the average value of the sliding standard deviation data in the last 40% of the sliding standard deviation sequence as the stationary threshold, and then find the first three consecutive sliding standard deviation data less than the stationary threshold from the first data, and take the median filtering window corresponding to the first data as the best initial window (for example, if the first three consecutive second-order difference results from the 12th to the 14th in the signal-to-noise ratio sequence meet the standard, then the window corresponding to the 12th second-order difference result is the best initial window of the signal-to-noise ratio). The window size of five continuous second-order difference results can avoid accidental error interference and accurately capture the fluctuation characteristics of the second-order difference sequence, ensuring that the sliding standard deviation sequence truly reflects the stability of the index change; taking twice the average value of the sliding standard deviation data in the last 40% as the stationary threshold and setting the condition that three consecutive sliding standard deviation data less than the stationary threshold meet the standard, can adapt to the fluctuation characteristics of different signals and exclude single accidental misjudgment, ensuring that the selected best initial window is stable and efficient, and providing a reliable basis for subsequent dynamic window adjustment. Therefore, through the standardized sliding standard deviation calculation and best initial window screening logic, the present application avoids subjective judgment errors and significantly improves the reliability of the filtering window screening and the scientificity of the method.

[0067] 4、The method for automatically calibrating seismic observation data based on median filtering, when processing the initial band-pass filtered signal point by point, the window needs to be adjusted according to the correlation coefficient: if the correlation coefficient of a data point is 0.995 (>0.99, determined as a noise-intensive area), the noise attenuation best initial window (assuming 21) is called, adjusted to 23 by 1.1 times, and the window is used to strengthen the elimination of intensive noise; if the correlation coefficient of a data point is 0.98 (<0.99, determined as an effective signal-intensive area), the effective signal retention best initial window (assuming 19) is called, adjusted to 17 by 0.9 times, and the window is used to reduce the loss of effective signal amplitude. This quantitative adjustment ratio realizes the fine filtering of "expanding the window to strongly reduce noise in the noise area and reducing the window to preserve details in the signal area", avoids incomplete noise reduction or signal loss caused by fixed length window or empirical window adjustment, and eliminates the randomness of manual adjustment through clear 1.1 times and 0.9 times parameters, ensuring the consistency of processing logic for the same type of data points in different scenarios and improving the stability of filtering effect. Therefore, through the quantitative adjustment of the median filtering window in different regions, the present application accurately balances the noise reduction effect and effective signal retention, further enhances the adaptability of the method and the stability of the filtering quality. BRIEF DESCRIPTION OF DRAWINGS

[0068] Figure 1 is a flowchart of the present application.

[0069] Figure 2 is an example graph of the best initial window selection in Example 1.

[0070] Figure 3 is a comparison graph of the signal-to-noise ratio in Example 1.

[0071] Figure 4 is a comparison graph of the noise signal in Example 1.

[0072] Figure 5 is a comparison graph of the effective signal in Example 1. DETAILED DESCRIPTION

[0073] The present application is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0074] Reference is made to Figure 1 — Figure 5 A median filter-based automatic calibration method for seismic observation data, the method comprising the following steps:

[0075] S1, first select a paragraph in the original signal with an amplitude variance below a preset threshold as a background noise paragraph to calculate the background noise level; then detect the time when the signal energy first continuously exceeds the background noise level by a first set threshold, and record it as the trigger time; then calibrate the trigger time based on the earthquake occurrence time and epicentral distance provided by the seismic network to determine the effective signal starting point, then start from the effective signal starting point, track the signal energy backward, when the signal energy falls within the background noise level range and continuously falls below a second set threshold for a specified time length, record the end time that meets the continuous time length condition as the effective signal termination point; finally, expand the effective signal starting point forward by a first set time length, and expand the effective signal termination point backward by a second set time length, to finally form a seismic signal time window containing pre-seismic background, complete seismic phase and post-seismic coda;

[0076] S2, first extract the signal segment within the seismic signal time window in S1, then perform band-pass filtering processing on the signal segment according to the type of the original signal in S1 to obtain an initial band-pass filtered signal;

[0077] S3, first, the initial band-pass filtered signal obtained in S2 is subjected to median filtering with median filtering windows of different window lengths in sequence to obtain filtered signals corresponding to each median filtering window; then, the signal-to-noise ratio, the effective signal retention rate and the noise attenuation rate corresponding to each median filtering window are calculated; the signal-to-noise ratio is obtained by calculating the root mean square amplitude ratio of the effective section and the noise section of the signal after median filtering; the effective signal retention rate is obtained by calculating the root mean square amplitude ratio of the energy of the effective signal section before and after median filtering; the noise attenuation rate is obtained by calculating the root mean square amplitude ratio of the energy of the noise section before and after median filtering;

[0078] S4, first, the signal-to-noise ratio, the effective signal retention rate and the noise attenuation rate corresponding to each median filtering window in S3 are calculated respectively to obtain the difference sequences thereof;

[0079] For the signal-to-noise ratio, the first-order difference sequence and the second-order difference sequence are calculated to obtain the first-order difference sequence of the signal-to-noise ratio and the second-order difference sequence of the signal-to-noise ratio;

[0080] For the effective signal retention rate and the noise attenuation rate, the second-order difference sequences are calculated respectively to obtain the second-order difference sequence of the effective signal retention rate and the second-order difference sequence of the noise attenuation rate;

[0081] Then, the second-order difference sequences of the signal-to-noise ratio, the effective signal retention rate and the noise attenuation rate are calculated to obtain the corresponding signal-to-noise ratio sliding standard deviation sequence, the effective signal retention rate sliding standard deviation sequence and the noise attenuation rate sliding standard deviation sequence; then, based on the sliding standard deviation sequences of the signal-to-noise ratio, the effective signal retention rate and the noise attenuation rate, the signal-to-noise ratio, the effective signal retention rate and the noise attenuation rate are screened respectively to preliminarily obtain the best initial window of the signal-to-noise ratio, the best initial window of the effective signal retention and the best initial window of the noise attenuation;

[0082] S5, first, the initial band-pass filtered signal obtained in S2 is subjected to median filtering with the best initial window of the signal-to-noise ratio obtained in S4 to obtain a preliminary median filtered signal, and then the difference between the initial band-pass filtered signal and the preliminary median filtered signal is calculated to obtain a residual signal;

[0083] S6, first, the best initial window of the signal-to-noise ratio obtained in S4 is taken as a sliding unit, and then the initial band-pass filtered signal obtained in S2 and the residual signal obtained in S5 are subjected to point-by-point correlation coefficient calculation to obtain the correlation coefficient corresponding to each data point in the initial band-pass filtered signal;

[0084] S7, according to the numerical value of the correlation coefficient of each data point obtained in S6, the initial band-pass filtered signal obtained in S2 is subjected to point-by-point median filtering again:

[0085] If the correlation coefficient of the data point is greater than 0.99, the data point is subjected to median filtering processing based on the median filtering window of the noise attenuation optimal initial window obtained in S4; if the correlation coefficient of the data point is less than or equal to 0.99, the data point is subjected to median filtering processing based on the median filtering window of the effective signal reservation optimal initial window obtained in S4; and after each data point in the initial band-pass filtered signal is processed, a secondary median filtered signal is obtained.

[0086] S8, the secondary median filtered signal obtained in S7 is subjected to band-pass filtering processing again to obtain a final calibration signal.

[0087] In S2, the specific manner of the band-pass filtering processing of the signal segment according to the type of the original signal in S1 is that: for different original signal types of propagation characteristics, corresponding band-pass filtering frequency bands are preset, and the original signal types at least include near earthquake signals, far earthquake signals and ordinary earthquake signals.

[0088] The preset manner of the corresponding band-pass filtering frequency band is:

[0089] When the original signal type is a near earthquake signal, the frequency band of the band-pass filtering is set to be low frequency 0.05 Hz to high frequency 30 Hz.

[0090] When the original signal type is a far earthquake signal, the frequency band of the band-pass filtering is set to be low frequency 0.025 Hz to high frequency 10 Hz.

[0091] When the original signal type is an ordinary earthquake signal, the frequency band of the band-pass filtering is set to be low frequency 0.05 Hz to high frequency 10 Hz.

[0092] In S3, the parameters of the median filtering windows of different lengths are determined by the following manner:

[0093] Based on the sampling rate of the seismic observation equipment and the high frequency cutoff frequency of the band-pass filtering frequency band in S2, the minimum window length, the window step and the maximum window length are calculated, and the minimum window length and the maximum window length are both odd numbers.

[0094] The calculation formula of the minimum window length is:

[0095]

[0096] In the formula, bw is the minimum window length, f is the sampling rate of the seismic observation equipment, h is the high frequency cutoff frequency of the band-pass filtering in S2, and round is the rounding function; if the calculation result of bw is an even number, the final value of bw is: the calculation result of bw plus 1; if the calculation result of bw is less than 3, the final value of bw is equal to 3. f ​

[0097] The formula for calculating the window step size is:

[0098]

[0099] In the formula, sw is the window step size, f is the sampling rate of the seismic observation equipment, and h f is the high-frequency cutoff frequency of the bandpass filter in S2, and round is the rounding function; if the calculated result of sw is odd, then the final value of sw is: the calculated result of sw plus 1; if the calculated result of sw is less than 2, then the final value of sw is set to 2.

[0100] The formula for calculating the maximum window length is:

[0101] ew = round(f × 1),

[0102] In the formula, ew is the maximum window length, and f is the sampling rate of the seismic observation equipment; if the calculated result of ew is even, the final value of ew is: the calculated result of ew plus 1, and the final value of ew must be greater than or equal to 20 times the final value of bw plus 2.

[0103] In step S3, the formula for calculating the signal-to-noise ratio corresponding to each median filtering window is as follows:

[0104]

[0105] In the formula, SNR i Represents the signal-to-noise ratio corresponding to the i-th median filter window, rms represents the root mean square function, and SG... i NS represents the effective signal component in the filtered signal corresponding to the i-th median filter window. i This represents the noise component in the filtered signal corresponding to the i-th median filter window;

[0106] The formula for calculating the effective signal retention rate for each median filter window is as follows:

[0107]

[0108] In the formula, SR i The median filter window represents the effective signal retention rate, rms represents the root mean square function, and SG represents the effective signal retention rate. i SG0 represents the effective signal component in the filtered signal corresponding to the i-th median filtering window, and SG0 represents the effective signal component in the initial bandpass filtered signal obtained in S2.

[0109] The formula for calculating the noise attenuation rate for each median filter window is as follows:

[0110]

[0111] wherein NR i denotes the noise attenuation rate corresponding to the i-th median filter window, rms denotes the root mean square calculation function, NS i denotes the noise component in the filtered signal corresponding to the i-th median filter window, NS0 denotes the noise component in the initial band-pass filtered signal obtained in S2.

[0112] In S4, the second-order difference sequences of the signal-to-noise ratio, the effective signal retention rate and the noise attenuation rate are respectively calculated to obtain the sliding standard deviation sequences of the signal-to-noise ratio, the effective signal retention rate and the noise attenuation rate in the following manner:

[0113] For the second-order difference sequences of the signal-to-noise ratio, the effective signal retention rate and the noise attenuation rate, each is independently processed in the following manner:

[0114] Step 1: define a sliding window, each sliding window is composed of 5 consecutive second-order difference results in the second-order difference sequence, specifically: take the k-th second-order difference result in the second-order difference sequence as the starting point, and include the k-th, k+1-th, k+2-th, k+3-th and k+4-th second-order difference results;

[0115] Step 2: start from the first second-order difference result in the second-order difference sequence, i.e. when k = 1, the window includes the 1st to 5th second-order difference results, calculate the standard deviation of the 5 second-order difference results in the sliding window to obtain 1 sliding standard deviation data; then slide the sliding window by 1 position each time, i.e. when k = 2, take the 2nd to 6th second-order difference results to form a new sliding window and calculate the standard deviation, and stop the calculation when the starting point k of the window satisfies k+4 is greater than the total length of the second-order difference sequence, i.e. the sliding window cannot completely include 5 second-order difference results.

[0116] Step 3: arrange all the calculated sliding standard deviation data in the order of their corresponding sliding window starting positions to obtain the sliding standard deviation sequence corresponding to the second-order difference sequence.

[0117] In S4, the sliding standard deviation sequences of the signal-to-noise ratio, the effective signal retention rate and the noise attenuation rate are respectively screened to obtain the signal-to-noise ratio optimal initial window, the effective signal retention optimal initial window and the noise attenuation optimal initial window in the following manner:

[0118] For the sliding standard deviation sequences of the signal-to-noise ratio, the effective signal retention rate and the noise attenuation rate, each is independently processed in the following manner:

[0119] In the first step, an average value of the 40% of the sliding standard deviation data in the sliding standard deviation sequence is calculated, and twice the average value is taken as the stationary threshold corresponding to the sliding standard deviation sequence.

[0120] In the second step, starting from the first sliding standard deviation data in the sliding standard deviation sequence, it is determined whether the three consecutive sliding standard deviation data are less than the stationary threshold obtained in the first step, that is, whether the jth, j+1th and j+2th sliding standard deviation data are less than the stationary threshold.

[0121] In the third step, the first sliding standard deviation data corresponding to the median filtering window in which the three consecutive sliding standard deviation data are less than the stationary threshold for the first time is determined as the optimal initial window of the index corresponding to the sliding standard deviation sequence.

[0122] In the S7, if the correlation coefficient of the data point is greater than 0.99, the median filtering window adjusted based on the noise attenuation optimal initial window obtained in the S4 is used for median filtering processing of the data point, and the length of the median filtering window is adjusted to 1.1 times the length of the noise attenuation optimal initial window.

[0123] If the correlation coefficient of the data point is less than or equal to 0.99, the median filtering window adjusted based on the effective signal retention optimal initial window obtained in the S4 is used for median filtering processing of the data point, and the length of the median filtering window is adjusted to 0.9 times the length of the effective signal retention optimal initial window.

[0124] In the S4, after the signal-to-noise ratio optimal initial window, the effective signal retention optimal initial window and the noise attenuation optimal initial window are obtained, the signal-to-noise ratio first-order difference sequence needs to be judged.

[0125] If there are 8 or more values less than 0 in the 10 consecutive first-order difference results in the signal-to-noise ratio first-order difference sequence, the previous median filtering window corresponding to the first first-order difference result less than 0 in the 10 consecutive first-order difference results is determined as the critical window at which the signal-to-noise ratio starts to decrease with the increase of the window length. If the window length of the signal-to-noise ratio optimal initial window obtained initially is greater than the window length of the critical window, the final value of the window length of the signal-to-noise ratio optimal initial window is corrected to the window length of the critical window. If the window length of the effective signal retention optimal initial window obtained initially is greater than the window length of the critical window, the final value of the window length of the effective signal retention optimal initial window is corrected to the window length of the critical window.

[0126] Example 1:

[0127] Referring to Figure 1 — Figure 5The application discloses a kind of median filter-based seismic observation data automatic calibration method, the method comprises the following steps: S1, first select the paragraph of original signal in amplitude variance below preset threshold, as background noise section, to calculate background noise level;Then detect the time when signal energy first continuously exceeds background noise level by first set threshold, record it as trigger time;Again, based on the earthquake time and epicentral distance provided by seismic network, the trigger time is calibrated to determine the effective signal starting point, then from the effective signal starting point, signal energy is tracked back, when signal energy falls within the background noise level range, and continuously lower than second set threshold reaches specified time length, after meeting the duration time condition, the end time that satisfies the duration time condition is recorded as effective signal termination point;Finally, the effective signal starting point is extended first set time length forward, and the effective signal termination point is extended second set time length backward, finally form the seismic signal time window containing pre-seismic background, complete seismic phase and post-seismic tail wave;S2, first extract the signal segment in the seismic signal time window in S1, then according to the type of original signal in S1, the signal segment is band-pass filtered to obtain initial band-pass filtered signal;S3, first, the initial band-pass filtered signal obtained in S2 is subjected to median filtering with different window lengths in sequence to obtain the filtered signal corresponding to each median filtering window;Then calculate the signal-to-noise ratio, effective signal retention rate and noise attenuation rate corresponding to each median filtering window;The signal-to-noise ratio is obtained by calculating the root mean square amplitude ratio of the effective section and noise section after median filtering;The effective signal retention rate is obtained by calculating the root mean square amplitude ratio of the energy of the effective signal section before and after median filtering;The noise attenuation rate is obtained by calculating the root mean square amplitude ratio of the energy of the noise section before and after median filtering;S4, first, the second-order difference of the signal-to-noise ratio, effective signal retention rate and noise attenuation rate corresponding to each median filtering window in S3 is calculated respectively to obtain the corresponding signal-to-noise ratio second-order difference sequence, effective signal retention rate second-order difference sequence and noise attenuation rate second-order difference sequence;Then, the second-order difference sequences of signal-to-noise ratio, effective signal retention rate and noise attenuation rate are calculated to obtain the corresponding signal-to-noise ratio sliding standard deviation sequence, effective signal retention rate sliding standard deviation sequence and noise attenuation rate sliding standard deviation sequence, and then the signal-to-noise ratio, effective signal retention rate and noise attenuation rate sliding standard deviation sequences are screened respectively to obtain the signal-to-noise ratio best initial window, effective signal retention best initial window and noise attenuation best initial window;S5, first, the signal-to-noise ratio best initial window obtained in S4 is used to perform median filtering on the initial band-pass filtered signal obtained in S2 to obtain a preliminary median filtered signal, and then the difference between the initial band-pass filtered signal and the preliminary median filtered signal is calculated to obtain a residual signal.S6. First, using the optimal initial window for signal-to-noise ratio obtained in S4 as the sliding unit, calculate the point-by-point correlation coefficient between the initial bandpass filtered signal obtained in S2 and the residual signal obtained in S5 to obtain the correlation coefficient corresponding to each data point in the initial bandpass filtered signal; S7. Based on the magnitude of the correlation coefficient of each data point obtained in S6, perform point-by-point median filtering on the initial bandpass filtered signal obtained in S2 again: if the correlation coefficient of the data point is greater than 0.99, then use the median filtering window adjusted based on the optimal initial window for noise attenuation obtained in S4 to perform median filtering on the data point; if the correlation coefficient of the data point is less than or equal to 0.99, then use the median filtering window adjusted based on the optimal initial window for effective signal retention obtained in S4 to perform median filtering on the data point. After processing each data point in the initial bandpass filtered signal, a secondary median filtered signal is obtained; S8. Perform bandpass filtering on the secondary median filtered signal obtained in S7 again to obtain the final calibration signal.

[0128] In application, taking the data recorded by two JCZ-360 ultra-wideband seismometers (wideband channel observation frequency band 50Hz-360s, data acquisition sampling rate 100Hz (sampling interval 0.01s)) installed at the Ezhuang Coal Mine as an example, the automatic calibration method for seismic observation data based on median filtering was applied to the two earthquake signals (M5.3 earthquake at 17:08:40 and M4.7 earthquake at 17:52:55) in the Hualien sea area on April 22, 2024. The specific process is as follows:

[0129] S1: Determine the time window (bt, et) for valid seismic signals

[0130] First, the surface and underground (-300 meters underground mine shaft) were examined. Figure 3 A preliminary analysis was conducted using the raw signals (denoted as data) from two seismometers (for comparison). Due to the large epicentral distance between the two earthquakes (belonging to teleseismic signals), the distance between the two instruments is negligible relative to the epicentral distance, and they are considered as co-located observations. This analysis was combined with the theoretical arrival times of the two earthquakes published by the China Earthquake Networks Center. Figure 3 (The vertical line marks the position), and the effective seismic signal time window in the original signal data is determined to be [bt, et] (where bt is set to 100 seconds before the theoretical arrival time of the M5.3 earthquake, and et is set to 100 seconds after the theoretical arrival time of the M4.7 earthquake, covering the pre-earthquake background noise, the entire propagation of the seismic phase, and the post-earthquake coma period of the two earthquakes).

[0131] S2: Extract the signal segment and perform bandpass filtering (to obtain bpdata).

[0132] Extract the signal segment in the [bt, et] time window from the original signal data. Since the target signal is a teleseismic signal, according to the frequency characteristics of the teleseismic signal (high proportion of low-frequency components, high-frequency interference is easy to remain), the frequency band of the band-pass filter is set to low frequency 40s (corresponding to frequency 0.025Hz) and high frequency 10Hz; if it is a local earthquake signal, the low frequency is set to 0.05Hz and the high frequency is set to 30Hz; if it is an ordinary earthquake signal, the low frequency is set to 0.05Hz and the high frequency is set to 10Hz. Perform band-pass filtering on the extracted signal segment to filter out low-frequency drift below 40s and high-frequency redundant noise above 10Hz, and obtain the initial band-pass filtered signal (denoted as bpdata, corresponding to Figure 3 “only band-pass filtered signal” in the middle).

[0133] S3: Multi-window median filtering and performance index calculation (SNR, SR, NR)

[0134] For bpdata, set the median filtering window parameters based on the seismometer sampling rate (100Hz): the initial filtering window is 3 sampling points (corresponding to 0.03s), the maximum filtering window is 1s (corresponding to 101 sampling points), and the length of all filtering windows is odd, and each time it is increased by 2 sampling points (i.e. the filtering window sequence is 3, 5, 7, …, 101 sampling points, a total of 50 different filtering windows).

[0135] Perform median filtering on bpdata with each filtering window in turn to obtain 50 filtered signals corresponding to the filtering windows; at the same time, calculate the performance index of each filtering window (the i-th window, i=1 corresponds to win=3) according to the following formula:

[0136] Signal-to-noise ratio: SNR(i)=rms(SG(i)) / rms(NS(i)), where SG(i) is the effective signal component after median filtering processing of the i-th filtering window, and NS(i) is the noise component after median filtering processing;

[0137] Effective signal retention rate: SR(i)=rms(SG(i)) / rms(SG0), where SG0 is the effective signal component of bpdata (initial signal before median filtering);

[0138] Noise attenuation rate: NR(i)=rms(NS(i)) / rms(NS0), where NS0 is the noise component of bpdata.

[0139] The calculation results are shown in Figure 2 from top to bottom, respectively, the signal-to-noise ratio change, the effective signal amplitude change, the noise amplitude change and the subsequent difference curve of the three under different filtering window lengths).

[0140] S4: Second-order difference, sliding standard deviation calculation and best initial window screening (winSNR, winSR, winNR)

[0141] Calculate the second-order difference: Calculate the second-order difference of adjacent data for the 50 groups of SNR(i), SR(i), NR(i) obtained in S3, to obtain the second-order difference sequence ddSNR of signal-to-noise ratio, the second-order difference sequence ddSR of effective signal retention rate, and the second-order difference sequence ddNR of noise attenuation rate, which are used to reflect the fluctuation trend of the three indicators.

[0142] Calculate the sliding standard deviation: Take "each second-order difference result and the next 4 second-order difference results" as the sliding window (a total of 5 consecutive difference results), and calculate the standard deviation of ddSNR, ddSR, and ddNR point by point to obtain the corresponding sliding standard deviation sequences σSNR, σSR, and σNR. The fluctuation degree of the second-order difference can be directly judged by the size of the standard deviation (the smaller the standard deviation, the more stable the indicator).

[0143] Screen the best initial window:

[0144] First, calculate the average value of the last 40% standard deviation results in σSNR, σSR, and σNR, denoted as mσSNR, mσSR, and mσNR, and then take twice their values as the stable threshold of the indicators changing from fluctuation to stability.

[0145] Take σSNR as an example, start searching from the first data in the sequence. If the consecutive three data (σSNR(i), σSNR(i+1), σSNR(i+2)) are all less than 2×mσSNR, then the median filtering window corresponding to σSNR(i) is the best initial window of signal-to-noise ratio winSNR. Similarly, the best initial windows of effective signal retention winSR and noise attenuation winNR are screened out.

[0146] Finally, winSNR=31 (corresponding to 31 sampling points, 0.31s), winSR=27 (0.27s), and winNR=37 (0.37s) are obtained (the results are derived from the stable segment of σSNR, σSR, and σNR). Figure 2

[0147] S5: Preliminary median filtering and residual signal calculation (mbpdata, deltdata)

[0148] ​The winSNR=31 is adopted to perform median filtering on the bpdata to obtain a preliminary median filtered signal mbpdata; and a residual signal deltdata is calculated according to the formula: deltdata=bpdata-mbpdata, the residual signal mainly contains noise which is not completely eliminated by the band-pass filtering and a small amount of effective signal details which are smoothed by the preliminary median filtering, and is used for judging the noise level of the original signal in the subsequent process.

[0149] S6: Point-by-point correlation coefficient calculation (CORR)

[0150] The winSNR=31 (31 sampling points) is adopted to perform point-by-point correlation coefficient calculation on the bpdata and the deltdata: around each data point in the bpdata, 31 continuous data points before and after the data point are taken to form a signal segment, and 31 continuous data points at the corresponding position of the deltdata are taken to form a residual segment, the correlation coefficient CORR(i) (i is the data point sequence number of the bpdata) of the two groups of segments is calculated, and finally a complete correlation coefficient sequence CORR is obtained.

[0151] According to the physical meaning, the deltdata is mainly noise, therefore, the greater the CORR(i) is, the stronger the correlation between the i-th data point and the noise is, that is, the greater the noise level of the data point is; otherwise, the noise level is small, and the proportion of the effective signal is high.

[0152] S7: Dynamic window adjustment and secondary median filtering (nmbpdata) in different regions

[0153] According to the numerical value of the CORR(i), different filtering windows are adopted to perform median filtering on each data point in the bpdata, and the specific rules are as follows:

[0154] If CORR(i)>0.99 (the noise level is high, and it is determined as a noise dense area): a filtering window of winNRx1.1 is adopted, that is, 37x1.1≈41 (an odd number of 41 sampling points, 0.41s), and the filtering window is expanded to strengthen the wrapping and elimination of the dense noise;

[0155] If CORR(i)≤0.99 (the noise level is low, and it is determined as an effective signal dense area): a filtering window of winSRx0.9 is adopted, that is, 27x0.9≈24 (an odd number of 25 sampling points, 0.25s), and the filtering window is reduced to reduce the amplitude loss of the effective signal.

[0156] After the above processing is completed on all data points in the bpdata, a secondary median filtered signal nmbpdata is obtained.

[0157] S8: Secondary band-pass filtering and generation of a final calibration signal (imdata)

[0158] The band-pass filtering (frequency band 40s-10Hz) consistent with S2 is performed on nmbpdata, and a small amount of redundant noise outside the frequency band (such as low-frequency fluctuation caused by filter window adjustment) possibly introduced in the secondary median filtering process is filtered out, and finally the calibration result of the surface seismic signal (denoted as imdata, corresponding to Figure 3 the "improved median filtered signal") is obtained.

[0159] From the processing effect: Figure 3 The signal-to-noise ratio of imdata (improved median filtered signal) is significantly higher than that of the original signal, the band-pass filtered signal and the conventional median filtered signal only, and the first arrival wave of the M5.3 earthquake can be clearly identified, and although the M4.7 earthquake is disturbed by the coda wave, the subsequent seismic phase can be effectively separated from the noise; Figure 4 The noise energy of imdata (improved median filtered signal) is significantly reduced compared with the conventional method, Figure 5 which shows that the effective signal amplitude retention rate is higher, and the waveform consistency with the underground control signal is good, which fully embodies the core advantages of the method "dynamic adaptive window, considering noise reduction and signal preservation".

[0160] The above only describes the preferred embodiments of the present application, and the protection scope of the present application is not limited to the above embodiments, but any equivalent modification or change made by those skilled in the art according to the disclosed content of the present application shall be included in the protection scope recited in the claims.

Claims

1. An automatic calibration method for seismic observation data based on median filtering, characterized in that: The method includes the following steps: S1. First, select the segment in the original signal where the amplitude variance is lower than a preset threshold as the background noise segment to calculate the background noise level. Then, detect the moment when the signal energy first continuously exceeds the background noise level to a first set threshold and record it as the trigger moment. Then, based on the earthquake occurrence time and epicentral distance provided by the seismic network, calibrate the trigger moment to determine it as the effective signal start point. Then, starting from the effective signal start point, track the signal energy backward. When the signal energy falls back to the background noise level range and continues to be lower than the second set threshold for a specified duration, record the end moment that meets the duration condition as the effective signal termination point. Finally, extend the effective signal start point forward by a first set duration and extend the effective signal termination point backward by a second set duration to finally form a seismic signal time window that includes the pre-earthquake background, complete seismic phase, and post-earthquake coma. S2. First, extract the signal segments within the time window of the seismic signal in S1, and then perform bandpass filtering on the signal segments according to the type of the original signal in S1 to obtain the initial bandpass filtered signal. S3. First, the initial bandpass filtered signal obtained in S2 is subjected to median filtering with median filtering windows of different window lengths in sequence to obtain the filtered signal corresponding to each median filtering window; then, the signal-to-noise ratio, effective signal retention rate, and noise attenuation rate corresponding to each median filtering window are calculated; the signal-to-noise ratio is obtained by calculating the root mean square amplitude ratio of the effective signal segment to the noise segment after median filtering; the effective signal retention rate is obtained by calculating the root mean square amplitude ratio of the energy of the effective signal segment before and after median filtering; the noise attenuation rate is obtained by calculating the root mean square amplitude ratio of the energy of the noise segment before and after median filtering. S4. First, calculate the difference sequence for each median filtering window in S3, considering its signal-to-noise ratio, effective signal retention rate, and noise attenuation rate: For the signal-to-noise ratio, calculate its first-order difference sequence and second-order difference sequence to obtain the first-order difference sequence and second-order difference sequence of the signal-to-noise ratio; For the effective signal retention rate and noise attenuation rate, the second-order difference sequences are calculated respectively to obtain the second-order difference sequence of the effective signal retention rate and the second-order difference sequence of the noise attenuation rate; Then, the moving standard deviations of the second-order difference sequences of signal-to-noise ratio, effective signal retention rate, and noise attenuation rate are calculated respectively to obtain the corresponding moving standard deviation sequences of signal-to-noise ratio, effective signal retention rate, and noise attenuation rate. Based on the moving standard deviation sequences of signal-to-noise ratio, effective signal retention rate, and noise attenuation rate, the optimal initial windows for signal-to-noise ratio, effective signal retention rate, and noise attenuation rate are preliminarily obtained. S5. First, using the optimal initial window for signal-to-noise ratio obtained in S4, perform median filtering on the initial bandpass filtered signal obtained in S2 to obtain a preliminary median filtered signal. Then, calculate the difference between the initial bandpass filtered signal and the preliminary median filtered signal to obtain the residual signal. S6. First, using the optimal initial window for signal-to-noise ratio obtained in S4 as the sliding unit, then perform point-by-point correlation coefficient calculation on the initial bandpass filtered signal obtained in S2 and the residual signal obtained in S5 to obtain the correlation coefficient corresponding to each data point in the initial bandpass filtered signal. S7. Based on the magnitude of the correlation coefficient of each data point obtained in S6, perform point-by-point median filtering on the initial bandpass filtered signal obtained in S2 again: If the correlation coefficient of the data point is greater than 0.99, then the median filter window adjusted based on the optimal initial window for noise attenuation obtained in S4 is used to perform median filtering on the data point; if the correlation coefficient of the data point is less than or equal to 0.99, then the median filter window adjusted based on the optimal initial window for effective signal retention obtained in S4 is used to perform median filtering on the data point. After processing each data point in the initial bandpass filtered signal, a secondary median filtered signal is obtained. S8. The secondary median filtered signal obtained in S7 is subjected to bandpass filtering again to obtain the final calibration signal.

2. The automatic calibration method for seismic observation data based on median filtering according to claim 1, characterized in that: In S2, the specific method for bandpass filtering the signal segment according to the type of the original signal in S1 is as follows: for the original signal types with different propagation characteristics, a corresponding bandpass filtering frequency band is preset, and the original signal types include at least near-field earthquake signals, far-field earthquake signals and ordinary earthquake signals.

3. The automatic calibration method for seismic observation data based on median filtering according to claim 2, characterized in that: The preset corresponding bandpass filter frequency band is set as follows: When the original signal type is a near-field signal, the bandpass filter frequency band is set from low frequency 0.05Hz to high frequency 30Hz; When the original signal type is a teleseismic signal, the bandpass filter frequency band is set from low frequency 0.025Hz to high frequency 10Hz; When the original signal type is a normal seismic signal, the bandpass filter frequency band is set from low frequency 0.05Hz to high frequency 10Hz.

4. The automatic calibration method for seismic observation data based on median filtering according to claim 1, characterized in that: In step S3, the parameters of the median filter windows of different lengths are determined in the following way: Based on the sampling rate of the seismic observation equipment and the high-frequency cutoff frequency of the bandpass filter band in S2, calculate the minimum window length, window step size, and maximum window length, where the minimum and maximum window lengths are both odd numbers.

5. The automatic calibration method for seismic observation data based on median filtering according to claim 4, characterized in that: The formula for calculating the minimum window length is: In the formula, bw is the minimum window length, f is the sampling rate of the seismic observation equipment, and h f is the high-frequency cutoff frequency of the bandpass filter in S2, and round is the rounding function; if the calculated result of bw is even, then the final value of bw is: the calculated result of bw plus 1; if the calculated result of bw is less than 3, then the final value of bw is set to 3. The formula for calculating the window step size is: In the formula, sw is the window step size, f is the sampling rate of the seismic observation equipment, and h f is the high-frequency cutoff frequency of the bandpass filter in S2, and round is the rounding function; if the calculated result of sw is odd, then the final value of sw is: the calculated result of sw plus 1; if the calculated result of sw is less than 2, then the final value of sw is set to 2. The formula for calculating the maximum window length is: ew = round(f × 1), In the formula, ew is the maximum window length, and f is the sampling rate of the seismic observation equipment; if the calculated result of ew is even, the final value of ew is: the calculated result of ew plus 1, and the final value of ew must be greater than or equal to 20 times the final value of bw plus 2.

6. The automatic calibration method for seismic observation data based on median filtering according to claim 1, characterized in that: In step S3, the formula for calculating the signal-to-noise ratio corresponding to each median filtering window is as follows: In the formula, SNR i Represents the signal-to-noise ratio corresponding to the i-th median filter window, rms represents the root mean square function, and SG... i NS represents the effective signal component in the filtered signal corresponding to the i-th median filter window. i This represents the noise component in the filtered signal corresponding to the i-th median filter window; The formula for calculating the effective signal retention rate for each median filter window is as follows: In the formula, SR i The median filter window represents the effective signal retention rate, rms represents the root mean square function, and SG represents the effective signal retention rate. i SG0 represents the effective signal component in the filtered signal corresponding to the i-th median filtering window, and SG0 represents the effective signal component in the initial bandpass filtered signal obtained in S2. The formula for calculating the noise attenuation rate for each median filter window is as follows: In the formula, NR i Represents the noise attenuation rate corresponding to the i-th median filter window, rms represents the root mean square function, and NS i NS0 represents the noise component in the filtered signal corresponding to the i-th median filtering window, and NS0 represents the noise component in the initial bandpass filtered signal obtained in S2.

7. The automatic calibration method for seismic observation data based on median filtering according to claim 1, characterized in that: In step S4, the specific method for calculating the moving standard deviation of the second-order difference sequences of signal-to-noise ratio, effective signal retention rate, and noise attenuation rate to obtain the moving standard deviation sequences of signal-to-noise ratio, effective signal retention rate, and noise attenuation rate is as follows: For the second-order difference sequence of signal-to-noise ratio, the second-order difference sequence of effective signal retention rate, and the second-order difference sequence of noise attenuation rate, each should be processed independently according to the following steps: Step 1: First, define a sliding window. Each sliding window consists of 5 consecutive second-order difference results in the second-order difference sequence. Specifically, it takes the k-th second-order difference result in the second-order difference sequence as the starting point and includes the k-th, k+1-th, k+2-th, k+3-th, and k+4-th second-order difference results. Step 2: Starting from the first second-order difference result of the second-order difference sequence, i.e., when k=1, the window contains the 1st to 5th second-order difference results. Calculate the standard deviation of the 5 second-order difference results within the sliding window to obtain one moving standard deviation data. Then, slide the sliding window forward one position at a time. When k=2, take the 2nd to 6th second-order difference results to form a new sliding window and calculate the standard deviation. Stop the calculation when the starting point k of the window satisfies that k+4 is greater than the total length of the second-order difference sequence, i.e., the sliding window cannot completely contain 5 second-order difference results. Step 3: Arrange all the calculated moving standard deviation data in order of their corresponding sliding window starting positions to obtain the moving standard deviation sequence corresponding to the second difference sequence.

8. The automatic calibration method for seismic observation data based on median filtering according to claim 7, characterized in that: In step S4, the specific method for filtering the moving standard deviation sequences based on the signal-to-noise ratio, effective signal retention rate, and noise attenuation rate to obtain the optimal initial window for signal-to-noise ratio, effective signal retention rate, and noise attenuation rate is as follows: For the signal-to-noise ratio moving standard deviation sequence, the effective signal retention rate moving standard deviation sequence, and the noise attenuation rate moving standard deviation sequence, each should be processed independently according to the following steps: The first step is to calculate the average of the last 40% of the moving standard deviation data in the moving standard deviation series, and then take twice the average as the stationarity threshold corresponding to the moving standard deviation series. The second step is to start from the first moving standard deviation data in the moving standard deviation sequence and check whether three consecutive moving standard deviation data are all less than the stationary threshold obtained in the first step, that is, to determine whether the j-th, j+1-th, and j+2-th moving standard deviation data are simultaneously less than the stationary threshold. The third step is to determine the median filter window corresponding to the first (j-th) moving standard deviation data among the three moving standard deviation data that first meet the condition that three consecutive moving standard deviation data are all less than the stationary threshold, which is found in the second step, as the optimal initial window for the index corresponding to the moving standard deviation sequence.

9. The automatic calibration method for seismic observation data based on median filtering according to claim 1, characterized in that: In step S7, if the correlation coefficient of the data point is greater than 0.99, then the median filter window adjusted based on the optimal initial window for noise attenuation obtained in step S4 is used to perform median filtering on the data point. Specifically, the length of the median filter window used to perform median filtering on the data point is adjusted to 1.1 times the length of the optimal initial window for noise attenuation. If the correlation coefficient of the data point is less than or equal to 0.99, then the median filter window adjusted based on the optimal initial window for retaining the effective signal obtained in S4 is used to perform median filtering on the data point. Specifically, the length of the median filter window used to perform median filtering on the data point is adjusted to 0.9 times the length of the optimal initial window for retaining the effective signal.

10. The automatic calibration method for seismic observation data based on median filtering according to claim 1, characterized in that: In step S4, after initially obtaining the optimal initial window for signal-to-noise ratio, the optimal initial window for retaining effective signal, and the optimal initial window for noise attenuation, it is also necessary to determine the first-order difference sequence of signal-to-noise ratio: If, in the first-order difference sequence of signal-to-noise ratio (SNR), there are 8 or more values ​​less than 0 out of 10 consecutive first-order difference results, then the median filtering window corresponding to the first first-order difference result less than 0 among these 10 consecutive first-order difference results is determined as the critical window where the SNR begins to decrease as the window length increases; if the window length of the initially obtained optimal initial window for SNR is greater than the window length of the critical window, then the final value of the window length of the optimal initial window for SNR is corrected to the window length of the critical window; if the window length of the initially obtained optimal initial window for effective signal retention is greater than the window length of the critical window, then the final value of the window length of the optimal initial window for effective signal retention is corrected to the window length of the critical window.

Images